PROPT Batch Fermentor
|
|
This page is part of the PROPT Manual. See PROPT Manual. |
Dynamic optimization of bioprocesses: efficient and robust numerical strategies 2003, Julio R. Banga, Eva Balsa-Cantro, Carmen G. Moles and Antonio A. Alonso
Case Study I: Optimal Control of a Fed-Batch Fermentor for Penicillin Production
Problem description
This problem considers a fed-batch reactor for the production of penicillin, as studied by Cuthrell and Biegler (1989). This problem has also been studied by many other authors (Dadebo & McAuley 1995, Banga & Seider 1996, Banga et al. 1997). We consider here the free terminal time version where the objective is to maximize the amount of penicillin using the feed rate as the control variable. It should be noted that the resulting NLP problem (after using CVP) does not seem to be multimodal, but it has been reported that local gradient methods do experience convergence problems if initialized with far-from-optimum profiles, or when a very refined solution is sought. Thus, this example will be excellent in order to illustrate the better robustness and efficiency of the alternative stochastic and hybrid approaches. The mathematical statement of the free terminal time problem is:
Find u(t) and t_f over t in [t0; t_f ] to maximize
subject to:
where x1, x2, and x3 are the biomass, penicillin and substrate concentrations (g=L), and x4 is the volume (L). The initial conditions are:
![{\displaystyle x(t_{0})=[1.5\ 0\ 0\ 7]'}](https://wikimedia.org/api/rest_v1/media/math/render/svg/209176216daee2b2bb340225ef66e004b288a791)
There are several path constraints (upper and lower bounds) for state variables (case III of Cuthrell and Biegler, 1989):
The upper and lower bounds on the only control variable (feed rate of substrate) are:
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.Solving the problem on multiple grids.
The problem is solved in two stages. First, a solution is computed for a small number of collocation points, then the number of collocation points is increased, and the problem is resolved. This saves time, compared to using the fine grid immediately.
toms t
toms t_f
nvec = [35 70 80 90 100];
for i=1:length(nvec)
n = nvec(i);
p = tomPhase('p', t, 0, t_f, n);
setPhase(p);
tomStates x1 x2 x3 x4
tomControls u
% Initial guess
% Note: The guess for t_f must appear in the list before
% expression involving t.
if i==1
x0 = {t_f == 126
icollocate(x1 == 1.5)
icollocate(x2 == 0)
icollocate(x3 == 0)
icollocate(x4 == 7)
collocate(u==11.25)};
else
% Copy the solution into the starting guess
x0 = {t_f == tf_init
icollocate(x1 == x1_init)
icollocate(x2 == x2_init)
icollocate(x3 == x3_init)
icollocate(x4 == x4_init)
collocate(u == u_init)};
end
% Box constraints
% Setting the lower limit for t, x1 and x4 to slightly more than zero
% ensures that division by zero is avoided during the optimization
% process.
cbox = {1 <= t_f <= 256
1e-8 <= mcollocate(x1) <= 40
0 <= mcollocate(x2) <= 50
0 <= mcollocate(x3) <= 25
1 <= mcollocate(x4) <= 10
0 <= collocate(u) <= 50};
% Various constants and expressions
h1 = 0.11*(x3./(0.006*x1+x3));
h2 = 0.0055*(x3./(0.0001+x3.*(1+10*x3)));
% Boundary constraints
cinit = initial({x1 == 1.5; x2 == 0
x3 == 0; x4 == 7});
% This final condition is not necesary, but helps convergence speed.
cfinal = final(h2.*x1-0.01*x2) == 0;
% ODEs and path constraints
ceq = collocate({
dot(x1) == h1.*x1-u.*(x1./500./x4)
dot(x2) == h2.*x1-0.01*x2-u.*(x2./500./x4)
dot(x3) == -h1.*x1/0.47-h2.*x1/1.2-x1.*...
(0.029*x3./(0.0001+x3))+u./x4.*(1-x3/500)
dot(x4) == u/500});
% Objective
objective = -final(x2)*final(x4);
options = struct;
options.name = 'Batch Fermentor';
%options.scale = 'auto';
%if i==1
% options.solver = 'multiMin';
% options.xInit = 20;
%end
solution = ezsolve(objective, {cbox, cinit, cfinal, ceq}, x0, options);Problem type appears to be: qpcon
Time for symbolic processing: 0.60504 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Batch Fermentor f_k -87.746072952326273000
sum(|constr|) 0.000000000153458576
f(x_k) + sum(|constr|) -87.746072952172810000
f(x_0) 0.000000000000000000
Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied
FuncEv 1 ConstrEv 1634 ConJacEv 1634 Iter 398 MinorIter 7624
CPU time: 4.180827 sec. Elapsed time: 4.193000 sec.
Problem type appears to be: qpcon
Time for symbolic processing: 0.58293 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Batch Fermentor f_k -87.965550967328838000
sum(|constr|) 0.000000635729659157
f(x_k) + sum(|constr|) -87.965550331599175000
f(x_0) -87.746072952325463000
Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied
FuncEv 1 ConstrEv 98 ConJacEv 98 Iter 58 MinorIter 1461
CPU time: 2.496016 sec. Elapsed time: 2.474000 sec.
Problem type appears to be: qpcon
Time for symbolic processing: 0.58713 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Batch Fermentor f_k -87.997927563703584000
sum(|constr|) 0.000000125045315595
f(x_k) + sum(|constr|) -87.997927438658266000
f(x_0) -87.966078734278753000
Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied
FuncEv 1 ConstrEv 118 ConJacEv 118 Iter 91 MinorIter 1433
CPU time: 5.350834 sec. Elapsed time: 5.354000 sec.
Problem type appears to be: qpcon
Time for symbolic processing: 0.58916 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Batch Fermentor f_k -87.808944255539416000
sum(|constr|) 0.468138598171738750
f(x_k) + sum(|constr|) -87.340805657367682000
f(x_0) -87.998155723318973000
Solver: snopt. EXIT=1. INFORM=31.
SNOPT 7.2-5 NLP code
Iteration limit reached
FuncEv 1 ConstrEv 14 ConJacEv 14 Iter 11 MinorIter 10000
CPU time: 8.299253 sec. Elapsed time: 8.336000 sec.
Warning: Solver returned ExitFlag = 1
The returned solution may be incorrect.
Problem type appears to be: qpcon
Time for symbolic processing: 0.59682 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Batch Fermentor f_k -88.039618526588725000
sum(|constr|) 0.000001401154355635
f(x_k) + sum(|constr|) -88.039617125434376000
f(x_0) -87.808944255540538000
Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied
FuncEv 1 ConstrEv 231 ConJacEv 231 Iter 148 MinorIter 5151
CPU time: 51.917133 sec. Elapsed time: 20.128000 sec.
Plot result
subplot(2,1,1);
ezplot([x1; x2; x3; x4]);
legend('x1','x2','x3','x4');
title('Batch Fermentor state variables');
subplot(2,1,2);
ezplot(u);
legend('u');
title('Batch Fermentor control');
drawnow
% Copy solution for initializing next round
x1_init = subs(x1,solution);
x2_init = subs(x2,solution);
x3_init = subs(x3,solution);
x4_init = subs(x4,solution);
u_init = subs(u,solution);
tf_init = subs(t_f,solution);
end