1 | function pass = helmersson |
---|
2 | |
---|
3 | % this example is from scherer, lpv control and full block |
---|
4 | % multipliers, automatica 37 (2001), pp. 361-375. |
---|
5 | |
---|
6 | clear all |
---|
7 | |
---|
8 | if ~exist ('a', 'var'), |
---|
9 | a = 1.0; |
---|
10 | elseif size (a) ~= [1 1], |
---|
11 | a = 1.0; |
---|
12 | end |
---|
13 | |
---|
14 | sz = 3; |
---|
15 | if ~exist ('ndd', 'var'), |
---|
16 | ndd = sz; |
---|
17 | elseif size (ndd) ~= [1 1], |
---|
18 | ndd = sz; |
---|
19 | end |
---|
20 | |
---|
21 | if ~exist ('nddr', 'var'), |
---|
22 | nddr = sz; |
---|
23 | elseif size (nddr) ~= [1 1], |
---|
24 | nddr = sz; |
---|
25 | end |
---|
26 | |
---|
27 | ddr = 2.0/nddr; |
---|
28 | dd = ddr/ndd; |
---|
29 | |
---|
30 | if ~exist('slack_factor', 'var'), |
---|
31 | slack_factor = 1.0; |
---|
32 | elseif size (slack_factor) ~= [1 1], |
---|
33 | slack_factor = 1.0; |
---|
34 | end |
---|
35 | |
---|
36 | A = -1; |
---|
37 | B = [1 1 1 1 0 1]; |
---|
38 | C = [0; 0; 0; 0; 0; 1]; |
---|
39 | D = [0 1 0 1 0 0; |
---|
40 | 0.5 0 0.5 0 1 0; |
---|
41 | 2*a 0 a 0 0 0; |
---|
42 | 0 -2*a 0 -a 0 1; |
---|
43 | 0 1 0 0 0 0; |
---|
44 | 1 0 0 0 0 0]; |
---|
45 | |
---|
46 | I = eye(size([A B],2)); |
---|
47 | |
---|
48 | blk = [1 4 1 1]; % states deltas perf control |
---|
49 | ABCD = [A B; C D]; |
---|
50 | |
---|
51 | |
---|
52 | deltar = 0.9*(-1+0.5*ddr:ddr:1-0.5*ddr); |
---|
53 | |
---|
54 | Deltars = zeros (4, 4, length(deltar)^2); |
---|
55 | m = 0; |
---|
56 | nd = length (deltar); |
---|
57 | for i1 = 1:nd, |
---|
58 | d1 = deltar(i1); |
---|
59 | for i2 = 1:nd, |
---|
60 | d2 = deltar(i2); |
---|
61 | m = m + 1; |
---|
62 | Deltars(:,:,m) = blkdiag (d1*eye(2), d2*eye(2)); |
---|
63 | end |
---|
64 | end |
---|
65 | |
---|
66 | nr = size(Deltars,3); |
---|
67 | |
---|
68 | delta = 0.9*(-0.5*ddr:dd:0.5*ddr); |
---|
69 | |
---|
70 | maxsize = max ([0 diff(delta); diff(delta) 0]); |
---|
71 | Deltas = zeros (4, 4, length(delta)^2, nr); |
---|
72 | slack = slack_factor * ones (length(delta)^2, 2); |
---|
73 | |
---|
74 | for ir = 1:nr |
---|
75 | maxsize = max ([0 diff(delta); diff(delta) 0]); |
---|
76 | slack = slack_factor * ones (length(delta)^2, 2); |
---|
77 | |
---|
78 | m = 1; |
---|
79 | nd = length (delta); |
---|
80 | for i1 = 1:nd, |
---|
81 | d1 = delta(i1); |
---|
82 | for i2 = 1:nd, |
---|
83 | d2 = delta(i2); |
---|
84 | Deltas(:,:,m,ir) = blkdiag (d1*eye(2), d2*eye(2)) + Deltars(:,:,ir); |
---|
85 | slack(m,:) = slack_factor * [maxsize(i1) maxsize(i2)]; |
---|
86 | m = m + 1; |
---|
87 | end |
---|
88 | end |
---|
89 | |
---|
90 | DDeltas = []; |
---|
91 | DDeltas(:,:,1) = blkdiag (eye(2), 0*eye(2)); |
---|
92 | DDeltas(:,:,2) = blkdiag (0*eye(2), eye(2)); |
---|
93 | end |
---|
94 | |
---|
95 | [P, R, gam, F] = synlpv (ABCD, blk, Deltas, DDeltas, slack); |
---|
96 | |
---|
97 | % evaluation |
---|
98 | |
---|
99 | Pd = P(3:10,3:10); |
---|
100 | Rd = R(3:10,3:10); |
---|
101 | |
---|
102 | pass = norm(gam-1.527)<1e-2; |
---|
103 | |
---|
104 | |
---|
105 | function [P_feas, R_feas, gam_feas, F] = synlpv (ABCD, blk, Deltas, ... |
---|
106 | DDeltas, slack) |
---|
107 | |
---|
108 | if nargin < 3, |
---|
109 | disp ('usage: [P, R, gam] = synlpv (A, blk, Deltas, DDeltas, slack)'); |
---|
110 | return; |
---|
111 | end |
---|
112 | |
---|
113 | if nargin < 4, DDeltas = []; end |
---|
114 | |
---|
115 | n = sum(blk(1:3)); |
---|
116 | A = ABCD(1:n,1:n); |
---|
117 | B = ABCD(1:n,n+1:end); |
---|
118 | C = ABCD(n+1:end,1:n); |
---|
119 | |
---|
120 | n1 = blk(1); % states |
---|
121 | n2 = blk(2); % deltas |
---|
122 | n3 = blk(3); % H-inf |
---|
123 | [n, k] = size (A); |
---|
124 | [kd, nd, q, r] = size (Deltas); |
---|
125 | [kdd, ndd, qd] = size (DDeltas); |
---|
126 | if isempty (DDeltas), qd = 0; end |
---|
127 | |
---|
128 | if n ~= n1+n2+n3, error ('blk not compatible with A'); end |
---|
129 | if n ~= n1+n2+n3, error ('blk not compatible with A'); end |
---|
130 | if k ~= n1+n2+n3, error ('blk not compatible with A'); end |
---|
131 | |
---|
132 | if n2 ~= nd | n2 ~= kd, error ('blk not compatible with Deltas'); end |
---|
133 | |
---|
134 | if qd > 0 & (kdd ~= kd | ndd ~= nd), |
---|
135 | error ('DDeltas not compatible with Deltas'); |
---|
136 | end |
---|
137 | |
---|
138 | if nargin < 5 & qd > 1, |
---|
139 | slack = ones(q,qd); |
---|
140 | end |
---|
141 | |
---|
142 | if length (slack) ~= q, |
---|
143 | error ('slack not compatible with Deltas'); |
---|
144 | end |
---|
145 | |
---|
146 | slackvar = any (slack(:) > 0); |
---|
147 | |
---|
148 | ns = n1; |
---|
149 | ks = n1; % state block is always square |
---|
150 | |
---|
151 | np = n3; % performance block (H-inf) |
---|
152 | kp = n3; |
---|
153 | |
---|
154 | tic |
---|
155 | yalmip ('clear'); |
---|
156 | |
---|
157 | % define the block diagonal structure of P and R |
---|
158 | |
---|
159 | Ps = sdpvar (ns, ns, 'symmetric'); |
---|
160 | Rs = sdpvar (ns, ns, 'symmetric'); |
---|
161 | zs = zeros(ns,ns); |
---|
162 | |
---|
163 | gam = sdpvar (1, 1); |
---|
164 | gami = sdpvar (1, 1); |
---|
165 | |
---|
166 | for ir = 1:r, |
---|
167 | Pd(ir).var = sdpvar (nd+kd, nd+kd, 'symmetric'); |
---|
168 | P(ir).var = blkdiag ([zs Ps; Ps zs], Pd(ir).var); |
---|
169 | Rd(ir).var = sdpvar (nd+kd, nd+kd, 'symmetric'); |
---|
170 | R(ir).var = blkdiag ([zs Rs; Rs zs], Rd(ir).var); |
---|
171 | |
---|
172 | % Hinf outputs and inputs |
---|
173 | P(ir).var = blkdiag (P(ir).var, gami*eye(np), -gam*eye(kp)); |
---|
174 | R(ir).var = blkdiag (R(ir).var, gam*eye(np), -gami*eye(kp)); |
---|
175 | |
---|
176 | if slackvar, |
---|
177 | % define slack variables of Delta inequalities |
---|
178 | for iq = 1:qd, |
---|
179 | DDP(iq,ir).var = sdpvar (nd, nd, 'symmetric'); |
---|
180 | DDR(iq,ir).var = sdpvar (nd, nd, 'symmetric'); |
---|
181 | end |
---|
182 | end |
---|
183 | end |
---|
184 | |
---|
185 | Pswap = [1:ns n+(1:ks) ns+(1:nd) n+ns+(1:kd) ... |
---|
186 | ns+nd+(1:np) n+ns+nd+(1:kp)]; |
---|
187 | |
---|
188 | Rswap = [1:ks k+(1:ns) ks+(1:kd) k+ns+(1:nd) ... |
---|
189 | ks+kd+(1:kp) k+ns+nd+(1:np)]; |
---|
190 | |
---|
191 | AI = [A; eye(k)]; |
---|
192 | IA = [eye(n); -A']; |
---|
193 | |
---|
194 | ACp = AI(Pswap,:)*null(C); |
---|
195 | |
---|
196 | ABp = IA(Rswap,:)*null(B'); |
---|
197 | |
---|
198 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
199 | % gam*gami > 1 |
---|
200 | F = lmi; |
---|
201 | |
---|
202 | % [gam 1; 1 gami] > 0, |
---|
203 | % which is equivalent to gam > 0 and gami > 1/gam |
---|
204 | % even if gami > 1/gam we can replace it by 1/gam and still |
---|
205 | % satisfying the quadratic constraints since gami terms appears as |
---|
206 | % in the P and R inequalities as |
---|
207 | % ... + U' gami U < 0 and ... - V' gami V > 0, respectively |
---|
208 | F = addlmi (F, [gam 1; 1 gami]); |
---|
209 | |
---|
210 | % state-space multipliers defined by [Ps I; I Rs] > 0 |
---|
211 | F = addlmi (F, [Ps eye(ns); eye(ns) Rs]); |
---|
212 | |
---|
213 | for ir = 1:r |
---|
214 | % Cp' [A;I]' P(ir) [A;I] Cp < 0 |
---|
215 | F = addlmi (F, -ACp'*P(ir).var*ACp); |
---|
216 | |
---|
217 | % Bp' [I;A']' R(ir) [I;A'] Cp > 0 |
---|
218 | F = addlmi (F, ABp'*R(ir).var*ABp); |
---|
219 | |
---|
220 | if slackvar, |
---|
221 | % slack variables |
---|
222 | for iq = 1:qd, |
---|
223 | % DDP(iq,ir) > [0; DD]' Pd(ir) [0; DD] |
---|
224 | ZDD = [zeros(nd,nd); DDeltas(:,:,iq)]; |
---|
225 | |
---|
226 | F = addlmi (F, DDP(iq,ir).var - ZDD'*Pd(ir).var*ZDD); |
---|
227 | |
---|
228 | % DDP(iq,ir) > 0 |
---|
229 | F = addlmi (F, DDP(iq,ir).var); |
---|
230 | |
---|
231 | % DDR(iq,ir) + [-DD'; 0]' Rd(ir) [-DD'; 0] > 0 |
---|
232 | DDZ = [-DDeltas(:,:,iq)'; zeros(kd,kd)]; |
---|
233 | F = addlmi (F, DDR(iq,ir).var + DDZ'*Rd(ir).var*DDZ); |
---|
234 | |
---|
235 | % DDR(iq,ir) > 0 |
---|
236 | F = addlmi (F, DDR(iq,ir).var); |
---|
237 | end |
---|
238 | end |
---|
239 | |
---|
240 | for iq = 1:q, |
---|
241 | % [I; D]' Pd(ir) [I; D] > 0.25 * sum DDP |
---|
242 | sumDDP = 0; |
---|
243 | if slackvar, |
---|
244 | for iqd = 1:qd, |
---|
245 | sumDDP = sumDDP + slack(iq,iqd)^2*DDP(iqd,ir).var; |
---|
246 | end |
---|
247 | end |
---|
248 | |
---|
249 | ID = [eye(nd); Deltas(:,:,iq,ir)]; |
---|
250 | F = addlmi (F, ID'*Pd(ir).var*ID - 0.25*sumDDP); |
---|
251 | |
---|
252 | % [-D'; I]' Rd(ir) [-D'; I] + 0.25 * sum DDR < 0 |
---|
253 | sumDDR = 0; |
---|
254 | if slackvar, |
---|
255 | for iqd = 1:qd, |
---|
256 | sumDDR = sumDDR + slack(iq,iqd)^2*DDR(iqd,ir).var; |
---|
257 | end |
---|
258 | end |
---|
259 | DI = [-Deltas(:,:,iq,ir)'; eye(kd)]; |
---|
260 | F = addlmi (F, -DI'*Rd(ir).var*DI - 0.25*sumDDR); |
---|
261 | end |
---|
262 | end |
---|
263 | |
---|
264 | toc |
---|
265 | |
---|
266 | % minimize gam |
---|
267 | %[c,A,b,K] = sedumidata(F,gam) |
---|
268 | %save dummy c A b K |
---|
269 | %break |
---|
270 | sol = solvesdp (F,gam) |
---|
271 | |
---|
272 | if isempty (sol), |
---|
273 | disp ('no solution found'); |
---|
274 | P = []; |
---|
275 | R = []; |
---|
276 | else |
---|
277 | gam_feas = sdp2mat (gam, sol); |
---|
278 | gami_feas = sdp2mat (gami, sol); |
---|
279 | disp (sprintf ('gam = %8.3e', gam_feas)) |
---|
280 | disp (sprintf ('gam*gami - 1 = %8.3e', gam_feas*gami_feas - 1)) |
---|
281 | P_feas = sdp2mat (P(1).var, sol); |
---|
282 | R_feas = sdp2mat (R(1).var, sol); |
---|
283 | end |
---|
284 | |
---|
285 | |
---|
286 | |
---|