|
3440
|
1 ## Copyright (C) 1996 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. |
|
3440
|
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. |
|
3440
|
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/>. |
|
3440
|
18 |
|
|
19 ## -*- texinfo -*- |
|
|
20 ## @deftypefn {Function File} {} dgkfdemo () |
|
5016
|
21 ## Octave Controls toolbox demo: |
|
|
22 ## @iftex |
|
|
23 ## @tex |
|
|
24 ## $ { \cal H }_2 $/$ { \cal H }_\infty $ |
|
|
25 ## @end tex |
|
|
26 ## @end iftex |
|
|
27 ## @ifinfo |
|
|
28 ## H-2/H-infinity |
|
|
29 ## @end ifinfo |
|
|
30 ## options demos. |
|
3440
|
31 ## @end deftypefn |
|
|
32 |
|
|
33 ## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> |
|
|
34 ## Created: June 1995 |
|
|
35 |
|
|
36 function dgkfdemo () |
|
|
37 |
|
6522
|
38 save_val = page_screen_output (); |
|
|
39 page_screen_output (0); |
|
3440
|
40 while (1) |
|
|
41 clc |
|
|
42 sel = 0; |
|
|
43 while (sel > 10 || sel < 1) |
|
|
44 sel = menu ("Octave H2/Hinfinity options demo", |
|
|
45 "LQ regulator", |
|
|
46 "LG state estimator", |
|
|
47 "LQG optimal control design", |
|
|
48 "H2 gain of a system", |
|
|
49 "H2 optimal controller of a system", |
|
|
50 "Hinf gain of a system", |
|
|
51 "Hinf optimal controller of a SISO system", |
|
|
52 "Hinf optimal controller of a MIMO system", |
|
|
53 "Discrete-time Hinf optimal control by bilinear transform", |
|
|
54 "Return to main demo menu"); |
|
|
55 endwhile |
|
|
56 if (sel == 1) |
|
|
57 disp("Linear/Quadratic regulator design:") |
|
|
58 disp("Compute optimal state feedback via the lqr command...") |
|
|
59 help lqr |
|
|
60 disp(" ") |
|
|
61 disp("Example:") |
|
|
62 A = [0, 1; -2, -1] |
|
|
63 B = [0; 1] |
|
|
64 Q = [1, 0; 0, 0] |
|
|
65 R = 1 |
|
|
66 disp("Q = state penalty matrix; R = input penalty matrix") |
|
|
67 prompt |
|
|
68 disp("Compute state feedback gain k, ARE solution P, and closed-loop") |
|
|
69 disp("poles as follows:"); |
|
|
70 cmd = "[k, p, e] = lqr(A,B,Q,R)"; |
|
|
71 run_cmd |
|
|
72 prompt |
|
|
73 disp("A similar approach can be used for LTI discrete-time systems") |
|
|
74 disp("by using the dlqr command in place of lqr (see LQG example).") |
|
|
75 elseif (sel == 2) |
|
|
76 disp("Linear/Gaussian estimator design:") |
|
|
77 disp("Compute optimal state estimator via the lqe command...") |
|
|
78 help lqe |
|
|
79 disp(" ") |
|
|
80 disp("Example:") |
|
|
81 A = [0, 1; -2, -1] |
|
|
82 disp("disturbance entry matrix G") |
|
|
83 G = eye(2) |
|
|
84 disp("Output measurement matrix C") |
|
|
85 C = [0, 1] |
|
|
86 SigW = [1, 0; 0, 1] |
|
|
87 SigV = 1 |
|
|
88 disp("SigW = input disturbance intensity matrix;") |
|
|
89 disp("SigV = measurement noise intensity matrix") |
|
|
90 prompt |
|
|
91 disp("Compute estimator feedback gain k, ARE solution P, and estimator") |
|
|
92 disp("poles via the command: ") |
|
|
93 cmd = "[k, p, e] = lqe(A,G,C,SigW,SigV)"; |
|
|
94 run_cmd |
|
|
95 disp("A similar approach can be used for LTI discrete-time systems") |
|
|
96 disp("by using the dlqe command in place of lqe (see LQG example).") |
|
|
97 elseif (sel == 3) |
|
|
98 disp("LQG optimal controller of a system:") |
|
|
99 disp("Input accepted as either A,B,C matrices or in system data structure form") |
|
|
100 disp("in both discrete and continuous time.") |
|
|
101 disp("Example 1: continuous time design:") |
|
|
102 prompt |
|
|
103 help lqg |
|
|
104 disp("Example system") |
|
|
105 A = [0, 1; .5, .5]; |
|
|
106 B = [0; 2]; |
|
|
107 G = eye(2) |
|
|
108 C = [1, 1]; |
|
4771
|
109 sys = ss(A, [B, G], C); |
|
3440
|
110 sys = syssetsignals(sys,"in", ... |
|
|
111 ["control input"; "disturbance 1"; "disturbance 2"]); |
|
|
112 sysout(sys) |
|
|
113 prompt |
|
|
114 disp("Filtering/estimator parameters:") |
|
|
115 SigW = eye(2) |
|
|
116 SigV = 1 |
|
|
117 prompt |
|
|
118 disp("State space (LQR) parameters Q and R are:") |
|
|
119 Q = eye(2) |
|
|
120 R = 1 |
|
|
121 cmd = "[K,Q1,P1,Ee,Er] = lqg(sys,SigW,SigV,Q,R,1);"; |
|
|
122 run_cmd |
|
|
123 disp("Check: closed loop system A-matrix is") |
|
|
124 disp(" [A, B*Cc]") |
|
|
125 disp(" [Bc*C, Ac ]") |
|
|
126 cmd = "[Ac, Bc, Cc] = sys2ss(K);"; |
|
|
127 run_cmd |
|
|
128 cmd = "Acl = [A, B*Cc; Bc*C, Ac]"; |
|
|
129 run_cmd |
|
|
130 disp("Check: poles of Acl:") |
|
|
131 Acl_poles = sortcom(eig(Acl)) |
|
|
132 disp("Predicted poles from design = union(Er,Ee)") |
|
|
133 cmd = "pred_poles = sortcom([Er; Ee])"; |
|
|
134 run_cmd |
|
|
135 disp("Example 2: discrete-time example") |
|
4771
|
136 cmd1 = "Dsys = ss(A, [G, B], C, [0, 0, 0], 1);"; |
|
3440
|
137 cmd2 = "[K,Q1,P1,Ee,Er] = lqg(Dsys,SigW, SigV,Q,R);"; |
|
|
138 disp("Run commands:") |
|
|
139 cmd = cmd1; |
|
|
140 run_cmd |
|
|
141 cmd = cmd2; |
|
|
142 run_cmd |
|
|
143 prompt |
|
|
144 disp("Check: closed loop system A-matrix is") |
|
|
145 disp(" [A, B*Cc]") |
|
|
146 disp(" [Bc*C, Ac ]") |
|
|
147 [Ac,Bc,Cc] = sys2ss(K); |
|
|
148 Acl = [A, B*Cc; Bc*C, Ac] |
|
|
149 prompt |
|
|
150 disp("Check: poles of Acl:") |
|
|
151 Acl_poles = sortcom(eig(Acl)) |
|
|
152 disp("Predicted poles from design = union(Er,Ee)") |
|
|
153 pred_poles = sortcom([Er;Ee]) |
|
|
154 elseif (sel == 4) |
|
|
155 disp("H2 gain of a system: (Energy in impulse response)") |
|
|
156 disp("Example 1: Stable plant:") |
|
|
157 cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)"; |
|
|
158 run_cmd |
|
|
159 disp("Put into Packed system form:") |
|
4771
|
160 cmd = "Asys = ss(A,B,C);"; |
|
3440
|
161 run_cmd |
|
|
162 disp("Evaluate system 2-norm (impulse response energy):"); |
|
|
163 cmd = "AsysH2 = h2norm(Asys)"; |
|
|
164 run_cmd |
|
|
165 disp("Compare with a plot of the system impulse response:") |
|
|
166 tt = 0:0.1:20; |
|
|
167 for ii=1:length(tt) |
|
|
168 ht(ii) = C*expm(A*tt(ii))*B; |
|
|
169 endfor |
|
|
170 plot(tt,ht) |
|
|
171 title("impulse response of example plant") |
|
|
172 prompt |
|
|
173 disp("Example 2: unstable plant") |
|
|
174 cmd = "A = [0, 1; 2, 1]"; |
|
|
175 eval(cmd); |
|
|
176 cmd = "B = [0; 1]"; |
|
|
177 eval(cmd); |
|
|
178 cmd = "C = [1, 0]"; |
|
|
179 eval(cmd); |
|
|
180 cmd = "sys_poles = eig(A)"; |
|
|
181 run_cmd |
|
|
182 prompt |
|
|
183 disp("Put into system data structure form:") |
|
4771
|
184 cmd="Bsys = ss(A,B,C);"; |
|
3440
|
185 run_cmd |
|
|
186 disp("Evaluate 2-norm:") |
|
|
187 cmd = "BsysH2 = h2norm(Bsys)"; |
|
|
188 run_cmd |
|
|
189 disp(" ") |
|
|
190 prompt("NOTICE: program returns a value without an error signal.") |
|
|
191 disp("") |
|
|
192 |
|
|
193 elseif (sel == 5) |
|
|
194 disp("H2 optimal controller of a system: command = h2syn:") |
|
|
195 prompt |
|
|
196 help h2syn |
|
|
197 prompt |
|
|
198 disp("Example system: double integrator with output noise and") |
|
|
199 disp("input disturbance:") |
|
|
200 disp(" "); |
|
|
201 disp(" -------------------->y2"); |
|
|
202 disp(" | _________"); |
|
|
203 disp("u(t)-->o-->| 1/s^2 |-->o-> y1"); |
|
|
204 disp(" ^ --------- ^"); |
|
|
205 disp(" | |"); |
|
|
206 disp(" w1(t) w2(t)"); |
|
|
207 disp(" ") |
|
|
208 disp("w enters the system through B1, u through B2") |
|
|
209 disp("z = [y1; y2] is obtained through C1, y=y1 through C2"); |
|
|
210 disp(" ") |
|
|
211 cmd = "A = [0, 1; 0, 0]; B1 = [0, 0; 1, 0]; B2 = [0; 1];"; |
|
|
212 disp(cmd) |
|
|
213 eval(cmd); |
|
|
214 cmd = "C1 = [1, 0; 0, 0]; C2 = [1, 0]; D11 = zeros(2);"; |
|
|
215 disp(cmd) |
|
|
216 eval(cmd); |
|
|
217 cmd = "D12 = [0; 1]; D21 = [0, 1]; D22 = 0; D = [D11, D12; D21, D22];"; |
|
|
218 disp(cmd) |
|
|
219 eval(cmd); |
|
|
220 disp("Design objective: compute U(s)=K(s)Y1(s) to minimize the closed") |
|
|
221 disp("loop impulse response from w(t) =[w1; w2] to z(t) = [y1; y2]"); |
|
|
222 prompt |
|
|
223 disp("First: pack system:") |
|
4771
|
224 cmd="Asys = ss(A, [B1, B2], [C1; C2], D);"; |
|
3440
|
225 run_cmd |
|
|
226 disp("Open loop multivariable Bode plot: (will take a moment)") |
|
|
227 cmd="bode(Asys);"; |
|
|
228 run_cmd |
|
|
229 prompt("Press a key to close plot and continue"); |
|
|
230 closeplot |
|
|
231 disp("Controller design command: (only need 1st two output arguments)") |
|
|
232 cmd="[K,gain, Kc, Kf, Pc, Pf] = h2syn(Asys,1,1);"; |
|
|
233 run_cmd |
|
|
234 disp("Controller is:") |
|
|
235 cmd = "sysout(K)"; |
|
|
236 run_cmd |
|
|
237 disp(["returned gain value is: ",num2str(gain)]); |
|
|
238 disp("Check: close the loop and then compute h2norm:") |
|
|
239 prompt |
|
|
240 cmd="K_loop = sysgroup(Asys,K);"; |
|
|
241 run_cmd |
|
|
242 cmd = "Kcl = sysconnect(K_loop,[3,4],[4,3]);"; |
|
|
243 run_cmd |
|
|
244 cmd = "Kcl = sysprune(Kcl,[1,2],[1,2]);"; |
|
|
245 run_cmd |
|
|
246 cmd="gain_Kcl = h2norm(Kcl)"; |
|
|
247 run_cmd |
|
|
248 cmd="gain_err = gain_Kcl - gain"; |
|
|
249 run_cmd |
|
|
250 disp("Check: multivarible bode plot:") |
|
|
251 cmd="bode(Kcl);"; |
|
|
252 run_cmd |
|
|
253 prompt |
|
|
254 disp("Related functions: is_dgkf, is_controllable, is_stabilizable,") |
|
|
255 disp(" is_observable, is_detectable") |
|
|
256 elseif (sel == 6) |
|
|
257 disp("Hinfinity gain of a system: (max gain over all j-omega)") |
|
|
258 disp("Example 1: Stable plant:") |
|
|
259 cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)"; |
|
|
260 run_cmd |
|
|
261 disp("Pack into system format:") |
|
4771
|
262 cmd = "Asys = ss(A,B,C);"; |
|
3440
|
263 run_cmd |
|
|
264 disp("The infinity norm must be computed iteratively by") |
|
|
265 disp("binary search. For this example, we select tolerance tol = 0.01, ") |
|
|
266 disp("min gain gmin = 1e-2, max gain gmax=1e4.") |
|
|
267 disp("Search quits when upper bound <= (1+tol)*lower bound.") |
|
|
268 cmd = "tol = 0.01; gmin = 1e-2; gmax = 1e+4;"; |
|
|
269 run_cmd |
|
|
270 cmd = "[AsysHinf,gmin,gmax] = hinfnorm(Asys,tol,gmin,gmax)" |
|
|
271 run_cmd |
|
|
272 disp("Check: look at max value of magntude Bode plot of Asys:"); |
|
|
273 [M,P,w] = bode(Asys); |
|
|
274 xlabel("Omega") |
|
|
275 ylabel("|Asys(j omega)| ") |
|
|
276 grid(); |
|
|
277 semilogx(w,M); |
|
|
278 disp(["Max magnitude is ",num2str(max(M)), ... |
|
|
279 ", compared with gmin=",num2str(gmin)," and gmax=", ... |
|
|
280 num2str(gmax),"."]) |
|
|
281 prompt |
|
|
282 disp("Example 2: unstable plant") |
|
|
283 cmd = "A = [0, 1; 2, 1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)"; |
|
|
284 run_cmd |
|
|
285 disp("Pack into system format:") |
|
4771
|
286 cmd = "Bsys = ss(A,B,C);"; |
|
3440
|
287 run_cmd |
|
|
288 disp("Evaluate with BsysH2 = hinfnorm(Bsys,tol,gmin,gmax)") |
|
|
289 BsysH2 = hinfnorm(Bsys,tol,gmin,gmax) |
|
|
290 disp(" ") |
|
|
291 disp("NOTICE: program returns a value without an error signal.") |
|
|
292 disp("") |
|
|
293 |
|
|
294 elseif (sel == 7) |
|
|
295 disp("Hinfinity optimal controller of a system: command = hinfsyn:") |
|
|
296 prompt |
|
|
297 help hinfsyn |
|
|
298 prompt |
|
|
299 disp("Example system: double integrator with output noise and") |
|
|
300 disp("input disturbance:") |
|
|
301 A = [0, 1; 0, 0] |
|
|
302 B1 = [0, 0; 1, 0] |
|
|
303 B2 = [0; 1] |
|
|
304 C1 = [1, 0; 0, 0] |
|
|
305 C2 = [1, 0] |
|
|
306 D11 = zeros(2); |
|
|
307 D12 = [0; 1]; |
|
|
308 D21 = [0, 1]; |
|
|
309 D22 = 0; |
|
|
310 D = [D11, D12; D21, D22] |
|
|
311 prompt |
|
|
312 disp("First: pack system:") |
|
4771
|
313 cmd="Asys = ss(A, [B1, B2], [C1; C2], D);"; |
|
3440
|
314 run_cmd |
|
|
315 prompt |
|
|
316 disp("Open loop multivariable Bode plot: (will take a moment)") |
|
|
317 cmd="bode(Asys);"; |
|
|
318 run_cmd |
|
|
319 prompt |
|
|
320 disp("Controller design command: (only need 1st two output arguments)") |
|
|
321 gmax = 1000 |
|
|
322 gmin = 0.1 |
|
|
323 gtol = 0.01 |
|
|
324 cmd="[K,gain] = hinfsyn(Asys,1,1,gmin,gmax,gtol);"; |
|
|
325 run_cmd |
|
|
326 disp("Check: close the loop and then compute h2norm:") |
|
|
327 prompt |
|
|
328 cmd="K_loop = sysgroup(Asys,K);"; |
|
|
329 run_cmd |
|
|
330 cmd = "Kcl = sysconnect(K_loop,[3,4],[4,3]);"; |
|
|
331 run_cmd |
|
|
332 cmd = "Kcl = sysprune(Kcl,[1,2],[1,2]);"; |
|
|
333 run_cmd |
|
|
334 cmd="gain_Kcl = hinfnorm(Kcl)"; |
|
|
335 run_cmd |
|
|
336 cmd="gain_err = gain_Kcl - gain"; |
|
|
337 run_cmd |
|
|
338 disp("Check: multivarible bode plot:") |
|
|
339 cmd="bode(Kcl);"; |
|
|
340 run_cmd |
|
|
341 prompt |
|
|
342 disp("Related functions: is_dgkf, is_controllable, is_stabilizable,") |
|
|
343 disp(" is_observable, is_detectable, buildssic") |
|
|
344 elseif (sel == 8) |
|
|
345 disp("Hinfinity optimal controller of MIMO system: command = hinfsyn:") |
|
|
346 prompt |
|
|
347 help hinfsyn |
|
|
348 prompt |
|
|
349 disp("Example system: Boeing 707-321 airspeed/pitch angle control") |
|
|
350 disp(" ") |
|
|
351 hinfdemo |
|
|
352 elseif (sel == 9) |
|
|
353 disp("Discrete time H-infinity control via bilinear transform"); |
|
|
354 prompt |
|
|
355 dhinfdemo |
|
|
356 elseif (sel == 10) |
|
|
357 return |
|
|
358 endif |
|
|
359 prompt |
|
|
360 endwhile |
|
6522
|
361 page_screen_output (save_val); |
|
3440
|
362 |
|
|
363 endfunction |
