% Stochastic constraint tightening for a scalar system.
A = 1;
B = 1;
Q = 0.5;
R = 1;
K = -dlqr(A, B, Q, R);
AK = A + B*K;
U = [-1, 1];
X = [-1, 1];
W = [-0.5, 0.5];
N = 10; % Number of time points over which to consider disturbance.
SKinf = W/(1 - AK);
Ubar = U - abs(K)*SKinf;
Xbar_robust = X - SKinf;
% Helper function for state evolution matrix.
E = evolution(AK, 1, N); % Converts from vector of w(k) to vector of e(k).
rand('state', 0);
sampleE = @(n) E*(W(1) + (W(2) - W(1))*rand(N, n));
Ntrials = 500; % How many times to perform the constraint tightening.
beta = 0.01;
epsilons = logspace(-3, 0, 16);
Ms = ceil(2./epsilons*(2 - log(beta)));
esim = sampleE(1000000); % Number of simulations to calculate epsilon_test.
Xbar = NaN(2, length(Ms), Ntrials);
eviolation = NaN(length(Ms), Ntrials); % Holds samples of epsilon_test.
for i = 1:length(Ms)
fprintf('%2d. epsilon = %10g (M = %d)\n', i, epsilons(i), Ms(i));
for j = 1:Ntrials
% Sample e sequence.
e = sampleE(Ms(i));
% Find min/max of e and check constraint violation.
emin = min([e(1,:), 0]);
emax = max([e(1,:), 0]);
eviolation(i, j) = mean(esim(1,:) < emin) + mean(esim(1,:) > emax);
% Tighten constraints on x.
Xbar(1, i, j) = min(X(1) - emin, 0);
Xbar(2, i, j) = max(X(2) - emax, 0);
end
end
% Plot tightened sets.
labels = {'Xbar_{lb}', 'Xbar_{ub}'};
locations = {'SouthWest', 'NorthWest'};
for i = 1:2
dum = violationquantiles(epsilons, Xbar(i,:,:), false());
plot(epsilons([1, end]), Xbar_robust(i)*[1, 1], '-k', ...
'DisplayName', 'Robust');
xlabel('\epsilon');
ylabel(labels{i});
legend('location', locations{i});
title('Tightened Constraints');
axis([min(epsilons), max(epsilons), i - 2.1, i - 0.9]);
end
% Plot violation probabilities.
q = violationquantiles(epsilons, eviolation, true());
plot(epsilons, epsilons, '--k');
legend('off');
legend('location', 'SouthEast');
ylabel('\epsilon_{test}');
xlabel('\epsilon');
axis([min(epsilons), max(epsilons), 1e-6, 1]);