Ind. Eng. Chem. Res. 2000, 39, 4287-4295
4287
Efficient Optimal Control of Bioprocesses Using Second-Order Information Eva Balsa-Canto and Julio R. Banga* Chemical Engineering Laboratory, Instituto de Investigaciones Marinas (CSIC), Eduardo Cabello 6, 36208 Vigo, Spain
Antonio A. Alonso Department of Chemical Engineering, Universidad de Vigo, Aptdo. 874, 36200 Vigo, Spain
Vassilios S. Vassiliadis Department of Chemical Engineering, Cambridge University, Pembroke Street, Cambridge CB2 3RA, United Kingdom
The dynamic optimization (open-loop optimal control) of bioprocesses is considered. It is shown how these problems can be solved using a recently developed method, based on the control vector parametrization concept, which makes use of second-order sensitivities to obtain exact gradients and Hessians for the objective function of the underlying dynamic process model. A further extension of this scheme, which makes use of restricted second-order information, is also presented. This extension results in an efficient way to solve general dynamic optimization problems, even for high levels of control discretization. This new approach allows the efficient and robust solution of two challenging case studies regarding the optimal control of fed-batch bioreactors taken from the open literature. 1. Introduction Most processes in the biotechnology and pharmaceutical industries operate in batch or in semicontinuous mode. To increase the productivity of these bioprocesses, many efforts have been devoted to their optimization and control. In this way, mathematical modeling and analysis have become fundamental tools in optimally designing and operating production facilities in the biotechnology industry.1,2 A class of problems that has received major attention is the dynamic optimization (open-loop optimal control) of fed-batch bioreactors. The use of fed-batch fermentation is very effective in overcoming undesired effects such as substrate inhibition and catabolite repression. An illustrative example of the importance of fed-batch culture is the large-scale production of monoclonal antibodies, with product revenues in the United States of $600 million in 1992 and an estimate of $4 billion by 1998.3 Many efforts have also been devoted to the related problems of dynamic model development and parameter estimation,4-6 fermentation process control,7 and robust operation.8,9 Recently, Conejeros and Vassiliadis10,11 have used a dynamic optimization framework to identify reaction bottlenecks for dynamic fermentation processes. This allowed a direct identification of enzymatic reactions for possible genetic manipulation by considering the impact of process time and operating conditions simultaneously. The dynamic optimization of fed-batch bioreactors is a challenging task, as most bioprocesses have highly nonlinear dynamics, and constraints are also frequently present on both the state and the control variables. Therefore, efficient and robust optimal control methods * To whom correspondence should be addressed. E-mail:
[email protected]. Phone: +34 986214473. Fax: +34 986292762.
are needed to obtain the optimal operating policies. Numerical methods for the solution of optimal control problems (OCP) are usually classified under two categories: indirect and direct approaches. Indirect (classical) approaches are based on the transformation of the original optimal control problem into a two-point boundary value problem using the necessary conditions of Pontryagin.12 Although several research studies have followed this approach,13-18 the resulting boundary value problems are usually very difficult to solve, especially when state constraints are present. Alternatively, direct approaches transform the OCP into a nonlinear programming (NLP) problem, either using control vector parametrization (CVP)19,20 or complete (control and state) parametrization.21 Luus22 proposed the use of iterative dynamic programming (IDP), but this method seems to be computationally expensive, especially for systems involving a large number of differential and algebraic equations (DAEs). Further, several search parameters of IDP must be adjusted by the user, so a number of exploratory runs are necessary. It should be noted that the NLPs arising from direct approaches are sometimes multimodal, so deterministic (gradient-based) optimization techniques may converge to local optima, especially if they are started far away from the global solution. Adaptive stochastic methods have been suggested as an alternative to surmount these difficulties.23 However, refined solutions are usually obtained at a large computational cost using these methods. Although there is always a tradeoff between convergence speed and robustness in both stochastic and deterministic methods, the latter usually have the opposite behavior; i.e., they converge very fast if they are started close to the global solution. Clearly, a convenient approach would be to combine both meth-
10.1021/ie990658p CCC: $19.00 © 2000 American Chemical Society Published on Web 09/29/2000
4288
Ind. Eng. Chem. Res., Vol. 39, No. 11, 2000
odologies to compensate for their weaknesses while enhancing their strengths. A hybrid (stochastic-deterministic) approach was considered by Balsa-Canto et al.24 with good results. This approach was developed by adequately combining the key elements of a stochastic and a deterministic method, taking advantage of their complementary features. Here, we present the application of a recently developed deterministic method based on the CVP concept, which makes use of second-order sensitivities to obtain exact gradients and Hessians. A further extension of this scheme, which makes use of restricted second-order information, is also presented. This extension results in an efficient way to solve general OCPs, even for high levels of control discretization. This is especially interesting in this context because most optimal control policies of bioreactors present a quite wild shape, with many abrupt changes. We show how this new approach allows the efficient solution of two challenging bioreactor case studies. In the case of multimodal problems, this new approach can be used to significantly improve the efficiency and robustness of stochastic-deterministic hybrids. Furthermore, because of its efficiency, this methodology can be extended and implemented on-line in closed loop, that is, in a nonlinear model predictive control (MPC) framework. 2. Statement of Optimal Control Problem Most dynamic optimization problems related to bioprocesses can be stated in the following general form, considering a fixed terminal time tf: Find the control vector u(t) over t ∈ [t0,tf] to minimize (or maximize) a performance index J(x,u),
J(x,u) ) Θ(x(tf))
(3)
An additional set of algebraic inequality constraints are the upper and lower bounds on the state and control variables, eqs 4 and 5:
xL e x(t) e xU
(4)
uL e u(t) e uU
(5)
3. Solution of Optimal Control Problems Using Second-Order Sensitivities The original system dynamics for the problem considered are the following,
f(x3 ,x,u,v)) 0n,
t ∈ [t0,tf]
(6)
with the following set of initial conditions,
x(t0)) x0(v)
(7)
(8)
which is common in control vector parametrization (CVP) approaches (a review of these may be found in Vassiliadis25). For example, in the case of a simple piecewise constant approximation defined by a uniform mesh of F time elements, v would be a vector with components the control values at each time element plus the time-invariant parameters originally present in the problem. Thus, its dimension will be F˜ ) σF + η. In Vassiliadis,25 the use of the first-order sensitivity equations was proposed as a means of obtaining gradients with respect to optimization parameters in conjuction with a CVP approach. These can be obtained by a chain rule differentiation applied to the system in eqs 6-8 with respect to the time-invariant parameters v,
∂f ∂x ∂f ∂u ∂f ∂f ∂x3 + + + ) 0n×F˜ , t ∈[t0,tf] ∂x3 ∂v ∂x ∂v ∂u ∂v ∂v
(9)
with the initial conditions
∂x0 ∂x (t ) ) (v) ∂v 0 ∂v
(10)
The system Jacobian in eq 9 is the partitioned matrix,
[∂x∂f3 , ∂x∂f ]
(2)
where x is the vector of state variables, with initial conditions x(t0) ) x0, and also subject to sets of inequality constraints (path, interior point, or end point, depending on the time values at which they are active), eq 3:
g(x(t),u(t)) e 0
u ) u(v)
(1)
subject to a set of differential-algebraic equality constraints, eq 2,
f(x3 ,x,u,v) ) 0
where x, x3 ∈Rn are the state (output) variables and their time derivatives, respectively, u ∈Rσ are control (input) variables, and v ∈Rη are time-invariant parameters. The functions f are such that f: Rn × Rn × Rσ × Rη f Rn. u indicates the control variable vector and v the timeindependent decision parameters. Depending on the implementation, the control variables may be approximated by some type of discretization (e.g., piecewise constant), which may involve some (or all) of the parameters in set v,
(11)
which is the same as the one in the resulting equations used in the integration of eq 6, as for example noted in Leis and Kramer.26 This allows efficient integration of eq 9 in a decoupled direct way, simultaneously with that of eq 6. A second differentiation (as presented by Vassiliadis et al.27) of eqs 9 and 10 with respect to v yields the following equation,
{[ {[
}
∂f ∂2x3 ∂f ∂2x X IF˜ X I + + F˜ ∂x3 ∂x ∂v2 ∂v2 ∂x3 T ∂2f ∂x3 ∂2f ∂x ∂2f ∂u ∂2f In X + + + + ∂v ∂x3 2 ∂v ∂x∂x3 ∂v ∂u∂x3 ∂v ∂v∂x3 ∂x T ∂2f ∂x3 ∂2f ∂x ∂2f ∂u ∂2f In X + + 2 + + ∂v ∂x3 ∂x ∂v ∂x ∂v ∂u∂x ∂v ∂v∂x
]
[
]
( ) ][ ] [ ( ) ][ ] ∂ f ∂x + [∂u∂f X I ] ∂∂vu + [I X (∂u∂v) ][∂x∂3 ∂uf ∂v∂x3 + ∂x∂u ∂v 2
F˜
T
2
n
2
2
]
∂ f ∂u ∂2f + + ∂u2 ∂v ∂v∂u ∂2f ∂x3 ∂2f ∂x ∂2f ∂u ∂2f + + + ∂x3 ∂v ∂v ∂x∂v ∂v ∂u∂v ∂v ∂v2
[
with initial conditions given by
2
]}
) 0n‚F˜ ×F˜ ,
t ∈[t0,tf] (12)
Ind. Eng. Chem. Res., Vol. 39, No. 11, 2000 4289
∂2x0 ∂2x (t ) ) (v) 0 ∂v2 ∂v2
(13)
The part of eq 12 in the second set of curly braces (“{}”) represents the right-hand side of that equation. It depends on second-derivative terms of the original system of equations (6) and it also involves the firstorder sensitivity terms. The expressions in (12) and in (13) comprise equations involving tensors of third order. Next, we consider the product of a time invariant vector p ∈ RG˜ with the second-order matrices (Hessians) for each state, that is, with ∂2xi/∂v2 for all t ∈[t0,tf]. The result of eq 12 is postmultiplied by p to obtain
{‚‚‚}‚p + {‚‚‚}‚p ) 0n‚F˜ ×1
(14)
By comparing terms the equivalent form is derived,
∂f T ∂f T Z4 + Z + A(x3 ,x,u,v)) 0n×F˜ ∂x3 ∂x
(15)
where Z, Z4 ∈ RF˜ ×n are the matrices whose columns are respectively given by the matrix vector products zi, z3 i ∈Rn, with
zi )
∂2xi ∂v2
p, z3 i )
∂2x˘ i p, ∂v2
i ) 1, 2, ..., n, ∀ t ∈ [t0,tf] (16)
( )
Finally, the set of initial conditions for these are
[ [
] ]
∂2x01 (v) ∂v2 · (Z(t0))T ) · 2 · ∂ x 0n pΤ (v) ∂v2 pΤ
(17)
It is clear that eq 15 is a first-order set of sensitivity equations, similar to eq 9. The Jacobian for these equations is the same as that in eq 11 used in the solution of the DAE system and the first-order sensitivity system. The only complication arises in the righthand side, which will be -A(x3 ,x,u,v); at first sight its calculation appears to require significant effort. Nonetheless, proper exploitation of sparsity and recurrent vector product terms can reduce this substantially, coupled with symbolic or automatic differentiation. Generation of the full second-order matrices for each xi, i ) 1, 2,..., n, can be viewed in terms of eq 15, in which the system is solved F˜ times, in each of which p ) ei is used, where ei are the unit vectors. Solution of the ith problem will yield as rows of ZT(t) the columns of ∂2xi/∂v2 (or its rows, as these matrices are symmetric) for all xi one column at a time. An interesting feature of the derivation in eq 15 is that for truncated Newton-type algorithms (e.g., Nash28), requiring only a product of the Hessian matrix with a search direction vector, the analysis presented provides a rigorous way to calculate these products exactly with comparable computational cost to that of the first-order sensitivity evaluation. The approach of using the Hessian search direction vector product evaluation scheme was tested with the TN (truncated Newton) optimization solver of Nash.28 This code was modified so that
the vector product can be obtained exactly as shown in eq 15. The TN approach uses an iterative scheme (conjugate-gradient-based) to solve the Newton equations of the optimization problem,
H(x)p ) -g(x)
(18)
where H(x) is the Hessian matrix, p is the search direction, and g(x) is the gradient vector. An iterative solution of the set of these equations requires the product γ ) H(x)p. In the case of a dynamic optimization problem, one clear way to obtain it is by using the exact product form in eq 15. Application of this scheme is an original contribution in this paper. To fully automate the application of this procedure, a Mathematica notebook was developed to derive symbolically the first-order sensitivities and the exact product form as presented in the equations (9) and (15). This notebook also builds automatically the source code (Fortran or C) of the resulting extended initial value problem (IVP) in a format ready to be used by standard solvers such as LSODA29 or DASSL. The modified TN solver is then used to solve the outer optimization problem, making use of the exact first- and (restricted) second-order information. In this way, optimal control problems can be solved very efficiently, even for high levels of control discretization. In the following, we will refer to this approach as one-run optimization. To achieve higher discretization values (F) with even more reduced computation times, we also considered a mesh-refining approach consisting of successive reoptimizations with increasing F values, where the initial and final control discretization values and integration tolerances are set by the user. The integration tolerances are taken at least as 2 orders of magnitude inferior to the optimization ones. The basic steps of this mesh-refining approach are outlined below: Step 1. Guess initial control profile (usually a flat control value u0(t) ) uguess) Step 2. Choose initial and final values for the discretization (Fi,Ff) and for the integration tolerances (int•toli,int•tolf). Optimization tolerances will be computed from the integration tolerances: opt•toli ) tol•ratio‚int•toli and opt•tolf ) tol•ratio‚int•tolf. Choose also a parameter rF used to increase the discretization level from one run to the next. Step 3. Compute the number of refining optimizations (NRO) and the tolerance-reduction step δtol using the expressions Ff/Fi ) rF(NRO-1) and δtol ) (1/NRO)(int•toli,int•tolf). Step 4. Initial optimization run (k ) 1): compute, using TN, the optimal control profile u1(t) corresponding to the initial control discretization Fi and the tolerances int•toli,opt•toli. Step 5. Optimization loop refining the control discretization, F, k)2 while (k
