home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / CS1.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  5.7 KB  |  210 lines
  1. if ~exist('decimate'),error('Sorry, this test case requires the Signal Processing Toolbox'),end
  2. echo on,clc
  3. %CASE STUDY NUMBER 1 : A GLASS TUBE DRAWING FACTORY
  4. %
  5. % In this Case Study we shall study data from a glass tube factory. The
  6. % experiments and the data are discussed in
  7. %
  8. % V. Wertz, G. Bastin and M. Heet: Identification of a glass tube
  9. % drawing bench. Proc. of the 10th IFAC Congress, Vol 10, pp 334-339
  10. % Paper number 14.5-5-2. Munich August 1987.
  11. %
  12. % The the output of the process is the thickness and the diameter of
  13. % the manufactured tube. The inputs are the air-pressure inside the
  14. % tube and the drawing speed.
  15. %
  16. % We shall in this tutorial study the process from the input speed
  17. % to the output thickness.
  18. %
  19. load thispe25.mat
  20. % 
  21. % First we look at the data. We split it into two halves, one for
  22. % estimation and one for validation:
  23. %
  24.  
  25. ze=thispe25(1001:1500,:);
  26. zv=thispe25(1501:2000,:);pause
  27. idplot(ze),pause
  28. % A close-up:
  29. pause
  30. idplot(ze,101:200),pause
  31. % Let us remove the mean values:
  32.  
  33. ze=dtrend(ze);
  34. zv=dtrend(zv);
  35.  
  36. % The sampling interval of the data is one second. We may perhaps
  37. % detect some rather high frequencies in the output. Let us therefore
  38. % first compute the input and output spectra:
  39.  
  40. sy=spa(ze(:,1));
  41. su=spa(ze(:,2));
  42. pause  % Red/solid is output spectrum and green/dashed is input
  43.  
  44. bodeplot([sy su])
  45. grid,pause
  46.  
  47. % We may note that the input has very little relative energy above 1 rad/sec
  48. % while a substantial part of the output's energy comes from frquencies above
  49. % 1 rad/sec. There are thus some high frequency disturbancies that may cause
  50. % some problem for the model building.
  51.  
  52. % We first compute the impulse response, using part of the data:
  53.  
  54. [ir,dum,cl]=cra(ze); pause
  55. plot(0:19,ir,'g',0:19,cl*ones(1,20),'-.r',0:19,-cl*ones(1,20),'-.r')
  56. title('The impulse response with confidence limits'),pause
  57.  
  58. % We see a quite substantial delay of about 12 samples in the impulse response.
  59.  
  60. % We also as a preliminary test compute the spectral analysis estimate:
  61.  
  62. [gs,phiv]=spa(ze);
  63. pause
  64. bodeplot(gs,1),pause
  65.  
  66. % We note, among other things, that the high frequency behavior is quite
  67. % uncertain.
  68.  
  69. % Let us do a quick check if we can pick up good dynamics by just computing
  70. % a fourth order ARX model using the estimation data and simulate that
  71. % model using the validation data. We know that the delay is about 12:
  72. pause  % Press a key to continue
  73.  
  74. m1=arx(ze,[4 4 12]);
  75. compare(zv,m1);pause
  76.  
  77. % A close up:
  78. pause
  79. compare(zv,m1,inf,101:200);pause
  80.  
  81. %
  82. % There are clear difficulties to deal with the high frequency components
  83. % of the output. That in conjunction with the long delay suggests that we
  84. % decimate the data  by four (i.e. low-pass filter it and pick every fourth 
  85. % value):
  86.  
  87. zd(:,1)=decimate(dtrend(thispe25(:,1)),4);
  88. zd(:,2)=decimate(dtrend(thispe25(:,2)),4);
  89.  
  90. zde=zd(1:500,:);
  91. zdv=zd(501:length(zd),:);
  92. pause
  93.  
  94. % Let us find a good structure for the decimated data.
  95.  
  96. % First compute the impulse response:
  97.  
  98. [ird,dum,cl]=cra(zde); pause
  99. plot(0:19,ird,'g',0:19,cl*ones(1,20),'-.r',0:19,-cl*ones(1,20),'-.r')
  100. title('The impulse response for decimated data with confidence limits'),pause
  101.  
  102. % We see that the delay is about 3 (which is consistent with what we saw
  103. % above.
  104.  
  105. % What does the "quick start" give?
  106.  
  107. pause
  108. compare(zdv,arx(zde,[4 4 3]));pause
  109.  
  110. % Seems much better that for the undecimated data. Let's go on to find
  111. % a good model structure. First we look for the delay:
  112.  
  113.  
  114. V = arxstruc(zde,zdv,struc(2,2,1:30));
  115. nn=selstruc(V,0)
  116.  
  117. % Then we test several different orders with and around this delay:
  118.  
  119. V = arxstruc(zde,zdv,struc(1:5,1:5,nn(3)-1:nn(3)+1));
  120.  
  121. % Make your own choice of orders based on the figure to follow.
  122. % (The "model quality" comments to be made below refer to the default choice)
  123.  
  124. pause,
  125. nn=selstruc(V)
  126.  
  127. % Let's compute and check the model for whatever model you just preferred:
  128.  
  129. m2=arx(zde,nn);
  130. m2=sett(m2,4); % The sampling interval
  131. pause
  132. compare(zdv,m2,inf,21:150);pause
  133.  
  134. % This seems reasonably good!
  135.  
  136. % We may also compare the bodeplots of the model m1 and m2:
  137.  
  138. g1=th2ff(m1);
  139. g2=th2ff(m2);
  140. pause
  141. bodeplot([g1 g2]),pause
  142.  
  143. % We see that the two models agree well up to the Nyquist frequency of
  144. % the slower sampled data.
  145.  
  146. % Let's test the residuals:
  147. pause
  148. resid(zdv,m2);pause
  149.  
  150. % This seems OK.
  151.  
  152. % What does the pole-zero diagram tell us?
  153.  
  154. zpplot(th2zp(m2),3,[],2),pause
  155.  
  156. % No clear indication of pole-zero cancellations!
  157.  
  158.  
  159. % Let us however also check if we can do better by modeling the noise
  160. % separately in a Box-Jenkins model:
  161.  
  162. mb=bj(zde,[nn(2) 2 2 nn(1) nn(3)]);
  163. mb=sett(mb,4); % The sampling interval
  164.  
  165. % Is this better than the ARX-model in simulation ?
  166. % (press key)
  167. pause
  168. clf reset
  169. subplot(211)
  170. compare(zdv,m2);
  171. subplot(212)
  172. compare(zdv,mb);pause
  173.  
  174. % The Box-Jenkins model is thus capable of describing the validation 
  175. % data set somewhat better. How about the Bode plots?
  176.  
  177. gb=th2ff(mb);
  178. pause
  179. clf reset
  180. bodeplot([g2 gb]),pause
  181.  
  182. % The models agree quite well, apart from the funny twist in the BJ model.
  183. % This twist, however, seems to have some physical relevance, since we
  184. % get better simulations in the latter case.
  185. % Finally we may compare the FPE:s of the two models:
  186.  
  187. [m2(2,1) mb(2,1)]
  188.  
  189. % To summarize: After decimation it is quite possible to build simple models
  190. % that are capable of reproducing the validation set in a good manner.
  191.  
  192. % We can however do quite well with a much simpler model; a [1 1 3] ARX model:
  193.  
  194. m3=arx(zde,[1 1 3]);
  195.  
  196. % Simulation:
  197.  
  198. pause
  199. clf reset,subplot(211),compare(zdv,m3);subplot(212),compare(zdv,mb);pause
  200.  
  201. % 5-step ahead prediction:
  202.  
  203. pause
  204. clf,subplot(211),compare(zdv,m3,5);subplot(212),compare(zdv,mb,5);pause
  205.  
  206. pause
  207. echo off
  208.  
  209.  
  210.