home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / BALDEMO.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  9.0 KB  |  254 lines
  1. %
  2. % BALDEMO Balanced series model reductin techniques demonstration
  3.  
  4. % R. Y. Chiang & M. G. Safonov 3/88
  5. % Copyright (c) 1988 by the MathWorks, Inc.
  6. % All Rights Reserved.
  7.  
  8. clc
  9. disp('    ')
  10. disp(' ')
  11. disp('           << Demo #1: Balanced Model Reduction Techniques >>')
  12. disp('  ')
  13. disp(' ')
  14. disp('            This demo will show:')
  15. disp(' ')
  16. disp('             1). Schur Balanced Model Reduction Technique')
  17. disp('     ')
  18. disp('             2). Square-root Balanced Model Reduction Technique ')
  19. disp('  ')
  20. disp(' ')
  21. disp('  ')
  22. disp('  ')
  23. disp(' ')
  24. disp(' ')
  25. disp('                                 (strike a key to continue ...)')
  26. pause
  27. clc
  28. disp(' ')
  29. disp(' ')
  30. disp('                   << Introduction of Balanced Methods >>')
  31. disp('  ')
  32. disp('  ')
  33. disp('    The calculation of the balancing transformation required in Moore`s')
  34. disp(' ')
  35. disp('    "balanced-transformation" model reduction procedure tends to be badly')
  36. disp(' ')
  37. disp('    conditioned, especially for non-minimal models that stands to benefit')
  38. disp(' ')
  39. disp('    the most from model reduction. The Schur procedure described here can')
  40. disp(' ')
  41. disp('    directly compute the Moore`s reduced model without balancing at all,')
  42. disp(' ')
  43. disp('    hence completely circumvents the original ill-conditioned approach.')
  44. disp(' ')
  45. disp('    The square-root method for some problem can be ill-conditioned as well.')
  46. disp('  ')
  47. disp('  ')
  48. disp('                                 (strike a key to continue ...)')
  49. pause
  50. clc
  51. disp(' ')
  52. disp(' ')
  53. disp(' State-space of the 10-state model to be reduced (stable & nonminimal):')
  54. a =[-6    -1     0     0     0     0     0     0     0     0;
  55.      1    -8     0     0     0     0     0     0     0     0;
  56.      0     0   -10     3     0     0     0     0     0     0;
  57.      0     0     1    -8     0     0     0     0     0     0;
  58.      0     0     0     0   -13    -3     9     0     0     0;
  59.      0     0     0     0     1    -8     0     0     0     0;
  60.      0     0     0     0     0     1    -8     0     0     0;
  61.      0     0     0     0     0     0     0   -14    -9     0;
  62.      0     0     0     0     0     0     0     1    -8     0;
  63.      0     0     0     0     0     0     0     0     0    -2];
  64. b =[1         0;
  65.     0         0;
  66.     0         1;
  67.     0         0;
  68.     1         0;
  69.     0         0;
  70.     0         0;
  71.     0         1;
  72.     0         0;
  73.     0.001     0.001];
  74. c =[0 1 0 1 0 0 0 0 0 500000;
  75.     0 0 0 0 0 0 -6 1 -2 500000];
  76. d = [0 0;0 0];
  77. a
  78. disp('                                 (strike a key to continue ...)')
  79. pause
  80. clc
  81. b
  82. disp('                                 (strike a key to continue ...)')
  83. pause
  84. clc
  85. c
  86. disp('                                 (strike a key to continue ...)')
  87. pause
  88. clc
  89. d
  90. disp('                                 (strike a key to continue ...)')
  91. pause
  92. format short e
  93. clc
  94. disp('               << Technique 1. Schur Balanced Model Reduction >>')
  95. disp('      -----------------------------------------------------------------')
  96. disp('       First let us compute the Hankel singular values of the model...')
  97. disp('  ')
  98. disp('       [hsv] = hksv(a,b,c);    % Computing Hankel singular values')
  99. disp('      -----------------------------------------------------------------')
  100. [hsv] = hksv(a,b,c);
  101. hsv = hsv';
  102. format compact
  103. hsv
  104. format loose
  105. disp(' ')
  106. disp('      ----------------------------------------------------------------')
  107. disp('       Obviously, the regular "balanced transformation" fails immedi-')
  108. disp('       ately due to the last few Hankel singular vlaues being "0" ...')
  109. disp('      ----------------------------------------------------------------')
  110. disp(' ')
  111. disp('        (strike a key to see a bar chart of the Hankel Singular Values ...)')
  112. pause
  113. bar(hsv)
  114. title('Hankel Singular Values')
  115. pause
  116. clc
  117. disp(' ')
  118. disp(' ')
  119. disp('      ---------------------------------------------------------------')
  120. disp('       If we truncate the model to 4-state, the infinity-norm error ') 
  121. disp('       bound can be computed beforehand ...')
  122. disp(' ')
  123. disp('       bndsch = 2*sum(hsv(1,5:10));  % Infinity-norm error bound')
  124. disp('      ---------------------------------------------------------------')
  125. bndsch = 2*sum(hsv(1,5:10))
  126. disp('                                 (strike a key to continue ...)')
  127. pause
  128. clc
  129. disp(' ')
  130. disp('  ')
  131. disp('     ---------------------------------------------------------------------')
  132. disp('      Let us start Schur Balanced Trunction (Schur-BT) to get a 4-state')
  133. disp('      reduced model ...')
  134. disp('  ')
  135. disp('      [as,bs,cs,ds,agsh,hsv,tlh,trh] = schbal(a,b,c,d,1,4); % Run SCHBAL.M')
  136. disp('     ---------------------------------------------------------------------')
  137. [as,bs,cs,ds,augsh,hsv,tlh,trh] = schbal(a,b,c,d,1,4);
  138. bndsch = augsh(1,2);
  139. disp('                                               ^')
  140. disp('        - - - Computing SV Bode plot of G(s) & G(s)- - - ')
  141. w = logspace(-4,4,50);
  142. svg = sigma(a,b,c,d,1,w); svg = 20*log10(svg);
  143. svsch = sigma(as,bs,cs,ds,1,w); svsch = 20*log10(svsch);
  144. [aser,bser,cser,dser] = addss(a,b,c,d,as,bs,-cs,-ds);
  145. disp(' ')
  146. disp('        - - - Computing SV Bode plot of the error system - - -')
  147. sversch = sigma(aser,bser,cser,dser,1,w); sversch = 20*log10(sversch);
  148. bndsch = 20*log10(bndsch)*ones(w);
  149. disp(' ')
  150. disp('                                                               ^')
  151. disp('                       (strike a key to see the plot of G(s) & G(s)...)')
  152. pause
  153. clg
  154. semilogx(w,svg,w,svsch)
  155. title('Schur Balanced Truncation Technique')
  156. xlabel('Frequency - Rad/Sec')
  157. ylabel('SV - db')
  158. text(0.001,0,'ORIGINAL & REDUCED MODEL (SAME SV PLOTS)')
  159. grid
  160. pause
  161. disp(' ')
  162. disp('                                                    ^')
  163. disp('            (strike a key to see the plot of G(s) - G(s) & its bound..)')
  164. pause
  165. clg
  166. semilogx(w,bndsch,w,sversch)
  167. title('Schur-BT Technique (Infinity-norm error bound)')
  168. xlabel('Frequency - Rad/Sec')
  169. ylabel('SV - db')
  170. grid
  171. text(100,-50,'Error Bound')
  172. text(1,-80,'Error')
  173. pause
  174. clc
  175. disp(' ')
  176. disp('  ')
  177. disp('                 << Technique 2. Square-Root Balanced Method >>')
  178. disp('     --------------------------------------------------------------------')
  179. disp('       Let us try Square-Root method and compare with Schur-BT method to')
  180. disp('       get a 4-state model...')
  181. disp(' ')
  182. disp('       [aq,bq,cq,dq,agq,hsv,tlq,trq] = balsq(a,b,c,d,1,4); % Run BALSQ.M')
  183. disp('     --------------------------------------------------------------------')
  184. [aq,bq,cq,dq,agq,hsv,tlq,trq] = balsq(a,b,c,d,1,4);
  185. disp('  ')
  186. disp('  ')
  187. disp('                                 (strike a key to continue ...)')
  188. pause
  189. clc
  190. disp(' ')
  191. disp(' ')
  192. disp('      The singular values Bode plots are theoretically identical to ')
  193. disp(' ')
  194. disp('      the Schur-BT technique, therefore we omit them, but the condition')
  195. disp('  ')
  196. disp('      numbers of the final transformations differ a lot ......')
  197. disp('  ')
  198. disp('  ')
  199. disp('                                 (strike a key to continue ...)')
  200. pause
  201. clc
  202. disp('   ')
  203. disp('  ')
  204. disp('            << Condition Number of Shcur-BT vs. Square-Root BT >>')
  205. disp(' ')
  206. disp('                    Schur-BT               Square-Root BT')
  207. disp('               -----------------         ------------------')
  208. disp('                Slbig     Srbig           Slbig      Srbig')
  209. disp('  ')
  210. disp('                 6.1       6.1            27216      4049.5')
  211. disp(' ')
  212. disp(' ')
  213. disp('      A smaller condition number means less sensitivity to numerical roundoff')
  214. disp(' ')
  215. disp('      errors, i.e., SCHMR is more numerically robust than BALMR.')
  216. disp(' ')
  217. disp('                                 (strike a key to continue ...)')
  218. pause
  219. clc
  220. disp(' ')
  221. disp(' ')
  222. disp('  << Comparison of Schur and Square-Root with Finite Precision Arithmetic >>')
  223. disp(' ---------------------------------------------------------------------------')
  224. disp('   Let`s truncate the transformation matrices SLBIG & SRBIG of the two')
  225. disp('   methods to 1/1000 th finite precision, then compare the results.')
  226. disp(' ')
  227. disp('   It shows that round-off error creates a "terrible" result for the ')
  228. disp('   Square-Root method as one would expect from its ill-conditioned ')
  229. disp('   transformation matrices !!')
  230. disp(' ---------------------------------------------------------------------------')
  231. tlhrd = round(tlh*1000)/1000;
  232. trhrd = round(trh*1000)/1000;
  233. tlqrd = round(tlq*1000)/1000;
  234. trqrd = round(trq*1000)/1000;
  235. ass = tlhrd'*a*trhrd; bss = tlhrd'*b; css = c*trhrd; dss = ds;
  236. aqq = tlqrd'*a*trqrd; bqq = tlqrd'*b; cqq = c*trqrd; dqq = dq;
  237. svss = sigma(ass,bss,css,dss,1,w); svss = 20*log10(svss);
  238. svqq = sigma(aqq,bqq,cqq,dqq,1,w); svqq = 20*log10(svqq);
  239. disp(' ')
  240. disp(' ')
  241. disp(' ')
  242. disp('                         (strike a key to see the comparison ..)')
  243. pause
  244. semilogx(w,svsch,w,svss,w,svqq)
  245. title('Finite Precision Arithmetic (Schur vs. Square-Root)')
  246. xlabel('Frequency - Rad/Sec')
  247. ylabel('SV - db')
  248. grid
  249. text(0.001,-60,'- & -- : Schur method')
  250. text(0.001,-70,'.. & .-: rounded Schur')
  251. text(0.001,-80,'- & -- : rounded Square-Root')
  252. pause
  253. %
  254. % ----- End of BALDEMO.M --- RYC/MGS %