home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / ACTDEMO.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  10.5 KB  |  313 lines
  1. %
  2. % ACTDEMO Digital Hydraulic Actuator 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('          << Demo #1: SISO Digital Hydraulic Actuator Design Example >>')
  14. disp('  ')
  15. disp(' ')
  16. disp('  ')
  17. disp('                         -------             -----        -----')
  18. disp('        +        /      | H-Inf |     /     |     |      | Hyd.|') 
  19. disp('      --->(X)----  -----| Cont- | ----  --->| ZOH |----->|     |--------->')
  20. disp('           ^ - Ts:0.01  | roller|   Ts:0.01 |     |      | ACT.|    |')
  21. disp('           |      sec    -------       sec   -----        -----     |')
  22. disp('           |                                                        |')
  23. disp('           |________________________________________________________|')
  24. disp('  ')
  25. disp('  ')
  26. disp('                                 (strike a key to continue ...)')
  27. pause
  28. format short e
  29. num = 9000;
  30. den = [1 30 700 1000];
  31. [a,b,c,d] = tf2ss(num,den);
  32. clc
  33. disp('  ')
  34. disp('  ')
  35. disp('   Hydraulic Actuator Open Loop:')
  36. disp(' ')
  37. disp('   ')
  38. disp('                                 9000')
  39. disp('             G(s) = -----------------------------')
  40. disp('                     s^3 + 30 s^2 + 700 s + 1000')
  41. disp('  ')
  42. disp(' ')
  43. disp('      Poles of the open loop plant:')
  44. poleg = roots(den)
  45. disp('  ')
  46. disp('  ')
  47. disp('  ')
  48. disp('                                 (strike a key to continue ...)')
  49. pause
  50. clc
  51. disp(' ')
  52. disp(' ')
  53. disp('                              << W-Plane Design >>')
  54. disp('  ')
  55. disp('  ')
  56. disp('           Let us attach the plant with a Z.O.H and convert it into ')
  57. disp(' ')
  58. disp('           W-plane (sampling period: 0.01 sec) ....')
  59. disp(' ')
  60. disp('  ')
  61. disp('         ------------------------------------------------------------')
  62. disp('           [az,bz] = c2d(a,b,0.01);   % Attach the plant with Z.O.H')
  63. disp('         ------------------------------------------------------------')
  64. [az,bz] = c2d(a,b,0.01);
  65. disp('           [ag,bg,cg,dg] = bilin(az,bz,c,d,-1,`Tustin`,0.01);')
  66. disp('                           % Convert the plant with Z.O.H to W-plane')
  67. disp('         ------------------------------------------------------------')
  68. [ag,bg,cg,dg] = bilin(az,bz,c,d,-1,'Tustin',0.01);
  69. disp('  ')
  70. disp('  ')
  71. disp('                                 (strike a key to continue ...)')
  72. pause
  73. clc
  74. disp(' ')
  75. disp(' ')
  76. disp('     - - - Computing Bode plot of the open loop plant (in s & w) - - -')
  77. w = logspace(-3,5,50);
  78. svg = bode(a,b,c,d,1,w); svg = 20*log10(svg);
  79. svw = bode(ag,bg,cg,dg,1,w); svw = 20*log10(svw);
  80. disp(' ')
  81. disp(' ')
  82. disp(' ')
  83. disp(' ')
  84. disp('                             (strike a key to see the plot ...)')
  85. pause
  86. clg
  87. plot(111)
  88. semilogx(w,svg,w,svw)
  89. title('SISO Hydraulic Actuator Open Loop (in s & w domain)')
  90. xlabel('Frequency - Rad/Sec')
  91. ylabel('SV - db')
  92. text(150,-30,'Nyquist Freq.: 100pi')
  93. grid
  94. pause
  95. clc
  96. disp('  ')
  97. disp('                       << Design Specifications >>  ')
  98. disp(' ')
  99. disp('      1). Robustness Spec. : closed loop bandwidth -- 30 r/s')
  100. disp('          Associated Weighting:')
  101. disp('     ')
  102. disp('                    -1     3.16(1 + s/300)')
  103. disp('                  W3(s) = -----------------')
  104. disp('                             (1 + s/10)')
  105. disp('   ')
  106. disp(' ')
  107. disp('      2). Performance Spec.: sensitivity reduction of at least 100:1')
  108. disp('                             up to approx. 1 r/s')
  109. disp('          Associated Weighting:')
  110. disp(' ')          
  111. disp('                   -1         -1   0.01(1 + s)^2')
  112. disp('                  W1(s) = Gam  * ---------------')
  113. disp('                                   (1 + s/30)^2')
  114. disp('   ')
  115. disp('          where "Gam" in our design goes from 1 --> 1.5')
  116. nuw3i = 3.16*[1/300 1]; dnw3i = [1/10 1];
  117. svw3i = bode(nuw3i,dnw3i,w); svw3i = 20*log10(svw3i);
  118. nuw1i = 0.01*[1 2 1]; dnw1i = conv([1/30 1],[1/30 1]);
  119. svw1i = bode(nuw1i,dnw1i,w); svw1i = 20*log10(svw1i);
  120. disp('  ')
  121. disp('  ')
  122. disp('                   (strike a key to see the plot of the weightings ...)')
  123. pause
  124. clg
  125. plot(111)
  126. semilogx(w,svw1i,w,svw3i)
  127. grid
  128. title('MIMO LSS Design Example -- Design Specifications')
  129. xlabel('Frequency - Rad/Sec')
  130. ylabel('1/W1 & 1/W3 - db')
  131. text(0.005,-20,'Sensitivity Spec.-- 1/W1(s)')
  132. text(100,0,'Robustness Spec.')
  133. text(1000,-10,'1/W3(s)')
  134. pause
  135. clc
  136. disp('                      << Problem Formulation >>')
  137. disp(' ')
  138. disp('      Form an augmented plant P(s) with these two weighting functions:')
  139. disp(' ')
  140. disp('                 1). W1 penalizing error signal "e"')
  141. disp('  ')
  142. disp('                 2). W3 penalizing plant output "y"')
  143. disp(' ')
  144. disp('      and find a stabilizing controller F(s) such that the Hinf-norm')
  145. disp('      of TF Ty1u1 is minimized and less than one, i.e.')
  146. disp('  ')
  147. disp('                         min |Ty1u1|   < 1,')
  148. disp('                         F(s)       inf')
  149. disp(' ')
  150. disp('      where ')
  151. disp('                       |               -1|')                         
  152. disp('               Ty1u1 = | Gam*W1*(I + GF) |  = | Gam * W1 * S  |')
  153. disp('                       |               -1|    |  W3 * (I - S) |')
  154. disp('                       |  W3*GF*(I + GF) |')
  155. disp('  ')
  156. disp('  ')
  157. disp('  ')
  158. disp('                                 (strike a key to continue ...)')
  159. pause
  160. clc
  161. disp('  ')
  162. disp(' ')
  163. disp('                              << DESIGN PROCEDURE >>')
  164. disp('  ')
  165. disp('            * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *')
  166. disp('            *    [Step 1]. Do plant augmentation (run AUGTF.M or      *')
  167. disp('            *              AUGSS.M)                                   *')
  168. disp('            *                                                         *')
  169. disp('            *    [Step 2]. Do H-Inf synthesis (run HINF.M)            *')
  170. disp('            *                                                         *')
  171. disp('            *    [Step 3]. Redo the plant augmentation by setting     *')
  172. disp('            *              new "Gam" --> 1.5 and rerun HINF.M         *')
  173. disp('            * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *')
  174. disp('  ')
  175. disp('  ')
  176. disp('                                 (strike a key to continue ...)')
  177. pause
  178. clc
  179. disp('  ')
  180. disp('           Assign the cost coefficient "Gam" --> 1 ')
  181. disp('       ')
  182. disp('           this will serve as the baseline design ....')
  183. disp('  ')
  184. disp('    ---------------------------------------------------------------')
  185. disp('     % Plant augmentation of the actuator:')
  186. disp('       w1 = [Gam*dnw1i;nuw1i];')
  187. disp('       w3 = [dnw3i;nuw3i];')
  188. disp('       w2 = [];')
  189. disp('       ss_g  = mksys(ag,bg,cg,dg);')
  190. disp('       [TSS_]=augtf(ss_g,w1,w2,w3);')
  191. disp('    ---------------------------------------------------------------')
  192. Gam = input('           Input cost coefficient "Gam" = ');
  193. w1 = [Gam*dnw1i;nuw1i];
  194. w3 = [dnw3i;nuw3i];
  195. w2 = [];
  196. ss_g = mksys(ag,bg,cg,dg);
  197. [TSS_] = augtf(ss_g,w1,w2,w3);
  198. disp('  ')
  199. disp('     - - - State-Space TSS_:=(A,B1,..,D21,D22) is ready for')
  200. disp('           the Small-Gain problem - - -')
  201. disp('  ')
  202. disp('  ')
  203. disp('                                 (strike a key to continue ...)')
  204. pause
  205. clc
  206. disp(' ')
  207. disp(' ')
  208. disp('   --------------------------------------------------------------- ')
  209. disp('    [ss_cp,ss_cl,hinfo] = hinf(TSS_);')
  210. disp('    [acp,bcp,ccp,dcp] = branch(ss_cp);')
  211. disp('    [acl,bcl,ccl,dcl] = branch(ss_cl);')
  212. disp('                    % Running file HINF.M for H-Inf optimization')
  213. disp('   ---------------------------------------------------------------')
  214. [ss_cp,ss_cl,hinfo] = hinf(TSS_);
  215. [acp,bcp,ccp,dcp] = branch(ss_cp);
  216. [acl,bcl,ccl,dcl] = branch(ss_cl);
  217. disp('  ')
  218. disp('  ')
  219. disp('                                 (strike a key to continue ...)')
  220. pause
  221. pltopt           % Preparing singular values for plotting
  222. svw1i1 = svw1i; hsvs1 = svs; hsvt1 = svt; hsvtt1 = svtt;
  223. disp('  ')
  224. disp('  ')
  225. disp('                                 (strike a key to continue ...)')
  226. pause
  227. clc      
  228. disp('  ')
  229. disp('  ')
  230. disp('  ')
  231. disp('     After a few iterations, we found a new Gam of 1.5 can push the')
  232. disp('  ')
  233. disp('     H-Inf cost function close to its limit. ')
  234. disp('  ')
  235. disp('  ')
  236. disp('            Input "Gam" --> 1.5, and try HINF again .....')
  237. disp('  ')
  238. disp('  ')
  239. disp('                                 (strike a key to continue ...)')
  240. pause
  241. Gam = input('           Input cost coefficient "Gam" = ');
  242. w1 = [Gam*dnw1i;nuw1i];
  243. [TSS_] = augtf(ss_g,w1,w2,w3);
  244. disp('  ')
  245. disp('  ')
  246. disp('                                 (strike a key to continue ...)')
  247. pause
  248. [ss_cp,ss_cl,hinfo] = hinf(TSS_);
  249. [acp,bcp,ccp,dcp] = branch(ss_cp);
  250. [acl,bcl,ccl,dcl] = branch(ss_cl);
  251. disp('  ')
  252. disp('  ')
  253. disp('                                 (strike a key to continue ...)')
  254. pause
  255. pltopt
  256. svw1i2 = svw1i; hsvs2 = svs; hsvt2 = svt; hsvtt2 = svtt;
  257. disp('  ')
  258. disp('  ')
  259. disp('                  (strike a key to see the plots of the comparison ...)')
  260. pause
  261. clg
  262. plot(111)
  263. semilogx(w,svw1i1,w,hsvs1,w,svw1i2,w,hsvs2)
  264. title('H-Inf W-Plane Actuator Design -- 1/W1 & Sensitivity Func.')
  265. xlabel('Frequency - Rad/Sec')
  266. ylabel('SV - db')
  267. grid
  268. text(0.002,10,'H-Inf (Gam = 1) ---> H-Inf (Gam = 1.5)')
  269. pause
  270. semilogx(w,svw3i,w,hsvt1,w,hsvt2)
  271. title('H-Inf W-Plane Actuator Design -- 1/W3 & Comp. Sens. Func.')
  272. xlabel('Frequency - Rad/Sec')
  273. ylabel('SV - db')
  274. grid
  275. text(0.002,-30,'H-Inf (Gam = 1) ---> H-Inf (Gam = 1.5)')
  276. pause
  277. semilogx(w,hsvtt1,w,hsvtt2)
  278. title('H-Inf W-Plane Actuator Design -- Cost function Ty1u1')
  279. xlabel('Frequency - Rad/Sec')
  280. ylabel('SV - db')
  281. grid
  282. text(0.002,-10,'H-Inf (Gam = 1) ---> H-Inf (Gam = 1.5)')
  283. pause
  284. clc
  285. disp('  ')
  286. disp('  ')
  287. disp('               << H-Inf Controller (Gam = 1.5) >>')
  288. disp('    ')
  289. disp('    Poles of Controller :')
  290. polecp = eig(acp)
  291. disp(' ')
  292. disp('                                 (strike a key to continue ...)')
  293. pause
  294. disp(' ')
  295. disp('    Transmission Zeros of the Controller :')
  296. tzcp = tzero(acp,bcp,ccp,dcp)
  297. disp(' ')
  298. disp('                                 (strike a key to continue ...)')
  299. pause
  300. svcp = sigma(acp,bcp,ccp,dcp,1,w); svcp = 20*log10(svcp);
  301. semilogx(w,svcp)
  302. title('Final 8-State H-Inf Controller')
  303. xlabel('Rad/Sec')
  304. ylabel('SV (db)')
  305. pause
  306. clc
  307. disp('    Poles of the Cost Function :')
  308. poletyu = eig(acl)
  309. disp('                                 (strike a key to continue ...)')
  310. pause
  311. %
  312. % ----- End of ACTDEMO.M --- RYC/MGS %
  313.