home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / OHKDEMO.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  12.2 KB  |  302 lines
  1. %
  2. % OHKDEMO Optimal Descriptor Hankel Model Reduction Technique & Anticausal
  3. %         Projection demonstration
  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 #2: Optimal Descriptor Hankel Approximation >>')
  14. disp('  ')
  15. disp(' ')
  16. disp('             This demo will show:')
  17. disp(' ')
  18. disp('               1). Optimal Descriptor Hankel Model Reduction')
  19. disp('                   (Optimal Hankel MDA)')
  20. disp('     ')
  21. disp('               2). Optimal Descriptor Anticausal Projection')
  22. disp('  ')
  23. disp(' ')
  24. disp('  ')
  25. disp('  ')
  26. disp(' ')
  27. disp(' ')
  28. disp('                                 (strike a key to continue ...)')
  29. pause
  30. clc
  31. disp(' ')
  32. disp(' ')
  33. disp('                   << Introduction of Hankel Methods >>')
  34. disp('  ')
  35. disp('  ')
  36. disp('    A basis-free descriptor system representation is shown to facilitate')
  37. disp(' ')
  38. disp('    the computation of all minimum-degree and optimal k-th order all-pass')
  39. disp(' ')
  40. disp('    extensions and Hankel-norm approximants. The descriptor representation')
  41. disp(' ')
  42. disp('    has the same simple form for both the optimal and suboptimal cases.')
  43. disp(' ')
  44. disp('    The method makes Hankel model reduction practical for non-minimal and')
  45. disp(' ')
  46. disp('    near non-minimal systems by eliminating the ill-conditioned calculation')
  47. disp(' ')
  48. disp('    of a minimal balanced realization. A simple, numerically sound method')
  49. disp(' ')
  50. disp('    based on singular-value demcompositon enables the results to be expressed')
  51. disp(' ')
  52. disp('    in state-space form.')
  53. disp('  ')
  54. disp('  ')
  55. disp('                                 (strike a key to continue ...)')
  56. pause
  57. clc
  58. format short e
  59. disp(' State-space of the 16-state model to be reduced (stable & nonminimal):')
  60. a1 = [
  61.   -2.0000e+03  -1.1493e+00   1.4858e+03  -6.0345e+02   7.7992e+02   4.5933e+02;
  62.    1.3549e-13  -2.0981e-02  -2.3740e+00   1.2214e+00  -2.4659e+00  -1.1860e+00;
  63.    3.0579e-14  -1.2737e-14  -5.6757e+00   2.2035e+00   2.0629e+01   1.0031e+01;
  64.   -1.5913e-14   5.1825e-15   8.1765e-06  -2.5785e-01   3.6547e+00   3.3445e+00;
  65.    7.7097e-15  -6.4280e-15  -1.0629e-05   7.4497e-05  -3.0000e+01  -1.5516e-04;
  66.    4.0211e-15  -3.9663e-15  -5.2035e-06   3.8799e-05   5.3245e-06  -3.0000e+01;
  67.    1.0512e-26  -1.3191e-26   7.4402e-18  -1.8651e-18  -2.4512e-16  -2.4600e-16;
  68.    1.5017e-26  -6.5896e-27   3.7201e-18  -9.3252e-19  -1.2256e-16  -1.2300e-16;
  69.             0            0            0            0            0            0;
  70.             0            0            0            0            0            0;
  71.             0            0            0            0            0            0;
  72.             0            0            0            0            0            0;
  73.             0            0            0            0            0            0;
  74.             0            0            0            0            0            0;
  75.             0            0            0            0            0            0;
  76.             0            0            0            0            0            0];
  77. a2 = [
  78.    1.1539e-09   5.7695e-10            0            0            0            0;
  79.    7.5319e-13   3.7660e-13            0            0            0            0;
  80.   -4.0627e-12   2.4833e-12            0            0            0            0;
  81.   -2.2146e-12  -2.9492e-12            0            0            0            0;
  82.    1.0331e-11   7.4681e-12            0            0            0            0;
  83.    6.3843e-12   4.5653e-12            0            0            0            0;
  84.   -1.0000e-02   3.9925e-24            0            0            0            0;
  85.    2.9281e-25  -1.0000e-02            0            0            0            0;
  86.             0            0  -1.6184e+01   1.5083e+02   4.8353e+02  -3.7230e+02;
  87.             0            0  -2.1463e+01  -1.1303e+02  -3.9936e+02   4.5213e+02;
  88.             0            0   1.0421e+01   7.2968e+01   2.5251e+02  -2.6766e+02;
  89.             0            0   4.2729e+00   5.0112e+01   1.6975e+02  -1.6710e+02;
  90.             0            0   1.3789e+01   1.0265e+02   3.5676e+02  -3.7530e+02;
  91.             0            0  -5.9738e+01  -3.3114e+02  -1.1675e+03   1.2905e+03;
  92.             0            0   2.0960e+03   1.1873e+04   4.1764e+04  -4.5988e+04;
  93.             0            0   1.6822e+03  -1.6399e+04  -5.2599e+04   4.0665e+04];
  94. a3 = [
  95.             0            0            0            0;
  96.             0            0            0            0;
  97.             0            0            0            0;
  98.             0            0            0            0;
  99.             0            0            0            0;
  100.             0            0            0            0;
  101.             0            0            0            0;
  102.             0            0            0            0;
  103.    2.3930e-01  -3.3054e-01  -6.5690e-03   1.2279e+00;
  104.   -6.8475e+00   6.3611e+00   1.5087e+00   4.1866e-02;
  105.    2.6085e+00  -2.4529e+00  -5.6231e-01   3.1315e-01;
  106.    5.9270e-01  -5.8463e-01  -1.1822e-01   3.5004e-01;
  107.    2.9069e+00  -5.8875e+00  -9.0978e-01   8.2380e-01;
  108.   -1.1292e+01   1.7963e+01   2.2323e+00   5.8339e-01;
  109.    3.6532e+02  -5.8425e+02  -4.9193e+01   1.9730e-01;
  110.   -9.0859e+00   9.8929e+00   1.0133e+00  -3.2254e+01];
  111. a = [a1 a2 a3];
  112. b =[
  113.   -9.3933e+02   1.1326e+03;
  114.    6.3293e+01  -1.5312e+01;
  115.   -2.7062e-04   4.0846e-03;
  116.    2.6371e-03   1.5393e-03;
  117.   -5.5859e-04   8.9453e-04;
  118.   -3.6171e-04   5.5788e-04;
  119.    7.5916e-16  -1.3986e-15;
  120.    3.7954e-16  -6.9930e-16;
  121.    1.5533e-02  -1.5516e-01;
  122.    1.7881e-01   5.2796e-03;
  123.   -8.6917e-02  -1.8393e-02;
  124.   -3.9949e-02  -2.3137e-02;
  125.   -1.2324e-01  -2.8050e-02;
  126.    4.7783e-01   3.8904e-02;
  127.   -1.6720e+01  -1.5469e+00;
  128.   -1.5476e+00   1.6728e+01];
  129. c1 = [
  130.   -9.1517e-08  -5.2591e-11   3.1223e-03  -1.2700e-03  -1.9162e-05  -4.5190e-04;
  131.    9.8922e-07   5.6846e-10  -1.7184e-03   6.9980e-04   9.5218e-04   5.0433e-04];
  132. c2 = [
  133.    3.6519e-16  -2.0788e-17  -2.0891e+01  -1.1933e+02  -4.1720e+02   4.6024e+02;
  134.    9.2135e-16   1.6474e-15  -1.6124e+01   1.6319e+02   5.2691e+02  -4.0617e+02];
  135. c3 = [
  136.   -4.8842e+00   5.2742e+00   1.4817e+00   7.4520e-03;
  137.   -2.8106e+00  -1.8255e+00  -4.2506e-02   1.3483e+00];
  138. c = [c1 c2 c3];
  139. d = [
  140.    1.6729e-01   1.5477e-02
  141.    1.5477e-02  -1.6729e-01];
  142. disp('     First 6 column of the A matrix:')
  143. a1
  144. disp('                                 (strike a key to continue ...)')
  145. pause
  146. clc
  147. disp('     Column 7 to 12 of the A matrix:')
  148. a2
  149. disp('                                 (strike a key to continue ...)')
  150. pause
  151. clc
  152. disp('     Column 13 to 16 of the A matrix:')
  153. a3
  154. disp('                                 (strike a key to continue ...)')
  155. pause
  156. clc
  157. b
  158. disp('                                 (strike a key to continue ...)')
  159. pause
  160. clc
  161. c
  162. disp('                                 (strike a key to continue ...)')
  163. pause
  164. clc
  165. d
  166. disp('                                 (strike a key to continue ...)')
  167. pause
  168. clc
  169. disp('                 << Example 1. Optimal Descriptor Hankel MDA >>')
  170. disp('       -----------------------------------------------------------------')
  171. disp('        First let us compute the Hankel singular values of the model...')
  172. disp('  ')
  173. disp('        [hsv] = hksv(a,b,c);    % Computing Hankel singular values')
  174. disp('       -----------------------------------------------------------------')
  175. [hsv] = hksv(a,b,c);
  176. hsv = hsv';
  177. format compact
  178. hsv
  179. format loose
  180. disp(' ')
  181. disp('      ----------------------------------------------------------------')
  182. disp('       Obviously, the "balanced" version of Hankel MDA fails immedi-')
  183. disp('       ately due to the last few Hankel singular vlaues being "0" ...')
  184. disp('      ----------------------------------------------------------------')
  185. disp(' ')
  186. disp('        (strike a key to see a bar chart of the Hankel Singular Values ...)')
  187. pause
  188. bar(hsv)
  189. title('Hankel Singular Values')
  190. pause
  191. clc
  192. disp(' ')
  193. disp(' ')
  194. disp('      ---------------------------------------------------------------')
  195. disp('       If we truncate the model to 6-state, the infinity-norm error ') 
  196. disp('       bound can be computed beforehand ...')
  197. disp(' ')
  198. disp('       bndsch = 2*sum(hsv(1,7:16));  % Infinity-norm error bound')
  199. disp('      ---------------------------------------------------------------')
  200. bndsch = 2*sum(hsv(1,7:16))
  201. disp('                                 (strike a key to continue ...)')
  202. pause
  203. clc
  204. disp(' ')
  205. disp('  ')
  206. disp('     ------------------------------------------------------------------')
  207. disp('      Let us start Optimal Descriptor Hankel Model Reduction to get a')
  208. disp('      6-state reduced order model ...')
  209. disp(' ')
  210. disp('      [ah,bh,ch,dh,bndhkl,hsv] = ohklmr(a,b,c,d,1,6); % Run OHKLMR.M')
  211. disp('     ------------------------------------------------------------------')
  212. [ah,bh,ch,dh,bndhkl,hsv] = ohklmr(a,b,c,d,1,6);
  213. disp('                                                 ^')
  214. disp('        - - - Computing SV Bode plot of G(s) and G(s) - - - ')
  215. w = logspace(-4,4,100);
  216. svg = sigma(a,b,c,d,1,w); svg = 20*log10(svg);
  217. svhkl = sigma(ah,bh,ch,dh,1,w); svhkl = 20*log10(svhkl);
  218. [aher,bher,cher,dher] = addss(a,b,c,d,ah,bh,-ch,-dh);
  219. disp('  ')
  220. disp('        - - - Computing SV Bode plot of the error system - - -')
  221. sverhkl = sigma(aher,bher,cher,dher,1,w); sverhkl = 20*log10(sverhkl);
  222. bndhkl = 20*log10(bndhkl)*ones(w);
  223. disp(' ')
  224. disp('                                                               ^')
  225. disp('                       (strike a key to see the plot of G(s) & G(s)...)')
  226. pause
  227. clg
  228. semilogx(w,svg,w,svhkl)
  229. title('Optimal Descriptor Hankel Model Reduction (16-state --> 6-state)')
  230. xlabel('Frequency - Rad/Sec')
  231. ylabel('SV - db')
  232. text(0.5,-120,'Original 16-State Model')
  233. text(0.001,-40,'6-State Approximation')
  234. grid
  235. pause
  236. disp(' ')
  237. disp('                                                 ^')
  238. disp('           (strike a key to see the plot of G(s)-G(s) and its bound...)')
  239. pause
  240. clg
  241. semilogx(w,bndhkl,w,sverhkl)
  242. title('Optimal Descriptor Hankel MDA (Infinity Error Bound)')
  243. xlabel('Frequency - Rad/Sec')
  244. ylabel('SV - db')
  245. text(1,-45,'Error Bound')
  246. text(1000,-100,'Error')
  247. grid
  248. pause
  249. clc
  250. disp(' ')
  251. disp('  ')
  252. disp('       << Example 2. Optimal Descriptor Hankel Anticausal Projection >>')
  253. disp('     -------------------------------------------------------------------')
  254. disp('       Let us try the same algorithm but do an anticausal projection')
  255. disp('       this time. This is also referred to the 0-th order Hankel MDA ..')
  256. disp(' ')
  257. disp('       [ahx,bhx,chx,dhx,ahy,bhy,chy,dhy,aug] = ohkapp(a,b,c,d,1,0); ')
  258. disp('                                                     % Run OHKAPP.M')
  259. disp('       State-space of the anticausal approximant: (ahy,bhy,chy,dhy)')
  260. disp('     -------------------------------------------------------------------')
  261. [ahx,bhx,chx,dhx,ahy,bhy,chy,dhy,aug] = ohkapp(a,b,c,d,1,0); 
  262. disp(' ')
  263. disp(' ')
  264. disp('                                 (strike a key to continue ...)')
  265. pause
  266. clc
  267. disp(' ')
  268. disp('    -----------------------------------------------------------------')
  269. disp('     Let us show that the maximum singular value of the total error ')
  270. disp('          ^')
  271. disp('     G(s)-G(s)  is an all-pass function ...')
  272. disp('             +')
  273. disp('           ^')
  274. disp('     where G(s)  is the anticausal Hankel approximant of G(s).')
  275. disp('              +')
  276. disp('    -----------------------------------------------------------------')
  277. disp('      ')
  278. disp(' ')
  279. disp('                                                         ^')
  280. disp('        - - - Computing singular value Bode plot of G(s)-G(s) - - -')
  281. disp('                                                            +')
  282. [aher,bher,cher,dher] = addss(a,b,c,d,ahy,bhy,-chy,-dhy);
  283. svher = sigma(aher,bher,cher,dher,1,w); svher = 20*log10(svher);
  284. disp(' ')
  285. disp(' ')
  286. disp('                (strike a key to see the "all-pass" error function ...)')
  287. pause
  288. subplot(211)
  289. semilogx(w,svher(1,:))
  290. title('Optimal Descriptor Hankel Anticausal Projection (Error System)')
  291. xlabel('Frequency - Rad/Sec')
  292. ylabel('max(SV) - db')
  293. subplot(212)
  294. semilogx(w,svher(2,:))
  295. xlabel('Frequency - Rad/Sec')
  296. ylabel('min(SV) - db')
  297. grid
  298. pause
  299. subplot(111)
  300. %
  301. % ----- End of OHKDEMO.M --- RYC/MGS 
  302.