home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / BHRDEMO.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  5.9 KB  |  205 lines
  1. %
  2. % BHRDEMO Compare model reduction techniques.
  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 #4: Comparison of Basis-Free Model Reduction Techniques >>')
  12. disp('  ')
  13. disp(' ')
  14. disp(' ')
  15. disp('        ')
  16. disp('              Technique 1). Schur Balanced Model Reduction')
  17. disp('     ')
  18. disp('              Technique 2). Optimal Descriptor Hankel MDA')
  19. disp(' ')
  20. disp('              Technique 3). Schur BST-REM')
  21. disp('  ')
  22. disp(' ')
  23. disp(' ')
  24. disp('                                 (strike a key to continue ...)')
  25. pause
  26. clc
  27. disp('  ')
  28. disp(' ')
  29. disp(' ')
  30. disp('   We will compare these three "basis-free" model reduction techniques')
  31. disp(' ')
  32. disp('   via a SISO example. This nonminimal model has 7 states and one of the')
  33. disp('  ')
  34. disp('   mode is low-damped (~ -45 db notch). After we truncate to 3-state model,')
  35. disp(' ')
  36. disp('   only Schur BST-REM technique matches well with the 50 db notch and')
  37. disp(' ')
  38. disp('   its associated phase angle. All the other methods failed.')
  39. disp('  ')
  40. disp('  ')
  41. disp(' ')
  42. disp(' ')
  43. disp('                                 (strike a key to continue ...)')
  44. pause
  45. clc
  46. disp(' ')
  47. disp(' ')
  48. disp('Transfer function of the 7-state model to be reduced (stable & nonminimal) :')
  49. format short e
  50. n1 = [
  51. 5.4121e-02   4.3366e+01   5.5430e+01   3.2436e+01   2.4399e+01   6.4144e+00];
  52. n2 = [   2.6329e+00   3.0046e-01];
  53. num = [n1 n2];
  54. d1 = [
  55. 1.0000e+00   1.2606e+01   5.3483e+01   9.0942e+01   7.1828e+01   2.7218e+01];
  56. d2 = [4.7546e+00   3.0046e-01];
  57. den = [d1 d2];
  58. disp(' ')
  59. disp(' ')
  60. disp(' ')
  61. disp(' ')
  62. disp(' ')
  63. disp('       0.05(s^7 + 801s^6 + 1024s^5 + 599s^4 + 451s^3 + 119s^2 + 49s + 5.55')
  64. disp(' (s) = -----------------------------------------------------------------------')
  65. disp('       s^7 + 12.6s^6 + 53.48s^5 + 90.94s^4 + 71.83s^3 + 27.22s^2 + 4.75s + 0.3')
  66. disp(' ')
  67. disp(' ')
  68. disp(' ')
  69. disp(' ss_ = mksys(a,b,c,d); % Turn it into new system matrix!')
  70. disp(' ')
  71. disp(' ')
  72. disp('                                 (strike a key to continue ...)')
  73. pause
  74. clc
  75. disp(' ')
  76. disp(' ')
  77. disp('   - - - Computing the Bode plot of the 7-state model - - -')
  78. w = logspace(-3,4,200);
  79. [ga,ph] = bode(num,den,w);
  80. [a,b,c,d] = tf2ss(num,den);
  81. ss_ = mksys(a,b,c,d);
  82. ga = 20*log10(ga);
  83. disp(' ')
  84. disp('  ')
  85. disp('                       (strike a key to see the Bode plot of G(s)...)')
  86. pause
  87. clg
  88. semilogx(w,ga)
  89. title('Original 7-State Model')
  90. xlabel('Frequency - Rad/Sec')
  91. ylabel('Gain - db')
  92. grid
  93. pause
  94. semilogx(w,ph)
  95. title('Original 7-State Model')
  96. xlabel('Frequency - Rad/Sec')
  97. ylabel('Phase - deg')
  98. grid
  99. pause
  100. clc
  101. disp('       -----------------------------------------------------------------')
  102. disp('        [hsv] = hksv(a,b,c);    % Computing Hankel singular values')
  103. disp('       -----------------------------------------------------------------')
  104. [hsv] = hksv(a,b,c);
  105. hsv = hsv';
  106. format compact
  107. hsv
  108. format loose
  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('    First let us try Schur Balanced model reduction to get a 3-state')
  121. disp('    reduced model ...')
  122. disp(' ')
  123. disp('    [ss_s,bnds,hsvs] = schmr(ss_,1,3);      % Run SCHMR.M')
  124. disp('   ------------------------------------------------------------------')
  125. [ss_s,bnds,hsvs] = schmr(ss_,1,3);
  126. [as,bs,cs,ds] = branch(ss_s);
  127. disp(' ')
  128. disp(' ')
  129. disp('                                 (strike a key to continue ...)')
  130. pause
  131. clc
  132. disp(' ')
  133. disp('  ')
  134. disp('   ------------------------------------------------------------------')
  135. disp('    Let us try optimal descriptor Hankel model reduction to get a')
  136. disp('    3-state reduced model ...')
  137. disp(' ')
  138. disp('    [ss_h,bndh,hsvh] = ohklmr(ss_,1,3);    % Run OHKLMR.M')
  139. disp('   ------------------------------------------------------------------')
  140. [ss_h,bndh,hsvh] = ohklmr(ss_,1,3);   
  141. [ah,bh,ch,dh] = branch(ss_h);
  142. disp(' ')
  143. disp(' ')
  144. disp('                                 (strike a key to continue ...)')
  145. pause
  146. clc
  147. disp(' ')
  148. disp('  ')
  149. disp('   ------------------------------------------------------------------')
  150. disp('    Let us try Schur BST-REM model reduction to get a 3-state model..')
  151. disp(' ')
  152. disp('    [ss_r,augr,hsvr] = bstschmr(ss_,1,3);  % Run BSTSCHMR.M')
  153. disp('   ------------------------------------------------------------------')
  154. [ss_r,augr,hsvr] = bstschmr(ss_,1,3);
  155. [ar,br,cr,dr] = branch(ss_r);
  156. disp(' ')
  157. hsvr = hsvr';
  158. format compact
  159. hsvr
  160. format loose
  161. disp(' ')
  162. disp('        (strike a key to see a bar chart of the Phase Matrix Hankel SV`s ...)')
  163. pause
  164. bar(hsvr)
  165. title('Hankel Singular Values of the Phase Matrix')
  166. pause
  167. clc
  168. disp(' ')
  169. disp(' ')
  170. disp('         - - - Computing the Bode Plots - - - ')
  171. [gs,phs] = bode(as,bs,cs,ds,1,w); gs = 20*log10(gs);
  172. [gh,phh] = bode(ah,bh,ch,dh,1,w); gh = 20*log10(gh);
  173. [gr,phr] = bode(ar,br,cr,dr,1,w); gr = 20*log10(gr);
  174. disp(' ')
  175. disp(' ')
  176. disp(' ')
  177. disp(' ')
  178. disp(' ')
  179. disp('                         (strike a key to see the Bode plots ..)')
  180. pause
  181. clg
  182. disp(' ')
  183. semilogx(w,ga,w,gs,w,gh,w,gr)
  184. title('Schur BST-REM vs. Schur BT and Des. Hankel (7-state --> 3-state)')
  185. xlabel('Frequency - Rad/Sec')
  186. ylabel('Gain - db')
  187. text(1,-30,'solid : original model')
  188. text(1,-33,'----- : Schur-BT')
  189. text(1,-36,'..... : Optimal Des. Hankel MDA')
  190. text(1,-39,'-.-.-.: Schur BST-REM')
  191. grid
  192. pause
  193. semilogx(w,ph,w,phs,w,phh,w,phr)
  194. title('Schur BST-REM vs. Schur BT and Des. Hankel (7-state --> 3-state)')
  195. xlabel('Frequency - Rad/Sec')
  196. ylabel('Phase - deg')
  197. text(0.002,-300,'solid : original model')
  198. text(0.002,-330,'----- : Schur-BT')
  199. text(0.002,-360,'..... : Optimal Des. Hankel MDA')
  200. text(0.002,-390,'-.-.-.: Schur BST-REM')
  201. grid
  202. pause
  203. %
  204. % ----- End of BHRDEMO.M --- RYC/MGS %
  205.