* Bug #11980 fixed - csim returned an error message in some cases: "Submatrix
[scilab.git] / scilab / modules / cacsd / macros / csim.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA - 
3 // 
4 // This file must be used under the terms of the CeCILL.
5 // This source file is licensed as described in the file COPYING, which
6 // you should have received as part of this distribution.  The terms
7 // are also available at    
8 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
9
10
11 function [y,x]=csim(u,dt,sl,x0,tol)
12     //Syntax:
13     //  [y [,x]]=csim(u,dt,sl,[x0]) 
14     // simulation of the controlled linear system sl.
15     // sl is assumed to be a continuous-time system.
16     // u is the control and x0 the initial state.
17     //
18     //u can be:
19     // - a function 
20     //    [inputs]=u(t)
21     // - a list
22     //    list(ut,parameter1,....,parametern) such that
23     //    inputs=ut(t,parameter1,....,parametern)
24     // - the character string 'impuls' for impulse response calculation
25     //    (here sl is assumed SISO without direct feedthrough and x0=0)
26     // - the character string 'step' for step response calculation 
27     //    (here sl is assumed SISO without direct feedthrough and x0=0)
28     //dt is a vector of instants with dt(1) = initial time
29     //                   that is:           x0=x
30     //                                          dt(1)
31     //
32     //y matrix such that:
33     //  y=[y       y  ...  y     ]
34     //      dt(1)   dt(2)   dt(n)
35     //x matrix such that:
36     //  x=[x       x  ...  x     ]
37     //      dt(1)   dt(2)   dt(n)
38     //
39     //See also:
40     // dsimul flts ltitr rtitr ode impl
41     //!
42
43     [lhs,rhs]=argn(0)
44     //
45     if rhs<3 then error(39),end
46     sltyp=typeof(sl)
47     if and(sltyp<>['state-space' 'rational']) then 
48         error(msprintf(_("%s: Wrong type for input argument #%d: %s data structure expected.\n"),"csim",3,"syslin"))
49     end
50     if sltyp=='rational' then sl=tf2ss(sl),end
51     if sl.dt<>'c' then 
52         warning(msprintf(gettext("%s: Input argument %d is assumed continuous time.\n"),"csim",1));
53     end
54     //
55     [a,b,c,d]=sl(2:5);
56     if type(d)==2&degree(d)>0 then 
57         d=coeff(d,0);
58         warning(msprintf(gettext("%s: Direct feedthrough set to its zero degree coefficient.\n"),"csim"));
59     end
60     [ma,mb]=size(b);
61     //
62     imp=0;step=0
63     text='if t==0 then y=0, else y=1,end'
64     //
65     select type(u)
66     case 10 then //input given by its type (step or impuls)
67         if mb<>1 then error(95,1);end;
68         if part(u,1)=='i' then
69             //impuse response
70             imp=1;
71             if norm(d,1)<>0 then
72                 warning(msprintf(gettext("%s: Direct feedthrough set to zero.\n"),"csim"));
73                 d=0*d;
74             end;
75         elseif part(u,1)=='s' then
76             step=1
77             if norm(d,1)<>0 then
78                 warning(msprintf(gettext("%s: Direct feedthrough set to zero.\n"),"csim"));
79                 d=0*d;
80             end;
81         else
82             error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),"csim",1,"""step"",""impuls"""))
83         end;
84         deff('[y]=u(t)',text);
85     case 11 then //input given by a function of time
86         comp(u)
87     case 13 then //input given by a function of time
88     case 1 then //input given by a vector of data
89         [mbu,ntu]=size(u);
90         if mbu<>mb | ntu<>size(dt,'*') then 
91             error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),"csim",1,2))
92         end
93     case 15 then  //input given by a list: function of time with parameters
94         uu=u(1),
95         if type(uu)==11 then 
96             comp(uu),
97             u(1)=uu,
98         end
99     else error(44,2)
100     end;
101     //
102     if rhs==3 then x0=sl(6),end
103     if imp==1|step==1 then x0=0*x0,end
104     nt=size(dt,'*');x=0*ones(ma,nt);
105     [a,v]=balanc(a);
106     v1=v;//save for backward transformation
107
108     //apply transformation u without matrix inversion
109     [k,l]=find(v<>0) //get the permutation
110
111     //apply right transformation 
112     v=v(k,l);//diagonal matrix
113     c=c(:,k)*v; 
114     //apply left transformation 
115     v=diag(1 ./diag(v));b=v*b(k,:);x0=v*x0(k)
116
117     [a,v2,bs]=bdiag(a,1);b=v2\b;c=c*v2;x0=v2\x0;
118     //form the differential equation function
119     if type(u)==1 then
120         //form a continuous time interpolation of the given data
121         ut=u;
122         if min(size(ut))==1 then ut=matrix(ut,1,-1),end
123         deff('[y]=u(t)',['ind=find(dt<=t);nn=ind($)'
124         'if (t==dt(nn)|nn==nt) then '
125         '   y=ut(:,nn)'
126         'else '
127         '   y=ut(:,nn)+(t-dt(nn))/(dt(nn+1)-dt(nn))*(ut(:,nn+1)-ut(:,nn))'
128         'end']);
129         deff('[ydot]=%sim2(%tt,%y)','ydot=ak*%y+bk*u(%tt)');
130     elseif type(u)<>15 then
131         deff('[ydot]=%sim2(%tt,%y)','ydot=ak*%y+bk*u(%tt)');
132         ut=ones(mb,nt);for k=1:nt, ut(:,k)=u(dt(k)),end
133     else
134         %sim2=u
135         tx=' ';for l=2:size(u), tx=tx+',%'+string(l-1);end;
136         deff('[ydot]=sk(%tt,%y,u'+tx+')','ydot=ak*%y+bk*u(%tt'+tx+')');
137         %sim2(0)=sk;u=u(1)
138         deff('[ut]=uu(t)',...
139         ['['+part(tx,3:length(tx))+']=%sim2(3:'+string(size(%sim2))+')';
140         'ut=ones(mb,nt);for k=1:nt, ut(:,k)=u(t(k)'+tx+'),end'])
141         ut=uu(dt);
142     end;
143
144     //simulation
145     k=1;
146     for n=bs', 
147         kk=k:k+n-1;
148         ak=a(kk,kk);
149         bk=b(kk,:);
150         nrmu=max([norm(bk*ut,1),norm(x0(kk))]);
151         if nrmu > 0 then
152             if rhs<5 then
153                 atol=1.d-12*nrmu;rtol=atol/100;
154             else
155                 atol=tol(1);rtol=tol(2);
156             end
157             xkk=ode('adams',x0(kk),dt(1),dt,rtol,atol,%sim2);
158             if size(xkk,2)<>size(x,2) then
159                 error(msprintf(_("%s: Simulation failed before final time is reached.\n"),"csim"))
160             end
161             x(kk,:)=xkk;
162             if imp==1 then x(kk,:)=ak*x(kk,:)+bk*ut,end
163         end;
164         k=k+n
165     end;
166     if imp==0 & step==0 then y=c*x+d*ut, else y=c*x,end
167     if lhs==2 then x=v1*v2*x,end
168 endfunction