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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
|
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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
|
+
-
+
-
+
+
+
-
+
+
+
+
-
-
+
+
-
+
-
-
-
-
+
+
+
+
-
-
-
+
+
-
+
-
-
+
+
-
-
-
-
-
-
+
-
-
+
+
-
-
-
-
-
+
-
-
-
-
+
+
+
-
-
+
-
-
-
-
+
-
-
+
-
-
-
-
+
-
-
-
-
+
+
-
-
-
-
+
+
-
-
-
+
-
-
-
-
+
+
-
-
-
+
+
-
-
-
+
+
-
-
-
-
-
-
-
+
+
+
+
+
+
-
-
-
-
+
-
-
-
-
-
-
+
+
-
-
-
+
-
-
+
-
-
-
-
-
-
-
-
-
-
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
-
-
-
+
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
+
-
-
-
+
-
-
-
+
+
-
-
-
-
-
+
+
-
-
-
-
-
+
+
+
+
+
-
-
-
+
+
+
+
+
+
+
+
+
+
-
+
-
+
+
+
-
+
-
+
+
+
+
+
|
echo ${mtt_options} > ${option_file}
fi
}
## Make the _ident.m file
make_ident() {
filename=${sys}_${rep}.m
date=`date`
echo Creating ${filename}
cat > ${filename} <<EOF
function [y,u,t] = ${sys}_ident (last, ppp_names, par_names, A_u, A_w, w, Q, extras)
function [epar,Y] = ${sys}_ident (y,u,t,par_names,Q,extras)
## usage: [y,u,t] = ${sys}_ident (last, ppp_names, par_names, A_u, A_w, w, Q, extras)
## usage: [epar,Y] = ${sys}_ident (y,u,t,par_names,Q,extras)
##
## last last time in run
## ppp_names Column vector of names of ppp params
## par_names Column vector of names of estimated params
## extras Structure containing additional info
##
## Created by MTT on ${date}
## Sensitivity system name
system_name = "s${sys}"
##Sanity check
if nargin<2
printf("Usage: [y,u,t] = ${sys}_ident(N, ppp_names[, par_names, extras])\n");
if nargin<3
printf("Usage: [y,u,t] = ${sys}_ident(y,u,t,par_names,Q,extras);");
return
endif
if nargin<4
if nargin<6
## Set up optional parameters
extras.criterion = 1e-3;
extras.emulate_timing = 0;
extras.max_iterations = 15;
extras.simulate = 1;
extras.v = 1e-6;
extras.verbose = 0;
extras.max_iterations = 10;
extras.simulate = 2;
extras.v = 1e-2;
extras.verbose = 1;
extras.visual = 1;
endif
## System info
[n_x,n_y,n_u,n_z,n_yz] = ${sys}_def;
sympar = ${sys}_sympar;
simpar = ${sys}_simpar;
sympars = s${sys}_sympar;
simpars = s${sys}_simpar;
t_ol = simpar.last;
## Number of intervals needed
N = ceil(last/t_ol);
## Setpoints
if extras.verbose
## Initial model state
x_0 = zeros(2*n_x,1);
printf("Open-loop interval %3.2f \n", simpar.last);
printf(" -- using info in ${sys}_simpar.txt\n");
printf("PPP optimisation from %3.2f to %3.2f\n", simpars.first, simpars.last);
printf(" -- using info in s${sys}_simpar.txt\n");
endif
t_horizon = [simpars.first+simpars.dt:simpars.dt:simpars.last]';
w_s = ones(length(t_horizon)-1,1)*w';
## Setup the indices of the adjustable stuff
if nargin<2
i_ppp = []
else
i_ppp = ppp_indices (ppp_names,sympar,sympars); # Parameters
endif
n_ppp = length(i_ppp(:,1));
if nargin<3
i_par = []
else
i_par = ppp_indices (par_names,sympar,sympars); # Parameters
endif
## Do some simulations to check things out
## System itself
par = ${sys}_numpar;
x_0_ol = ${sys}_state(par);
[y_ol,x_ol, t_ol] = ${sys}_sim(x_0_ol, par, simpar, ones(1,n_u));
simpar_OL = simpar;
simpar_OL.last = simpars.last;
[y_OL,x_OL, t_OL] = ${sys}_sim(x_0_ol, par, simpar_OL, ones(1,n_u));
pars = s${sys}_numpar;
x_0_ppp = s${sys}_state(pars);
[y_ppp,y_par,x_ppp, t_ppp] = s${sys}_ssim(x_0_ppp, pars, simpars, ones(1,n_u));
simpar_PPP = simpars;
simpar_PPP.first = simpar.first;
[y_PPP,y_par,x_PPP, t_PPP] = s${sys}_ssim(x_0_ppp, pars, simpar_PPP, ones(1,n_u));
## Basis functions
Us = ppp_ustar(A_u,1,t_OL',0,0)';
if extras.visual #Show some graphs
figure(2);
grid; title("Outputs of ${sys}_sim and s${sys}_ssim");
plot(t_ol,y_ol, '*', t_ppp, y_ppp, '+', t_OL, y_OL, t_PPP, y_PPP);
figure(3);
## Create a table of the parameters
[n_par,m_par] = size(i_par);
grid; title("Basis functions");
plot(t_OL, Us);
endif
## Do it
[y,u,t,P,U,t_open,t_ppp,t_est,its_ppp,its_est] \
= ppp_nlin_run ("${sys}",i_ppp,i_par,A_u,w_s,N,Q,extras);
## Compute values at ends of ol intervals
T_open = cumsum(t_open);
T_open = T_open(1:length(T_open)-1); # Last point not in t
j=[];
for i = 1:length(T_open)
j = [j; find(T_open(i)*ones(size(t))==t)];
y_open = y(j,:);
u_open = u(j,:);
## Plots
gset nokey
gset nogrid
#eval(sprintf("gset xtics %g", simpar.last));
#gset noytics
gset format x "%.1f"
gset format y "%.2f"
gset term fig monochrome portrait fontsize 20 size 20 20 metric \
thickness 4
gset output "${sys}_ident.fig"
title("");
xlabel("Time (s)");
ylabel("u");
subplot(2,1,2); plot(t,u,'-', T_open, u_open,"+");
#subplot(2,1,2); plot(t,u);
ylabel("y");
subplot(2,1,1); plot(t,y,'-', T_open, y_open,"+");
#subplot(2,1,1); plot(t,y);
oneplot;
gset term fig monochrome portrait fontsize 20 size 20 10 metric \
thickness 4
gset output "${sys}_ident.basis.fig"
title("");
xlabel("Time (s)");
ylabel("Basis functions");
plot(t_OL, Us);
## Create plot against time
TTT = [ [0;T_open] [T_open; last] ]';
TT = TTT(:);
[n,m] = size(P);
if m>0
P = P(1:n-1,:); # Loose last point
PP = [];
for j=1:m
pp = [P(:,j) P(:,j)]';
PP = [PP pp(:)];
endfor
oneplot;
gset output "${sys}_ident.par.fig"
title("");
xlabel("Time (s)");
ylabel("Parameters");
plot(TT,PP);
endif
[n,m] = size(U);
if m>0 oneplot;
gset output "${sys}_ident.U.fig"
title("");
xlabel("Time (s)");
ylabel("U");
[n,m] = size(U);
U = U(1:n-1,:); # Loose last point
UU = [];
for j=1:m
uu = [U(:,j) U(:,j)]';
UU = [UU uu(:)];
endfor
plot(TT,UU);
endif
function [last, ppp_names, par_names, A_u, A_w, w, Q, extras] = ${sys}_ident_numpar
function [y,u,t,par_names,Q,extras] = ${sys}_ident_numpar;
## usage: [last, ppp_names, par_names, A_u, A_w, w, Q, extras] = ${sys}_ident_numpar ()
##
##
## usage: [y,u,t,par_names,Q,extras] = ${sys}_ident_numpar;
## Edit for your own requirements
## Created by MTT on ${date}
## Last time of run
last = 10;
## Specify basis functions
A_w = zeros(1,1);
n_ppp = ${nu};
A_u = ppp_aug(A_w,laguerre_matrix(n_ppp-1,2.0));
## Names of ppp parameters
ppp_names = "";
for i=1:n_ppp
name = sprintf("ppp_%i", i);
ppp_names = [ppp_names; name];
## This section sets up the data source
## simulate = 0 Real data (you supply ${sys}_ident_data.dat)
## simulate = 1 Real data input, simulated output
## simulate = 2 Unit step input, simulated output
simulate = 2;
## System info
[n_x,n_y,n_u,n_z,n_yz] = ${sys}_def;
simpars = s${sys}_simpar;
## Access or create data
if (simulate<2) # Get the real data
if (exist("${sys}_ident_data.dat")==2)
printf("Loading ${sys}_ident_data.dat\n");
load ${sys}_ident_data.dat
else
printf("Please create a loadable file ${sys}_ident_data.dat containing y,u and t\n");
return
endif
else
switch simulate
case 2 # Step simulation
t = [0:simpars.dt:simpars.last]';
u = ones(size(t));
otherwise
error(sprintf("simulate = %i not implemented", simulate));
endswitch
endif
if (simulate>0)
par = ${sys}_numpar();
x_0 = ${sys}_state(par);
dt = t(2)-t(1);
simpars.dt = dt;
simpars.last = t(length(t));
y = ${sys}_sim(zeros(n_x,1), par, simpars, u);
endif
## Default parameter names - Put in your own here
sympar = ${sys}_sympar; # Symbolic params as structure
par_names = struct_elements (sympar); # Symbolic params as strings
[n,m] = size(par_names); # Size the string list
## Sort by index
for [i,name] = sympar
par_names(i,:) = sprintf("%s%s",name, blanks(m-length(name)));
endfor
## Estimated parameters
par_names = [];
## Weights
Q = ones(${ny},1);
## Setpoint
w = zeros(${ny},1); w(1) = 1;
Set up optional parameters
extras.criterion = 1e-3;
extras.max_iterations = 15;
extras.simulate = 1;
extras.v = 1e-6;
extras.verbose = 0;
extras.visual = 0;
extras.max_iterations = 10;
extras.simulate = simulate;
extras.v = 1e-2;
extras.verbose = 1;
extras.visual = 1;
endfunction
EOF
}
make_dat2() {
## Inform user
echo Creating ${dat2_file}
## Use octave to generate the data
octave -q <<EOF
[last, ppp_names, par_names, A_u, A_w, w, Q, extras] = ${sys}_ident_numpar;
[y,u,t] = ${sys}_ident(last, ppp_names, par_names, A_u, A_w, w, Q, extras);
data = [t,y,u];
[y,u,t,par_names,Q,extras] = ${sys}_ident_numpar;
[epar,Y] = ${sys}_ident (y,u,t,par_names,Q,extras);
[N,M] = size(Y);
y_est = Y(:,M);
data = [t,y_est,u];
save -ascii ${dat2_file} data
EOF
## Tidy up the latex stuff - convert foo_123 to foo_{123}
cat ${sys}_ident_par.tex > mtt_junk
sed -e "s/_\([a-z0-9,]*\)/_{\1}/g" < mtt_junk >${sys}_ident_par.tex
rm mtt_junk
}
case ${lang} in
numpar.m)
## Make the numpar stuff
make_ident_numpar;
;;
m)
## Make the code
make_ident;
;;
dat2|fig|basis.fig|par.fig|U.fig)
dat2)
## The dat2 language (output data) & fig file
rm ${sys}_ident*.fig
ps|basis.ps|par.ps|U.ps)
ps)
figs=`ls ${sys}_ident*.fig | sed -e 's/\.fig//'`
echo $figs
|