home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / HMATDEMO.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  14.0 KB  |  401 lines
  1. %
  2. % HMATDEMO Fighter Design Demonstration
  3. %
  4.  
  5. % R. Y. Chiang & M. G. Safonov 3/88
  6. % Copyright (c) 1988 by the MathWorks, Inc.
  7. % All Rights Reserved.
  8. % ---------------------------------------------------------------
  9. clear
  10. clc
  11. disp(' ')
  12. disp('  ')
  13. disp('  ')
  14. disp('                << Demo #2: MIMO Fighter Design Example >>')
  15. disp(' ')
  16. disp(' ')
  17. disp('                        ----')
  18. disp('                         \  \          ----')
  19. disp('                       == \ = \         \  \')
  20. disp('                          /     \        \  \')
  21. disp('                    / ---------------------------- \')
  22. disp('                   |         FIGHTER OF 1990        >>>> -----')
  23. disp('                    \ ---------------------------- /')
  24. disp('                          \     /        /  /')
  25. disp('                       == / = /         /  /')
  26. disp('                         /  /          ----')
  27. disp('                        ----')
  28. format short e
  29. ag =[
  30. -2.2567e-02  -3.6617e+01  -1.8897e+01  -3.2090e+01   3.2509e+00  -7.6257e-01;
  31. 9.2572e-05  -1.8997e+00   9.8312e-01  -7.2562e-04  -1.7080e-01  -4.9652e-03;
  32. 1.2338e-02   1.1720e+01  -2.6316e+00   8.7582e-04  -3.1604e+01   2.2396e+01;
  33. 0            0   1.0000e+00            0            0            0;
  34. 0            0            0            0  -3.0000e+01            0;
  35. 0            0            0            0            0  -3.0000e+01];
  36. bg = [0     0;
  37.       0     0;
  38.       0     0;
  39.       0     0;
  40.      30     0;
  41.       0    30];
  42. cg = [0     1     0     0     0     0;
  43.       0     0     0     1     0     0];
  44. dg = [0     0;
  45.       0     0];
  46. disp('  ')
  47. disp('  ')
  48. disp('                                 (strike a key to continue ...)')
  49. pause
  50. clc
  51. disp('  ')
  52. disp('   State-Space of the fighter pitch axis (trimed @ 25000 ft, 0.9 mach):')
  53. ag 
  54. bg
  55. disp('                                 (strike a key to continue ...)')
  56. pause
  57. clc
  58. cg
  59. dg
  60. disp('                                 (strike a key to continue ...)')
  61. pause
  62. clc
  63. disp('  ')
  64. disp('  ')
  65. disp('    Poles of the Plant (unstable Phugoid modes):')
  66. disp('  ')
  67. disp('    ---------------------------------------------------------------')
  68. disp('       poleg = eig(ag)      % Computing the poles of the plant')
  69. disp('    ---------------------------------------------------------------')
  70. poleg  = eig(ag)       
  71. disp(' ')
  72. disp(' ')
  73. disp('                                 (strike a key to continue ...)')
  74. pause
  75. clc
  76. disp(' ')
  77. disp('  ')
  78. disp('    Transmission Zeros of the Plant (minimum phase):')
  79. disp('  ')
  80. disp('    -------------------------------------------------------------------')
  81. disp('      tzerog = tzero(ag,bg,cg,dg)   % Computing the transmission zeros')
  82. disp('    -------------------------------------------------------------------')
  83. tzerog = tzero(ag,bg,cg,dg)      % Computing the transmission zeros
  84. disp('  ')
  85. disp(' ')
  86. disp(' ')
  87. disp('                                 (strike a key to continue ...)')
  88. pause
  89. clc
  90. disp(' ')
  91. disp('      - - - Computing SV Bode plot of the open loop plant - - -')
  92. w = logspace(-3,5,50);
  93. svg = sigma(ag,bg,cg,dg,1,w); svg = 20*log10(svg);
  94. disp(' ')
  95. disp(' ')
  96. disp(' ')
  97. disp(' ')
  98. disp('                (strike a key to see the SV Bode plot of G(s) ...)')
  99. pause
  100. clg
  101. semilogx(w,svg)
  102. title('MIMO Fighter Pitch Axis Open Loop')
  103. xlabel('Frequency - Rad/Sec')
  104. ylabel('SV - db')
  105. grid
  106. pause
  107. clc
  108. disp('  ')
  109. disp('                    << Design Specifications >>  ')
  110. disp(' ')
  111. disp('      1). Robustness Spec. : -2 roll-off, -20 db @ 100 Rad/Sec.')
  112. disp('          Associated Weighting:')
  113. disp('     ')
  114. disp('                    -1     1000')
  115. disp('                  W3(s) = ------ * I       (fixd)')
  116. disp('                            s^2     2x2')
  117. disp('   ')
  118. disp(' ')
  119. disp('      2). Performance Spec.: minimizing the sensitivity function')
  120. disp('                             as much as possible.')
  121. disp('          Associated Weighting:')
  122. disp(' ') 
  123. disp('                    -1       -1   s + 0.01')
  124. disp('                  W1(s) = Gam  * ----------- *  I')
  125. disp('                                     1           2x2')
  126. disp('   ')
  127. disp('          where "Gam" in our design goes from 1 --> 8.4 --> 16.8.')
  128. k = 1000; tau = 5.0000e-04; mn = [2 2];
  129. nuw3i = [0 0 k]; dnw3i = [1 0 0]; 
  130. svw3i = bode(nuw3i,dnw3i,w); svw3i = 20*log10(svw3i);
  131. nuw1i = [1.0000e+00   1.0000e-02]; dnw1i =[0 1];
  132. svw1i = bode(nuw1i,dnw1i,w); svw1i = 20*log10(svw1i);
  133. disp('  ')
  134. disp('  ')
  135. disp('                   (strike a key to see the plot of the weightings ...)')
  136. pause
  137. semilogx(w,svw1i,w,svw3i)
  138. grid
  139. title('MIMO Fighter Design Example -- Design Specifications')
  140. xlabel('Frequency - Rad/Sec')
  141. ylabel('1/W1 & 1/W3 - db')
  142. text(0.003,-70,'Sensitivity Spec.')
  143. text(0.003,-100,'1/W1(s)')
  144. text(200,-20,'Robustness Spec.')
  145. text(1000,-50,'1/W3(s)')
  146. pause
  147. clc
  148. disp('                      << Problem Formulation >>')
  149. disp(' ')
  150. disp('    Form an augmented plant P(s) with these two weighting functions:')
  151. disp(' ')
  152. disp('               1). W1 penalizing error signal "e"')
  153. disp('  ')
  154. disp('               2). W3 penalizing plant output "y"')
  155. disp('  ')
  156. disp('    and find a stabilizing controller F(s) such that the H2 or H-Inf norm')
  157. disp('    of TF Ty1u1 is minimized and its H-Inf norm is less than one, i.e.')
  158. disp('  ')
  159. disp('                        min |Ty1u1|   < 1,')
  160. disp('                        F(s)       inf')
  161. disp(' ')
  162. disp('      where ')
  163. disp('                       |               -1|')                         
  164. disp('               Ty1u1 = | Gam*W1*(I + GF) |  = | Gam * W1 * S  |')
  165. disp('                       |               -1|    |  W3 * (I - S) |')
  166. disp('                       |  W3*GF*(I + GF) |')
  167. disp('  ')
  168. disp('  ')
  169. disp('  ')
  170. disp('                                 (strike a key to continue ...)')
  171. pause
  172. clc
  173. disp('  ')
  174. disp(' ')
  175. disp('                              << DESIGN PROCEDURE >>')
  176. disp('  ')
  177. disp('            * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *')
  178. disp('            *    [Step 1]. Do plant augmentation (run AUGTF.M or      *')
  179. disp('            *              AUGSS.M)                                   *')
  180. disp('            *                                                         *')
  181. disp('            *    [Step 2]. Do H2 synthesis (run H2LQG.M)              *')
  182. disp('            *                                                         *')
  183. disp('            *    [Step 3]. Redo the plant augmentation for a          *')
  184. disp('            *              new "Gam" --> 8.4 and rerun H2LQG.M        *')
  185. disp('            *                                                         *')
  186. disp('            *    [Step 4]. Redo the plant augmentation for a          *')
  187. disp('            *              higher "Gam" --> 16.8, then run HINF.M     *')
  188. disp('            * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *')
  189. disp('  ')
  190. disp('  ')
  191. disp('                                 (strike a key to continue ...)')
  192. pause
  193. clc
  194. disp('  ')
  195. disp('  ')
  196. disp('           Assign the cost coefficient "Gam" --> 1 ')
  197. disp('       ')
  198. disp('           this will serve as the baseline design ....')
  199. disp('  ')
  200. disp('            ----------------------------------------------------------')
  201. disp('             % Plant augmentation of the fighter:')
  202. disp('             ss_g = mksys(ag,bg,cg,dg);')
  203. disp('             w1 = [Gam*dnw1i;nuw1i;Gam*dnw1i;nuw1i];')
  204. disp('             w2 = [];')
  205. disp('             w3 = [0 1 0 0;0 0 0 k;tau 1 0 0;0 0 0 k];')
  206. disp('             [TSS_] = augtf(ss_g,w1,w2,w3);')
  207. disp('            ----------------------------------------------------------')
  208. disp('   ')
  209. disp(' ')
  210. Gam = input('           Input the cost coefficient "Gam" = ');
  211. ss_g = mksys(ag,bg,cg,dg);
  212. w1 = [Gam*dnw1i;nuw1i;Gam*dnw1i;nuw1i];
  213. w2 = [];
  214. w3 = [0 1 0 0;0 0 0 k;tau 1 0 0;0 0 0 k];
  215. [TSS_]=augtf(ss_g,w1,w2,w3);
  216. disp(' ')
  217. disp('     - - - State-Space TSS_:=(A,B1,B2,..,D21,D22) is ready for')
  218. disp('           the Small-Gain problem - - -')
  219. disp('  ')
  220. disp('  ')
  221. disp('                                 (strike a key to continue ...)')
  222. pause
  223. clc
  224. disp(' ')
  225. disp(' ')
  226. disp('    --------------------------------------------------------------')
  227. disp('     [ss_cp,ss_cl] = h2lqg(TSS_);  % Running H2LQG.M')
  228. disp('     [acp,bcp,ccp,dcp] = branch(ss_cp);')
  229. disp('     [acl,bcl,ccl,dcl] = branch(ss_cl);')
  230. disp('    --------------------------------------------------------------')
  231. [ss_cp,ss_cl] = h2lqg(TSS_);
  232. [acp,bcp,ccp,dcp] = branch(ss_cp);
  233. [acl,bcl,ccl,dcl] = branch(ss_cl);
  234. disp('  ')
  235. disp('  ')
  236. disp('                                 (strike a key to continue ...)')
  237. pause
  238. pltopt           % Preparing singular values for plotting
  239. svw1i1 = svw1i; h2svs1 = svs; h2svt1 = svt; h2svtt1 = svtt;
  240. disp('  ')
  241. disp('  ')
  242. disp('                                 (strike a key to continue ...)')
  243. pause
  244. clc      
  245. disp('  ')
  246. disp('  ')
  247. disp('    After a few iterations without changing the dynamics of the weights,')
  248. disp(' ')
  249. disp('    we found a new Gam of 8.4 can push the H2 cost function close to its')
  250. disp('  ')
  251. disp('    "shaping limit". ')
  252. disp('  ')
  253. disp('  ')
  254. disp('            Input "Gam" --> 8.4, and try H2LQG again .....')
  255. disp('  ')
  256. disp('  ')
  257. Gam = input('           Input the cost coefficient "Gam" = ');
  258. w1 = [Gam*dnw1i;nuw1i;Gam*dnw1i;nuw1i];
  259. [TSS_]=augtf(ss_g,w1,w2,w3);
  260. disp(' ')
  261. disp('     - - - New state-space TSS_ is ready for')
  262. disp('           the Small-Gain problem - - -')
  263. disp('  ')
  264. disp('                                 (strike a key to continue ...)')
  265. pause
  266. clc
  267. disp(' ')
  268. disp(' ')
  269. disp('    --------------------------------------------------------------')
  270. disp('     [ss_cp,ss_cl] = h2lqg(TSS_);  % Running H2LQG.M')
  271. disp('     [acp,bcp,ccp,dcp] = branch(ss_cp);')
  272. disp('     [acl,bcl,ccl,dcl] = branch(ss_cl);')
  273. disp('    --------------------------------------------------------------')
  274. [ss_cp,ss_cl] = h2lqg(TSS_);
  275. [acp,bcp,ccp,dcp] = branch(ss_cp);
  276. [acl,bcl,ccl,dcl] = branch(ss_cl);
  277. disp('  ')
  278. disp('  ')
  279. disp('                                 (strike a key to continue ...)')
  280. pause
  281. pltopt
  282. svw1i2 = svw1i; h2svs2 = svs; h2svt2 = svt; h2svtt2 = svtt;
  283. disp('  ')
  284. disp('  ')
  285. disp('                                 (strike a key to continue ...)')
  286. pause
  287. clc
  288. disp('  ')
  289. disp('  ')
  290. disp('               Now let us try H-Inf synthesis .....')
  291. disp('   ')
  292. disp('    After a few iterations on the parameter "Gam", we found "Gam" can be')
  293. disp(' ')
  294. disp('    increased to 16.8, which is twice of the H2 approach !!  That is')
  295. disp('   ')
  296. disp('    using "Gam-iteration" one can squeeze more for a particular design')
  297. disp(' ')
  298. disp('    spec. than H2, hence, shape the singular values "exactly" to the pre-')
  299. disp('  ')
  300. disp('    specified frequency domain specifications.')
  301. disp('  ')
  302. disp('    Note that for H-Inf synthesis, the D11 matrix (i.e. W1(inf)) can')
  303. disp('  ')
  304. disp('    be NONZERO now (--> denominator of W1 = [0.01 1];).')
  305. disp(' ')
  306. disp('  ')
  307. disp('                                 (strike a key to continue ...)')
  308. pause
  309. clc
  310. dnw1i = [0.01 1];
  311. disp('            ----------------------------------------------------------')
  312. disp('             % Plant augmentation of the fighter:')
  313. disp('             w1 = [Gam*dnw1i;nuw1i;Gam*dnw1i;nuw1i];')
  314. disp('             [TSS_]=augtf(ss_g,w1,w2,w3);')
  315. disp('            ----------------------------------------------------------')
  316. disp('   ')
  317. disp('             Input "Gam" --> 16.8, and try HINF.M ....')
  318. disp(' ')
  319. Gam = input('           Input the cost coefficient "Gam" = ');
  320. w1 = [Gam*dnw1i;nuw1i;Gam*dnw1i;nuw1i];
  321. [TSS_]=augtf(ss_g,w1,w2,w3);
  322. disp(' ')
  323. disp('     - - - State-Space TSS_ is ready for')
  324. disp('           the Small-Gain problem - - -')
  325. disp('  ')
  326. disp('  ')
  327. disp('                                 (strike a key to continue ...)')
  328. pause
  329. clc
  330. disp(' ')
  331. disp(' ')
  332. disp('   --------------------------------------------------------------- ')
  333. disp('    [ss_cp,ss_cl,hinfo] = hinf(TSS_);')
  334. disp('    [acp,bcp,ccp,dcp] = branch(ss_cp);')
  335. disp('    [acl,bcl,ccl,dcl] = branch(ss_cl);')
  336. disp('                    % Running file HINF.M for H-Inf optimization')
  337. disp('   ---------------------------------------------------------------')
  338. [ss_cp,ss_cl,hinfo] = hinf(TSS_);
  339. [acp,bcp,ccp,dcp] = branch(ss_cp);
  340. [acl,bcl,ccl,dcl] = branch(ss_cl);
  341. disp('  ')
  342. disp('  ')
  343. disp('                                 (strike a key to continue ...)')
  344. pause
  345. pltopt
  346. svw1i3 = svw1i; hinfsvs = svs; hinfsvt = svt; hinfsvtt = svtt;
  347. disp('  ')
  348. disp('  ')
  349. disp('                  (strike a key to see the plots of the comparison ...)')
  350. pause
  351. clg
  352. semilogx(w,svw1i1,w,h2svs1,w,svw1i2,w,h2svs2,w,svw1i3,w,hinfsvs)
  353. title('H2 & H-Inf Fighter Design -- 1/W1 & Sensitivity Func.')
  354. xlabel('Frequency - Rad/Sec')
  355. ylabel('SV - db')
  356. grid
  357. text(0.002,30,'H2 (Gam = 1) ---> H2 (Gam = 8.4) ---> H-Inf (Gam = 16.8)')
  358. pause
  359. semilogx(w,svw3i,w,h2svt1,w,h2svt2,w,hinfsvt)
  360. title('H2 & H-Inf Fighter Design -- 1/W3 & Comp. Sens. Func.')
  361. xlabel('Frequency - Rad/Sec')
  362. ylabel('SV - db')
  363. grid
  364. text(0.005,100,'H2 (Gam = 1) ---> H2 (Gam = 8.4) ---> H-Inf (Gam = 16.8)')
  365. pause
  366. semilogx(w,h2svtt1,w,h2svtt2,w,hinfsvtt)
  367. title('H2 & H-Inf Fighter Design -- Cost function Ty1u1')
  368. xlabel('Frequency - Rad/Sec')
  369. ylabel('SV - db')
  370. grid
  371. text(0.002,-25,'H2 (Gam = 1) ---> H2 (Gam = 8.4) ---> H-Inf (Gam = 16.8)')
  372. pause
  373. clc
  374. disp('  ')
  375. disp('  ')
  376. disp('             << 8-State H-Inf Controller (Gam = 16.8) >>')
  377. disp('    ')
  378. disp('    Poles of the Controller :')
  379. polecp = eig(acp)
  380. disp(' ')
  381. disp('                                 (strike a key to continue ...)')
  382. pause
  383. disp(' ')
  384. disp('    Transmission Zeros of the Controller :')
  385. tzcp = tzero(acp,bcp,ccp,dcp)
  386. disp(' ')
  387. disp('                                 (strike a key to continue ...)')
  388. pause
  389. svcp = sigma(acp,bcp,ccp,dcp,1,w); svcp = 20*log10(svcp);
  390. semilogx(w,svcp)
  391. title('8-State H-Inf Controller (Gam = 16.8)')
  392. xlabel('Rad/Sec')
  393. ylabel('SV (db)')
  394. pause
  395. clc
  396. disp('    Poles of the Cost Function :')
  397. poletyu = eig(acl)
  398. disp('                                 (strike a key to continue ...)')
  399. pause
  400. %
  401. % ------ End of HIMATDEMO.