Figure 3.4:

Minimum feasible \alpha for varying N. Note that we require \alpha \in [0, 1).

Code for Figure 3.4

Text of the GNU GPL.

main.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
% Example Constraint tightening for a linear system.
A = [1, 1; 0, 1];
B = [0; 1];
K = [-0.4, -1.2];

% State constraints. u is unconstrained.
xmax = 1;
umax = 1;
Z = struct();
[GH, Z.psi] = hyperrectangle([xmax*[-1; -1]; -umax], [xmax*[1; 1]; umax]);
Z.G = GH(:,1:2);
Z.H = GH(:,3);

% Disturbance set. Because all sets are convex, we only have to keep track of
% the vertices.
infnorm = @(x) max(abs(x(:)));

wmax = 0.1;
W = [wmax, wmax, -wmax, -wmax;
     wmax, -wmax, -wmax, wmax];
normW = infnorm(W);

KW = K*W;
normKW = infnorm(KW);

% For each value of N, see how far you need to squeeze.
Nmax = 25;
AkW = cell(Nmax + 1, 1);
alphaW = cell(Nmax + 1, 1);
normZbars = NaN(3, Nmax + 1);
Nkeep = [2, 3, 4, 5, 6, 7, 8, 10, 15];
Ak = eye(size(A));
alphas = NaN(Nmax + 1, 1);
thetas = cell(Nmax + 1, 1);
for n = 1:(Nmax + 1)
    AkW{n} = Ak*W;

    normAkW = infnorm(AkW{n});
    normKAkW = infnorm(K*AkW{n});

    alphas(n) = max(normAkW/normW, normKAkW/normKW);
    alphaW{n} = alphas(n)*W;

    if alphas(n) < 1
        thetas{n} = tightenconstraints(Z, A, B, K, wmax*[1; 1], n);
    else
        thetas{n} = NaN(size(Z.psi));
    end
    psibar = Z.psi - thetas{n}/(1 - alphas(n) + eps());
    normZbars(:,n) = psibar([1; 3; 5]);

    Ak = (A + B*K)*Ak;
end

figure();
semilogy(0:Nmax, alphas, '-ok', [1, Nmax], [1, 1], '--r');
xlabel('N')
ylabel('Minimum \alpha');

figure();
hold('on');
Nplots = ceil(sqrt(length(Nkeep)));
for i = 1:length(Nkeep)
    subplot(Nplots, Nplots, i);
    n = Nkeep(i) + 1;
    Vin = AkW{n};
    Vout = alphaW{n};
    plot(Vin(1,[1:end,1]), Vin(2,[1:end,1]), '-or', ...
         Vout(1,[1:end,1]), Vout(2,[1:end,1]), '-sb');
    title(sprintf('N = %d: \\alpha = %g', n, alphas(n)));
end

figure();
N = 0:Nmax;
plot(N, normZbars(1,:), '-or', N, normZbars(2,:), '-og', ...
     N, normZbars(3,:), '-ob');
xlabel('N')
ylabel('Bound');
legend('x_1', 'x_2', 'u');

% Save data.
data = struct('N', 0:Nmax, 'alpha', alphas, 'normZbar', normZbars);
save('-v7', 'calcalpha.mat', '-struct', 'data');

hyperrectangle.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
function [A, b] = hyperrectangle(lb, ub)
% [A, b] = hyperrectangle(lb, ub)
%
% Returns halfspace representation of hyperrectangle with bounds lb and ub.
% Any infinite or NaN bounds are ignored.

narginchk(2, 2);
if ~isvector(lb) || ~isvector(ub) || length(lb) ~= length(ub)
    error('Inputs must be vectors of the same size!');
end

A = kron(eye(length(lb)), [1; -1]);
b = reshape([ub(:)'; -lb(:)'], 2*length(lb), 1);

goodrows = ~isinf(b) & ~isnan(b);
A = A(goodrows,:);
b = b(goodrows);

end%function

tightenconstraints.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
function theta = tightenconstraints(Z, A, B, K, W, N)
% theta = tightenconstraints(Z, A, B, K, W, N)
%
% Tightens state and input constraints for a box W.
%
% Z should be a struct with fields G, H, and psi to define the feasible set Z as
%
%    Gx + Hu <= psi
%
% A, B, and K should be the system model and controller gain. W should be a
% vector defining the maximum absolute values for w. The system evolves as
%
%    x^+ = Ax + Bu + w,     u = v + Kx,     -W <= w <= W
%
% This means that W must be a symmetric box in R^n (we make this restriction
% because optimization becomes particularly easy).
%
% Returns theta such that the tighter constraints
%
%    Gx + H(v + Kx) <= psi - theta
%
% are satisfied for any N realizations of W.
narginchk(6, 6);
if ~isstruct(Z) || ~all(isfield(Z, {'G', 'H'}))
    error('Invalid input for Z!');
end

C = Z.G + Z.H*K;
A = A + B*K;
theta = zeros(size(size(C, 1)));
Ak = eye(size(A));
W = abs(W);
for i = 0:N
    theta = theta + abs(C*Ak)*W;
    Ak = A*Ak;
end

end%function