diff options
Diffstat (limited to 'modules/cacsd/macros/csim.sci')
| -rwxr-xr-x | modules/cacsd/macros/csim.sci | 165 |
1 files changed, 165 insertions, 0 deletions
diff --git a/modules/cacsd/macros/csim.sci b/modules/cacsd/macros/csim.sci new file mode 100755 index 000000000..d14fe4e3e --- /dev/null +++ b/modules/cacsd/macros/csim.sci @@ -0,0 +1,165 @@ +// 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 [y,x]=csim(u,dt,sl,x0,tol) + //Syntax: + // [y [,x]]=csim(u,dt,sl,[x0]) + // simulation of the controlled linear system sl. + // sl is assumed to be a continuous-time system. + // u is the control and x0 the initial state. + // + //u can be: + // - a function + // [inputs]=u(t) + // - a list + // list(ut,parameter1,....,parametern) such that + // inputs=ut(t,parameter1,....,parametern) + // - the character string 'impuls' for impulse response calculation + // (here sl is assumed SISO without direct feedthrough and x0=0) + // - the character string 'step' for step response calculation + // (here sl is assumed SISO without direct feedthrough and x0=0) + //dt is a vector of instants with dt(1) = initial time + // that is: x0=x + // dt(1) + // + //y matrix such that: + // y=[y y ... y ] + // dt(1) dt(2) dt(n) + //x matrix such that: + // x=[x x ... x ] + // dt(1) dt(2) dt(n) + // + //See also: + // dsimul flts ltitr rtitr ode impl + //! + + [lhs,rhs]=argn(0) + // + if rhs<3 then error(39),end + sltyp=typeof(sl) + if and(sltyp<>["state-space" "rational"]) then + error(msprintf(_("%s: Wrong type for input argument #%d: %s data structure expected.\n"),"csim",3,"syslin")) + end + if sltyp=="rational" then sl=tf2ss(sl),end + if sl.dt<>"c" then + warning(msprintf(gettext("%s: Input argument #%d is assumed continuous time.\n"),"csim",1)); + end + // + [a,b,c,d]=sl(2:5); + if degree(d)>0 then + error(msprintf(gettext("%s: Wrong type for input argument #%d: A proper system expected\n"),"csim",1)); + end + [ma,mb]=size(b); + // + imp=0;step=0 + text="if t==0 then y=0, else y=1,end" + // + select type(u) + case 10 then //input given by its type (step or impuls) + if mb<>1 then + error(msprintf(gettext("%s: Wrong type for input argument #%d: A SIMO expected.\n"),"csim",1)); + end; + if part(u,1)=="i" then + //impulse response + imp=1; + dt(dt==0)=%eps^2; + elseif part(u,1)=="s" then + step=1 + if norm(d,1)<>0 then + dt(dt==0)=%eps^2; + end; + else + error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),"csim",1,"""step"",""impuls""")) + end; + deff("[y]=u(t)",text); + case 11 then //input given by a function of time + comp(u) + case 13 then //input given by a function of time + case 1 then //input given by a vector of data + [mbu,ntu]=size(u); + if mbu<>mb | ntu<>size(dt,"*") then + error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),"csim",1,2)) + end + case 15 then //input given by a list: function of time with parameters + uu=u(1), + if type(uu)==11 then + comp(uu), + u(1)=uu, + end + else error(44,2) + end; + // + if rhs==3 then x0=sl(6),end + if imp==1|step==1 then x0=0*x0,end + nt=size(dt,"*");x=0*ones(ma,nt); + [a,v]=balanc(a); + v1=v;//save for backward transformation + + //apply transformation u without matrix inversion + [k,l]=find(v<>0) //get the permutation + + //apply right transformation + v=v(k,l);//diagonal matrix + c=c(:,k)*v; + //apply left transformation + v=diag(1 ./diag(v));b=v*b(k,:);x0=v*x0(k) + + [a,v2,bs]=bdiag(a,1);b=v2\b;c=c*v2;x0=v2\x0; + //form the differential equation function + if type(u)==1 then + //form a continuous time interpolation of the given data + ut=u; + if min(size(ut))==1 then ut=matrix(ut,1,-1),end + deff("[y]=u(t)",["ind=find(dt<=t);nn=ind($)" + "if (t==dt(nn)|nn==nt) then " + " y=ut(:,nn)" + "else " + " y=ut(:,nn)+(t-dt(nn))/(dt(nn+1)-dt(nn))*(ut(:,nn+1)-ut(:,nn))" + "end"]); + deff("[ydot]=%sim2(%tt,%y)","ydot=ak*%y+bk*u(%tt)"); + elseif type(u)<>15 then + deff("[ydot]=%sim2(%tt,%y)","ydot=ak*%y+bk*u(%tt)"); + ut=ones(mb,nt);for k=1:nt, ut(:,k)=u(dt(k)),end + else + %sim2=u + tx=" ";for l=2:size(u), tx=tx+",%"+string(l-1);end; + deff("[ydot]=sk(%tt,%y,u"+tx+")","ydot=ak*%y+bk*u(%tt"+tx+")"); + %sim2(0)=sk;u=u(1) + deff("[ut]=uu(t)",... + ["["+part(tx,3:length(tx))+"]=%sim2(3:"+string(size(%sim2))+")"; + "ut=ones(mb,nt);for k=1:nt, ut(:,k)=u(t(k)"+tx+"),end"]) + ut=uu(dt); + end; + + //simulation + k=1; + for n=bs', + kk=k:k+n-1; + ak=a(kk,kk); + bk=b(kk,:); + nrmu=max([norm(bk*ut,1),norm(x0(kk))]); + if nrmu > 0 then + if rhs<5 then + atol=1.d-12*nrmu;rtol=atol/100; + else + atol=tol(1);rtol=tol(2); + end + xkk=ode("adams",x0(kk),0,dt,rtol,atol,%sim2); + if size(xkk,2)<>size(x,2) then + error(msprintf(_("%s: Simulation failed before final time is reached.\n"),"csim")) + end + x(kk,:)=xkk; + if imp==1 then x(kk,:)=ak*x(kk,:)+bk*ut,end + end; + k=k+n + end; + y=c*x+d*ut + if lhs==2 then x=v1*v2*x,end +endfunction |