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
|
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
|
-
+
+
+
-
-
+
+
+
|
simpar_model = simpar; # Initial internal model simulation params
simpar_pred = simpar; # Initial prediction simulation params
T_ol_0 = simpar.last; # The initial specified interval
n_t = round(simpar.last/simpar.dt); # Corresponding length
x_0 = eval(sprintf("%s_state(par);", system_name));
x_0_model = x_0;
[n_x,n_y,n_u] = eval(sprintf("%s_def;", system_name));
[n_par,m_par] = size(i_par);
if nargin<8
Q = ones(n_y,1);
endif
## Sensitivity system details -- defines moving horizon simulation
simpars = eval(sprintf("%s_simpar;", s_system_name));
pars = eval(sprintf("%s_numpar;", s_system_name));
x_0s = eval(sprintf("%s_state(pars);", s_system_name));
x_0_models = x_0s;
p = []; # Initialise saved parameters
## Times
## -- within opt horizon
n_Tau = round(simpars.last/simpars.dt);
dtau = simpars.dt;
Tau = [0:n_Tau-1]'*dtau;
[n_tau,n_w] = size(w_s);
tau = Tau(n_Tau-n_tau+1:n_Tau);
w = w_s(n_tau,:); # Final value of setpoint
## Main simulation loop
y = [];
x = [];
u = [];
t = [];
p = [];
t_last = 0;
UU = [];
UU_l =[];
UU_c =[];
t_ppp = [];
t_est = [];
its_ppp = [];
its_est = [];
t_open = [];
x_nexts = zeros(2*n_x,1);
## Initial U is zero
[n_U,junk] = size(A_u);
U = zeros(n_U,1);
## Initialise saved U
UU = [];
## Create input basis functions
u_star_tau = ppp_ustar(A_u,1,Tau',0,0,n_u-n_U);
## Reverse time to get "previous" U
U_old = U;
|
121
122
123
124
125
126
127
128
129
130
131
132
133
134
|
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
|
+
+
+
+
+
+
+
+
|
A_u_sim = A_u;
simpar_sim = simpar;
T_total = simpar.last;
endif
for i = 0:N_ol # Main loop
printf("%i",i);
UU = [UU; U']; # Save control U
if n_par>0
par_est = pars(i_par(:,1));
p = [p; par_est']; # Save up the estimated parameters
endif
if (extras.simulate==1)
[y_ol,u_ol] = ppp_RT_sim(U); # Simulate
else
[y_ol,u_ol] = ppp_RT(U); # Real thing
endif
t_start = time; # Initialise time
|
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
|
154
155
156
157
158
159
160
161
162
163
164
165
166
167
|
-
-
-
|
t_open = [t_open;T_ol];
## Generate input to actual system
u_star_t = ppp_ustar(A_u,1,t_ol',0,0,n_u-n_U);
## Tune parameters/states
if (estimating_parameters==1)
## Save up the estimated parameters
par_est = pars(i_par(:,1));
p = [p; par_est'];
## Set up according to interval length
if (T_ol>T_ol_0) ## Truncate data
simpar_est.last = T_ol_0;
y_est = y_ol(1:n_t+1,:);
else
simpar_est.last = T_ol;
|
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
|
249
250
251
252
253
254
255
256
257
258
259
260
261
262
|
-
|
pars(i_ppp(:,1)) = U; # Put final value of U into the parameter vector
## Save up data
y_ol = y_ol(1:n_ol,:);
y = [y; y_ol];
u = [u; u_ol];
UU = [UU; U'];
t = [t; t_ol+t_last*ones(n_ol,1) ];
t_last = t_last + T_ol;
endif
if (extras.simulate==1) # Do the actual simulation
|