home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / LINFDEMO.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  12.5 KB  |  361 lines
  1. %
  2. % LINFDEMO Large space structure design demonstration using LINF
  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. clc
  10. disp('      << Linfdemo: MIMO Large Space Structure Design Example >>')
  11. disp('                                          __________')
  12. disp('         Secondary Mirror ---->   -ooo-         ^')
  13. disp('                                 / \ / \        |')
  14. disp('                                 |  X  |        |')
  15. disp('                                 | / \ |        |')
  16. disp('                                 -------        |')
  17. disp('                                 | \ / |        |')
  18. disp('                                 |  X  |        |')
  19. disp('                                 | / \ |         ')
  20. disp('                  Lense ------>  ---O---     7.4 Meters')
  21. disp('                                |       |        ')
  22. disp('                               |  \   /  |      |')
  23. disp('                               |   \ /   |      |')
  24. disp('                               |    X    |      |')
  25. disp('                              |    / \    |     |')
  26. disp('                              |   /   \   |     |')
  27. disp('                             /   /     \   \    |')
  28. disp('        Primary Mirror -->  (_OOOOOO____\___)   |')
  29. disp('                             \             /    |')
  30. disp('                              \___________/ ____v___')
  31. disp(' ')
  32. disp('                                 (strike a key to continue ...)')
  33. pause
  34. format short e
  35. ag =[-9.9005e-01   4.7491e-04   4.8994e-01   1.9219e+00;
  36.       9.0643e-04  -9.8765e-01   1.9010e+00  -4.9180e-01;
  37.      -4.9605e-01  -1.9005e+00  -3.1170e+02   4.9716e+00;
  38.      -1.9215e+00   4.9069e-01  -7.7879e+00  -3.9831e+02];
  39. bg = [7.8273e-01  -6.1398e-01;
  40.       6.1298e-01   7.8260e-01;
  41.       7.8349e-01   5.9597e-01;
  42.       6.0685e-01  -7.8784e-01];
  43. cg = [7.8291e-01   6.1279e-01  -7.8163e-01  -6.0610e-01;
  44.      -6.1438e-01   7.8200e-01  -5.9841e-01   7.8842e-01];
  45. dg = [0     0;
  46.       0     0];
  47. clc
  48. disp('  ')
  49. disp('   State-Space of the Large Space Structure:')
  50. disp('   (after colocated rate feedback and model reduction)')
  51. ag 
  52. bg
  53. disp('                                 (strike a key to continue ...)')
  54. pause
  55. clc
  56. cg
  57. dg
  58. disp('                                 (strike a key to continue ...)')
  59. pause
  60. clc
  61. disp('  ')
  62. disp('  ')
  63. disp('    Poles of the Plant (stable):')
  64. disp('  ')
  65. disp('    ---------------------------------------------------------------')
  66. disp('       poleg = eig(ag)      % Computing the poles of the plant')
  67. disp('    ---------------------------------------------------------------')
  68. poleg  = eig(ag)       
  69. disp(' ')
  70. disp(' ')
  71. disp('                                 (strike a key to continue ...)')
  72. pause
  73. clc
  74. disp('  ')
  75. disp('    ')
  76. disp('    Transmission Zeros of the Plant (minimum phase):')
  77. disp('  ')
  78. disp('    -------------------------------------------------------------------')
  79. disp('      tzerog = tzero(ag,bg,cg,dg)   % Computing the transmission zeros')
  80. disp('    -------------------------------------------------------------------')
  81. tzerog = tzero(ag,bg,cg,dg)      % Computing the transmission zeros
  82. disp('  ')
  83. disp(' ')
  84. disp(' ')
  85. disp('                                 (strike a key to continue ...)')
  86. pause
  87. clc
  88. disp(' ')
  89. disp(' - - - Computing singular value Bode plot of the open loop plant - - -')
  90. w = logspace(-3,5,100);
  91. svg = sigma(ag,bg,cg,dg,1,w); svg = 20*log10(svg);
  92. disp('  ')
  93. disp(' ')
  94. disp(' ')
  95. disp('                       (strike a key to see the open loop plant ...)')
  96. pause
  97. subplot
  98. semilogx(w,svg)
  99. title('MIMO Large Space Structure Open Loop')
  100. xlabel('Frequency - Rad/Sec')
  101. ylabel('SV - db')
  102. grid
  103. pause
  104. clc
  105. disp('  ')
  106. disp('                         << Design Specifications >>  ')
  107. disp(' ')
  108. disp('      1). Robustness Spec. : closed loop bandwidth -- 2000 r/s (300 hz)')
  109. disp('          Associated Weighting:')
  110. disp('     ')
  111. disp('                    -1     2000')
  112. disp('                  W3(s) = ------ * I       (fixd)')
  113. disp('                             s      2x2')
  114. disp('   ')
  115. disp(' ')
  116. disp('      2). Performance Spec.: sensitivity reduction of at least 100:1')
  117. disp('                             up to approx. 100 r/s')
  118. disp('          Associated Weighting:')
  119. disp(' ')          
  120. disp('                    -1       -1   0.01(1 + s/100)^2')
  121. disp('                  W1(s) = Gam  * --------------------- *  I')
  122. disp('                                     (1 + s/5000)^2        2x2')
  123. disp('   ')
  124. disp('          where "Gam" in our design goes from 1 --> 1.5.')
  125. k = 2000; 
  126. nuw3i = [0 k]; dnw3i = [1 0];
  127. svw3i = bode(nuw3i,dnw3i,w); svw3i = 20*log10(svw3i);
  128. nuw1i = 0.01*conv([1/100 1],[1/100 1]); dnw1i = conv([1/5000 1],[1/5000 1]);
  129. svw1i = bode(nuw1i,dnw1i,w); svw1i = 20*log10(svw1i);
  130. disp('  ')
  131. disp('  ')
  132. disp('         (strike a key to see the plot of weighting functions ...)')
  133. pause
  134. subplot
  135. axis([0 5 -40 40])
  136. semilogx(w,svw1i,w,svw3i)
  137. grid
  138. title('MIMO LSS Design Example -- Design Specifications')
  139. xlabel('Frequency - Rad/Sec')
  140. ylabel('1/W1 & 1/W3 - db')
  141. text(2,-20,'Sensitivity Spec.-- 1/W1(s)')
  142. text(2,10,'Robustness Spec.-- 1/W3(s)')
  143. pause
  144. axis
  145. clc
  146. disp('                      << Problem Formulation >>')
  147. disp(' ')
  148. disp('      Form an augmented plant P(s) with these two weighting functions:')
  149. disp(' ')
  150. disp('                 1). W1 penalizing error signal "e"')
  151. disp('  ')
  152. disp('                 2). W3 penalizing plant output "y"')
  153. disp(' ')
  154. disp('      and find a stabilizing controller F(s) such that the Hinf-norm')
  155. disp('      of TF Ty1u1 is minimized and less than one, i.e.')
  156. disp('  ')
  157. disp('                         min |Ty1u1|   < 1,')
  158. disp('                         F(s)       inf')
  159. disp(' ')
  160. disp('      where ')
  161. disp('                       |               -1|')
  162. disp('               Ty1u1 = | Gam*W1*(I + GF) |  = | Gam * W1 * S  |')
  163. disp('                       |               -1|    |               |')  
  164. disp('                       |  W3*GF*(I + GF) |    |  W3 * (I - S) |')
  165. disp('  ')
  166. disp('  ')
  167. disp('  ')
  168. disp('                                 (strike a key to continue ...)')
  169. pause
  170. clc
  171. disp('  ')
  172. disp(' ')
  173. disp('                              << DESIGN PROCEDURE >>')
  174. disp('  ')
  175. disp('            * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *')
  176. disp('            *    [Step 1]. Do plant augmentation (run AUGSS.M or      *')
  177. disp('            *              AUGTF.M)                                   *')
  178. disp('            *                                                         *')
  179. disp('            *    [Step 2]. Do H-inf synthesis (run LINF.M)            *')
  180. disp('            *                                                         *')
  181. disp('            *    [Step 3]. Redo the plant augmentation by setting     *')
  182. disp('            *              a new "Gam" --> 1.5 and rerun LINF.M       *')
  183. disp('            * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *')
  184. disp('  ')
  185. disp('  ')
  186. disp('                                 (strike a key to continue ...)')
  187. pause
  188. clc
  189. disp('  ')
  190. disp('  ')
  191. disp('  ')
  192. disp('           Assign the cost coefficient "Gam" --> 1 ')
  193. disp('       ')
  194. disp('           this will serve as the baseline design ....')
  195. disp('  ')
  196. disp(' ')
  197. Gam = input('             Input the cost coefficient "Gam" = ');
  198. disp('  ')
  199. disp('     ------------------------------------------------------------------')
  200. disp('         % Plant augmentation of the LSS:')
  201. disp('           ss_g = mksys(ag,bg,cg,dg);')
  202. disp('           w1 = [Gam*dnw1i;nuw1i;Gam*dnw1i;nuw1i];')
  203. disp('           w2 = [];')
  204. disp('           w3 = [dnw3i;nuw3i;dnw3i;nuw3i];')
  205. disp('           TSS_  = augtf(ss_g,w1,w2,w3);')
  206. disp('           [A,B1,B2,C1,C2,D11,D12,D21,D22] = branch(TSS_);')
  207. disp('     ------------------------------------------------------------------')
  208. ss_g = mksys(ag,bg,cg,dg); 
  209. w1 = [Gam*dnw1i;nuw1i;Gam*dnw1i;nuw1i];
  210. w2 = [];
  211. w3 = [dnw3i;nuw3i;dnw3i;nuw3i];
  212. TSS_ = augtf(ss_g,w1,w2,w3);
  213. [A,B1,B2,C1,C2,D11,D12,D21,D22] = branch(TSS_);
  214. [aa,bb,cc,mm,tt] = obalreal(A,[B1,B2],[C1;C2]);
  215. A = aa; B1 = bb(:,1:2); B2 = bb(:,3:4); C1 = cc(1:4,:); C2 = cc(5:6,:);
  216. disp('  ')
  217. disp('     - - - State-Space (A,B1,B2,C1,C2,D11,D12,D21,D22) is ready for')
  218. disp('           the Small-Gain problem - - -')
  219. disp(' ')
  220. disp(' ')
  221. disp('  ')
  222. disp('                                 (strike a key to continue ...)')
  223. pause
  224. clc
  225. disp(' ')
  226. disp(' ')
  227. disp('  ')
  228. disp('    ---------------------------------------------------------------')
  229. disp('     linf      % Running script file LINF.M for H-inf optimization')
  230. disp('    ---------------------------------------------------------------')
  231. linf
  232. syscp = [acp,bcp;ccp dcp]; xcp = size(acp)*[1;0];
  233. syscl = [acl,bcl;ccl dcl]; xcl = size(acl)*[1;0];
  234. disp('  ')
  235. disp('  ')
  236. disp('                                 (strike a key to continue ...)')
  237. pause
  238. pltopt           % Preparing singular values for plotting
  239. svw1i1 = svw1i; hsvs1 = svs; hsvt1 = svt; hsvtt1 = svtt;
  240. disp('  ')
  241. disp('  ')
  242. disp('                                 (strike a key to continue ...)')
  243. pause
  244. clc      
  245. disp('  ')
  246. disp('  ')
  247. disp('  ')
  248. disp('     After a few iterations, we found a new Gam of 1.5 can push the')
  249. disp('  ')
  250. disp('     H-inf cost function close to its limit. ')
  251. disp('  ')
  252. disp('  ')
  253. disp('            Input "Gam" --> 1.5, and try LINF again .....')
  254. disp('  ')
  255. disp('  ')
  256. Gam = input('             Input the cost coefficient "Gam" = ');
  257. disp('  ')
  258. disp('     ------------------------------------------------------------------')
  259. disp('         % Adjust plant augmentation:')
  260. disp('           w1 = [Gam*dnw1i;nuw1i;Gam*dnw1i;nuw1i];')
  261. disp('           TSS_ = augtf(ss_g,w1,w2,w3);')
  262. disp('           [A,B1,B2,C1,C2,D11,D12,D21,D22] = branch(TSS_);')
  263. disp('     ------------------------------------------------------------------')
  264. w1 = [Gam*dnw1i;nuw1i;Gam*dnw1i;nuw1i];
  265. TSS_ = augtf(ss_g,w1,w2,w3);
  266. [A,B1,B2,C1,C2,D11,D12,D21,D22] = branch(TSS_);
  267. [aa,bb,cc,mm,tt] = obalreal(A,[B1,B2],[C1;C2]);
  268. A = aa; B1 = bb(:,1:2); B2 = bb(:,3:4); C1 = cc(1:4,:); C2 = cc(5:6,:);
  269. disp('  ')
  270. disp('     - - - State-Space (A,B1,B2,C1,C2,D11,D12,D21,D22) is ready for')
  271. disp('           the Small-Gain problem - - -')
  272. disp('  ')
  273. disp('  ')
  274. disp('                                 (strike a key to continue ...)')
  275. pause
  276. linf
  277. syscp = [acp,bcp;ccp dcp]; xcp = size(acp)*[1;0];
  278. syscl = [acl,bcl;ccl dcl]; xcl = size(acl)*[1;0];
  279. disp('  ')
  280. disp('  ')
  281. disp('                                 (strike a key to continue ...)')
  282. pause
  283. pltopt
  284. svw1i2 = svw1i; hsvs2 = svs; hsvt2 = svt; hsvtt2 = svtt;
  285. disp('  ')
  286. disp('  ')
  287. disp('                             (strike a key to see the plots ...)')
  288. pause
  289. semilogx(w,svw1i1,w,hsvs1,w,svw1i2,w,hsvs2)
  290. title('H-inf LSS Design -- W1 & Sensitivity Func.')
  291. xlabel('Frequency - Rad/Sec')
  292. ylabel('SV - db')
  293. grid
  294. text(0.01,0,'H-inf (Gam = 1) ---> H-inf (Gam = 1.5)')
  295. pause
  296. axis([0 5 -40 40])
  297. subplot
  298. semilogx(w,svw3i,w,hsvt1,w,hsvt2)
  299. title('H-inf LSS Design -- W3 & Comp. Sens. Func.')
  300. xlabel('Frequency - Rad/Sec')
  301. ylabel('SV - db')
  302. grid
  303. text(2,20,'H-inf (Gam = 1) ---> H-inf (Gam = 1.5)')
  304. pause
  305. axis
  306. semilogx(w,hsvtt1(1,:),w,hsvtt2(1,:))
  307. title('H-inf LSS Design -- Cost function Ty1u1')
  308. xlabel('Frequency - Rad/Sec')
  309. ylabel('SV - db')
  310. grid
  311. text(0.001,-1.5,'H-inf (Gam = 1) ---> H-inf (Gam = 1.5)')
  312. pause
  313. clc
  314. disp('  ')
  315. disp('  ')
  316. disp('               << H-inf Controller (Gam = 1.5) >>')
  317. disp('    ')
  318. disp('    Poles of Controller :')
  319. polecp = eig(acp)
  320. disp(' ')
  321. disp('                                 (strike a key to continue ...)')
  322. pause
  323. clc
  324. disp('  ')
  325. disp('    State-Space of the final H-inf Controller:')
  326. disp('   ')
  327. acp
  328. disp('                                 (strike a key to continue ...)')
  329. pause
  330. clc
  331. bcp
  332. disp('                                 (strike a key to continue ...)')
  333. pause
  334. clc
  335. ccp
  336. dcp   
  337. disp('                                 (strike a key to continue ...)')
  338. pause
  339. clc
  340. disp('   ')
  341. disp('    Poles of closed-loop TF matrix Ty1u1:')
  342. poletyu = eig(acl)
  343. disp('                                 (strike a key to continue ...)')
  344. pause
  345. clc
  346. disp(' ')
  347. disp(' ')
  348. disp('  ')
  349. disp('  ')
  350. disp('                   * * * * * * * * * * * * * * * * *')
  351. disp('                   *  End of LINFDEMO  ......      *')
  352. disp('                   *                               *')
  353. disp('                   * Good luck with your design !! *')
  354. disp('                   * * * * * * * * * * * * * * * * *')
  355. %
  356. % ----- End of LINFDEMO.M --- RYC/MGS %
  357.  
  358.  
  359.  
  360.  
  361.