home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / REMDEMO.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  9.1 KB  |  229 lines
  1. %
  2. % REMDEMO Balanced Stochastic Truncation and Relative Error Model Reduction
  3. %         via Schur method (Schur BST-REM).
  4. %
  5.  
  6. % R. Y. Chiang & M. G. Safonov 3/88
  7. % Copyright (c) 1988 by the MathWorks, Inc.
  8. % All Rights Reserved.
  9. % -------------------------------------------------------------------------
  10. clc
  11. disp('    ')
  12. disp(' ')
  13. disp('            << Demo #3: Schur BST-REM Model Reduction Technique >>')
  14. disp('  ')
  15. disp(' ')
  16. disp(' ')
  17. disp(' ')
  18. disp('          Example: 16-State Nonminimal MIMO Model reduced to 6-State')
  19. disp('     ')
  20. disp('  ')
  21. disp(' ')
  22. disp('  ')
  23. disp('  ')
  24. disp(' ')
  25. disp(' ')
  26. disp('                                 (strike a key to continue ...)')
  27. pause
  28. clc
  29. disp(' ')
  30. disp(' ')
  31. disp('                   << Introduction of Schur BST-REM >>')
  32. disp('  ')
  33. disp('  ')
  34. disp('    A Schur Relative Error Method (REM) for state-space model order reduction')
  35. disp(' ')
  36. disp('    is described. Our algorithm is based on Desai`s Balanced Stochastic')
  37. disp(' ')
  38. disp('    Truncation (BST) technique for which Wang and Safonov have obtained')
  39. disp(' ')
  40. disp('    an L-inf norm relative-error bound. However, unlike previous methods,')
  41. disp(' ')
  42. disp('    our Schur method completely circumvents the numerically delicate initial')
  43. disp(' ')
  44. disp('    step of obtaining a minimal balanced stochastic realization (BSR) of the')
  45. disp(' ')
  46. disp('    power spectrum matrix G(s)G`(-s).')
  47. disp(' ')
  48. disp('  ')
  49. disp('  ')
  50. disp('                                 (strike a key to continue ...)')
  51. pause
  52. clc
  53. format short e
  54. disp(' State-space of the 16-state model to be reduced (stable & nonminimal):')
  55. a1 = [
  56.   -2.0000e+03  -1.1493e+00   1.4858e+03  -6.0345e+02   7.7992e+02   4.5933e+02;
  57.    1.3549e-13  -2.0981e-02  -2.3740e+00   1.2214e+00  -2.4659e+00  -1.1860e+00;
  58.    3.0579e-14  -1.2737e-14  -5.6757e+00   2.2035e+00   2.0629e+01   1.0031e+01;
  59.   -1.5913e-14   5.1825e-15   8.1765e-06  -2.5785e-01   3.6547e+00   3.3445e+00;
  60.    7.7097e-15  -6.4280e-15  -1.0629e-05   7.4497e-05  -3.0000e+01  -1.5516e-04;
  61.    4.0211e-15  -3.9663e-15  -5.2035e-06   3.8799e-05   5.3245e-06  -3.0000e+01;
  62.    1.0512e-26  -1.3191e-26   7.4402e-18  -1.8651e-18  -2.4512e-16  -2.4600e-16;
  63.    1.5017e-26  -6.5896e-27   3.7201e-18  -9.3252e-19  -1.2256e-16  -1.2300e-16;
  64.             0            0            0            0            0            0;
  65.             0            0            0            0            0            0;
  66.             0            0            0            0            0            0;
  67.             0            0            0            0            0            0;
  68.             0            0            0            0            0            0;
  69.             0            0            0            0            0            0;
  70.             0            0            0            0            0            0;
  71.             0            0            0            0            0            0];
  72. a2 = [
  73.    1.1539e-09   5.7695e-10            0            0            0            0;
  74.    7.5319e-13   3.7660e-13            0            0            0            0;
  75.   -4.0627e-12   2.4833e-12            0            0            0            0;
  76.   -2.2146e-12  -2.9492e-12            0            0            0            0;
  77.    1.0331e-11   7.4681e-12            0            0            0            0;
  78.    6.3843e-12   4.5653e-12            0            0            0            0;
  79.   -1.0000e-02   3.9925e-24            0            0            0            0;
  80.    2.9281e-25  -1.0000e-02            0            0            0            0;
  81.             0            0  -1.6184e+01   1.5083e+02   4.8353e+02  -3.7230e+02;
  82.             0            0  -2.1463e+01  -1.1303e+02  -3.9936e+02   4.5213e+02;
  83.             0            0   1.0421e+01   7.2968e+01   2.5251e+02  -2.6766e+02;
  84.             0            0   4.2729e+00   5.0112e+01   1.6975e+02  -1.6710e+02;
  85.             0            0   1.3789e+01   1.0265e+02   3.5676e+02  -3.7530e+02;
  86.             0            0  -5.9738e+01  -3.3114e+02  -1.1675e+03   1.2905e+03;
  87.             0            0   2.0960e+03   1.1873e+04   4.1764e+04  -4.5988e+04;
  88.             0            0   1.6822e+03  -1.6399e+04  -5.2599e+04   4.0665e+04];
  89. a3 = [
  90.             0            0            0            0;
  91.             0            0            0            0;
  92.             0            0            0            0;
  93.             0            0            0            0;
  94.             0            0            0            0;
  95.             0            0            0            0;
  96.             0            0            0            0;
  97.             0            0            0            0;
  98.    2.3930e-01  -3.3054e-01  -6.5690e-03   1.2279e+00;
  99.   -6.8475e+00   6.3611e+00   1.5087e+00   4.1866e-02;
  100.    2.6085e+00  -2.4529e+00  -5.6231e-01   3.1315e-01;
  101.    5.9270e-01  -5.8463e-01  -1.1822e-01   3.5004e-01;
  102.    2.9069e+00  -5.8875e+00  -9.0978e-01   8.2380e-01;
  103.   -1.1292e+01   1.7963e+01   2.2323e+00   5.8339e-01;
  104.    3.6532e+02  -5.8425e+02  -4.9193e+01   1.9730e-01;
  105.   -9.0859e+00   9.8929e+00   1.0133e+00  -3.2254e+01];
  106. a = [a1 a2 a3];
  107. b =[
  108.   -9.3933e+02   1.1326e+03;
  109.    6.3293e+01  -1.5312e+01;
  110.   -2.7062e-04   4.0846e-03;
  111.    2.6371e-03   1.5393e-03;
  112.   -5.5859e-04   8.9453e-04;
  113.   -3.6171e-04   5.5788e-04;
  114.    7.5916e-16  -1.3986e-15;
  115.    3.7954e-16  -6.9930e-16;
  116.    1.5533e-02  -1.5516e-01;
  117.    1.7881e-01   5.2796e-03;
  118.   -8.6917e-02  -1.8393e-02;
  119.   -3.9949e-02  -2.3137e-02;
  120.   -1.2324e-01  -2.8050e-02;
  121.    4.7783e-01   3.8904e-02;
  122.   -1.6720e+01  -1.5469e+00;
  123.   -1.5476e+00   1.6728e+01];
  124. c1 = [
  125.   -9.1517e-08  -5.2591e-11   3.1223e-03  -1.2700e-03  -1.9162e-05  -4.5190e-04;
  126.    9.8922e-07   5.6846e-10  -1.7184e-03   6.9980e-04   9.5218e-04   5.0433e-04];
  127. c2 = [
  128.    3.6519e-16  -2.0788e-17  -2.0891e+01  -1.1933e+02  -4.1720e+02   4.6024e+02;
  129.    9.2135e-16   1.6474e-15  -1.6124e+01   1.6319e+02   5.2691e+02  -4.0617e+02];
  130. c3 = [
  131.   -4.8842e+00   5.2742e+00   1.4817e+00   7.4520e-03;
  132.   -2.8106e+00  -1.8255e+00  -4.2506e-02   1.3483e+00];
  133. c = [c1 c2 c3];
  134. d = [
  135.    1.6729e-01   1.5477e-02
  136.    1.5477e-02  -1.6729e-01];
  137. disp('     First 6 column of the A matrix:')
  138. a1
  139. disp('                                 (strike a key to continue ...)')
  140. pause
  141. clc
  142. disp('     Column 7 to 12 of the A matrix:')
  143. a2
  144. disp('                                 (strike a key to continue ...)')
  145. pause
  146. clc
  147. disp('     Column 13 to 16 of the A matrix:')
  148. a3
  149. disp('                                 (strike a key to continue ...)')
  150. pause
  151. clc
  152. b
  153. disp('                                 (strike a key to continue ...)')
  154. pause
  155. clc
  156. c
  157. disp('                                 (strike a key to continue ...)')
  158. pause
  159. clc
  160. d
  161. disp('                                 (strike a key to continue ...)')
  162. pause
  163. clc
  164. disp('   ---------------------------------------------------------------------')
  165. disp('    Let us start Schur BST-REM model reduction to get a 6-state reduced')
  166. disp('    order model ...')
  167. disp(' ')
  168. disp('    [ar,br,cr,dr,aug,hsv] = bstschmr(a,b,c,d,1,6); % Run BSTSCHMR.M')
  169. disp('   ---------------------------------------------------------------------')
  170. [ar,br,cr,dr,aug,hsv] = bstschmr(a,b,c,d,1,6);
  171. hsv = hsv';
  172. format compact
  173. hsv
  174. format loose
  175. disp(' ')
  176. disp('      ----------------------------------------------------------------')
  177. disp('       Obviously, the regular "Balanced" Stochastic Trunciton/Relative')
  178. disp('       Error Method fails immediately due to the last Hankel singular ')
  179. disp('       vlaue of the "phase matrix" being "0" ...')
  180. disp('      ----------------------------------------------------------------')
  181. disp(' ')
  182. disp('        (strike a key to see a bar chart of the Hankel Singular Values ...)')
  183. pause
  184. bar(hsv)
  185. title('Hankel Singular Values of the Phase Matrix')
  186. pause
  187. clc
  188. disp(' ')
  189. disp(' ')
  190. disp('                                                 ^')
  191. disp('        - - - Computing SV Bode plot of G(s) and G(s) - - - ')
  192. w = logspace(-4,4,100);
  193. svg = sigma(a,b,c,d,1,w); svg = 20*log10(svg);
  194. svrem = sigma(ar,br,cr,dr,1,w); svrem = 20*log10(svrem); 
  195. disp(' ')
  196. disp('        - - - Computing SV Bode plot of the relative error - - -')
  197. [arer,brer,crer,drer] = addss(a,b,c,d,ar,br,-cr,-dr);
  198. [ai,bi,ci,di] = ssinv(a,b,c,d);
  199. [arer,brer,crer,drer] = series(arer,brer,crer,drer,ai,bi,ci,di);
  200. sverrem = sigma(arer,brer,crer,drer,1,w); sverrem = 20*log10(sverrem);
  201. remdb = 20*log10(aug(1,2))*ones(w);
  202. disp('  ')
  203. disp('                                                              ^')
  204. disp('                    (strike a key to see the plot of G(s) and G(s) ...)')
  205. pause
  206. clg
  207. semilogx(w,svg,w,svrem)
  208. title('Schur BST-REM Model Reduction (16-state --> 6-state)')
  209. xlabel('Frequency - Rad/Sec')
  210. ylabel('SV - db')
  211. text(1,-100,'Original 16-State')
  212. text(0.001,-60,'6-State Approximation')
  213. grid
  214. pause
  215. disp(' ')
  216. disp(' ')
  217. disp('  (strike a key to see the plot of the relative error and its bound ..)')
  218. pause
  219. clg
  220. semilogx(w,remdb,w,sverrem)
  221. title('Schur BST-REM Model Reduction (16-State ---> 6-State)')
  222. xlabel('Frequency - Rad/Sec')
  223. ylabel('SV - db')
  224. grid
  225. text(10,40,'Relative Error Bound')
  226. text(100,-20,'Relative Error')
  227. pause
  228. %
  229. % ----- End of REMDEMO.M --- RYC/MGS %