diff options
Diffstat (limited to 'modules/cacsd/macros/armax1.sci')
| -rwxr-xr-x | modules/cacsd/macros/armax1.sci | 124 |
1 files changed, 124 insertions, 0 deletions
diff --git a/modules/cacsd/macros/armax1.sci b/modules/cacsd/macros/armax1.sci new file mode 100755 index 000000000..95eafcf74 --- /dev/null +++ b/modules/cacsd/macros/armax1.sci @@ -0,0 +1,124 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) INRIA +// Copyright (C) ENPC - J-Ph. Chancelier +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + +function [arc,resid]=armax1(r,s,q,y,u,b0f) + // [arc, resid] = armax1(r, s, q, y, u, [b0f]) + // + // Compute the coefficient of monodimensional ARMAX + // A(z^-1)y= B(z^-1)u + D(z^-1)sig*e(t) + // e(t) is a white noise of variance 1 + // A(z)= 1+a1*z+...+a_r*z^r; ( r=0 => A(z)=1) + // B(z)= b0+b1*z+...+b_s z^s ( s=-1 => B(z)=0) + // D(z)= 1+d1*z+...+d_q*z^q ( q=0 => D(z)=1) + // + // Intput: + // y: output process + // u: input process + // r, s and q: auto-regression orders r >= 0, s >= 1 and moving average q. + // b0f is a optional parameter. By default, b0f is 0 and it means that + // this parameter must be identified. If b0f is 1, then b0f is supposed + // to be zero and is not identified. + // + // Output: + // arc is the tlist with fields + // a is the vector [1,a1,...,a_r] + // b is the vector [b0, ..., b_s] + // d is the vector [1, d1, ..., d_q] + // sig is the estimated standard deviation + // resid = [sig*echap(1), ..., ] + // + // Method: + // Cfre : Eykhoff (trends and progress in system identification) page 91 + // Introducing + // z(t)=[y(t-1),..,y(t-r),u(t),...,u(t-s),e(t-1),...,e(t-q)] and + // coef= [-a1,..,-ar,b0,...,b_s,d1,...,d_q]', we obtain + // y(t)= coef'* z(t) + sig*e(t). + // We use the sequential version of the AR estimation where e(t-i) is + // replaced by its estimated (Method RLLS). + // If q=0, it is a sequential version of the least squares algorithm given + // in armax function + + [lhs,rhs]=argn(0) + if rhs<=5,b0f=0;end + if s==-1,b0f=0;end // Seems not natural, but makes things work + u=matrix(u,1,-1);y=matrix(y,1,-1); //make u and y row vectors + [n1,n2]=size(y) + if size(y,"*")<>size(u,"*") then + error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same numbers of elements expected.\n"),.. + "armax1",4,5)); + end + // + t0=max(max(r,s+1),1)+1; + if r<>0;XTM1=y((t0-1):-1:(t0-r));else XTM1=[];end + if s<>-1;UTM1=u(t0-b0f:-1:(t0-s));else UTM1=[];end + if q<>0;ETM1=0*ones(1,q);else ETM1=[];end + npar=r+s+1-b0f+q + CTM1=0*ones(npar,1); + ZTM1=[XTM1,UTM1,ETM1]'; + PTM1=10.0*eye(npar,npar); + + for t=t0+1:n2, + if r<>0;XT=[ y(t-1), XTM1(1:(r-1))];else XT=[];end + if s<>-1;UT=[ u(t-b0f), UTM1(1:(s-b0f))];else UT=[];end + eeTM1=y(t-1)- CTM1'*ZTM1; + if q<>0;ET=[ eeTM1, ETM1(1:(q-1))];else ET=[];end + ZT=[XT,UT,ET]'; + // + KT=PTM1*ZT*(1/(1+ ZT'*PTM1*ZT)) + CT=CTM1+KT*(y(t)-ZT'*CTM1) + PT=PTM1-KT*ZT'*PTM1 + XTM1=XT;UTM1=UT;CTM1=CT;ETM1=ET;ZTM1=ZT;PTM1=PT; + end + // The coefficient a, b and d are extracted + // + if r<>0;a=[1;-CT(1:r)]';else a=1;end + if s<>-1; + if b0f==1,b=[0;CT(r+1:(r+s+1-b0f))]';else + b=[CT(r+1:(r+s+1-b0f))]';end + if q<>0;d=[1;CT(r+s+2-b0f:(r+s+q+1-b0f))]';else d=[1];end + else + b=0; + if q<>0;d=[1;CT(r+s+1-b0f:(r+s+q-b0f))]';else d=[1];end + end + // Simulation to get the prediction error + // + [sig,resid]=epred(r,s,q,CTM1,y,u,b0f); + arc=armac(a,b,d,1,1,sig); + + + +endfunction + +function [sig,resid]=epred(r,s,q,coef,y,u,b0f) + //============================================= + // [sig,resid] = epred(r,s,q,coef,y,u,b0f) + // Used by armax1 function to compute the prediction error + // coef= [-a1,..,-ar,b0,...,b_s,d1,...,d_q]' + // or + // coef= [-a1,..,-ar,b1,...,b_s,d1,...,d_q]' si b0f=1 + //! + [n1,n2]=size(y); + t0=max(max(r,s+1),1)+1; + if r<>0;XTM1=y((t0-1):-1:(t0-r));else XTM1=[];end + if s<>-1;UTM1=u(t0-b0f:-1:(t0-s));else UTM1=[];end + if q<>0;ETM1=0*ones(1,q);else ETM1=[];end + npar=r+s+1-b0f+q + ZTM1=[XTM1,UTM1,ETM1]'; + resid=0*ones(1,n2); + for t=t0+1:n2, + if r<>0;XT=[ y(t-1), XTM1(1:(r-1))];else XT=[];end + if s<>-1;UT=[ u(t-b0f), UTM1(1:(s-b0f))];else UT=[];end + resid(t)=y(t-1)- coef'*ZTM1; + if q<>0;ET=[ resid(t), ETM1(1:(q-1))];else ET=[];end + ZT=[XT,UT,ET]'; + XTM1=XT;UTM1=UT;ETM1=ET;ZTM1=ZT; + end + sig=sqrt(sum(resid.*resid)/size(resid,"*")) +endfunction |