home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / THINIT.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.9 KB  |  56 lines
  1. function eta=thinit(eta0,R,par,sp);
  2. %THINIT gives (random) initial values to the parameters, with guaranteed 
  3. %    predictor and/or system stability
  4. %
  5. %    TH = thinit(TH0)
  6. %or
  7. %    TH = thinit(TH0,R,PAR,SP)
  8. %
  9. %    TH0:Original theta corresponding to state-space model (See help th(ss)) 
  10. %    TH: Modified TH0 with new initial parameter values
  11. %    R: variance of random initial parameters 
  12. %       (row vector with dim = dimension of parameter vector)
  13. %    PAR: mean of initial parameters (PAR =[] gives default)
  14. %    AP: SP='p' :stability of predictor assured
  15. %        SP='s' :stability of system is assured
  16. %        SP='b' :both predictor and system will be stable
  17. %
  18. %    Defaults: SP='p', PAR=the parameter values of TH0, R=[1 ... 1]
  19.  
  20. %    L. Ljung 10-2-1990
  21. %    Copyright (c) 1990 by the MathWorks, Inc.
  22. %    All Rights Reserved
  23.  
  24. if ~isthss(eta0),error('This routine is for state-space model structures only!'),end
  25. [nr,nc]=size(eta0);
  26. sspm=getmfth(eta0);
  27. T=eta0(1,2); Tmod=abs(T);
  28. nd=eta0(1,5);
  29. ms=getargth(eta0);
  30. if nargin<2;R=[];end
  31. if nargin<3;par=[];end
  32. if nargin<4;sp=[];end
  33.  
  34. if isempty(R),R=ones(nd,1);end,if isempty(sp),sp='p';end
  35. if isempty(par),par=eta0(3,1:nd);end
  36.  
  37. if length(R)~=nd,error('The length of the vector R must equal the number of free parameters in th!'),end
  38. if length(par)~=nd,error('The length of the parameter vector must equal the number of free parameters in th!'),end
  39. test=2;nr=1;if norm(R)==0,noR=1;else noR=0;end
  40. while test>1 & nr<50
  41.     parval=rand(1,nd)*diag(sqrt(R))+par;
  42.     if any(eta0(2,8)==[2 3]),Tmod=-1;end
  43.     [a,b,c,d,k]=feval(sspm,parval,Tmod,ms);
  44.  
  45.     if sp=='s',test=max(abs(eig(a)));
  46.     elseif sp=='p',test=max(abs(eig(a-k*c)));
  47.     else test=max([max(abs(eig(a))), max(abs(eig(a-k*c)))]);
  48.     end
  49. nr=nr+1; if nr==50,error('50 samples drawn without satisfying the stability conditions'),end
  50. if noR,nr=50;if test>1,disp('WARNING: System/predictor unstable'),end,end
  51. end
  52. eta=eta0;
  53. eta(3,1:nd)=parval;
  54. eta(2,7)=27;
  55. if eta(2,8)==3,eta(2,7)=37;end
  56.