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