|
3430
|
1 ## Copyright (C) 1996, 1998 Auburn University. All rights reserved. |
|
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
7016
|
6 ## under the terms of the GNU General Public License as published by |
|
|
7 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
8 ## your option) any later version. |
|
3430
|
9 ## |
|
7016
|
10 ## Octave is distributed in the hope that it will be useful, but |
|
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
13 ## General Public License for more details. |
|
3430
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
7016
|
16 ## along with Octave; see the file COPYING. If not, see |
|
|
17 ## <http://www.gnu.org/licenses/>. |
|
3430
|
18 |
|
|
19 ## -*- texinfo -*- |
|
3502
|
20 ## @deftypefn {Function File} {[@var{k}, @var{g}, @var{gw}, @var{xinf}, @var{yinf}] =} hinfsyn (@var{asys}, @var{nu}, @var{ny}, @var{gmin}, @var{gmax}, @var{gtol}, @var{ptol}, @var{tol}) |
|
3430
|
21 ## |
|
|
22 ## @strong{Inputs} input system is passed as either |
|
|
23 ## @table @var |
|
3502
|
24 ## @item asys |
|
5016
|
25 ## system data structure (see @command{ss}, @command{sys2ss}) |
|
3430
|
26 ## @itemize @bullet |
|
|
27 ## @item controller is implemented for continuous time systems |
|
5016
|
28 ## @item controller is @strong{not} implemented for discrete time systems (see |
|
|
29 ## bilinear transforms in @command{c2d}, @command{d2c}) |
|
3430
|
30 ## @end itemize |
|
|
31 ## @item nu |
|
|
32 ## number of controlled inputs |
|
|
33 ## @item ny |
|
|
34 ## number of measured outputs |
|
|
35 ## @item gmin |
|
5016
|
36 ## initial lower bound on |
|
|
37 ## @iftex |
|
|
38 ## @tex |
|
|
39 ## $ { \cal H }_\infty $ |
|
|
40 ## @end tex |
|
|
41 ## @end iftex |
|
|
42 ## @ifinfo |
|
|
43 ## H-infinity |
|
|
44 ## @end ifinfo |
|
|
45 ## optimal gain |
|
3430
|
46 ## @item gmax |
|
5016
|
47 ## initial upper bound on |
|
|
48 ## @iftex |
|
|
49 ## @tex |
|
|
50 ## $ { \cal H }_\infty $ |
|
|
51 ## @end tex |
|
|
52 ## @end iftex |
|
|
53 ## @ifinfo |
|
|
54 ## H-infinity |
|
|
55 ## @end ifinfo |
|
|
56 ## Optimal gain. |
|
3430
|
57 ## @item gtol |
|
5016
|
58 ## Gain threshold. Routine quits when @var{gmax}/@var{gmin} < 1+tol. |
|
3430
|
59 ## @item ptol |
|
5016
|
60 ## poles with @code{abs(real(pole))} |
|
|
61 ## @iftex |
|
|
62 ## @tex |
|
|
63 ## $ < ptol \Vert H \Vert $ |
|
|
64 ## @end tex |
|
|
65 ## @end iftex |
|
|
66 ## @ifinfo |
|
|
67 ## < ptol*||H|| |
|
|
68 ## @end ifinfo |
|
|
69 ## (@var{H} is appropriate |
|
3430
|
70 ## Hamiltonian) are considered to be on the imaginary axis. |
|
5016
|
71 ## Default: 1e-9. |
|
3430
|
72 ## @item tol |
|
5016
|
73 ## threshold for 0. Default: 200*@code{eps}. |
|
3430
|
74 ## |
|
7001
|
75 ## @var{gmax}, @var{min}, @var{tol}, and @var{tol} must all be positive scalars. |
|
3430
|
76 ## @end table |
|
|
77 ## @strong{Outputs} |
|
|
78 ## @table @var |
|
3502
|
79 ## @item k |
|
5016
|
80 ## System controller. |
|
3430
|
81 ## @item g |
|
5016
|
82 ## Designed gain value. |
|
3502
|
83 ## @item gw |
|
5016
|
84 ## Closed loop system. |
|
3502
|
85 ## @item xinf |
|
5016
|
86 ## @acronym{ARE} solution matrix for regulator subproblem. |
|
3502
|
87 ## @item yinf |
|
5016
|
88 ## @acronym{ARE} solution matrix for filter subproblem. |
|
3430
|
89 ## @end table |
|
|
90 ## |
|
5016
|
91 ## References: |
|
3430
|
92 ## @enumerate |
|
5016
|
93 ## @item Doyle, Glover, Khargonekar, Francis, @cite{State-Space Solutions |
|
|
94 ## to Standard} |
|
|
95 ## @iftex |
|
|
96 ## @tex |
|
|
97 ## $ { \cal H }_2 $ @cite{and} $ { \cal H }_\infty $ |
|
|
98 ## @end tex |
|
|
99 ## @end iftex |
|
|
100 ## @ifinfo |
|
|
101 ## @cite{H-2 and H-infinity} |
|
|
102 ## @end ifinfo |
|
|
103 ## @cite{Control Problems}, @acronym{IEEE} @acronym{TAC} August 1989. |
|
3430
|
104 ## |
|
5016
|
105 ## @item Maciejowksi, J.M., @cite{Multivariable feedback design}, |
|
|
106 ## Addison-Wesley, 1989, @acronym{ISBN} 0-201-18243-2. |
|
3430
|
107 ## |
|
5016
|
108 ## @item Keith Glover and John C. Doyle, @cite{State-space formulae for all |
|
|
109 ## stabilizing controllers that satisfy an} |
|
|
110 ## @iftex |
|
|
111 ## @tex |
|
|
112 ## $ { \cal H }_\infty $@cite{norm} |
|
|
113 ## @end tex |
|
|
114 ## @end iftex |
|
|
115 ## @ifinfo |
|
|
116 ## @cite{H-infinity-norm} |
|
|
117 ## @end ifinfo |
|
|
118 ## @cite{bound and relations to risk sensitivity}, |
|
|
119 ## Systems & Control Letters 11, Oct. 1988, pp 167--172. |
|
3430
|
120 ## @end enumerate |
|
|
121 ## @end deftypefn |
|
|
122 |
|
|
123 ## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> |
|
|
124 ## Created: August 1995 |
|
|
125 ## Updated for Packed system structures December 1996 by John Ingram |
|
|
126 ## |
|
|
127 ## Revised by Kai P. Mueller April 1998 to solve the general H_infinity |
|
|
128 ## problem using unitary transformations Q (on w and z) |
|
|
129 ## and non-singular transformations R (on u and y). |
|
|
130 |
|
|
131 function [K, g, GW, Xinf, Yinf] = hinfsyn (Asys, nu, ny, gmin, gmax, gtol, ptol, tol) |
|
|
132 |
|
|
133 if( (nargin < 1) | (nargin > 8) ) |
|
6046
|
134 print_usage (); |
|
3430
|
135 endif |
|
|
136 ## set default arguments |
|
|
137 if(nargin < 8) |
|
|
138 tol = 200*eps; |
|
|
139 elseif(!is_sample(tol)) |
|
|
140 error("tol must be a positive scalar.") |
|
|
141 endif |
|
|
142 if(nargin < 7) |
|
|
143 ptol = 1e-9; |
|
|
144 elseif(!is_sample(ptol)) |
|
|
145 error("hinfsyn: ptol must be a positive scalar"); |
|
|
146 endif |
|
|
147 |
|
|
148 if(!is_sample(gmax) | !is_sample(gmin) | !is_sample(gtol) ) |
|
|
149 error(["hinfsyn: gmax=",num2str(gmax),", gmin=",num2str(gmin), ... |
|
|
150 "gtol=",num2str(gtol), " must be positive scalars."]) |
|
|
151 endif |
|
|
152 |
|
|
153 [chkdgkf,dgs] = is_dgkf(Asys,nu,ny,tol); |
|
|
154 |
|
|
155 if (! chkdgkf ) |
|
|
156 disp("hinfsyn: system does not meet required assumptions") |
|
|
157 help is_dgkf |
|
|
158 error("hinfsyn: exit"); |
|
|
159 endif |
|
|
160 |
|
|
161 ## extract dgs information |
|
|
162 nw = dgs.nw; nu = dgs.nu; |
|
|
163 A = dgs.A; B1 = dgs.Bw; B2 = dgs.Bu; |
|
|
164 C1 = dgs.Cz; D11 = dgs.Dzw; D12 = dgs.Dzu; nz = dgs.nz; |
|
|
165 C2 = dgs.Cy; D21 = dgs.Dyw; D22 = dgs.Dyu; ny = dgs.ny; |
|
|
166 d22nz = dgs.Dyu_nz; |
|
|
167 dflg = dgs.dflg; |
|
|
168 |
|
|
169 ## recover i/o transformations |
|
|
170 R12 = dgs.Ru; R21 = dgs.Ry; |
|
|
171 [ncstates, ndstates, nin, nout] = sysdimensions(Asys); |
|
|
172 Atsam = sysgettsam(Asys); |
|
|
173 [Ast, Ain, Aout] = sysgetsignals(Asys); |
|
|
174 |
|
|
175 BB = [B1, B2]; |
|
|
176 CC = [C1 ; C2]; |
|
|
177 DD = [D11, D12 ; D21, D22]; |
|
|
178 |
|
|
179 if (dflg == 0) |
|
|
180 n = ncstates; |
|
|
181 ## perform binary search to find gamma min |
|
|
182 ghi = gmax; |
|
|
183 ## derive a lower lower bound for gamma from D matrices |
|
|
184 xx1 = norm((eye(nz) - (D12/(D12'*D12))*D12')*D11); |
|
|
185 xx2 = norm(D11*(eye(nw)-(D21'/(D21*D21'))*D21)); |
|
|
186 glo = max(xx1, xx2); |
|
|
187 if (glo > gmin) |
|
|
188 disp(" *** D matrices indicate a greater value of gamma min."); |
|
|
189 fprintf(" gamma min (%f) superseeded by %f.", gmin, glo); |
|
|
190 glo = xx1; |
|
|
191 else |
|
|
192 glo = gmin; |
|
|
193 endif |
|
|
194 if (glo > ghi) |
|
|
195 fprintf(" *** lower bound of gamma greater than upper bound(%f)", ... |
|
|
196 glo, ghi); |
|
|
197 disp(" *** unable to continue, Goodbye."); |
|
|
198 return; |
|
|
199 endif |
|
|
200 |
|
|
201 de = ghi - glo; |
|
|
202 g = glo; |
|
|
203 search_state = 0; |
|
|
204 iteration_finished = 0; |
|
|
205 disp(" o structural tests passed, start of iteration..."); |
|
|
206 disp(" o <-> test passed # <-> test failed - <-> cannot test"); |
|
|
207 printf("----------------------------------------"); |
|
|
208 printf("--------------------------------------\n"); |
|
|
209 |
|
|
210 ## ------123456789012345678901234567890123456789012345678901234567890 |
|
|
211 printf(" .........X......... .........Y......... "); |
|
|
212 printf(".Z. PASS REMARKS\n"); |
|
|
213 printf(" ga iax nev ene sym pos iax nev ene sym pos "); |
|
|
214 printf("rho y/n ======>\n"); |
|
|
215 printf("----------------------------------------"); |
|
|
216 printf("--------------------------------------\n"); |
|
|
217 |
|
|
218 ## set up error messages |
|
4771
|
219 errmesg = {" o o o o o ", ... |
|
3430
|
220 " # - - - - ", ... |
|
|
221 " o # - - - ", ... |
|
|
222 " o o # - - ", ... |
|
|
223 " o o o # - ", ... |
|
|
224 " o o o o # ", ... |
|
4771
|
225 " - - - - - "}; |
|
|
226 errdesx = {"", ... |
|
3430
|
227 "X im eig.", ... |
|
|
228 "Hx not Ham.", ... |
|
|
229 "X inf.eig", ... |
|
|
230 "X not symm.", ... |
|
|
231 "X not pos", ... |
|
4771
|
232 "R singular"}; |
|
3430
|
233 |
|
4771
|
234 errdesy = {" ", ... |
|
3430
|
235 "Y im eig.", ... |
|
|
236 "Hy not Ham.", ... |
|
|
237 "Y inf.eig", ... |
|
|
238 "Y not symm.", ... |
|
|
239 "Y not pos", ... |
|
4771
|
240 "Rtilde singular"}; |
|
3430
|
241 |
|
|
242 |
|
|
243 ## now do the search |
|
|
244 while (!iteration_finished) |
|
|
245 switch (search_state) |
|
|
246 case (0) g = ghi; |
|
|
247 case (1) g = glo; |
|
|
248 case (2) g = 0.5 * (ghi + glo); |
|
|
249 otherwise error(" *** This should never happen!"); |
|
|
250 endswitch |
|
|
251 printf("%10.4f ", g); |
|
|
252 |
|
|
253 ## computing R and R~ |
|
|
254 d1dot = [D11, D12]; |
|
|
255 R = zeros(nin, nin); |
|
|
256 R(1:nw,1:nw) = -g*g*eye(nw); |
|
|
257 R = R + d1dot' * d1dot; |
|
|
258 ddot1 = [D11; D21]; |
|
|
259 Rtilde = zeros(nout, nout); |
|
|
260 Rtilde(1:nz,1:nz) = -g*g*eye(nz); |
|
|
261 Rtilde = Rtilde + ddot1 * ddot1'; |
|
|
262 |
|
|
263 [Xinf,x_ha_err] = hinfsyn_ric(A,BB,C1,d1dot,R,ptol); |
|
|
264 [Yinf,y_ha_err] = hinfsyn_ric(A',CC',B1',ddot1',Rtilde,ptol); |
|
|
265 |
|
|
266 ## assume failure for this gamma |
|
|
267 passed = 0; |
|
|
268 rerr=""; |
|
|
269 if (!x_ha_err && !y_ha_err) |
|
|
270 ## test spectral radius condition |
|
|
271 rho = max(abs(eig(Xinf * Yinf))); |
|
|
272 if (rho < g*g) |
|
|
273 ## spectral radius condition passed |
|
|
274 passed = 1; |
|
|
275 else |
|
|
276 rerr = sprintf("rho=%f",rho); |
|
|
277 endif |
|
|
278 endif |
|
|
279 |
|
|
280 if(x_ha_err >= 0 & x_ha_err <= 6) |
|
4771
|
281 printf("%s",errmesg{x_ha_err+1}); |
|
|
282 xerr = errdesx{x_ha_err+1}; |
|
3430
|
283 else |
|
|
284 error(" *** Xinf fail: this should never happen!"); |
|
|
285 endif |
|
|
286 |
|
|
287 if(y_ha_err >= 0 & y_ha_err <= 6) |
|
4771
|
288 printf("%s",errmesg{y_ha_err+1}); |
|
|
289 yerr = errdesy{y_ha_err+1}; |
|
3430
|
290 else |
|
|
291 error(" *** Yinf fail: this should never happen!"); |
|
|
292 endif |
|
|
293 |
|
|
294 if(passed) printf(" y all tests passed.\n"); |
|
|
295 else printf(" n %s/%s%s\n",xerr,yerr,rerr); endif |
|
|
296 |
|
|
297 if (passed && (de/g < gtol)) |
|
|
298 search_state = 3; |
|
|
299 endif |
|
|
300 |
|
|
301 switch (search_state) |
|
|
302 case (0) |
|
|
303 if (!passed) |
|
|
304 ## upper bound must pass but did not |
|
|
305 fprintf(" *** the upper bound of gamma (%f) is too small.\n", g); |
|
|
306 iteration_finished = 2; |
|
|
307 else |
|
|
308 search_state = 1; |
|
|
309 endif |
|
|
310 case (1) |
|
|
311 if (!passed) search_state = 2; |
|
|
312 else |
|
|
313 ## lower bound must not pass but passed |
|
|
314 fprintf(" *** the lower bound of gamma (%f) passed.\n", g); |
|
|
315 iteration_finished = 3; |
|
|
316 endif |
|
|
317 case (2) |
|
|
318 ## Normal case; must check that singular R, Rtilde wasn't the problem. |
|
|
319 if ((!passed) & (x_ha_err != 6) & (y_ha_err != 6) ) glo = g; |
|
|
320 else ghi = g; endif |
|
|
321 de = ghi - glo; |
|
|
322 case (3) iteration_finished = 1; # done |
|
|
323 otherwise error(" *** This should never happen!"); |
|
|
324 endswitch |
|
|
325 endwhile |
|
|
326 |
|
|
327 printf("----------------------------------------"); |
|
|
328 printf("--------------------------------------\n"); |
|
|
329 if (iteration_finished != 1) |
|
|
330 K = []; |
|
|
331 else |
|
|
332 ## success: compute controller |
|
|
333 fprintf(" hinfsyn final: glo=%f ghi=%f, test gain g=%f\n", \ |
|
|
334 glo, ghi, g); |
|
|
335 printf("----------------------------------------"); |
|
|
336 printf("--------------------------------------\n"); |
|
|
337 Z = inv(eye(ncstates) - Yinf*Xinf/g/g); |
|
|
338 F = -R \ (d1dot'*C1 + BB'*Xinf); |
|
|
339 H = -(B1*ddot1' + Yinf*CC') / Rtilde; |
|
|
340 K = hinf_ctr(dgs,F,H,Z,g); |
|
|
341 |
|
|
342 Kst = strappend(Ast,"_K"); |
|
|
343 Kin = strappend(Aout((nout-ny+1):(nout)),"_K"); |
|
|
344 Kout = strappend(Ain((nin-nu+1):(nin)),"_K"); |
|
|
345 [Ac, Bc, Cc, Dc] = sys2ss(K); |
|
4771
|
346 K = ss(Ac,Bc,Cc,Dc,Atsam,ncstates,ndstates,Kst,Kin,Kout); |
|
3430
|
347 if (nargout >= 3) |
|
|
348 GW = starp(Asys, K); |
|
|
349 endif |
|
|
350 endif |
|
|
351 |
|
|
352 elseif(ndstates) |
|
|
353 |
|
|
354 ## discrete time solution |
|
|
355 error("hinfsyn: discrete-time case not yet implemented") |
|
|
356 |
|
|
357 endif |
|
|
358 |
|
|
359 endfunction |
