diff options
Diffstat (limited to 'modules/cacsd/macros/arsimul.sci')
| -rwxr-xr-x | modules/cacsd/macros/arsimul.sci | 130 |
1 files changed, 130 insertions, 0 deletions
diff --git a/modules/cacsd/macros/arsimul.sci b/modules/cacsd/macros/arsimul.sci new file mode 100755 index 000000000..c18f46fc1 --- /dev/null +++ b/modules/cacsd/macros/arsimul.sci @@ -0,0 +1,130 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) INRIA - +// +// 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 z=arsimul(x1,x2,x3,x4,x5,x6,x7,x8) + // function z=arsimul(a,b,d,sig,u,up,yp,ep) + + [lhs,rhs]=argn(0) + // switch to ar representation + if type(x1)<>15&type(x1)<>16 then + if rhs < 5, + error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"arsimul",5,8)); + end; + ar=armac(x1,x2,x3,size(x1,"r"),size(x5,"r"),x4); + if rhs==5,z=arsimul(ar,x5);return;end + if rhs==6,z=arsimul(ar,x5,x6);return;end + if rhs==7,z=arsimul(ar,x5,x6,x7);return;end; + if rhs==8,z=arsimul(ar,x5,x6,x7,x8);return;end; + end + // Here the call is always arsimul(ar,....) + a=x1("a");b=x1("b");d=x1("d");sig=x1("sig"); + u=x2; + [bl,bc]=size(b);[al,ac]=size(a);[dl,dc]=size(d); + [mmu,Nu]=size(u); + if mmu<>x1("nu") then + error(msprintf(gettext("%s: Number of rows of %s are incompatible with %s object.\n"),"arsimul","u","arma")); + return; + end; + // X = [y_{n-1},y_{n-2},...y_{n-ka},u_{n-1},....,u_{n-kb},e_{n-1},....,e_{n-kd} + a1=a(:,al+1:$);al1=size(a1,"c"); + b1=b(:,mmu+1:$);bl1=size(b1,"c"); + d1=d(:,al+1:$);dl1=size(d1,"c"); + A=[-a1,b1,d1]; + // y_{n} = A*X + b(:,1:mmu)*u_{n}+d(:,1:al)*e_{n} + // in the system fff x=[y_n;X]; + if A==[] then + deff("[xkp1]=fff(k,x)",... + ["ukp1=u(:,k+1);dkp1=br(:,k+1);"; + "xkp1= b(:,1:mmu)*ukp1+d(:,1:al)*dkp1"]); + else + deff("[xkp1]=fff(k,x)",... + ["x=x(al+1:$);ukp1=u(:,k+1);dkp1=br(:,k+1);"; + "ykp1= A*x + b(:,1:mmu)*ukp1+d(:,1:al)*dkp1"; + "xkp1=[];" + "if al1>0; xkp1=[ykp1;x(1:al1-al)];end;"; + "if bl1>0; xkp1=[xkp1;ukp1;x(al1+1:al1+bl1-mmu)];end;"; + "if dl1>0; xkp1=[xkp1;dkp1;x(al1+bl1+1:al1+bl1+dl1-al)];end;"; + "xkp1=[ykp1;xkp1];" ]); + end + // Noise simulation. + br=sig*rand(al,Nu,"normal"); + //br=[-2,1,0.5] + // Initial condition + // the first call to fff will be fff(0,x) + // x must be set to + // [ y_{0},...y{-ak},u_{0},...u_{-bk},d_{0},...d_{-dk} + // where ak= al1/al -1; bk= bl1/mmu -1 ; dk = dl1/al-1 + // past conditions for up + //-------------------------- + if rhs <=2, + up=0*ones(bl1,1); + else + up=x3; + if bl1==0 then + if up<>[] then error(msprintf(gettext("%s: Wrong size for input argument #%d: An empty matrix expected.\n"),"arsimul",3)) + return ; + end; + else + up_s=size(up) + if up_s(1)<>mmu|up_s(2)<>(bl1/mmu) then + error(msprintf(gettext("%s: %s must be of dimension (%s, %s).\n"),"arsimul","up=[u(0),u(-1),..,]",string(mmu),string(bl1/mmu))); + return + end + up=matrix(up,bl1,1); + end + end + // past conditions for yp + //-------------------------- + if rhs <=3, + yp=0*ones(al1,1) + else + yp=x4; + if al1==0 then + if yp<>[] then error(msprintf(gettext("%s: Wrong size for input argument #%d: An empty matrix expected.\n"),"arsimul",4)) + return ; + end; + else + yp_s=size(yp); + if yp_s(1)<>al|yp_s(2)<>(al1/al) then + error(msprintf(gettext("%s: %s must be of dimension (%s, %s).\n"), "arsimul","yp=[y(0),y(-1),..,]",string(al), string(al1/al))); + return; + end + yp=matrix(yp,al1,1); + end + end + // past conditions for ep + //-------------------------- + if rhs <=4, + ep=0*ones(dl1,1); + else + ep=x5 + if dl1==0 then + if ep<>[] then error(msprintf(gettext("%s: Wrong size for input argument #%d: An empty matrix expected.\n"),"arsimul",5)) + return ; + end; + else + ep_s=size(ep); + if ep_s(1)<>al|ep_s(2)<>(dl1/al) then + error(msprintf(gettext("%s: %s must be of dimension (%s, %s).\n"), "arsimul","ep=[e(0),e(-1),..,]",string(al), string(dl1/al))); + return; + end + ep=matrix(ep,dl1,1); + end; + end; + xi=[yp;up;ep]; + // If A=[] it is a degenerate case which also work + // but xi must be set to a scalar value to provide proper + // result dimensions. + // + xi=[0*ones(al,1);xi]; + z=ode("discrete",xi,0,1:Nu,fff); + // Now z contains y_{1},.....y_{Nu}; + z=z(1:al,:) +endfunction |