home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / IDDEMO7.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  8.0 KB  |  268 lines
  1. echo on
  2. %     In this demo we shall demonstrate how to use several the SITB to 
  3. %    estimate paramaters in state-space models. We shall investigate 
  4. %    data produced by a (simulated) dc-motor. We first load the data:
  5.  
  6. load dcm-data
  7.  
  8. %     The matrix y contains the two outputs: y1 is the angular position of 
  9. %    the motor shaft and y2 is the angular velocity. There are 400 data 
  10. %    points and the sampling interval is 0.1 seconds. The input is 
  11. %     contained in the vector u. It is the input voltage to the motor.
  12. %    Press a key for plots of the input-output data.
  13.  
  14. z=[y u];
  15.  
  16. pause,idplot(z,[],0.1,2),pause
  17.  
  18. clc
  19.  
  20. %     We shall build a model of the dc-motor. The dynamics of the motor is 
  21. %    well known. If we choose x1 as the angular position and x2 as the 
  22. %    angular velocity it is easy to set up a state-space model of the 
  23. %    following character neglecting disturbances:
  24. %    (see page 84 in Ljung(1987):
  25.  
  26. %           | 0     1  |      | 0   |
  27. %  d/dt x = |          | x  + |     | u
  28. %           | 0   -th1 |      | th2 |
  29. %
  30. %
  31. %          | 1  0 |
  32. %     y =  |      | x
  33. %          | 0  1 |  
  34. %
  35. %
  36. %     The parameter th1 is here the inverse time-constant of the motor and 
  37. %    th2 is such that th2/th1 is the static gain from input to the angular
  38. %    velocity. (See Ljung(1987) for how th1 and th2 relate to the physical 
  39. %    parameters of the motor.)
  40. %     We shall estimate these two parameters from the observed data. Strike 
  41. %    a key to continue.
  42.  
  43. pause
  44.  
  45. %     1. FREE PARAMETERS
  46. %
  47. %     We first have to define the structure in terms of the general 
  48. %    description
  49. %
  50. %     d/dt x = A x + B u + K e
  51. %
  52. %        y   = C x + D u + e
  53. %
  54. %     Strike a key to continue.
  55.  
  56. pause
  57.  
  58. %
  59. %     Any parameter to be estimated is entered as NaN (Not a Number). 
  60. %    Thus we have
  61. A=[0 1
  62.    0 NaN];
  63. B=[0
  64.    NaN];
  65. C=[1 0
  66.    0 1];
  67. D=[0
  68.    0];
  69. K=[0 0
  70.    0 0];
  71. X0=[0
  72.     0];   %X0 is the initial value of the state vector; it could also be
  73.           %entered as parameters to be identified.
  74.  
  75. %     The model structure is now defined by
  76.  
  77. ms = modstruc(A,B,C,D,K,X0);pause     %    Strike a key to continue
  78.  
  79. %     We shall produce an initial guess for the parameters. Let us guess 
  80. %    that the time constant is one second and that the static gain is 0.28.
  81. %    This gives:
  82.  
  83. th_guess = [-1 0.28];
  84.  
  85. %     We now collect all this information into the "theta-format" by
  86.  
  87. dcmodel = ms2th(ms,'c',th_guess,[],0.1);
  88.  
  89. %     'c' denotes that the parametrization refers to a continuous-time model
  90. %     The last argument, 0.1,is the sampling interval for the collected data
  91.  
  92. %     The prediction error (maximum likelihood) estimate of the parameters 
  93. %    is now computed by (Press a key to continue)
  94.  
  95. pause, dcmodel = pem(z,dcmodel); 
  96.  
  97. %     We can display the result by the 'present' command.
  98. %     The imaginary parts denote standard deviation. 
  99. %    Strike a key to continue.
  100.  
  101. pause, present(dcmodel), pause
  102.  
  103. %     The estimated values of the parameters are quite close to those used 
  104. %    when the data were simulated (-4 and 1).
  105.  
  106. %     To evaluate the model's quality we can simulate the model with the 
  107. %    actual input by and compare it with the actual output.
  108. %    Strike a key to continue.
  109.  
  110. pause, compare(z,dcmodel); pause
  111.  
  112. %     We can now, for example plot zeros and poles and their uncer-
  113. %     tainty regions. We will draw the regions corresponding to 10 standard
  114. %     deviations, since the model is quite accurate. Note that the pole at 
  115. %    the origin is absolutely certain, since it is part of the model 
  116. %    structure; the integrator from angular velocity to position. 
  117. %    Strike a key to continue.
  118.  
  119. pause, zpplot(th2zp(dcmodel),10), pause
  120.  
  121. %     2. COUPLED PARAMETERS
  122. %
  123. %     Suppose that we accurately know the static gain of the dc-motor (from
  124. %     input voltage to angular velocity, e.g. from a previous step-response
  125. %     experiment. If the static gain is G, and the time constant of the 
  126. %    motor is t, then the state-space model becomes
  127. %
  128. %            |0     1|    |  0  |
  129. %  d/dt x =  |       |x + |     | u
  130. %            |0  -1/t|    | G/t | 
  131. %
  132. %            |1   0|
  133. %     y   =  |     | x
  134. %            |0   1|
  135. %
  136.  
  137. pause % Press a key to continue
  138.  
  139. %     With G known, there is a dependance between the entries in the 
  140. %    different matrices. In order to describe that, the earlier used way 
  141. %    witn "NaN" will not be sufficient. We thus have to write an m-file 
  142. %    which produces the A,B,C,D,K and X0 matrices as outputs, for each 
  143. %    given parameter vector as input. It also takes auxiliary arguments as 
  144. %    inputs, so that the user can change certain things in the model 
  145. %    structure, without having to edit the m-file. In this case we let the 
  146. %    known static gain G be entered as such an argument. The m-file that 
  147. %    has been written has the name 'motor'. Press a key for a listing.
  148.  
  149. pause, type motor
  150.  
  151. pause % Press a key to continue.
  152.  
  153. %     We now create a "THETA-structure" corresponding to this model 
  154. %    structure:The assumed time constant will be
  155.  
  156. tc_guess = 1;
  157.  
  158. %     We also give the value 0.25 to the auxiliary variable G (gain)
  159. %     and sampling interval
  160.  
  161. aux = 0.25;
  162.  
  163. dcmm = mf2th('motor','c',tc_guess,aux,[],0.1);
  164.  
  165. pause   % Press a key to continue
  166.  
  167. %     The time constant is now estimated by
  168.  
  169. dcmm = pem([y u],dcmm);
  170.  
  171. %     We have thus now estimated the time constant of the motor directly. 
  172. %    Its value is in good agreement with the previous estimate. 
  173. %    Strike a key to continue.
  174.  
  175. pause, present(dcmm),  pause
  176.  
  177. %     With this model we may now proceed to test various aspects as before. 
  178. %    The syntax of all the commands is identical to the previous case.
  179.  
  180.  
  181. %     3. MULTIVARIATE ARX-MODELS.
  182. %
  183. %     The state-space part of the toolbox also handles multivariate (several
  184. %     outputs) ARX models. By a multivariate ARX-model we mean the 
  185. %    following: Strike a key to continue.
  186.  
  187. pause
  188.  
  189. %      A(q) y(t) = B(q) u(t) + e(t)
  190. %
  191. %     Here A(q) is a ny | ny matrix whose entries are polynomials in the 
  192. %    delay operator 1/q. The k-l element is denoted by a_kl(q):
  193. %
  194. %                          -1                   -nakl
  195. %     a_kl(q) = 1 + a-1 q    + .... + a-nakl q    
  196. %
  197. %     It is thus a polynomial in 1/q of degree nakl.
  198.  
  199. %     Similarily B(q) is a ny | nu matrix, whose kj-element is
  200.  
  201. %                      -nkkj         -nkkj-1                 -nkkj-nbkj
  202. %     b_kj(q) = b-0 q      +  b-1 q        + ... + b-nbkj q
  203.  
  204. %     There is thus a delay of nkkj from input number j to output number k
  205.  
  206. %     The most common way to create those would be to use the ARX-command.
  207. %     The orders are specified as
  208. %     nn = [na nb nk]
  209. %     with na being a ny | ny matrix whose kj-entry is nakj; nb and nk are
  210. %     defined anlogously.  Strike a key to continue.
  211.  
  212. pause
  213.  
  214. %     Let's test some ARX-models on the dc-data. First we could simply build
  215. %     a general second order model:
  216.  
  217. na = [ 2 2; 2 2];
  218. nb = [2; 2];
  219. nk = [1;1];
  220.  
  221. dcarx1 = arx([y u],[na nb nk]);
  222.  
  223. %     The result, dcarx1, is stored in the THETA-format, and all previous 
  224. %    commands apply. We could for example explicitly determine the 
  225. %    ARX-polynomials by 
  226.  
  227. [A,B] = th2arx(dcarx1);
  228.  
  229. %      Strike a key to continue.
  230.  
  231. pause, A, B, pause
  232.  
  233. %     We could also test a structure, where we know that y1 is obtained by
  234. %     filtering y2 through a first order filter. (The angle is the integral
  235. %     of the angular velocity). We could then also postulate a first order
  236. %     dynamics from input to output number 2:
  237.  
  238. na = [1 1; 0 1];
  239. nb = [0 ; 1];
  240. nk = [1 ; 1];
  241.  
  242. dcarx2 = arx(z,[na nb nk]);
  243. dcarx2 = sett(dcarx2,0.1);   % Setting the sampling interval to 0.1 seconds
  244. [Am, Bm] = th2arx(dcarx2);
  245.  
  246. Am, Bm
  247.  
  248. %     Strike e key to continue
  249.  
  250. pause
  251.  
  252. %     Finally, we could compare the bodeplots obtained from the input to
  253. %     output one for the different models by
  254.  
  255. g1 = th2ff(dcmodel); % The first calculated model
  256. g2 = th2ff(dcmm);    % The second model (known dc-gain)
  257. g3 = th2ff(dcarx2);  % The second arx-model
  258.  
  259. pause,bodeplot([g1 g2 g3]),pause
  260.  
  261. %     The two first models are in more or less exact agreement. The 
  262. %    arx-models are not so good, due to the bias caused by the non-white 
  263. %    equation error noise. (We had white measurement noise in the 
  264. %    simulations).
  265. pause
  266. set(gcf,'NextPlot','replace');
  267.  
  268.