Mercurial > hg > octave-thorsten

comparison scripts/control/dgkfdemo.m @ 3400:18366d37e7dd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1999-12-22 23:36:09 by jwe]
author jwe
date Wed, 22 Dec 1999 23:36:14 +0000
parents 0f515bc98460
children 8625164a0a39
comparison
equal deleted inserted replaced
3399:e665633c76af 3400:18366d37e7dd
28 28
29 save_val = page_screen_output; 29 save_val = page_screen_output;
30 page_screen_output = 1; 30 page_screen_output = 1;
31 while (1) 31 while (1)
32 clc 32 clc
33 menuopt=0; 33 sel = 0;
34 while(menuopt > 10 || menuopt < 1) 34 while (sel > 10 || sel < 1)
35 menuopt = menu('Octave H2/Hinfinity options demo', ... 35 sel = menu ("Octave H2/Hinfinity options demo",
36 'LQ regulator', ... 36 "LQ regulator",
37 'LG state estimator', ... 37 "LG state estimator",
38 'LQG optimal control design', ... 38 "LQG optimal control design",
39 'H2 gain of a system', ... 39 "H2 gain of a system",
40 'H2 optimal controller of a system', ... 40 "H2 optimal controller of a system",
41 'Hinf gain of a system', ... 41 "Hinf gain of a system",
42 'Hinf optimal controller of a SISO system', ... 42 "Hinf optimal controller of a SISO system",
43 'Hinf optimal controller of a MIMO system', ... 43 "Hinf optimal controller of a MIMO system",
44 'Discrete-time Hinf optimal control by bilinear transform', ... 44 "Discrete-time Hinf optimal control by bilinear transform",
45 'Return to main demo menu'); 45 "Return to main demo menu");
46 endwhile 46 endwhile
47 if (menuopt == 1) 47 if (sel == 1)
48 disp('Linear/Quadratic regulator design:') 48 disp("Linear/Quadratic regulator design:")
49 disp('Compute optimal state feedback via the lqr command...') 49 disp("Compute optimal state feedback via the lqr command...")
50 help lqr 50 help lqr
51 disp(' ') 51 disp(" ")
52 disp('Example:') 52 disp("Example:")
53 A = [0, 1; -2, -1] 53 A = [0, 1; -2, -1]
54 B = [0; 1] 54 B = [0; 1]
55 Q = [1, 0; 0, 0] 55 Q = [1, 0; 0, 0]
56 R = 1 56 R = 1
57 disp("Q = state penalty matrix; R = input penalty matrix") 57 disp("Q = state penalty matrix; R = input penalty matrix")
58 prompt 58 prompt
59 disp('Compute state feedback gain k, ARE solution P, and closed-loop') 59 disp("Compute state feedback gain k, ARE solution P, and closed-loop")
60 disp('poles as follows:'); 60 disp("poles as follows:");
61 cmd = "[k, p, e] = lqr(A,B,Q,R)"; 61 cmd = "[k, p, e] = lqr(A,B,Q,R)";
62 run_cmd 62 run_cmd
63 prompt 63 prompt
64 disp("A similar approach can be used for LTI discrete-time systems") 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).") 65 disp("by using the dlqr command in place of lqr (see LQG example).")
66 elseif (menuopt == 2) 66 elseif (sel == 2)
67 disp('Linear/Gaussian estimator design:') 67 disp("Linear/Gaussian estimator design:")
68 disp('Compute optimal state estimator via the lqe command...') 68 disp("Compute optimal state estimator via the lqe command...")
69 help lqe 69 help lqe
70 disp(' ') 70 disp(" ")
71 disp('Example:') 71 disp("Example:")
72 A = [0, 1; -2, -1] 72 A = [0, 1; -2, -1]
73 disp("disturbance entry matrix G") 73 disp("disturbance entry matrix G")
74 G = eye(2) 74 G = eye(2)
75 disp("Output measurement matrix C") 75 disp("Output measurement matrix C")
76 C = [0, 1] 76 C = [0, 1]
77 SigW = [1, 0; 0, 1] 77 SigW = [1, 0; 0, 1]
78 SigV = 1 78 SigV = 1
79 disp("SigW = input disturbance intensity matrix;") 79 disp("SigW = input disturbance intensity matrix;")
80 disp("SigV = measurement noise intensity matrix") 80 disp("SigV = measurement noise intensity matrix")
81 prompt 81 prompt
82 disp('Compute estimator feedback gain k, ARE solution P, and estimator') 82 disp("Compute estimator feedback gain k, ARE solution P, and estimator")
83 disp('poles via the command: ') 83 disp("poles via the command: ")
84 cmd = "[k, p, e] = lqe(A,G,C,SigW,SigV)"; 84 cmd = "[k, p, e] = lqe(A,G,C,SigW,SigV)";
85 run_cmd 85 run_cmd
86 disp("A similar approach can be used for LTI discrete-time systems") 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).") 87 disp("by using the dlqe command in place of lqe (see LQG example).")
88 elseif (menuopt == 3) 88 elseif (sel == 3)
89 disp('LQG optimal controller of a system:') 89 disp("LQG optimal controller of a system:")
90 disp('Input accepted as either A,B,C matrices or in system data structure form') 90 disp("Input accepted as either A,B,C matrices or in system data structure form")
91 disp('in both discrete and continuous time.') 91 disp("in both discrete and continuous time.")
92 disp("Example 1: continuous time design:") 92 disp("Example 1: continuous time design:")
93 prompt 93 prompt
94 help lqg 94 help lqg
95 disp("Example system") 95 disp("Example system")
96 A = [0, 1; .5, .5]; 96 A = [0, 1; .5, .5];
140 prompt 140 prompt
141 disp("Check: poles of Acl:") 141 disp("Check: poles of Acl:")
142 Acl_poles = sortcom(eig(Acl)) 142 Acl_poles = sortcom(eig(Acl))
143 disp("Predicted poles from design = union(Er,Ee)") 143 disp("Predicted poles from design = union(Er,Ee)")
144 pred_poles = sortcom([Er;Ee]) 144 pred_poles = sortcom([Er;Ee])
145 elseif (menuopt == 4) 145 elseif (sel == 4)
146 disp('H2 gain of a system: (Energy in impulse response)') 146 disp("H2 gain of a system: (Energy in impulse response)")
147 disp('Example 1: Stable plant:') 147 disp("Example 1: Stable plant:")
148 cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)"; 148 cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)";
149 run_cmd 149 run_cmd
150 disp("Put into Packed system form:") 150 disp("Put into Packed system form:")
151 cmd = "Asys = ss2sys(A,B,C);"; 151 cmd = "Asys = ss2sys(A,B,C);";
152 run_cmd 152 run_cmd
159 ht(ii) = C*expm(A*tt(ii))*B; 159 ht(ii) = C*expm(A*tt(ii))*B;
160 endfor 160 endfor
161 plot(tt,ht) 161 plot(tt,ht)
162 title("impulse response of example plant") 162 title("impulse response of example plant")
163 prompt 163 prompt
164 disp('Example 2: unstable plant') 164 disp("Example 2: unstable plant")
165 cmd = "A = [0, 1; 2, 1]"; 165 cmd = "A = [0, 1; 2, 1]";
166 eval(cmd); 166 eval(cmd);
167 cmd = "B = [0; 1]"; 167 cmd = "B = [0; 1]";
168 eval(cmd); 168 eval(cmd);
169 cmd = "C = [1, 0]"; 169 cmd = "C = [1, 0]";
170 eval(cmd); 170 eval(cmd);
171 cmd = "sys_poles = eig(A)"; 171 cmd = "sys_poles = eig(A)";
172 run_cmd 172 run_cmd
173 prompt 173 prompt
174 disp('Put into system data structure form:') 174 disp("Put into system data structure form:")
175 cmd="Bsys = ss2sys(A,B,C);"; 175 cmd="Bsys = ss2sys(A,B,C);";
176 run_cmd 176 run_cmd
177 disp('Evaluate 2-norm:') 177 disp("Evaluate 2-norm:")
178 cmd = "BsysH2 = h2norm(Bsys)"; 178 cmd = "BsysH2 = h2norm(Bsys)";
179 run_cmd 179 run_cmd
180 disp(' ') 180 disp(" ")
181 prompt('NOTICE: program returns a value without an error signal.') 181 prompt("NOTICE: program returns a value without an error signal.")
182 disp('') 182 disp("")
183 183
184 elseif (menuopt == 5) 184 elseif (sel == 5)
185 disp('H2 optimal controller of a system: command = h2syn:') 185 disp("H2 optimal controller of a system: command = h2syn:")
186 prompt 186 prompt
187 help h2syn 187 help h2syn
188 prompt 188 prompt
189 disp("Example system: double integrator with output noise and") 189 disp("Example system: double integrator with output noise and")
190 disp("input disturbance:") 190 disp("input disturbance:")
242 cmd="bode(Kcl);"; 242 cmd="bode(Kcl);";
243 run_cmd 243 run_cmd
244 prompt 244 prompt
245 disp("Related functions: is_dgkf, is_controllable, is_stabilizable,") 245 disp("Related functions: is_dgkf, is_controllable, is_stabilizable,")
246 disp(" is_observable, is_detectable") 246 disp(" is_observable, is_detectable")
247 elseif (menuopt == 6) 247 elseif (sel == 6)
248 disp('Hinfinity gain of a system: (max gain over all j-omega)') 248 disp("Hinfinity gain of a system: (max gain over all j-omega)")
249 disp('Example 1: Stable plant:') 249 disp("Example 1: Stable plant:")
250 cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)"; 250 cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)";
251 run_cmd 251 run_cmd
252 disp('Pack into system format:') 252 disp("Pack into system format:")
253 cmd = "Asys = ss2sys(A,B,C);"; 253 cmd = "Asys = ss2sys(A,B,C);";
254 run_cmd 254 run_cmd
255 disp('The infinity norm must be computed iteratively by') 255 disp("The infinity norm must be computed iteratively by")
256 disp('binary search. For this example, we select tolerance tol = 0.01, ') 256 disp("binary search. For this example, we select tolerance tol = 0.01, ")
257 disp('min gain gmin = 1e-2, max gain gmax=1e4.') 257 disp("min gain gmin = 1e-2, max gain gmax=1e4.")
258 disp('Search quits when upper bound <= (1+tol)*lower bound.') 258 disp("Search quits when upper bound <= (1+tol)*lower bound.")
259 cmd = "tol = 0.01; gmin = 1e-2; gmax = 1e+4;"; 259 cmd = "tol = 0.01; gmin = 1e-2; gmax = 1e+4;";
260 run_cmd 260 run_cmd
261 cmd = "[AsysHinf,gmin,gmax] = hinfnorm(Asys,tol,gmin,gmax)" 261 cmd = "[AsysHinf,gmin,gmax] = hinfnorm(Asys,tol,gmin,gmax)"
262 run_cmd 262 run_cmd
263 disp("Check: look at max value of magntude Bode plot of Asys:"); 263 disp("Check: look at max value of magntude Bode plot of Asys:");
264 [M,P,w] = bode(Asys); 264 [M,P,w] = bode(Asys);
265 xlabel('Omega') 265 xlabel("Omega")
266 ylabel('|Asys(j omega)| ') 266 ylabel("|Asys(j omega)| ")
267 grid(); 267 grid();
268 semilogx(w,M); 268 semilogx(w,M);
269 disp(["Max magnitude is ",num2str(max(M)), ... 269 disp(["Max magnitude is ",num2str(max(M)), ...
270 ", compared with gmin=",num2str(gmin)," and gmax=", ... 270 ", compared with gmin=",num2str(gmin)," and gmax=", ...
271 num2str(gmax),"."]) 271 num2str(gmax),"."])
272 prompt 272 prompt
273 disp('Example 2: unstable plant') 273 disp("Example 2: unstable plant")
274 cmd = "A = [0, 1; 2, 1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)"; 274 cmd = "A = [0, 1; 2, 1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)";
275 run_cmd 275 run_cmd
276 disp("Pack into system format:") 276 disp("Pack into system format:")
277 cmd = "Bsys = ss2sys(A,B,C);"; 277 cmd = "Bsys = ss2sys(A,B,C);";
278 run_cmd 278 run_cmd
279 disp('Evaluate with BsysH2 = hinfnorm(Bsys,tol,gmin,gmax)') 279 disp("Evaluate with BsysH2 = hinfnorm(Bsys,tol,gmin,gmax)")
280 BsysH2 = hinfnorm(Bsys,tol,gmin,gmax) 280 BsysH2 = hinfnorm(Bsys,tol,gmin,gmax)
281 disp(' ') 281 disp(" ")
282 disp('NOTICE: program returns a value without an error signal.') 282 disp("NOTICE: program returns a value without an error signal.")
283 disp('') 283 disp("")
284 284
285 elseif (menuopt == 7) 285 elseif (sel == 7)
286 disp('Hinfinity optimal controller of a system: command = hinfsyn:') 286 disp("Hinfinity optimal controller of a system: command = hinfsyn:")
287 prompt 287 prompt
288 help hinfsyn 288 help hinfsyn
289 prompt 289 prompt
290 disp("Example system: double integrator with output noise and") 290 disp("Example system: double integrator with output noise and")
291 disp("input disturbance:") 291 disp("input disturbance:")
330 cmd="bode(Kcl);"; 330 cmd="bode(Kcl);";
331 run_cmd 331 run_cmd
332 prompt 332 prompt
333 disp("Related functions: is_dgkf, is_controllable, is_stabilizable,") 333 disp("Related functions: is_dgkf, is_controllable, is_stabilizable,")
334 disp(" is_observable, is_detectable, buildssic") 334 disp(" is_observable, is_detectable, buildssic")
335 elseif (menuopt == 8) 335 elseif (sel == 8)
336 disp('Hinfinity optimal controller of MIMO system: command = hinfsyn:') 336 disp("Hinfinity optimal controller of MIMO system: command = hinfsyn:")
337 prompt 337 prompt
338 help hinfsyn 338 help hinfsyn
339 prompt 339 prompt
340 disp("Example system: Boeing 707-321 airspeed/pitch angle control") 340 disp("Example system: Boeing 707-321 airspeed/pitch angle control")
341 disp(" ") 341 disp(" ")
342 hinfdemo 342 hinfdemo
343 elseif (menuopt == 9) 343 elseif (sel == 9)
344 disp("Discrete time H-infinity control via bilinear transform"); 344 disp("Discrete time H-infinity control via bilinear transform");
345 prompt 345 prompt
346 dhinfdemo 346 dhinfdemo
347 elseif (menuopt == 10) 347 elseif (sel == 10)
348 return 348 return
349 endif 349 endif
350 prompt 350 prompt
351 endwhile 351 endwhile
352 page_screen_output = save_val; 352 page_screen_output = save_val;