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