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
|
# Copyright (C) 2000 by Peter J. Gawthrop
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
###############################################################
# Tell user
Sys=$1
method=$2
if [ -z "$method" ]; then
method=implicit
fi
#
echo "Creating $1_sim.m with $method integration method"
if [ $method = "implicit" ]; then
ode=cse
odeo=cseo
else
ode=ode
odeo=odeo
fi
# Find system constants
Nx=`mtt_getsize $Sys x` # States
Nu=`mtt_getsize $Sys u` # Inputs
Ny=`mtt_getsize $Sys y` # Inputs
# Header
lang_header -noglobal $1 sim m 'x0,u,t,par' '[y,x]' > $1_sim.m
cat >> $1_sim.m <<EOF
## Initialise
ui = zero_input($Nu); # Zero the input
xi = x0; # Read in initial state
yi = $1_cseo(xi,ui,t(1),par); # First output
##Sizes
N = length(t);
## Initialise arrays
x = zeros($Nx,N);
y = zeros($Ny,N);
y(:,1) = yi(:);
for i = 1:N
ti = t(i);
ui = u(i);
y(:,i) = yi;
x(:,i) = xi;
dxi = $1_cse(xi,ui,ti,par); # State derivative
A = $1_smxa(xi,ui,dt,par); # (I-Adt)
A = reshape(A,$Nx,$Nx);
Ax = $1_smxax(xi,ui,dt,par); # (I-Adt)x
#open = eval(sprintf("%s_switchopen(x);", system_name)); # Open switches
#x = mtt_implicit(x,dx,A,Ax,dt,$Nx,zeros(20,1)); # Implicit update
xi = A\(Ax + dxi*dt); # Implicit update
yi = $1_cseo(xi,ui,ti,par); # Output
|
>
>
>
>
|
>
>
>
>
|
>
|
>
>
>
>
|
|
<
<
<
|
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
|
|
>
>
<
>
|
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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
|
# Copyright (C) 2000 by Peter J. Gawthrop
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.1 2000/04/08 10:43:26 peterg
## Initial revision
##
###############################################################
# Tell user
Sys=$1
method=$2
computation=$3
if [ -z "$method" ]; then
method=implicit
fi
if [ -n "$computation" ]; then
blurb="for language $computation"
fi
echo "Creating $1_sim.m with $method integration method $blurb"
if [ $method = "implicit" ]; then
ode=cse
odeo=cseo
else
ode=ode
odeo=odeo
fi
# Find system constants
Nx=`mtt_getsize $Sys x` # States
Nu=`mtt_getsize $Sys u` # Inputs
Ny=`mtt_getsize $Sys y` # Inputs
Npar=`wc -l $Sys\_sympar.txt | awk '{print $1}'`
# Header
lang_header -noglobals $1 sim m 'x0,u,t,par,sensitivities' '[y,x]' > $1_sim.m
cat >> $1_sim.m <<EOF
if nargin<5
sensitivities = [];
endif;
## Initialise
[ui] = zero_input($Nu); # Zero the input
[xi] = x0; # Read in initial state
##Sizes
N = length(t);
## Initialise arrays
x = zeros($Nx,N);
y = zeros($Ny,N);
## Timing parameters
first = t(1);
dt = t(2) - t(1);
last = t(length(t))+dt;
n_sens = length(sensitivities);
EOF
if [ "$computation" = "c" ]; then
cat >> $1_sim.m <<EOF
## Create the command string
command = "";
for sensitivity_index = [0 sensitivities]
args = "";
for i=1:$Nx
args = sprintf("%s %g", args, x0(i));
endfor
par(sensitivities) = zeros(1,n_sens);
if sensitivity_index>0
par(sensitivity_index) = 1;
i_1 = $Ny/2 + 2;
else
i_1 = 2;
endif;
i_2 = i_1 + $Ny/2 -1;
for i=1:$Npar
args = sprintf("%s %g", args, par(i));
endfor
args = sprintf("%s %g %g %g", args, first, dt, last);
command = sprintf("%s ./$1_ode2odes.out %s | cut -f %i-%i;", command, args, i_1, i_2);
endfor;
command
## And execute it ...
yy=str2num(system(command));
[Ny,My] = size(yy);
## Reshape
Ny/N,N
y = reshape(yy',N,Ny/N)';
endfunction
EOF
else
cat >> $1_sim.m <<EOF
A = zeros($Nx,$Nx);
Ax = zeros($Nx,1);
dx = zeros($Nx,1);
## Step size
dt = t(2)-t(1);
iFirst = first/dt;
for i = 1:N
ti = t(i);
ui = u(:,i);
yi = $1_cseo(xi,ui,ti,par); # Output
if i>
y(:,i) = yi;
x(:,i) = xi;
dxi = $1_cse(xi,ui,ti,par); # State derivative
A = $1_smxa(xi,ui,dt,par); # (I-Adt)
A = reshape(A,$Nx,$Nx);
Ax = $1_smxax(xi,ui,dt,par); # (I-Adt)x
#open = eval(sprintf("%s_switchopen(x);", system_name)); # Open switches
#x = mtt_implicit(x,dx,A,Ax,dt,$Nx,zeros(20,1)); # Implicit update
xi = A\(Ax + dxi*dt); # Implicit update
endfor;
endfunction
EOF
fi
|