home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / IDDEMO5.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.8 KB  |  117 lines
  1. hold off
  2. %IDDEMO5
  3. echo on
  4. clc
  5. %    This demo describes the group of RECURSIVE (ON-LINE)
  6. %    algorithms.  The recursive M-files include: RPEM, RPLR,
  7. %    RARMAX, RARX, ROE, and RBJ.   These algorithms implement
  8. %    all the recursive algorithms described in Chapter 11 of
  9. %    Ljung(1987).
  10. %
  11. %    RPEM is the general Recursive Prediction Error Algorithm
  12. %    for arbitrary multiple-input-single-output models
  13. %    (the same models as PEM works for).
  14. %
  15. %    PRLR is the general Recursive PseudoLinear Regression method
  16. %    for the same family of models.
  17. %
  18. %    RARX is a more efficient version of RPEM (and RPLR) for the
  19. %    ARX-case.
  20. %
  21. %    ROE, RARMAX and RBJ are more efficient versions of RPEM for
  22. %    the OE, ARMAX, and BJ cases (compare these functions to the
  23. %    off-line methods).
  24.  
  25. pause    % Press any key to continue.
  26. clc
  27. %    ADAPTATION MECHANISMS:  Each of the algorithms implement the
  28. %    four most common adaptation principles:
  29. %
  30. %    KALMAN FILTER approach: The true parameters are supposed to
  31. %    vary like a random walk with incremental covariance matrix
  32. %    R1. 
  33. %
  34. %    FORGETTING FACTOR approach: Old measurements are discounted
  35. %    exponentially. The base of the decay is the forgetting factor
  36. %    lambda.
  37. %
  38. %    GRADIENT method: The update step is taken as a gradient step
  39. %    of length gamma (th_new=th_old + gamma*psi*epsilon).
  40. %
  41. %    NORMALIZED GRADIENT method: As above, but gamma is replaced by
  42. %    gamma/(psi'*psi). The Gradient methods are also
  43. %    known as LMS (least mean squares) for the ARX case.
  44.  
  45. pause    % Press any key to continue.
  46. clc
  47. %    Let's pick a model and generate some input-output data:
  48.  
  49. u = sign(randn(50,1)); e = 0.2*randn(50,1);
  50. th0 = poly2th([1 -1.5 0.7],[0 1 0.5],[1 -1 0.2]);
  51. y = idsim([u e],th0); z = [y u];
  52.  
  53. pause, idplot(z), pause   % Press a key to see data.
  54. clc
  55. %    First we build an Output-Error model of the data we just
  56. %    plotted.  Use a second order model with one delay, and apply
  57. %    the forgetting factor algorithm with lambda = 0.98:
  58.  
  59. thm1 = roe(z,[2 2 1],'ff',0.98);  % It may take a while ....
  60.  
  61. %    The four parameters can now be plotted as functions of time.
  62.  
  63. pause, plot(thm1), pause  % Press a key to see plot.
  64. clc
  65. %    The true values are as follows: (Press a key for plot)
  66.  
  67. pause, hold on, plot(ones(50,1)*[1 0.5 -1.5 0.7]), pause
  68. hold off
  69. clc
  70. %    Now let's try a second order ARMAX model, using the RPLR
  71. %    approach (i.e ELS) with Kalman filter adaptation, assuming
  72. %    a parameter variance of 0.001:
  73.  
  74. thm2 = rplr(z,[2 2 2 0 0 1],'kf',0.001*eye(6));
  75.  
  76. pause    % Press any key for plot.
  77. plot(thm2), axis([0 50 -2 2]), pause
  78. clc
  79. % The true values are as follows: (Press any key for plot)
  80.  
  81. pause, hold on, plot(ones(50,1)*[-1.5 0.7 1 0.5 -1 0.2]), pause
  82. hold off
  83. clc
  84. %    So far we have assumed that all data are available at once.
  85. %    We are thus studying the variability of the system rather
  86. %    than doing real on-line calculations. The algorithms are
  87. %    also prepared for such applications, but they must then
  88. %    store more update information. The conceptual update then 
  89. %    becomes:
  90.  
  91. %  1. Wait for measurements y and u.
  92. %  2. Update: [th,yh,p,phi] = rarx([y u],[na nb nk],'ff',0.98,th',p,phi)
  93. %  3. Use th for whatever on-line application required.
  94. %  4. Go to 1.
  95.  
  96. %    Thus the previous estimate th is fed back into the algorithm
  97. %    along with the previous value of the "P-matrix" and the data
  98. %    vector phi.
  99.  
  100. %    We now do an example of this where we plot just the current
  101. %    value of th. You can look at the code after the plot is done.
  102.  
  103. pause    % Press any key to continue
  104.  
  105. [th,yh,p,phi] = rarx(z(1,:),[2 2 1],'ff',0.98);
  106. plot(1,th(1),'*',1,th(2),'+',1,th(3),'o',1,th(4),'*'),axis([1 50 -2 2]),hold on;
  107.  
  108. for k = 2:50
  109.     [th,yh,p,phi] = rarx(z(k,:),[2 2 1],'ff',0.98,th',p,phi);
  110.     plot(k,th(1),'*',k,th(2),'+',k,th(3),'o',k,th(4),'*')
  111. end,pause
  112. pause
  113. hold off
  114. echo off
  115. set(gcf,'NextPlot','replace');
  116.  
  117.