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
|
function [y,x,u,t,UU,UU_c,UU_l] = ppp_nlin_sim (system_name,i_ppp,i_par,A_u,w_s,N_ol,extras)
## usage: [y,x,u,t,U,U_c,U_l] = ppp_nlin_sim (system_name,A_u,tau,t_ol,N,w)
##
##
## Simulate nonlinear PPP
## Copyright (C) 2000 by Peter J. Gawthrop
## Defaults
if nargin<7
extras.U_initial = "zero";
extras.U_next = "continuation";
extras.criterion = 1e-5;
extras.max_iterations = 10;
extras.v = 0.1;
extras.verbose = 0;
endif
## Names
s_system_name = sprintf("s%s", system_name);
## System details -- defines simulation within ol interval
par = eval(sprintf("%s_numpar;", system_name));
simpar = eval(sprintf("%s_simpar;", system_name));
x_0 = eval(sprintf("%s_state(par);", system_name));
[n_x,n_y,n_u] = eval(sprintf("%s_def;", system_name));
## Sensitivity system details -- defines moving horizon simulation
simpars = eval(sprintf("%s_simpar;", s_system_name));
sympars = eval(sprintf("%s_sympar;", s_system_name));
|
|
|
>
<
|
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
|
function [y,x,u,t,p,UU,UU_c,UU_l,t_ppp,t_est] = ppp_nlin_sim (system_name,i_ppp,i_par,A_u,w_s,N_ol,extras)
## usage: [y,x,u,t,p,UU,UU_c,UU_l,t_ppp,t_est] = ppp_nlin_sim (system_name,i_ppp,i_par,A_u,w_s,N_ol,extras)
##
##
## Simulate nonlinear PPP
## Copyright (C) 2000 by Peter J. Gawthrop
## Defaults
if nargin<7
extras.U_initial = "zero";
extras.U_next = "continuation";
extras.criterion = 1e-5;
extras.max_iterations = 10;
extras.v = 0.1;
extras.verbose = 0;
extras.estimate = 1;
endif
## Names
s_system_name = sprintf("s%s", system_name);
## System details -- defines simulation within ol interval
par = eval(sprintf("%s_numpar;", system_name));
simpar = eval(sprintf("%s_simpar;", system_name));
x_0 = eval(sprintf("%s_state(par);", system_name));
[n_x,n_y,n_u] = eval(sprintf("%s_def;", system_name));
## Sensitivity system details -- defines moving horizon simulation
simpars = eval(sprintf("%s_simpar;", s_system_name));
pars = eval(sprintf("%s_numpar;", s_system_name));
## Times
## -- within opt horizon
n_Tau = round(simpars.last/simpars.dt);
dtau = simpars.dt;
Tau = [0:n_Tau-1]'*dtau;
|
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
|
[k_x,k_w,K_x,K_w] = ppp_lin(A,B,C,D,A_u,0,tau');
## Main simulation loop
y = [];
x = [];
u = [];
t = [];
t_last = 0;
UU = [];
UU_l =[];
UU_c =[];
x_0s = zeros(2*n_x,1);
if strcmp(extras.U_initial,"linear")
U = K_w*w - K_x*x_0;
elseif strcmp(extras.U_initial,"zero")
U = zeros(n_U,1);
else
|
>
>
>
>
>
>
|
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
|
[k_x,k_w,K_x,K_w] = ppp_lin(A,B,C,D,A_u,0,tau');
## Main simulation loop
y = [];
x = [];
u = [];
t = [];
p = [];
t_last = 0;
UU = [];
UU_l =[];
UU_c =[];
t_ppp = [];
t_est = [];
x_0s = zeros(2*n_x,1);
if strcmp(extras.U_initial,"linear")
U = K_w*w - K_x*x_0;
elseif strcmp(extras.U_initial,"zero")
U = zeros(n_U,1);
else
|
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
|
tick = time;
if extras.max_iterations>0
[U, U_all, Error, Y] = ppp_nlin(system_name,x_0s,pars,simpars,u_star_tau,w_s,i_ppp,extras);
pars(i_ppp(:,1)) = U; # Put final value of U into the parameter vector
else
Error = [];
endif
opt_time = time-tick;
printf("Optimisation %i took %i iterations and %2.2f sec\n", i, \
length(Error), opt_time);
## Generate control
u_ol = u_star_t*U; # Not used - just for show
## Simulate system over one ol interval
[y_ol,ys_ol,x_ol] = eval(sprintf("%s_ssim(x_0s, pars, simpar, u_star_t);", s_system_name));
x_0 = x_ol(n_t+1,:)'; # Extract state for next time
y_ol = y_ol(1:n_t,:); # Avoid extra points due to rounding error
x_ol = x_ol(1:n_t,:); # Avoid extra points due to rounding error
y = [y; y_ol];
x = [x; x_ol];
u = [u; u_ol];
UU = [UU; U'];
UU_l = [UU_l; U_l'];
UU_c = [UU_c; U_c'];
t = [t; t_ol+t_last*ones(n_t,1) ];
t_last = t_last + T_ol;
endfor
endfunction
|
|
<
|
>
|
>
>
>
>
>
>
>
>
>
>
>
|
|
>
>
<
|
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
|
tick = time;
if extras.max_iterations>0
[U, U_all, Error, Y] = ppp_nlin(system_name,x_0s,pars,simpars,u_star_tau,w_s,i_ppp,extras);
pars(i_ppp(:,1)) = U; # Put final value of U into the parameter vector
else
Error = [];
endif
ppp_time = time-tick;
t_ppp = [t_ppp;ppp_time];
## Generate control
u_ol = u_star_t*U; # Not used - just for show
## Simulate system over one ol interval
par(i_ppp(:,3)) = pars(i_ppp(:,1)); # Update the simulation ppp weights
[y_ol,x_ol] = eval(sprintf("%s_sim(x_0, par, simpar, u_star_t);", system_name));
## Tune parameters/states
if (extras.estimate==1)
tick = time;
par_est = pars(i_par(:,1));
p = [p; par_est'];
pars = ppp_optimise(s_system_name,x_0s,pars,simpar,u_star_t,y_ol,i_par,extras);
est_time = time-tick;
t_est = [t_est;est_time];
endif
x_0 = x_ol(n_t+1,:)'; # Extract state for next time
y_ol = y_ol(1:n_t,:); # Avoid extra points due to rounding error
x_ol = x_ol(1:n_t,:); # Avoid extra points due to rounding error
y = [y; y_ol];
x = [x; x_ol];
u = [u; u_ol];
UU = [UU; U'];
UU_l = [UU_l; U_l'];
UU_c = [UU_c; U_c'];
t = [t; t_ol+t_last*ones(n_t,1) ];
t_last = t_last + T_ol;
endfor
endfunction
|