-
Notifications
You must be signed in to change notification settings - Fork 2
/
test_perf.m
163 lines (134 loc) · 5.99 KB
/
test_perf.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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% File : test_perf.m %
% %
% Author : Tobias Holicki %
% Date : 13.04.2023 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This script realizes the example in Section 2.5.2 of [1].
% It deals with the H_infty performance analysis of some interconnection
% inspired by a weighted tracking configuration where only a small part is
% affected by impulses. Note that the whole IQC approach in this context is
% motivated by the idea that usually not all of the system matrices are
% affected by impulses.
% We compare the numerically determined upper bounds on the robust energy
% gain. These are computed with three different tests:
% - A variation of the clock based approach from [2].
% - Another variation of the clock based approach from [2] involving slack
% variables that are restricted to be constant. Such conditions appear
% for example when designing clock-independent static state-feedback
% controllers.
% - A variation of the IQC theorem from [3] where we view the impulse as
% generated by an uncertain operator and where we constructed an IQC
% for this operator based on the lifting approach.
%
% [1] T. Holicki, C. W. Scherer, IQC Based Analysis and Estimator Design
% for Discrete-Time Systems Affected by Impulsive Uncertainties, 2023
% [2] C. Briat, Convex conditions for robust stability analysis and
% stabilization of linear aperiodic impulsive and sampled-data systems
% under dwell-time constraints, 2013
% [3] C. W. Scherer, J. Veenman, Stability analysis by dynamic dissipation
% inequalities: On merging frequency-domain techniques with time domain
% conditions, 2018
% Clean up
clc
clear
% Addpath for auxiliary functions
addpath(genpath('AuxiliaryFunctions'));
%% System data
% Parameters involved in the satellite model
J1 = 1;
J2 = 0.1;
k = 0.091;
b = 0.04;
% Two-degree-of-freedom design configuration
A = [0, 1, 0, 0; [-k, -b, k, b]/J2; 0, 0, 0, 1; [k, b, -k, -b]/J1];
Bnru = [[0; 1; 0; 0], [0; 0; 0; 0], [0; 0; 0; 1/J1]];
Ce = [1, 0, 0, 0; 0, 0, 0, 0];
Denru = [[0; 0], [-1; 0], [0; 1]];
Cy = [1, 0, 0, 0; 0, 0, 0, 0];
Dynru = [[0; 0], [0; 1], [0; 0]];
sys = ss(A, Bnru, [Ce; Cy], [Denru; Dynru], ...
'StateName', {'theta2', 'd theta2', 'theta1', 'd theta1'}, ...
'InputName', {'n', 'r', 'u'}, ...
'OutputName', {'y-r', 'u', 'y', 'r'});
inp = [2, 1]; % [gen. dist, control]
out = [2, 2]; % [error, measurements]
% Discretize
Ts = 0.01; % Sampling time
sysd = c2d(sys, Ts);
% Performance weights
s = zpk('s');
w1 = c2d((0.5*s +0.433) / (s + 0.00433), Ts);
w2 = 0.1;
% Weighted open-loop configuration
sysw = blkdiag(w1, w2, eye(2)) * sysd * blkdiag(0.01, 1, 1);
% Design an Hinfty controller for the nominal weighted generalized plant
[K, ~, ga] = hinfsyn(sysw, out(2), inp(2));
disp(['Achieved performance for nominal plant: ', num2str(ga)]);
%% Limited information setting
% We assume that the output y is no longer available at all times, but r
% still is. The previously obtained controller is then driven by a
% piecewise constant version of y. This is achieved by incorporating a hold
% device which can be modeled as an impulsive system.
% Naturally, the controller will have a poor performance in this case, but
% we are only interested in analyzing it. See for example the work by
% Leonid Mirkin for reasonable approaches to design controllers or modify
% the given one in this context.
% Unweighted open loop with impulsive operator (discretized)
A = blkdiag(sysd.a, 1);
Bw = [zeros(4, 1); 1];
Bnru = [sysd.b; zeros(1, 3)];
Cz = [Cy(1, :), -1];
Dzwnru = zeros(1, 1+3);
Ce = [Ce, [0; 0]];
Dewnru = [zeros(2, 1), Denru];
Cy = blkodiag(1, Cy(2, :));
Dywnru = [zeros(2, 1), Dynru];
sysI = ss(A, [Bw, Bnru], [Cz; Ce; Cy], [Dzwnru; Dewnru; Dywnru], Ts, ...
'StateName', {'theta2', 'd theta2', 'theta1', 'd theta1', 'y'},...
'InputName', {'w'; 'n'; 'r'; 'u'}, ...
'OutputName', {'z'; 'y-r'; 'u'; 'ty'; 'r'});
% Weighted open-loop configuration with same weights as before
sysw = blkdiag(1, w1, w2, eye(2)) * sysI * blkdiag(1, 0.01, 1, 1);
% Close loop with previsouly obtained controller
clw = lft(sysw, K);
% Standard impulsive description
sysF = lft(0, clw);
sysJ = lft(1, clw);
% Dwell-time samples
T = {[5, 5], [5, 6], [5, 7], [5, 8], [5, 9], [6, 6], [6, 7], [6, 8], ...
[6, 9], [7, 7], [7, 8], [7, 9]};
% Initialization of upper bounds on energy gain
ga = zeros(6, length(T));
% Compute upper bounds on the energy gain for different dwell-times with
% different algorithms
for t = 1 : length(T)
% IQC based analysis
for nu = 1 : 4
psi = basis_filter(nu, 2, sampling_time=sysI.Ts, pole=0, type=2);
ga(nu, t) = ana_perf_iqc(clw, T{t}, 1, psi=psi);
end
% Clock based analysis
ga(5, t) = ana_perf_clock(sysF, sysJ, T{t});
% Clock and slack variable analysis
ga(6, t) = ana_perf_clock_sv(sysF, sysJ, T{t});
end
%% Plot results
figure
hold on; box on
plot(1:length(T), ga(1, :), 'x', 'LineWidth', 1.8, 'MarkerSize', 10)
plot(1:length(T), ga(2, :), 'x', 'LineWidth', 1.8, 'MarkerSize', 10)
plot(1:length(T), ga(3, :), 'x', 'LineWidth', 1.8, 'MarkerSize', 10)
plot(1:length(T), ga(4, :), 'x', 'LineWidth', 1.8, 'MarkerSize', 10)
plot(1:length(T), ga(5, :), 'o', 'LineWidth', 1.5, 'MarkerSize', 8)
plot(1:length(T), ga(6, :), 'd', 'LineWidth', 1.5, 'MarkerSize', 8)
xticks(1:length(T));
tostr = @(T) sprintf('(%i, %i)', T(1), T(2));
xticklabels(cellfun(tostr, T, 'UniformOutput', false))
xlim([0.5, length(T)+0.5])
xlabel('(T_{min}, T_{max})')
ylabel('Energy gain upper bound')
legend('Thm 2.16, \nu=1', 'Thm 2.16, \nu=2', 'Thm 2.16, \nu=3', ...
'Thm 2.16, \nu=4', 'Theorem 2.6', 'Theorem 2.8', 'Location', ...
'northwest')