|
3 | 3 | %
|
4 | 4 | % You'll learn:
|
5 | 5 | % +: how to solve optimization problems
|
6 |
| -% +: How to apply a model predictice control with non linear models |
| 6 | +% +: How to apply a model predictive control with nonlinear models |
7 | 7 | %
|
8 | 8 | %% The problem
|
9 | 9 | %
|
|
12 | 12 | % lb < u < ub
|
13 | 13 | % 0 < |du| < du_max
|
14 | 14 | %
|
| 15 | +% About the process: |
| 16 | +% CSTR with van de vusse reaction system |
| 17 | +% A -> B |
| 18 | +% B -> C |
| 19 | +% 2A -> D |
| 20 | +% |
15 | 21 | % ============================================================
|
16 | 22 | % Author: ataide@peq.coppe.ufrj.br
|
17 | 23 | % homepage: github.com/asanet
|
|
33 | 39 | % The objective function parameters
|
34 | 40 | Hp = 1000; % Prediction horizon
|
35 | 41 | nt = 50; % Time sampling points
|
36 |
| -P_sp = 1; % Set point violetion penalty |
| 42 | +P_sp = 1; % Set point violation penalty |
37 | 43 | P_u = 1e6; % Saturation penalty
|
38 | 44 | P_du = 1e6; % Delta u penalty
|
39 | 45 | lb = 273.15; % lower bound to u
|
|
60 | 66 | % Simulation span
|
61 | 67 | tspan = 1:3*Hp;
|
62 | 68 |
|
63 |
| -% The intial u profile (must have Hc values) |
| 69 | +% The initial u profile (must have Hc values) |
64 | 70 | u = [373.15 0 0 0]';
|
65 | 71 | umi = u(1);
|
66 | 72 |
|
|
109 | 115 | um(i+1) = u(1);
|
110 | 116 | umi = u(1);
|
111 | 117 |
|
112 |
| - % Disturbance in tal |
| 118 | + % Disturbance in tau |
113 | 119 | if i == fix(nsamples/2)
|
114 | 120 | Tf = 420;
|
115 | 121 | end
|
116 | 122 | end
|
117 | 123 |
|
118 |
| -%% plot the data |
| 124 | +% Plot the data |
119 | 125 | figure;
|
120 | 126 | subplot(211)
|
121 | 127 | h = plot(tm,ym(:,3),tm,spval*ones(nsamples,1),'LineWidth',1.5);
|
|
136 | 142 |
|
137 | 143 | function f = mpc_cost_function(u,ynow)
|
138 | 144 |
|
139 |
| - % Solve the model until the preditction horizon |
| 145 | + % Solve the model until the prediction horizon |
140 | 146 | ts1 = linspace(0,Hp,nt);
|
141 | 147 | [t,y] = ode15s(@model,ts1,ynow,simulationOpt,u);
|
142 | 148 |
|
|
151 | 157 | % Delta u penalty
|
152 | 158 | f3 = -min(0, du_max - max(abs([umi- u(1); u(2:end)])));
|
153 | 159 |
|
154 |
| - % THe objective function |
| 160 | + % The objective function |
155 | 161 | f = P_sp*f1 + P_u*f2 + P_du*f3;
|
156 | 162 | end
|
157 | 163 |
|
158 | 164 | function dy = model(t,y,u)
|
159 | 165 | % Van de vusse reaction in a CSTR
|
160 | 166 | Ca = y(1); Cb = y(2); T = y(3);
|
161 | 167 |
|
162 |
| - % the regularization function |
| 168 | + % The regularization function |
163 | 169 | reg = 1/2 + tanh(p*(t - td))/2;
|
164 | 170 |
|
165 |
| - % control variable |
| 171 | + % Control variable |
166 | 172 | Tc = sum(u.*reg);
|
167 | 173 |
|
168 |
| - % reaction rates |
| 174 | + % Reaction rates |
169 | 175 | r1 = k1*exp(-E1/R/T)*Ca;
|
170 | 176 | r2 = k2*exp(-E2/R/T)*Cb;
|
171 | 177 | r3 = k3*exp(-E3/R/T)*Ca^2;
|
172 | 178 |
|
173 |
| - % mass balances |
| 179 | + % Mass balances |
174 | 180 | dCa = (Caf - Ca)/tau -r1 -2*r3;
|
175 | 181 | dCb = -Cb/tau + r1 - r2;
|
176 | 182 | dT = (Tf - T)/tau -1/rho/cp*( UA/V*(T - Tc) + H1*r1 + H2*r2 +H3*r3 );
|
|
0 commit comments