## Code for Figure 3.5

### 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
19function [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
38function 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