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