|
3430
|
1 ## Copyright (C) 1996, 1998 Kai P. Mueller |
|
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
|
6 ## under the terms of the GNU General Public License as published by the |
|
|
7 ## Free Software Foundation; either version 2, or (at your option) any |
|
|
8 ## later version. |
|
|
9 ## |
|
|
10 ## Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
13 ## for more details. |
|
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
|
16 ## along with Octave; see the file COPYING. If not, write to the Free |
|
|
17 ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. |
|
|
18 |
|
|
19 ## -*- texinfo -*- |
|
|
20 ## @deftypefn {Function File} {} hinfdemo () |
|
|
21 ## |
|
5016
|
22 ## @iftex |
|
|
23 ## @tex |
|
|
24 ## $ { \cal H }_\infty $ |
|
|
25 ## @end tex |
|
|
26 ## @end iftex |
|
|
27 ## @ifinfo |
|
|
28 ## H-infinity |
|
|
29 ## @end ifinfo |
|
|
30 ## design demos for continuous @acronym{SISO} and @acronym{MIMO} systems and a |
|
|
31 ## discrete system. The @acronym{SISO} system is difficult to control because |
|
|
32 ## it is non-minimum-phase and unstable. The second design example |
|
|
33 ## controls the @command{jet707} plant, the linearized state space model of a |
|
|
34 ## Boeing 707-321 aircraft at @var{v}=80 m/s |
|
|
35 ## @iftex |
|
|
36 ## @tex |
|
|
37 ## ($M = 0.26$, $G_{a0} = -3^{\circ}$, ${\alpha}_0 = 4^{\circ}$, ${\kappa}= 50^{\circ}$). |
|
|
38 ## @end tex |
|
|
39 ## @end iftex |
|
|
40 ## @ifinfo |
|
|
41 ## (@var{M} = 0.26, @var{Ga0} = -3 deg, @var{alpha0} = 4 deg, @var{kappa} = 50 deg). |
|
|
42 ## @end ifinfo |
|
|
43 ## Inputs: (1) thrust and (2) elevator angle |
|
|
44 ## Outputs: (1) airspeed and (2) pitch angle. The discrete system is a |
|
3430
|
45 ## stable and second order. |
|
|
46 ## |
|
|
47 ## @table @asis |
|
5016
|
48 ## @item @acronym{SISO} plant: |
|
|
49 ## |
|
|
50 ## @iftex |
|
|
51 ## @tex |
|
|
52 ## $$ G(s) = { s-2 \over (s+2) (s-1) } $$ |
|
|
53 ## @end tex |
|
|
54 ## @end iftex |
|
|
55 ## @ifinfo |
|
|
56 ## @example |
|
3430
|
57 ## @group |
|
|
58 ## s - 2 |
|
|
59 ## G(s) = -------------- |
|
|
60 ## (s + 2)(s - 1) |
|
5016
|
61 ## @end group |
|
|
62 ## @end example |
|
|
63 ## @end ifinfo |
|
|
64 ## |
|
|
65 ## @example |
|
|
66 ## @group |
|
3430
|
67 ## |
|
|
68 ## +----+ |
|
|
69 ## -------------------->| W1 |---> v1 |
|
|
70 ## z | +----+ |
|
5016
|
71 ## ----|-------------+ |
|
|
72 ## | | |
|
3430
|
73 ## | +---+ v y +----+ |
|
|
74 ## u *--->| G |--->O--*-->| W2 |---> v2 |
|
|
75 ## | +---+ | +----+ |
|
|
76 ## | | |
|
|
77 ## | +---+ | |
|
|
78 ## -----| K |<------- |
|
|
79 ## +---+ |
|
|
80 ## @end group |
|
5016
|
81 ## @end example |
|
|
82 ## |
|
|
83 ## @iftex |
|
|
84 ## @tex |
|
|
85 ## $$ { \rm min } \Vert T_{vz} \Vert _\infty $$ |
|
|
86 ## @end tex |
|
|
87 ## @end iftex |
|
|
88 ## @ifinfo |
|
|
89 ## @example |
|
|
90 ## min || T || |
|
|
91 ## vz infty |
|
|
92 ## @end example |
|
|
93 ## @end ifinfo |
|
3430
|
94 ## |
|
5016
|
95 ## @var{W1} und @var{W2} are the robustness and performance weighting |
|
|
96 ## functions. |
|
|
97 ## |
|
|
98 ## @item @acronym{MIMO} plant: |
|
|
99 ## The optimal controller minimizes the |
|
|
100 ## @iftex |
|
|
101 ## @tex |
|
|
102 ## $ { \cal H }_\infty $ |
|
|
103 ## @end tex |
|
|
104 ## @end iftex |
|
|
105 ## @ifinfo |
|
|
106 ## H-infinity |
|
|
107 ## @end ifinfo |
|
|
108 ## norm of the |
|
|
109 ## augmented plant @var{P} (mixed-sensitivity problem): |
|
|
110 ## @example |
|
3430
|
111 ## @group |
|
|
112 ## w |
|
|
113 ## 1 -----------+ |
|
|
114 ## | +----+ |
|
|
115 ## +---------------------->| W1 |----> z1 |
|
|
116 ## w | | +----+ |
|
|
117 ## 2 ------------------------+ |
|
|
118 ## | | | |
|
|
119 ## | v +----+ v +----+ |
|
|
120 ## +--*-->o-->| G |-->o--*-->| W2 |---> z2 |
|
|
121 ## | +----+ | +----+ |
|
|
122 ## | | |
|
|
123 ## ^ v |
|
5016
|
124 ## u y (to K) |
|
|
125 ## (from controller K) |
|
|
126 ## @end group |
|
|
127 ## @end example |
|
3430
|
128 ## |
|
5016
|
129 ## @iftex |
|
|
130 ## @tex |
|
|
131 ## $$ \left [ \matrix{ z_1 \cr |
|
|
132 ## z_2 \cr |
|
|
133 ## y } \right ] = |
|
|
134 ## P \left [ \matrix{ w_1 \cr |
|
|
135 ## w_2 \cr |
|
|
136 ## u } \right ] $$ |
|
|
137 ## @end tex |
|
|
138 ## @end iftex |
|
|
139 ## @ifinfo |
|
|
140 ## @example |
|
|
141 ## @group |
|
3430
|
142 ## + + + + |
|
|
143 ## | z | | w | |
|
|
144 ## | 1 | | 1 | |
|
|
145 ## | z | = [ P ] * | w | |
|
|
146 ## | 2 | | 2 | |
|
|
147 ## | y | | u | |
|
|
148 ## + + + + |
|
|
149 ## @end group |
|
5016
|
150 ## @end example |
|
|
151 ## @end ifinfo |
|
3430
|
152 ## |
|
5016
|
153 ## @item Discrete system: |
|
3430
|
154 ## This is not a true discrete design. The design is carried out |
|
|
155 ## in continuous time while the effect of sampling is described by |
|
|
156 ## a bilinear transformation of the sampled system. |
|
5016
|
157 ## This method works quite well if the sampling period is ``small'' |
|
3430
|
158 ## compared to the plant time constants. |
|
|
159 ## |
|
5016
|
160 ## @item The continuous plant: |
|
|
161 ## @iftex |
|
|
162 ## @tex |
|
|
163 ## $$ G(s) = { 1 \over (s+2)(s+1) } $$ |
|
|
164 ## @end tex |
|
|
165 ## @end iftex |
|
|
166 ## |
|
|
167 ## @ifinfo |
|
|
168 ## @example |
|
3430
|
169 ## @group |
|
|
170 ## 1 |
|
|
171 ## G (s) = -------------- |
|
|
172 ## k (s + 2)(s + 1) |
|
|
173 ## |
|
|
174 ## @end group |
|
5016
|
175 ## @end example |
|
|
176 ## @end ifinfo |
|
|
177 ## |
|
|
178 ## is discretised with a @acronym{ZOH} (Sampling period = @var{Ts} = 1 second): |
|
|
179 ## @iftex |
|
|
180 ## @tex |
|
|
181 ## $$ G(z) = { 0.199788z + 0.073498 \over (z - 0.36788) (z - 0.13534) } $$ |
|
|
182 ## @end tex |
|
|
183 ## @end iftex |
|
|
184 ## @ifinfo |
|
|
185 ## @example |
|
3430
|
186 ## @group |
|
|
187 ## |
|
|
188 ## 0.199788z + 0.073498 |
|
5016
|
189 ## G(z) = -------------------------- |
|
3430
|
190 ## (z - 0.36788)(z - 0.13534) |
|
5016
|
191 ## @end group |
|
|
192 ## @end example |
|
|
193 ## @end ifinfo |
|
|
194 ## |
|
|
195 ## @example |
|
|
196 ## @group |
|
3430
|
197 ## |
|
|
198 ## +----+ |
|
|
199 ## -------------------->| W1 |---> v1 |
|
|
200 ## z | +----+ |
|
5016
|
201 ## ----|-------------+ |
|
|
202 ## | | |
|
3430
|
203 ## | +---+ v +----+ |
|
|
204 ## *--->| G |--->O--*-->| W2 |---> v2 |
|
|
205 ## | +---+ | +----+ |
|
|
206 ## | | |
|
|
207 ## | +---+ | |
|
|
208 ## -----| K |<------- |
|
|
209 ## +---+ |
|
|
210 ## @end group |
|
5016
|
211 ## @end example |
|
|
212 ## @iftex |
|
|
213 ## @tex |
|
|
214 ## $$ { \rm min } \Vert T_{vz} \Vert _\infty $$ |
|
|
215 ## @end tex |
|
|
216 ## @end iftex |
|
|
217 ## @ifinfo |
|
|
218 ## @example |
|
|
219 ## min || T || |
|
|
220 ## vz infty |
|
|
221 ## @end example |
|
|
222 ## @end ifinfo |
|
|
223 ## @var{W1} and @var{W2} are the robustness and performance weighting |
|
|
224 ## functions. |
|
3430
|
225 ## @end table |
|
|
226 ## @end deftypefn |
|
|
227 |
|
|
228 ## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> |
|
|
229 ## Created: April 30, 1998 |
|
|
230 |
|
|
231 yn = []; |
|
|
232 while (length(yn) < 1) |
|
|
233 yn = input(" * [s]iso, [m]imo, or [d]iscrete design? [no default]: ","S"); |
|
|
234 endwhile |
|
|
235 if ((yn(1) == "s") | (yn(1) == 'S')) |
|
|
236 sys_type = 1; |
|
|
237 elseif ((yn(1) == "m") | (yn(1) == 'M')) |
|
|
238 sys_type = 2; |
|
|
239 elseif ((yn(1) == "d") | (yn(1) == 'D')) |
|
|
240 sys_type = 3; |
|
|
241 else |
|
|
242 disp(" *** no system type specified, hinfdemo terminated."); |
|
|
243 return; |
|
|
244 endif |
|
|
245 |
|
|
246 echo off |
|
|
247 switch (sys_type) |
|
|
248 |
|
|
249 case (1) |
|
|
250 ## siso |
|
|
251 disp(" "); |
|
|
252 disp(" ----------------------------------------------"); |
|
|
253 disp(" H_infinity optimal control for the SISO plant:"); |
|
|
254 disp(" "); |
|
|
255 disp(" s - 2"); |
|
|
256 disp(" G(s) = --------------"); |
|
|
257 disp(" (s + 2)(s - 1)"); |
|
|
258 disp(" "); |
|
|
259 disp(" ----------------------------------------------"); |
|
|
260 disp(" "); |
|
|
261 |
|
|
262 ## weighting on actuator u |
|
|
263 W1 = wgt1o(0.05, 100.0, 425.0); |
|
|
264 ## weighting on controlled variable y |
|
|
265 W2 = wgt1o(10.0, 0.05, 0.001); |
|
|
266 ## plant |
|
|
267 G = tf2sys([1 -2],[1 1 -2]); |
|
|
268 |
|
|
269 ## need One as the pseudo transfer function One = 1 |
|
|
270 One = ugain(1); |
|
|
271 disp(" o forming P..."); |
|
|
272 psys = buildssic([1 4;2 4;3 1],[3],[2 3 5],[3 4],G,W1,W2,One); |
|
|
273 disp(" "); |
|
|
274 disp(" o controller design..."); |
|
|
275 [K, gfin, GW]=hinfsyn(psys, 1, 1, 0.1, 10.0, 0.02); |
|
|
276 disp(" "); |
|
|
277 disp("-- OK ----------------------------------------------"); |
|
|
278 |
|
|
279 disp(" Closed loop poles:"); |
|
|
280 damp(GW); |
|
|
281 ## disp(" o Testing H_infinity norm: (hinfnorm does not work)"); |
|
|
282 ## hinfnorm(GW); |
|
|
283 |
|
|
284 disp(" "); |
|
|
285 yn = input(" * Plot closed loop step response? [n]: ","S"); |
|
|
286 if (length(yn) >= 1) |
|
|
287 if ((yn(1) == "y") || (yn(1) == 'Y')) |
|
|
288 disp(" o step responses of T and KS..."); |
|
|
289 GW = buildssic([1 2; 2 1], [], [1 2], [-2], G, K); |
|
|
290 figure(1); |
|
|
291 step(GW, 1, 10); |
|
|
292 endif |
|
|
293 endif |
|
|
294 |
|
|
295 case (2) |
|
|
296 ## mimo |
|
|
297 disp(" "); |
|
|
298 disp(" -----------------------------------------------"); |
|
|
299 disp(" H_inf optimal control for the jet707 plant"); |
|
|
300 disp(" -----------------------------------------------"); |
|
|
301 disp(" "); |
|
|
302 |
|
|
303 ## Weighting function on u (robustness weight) |
|
|
304 ww1 = wgt1o(0.01,5,0.9); |
|
|
305 ww2 = wgt1o(0.01,5,2.2); |
|
|
306 W1 = buildssic([1 0;2 0],[],[1 2],[1 2],ww1,ww2); |
|
|
307 ## Weighting function on y (performance weight) |
|
|
308 ww1 = wgt1o(250,0.1,0.0001); |
|
|
309 ww2 = wgt1o(250,0.1,0.0002); |
|
|
310 W2 = buildssic([1 0;2 0],[],[1 2],[1 2],ww1,ww2); |
|
|
311 ## plant (2 x 2 system) |
|
|
312 G = jet707; |
|
|
313 |
|
|
314 disp(" o forming P..."); |
|
|
315 One = ugain(2); |
|
|
316 Clst = [1 7; 2 8; 3 7; 4 8; 5 1; 6 2]; |
|
|
317 P = buildssic(Clst,[5 6],[3:6 9 10],[1 2 5:8],G,W1,W2,One); |
|
|
318 |
|
|
319 disp(" "); |
|
|
320 disp(" o controller design..."); |
|
|
321 K = hinfsyn(P, 2, 2, 0.25, 10.0, 0.005); |
|
|
322 |
|
|
323 disp(" "); |
|
|
324 yn = input(" * Plot closed loop step responses? [n]: ","S"); |
|
|
325 if (length(yn) >= 1) |
|
|
326 if ((yn(1) == "y") || (yn(1) == 'Y')) |
|
|
327 disp(" o step responses of T and KS..."); |
|
|
328 GW = buildssic([1 3;2 4;3 1;4 2],[],[1 2 3 4],[-3 -4],G,K); |
|
|
329 |
|
|
330 disp(" "); |
|
|
331 disp(" FIGURE 1: speed refence => 1, pitch angle ref. => 0"); |
|
|
332 disp(" ==================================================="); |
|
|
333 disp(" y1: speed (should be 1)"); |
|
|
334 disp(" y2: pitch angle (should remain 0)"); |
|
|
335 disp(" y3: thrust (should be a slow transient)"); |
|
|
336 disp(" y6: elevator (should be a faster transient)"); |
|
|
337 disp(" "); |
|
|
338 disp(" FIGURE 2: speed refence => 0, pitch angle ref. => 1"); |
|
|
339 disp(" ==================================================="); |
|
|
340 disp(" y1: speed (should remain 0)"); |
|
|
341 disp(" y2: pitch angle (should be 1)"); |
|
|
342 disp(" y3: thrust (should be a slow transient)"); |
|
|
343 disp(" y6: elevator (should be a faster transient)"); |
|
|
344 disp(" "); |
|
|
345 figure(1) |
|
|
346 step(GW); |
|
|
347 figure(2) |
|
|
348 step(GW,2); |
|
|
349 endif |
|
|
350 endif |
|
|
351 |
|
|
352 case (3) |
|
|
353 ## discrete |
|
|
354 disp(" "); |
|
|
355 disp(" --------------------------------------------------"); |
|
|
356 disp(" Discrete H_infinity optimal control for the plant:"); |
|
|
357 disp(" "); |
|
|
358 disp(" 0.199788z + 0.073498"); |
|
|
359 disp(" G(s) = --------------------------"); |
|
|
360 disp(" (z - 0.36788)(z - 0.13533)"); |
|
|
361 disp(" --------------------------------------------------"); |
|
|
362 disp(" "); |
|
|
363 |
|
|
364 ## sampling time |
|
|
365 Ts = 1.0; |
|
|
366 ## weighting on actuator value u |
|
|
367 W1 = wgt1o(0.1, 200.0, 50.0); |
|
|
368 ## weighting on controlled variable y |
|
|
369 W2 = wgt1o(350.0, 0.05, 0.0002); |
|
|
370 ## omega axis |
|
|
371 ww = logspace(-4.99, 3.99, 100); |
|
|
372 if (columns(ww) > 1); ww = ww'; endif |
|
|
373 |
|
|
374 ## continuous plant |
|
|
375 G = tf2sys(2,[1 3 2]); |
|
|
376 ## discrete plant with zoh |
|
|
377 Gd = c2d(G, Ts); |
|
|
378 ## w-plane (continuous representation of the sampled system) |
|
|
379 Gw = d2c(Gd, "bi"); |
|
|
380 |
|
|
381 disp(" "); |
|
|
382 disp(" o building P..."); |
|
|
383 ## need One as the pseudo transfer function One = 1 |
|
|
384 One = ugain(1); |
|
|
385 psys = buildssic([1 4;2 4;3 1],[3],[2 3 5],[3 4],Gw,W1,W2,One); |
|
|
386 disp(" o controller design..."); |
|
|
387 [K, gfin, GWC] = hinfsyn(psys, 1, 1, 0.1, 10.0, 0.02); |
|
|
388 |
|
|
389 disp(" "); |
|
|
390 fig_n = 1; |
|
|
391 yn = input(" * Plot magnitudes of W1KS and W2S? [n]: ","S"); |
|
|
392 if (length(yn) >= 1) |
|
|
393 if ((yn(1) == "y") || (yn(1) == 'Y')) |
|
|
394 disp(" o magnitudes of W1KS and W2S..."); |
|
|
395 gwx = sysprune(GWC, 1, 1); |
|
|
396 mag1 = bode(gwx, ww); |
|
|
397 if (columns(mag1) > 1); mag1 = mag1'; endif |
|
|
398 gwx = sysprune(GWC, 2, 1); |
|
|
399 mag2 = bode(gwx, ww); |
|
|
400 if (columns(mag2) > 1); mag2 = mag2'; endif |
|
|
401 figure(fig_n) |
|
|
402 fig_n = fig_n + 1; |
|
5215
|
403 __gnuplot_set__ grid |
|
3430
|
404 loglog(ww, [mag1 mag2]); |
|
|
405 endif |
|
|
406 endif |
|
|
407 |
|
|
408 Kd = c2d(K, "bi", Ts); |
|
|
409 GG = buildssic([1 2; 2 1], [], [1 2], [-2], Gd, Kd); |
|
|
410 disp(" o closed loop poles..."); |
|
|
411 damp(GG); |
|
|
412 |
|
|
413 disp(" "); |
|
|
414 yn = input(" * Plot closed loop step responses? [n]: ","S"); |
|
|
415 if (length(yn) >= 1) |
|
|
416 if ((yn(1) == "y") || (yn(1) == 'Y')) |
|
|
417 disp(" o step responses of T and KS..."); |
|
|
418 figure(fig_n) |
|
|
419 step(GG, 1, 10); |
|
|
420 endif |
|
|
421 endif |
|
|
422 |
|
|
423 endswitch |
|
|
424 |
|
|
425 disp(" o hinfdemo terminated successfully."); |
|
|
426 |
|
|
427 ## KPM-hinfdemo/End |
