Skip to content

Commit e3dbd6b

Browse files
addat10addat10
committed
Adds the pointer to the submitted paper and some other comments in the analysis scripts.
1 parent 8eaa85b commit e3dbd6b

6 files changed

Lines changed: 169 additions & 85 deletions

‎analysis_scripts/bisection_exponent.m‎

Lines changed: 47 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,23 @@
1+
%---------------------------------------------------------------------------------------------------
2+
% For Paper
3+
% "Robust Performance Analysis of Source-Seeking Dynamics with Integral Quadratic Constraints"
4+
% by Adwait Datar and Herbert Werner
5+
% Copyright (c) Institute of Control Systems, Hamburg University of Technology. All rights reserved.
6+
% Licensed under the GPLv3. See LICENSE in the project root for license information.
7+
% Author(s): Adwait Datar
8+
%---------------------------------------------------------------------------------------------------
9+
110
function [alpha_best,P_ret]=bisection_exponent(G_veh,m,L,alpha_lims,tolerences,multiplier_class)
211
% This function runs a bisection algorithm for obtaining the best
312
% exponent alpha by calling the analysis routine for a fixed alpha
413
% multiple times.
5-
14+
% G_veh -> Plant model (typically of a vehicle)
15+
% m -> strong convexity parameter (lower bound)
16+
% L -> Lipschitz constant for gradients (upper bound)
17+
% alpha_lims -> Initial search space for bisection
18+
% tolerences -> tolerences for the optimization: bisection, LMI_tol, condition no tol
19+
% multiplier_class -> IQC multipliers such as Circle Criterion, Zames Falb, etc.
20+
621
% Get tolernces
722
cond_tol=tolerences.cond_tol;
823
cvx_tol=tolerences.cvx_tol;
@@ -16,21 +31,21 @@
1631
[status,P_ret]=verify_exp_stab_FBCC(G_veh,alpha_lims(1),m,L,cond_tol,cvx_tol);
1732
case 3 % Zames Falb Multipliers with CC
1833
[status,P_ret]=verify_exp_stab_ZF(G_veh,alpha_lims(1),m,L,cond_tol,cvx_tol);
19-
case 6
34+
case 6 % Zames Falb Multipliers with a basis for LTI systems
2035
rho=multiplier_class.rho;
2136
psi_order=multiplier_class.psi_order;
2237
odd_flag=multiplier_class.odd_flag;
2338
causal_flag=multiplier_class.causal_flag; % 1: causal, -1:anti-causal, 0:non-causal
2439
[status,P_ret]=verify_exp_stab_ZF_basis(G_veh,alpha_lims(1),m,L,odd_flag,causal_flag,rho,psi_order,cond_tol,cvx_tol);
25-
case 7
40+
case 7 % Zames Falb Multipliers with a basis for LPV systems
2641
rho=multiplier_class.rho;
2742
psi_order=multiplier_class.psi_order;
2843
odd_flag=multiplier_class.odd_flag;
2944
causal_flag=multiplier_class.causal_flag; % 1: causal, -1:anti-causal, 0:non-causal
3045
[status,P_ret]=verify_exp_stab_ZF_basis_LPV_example(G_veh,alpha_lims(1),m,L,odd_flag,causal_flag,rho,psi_order,cond_tol,cvx_tol);
3146
end
32-
33-
if ~status
47+
48+
if ~status % return -1 and stop if Infeasible for alpha_low
3449
%error('Infeasible for alpha_low. Choose a smaller alpha_low')
3550
alpha_best=-1;
3651
return
@@ -44,38 +59,37 @@
4459
[status,P]=verify_exp_stab_FBCC(G_veh,alpha_lims(2),m,L,cond_tol,cvx_tol);
4560
case 3 % Zames Falb Multipliers with CC
4661
[status,P]=verify_exp_stab_ZF(G_veh,alpha_lims(2),m,L,cond_tol,cvx_tol);
47-
case 6
62+
case 6 % Zames Falb Multipliers with a basis for LTI systems
4863
[status,P]=verify_exp_stab_ZF_basis(G_veh,alpha_lims(2),m,L,odd_flag,causal_flag,rho,psi_order,cond_tol,cvx_tol);
49-
case 7
64+
case 7 % Zames Falb Multipliers with a basis for LPV systems
5065
[status,P]=verify_exp_stab_ZF_basis_LPV_example(G_veh,alpha_lims(2),m,L,odd_flag,causal_flag,rho,psi_order,cond_tol,cvx_tol);
51-
end
52-
if status
66+
end
67+
if status
5368
alpha_best=alpha_lims(2); % Return alpha_high if feasible
5469
P_ret=P;
55-
else
56-
57-
% Start the bisection
58-
while alpha_lims(2)-alpha_lims(1)>bisect_tol
59-
alpha_mid=mean(alpha_lims);
60-
switch multiplier_class.id
61-
case 1 % Circle Criterion
62-
[status,P]=verify_exp_stab_CC(G_veh,alpha_mid,m,L,cond_tol,cvx_tol);
63-
case 2 % Full block circle criterion
64-
[status,P]=verify_exp_stab_FBCC(G_veh,alpha_mid,m,L,cond_tol,cvx_tol);
65-
case 3 % Zames Falb Multipliers
66-
[status,P]=verify_exp_stab_ZF(G_veh,alpha_mid,m,L,cond_tol,cvx_tol);
67-
case 6
68-
[status,P]=verify_exp_stab_ZF_basis(G_veh,alpha_mid,m,L,odd_flag,causal_flag,rho,psi_order,cond_tol,cvx_tol);
69-
case 7
70-
[status,P]=verify_exp_stab_ZF_basis_LPV_example(G_veh,alpha_mid,m,L,odd_flag,causal_flag,rho,psi_order,cond_tol,cvx_tol);
71-
end
72-
if status
73-
alpha_lims(1)=alpha_mid;
74-
P_ret=P;
75-
else
76-
alpha_lims(2)=alpha_mid;
70+
else
71+
% Start the bisection and repeat until the limits are within bisect_tol
72+
while alpha_lims(2)-alpha_lims(1)>bisect_tol
73+
alpha_mid=mean(alpha_lims); % bisect the current limits
74+
switch multiplier_class.id
75+
case 1 % Circle Criterion
76+
[status,P]=verify_exp_stab_CC(G_veh,alpha_mid,m,L,cond_tol,cvx_tol);
77+
case 2 % Full block circle criterion
78+
[status,P]=verify_exp_stab_FBCC(G_veh,alpha_mid,m,L,cond_tol,cvx_tol);
79+
case 3 % Zames Falb Multipliers
80+
[status,P]=verify_exp_stab_ZF(G_veh,alpha_mid,m,L,cond_tol,cvx_tol);
81+
case 6 % Zames Falb Multipliers with a basis for LTI systems
82+
[status,P]=verify_exp_stab_ZF_basis(G_veh,alpha_mid,m,L,odd_flag,causal_flag,rho,psi_order,cond_tol,cvx_tol);
83+
case 7 % Zames Falb Multipliers with a basis for LPV systems
84+
[status,P]=verify_exp_stab_ZF_basis_LPV_example(G_veh,alpha_mid,m,L,odd_flag,causal_flag,rho,psi_order,cond_tol,cvx_tol);
85+
end
86+
if status % Choose the upper half as the new limits if feasible
87+
alpha_lims(1)=alpha_mid;
88+
P_ret=P;
89+
else % Choose the lower half as the new limits if infeasible
90+
alpha_lims(2)=alpha_mid;
91+
end
7792
end
78-
end
79-
alpha_best=alpha_lims(1);
93+
alpha_best=alpha_lims(1);
8094
end
8195
end

‎analysis_scripts/verify_exp_stab_CC.m‎

Lines changed: 13 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,14 @@
1+
%---------------------------------------------------------------------------------------------------
2+
% For Paper
3+
% "Robust Performance Analysis of Source-Seeking Dynamics with Integral Quadratic Constraints"
4+
% by Adwait Datar and Herbert Werner
5+
% Copyright (c) Institute of Control Systems, Hamburg University of Technology. All rights reserved.
6+
% Licensed under the GPLv3. See LICENSE in the project root for license information.
7+
% Author(s): Adwait Datar
8+
%---------------------------------------------------------------------------------------------------
19
function [status,P]=verify_exp_stab_CC(G_veh,alpha,sec_1,sec_2,cond_tol,tol)
2-
% This function runs the analysis LMI with cvx and return the status and
3-
% the storage function matrix P.
10+
% This function runs the analysis LMI with cvx using the circle criterion and return the status and
11+
% the storage matrix P.
412

513
[A_G,B_G,C_G,D_G]=ssdata(G_veh);
614

@@ -29,7 +37,9 @@
2937

3038
end
3139
function [status,P]=verify_exp_stab_alpha_CC_LMI(PSI_GI,M,alpha,cond_tol,tol)
32-
% This function runs the exp-stab analysis KYP Lemma LMI
40+
% This function runs the exp-stab analysis KYP Lemma LMI for the static
41+
% multiplier M corresponding to the circle criterion
42+
3343
status=false;
3444
% Get state-space matrics
3545
[A,B,C,D]=ssdata(PSI_GI);

‎analysis_scripts/verify_exp_stab_FBCC.m‎

Lines changed: 18 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,14 @@
1+
%---------------------------------------------------------------------------------------------------
2+
% For Paper
3+
% "Robust Performance Analysis of Source-Seeking Dynamics with Integral Quadratic Constraints"
4+
% by Adwait Datar and Herbert Werner
5+
% Copyright (c) Institute of Control Systems, Hamburg University of Technology. All rights reserved.
6+
% Licensed under the GPLv3. See LICENSE in the project root for license information.
7+
% Author(s): Adwait Datar
8+
%---------------------------------------------------------------------------------------------------
19
function [status,X]=verify_exp_stab_FBCC(G_veh,alpha,sec_1,sec_2,cond_tol,tol)
2-
% This function runs the analysis LMI with cvx and returns the status and
3-
% the storage function matrix P.
10+
% This function runs the analysis LMI with cvx using the full block circle criterion and returns
11+
% the status and the storage function matrix P.
412

513
[A,B,C,D]=ssdata(G_veh);
614
[~,nu]=size(D);
@@ -15,7 +23,12 @@
1523
[status,X]=verify_exp_stab_FBCC_LMI(Psi_GI,alpha,sec_1,sec_2,cond_tol,tol);
1624
end
1725
function [status,X]=verify_exp_stab_FBCC_LMI(Psi_GI,alpha,sec_1,sec_2,cond_tol,tol)
18-
% This function runs the exp-stab analysis KYP Lemma LMI
26+
% This function runs the analysis LMI with cvx using the full block circle criterion and return
27+
% the status and the storage matrix P.
28+
% Refer to the following paper for details:
29+
% Scherer, C., 2021. Dissipativity and Integral Quadratic Constraints, Tailored computational
30+
% robustness tests for complex interconnections. arXiv preprint arXiv:2105.07401.
31+
1932
status=false;
2033
[A,B,C,D]=ssdata(Psi_GI);
2134
n=size(A,1);
@@ -31,12 +44,11 @@ variable S(dim,dim)
3144
variable R(dim,dim) symmetric
3245

3346
% Multiplier class
34-
% Create the multiplier class
47+
% Create the vertices of the convex polytope and verify the LMI there
3548
dec=0:1:(2^dim)-1;
3649
vec=dec2bin(dec)-'0';
3750
d=sec_1*(~vec)+sec_2*(vec);
38-
E=@(delta)[eye(dim);diag(delta)];
39-
51+
E=@(delta)[eye(dim);diag(delta)];
4052
C1=E(d(1,:))'*[Q,S;S',R]*E(d(1,:));
4153
C2=E(d(2,:))'*[Q,S;S',R]*E(d(2,:));
4254
C3=E(d(3,:))'*[Q,S;S',R]*E(d(3,:));

‎analysis_scripts/verify_exp_stab_ZF.m‎

Lines changed: 10 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,14 @@
1+
%---------------------------------------------------------------------------------------------------
2+
% For Paper
3+
% "Robust Performance Analysis of Source-Seeking Dynamics with Integral Quadratic Constraints"
4+
% by Adwait Datar and Herbert Werner
5+
% Copyright (c) Institute of Control Systems, Hamburg University of Technology. All rights reserved.
6+
% Licensed under the GPLv3. See LICENSE in the project root for license information.
7+
% Author(s): Adwait Datar
8+
%---------------------------------------------------------------------------------------------------
19
function [status,P]=verify_exp_stab_ZF(G_veh,alpha,sec_1,sec_2,cond_tol,tol)
2-
% This function runs the analysis LMI with cvx and returns the status and
3-
% the storage function matrix P.
10+
% This function runs the analysis LMI with cvx using causal Zames Falb multipliers and returns the
11+
% status and the storage function matrix P.
412

513
[A_G,B_G,C_G,D_G]=ssdata(G_veh);
614

‎analysis_scripts/verify_exp_stab_ZF_basis.m‎

Lines changed: 42 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,17 @@
1+
%---------------------------------------------------------------------------------------------------
2+
% For Paper
3+
% "Robust Performance Analysis of Source-Seeking Dynamics with Integral Quadratic Constraints"
4+
% by Adwait Datar and Herbert Werner
5+
% Copyright (c) Institute of Control Systems, Hamburg University of Technology. All rights reserved.
6+
% Licensed under the GPLv3. See LICENSE in the project root for license information.
7+
% Author(s): Adwait Datar
8+
%---------------------------------------------------------------------------------------------------
19
function [status,P]=verify_exp_stab_ZF_basis(G_veh,alpha,sec_1,sec_2,odd_flag,causal_flag,rho,n_psi,cond_tol,tol)
2-
% This function runs the analysis LMI with cvx and returns the status and
3-
% the storage function matrix P with Zames Falb Multipliers
10+
% This function runs the analysis LMI with cvx using general Zames Falb multipliers and returns the
11+
% status and the storage function matrix P
12+
% Refer to the following paper for details on the parameterization:
13+
% [1] Veenman, J., Scherer, C.W. and Köroğlu, H., 2016. Robust stability and performance analysis based on
14+
% integral quadratic constraints. European Journal of Control, 31, pp.1-32.
415

516
[A_G,B_G,C_G,D_G]=ssdata(G_veh);
617

@@ -14,7 +25,7 @@
1425
end
1526
n_delta=nu;
1627

17-
% Define the multiplier psi as in Scherer and Veenman
28+
% Define the multiplier psi similar to [1]
1829
if n_psi>=1
1930
A_psi_alpha=(rho-2*alpha)*diag(ones(1,n_psi))+diag(ones(1,n_psi-1),-1);
2031
B_psi_alpha=[1;zeros(n_psi-1,1)];
@@ -27,7 +38,7 @@
2738
end
2839

2940

30-
% Define the system psi_tilde for the positivity constraint
41+
% Define the system psi_tilde for the positivity constraint as in [1]
3142
psi_b_til=[];
3243
if n_psi>=1
3344
s=tf('s');
@@ -60,7 +71,7 @@
6071
[n,nu]=size(B);
6172
n_psi=size(Av,1);
6273
n_psi_ch=n_psi/nu;
63-
% Get R_nu as in Scherer and Veenman
74+
% Get R_nu as in [1]
6475
Rv=diag(sqrt(factorial((0:n_psi_ch-1))));
6576

6677
if n_psi_ch>=1
@@ -72,30 +83,34 @@
7283
cvx_begin sdp
7384
cvx_precision high
7485

86+
% The variables P2 and P4 are not required if the uncertainty is not
87+
% odd but are declared here to treat the more general case
7588
variable X(n,n) symmetric
7689
variable P0
7790
variable P1(1,n_psi_ch)
7891
variable P2(1,n_psi_ch)
7992
variable P3(1,n_psi_ch)
8093
variable P4(1,n_psi_ch)
81-
82-
if n_psi_ch>=2
94+
95+
% Need more variables for the positivity constraint for multipliers of
96+
% order higher than 1
97+
if n_psi_ch>=2
8398
variable X1(n_psi_til,n_psi_til) symmetric
8499
variable X2(n_psi_til,n_psi_til) symmetric
85100
variable X3(n_psi_til,n_psi_til) symmetric
86101
variable X4(n_psi_til,n_psi_til) symmetric
87102
end
88-
103+
89104
P_12=[ P0, P4-P3;
90105
P2'-P1', zeros(n_psi_ch)];
91106

107+
% Form the KYP Lemma LMI
92108
L1=[A'*X + X*A + 2*alpha*X, X*B;
93109
B'*X, zeros(nu)];
94110
L2=[C';D']*[zeros(n_psi+nu),kron(P_12,eye(nu));kron(P_12',eye(nu)),zeros(n_psi+nu)]*[C,D];
95-
96111
LMI=L1+L2;
97112

98-
if n_psi_ch>=2
113+
if n_psi_ch>=2 % Positivity constraint LMIs for multipliers of order higher than 1
99114
MAT=[eye(n_psi_til), zeros(n_psi_til,1);...
100115
Av_til, Bv_til;...
101116
Rv*Cv_til, Rv*Dv_til];
@@ -105,7 +120,7 @@ variable X4(n_psi_til,n_psi_til) symmetric
105120
X, zeros(n_psi_til), zeros(n_psi_til,n_psi_ch);...
106121
zeros(n_psi_ch,n_psi_til), zeros(n_psi_ch,n_psi_til), diag(P)]*...
107122
MAT);
108-
elseif n_psi_ch==1
123+
elseif n_psi_ch==1 % Positivity constraint LMIs for multipliers of order 1
109124
MAT=[Rv*Dv_til];
110125
pos_Lemma=@(P)(MAT'*diag(P)*MAT);
111126
end
@@ -117,7 +132,7 @@ variable X4(n_psi_til,n_psi_til) symmetric
117132
X <= cond_tol*tol*eye(n);
118133
LMI<=-tol*eye(n+nu);
119134

120-
% Norm constraint
135+
% Norm constraint on the multiplier
121136
if n_psi_ch>=1
122137
kron(eye(nu),P0) + kron((P1+P2+P3+P4),eye(nu))*inv(Av)*Bv >=tol*eye(nu);
123138
elseif n_psi_ch==0
@@ -127,12 +142,12 @@ variable X4(n_psi_til,n_psi_til) symmetric
127142
% Positivity constraint for the two cases of odd and non-odd uncertainty
128143
% Case 1: If the uncertainty is odd
129144
if odd_flag==1
130-
if n_psi_ch==1
145+
if n_psi_ch==1 % Add positivity constraint for multipliers of order 1
131146
pos_Lemma(P1)>=tol*eye(n_psi_ch);
132147
pos_Lemma(P3)>=tol*eye(n_psi_ch);
133148
pos_Lemma(P2)>=tol*eye(n_psi_ch);
134149
pos_Lemma(P4)>=tol*eye(n_psi_ch);
135-
elseif n_psi_ch>=2
150+
elseif n_psi_ch>=2 % Add positivity constraints for multipliers of order higher than 1
136151
X1>=tol*eye(n_psi_til)
137152
pos_Lemma(X1,P1)>=tol*eye(n_psi_til+1);
138153
X3>=tol*eye(n_psi_til)
@@ -142,32 +157,38 @@ variable X4(n_psi_til,n_psi_til) symmetric
142157
X4>=tol*eye(n_psi_til)
143158
pos_Lemma(X4,P4)>=tol*eye(n_psi_til+1);
144159
end
145-
if causal_flag==1
146-
% Impose these constraints if only searching for causal multipliers
160+
if causal_flag==1 % Causality of multiplier is enforced by setting P1=P2
147161
P1==P2;
148162
end
149-
if causal_flag==-1
150-
% Impose these constraints if only searching for anti-causal multipliers
163+
if causal_flag==-1 % Anti-causality of multiplier is enforced by setting P3=P4
151164
P3==P4;
152165
end
153166
end
154167
% Case 2: if the uncertainty is not odd
155-
if odd_flag==0
156-
if causal_flag==1
168+
if odd_flag==0 % Set the variables P2,P4 to zero
169+
% Since Causaility is enforced by by setting P1=P2 and since
170+
% P2=P4=0, the only variable is P3
171+
if causal_flag==1
157172
P2==zeros(1,n_psi_ch);
158173
P4==zeros(1,n_psi_ch);
159174
P1==P2;
160-
if n_psi_ch==1
175+
% Constraint on P3 to ensure positivity of multipliers of
176+
% order 1 or higher
177+
if n_psi_ch==1
161178
pos_Lemma(P3)>=tol*eye(n_psi_ch);
162179
elseif n_psi_ch>=2
163180
X3>=tol*eye(n_psi_til)
164181
pos_Lemma(X3,P3)>=tol*eye(n_psi_til+1);
165182
end
166183
end
184+
% Since Anti-causaility is enforced by by setting P3=P4 and since
185+
% P2=P4=0, the only variable is P1
167186
if causal_flag==-1
168187
P2==zeros(1,n_psi_ch);
169188
P4==zeros(1,n_psi_ch);
170189
P3==P4;
190+
% Constraint on P1 to ensure positivity of multipliers of
191+
% order 1 or higher
171192
if n_psi_ch==1
172193
pos_Lemma(P1)>=tol*eye(n_psi_ch);
173194
elseif n_psi_ch>=2

0 commit comments

Comments
 (0)