Ind. Eng. Chem. Res. 2006, 45, 4851-4854
4851
Comment on “Optimization of Biochemical Fed-Batch ReactorsTransition from a Nonsingular to a Singular Problem” Hwa Sung Shin* and Henry C. Lim Department of Chemical Engineering and Materials Science, UniVersity of California, IrVine, California 92697 Sir: Jayant and Pushpavanam1 treated singular fed-batch operations by embedding them into a sequence of nonsingular problems that are forced to the singular solution in the limit as F f 0, which is a scheme that was originally proposed by Jacobson et al.2 They were not able to obtain singular solutions, because of numerical difficulty of the proposed approach, and, therefore, they advocated the use of nonsingular solutions to infer the structure for singular solutions. The purpose of this correspondence is to provide complete solutions to the problem using Pontryagin’s Minimum Principle and singular control and analysis of the Hamiltonian function, as well as to highlight several concerns, including the following: (i) the numerical technique used by Jayant and Pushpavanam1 had been used earlier by Menawat et al.,3 who reported numerical instability of such an approach; (ii) among the four cases considered, three cases (cases II through IV) are not properly defined, leading to unbound final times and/or volumes; (iii) they have chosen extremely small final times in their numerical studies so that the cell yields they reported amount to only 0.1%, versus the potential yield of 80%, and, thus, the results are physically impractical; and (iv) their proposed inference to a singular solution from their numerical solutions to the nonsingular problem, the case referenced as Figure 3a by Jayant and Pushpavanam,1 is incorrect.
However, Jayant and Pushpavanam1 chose to have free and fixed final volumes, and, therefore, we proceed without the volume constraint. This problem was imbedded into a nonsingular problem by modifying the objective function and allowing F f 0, so that, in the limit, one can obtain the solution to the original singular problem.
∫0t (F - Fh )2 dt]
MinJ ) lim[J′ ) -XfVf + F F*(t)
Ff0
f
(3)
Because the yield is constant, a stoichiometric relation can be obtained from eq 1:
S ) SF +
(
)
V0 X0 X + S0 - SF - ) g(X,V) V Y Y
(4)
Thus, we can work with two mass balance equations instead of the three:
X˙ ) µ(g)X -
FX V
V˙ ) F
(for X(0) ) X0)
(for V(0) ) V0)
(5a) (5b)
Solution via Pontryagin’s Minimum Principle Mathematical Model The problems treated are maximization cell mass produced at free and fixed final times for the process described by the following mass balances:
d(XV) 5.2S XV ) µXV ) 2 dt S + 25S + 62.9 d(SV) µXV ) FSF ) dt Y 5.2S XV FSF 2 0.8(S + 25S + 62.9) dV )F dt
(
MinH ) λX µ[X,V] F*(t)
(for X(0) ) X0) (1a)
(for S(0) ) S0) (1b)
(for V(0) ) V0)
MaxJ ) XfVf f Min(- J ) -XfVf)
((
FX + λVF V
)
(
( )
[ ] [
)
)
It is apparent that the final volume must be fixed. Otherwise, an infinite final volume would minimize the performance index. * To whom all correspondence should be addressed. E-mail address:
[email protected].
]
(8) The singular feed rate is determined by recognizing that the switching function is identically zero over a finite time interval, φ ) λV ) 0, and, therefore, its derivatives of all orders must also vanish:
φ ) λV ) 0
(2)
F*(t)
)
(6)
∂µ F λX µ + X ∂X V dλ d λX ) )(7) ∂µ FX dt dt λV λX X + 2 ∂x3 V -Vf λ (t ) λ(tf) ) X f ) (free if Vf is fixed, -Xf if Vf is free) λV(tf)
(1c)
The objective is to maximize the total amount of cell mass produced at the final time by applying the optimal feed rate profile, F*(t)
For the previously stated problem (eqs 1 and 2), Pontryagin’s minimum principle (PMP) provides the Hamiltonian to be minimized and the adjoint variables,
(9)
dφ dλV dµ dS ) ) -λXX dt dt dS dV V0X0 dµ 1 ) λXX 2 SF - S0 ) 0 (10) Yx/s dS V
( )(
From eq 10,
10.1021/ie060127i CCC: $33.50 © 2006 American Chemical Society Published on Web 05/24/2006
)
4852
Ind. Eng. Chem. Res., Vol. 45, No. 13, 2006
(
)
dµ 5.2S d ) )0 dS dS S2 + 25S + 62.9
(11a)
S ) Sm ) 7.931
(11b)
or
Equation 11 implies that the substrate concentration must be chosen properly to maximize the specific growth rate. The second time derivative is obtained by differentiating eq 10:
(
)
X0V0 d µ dS dφ ) SF - S0 )0 2 Y dS2 dt dt 2
2
Figure 1. Feed rate profiles for low, appropriate, and high initial substrate concentration (S(0)), small tf, and small Fmax, where Ss is the appropriate initial substrate concentration; the abscissa is time (t) and the ordinate is feed rate (F).
The nonsingular formulation leads to eq 7 from the work of Jayant and Pushpavanam:1
(12) F)F h+
Equation 12 implies that
dS )0 dt
(13)
According to eq 13, the substrate concentration must be constant during the entire singular period. From eqs 11 and 13, we conclude that the substrate concentration is maintained constant at the singular arc Sm (Sm ) 7.931) corresponding to the maximum specific growth rate, µ(Sm) ) µmax. Therefore, the singular feed rate is obtained from eq 2:
σmXV Fsin ) SF - Sm
(
1 λXX - λV 2F V
)
The proposed solution for the singular case is obtained by solving the nonsingular problem repeatedly with decreasing values of the weighting factor F until F f 0:
lim F ) F h+ Ff0
(
1 λXX - λV 2F V
)
(16)
We note in eq 16 that, to obtain a bounded solution for singular case, the term in the parenthesis (λXX/V - λV) must also vanish as F f 0:
λ XX - λV f 0 V
(14)
(17)
Therefore, eq 16 becomes Therefore, the optimal feed rate profile should maximize the specific growth rate to produce the maximum amount of cells. Thus, the optimal feed rate profile should first force the substrate concentration to Sm as fast as possible and then remain at that point until the reactor is full. Therefore, the first portion of the optimal feed rate should be Fmax if S(0) < Sm and Fmin ) 0 if S(0) > Sm, followed by the singular feed rate Fsin until the reactor is full (the time that corresponds to this point is denoted as tfull) and, finally, a batch period if the given final time is greater than the time need to fill the reactor (tf > tfull). The optimal sequence is
Fmax f Fsin f Fmin ) 0
(if S(0) < Sm)
Fmin ) 0 f Fsin f Fmin ) 0 Fsin f Fmin ) 0
(15a)
(if S(0) > Sm) (15b)
(if S(0) ) Sm)
(15c)
If the final time tf and the maximum feed rate Fmax are not sufficiently large enough, the aforementioned profiles face constraints and can lead to various altered profile forms, as shown in Figure 1. With the above treatment, we proceed to consider the work of Jayant and Pushpavanam.1 Numerical Technique The method of embedding the singular problem in a sequence of nonsingular problems that converge to the singular solution in the limit as F f 0 (actually, 1/ f 0) was originally proposed by Jacobson et al.2 and was used by Menawat et al.3 for baker’s yeast and obtained the singular feed rate portion only, rather than the entire feed rate profile. They reported that numerical instabilities are likely to occur as F becomes small and indicated the difficulty of such a scheme.
h+ lim F ) F Ff0
)F h+
(
1 λXX - λV 2F V
)
(λXX/V - λV) f 0 2F f 0
(18)
Because of the sensitivity of the last term, numerical conversion may be difficult to realize. Equation 17 is, precisely, the socalled switching function φ(t) ) 0 over a finite time interval. It is not certain how one would force numerically the switching function to be zero over a finite time interval as F f 0. We know from the singular solution that the switching function can be positive, negative, or zero over finite time intervals. It is unclear how one would overcome the numerically sensitive limiting process on a nonsingular problem to obtain the solution to a singular problem. Cases I, II, III, and IV Jayant and Pushpavanam1 considered four cases: case I (fixed tf and fixed Vf), case II (fixed tf and free Vf), case III (free tf and fixed Vf), and case IV (free tf and free Vf). Let us look at these cases individually. The intent is to obtain singular solutions from their nonsingular solutions. Case I: Fixed Final Time tf and Fixed Final Volume Vf. For this case, there are four boundary conditions: three on the state variables and one on the adjoint variable.
X(0) ) X0 V(0) ) V0 V(tf) ) Vf and
Ind. Eng. Chem. Res., Vol. 45, No. 13, 2006 4853
λX(tf) ) -Vf
which reduces to
Because we have four differential equationsstwo each for the state and adjoint variablesswe need four boundary conditions, and, therefore, this case is properly defined. Case II: Fixed Final Time tf and Free Final Volume Vf. For this case, there are only three boundary conditions: two on the state variables and one on the adjoint variable.
H(tf) ) λX(tf)X(tf)µ(tf) + FF h2 ) 0 which, under the limiting process, as F f 0, yields
H(tf) ) λX(tf)X(tf)µ(tf) ) 0
µ(tf) ) 0
V(0) ) V0
λV(tf) ) -Xf Because we have one less boundary condition than necessary to solve four differential equations, this case is ill-posed. The performance index is linear in Vf; therefore, it is obvious that one should pick the largest final volume possible to minimize the objective function, Vf ) ∞. Another way to look at this case is through the Hamiltonian. The Hamiltonian is given as
[
H ) λX Xµ -
FX + λVF + F(F - F h )2 ) constant (19) V
]
Assume that there is a short batch period after the bioreactor volume is full. The feed rate then must be zero at the final time. Thus, at the final time, and as F f 0, the Hamiltonian reduces to
H(tf) ) λX(tf)X(tf)µ(tf) ) -V(tf)X(tf)µ(tf) < 0, constant (20) To minimize the Hamiltonian, one would pick Vf to be as large as possible, Vf ) ∞. Case III: Free Final Time tf and Fixed Final Volume Vf. For this case, there are four boundary conditions: three on the state variables and one on the adjoint variable.
Thus, the specific growth rate must vanish at the final time tf. This implies that the substrate concentration must be either zero or very large (inhibition), so that the specific growth rate vanishes. The substrate concentration cannot exceed the feed concentration SF. Therefore, the substrate concentration must be zero. Theoretically, the substrate concentration would require an infinite time to reach the zero value, and, therefore, this implies that the final time is also infinity: tf f ∞. Essentially, the minimum principle yields an infinite final time if the final time does not appear explicitly in the objective function and is free and the cells do not decay. Because the final time is infinity, all of the substrate supplied to the reactor would be converted to cell mass at the specified 80% yield: 1500 g × 0.8 ) 1200 g. Thus, any mode of feeding would lead to an 80% conversion. The solution is nonunique and there is no unique optimal feed rate profile. The problem is thus ill-posed. Case IV: Free Final Time tf and Free Final Volume Vf. This case is a combination of cases II and III. Basically, one must iterate on both Vf and tf. As we have previously observed, the optimal solution demands both an infinite reactor volume (Vf f ∞) and infinite final time (tf f ∞). This also is an illposed problem. From the aforementioned argument, we conclude that the only valid case is case I (fixed final time and fixed final volume) and cases II, III, and IV are ill-posed.
The kinetic parameters that Jayant and Pushpavanam1 used correspond to the maximum specific growth rate of 0.127 that was occurring at the substrate concentration of 7.93: µ(S ) 7.93) ) 0.127. Thus, the minimum doubling time is given as
V(0) ) V0 V(tf) ) Vf and
td,min ) λX(tf) ) -Vf
Because we have four differential equations, two each for the state and the adjoint variables, we need four boundary conditions, and, therefore, this case seems to be properly defined. However, the free final time causes a problem. It is apparent that there is no penalty for a large final time (tf f ∞), and, therefore, it is possible to allow the cell growth to proceed over an infinite time period, yielding the theoretical yield of 80%, regardless of any feed rate profile used, including a batch. We can also show this more rigorously by analyzing the Hamiltonian. Assume that there is a short batch period after the reactor volume is full. The feed rate then must be zero at the final time. The Hamiltonian is identically zero, because tf is free:
[
FX + λVF + F(F - F h )2 ) 0 V
]
(24)
Numerical Examples
X(0) ) X0
H ) λX Xµ -
(23)
Because λX(tf)X(tf) * 0, eq 10 implies that
X(0) ) X0
and
(22)
(21)
ln 2 ) 5.444 0.1273
Jayant and Pushpavanam1 used final times of 10, 11.34, and 32.2 h for cases II, III, and IV, respectively. In fact, the actual average doubling time (due to varying substrate concentration) is >5.444. Thus, the final times chosen represent, at most, 1.8, 2.08, and 5.91 times the minimum doubling time. These final times are far too short. In fact, the cell mass obtained by fedbatch operation is insignificant, in comparison to the potential cell mass that could have been obtained. Table 1 shows the amount of cells obtained and the maximum possible cell mass. Because the feed substrate concentration is 500 g/L and the available reactor volume is (Vf - V0 ) Vf - 7), the maximum cell mass that can be obtained is [500 × (Vf - 7) × 0.8] g. Table 1 clearly indicates that the conditions chosensin particular, the final timeslead to negligible conversions for all cases (∼0.1%). The final time must be chosen carefully by considering the value of the specific growth rate and the amount of substrate to be used. In these examples, the final time is at
4854
Ind. Eng. Chem. Res., Vol. 45, No. 13, 2006
Table 1. Numerical Results Reported by Jayant and Pushpavanam1
case I II III IV
I
II
III
IV
V
VI
tf (h)
Vf (L)
X0V0 (g)
XfVf, obtained (g)
potentially obtainable maximum cell mass, SF(Vf - V0)YX/S ) 400(Vf - 7) (g)
actual cell mass obtained, as a percentage of the potential maximum, IV/V × 100
10.0 10.0 11.34 32.20
10.0 9.84 10.0 15.08
0.001 0.001 0.001 0.001
1.243 1.264 1.352 2.517
1200 1136 1200 3232
0.104 0.111 0.113 0.078
least 2 orders of magnitude smaller than practical. Thus, the examples treated have very little resemblance to practical situations. Inference from a Nonsingular Solution to a Singular Solution Jayant and Pushpavanam1 claimed that the solution to a nonsingular problem may be extrapolated to infer the solution to the singular problem. Although this is a feasible approach, one should be cautious with this approach, because it is difficult to guess the sequence of feed rates for the singular problem from the numerical results of the nonsingular problem. The proposed sequence proposed by those authors1 for the singular problem from their numerical results (Figure 3a in their paper1), Fmax f Fmin f Fsin, is incorrect. There is no period of minimum feed rate that follows the period of maximum feed rate. The optimal sequence for the singular case is Fmax f Fsin f Fmin ) 0. Numerical results for the nonsingular case seem to average out the feed sequence of the singular case, which is very similar to a Fourier series representation of a sequence of piecewise continuous functions. However, it is difficult to assess exactly whether the averaging is done over a sequence of Fmax f Fmin f Fsin or Fmax f Fsin. Even experts such as Jayant and
Pushpavanam1 missed this point. If the objective is to gain qualitative insight to the optimal feeding policy, a straightforward application to the singular problem may yield much more insight than the numerical solutions of the nonsingular problem. Some theoretical analyses prior to a purely numerical approach is both instructive and beneficial in obtaining solutions to singular feed rate problems in fed-batch cultures. Acknowledgment One of the authors (H.S.S.) was supported in part by Korea Science and Engineering Foundation (KOSEF) and University of California, Irvine. Literature Cited (1) Jayant, A.; Pushpavanam, S. Optimization of a Biochemical FedBatch ReactorsTransition from a Nonsingular to a Singular Problem. Ind. Eng. Chem. Res. 1998, 37, 4314-4321. (2) Jacobson, D. H.; Gershwin, S. B.; Lele, M. M. Computation of Optimal Singular Controls. IEEE Trans. Autom. Control 1970, AC-15 (1), 67-73. (3) Menawat, A.; Mutharasan, R.; Coughanowr, D. R. Singular optimal control strategy for a fed-batch bioreactor: Numerical approach. AIChE J. 1987, 33, 776-783.
IE060127I
