home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / IDDEMO3.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.0 KB  |  128 lines
  1. echo on
  2. clc
  3. %    This case study is based on the same "hairdyer" data as study # 1.
  4. %    The process consists of air being fanned through a tube. The air
  5. %    is heated at the inlet of the tube, and the input is the voltage
  6. %    applied to the heater. The output is the temperature at the outlet
  7. %    of the tube.
  8.  
  9. %    First, load the data:
  10.  
  11. load dryer2
  12.  
  13. %    Form a data set for estimation of the first half, and a reference
  14. %    set for validation purposes of the second half:
  15.  
  16. ze = [y2(1:500) u2(1:500)];
  17. zr = [y2(501:1000) u2(501:1000)];
  18.  
  19. %    Detrend each of the sets:
  20.  
  21. ze = dtrend(ze); zr = dtrend(zr);
  22.  
  23. %    Take a look at the data:
  24.  
  25. pause, idplot(ze,200:350), pause    % Press any key to continue
  26. clc
  27. %    We shall work with models of the following kind:
  28.  
  29. % y(t) + a1*y(t-1) + ... + ana*y(t-na) = b1*u(t-nk) + ... + bnb*u(t-nb-nk+1)
  30.  
  31. %    First we try and determine the time delay (nk). We then
  32. %    select a second order model (na=nb=2), and try out every time delay
  33. %    between 1 and 10. The loss function for the different models
  34. %    are computed using the reference data set:
  35.  
  36. V = arxstruc(ze,zr,struc(2,2,1:10));
  37.  
  38. %    We now select that delay that gives the best fit for the
  39. %    reference set:
  40.  
  41. [nn,Vm] = selstruc(V,0);
  42.  
  43. %    The chosen structure was
  44. nn
  45. pause % Press any key to continue.
  46. clc
  47. %    We can also check how the fit depends on the delay. The logarithms
  48. %    of a quadratic loss function are given as the first row, while
  49. %    the indexes na, nb and nk are given as a column below the correspon-
  50. %    ding loss function. 
  51.  
  52. pause % Press any key to continue.
  53. clc
  54. Vm
  55. %    The choice of three delays is thus rather clear. Now test the
  56. %    orders. We check the fit for all 25 combinations of up to 5
  57. %    b-parameters and up to 5 a-parameters, all with delay 3:
  58.  
  59. V = arxstruc(ze,zr,struc(1:5,1:5,3));
  60.  
  61. pause % Press any key to continue.
  62. clc
  63. %    The best fit for the reference data set is obtained for
  64.  
  65. nn = selstruc(V,0)
  66.  
  67. pause % Press any key to continue.
  68. clc
  69. %    It may be advisable to check yourself how much the fit is
  70. %    improved for the higher order models. The following plot
  71. %    shows the fit as a function of the number of parameters used.
  72. %    (Note that several different model structures use the same
  73. %    number of parameters.) You will then be asked to enter the
  74. %    number of parameters you think gives a sufficiently good fit.
  75. %    The routine then selects that structure with this many parameters
  76. %    that gives the best fit. 
  77.  
  78. pause, nns = selstruc(V)      % Press a key to continue.
  79.  
  80. %    The best fit is thus obtained for nn = [4 4 3], while we
  81. %    see that the improved fit compared to nn = [2 2 3] is rather
  82. %    marginal. 
  83.  
  84. pause % Press any key to continue.
  85. clc
  86. %    Let's compute the fourth order model:
  87.  
  88. th4 = arx(ze,[4 4 3]);
  89.  
  90. %    We now check the pole-zero configuration for the fourth order
  91. %    model. We also include confidence regions for the poles and zeros
  92. %    corresponding to 3 standard deviations. 
  93.  
  94. pause, zpplot(th2zp(th4),3), pause      % Press any key to for plot.
  95. clc
  96. %    It is thus quite clear that the two complex conjugated pole-zero
  97. %    pairs are likely to cancel, so that a second order model would be
  98. %    adequate. We thus compute
  99.  
  100. th2 = arx(ze,[2 2 3]);
  101.  
  102. %    We can test how well this model is capable of reproducing the
  103. %    reference data set. To compare the simulated output from the
  104. %    model with the actual output (plotting the mid 200 data points)
  105. %    we invoke (press key for plot)
  106.  
  107. pause,compare(zr,th2,150:350);pause
  108.  
  109. clc
  110. %    We may also check the residuals ("leftovers") of this model,
  111. %    i.e., what is left unexplained by the model.
  112.  
  113. pause % Press any key to continue.
  114.  
  115. e = resid(ze,th2); pause
  116.  
  117. plot(e), title('The residuals'),pause
  118. clc
  119. %    We see that the residual are quite small compared to the signal
  120. %    level of the output, that they are reasonably (although not
  121. %    perfectly) well uncorrelated with the input and between themselves.
  122. %    We can thus be (provisionally) satisfied with the model th2.
  123.  
  124. pause % Press any key to return to main menu.
  125. echo off
  126. set(gcf,'NextPlot','replace');
  127.  
  128.