home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / IDDEMO2.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.9 KB  |  150 lines
  1. echo on
  2. clc
  3. %    In this example we compare several different methods of 
  4. %    identification.
  5.  
  6. %    We start by forming some simulated data.  Here are the
  7. %    numerator and denominator coefficients of the transfer
  8. %    function model of a system:
  9.  
  10. B = [0 1 0.5];
  11. A = [1 -1.5 0.7];
  12.  
  13. %    These can be put into the special "theta" structure used by
  14. %    the toolbox:
  15.  
  16. th0 = poly2th(A,B,[1 -1 0.2]);
  17.  
  18. %    The third argument, [1 -1 0.2], gives a characterization
  19. %    of the disturbances that act on the system.
  20.  
  21. pause    % Press any key to continue.
  22. clc
  23.  
  24. %    Now we generate an input signal U, a disturbance signal E,
  25. %    and simulate the response of the model to these inputs:
  26.  
  27. u = sign(randn(350,1));
  28. e = randn(350,1);
  29. y = idsim([u e],th0);
  30. z = [y u]; % A matrix with the output and inputs as column vectors.
  31.  
  32. %    Plot first 100 values of the input U and the output Y.
  33.  
  34. pause, idplot(z,1:100), pause  % Press any key for plot.
  35. clc
  36. %    Now that we have corrupted, simulated data, we can estimate models
  37. %    and make comparisons.  Let's start with spectral analysis.  We
  38. %    invoke SPA which returns a spectral analysis estimate of
  39. %    the transfer function GS and the noise spectrum NSS:
  40.  
  41. [GS,NSS] = spa(z);
  42.  
  43. pause, bodeplot(GS), pause   % Press any key for plot. 
  44. clc
  45.  
  46. %    Now use the Least Squares method to find an ARX-model with
  47. %    two poles, one zero, and a single delay on the input:
  48.  
  49. a2 = arx(z,[2 2 1]);
  50.  
  51. pause, present(a2), pause   % Press any key to present results.
  52. clc
  53.  
  54. %    Now use the instrumental variable method (IV4) to find a model
  55. %    with two poles, one zero, and a single delay on the input:
  56.  
  57. i2 = iv4(z,[2 2 1]);
  58.  
  59.  
  60. pause, present(i2), pause     % Press any key to present results.
  61. clc
  62. %    Calculate transfer functions for the two models obtained by
  63. %    ARX and IV4.  Show Bode plots comparing the Spectral Analysis
  64. %    estimate with the ARX and IV4 estimates:
  65.  
  66. [Ga2,NSa2] = th2ff(a2);        % From ARX
  67. Gi2 = th2ff(i2);            % From IV4
  68.  
  69. pause, bodeplot([GS Ga2 Gi2]), pause  % Press any key for plots.
  70.  
  71. clc
  72. %    Calculate and plot residuals for the model obtained by IV4:
  73.  
  74. e = resid(z,i2); pause
  75. plot(e), title('Residuals, IV4 method'), pause
  76. clc
  77.  
  78. %    Now compute second order ARMAX model:
  79.  
  80. am2 = armax(z,[2 2 2 1]);
  81.  
  82. pause    % Press any key to continue.
  83. clc
  84. %    Now compute second order BOX-JENKINS model:
  85.  
  86. bj2 = bj(z,[2 2 2 2 1]);
  87.  
  88. %    Compute transfer functions and noise spectra for these models:
  89.  
  90. [Gam2,NSam2] = th2ff(am2);
  91. [Gbj2,NSbj2] = th2ff(bj2);
  92.  
  93. %    Compare these transfer functions with those obtained by the
  94. %    IV4 method and the spectral analysis method:
  95.  
  96. pause    % Press any key for plots.
  97.  
  98. bodeplot([Gam2 Gbj2 Gi2 GS]), pause
  99. bodeplot([NSam2 NSbj2 NSS]), pause, clc
  100.  
  101. %    Plot residuals for the ARMAX and BOX-JENKINS models
  102.  
  103. e1 = resid(z,am2); pause
  104. e2 = resid(z,bj2); pause
  105. plot([e1 e2]), title('ARMAX and BOX-JENKINS residuals'), pause
  106. clc
  107. %    Simulate the ARMAX and BOX-JENKINS models with the real input
  108.  
  109. yham2 = idsim(u,am2);
  110. yhbj2 = idsim(u,bj2);
  111.  
  112. pause    % Press any key for plot.
  113.  
  114. plot([y yham2 yhbj2]), ...
  115.  title('Real output and simulated model outputs '), pause
  116. clc
  117. %    Now we'll compute the poles and zeros of the ARMAX and
  118. %    BOX-JENKINS models
  119.  
  120. zpam2 = th2zp(am2);
  121. zpbj2 = th2zp(bj2);
  122.  
  123.     % Press any key for plot.
  124. pause    % Press keys to advance poles and zeros on plot, from the ARMAX
  125.     % model to the Box-Jenkins model.
  126.  
  127. zpplot([zpam2 zpbj2]), pause, clc
  128.  
  129. % Clearly, AM2 and BJ2 both give good residuals and simulated outputs.
  130. % They are also in good mutual agreement. Let's check which has the
  131. % smaller AIC:  
  132.  
  133. [am2(2,1) bj2(2,1)]
  134. pause    % Press any key to continue.
  135. clc
  136. %    Since we generated the data,  we enjoy the luxury of comparing
  137. %    the model to the true system:
  138.  
  139. [G0,NS0] = th2ff(th0);
  140. zp0 = th2zp(th0);
  141.  
  142. pause    % Press any key for plots of comparisons.
  143.  
  144. bodeplot([G0 Gam2]), pause
  145. bodeplot([NS0 NSam2]), pause
  146. zpplot([zp0 zpam2]), pause
  147. echo off
  148. set(gcf,'NextPlot','replace');
  149.  
  150.