1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
|
function [y,x]=csim(u,dt,sl,x0)
//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 chracter 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
//!
// Copyright INRIA
[lhs,rhs]=argn(0)
//
if rhs<3 then error(39),end
//-compat type(sl)<>15 retained for list/tlist compatibility
if type(sl)<>15&type(sl)<>16 then error(56,1),end
flag=sl(1)
select flag(1)
case 'lss' then ,
case 'r' then sl=tf2ss(sl)
else error(97,1),
end;
if sl(7)<>'c' then warning('csim: time domain is assumed continuous'),end
//
[a,b,c,d]=sl(2:5);
if type(d)==2°ree(d)>0 then d=coeff(d,0);warning('D set to constant');end
[ma,mb]=size(b);
//
imp=0;text='if t==0 then y=0, else y=1,end'
//
select type(u)
case 10 then //
if mb<>1 then error(95,1);end;
if part(u,1)=='i' then
imp=1;
if norm(d,1)<>0 then
warning('direct feedthrough (d) <> 0;set to zero');
d=0*d;
end;
end;
deff('[y]=u(t)',text);
case 11,comp(u)
case 13,
case 15 then
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 then x0=0*x0,end
nt=size(dt,'*');x=0*ones(ma,nt)
[a,v,bs]=bdiag(a,1);b=v\b;c=c*v;x0=v\x0;
//
if 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=maxi([norm(bk*ut,1),norm(x0(kk))])
if nrmu > 0 then
atol=1.d-7*nrmu;rtol=atol/100
x(kk,:)=ode('adams',x0(kk),dt(1),dt,rtol,atol,%sim2)
if imp==1 then x(kk,:)=ak*x(kk,:)+bk*ut,end
end;
k=k+n
end;
if imp==0 then y=c*x+d*ut,else y=c*x,end
if lhs==2 then x=v*x,end
|
