% 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');