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