2526
Ind. Eng. Chem. Res. 2007, 46, 2526-2534
Optimization of Metabolite Production in Fed-Batch Cultures: Use of Sufficiency and Characteristics of Singular Arc and Properties of Adjoint Vector in Numerical Computation Hwa Sung Shin and Henry C. Lim* Department of Chemical Engineering and Materials Science, UniVersity of California, IrVine, California 92697
Fed-batch fermentation for maximization of metabolite production at a final time is normally a singular control problem. The singular optimal feed rate strategy is dependent on the relative position of peaks in the specific rates of cell growth (µ), substrate consumption (σ), and product formation (π). Analyses and application of Pontryagin’s minimum principle, singular control theory, the generalized Legendre Clebsch conditions and the properties of adjoint variables that represent the gradients of the performance index, with respect to state variables λi(t) ) ∂J/∂xi, yield sufficient conditions for the existence of singular arcs and their characteristics for processes with various forms of µ, σ, and π. From the analysis, we can infer the optimal feed rate structure, the region of singular arc, the sign of adjoint variables, and the value of state variables on a singular arc, and we can apply them in numerical optimization techniques such as the multiple shooting method, in which proper selection of the initial values at each switching time is essential and critical for effective and efficient computation. Introduction Metabolite maximization in the fed-batch fermentation is one of the most important problems in biochemical engineering and, thus, has received much attention.1-3 Typically, a manipulated variable is the feed rate of limiting nutrients. According to Pontryagin’s minimum principle (PMP),4 the problem is singular, because the manipulated variable appears linearly in the mass-balance equations and, therefore, in the Hamiltonian function. The singular control deduced by PMP is limited to a region called singular arc, and the optimal feed rate is a concatenation of singular, maximum, and minimum feed rates. The optimization problems are reduced to two-point boundary value problems via PMP. Usually, a multiple shooting method has been used to obtain numerical solutions.5 However, numerical methods based on PMP may give nonunique solutions and do not distinguish between maximum and minimum trajectories, especially in singular control problems.6 For this reason, identification of a singular region or the absence of a singular region and the characteristics of singular arc can simplify the optimization of fed-batch processes considerably. Unstructured and nonsegregational models with specific rates of cell growth rate (µ), metabolite production rate (π), and substrate consumption rate (σ) have been used to model fedbatch fermentation. The specific ratessµ, σ, and πsare allowed to be dependent on the limiting nutrient concentration, normally, the substrate concentration. However, the specific rates for morecomplex models are known to be dependent not only on the substrate concentration but also on metabolite concentration.7-9 In rare cases, the specific growth rate is known to be even dependent on cell mass concentration. Because of the complexity of a general model, sufficient conditions for the existence of a singular arc or the analysis of singular arcs have not been reported. Thus, former numerical optimizations had assumed the existence of a singular arc and incorporated it into the numerical schemes. In this paper, the conditions for the existence of singular arcs are developed through PMP, generalized * To whom correspondence should be addressed. Tel.: 949 8249169. Fax: 949 824-2541. E-mail:
[email protected].
Legendre Clebsch conditions, singular control theory, and the properties of adjoint variables. Therefore, we can a priori determine the presence or absence of a singular arc using the sufficient conditions. The treatment is limited to the determination of an optimal feed rate profile. Monitoring and control of the fed-batch culture are other aspects of fed-batch culture optimization and are not covered in this paper. Problem Formulation The performance index to be maximized is the amount of product at the final time; it is expressed as
MinF(t)[PI ) -P(tf)V(tf)] The dynamic behavior of the system is represented by the cell mass, substrate, product, and overall mass balances; this is a system of ordinary differential equations (ODEs), as given below.
[] [ ] []
x˘ 1 µx1 0 x˘ 2 S -σx1 + F F, x˘ 3 ) πx1 0 x˘ 4 1 0
x3 ) a + bF
(1)
where x represents the total amounts of cell mass, substrate, and product and bioreactor volume, respectively; x1 ) XV, x2 ) SV, x3 ) PV, and x4 ) V. In most cases, the specific rates µ, σ, and π are functions of the substrate concentration alone or both the substrate and product concentrations. Therefore, we consider the specific rates to be arbitrary functions of the substrate and product concentrations; µ(S,P), σ(S,P), and π(S,P). The performance index to be minimized is determined by maximizing the total amount of metabolite produced [P(tf)V(tf)] at the final process time, tf:
Max [P(tf)V(tf)] ) Max [x3(tf)] ≡ Min [J ) -x3(tf)] (2) F(t)
F(t)
F(t)
According to Pontryagin’s Minimum Principle (PMP), this is equivalent to minimizing the Hamiltonian H, which is the linear
10.1021/ie061001o CCC: $37.00 © 2007 American Chemical Society Published on Web 03/23/2007
Ind. Eng. Chem. Res., Vol. 46, No. 8, 2007 2527
combination of the right-hand side of the mass-balance equations (eq 1) with the so-called adjoint variables (λ):
the singular feed rate, we take sequential time derivatives of the switching function until the feed rate F(t) appears explicitly:
H ) λTa + λTbF ) (λ1 µ - λ2σ + λ3π)x1 + (SF λ2 + λ4)F ) H1 + φF (3)
φ ) λTb ) SF λ2 + λ4 ) 0
where H1 denotes the terms that do not contain F and φ is the coefficient associated with F, the so-called switching function. PMP yields a system of ODEs for the adjoint variables, subject to the final conditions (transversality conditions) given below:
[] λ˙ 1
λ˙ 2 λ˙ 3
[
φ˙ ) -λTaxb ) -λTc ) -λT[c1,c2,c3]T ) -{(λ1 µs - λ2σs + λ3 πs)(SF - S) (λ1 µp - λ2σp + λ3 πp)P}
λ1 µ - λ2σ + λ3 π (λ1 µs - λ2σs + λ3 πs)x1/x4
) - (λ1 µp - λ2σp + λ3 πp)x1/x4 {(λ1 µs - λ2σs + λ3 πs)(-S) + (λ1 µp - λ2σp + λ3 πp)(-P)}x1/x4
[λ1(tf) λ2(tf) λ3(tf) λ4(tf) ] ) [0 0 -1 η ]
]
[{
) λ1
∂µ
∂ ln S˜
(5)
when φ ) SF λ2 + λ4 < 0
Fmin
when φ ) SF λ2 + λ4 > 0
} {
∂σ
∂ ln S˜
∂σ
+
∂ ln P
{
∂π
∂ ln S˜
}
+
+
∂π
}]
x1 ∂ ln P x4
()
where
S˜ = SF - S
(10a)
x1 [c1,c2,c3] ) [- µˆ , σˆ , - πˆ ] x4
(10b)
χˆ =
∂χ ∂ ln S˜
+
∂χ
, χ ) µ, σ, π
(10c)
∂ ln P
Aˆ = λ1 µˆ - λ2σˆ + λ3 πˆ
(10d)
Note that the newly defined gradients µˆ , σˆ , and πˆ are logarithmic gradients, which reduce to the usual gradient when the specific rates are functions of the substrate concentration only. Further differentiating eq 9 results in
F ) Fsingular when φ ) SF λ2 + λ4 ) 0 over finite time interval(s), tq < t < tq+1 when x4(t) ) Vmax
∂ ln P
- λ2
x1 x1 ≡ [λ1 µˆ - λ2σˆ + λ3 πˆ ] ) -λT[-µˆ , σˆ , -πˆ ]T ) x4 x4 x1 Aˆ ) 0 (9) x4
Because the Hamiltonian H is linear, with respect to the feed rate F, the Hamiltonian is minimized according to the sign of the switching function φ(t):
Fmax
∂µ
+
λ3
H (t) ) H (tf) ) λ3(tf)π(tf)x1(tf) ) -π(tf)x1(tf) < 0 (6)
0
x1 x4
x1 λ3{πs(SF - S) - πpP}] x4
(4)
where the subscripts s and p are used to denote partial derivatives µs ) ∂µ/∂S, σs ) ∂σ/∂S, πs ) ∂π/∂S, µp ) ∂µ/∂P, σp ) ∂σ/∂P, πp ) ∂π/∂P, and φ(tf) ) η > 0 because F ) Fmin ) 0 at the final time tf. The Hamiltonian H is constant over the entire time interval and F ) 0 at the final time. Therefore, the constant value of the Hamiltonian obtained from eqs 3 and 5 is negative.
{
()
) -[λ1{µs(SF - S) - µpP} - λ2{σs(SF - S) - σpP} +
∂H )∂x
λ˙ 4
(8)
(7)
The optimal feed rate profile is any concatenation of Fmax, Fmin, and Fsingular. When the initial conditions (X0, S0, P0, V0) lie on the singular hyperspace, it is possible for the optimal feed rate to be singular from the start until the reactor volume is full. It is also possible, if Fmax is not sufficiently large, that the singular feed rate can reach the upper limit and, therefore, the singular feed rate must be followed by Fmax before the reactor volume is full. If the chosen final time tf is very small, then it is possible for the feed rate to be the maximum to fill the reactor volume completely within the given final time. Thus, the shape of the concatenation of feed rates is affected by upper and lower limits on the feed rate, the initial values, and the final time. Singular Feed Rate On the singular arc, φ(t) )0 over the finite time interval and, therefore, its higher derivatives must also vanish. Thus, to obtain
φ¨ ) λT[(axc - cxa) - cxbF] ) λT[e - hF] ) λTg ) [(λ˙ 1 µˆ - λ˙ 2σˆ + λ˙ 3 πˆ ) + (λ1 µˆ s - λ2σˆ s + λ3 πˆ s)S˙ + (λ1 µˆ p - λ2σˆ p + λ3 πˆ p)P˙ ]X ) [-A µˆ + (Asσˆ X + Aˆ s S˙ ) + (-Ap πˆ X + Aˆ pP˙ )]X
[
) -A µˆ + (Asσˆ - Aˆ sσ)X + (Aˆ p π - Ap πˆ )X + (Aˆ s(SF - S) - Aˆ pP) w -Aµˆ +
[( )
( )]
Ap ∧ 2 As Yp/s π σ
∧
σ2X - Aˆˆ
()
( )]
F X (11) x4
F ) x4
[λ1, λ2, λ3][g1, g2, g3]T ) 0 Because of the fact that the feed rate appears in eq 11, no higher-
2528
Ind. Eng. Chem. Res., Vol. 46, No. 8, 2007
order derivative is required to be taken, and the singular feed rate is obtained from eq 11:
-Aµˆ + Fsin )
[( )
( )]
Ap ∧ 2 As Yp/s π σ
()
∧
σ2X
Aˆˆ x4
(13)
The Generalized Legendre Clebsch (GLC) condition, obtained from eq 11, is
(-1)
( )
Aˆˆ ∂ d2φ ) g 0 w (λ1 µˆˆ - λ2σˆˆ + λ3 πˆˆ ) g 0 ∂F dt2 x4 (14)
Equation 14 being an inequality, there are two cases: one being equal to zero and the other being greater than zero. The equality case trivially satisfies the GLC condition. The equation λTh ) 0 leads to λTe ) 0 in eq 11. These two relationships and eq 9 then can be rearranged as follows:
[ ][ ] e1 e2
e3
µˆ -σˆ πˆ
µˆˆ -σˆˆ πˆˆ
λ2 ) Mλ ) 0 λ3
(15)
e1(σˆˆ πˆ - σˆ πˆˆ ) - e2( µˆ πˆˆ - µˆˆ πˆ ) + e3( µˆˆ σˆ - µˆ σˆˆ ) ) K(x1,x2,x3,x4) ) 0 (16) K is a singular arc, which is a function of the state variables over the singular interval, and the corresponding feed rate is obtained by differentiating it with respect to time:
( )( ) ( )( ) ( )( ) ( )( ) ∂K dx1 ∂K dx2 ∂K dx3 ∂K dx4 + + + ∂x1 dt ∂x2 dt ∂x3 dt ∂x4 dt
) (Kx1 µ - Kx2σ + Kx3 π)x1 + (Kx2SF + Kx4)F ) 0 wF)
(17)
-(Kx1 µ - Kx2σ + Kx3 π)x1 Kx2SF + Kx4
where the subscripts refer to partial derivatives, Kxj = ∂K/∂xj. Thus, it is clear that the equality λ1 µˆˆ - λ2σˆˆ + λ3 πˆˆ ) 0 leads to a singular arc and a singular feed rate. However, this is a special case. Therefore, we can generally assume the case of inequality
λ1 µˆˆ - λ2σˆˆ + λ3 πˆˆ > 0
(18)
Existence of Singular Arcs The optimal feed rate structure of this problem may have degenerate cases for certain initial conditions, functional forms of µ, σ, and π, and constraints (such as feed rate limits). However, we assume that the feed rate limit is large, so that the singular feed rate does not exceed the limit, Fmax.
(19)
Generally, it is not possible to assess the time behavior of adjoint variables a priori. However, such an assessment can lead to a very efficient way to predict the optimal feed rate strategy, because the adjoint variables represent the gradients of the performance index PI, with respect to state variables,10 λi(t) ) ∂J/∂xi. However, despite a vast amount of literature on fedbatch optimization, this fact has not been utilized. We will utilize this information to elucidate the existence of singular arcs by investigating behaviors of adjoint variables over a singular region. From eqs 9 and 4, over the singular arc, φ˙ ) λ˙ 2(SF - S) - λ˙ 3P ) 0, but SF - S > 0 and P > 0. Therefore,
sign(λ˙ 2) ) sign(λ˙ 3)
(20)
Equations 9 and 11 are both linear in three adjoint variables. Therefore, these two adjoint variables can be expressed linearly in terms of the one remaining adjoint variable. For example,
[][ ][ ] λ2 c2 c3 λ3 ) g2 g3
λ1
To avoid the triviality of the adjoint variables of eq 15, the determination of M must be zero: det(M) ) 0. This gives
K˙ )
λ1(t) < 0
(12)
On the singular arc (φ ) 0), eqs 3 and 6 imply
λ1 µ - λ2σ + λ3 π < 0
First, we investigate the sign of λ1. When the feed rate is either singular (φ ) 0) or zero (Fmin ) 0), eqs 4 and 13 show that λ˙ 1(t) > 0. However, λ1(tf) ) 0 at the final time. Therefore, λ1 is negative on the singular arc:
-1
-c1 -g1 λ1
(21)
Substituting λ2 and λ3 from eq 21 into eq 4 and integrating from the time of entry to the singular arc t s0 to any time on the singular arc ts yields
λ˙ 1 ) -f(x,F)λ1 w λ1(ts) ) λ1(t s0) exp(-
∫t t f dt) s
s 0
(22)
This implies that the sign of λ1 over the singular arc is determined by λ1 at the entry time to the singular arc, i.e., the sign of λ1 does not change over the singular arc. It is simple to show that other adjoint variables behave similarly, and, thus, we can conclude that all adjoint variables, λ1 through λ4, do not change in sign on the singular arc. Now, we must determine the definiteness of the signs. The adjoint variables are the gradient of a performance index,10 with respect to corresponding state variable, for example,
λ3(t) )
∂J ∂J ∂x1(tf) ∂J ∂x2(tf) ) + + ∂x3 ∂x1(tf) ∂x3(t) ∂x2(tf) ∂x3(t) ∂J ∂x4(tf) ∂J ∂x3(tf) + ∂x3(tf) ∂x3(t) ∂x4(tf) ∂x3(t)
)-
∂x3(tf) ∂x3(t)
+η
∂x4(tf) ∂x3(t)
)-
∂x3(tf) ∂x3(t)
(23)
where eq 2 and x4(tf) ) Vmax are used. At a time t s/ over the singular arc, we assume that an infinitesimal increase in x3(t s/) causes x3(tf) to decrease, then ∂x3(tf)/∂x3(t s/) < 0 and therefore, λ3(t s/) > 0. However, as previously concluded, the adjoint variables do not change in sign and therefore, λ3(ts) > 0 at all times ts over the singular arc. However, this implies that x3(tf) decreases, despite a concurrent infinitesimal increase in x3(ts) at all singular times. Therefore, the aforementioned assumption is incorrect, because the increase in the amount of metabolite x3(t s/) over the singular arc should result in increasing the amount of metabolite at the final time x3(tf). Consequently, the sign of λ3 is negative over the singular arc. Now, using eqs 8
Ind. Eng. Chem. Res., Vol. 46, No. 8, 2007 2529
and 9, the adjoint variables are expressed in terms of λ1 and λ3:
[
]
λ µˆ + λ3 πˆ λ µˆ + λ3 πˆ λ3 -SF 1 [λ1 λ2 λ3 λ4 ]) λ1 1 σˆ σˆ (24) Substituting the adjoint variables of eq 24 into the corresponding adjoint variables of eq 13, we obtain
{ ()
( )}
{ ()
( )}
σ2 1 µ∧ π∧ µ∧ π∧ λ1 + λ3 >0w λ1 + λ3 > 0 (25) σˆ σ σ σ′ σ σ Inequality 25 gives conditions that narrow the singular arc region considerably. However, the utility of inequality 25 is limited, because it contains λ1 and λ3, which are not known a priori. At least the knowledge of the signs of λ1 and λ3 would make it more useful. This will be done later. Inequality 25 is not the only sufficient condition for the existence of singular arcs. Other conditions for the existence of a singular arc can be deduced by considering the sign of λ2, which is not known a priori. We will investigate singular arcs, with respect to the sign of λ2. We explore the consequence of assuming λ2(t) is negative. With the assumption that λ2 is negative on the singular arc, we seek the conditions for which the singular arc exists. Because λ3 is negative, and because of eq 9, both λ1 and λ2 are negative on the singular arc:
λ1 )
λ2σˆ - λ3 πˆ 0 (28)
Criteria, in terms of specific rates only, independent of adjoint variables, can be readily applied to assess the regions in which the singular arc can exist and, therefore, are more desirable. These are obtained from inequalities 25 and 28. Let µˆ and σˆ have the same sign. If (µ/σˆ ) ) 0, the inequality µˆ (π/σˆ ) < 0 guarantees inequalities 25 and 28. However, if µˆ (µ/σˆ ) < 0, µˆ (π/σ/µ/σ)∧ < 0 must be satisfied and these two cases can be rewritten as follows:
(σµ)
∧
) 0, µˆ
∧
(πσ)
< 0, µˆ σˆ > 0
(29a)
or
π µˆ σ
∧
()
e µˆ
∧
(σµ) (πµ) < 0,
µˆ σˆ > 0
(29b)
Other conditions can be obtained by substituting inequalities 26 or 27 into inequality 18.
{ () ( )} ( ) { () ( ) } ( ) µˆ σˆ λ1 σˆ
∧
+ λ3
πˆ σˆ
∧
(λ1 µ + λ3 π)∧ σˆ
) σˆ
∧
> 0 (30)
(λ2σ - λ3 π)∧ µˆ ∧ ˆ 2 πˆ ∧ ˆ 2 1 λ2 (σ) + λ3 ( µ) ) (- µˆ ) µˆ σˆ µˆ µˆ
∧
>0 (31)
Once again, inequalities 30 and 31 contain adjoint variables, and therefore, their utilities are limited. Criteria in terms of specific rates only are obtained from inequalities 30 and 31. Following the procedure used to develop inequality 29, we obtain the criteria for the existence of a singular arc that satisfy inequalities 30 and 31:
() µˆ σˆ
∧
()
πˆ ) 0, σˆ σˆ
∧
< 0, µˆ σˆ > 0
(32a)
or
()
µˆ σˆ σˆ
∧
( ) ()( )
πˆ < 0, σˆ σˆ
∧
e σˆ
µˆ σˆ
∧
πˆ , µˆ σˆ > 0 µˆ
(32b)
In summary, the sufficient conditions for the existence of a singular arc are those satisfying inequalities 25, 27, 28, 30, and 31. Utilities of these conditions are limited, because they contain unknown adjoint variables; however, they may be used to ascertain computational results a posteriori. Sufficient conditions, in terms of specific rates only, are inequalities 29 and 32. These can be applied readily to determine a priori the admissibility of a singular feed rate. The sufficient conditions for the existence of a singular arc, above the inequalities, are classified according to the signs of µˆ , σˆ , and πˆ in Tables 1-3. Table 1 lists the conditions that are sufficient for the existence of a singular arc, in terms of specific rates. These conditions can be used to determine a priori the existence of a singular feed rate and incorporate it into a computational scheme to determine the optimal feed rate profile. With given specific rates and a known substrate and metabolite concentration at a specific time, we can determine what the feed rate should be at that time: bang-bang or singular. This information is extremely useful in a numerical scheme to obtain the optimal feed rate profile. Table 2 lists conditions that are sufficient for the nonexistence of a singular arc, in terms of specific rates. These conditions can be used to eliminate a priori the existence of a singular feed rate and incorporate it into a computational scheme to determine the optimal feed rate profile. Broader sufficient conditions that satisfy inequalities 25-28 are listed in Table 3. Most of these conditions, except cases 3 and 8, cannot be used to predict a singular feed rate a priori, because of unknown adjoint variables, however they may be used to ascertain computational results a posteriori, to prove optimality. Generally, an inclusion of a singular arc in the optimal feed rate structure is preferable for the performance index, and, therefore, an optimal feed rate should first force the substrate and metabolite concentration from the initial value to be on the singular region defined by the inequalities given in Table 3. So far, we have used the necessary conditions from PMP and GLC to identify regions where a singular arc exists for processes with various forms of µ, σ, and π. As a result, we have identified singular arc regions that are defined by these conditions. This is very significant, in that we could infer the optimal feed rate structure, the sign of adjoint variables, and the value of state variables on the singular arc and apply them
2530
Ind. Eng. Chem. Res., Vol. 46, No. 8, 2007
Table 1. Sufficient Conditions for the Existence of a Singular Arc, in Terms of Specific Rates σˆ
µˆ
πˆ
+
+
+
-
-
-
constraints
σˆ µˆ πˆ g g > 0, σ µ π 0>
πˆ µˆ σˆ g g π µ σ
σˆˆ µˆˆ πˆˆ g g σˆ µˆ πˆ πˆˆ µˆˆ σˆˆ g g πˆ µˆ σˆ
index 1
σˆˆ µˆˆ πˆˆ except ) ) σˆ µˆ πˆ
and
σˆ µˆ πˆ ) ) σ µ π
2
Table 2. Sufficient Conditions for Nonexistence of Singular Arc in Terms of Specific Rates σˆ
µˆ
πˆ
constraints
σˆ µˆ πˆ e e σ µ π
or
σˆ πˆ µˆ 0< e e σ π µ
or
πˆ µˆ σˆ e e µˆ /µ ) σˆ /σ, which leads to (π/µ)∧ > 0 and (πˆ /µˆ )∧ > 0. The middle panel of Figure 1 is a contour plot of (π/µ)∧. It shows that the substrate concentration should be S > 10 if a singular arc is to be a part of the optimal trajectory, i.e., (π/µ)∧ > 0. Because of the fact that the initial substrate concentration is 0.001 g/L, the first feed rate is Fmax, to increase the substrate concentration so that a singular feed rate is a part of the optimal feed rate profile. The substrate is transformed to an alcohol or cell mass and their rates of transformation are dictated by the
∂µ ∂S˜
(t1)(S˜ 2 - S˜ 1) +
(t1)S˜ 2 +
∂µ ∂P
∂µ ∂P
(t1)(P2 - P1) )
(t1)P2 - µˆ (t1) ≈ µˆ (t2) - µˆ (t1) (33)
where we assume
∂µ ∂µ ∂µ ∂µ (t1)S˜ 2 + (t1)P2 ≈ (t2)S˜ 2 + (t2)P2 ) µˆ (t2) ∂P ∂P ∂S˜ ∂S˜ because ∆t is infinitesimal. The same argument applies to (π/µ)∧. Consequently, there are linear relationships between ∆µ and ∆µˆ and between ∆(π/µ) and ∆(π/µ)∧:
∆µ ∝ ∆µˆ
(34a)
and
∆
( πµ) ∝ ∆( πµ)
∧
(34b)
This equation implies that a change in the logarithmic gradient of the specific growth rate ∆µ∧ is proportional to a change in the specific growth rate ∆µ, and a change in the logarithmic gradient of the ratio of specific rates ∆(π/µ)∧ is also proportional to a change in the ratio of specific rates ∆(π/µ). In other words,
Ind. Eng. Chem. Res., Vol. 46, No. 8, 2007 2531 Table 3. Conditions Imposed on the Singular Arc by Inequalities 25-28 σˆ
µˆ
πˆ
+ +
-
+
+
+
-
constraints
+ + + +
+ + + +
+ + + +
-
+ +
+ -
-
-
+
-
-
+
-
-
-
πˆ > λ1 > ∧ µˆ (µ/σ)
( )
(π/σ)∧ ∧
> λ1 ,
(µ/σ) ˆσ µˆ πˆ ) > >0 σ µ π
σˆ πˆ µˆ g > > 0, σ π µ µˆ σˆ πˆ > > > 0, µ σ π πˆ σˆ µˆ > > > 0, π σ µ
0> or
σˆ µˆ πˆ > g >0 σ µ π
()
0>
πˆ µˆ σˆ g > π µ σ
0>
µˆ πˆ σˆ > g , µ π σ
0>
πˆ σˆ µˆ > > , π σ µ
0>
µˆ σˆ πˆ > > , µ σ π
-
µ ∧ , σ
we can interpret the logarithmic gradient of specific rates in terms of plain specific rates. With these interpretations, we proceed to analyze the case study. The first feed rate Fmax brings the substrate concentration to a point that is favorable to µˆ but not favorable to (π/µ)∧ (see the contour plots in Figure 1). In other words, the substrate is consumed to favor, first, production of the cell mass and then, to maximize the alcohol production, the culture environment must be changed to favor alcohol production. Thus, the cell growth rate (µ) must not increase during the second feed period. From eq 34, this means that µˆ must not increase during the alcohol production stage and the arrow in the top panel of Figure 1 (µˆ ) illustrates this phenomena. Therefore, the second feed rate must be Fmin until the state moves onto a singular arc. The third feed rate, Fsingular, continues until the reactor volume is full, not increasing µˆ abruptly, as represented by the arrow in the middle panel of Figure 1 ((π/µ)∧), and then a batch operation, Fmin ) 0, occurs until the final time tf. The simulation results are consistent with the predictions. The optimal feed rate sequence is Fmax f Fmin f Fsingularf 0, which is consistent with the profile of the switching function, φ < 0 f φ > 0 f φ ) 0 f φ > 0. The substrate concentration increases and then decreases gradually as expected (see Figure 2a). Figure 2b shows the time variations of µˆ , σˆ , and πˆ , as well as those of µ/µˆ and π/πˆ . These variations are consistent with the existence conditions of a singular arc given in Table 3. These results are for MinF(t)[PI ) -P(tf)V(tf)]. This case study
3
∧
∧
2
2
∧
∧
or
2
2
∧
0 > λ2 >
4
2
∧
(µ/σ) >λ (π/σ µ/σ) (µ/σ) π/σ (µ/σ) , 0>λ >( ) µ/σ (µ/σ) π/σ (µ/σ) >λ , 0>( ) µ/σ (µ/σ) 0>
( )
0>
2
πˆ > λ1 µˆ
(π/σ) ∧ (µ/σ)2 0 > λ2 > , (µ/σ) (µ/σ)∧ πˆ µ∧ g 0, 0 > > λ1 σ µˆ
()
no 1
(π/σ)∧
2 π/σ ∧ (µ/σ) 0 > λ2 > , µ/σ (µ/σ)∧
0>
index
2
(πσ) π 0>( ) σ
0>
∧
∧
5
∧
6
(σµ) µ -λ ( ) σ - λ1
1
no 7
(π/σ)∧
πˆ > λ1 > µˆ (µ/σ)∧
8
( )
(π/σ) ∧ (µ/σ)2 , (µ/σ) (µ/σ)∧ πˆ µˆ σˆ 0> > ) π µ σ
(π/σ)∧ πˆ > λ1 > µ (µ/σ)∧
9
10 11
2 π/σ ∧ (µ/σ) > λ2 µ/σ (µ/σ)∧
( ) π/σ (µ/σ) π 0>λ >( ) , ( ) µ/σ (µ/σ) σ π/σ (µ/σ) π 0>( ) >λ , ( ) µ/σ (µ/σ) σ 0>
∧
2
∧
∧
∧
2
∧
2
∧
∧
2
∧
(σµ) µ -λ ( ) σ - λ1
∧
1
>0 >0
12
13
illustrates the importance of the newly defined special variations of specific rates µˆ , σˆ , πˆ , (π/µ)∧, and (πˆ /µˆ )∧, in regard to predicting the optimal feed rate strategy. It is well-known that the initial conditions greatly affect the performance index PI. Because the singular feed rate is superior to the bang-bang-only feed rates (Fmax f Fmin ) 0 or Fmin ) 0 f Fmax) for the performance index, if the initial conditions are chosen properly to be on the singular arc, the first bangbang feed rate is not necessary and, therefore, the performance index would improve. In this case, the sufficient conditions for singular arcs give an insight to the substrate concentration. In this case study, the initial substrate concentration, which is expected to have a value of >10 g/L, following the singular arc condition previously described, will increase and then decrease, because the process must follow the direction favorable to the cell mass and then the metabolite. A proper set of initial conditions that places the process on the singular arc is obtained through a parameter optimization of the two-point boundary value problem. Allowing the initial amount of substrate to be free, and keeping intact the other initial conditions (inoculum size and volume), [x1, x3, x4](t0) ) [0.1 g, 0 g, 1 L], the initial amount of substrate was determined to be 17.251 g. Thus, the new initial conditions necessary to make Fsingular f 0 the optimal feed rate sequence are [x1, x2, x3, x4](t0) ) [0.1 g, 17.251 g, 0 g, 1 L]. In other words, the performance index is minimized by proper selections of the optimal feed-rate profile F(t) and the initial amount of substrate S0V0 (or the initial substrate
2532
Ind. Eng. Chem. Res., Vol. 46, No. 8, 2007
Figure 1. Contour plots of µˆ , (π/µ)∧, and (πˆ /µˆ )∧ of the case study (ethanol production by S. cereVisiae), µ ) 0.408S/[(1 + P/16)(0.22 + S)], σ ) 10µ, and π ) S/[(1 + P/71.5)(0.44 + S)]. Dark arrow indicates the state variable direction with time.
concentration): MinF(t), V0S0[PI ) -P(tf)V(tf)]. The final amount of alcohol produced is 767.41 g, using Fsingular f Fmin ) 0, which is comparable to the value of 764.47 g, which is obtained with the initial conditions not on the singular arc: Fmax f Fmin f Fsingular f 0. Figure 3 compares the state variable profiles, the performance indices, and the feed rates for these two sets of initial conditions. As expected, through the analysis of the sufficient conditions for a singular arc, the optimum feed rate profile should place the culture onto the singular arc as soon as possible and apply the singular feed rate until the culture volume reaches the maximum. The first portion of the feed rate to reach the singular arc then should be bang-bang, the maximum feed rate followed by a batch period, Fmax f Fmin ) 0. The aforementioned results were obtained by picking the optimum initial amount of substrate while keeping the remainder intact. A further improvement can be obtained by picking the optimum initial volume, in addition to the optimum initial amount of substrate. As far as the inoculum is concerned, this being an autocatalytic reaction, the larger the better, only limited by a practical constraint in implementing the startup procedure. Therefore, we did not optimize the inoculum size and used a fixed value (0.1 g). By allowing the initial amount of substrate and volume to be free, we seek, for Fsingular, to satisfy the degrees of freedom associated with the two-point boundary value problem.12 In other words, we are searching the singular arc, with respect to both the optimum initial amount of substrate and the optimum initial volume. These values were determined to be 4.4768 g and 0.24396 L, respectively. Thus, the optimum initial conditions used are [x1, x2, x3, x4](t0) ) [0.1 g, 4.49 g, 0 g, 0.244 L]. In this case, the performance index is minimized by proper selections of the initial amount of substrate S0V0 and the initial volume, V0 and determining the singular feed rate:
Figure 2. Time profiles of state variables, switching function, and feed rate of the case study (ethanol production by S. cereVisiae), µ ) 0.408S/[(1 + P/16)(0.22 + S)], σ ) 10µ, and π ) S/[(1 + P/71.5)(0.44 + S)]. [x1, x2, x3, x4](t0) ) [0.1 g, 0.001 g, 0 g, 1 L], [SF, tf, Fmax, Fmin, Vmax ] ) [150 (g/L), 39.5 (h), 2 (L/h), 0 (L/h), 10 L]. Panel a shows state variables, feed rate, and switching function behaviors; panel b shows the characteristics of the singular arc.
Figure 3. Comparative time profiles of state variables, switching function, and feed rate of the case study (ethanol production by S. cereVisiae), M = Fmax, m = Fmin, S = Fsingular, MmSm = Fmax f Fmin f Fsingular f 0(Fmin ), Sm = Fsingular f Fmin, µ ) {0.408S}/[(1 + P/16)(0.22 + S)], σ ) 10µ, and π ) S/[(1 + P/71.5)(0.44 + S)]. [SF, tf, Fmax, Fmin, Vmax ] ) [150 (g/L), 39.5 (h), 2 (L/h), 0 (L/h), 10 L].
Min [PI ) -P(tf)V(tf)]
F(t),V0S0,V0
The amount of alcohol resulting from these initial conditions using Fsingular was 804.11 g, as compared to 767.41 g using
Ind. Eng. Chem. Res., Vol. 46, No. 8, 2007 2533
Fsingular f Fmin ) 0 and with 764.47 g obtained with the initial conditions not on the singular arc, Fmax f Fmin f Fsingular f 0. The substrate concentration behavior is similar to the case of Fsingular - Fmin, but the concentration increases higher and then decreases. The sufficient conditions of a singular arc can be used to obtain an excellent value of the substrate concentration at the start of the singular feed rate. This value can be most effectively used to provide the initial estimates necessary in a numerical scheme, such as the multiple shooting method. Discussion The maximization of metabolite at the final time is a class of very typical and important problems. This paper examines processes whose specific rates µ, σ, and π, are dependent on the concentrations of both the substrate and metabolite (S and P, respectively), not just S alone. The sufficient conditions for the existence of a singular arc are developed for processes with various forms of µ, σ, and π through the use of PMP, singular control theory, GLC conditions, and the properties of adjoint variables through the use of newly defined special variations µˆ , σˆ , πˆ , (π/µ)∧, and (πˆ /µˆ )∧. These conditions allow us to quickly determine if a singular control(s) is a part of the optimal feed rate sequence. A multiple shooting method is a well-accepted optimization technique for numerical solution of optimization problems but requires proper selection of the initial values at each switching time for effective and efficient computation. In this sense, information elucidated in this paper is very useful, in that we could not only provide proper initial values for the shooting method but also infer from it the optimal feed rate structure, the region of singular arc, the sign of adjoint variables, and the value of state variables on a singular arc, and then apply them in numerical optimization techniques such as the multiple shooting method. The aforementioned development is based on a class of fedbatch fermentation processes that can be described by four dynamic mass-balance equations with an arbitrary form of the specific rates of cell growth, substrate consumption, and product formation (µ(S,P), σ(S,P), and µ(S,P), respectively). In this development, it is assumed that the model is time-invariant; that is, the form and the values of the parameters in the model remain unchanged throughout the course of fed-batch operation. However, in practice, a model is, at best, an approximation of the complex biological reactions that occur in the cells, and, therefore, it may not accurately represent the culture during the entire period, which may include a growth phase and a production phase, which may be distinct in some organisms. Thus, there may be substantially variability from run to run or even within a run. To address this variability, it may be necessary to incorporate adaptive features into the scheme so that the variability within a run or from run to run can be addressed through an adaptive optimization scheme. Therefore, the ability to compute the optimum feed rate profile rapidly is essential in implementing the optimization scheme online or offline. When the singular feed rate can be obtained in a feedback mode (as a function of the state variables alone), the optimum feed rate can be implemented in a feedback control mode. This, in return, leads to the problem of monitoring fed-batch process variables (state variables), including the problem of measuring all state variables. Thus, there is a measurement problem that may require an estimation scheme for variables that are not readily measurable. In addition to determining the optimal feed rate profile, there is a vast area of monitoring and control
problems of the fed-batch process to implement adaptively the optimum feed rate through the course of one run or from run to run. This is beyond the scope of the present work and the interested reader may refer elsewhere13,14 for a detailed treatment of this topic. Acknowledgment One (HSS) of the authors was supported in part by a fellowship from Korea Science and Engineering Foundation (KOSEF) and a grant from the University of California, Irvine. Nomenclature F ) feed rate (cm3/s) Fmax ) maximum feed rate (cm3/s) Fmin ) minimum feed rate (cm3/s) Fsingular ) singular feed rate (cm3/s) H ) Hamiltonian function H1 ) part of Hamiltonian function k ) metabolite degradation rate P ) metabolite concentration (g/cm3) P0 ) initial metabolite concentration (g/cm3) S ) substrate concentration (g/cm3) SF ) feed substrate concentration (g/cm3) S0 ) initial substrate concentration (g/cm3) t ) time (s) tf ) final time (s) V ) bioreactor volume (cm3) Vmax ) maximum bioreactor volume (cm3) V0 ) initial volume (cm3) x1 ) total cell mass (g) x2 ) total substrate (g) x3 ) total metabolite (cm3) x4 ) bioreactor volume (cm3) X ) cell mass concentration (g/cm3) X0 ) initial cell mass concentration (g/cm3) Yx/s ) cell mass yield (g/g) Yp/s ) metabolite yield (g/g) Greek Symbols φ ) switching function µ ) specific growth rate (s-1) π ) specific metabolite production rate (s-1) σ ) specific substrate consumption rate (s-1) Accents and Superscripts ∧ ∧ ∧
) derivative with respect to ∂/∂(ln Sˆ ) + ∂/∂(ln P) ) second derivative with respect to ∂/∂(ln Sˆ ) + ∂/∂(ln P)
Literature Cited (1) Rani, K. Y.; Rao, V. S. R. Control of FermenterssA Review. Bioprocess. Eng. 1999, 21, 77-88. (2) Modak, J. M.; Lim, H. C.; Tayeb, Y. J. General Characteristics of Optimal Feed Rate Profiles for Various Fed-batch Fermentation Processes. Biotechnol. Bioeng. 1986, 28, 1396-1407. (3) Chevalot, I.; Marc, A. Interest of Fed-batch Culture for the Production of a Membrane Bound Protein by an Adherent Animal Cell. Biotechnol. Lett. 1993, 15 (8), 791-796. (4) . Pontryagin, L. S.; Boltyanskii, V. G.; Gamkrelidze, R. V.; Mishchenko, E. F. The Mathematical Theory of Optimal Processes; Interscience: New York, 1962. (5) Bryson, A. E.; Ho, Y. C. Applied Optimal Control; Wiley: New York, 1975. (6) Krener, A. J. The High Order Maximum Principle and its Application to Singular Extremals. SIAM. J. Control. Optim. 1977, 15 (2), 256-293.
2534
Ind. Eng. Chem. Res., Vol. 46, No. 8, 2007
(7) Aiba, S.; Shoda, M; Nagatani, M. Kinetics of Product Inhibition in Alcohol Fermentation. Biotechnol. Bioeng. 1968, 10, 845-864. (8) BajpaI R. K.; Reub, M. Evaluation of Feeding Strategies in Carbonregulated Secondary Metabolite Production Through Mathematical Modeling. Biotechnol. Bioeng. 1981, 23, 717-738. (9) Ensari, S.; Kim, J. H.; Lim, H. C. Unstructured Model for L-lysine Fermentation Under Controlled Dissolved Oxygen. Biotechnol. Prog. 2003, 19, 1387-1390. (10) Jacobson, D. H. A New Necessary Condition of Optimality for Singular Control Problems. Siam J. Control 1969, 7, 578-595. (11) Hong, J. Optimal Substrate Feeding Policy for a Fed Batch Fermentation with Substrate and Product Inhibition Kinetics. Biotechnol. Bioeng. 1986, 28, 1421-1431. (12) Shin, H. S.; Lim, H. C. Degrees of Freedom Approach to
Admissible Feed-Rate Profiles for Optimal Fed-Batch Fermentation. Ind. Eng. Chem. Res. 2006, 45 (18), 6227-6235. (13) Flores-Cerrillo, J.; MacGregor, J. F. Control of Batch Product Quality by Trajectory Manipulation using Latent Variable Models. J. Process Control 2004, 14, 539-553. (14) Zhang, H.; Lennox, B. Integrated condition monitoring and control of fed-batch fermentation processes. J. Process Control 2004, 14, 41-50.
ReceiVed for reView July 31, 2006 ReVised manuscript receiVed December 24, 2006 Accepted February 12, 2007 IE061001O
