← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to macros/auto/csim.sci

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
macros/auto/csim.sci
 
1
function [y,x]=csim(u,dt,sl,x0,tol)
 
2
//Syntax:
 
3
//  [y [,x]]=csim(u,dt,sl,[x0]) 
 
4
// simulation of the controlled linear system sl.
 
5
// sl is assumed to be a continuous-time system.
 
6
// u is the control and x0 the initial state.
 
7
//
 
8
//u can be:
 
9
// - a function 
 
10
//    [inputs]=u(t)
 
11
// - a list
 
12
//    list(ut,parameter1,....,parametern) such that
 
13
//    inputs=ut(t,parameter1,....,parametern)
 
14
// - the character string 'impuls' for impulse response calculation
 
15
//    (here sl is assumed SISO without direct feedthrough and x0=0)
 
16
// - the character string 'step' for step response calculation 
 
17
//    (here sl is assumed SISO without direct feedthrough and x0=0)
 
18
//dt is a vector of instants with dt(1) = initial time
 
19
//                   that is:           x0=x
 
20
//                                          dt(1)
 
21
//
 
22
//y matrix such that:
 
23
//  y=[y       y  ...  y     ]
 
24
//      dt(1)   dt(2)   dt(n)
 
25
//x matrix such that:
 
26
//  x=[x       x  ...  x     ]
 
27
//      dt(1)   dt(2)   dt(n)
 
28
//
 
29
//See also:
 
30
// dsimul flts ltitr rtitr ode impl
 
31
//!
 
32
// Copyright INRIA
 
33
[lhs,rhs]=argn(0)
 
34
//
 
35
if rhs<3 then error(39),end
 
36
 
 
37
if type(sl)<>16 then error(56,1),end
 
38
flag=sl(1)
 
39
select flag(1)
 
40
  case 'lss' then ,
 
41
  case  'r'  then sl=tf2ss(sl)
 
42
  else  error(97,1),
 
43
end;
 
44
if sl(7)<>'c' then warning('csim: time domain is assumed continuous'),end
 
45
//
 
46
[a,b,c,d]=sl(2:5);
 
47
if type(d)==2&degree(d)>0 then d=coeff(d,0);warning('D set to constant');end
 
48
[ma,mb]=size(b);
 
49
//
 
50
imp=0;text='if t==0 then y=0, else y=1,end'
 
51
//
 
52
select type(u)
 
53
 case 10 then //
 
54
    if mb<>1 then error(95,1);end;
 
55
    if part(u,1)=='i' then
 
56
           imp=1;
 
57
           if norm(d,1)<>0 then
 
58
               warning('direct feedthrough (d) <> 0;set to zero');
 
59
               d=0*d;
 
60
           end;
 
61
    end;
 
62
    deff('[y]=u(t)',text);
 
63
 case 11,comp(u)
 
64
 case 13,
 
65
 case 1 then
 
66
   [mbu,ntu]=size(u);
 
67
   if mbu<>mb | ntu<>size(dt,'*') then error('wrong size of u'), end
 
68
 case 15 then
 
69
    uu=u(1),
 
70
    if type(uu)==11 then 
 
71
      comp(uu),
 
72
      u(1)=uu,
 
73
    end
 
74
 else error(44,2)
 
75
end;
 
76
//
 
77
if rhs==3 then x0=sl(6),end
 
78
if imp==1 then x0=0*x0,end
 
79
nt=size(dt,'*');x=0*ones(ma,nt)
 
80
[a,v,bs]=bdiag(a,1);b=v\b;c=c*v;x0=v\x0;
 
81
//
 
82
if type(u)==1 then
 
83
   ut=u;
 
84
   deff('[y]=u(t)',['ind=find(dt<=t)','nn=ind($)',...
 
85
        'if (t==dt(nn)|nn==nt) then y=ut(nn),...
 
86
         else y=ut(nn)+(t-dt(nn))/(dt(nn+1)-dt(nn))*(ut(nn+1)-ut(nn)), end']);
 
87
   deff('[ydot]=%sim2(%tt,%y)','ydot=ak*%y+bk*u(%tt)');
 
88
elseif type(u)<>15 then
 
89
  deff('[ydot]=%sim2(%tt,%y)','ydot=ak*%y+bk*u(%tt)');
 
90
  ut=ones(mb,nt);for k=1:nt, ut(:,k)=u(dt(k)),end
 
91
else
 
92
  %sim2=u
 
93
  tx=' ';for l=2:size(u), tx=tx+',%'+string(l-1);end;
 
94
  deff('[ydot]=sk(%tt,%y,u'+tx+')','ydot=ak*%y+bk*u(%tt'+tx+')');
 
95
  %sim2(0)=sk;u=u(1)
 
96
  deff('[ut]=uu(t)',...
 
97
     ['['+part(tx,3:length(tx))+']=%sim2(3:'+string(size(%sim2))+')';
 
98
      'ut=ones(mb,nt);for k=1:nt, ut(:,k)=u(t(k)'+tx+'),end'])
 
99
ut=uu(dt);
 
100
end;
 
101
//simulation
 
102
k=1;
 
103
for n=bs',
 
104
  kk=k:k+n-1
 
105
  ak=a(kk,kk)
 
106
  bk=b(kk,:)
 
107
  nrmu=maxi([norm(bk*ut,1),norm(x0(kk))])
 
108
  if nrmu > 0 then
 
109
    if rhs<5 then
 
110
      atol=1.d-10*nrmu;rtol=atol/100
 
111
    else
 
112
      atol=tol(1);rtol=tol(2)
 
113
    end
 
114
    x(kk,:)=ode('adams',x0(kk),dt(1),dt,rtol,atol,%sim2)
 
115
    if imp==1 then x(kk,:)=ak*x(kk,:)+bk*ut,end
 
116
  end;
 
117
  k=k+n
 
118
end;
 
119
if imp==0 then y=c*x+d*ut,else y=c*x,end
 
120
if lhs==2 then x=v*x,end
 
121