1998
Ind. Eng. Chem. Res. 1999, 38, 1998-2004
Optimization of a Biochemical Fed-Batch Reactor Using Sequential Quadratic Programming S. Pushpavanam,* Sanyasi Rao, and Imran Khan Department of Chemical Engineering, Indian Institute of Technology, Kanpur 208016, India
In this paper, we discuss the use of sequential quadratic programming to obtain the optimum feed rate strategy for the fed-batch operation of two biochemical systems: (I) lysine and (II) alcohol fermentation. We consider the constrained optimization problem where the reactor operation is such that the final volume is fixed. Both fixed and free time of operations are considered. In this approach, the entire duration of the reactor operation is divided into equal subintervals. The feed rate is approximated as discretized charges which are added to the reactor at the beginning of each subinterval. The effect of the number of stages on the optimum performance for both systems is investigated. The results are compared with those obtained using Pontryagin’s maximum principle in the literature. The advantages and the applications of using this approach are discussed. Introduction The rates of many biochemical reactions are characterized by substrate inhibition and or product inhibition. Such systems are known to exhibit an optimum performance when the reactor is operated in the fed-batch mode. Here the substrate is fed to the reactor containing an initial amount of biomass over a period of time. Pontryagin’s maximum principle has been used to solve such optimal control problems in the past.1 Because the control variable, i.e., the substrate feed rate, occurs linearly in the system, the optimal policy is bang-bang or follows a singular trajectory or both. The control variable is maintained at its upper or lower bound, in bang-bang control. Along the singular arc, the control variable varies continuously between the lower and the upper bound. The optimal behavior of the reactor in this approach has to be determined using different systems of equations in the different regions, depending on certain conditions being satisfied. The optimal profile for such singular problems can be discontinuous. The numerical computation of such profiles usually requires a priori information on the sequence of control actions which render the profile optimal. This is usually based on an understanding of the physics of the interactions between the different variables. Alternatively, it may be possible to determine the sequence of control actions on the basis of some mathematical analysis, as carried out by Modak and Lim.1 However, these methods are restrictive and may be applicable only to a limited set of problems. Moreover, the experimental implementation of a substrate feed rate which varies continuously with time will require very intricate control instrumentation and controllers. The optimal performance of fermentation reactors in the fed-batch mode has been investigated with the objective of maximizing either (i) the biomass production or (ii) the product formed.2,3 Two classical methods that have been used to study the problem are (i) Pontryagin’s * Current address: Department of Chemical Engineering, I. I. T Madras 600036, India. E-mail:
[email protected]. Fax: (91) (044)-2350509. Telephone: (91) (044)-235-1365 extn. 3934.
continuous maximum principle and (ii) nonlinear optimization techniques such as sequential quadratic programming (SQP).4,5 The model system most extensively investigated theoretically is characterized by a single reaction. Here, reaction rate varies linearly with biomass and is assumed to have a Haldane dependency on the substrate (the system thus exhibits substrate inhibition). The manipulated variable, i.e., the substrate feed rate (used to determine the optimal conditions), occurs linearly in the governing ordinary differential equations. This renders the problem singular. The optimal feed policy as determined using Pontryagin’s maximum principle is hence either bang-bang or follows a singular arc. The optimal profile of the manipulated variable hence can be discontinuous. The determination of the optimal policy requires solving different sets of equations along different regions as determined by the necessary and sufficient conditions of the maximum principle. In addition, the numerical scheme is very sensitive to the initial guesses. In most problems, one almost needs to know the values of the variables being iterated upon (i.e., the solution) before one can find them while using this procedure. Nonlinear programming techniques such as SQP can be used to determine a discrete set of variables. Hence, this method is based on discrete approximations of a continuous variable. Physical intuition is used to formulate the problem and to render it amenable for solution using this technique. Here, the numerical method is computationally intensive, and the problem with sensitivity to the initial guesses is eliminated. This technique is hence very robust. Weigand6 considered the problem of maximizing cell production in a repeated fed-batch mode of operation. He obtained an analytical expression relating the total reaction time with biomass concentration under optimal conditions. In his analysis, however, he did not assume an upper bound on the permissible flow rate. This has led to his results being of only theoretical significance. Bonte et al.7 and Modak et al.8 determined the optimal feeding profiles for different kinetic models. They developed a computational scheme to obtain the feeding policy, which maximized (a) penicillin produc-
10.1021/ie9805483 CCC: $18.00 © 1999 American Chemical Society Published on Web 03/26/1999
Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999 1999
tion and (b) cell mass production. The optimal control problem considered in these works is singular in nature. The optimal policy consists of different segments: (i) when the feed rate is held at its minimum, (ii) when the feed rate is held at its maximum, and (iii) when the feed rate follows a singular arc. Cazzador2 qualitatively analyzed the effect of initial conditions and substrate feed concentrations on the different profiles. His method of analysis can be applied only for one-dimensional problems where the system is described by one state variable. Here, the governing equations describe the evolution of a single scalar dependent variable. The numerical computation of the optimal profile for most real problems, however, requires physical insight into the problem. This is necessary to determine the exact sequence of the different portions which constitute the optimal profile.8 Such an approach may not be applicable to multivariable systems where the interactions are very complex or are properly understood. In a nonbiochemical engineering context, optimal feed rate policies were obtained using pulse feeds at discrete intervals of time by Levin,5 Biegler,9 Morrison and Sargent,10 and Vasiliadis et al.11 These researchers have used sequential quadratic programming to determine the optimal profile of the manipulated variable. Shukla and Pushpavanam12 used a sequential quadratic programming technique to obtain the optimum feed rate profile for a repeated fed-batch mode of operation. They divided the entire interval of operation into different segments. The feed rate in each interval was assumed to be constant. The optimal feed rate policy thus determined does not use any information from the physics and is based on mathematically approximating a continuous function as a sequence of constants over different subintervals. This approach provides a general mathematical basis for obtaining the optimal solution. It is useful when the interactions between the different variables is complex and when it is not possible to infer the sequence of control policies which will give rise to optimal behavior. The disadvantage of this technique lies in the fact that it is computationally intensive. In this paper, we discuss how sequential quadratic programming can be used to obtain a practical feeding policy which optimizes the performance of a multivariable biochemical reactor. In this approach, we approximate the feed rate over the reactor using instantaneous pulse charges. Here, the entire cycle of reactor operation is divided into different subintervals. The substrate is assumed to be added instantaneously at discrete points of time. These time instants are at the beginning of each subinterval. Consequently, the various concentrations exhibit discontinuities here. In each subinterval, the reactor is operated in a batch mode. The pulses which are added at the beginning of a subinterval instantaneously change the reactor state. This new state determines the initial conditions for that subinterval. The optimization problem is now reduced to determining the different pulses which maximizes the total amount of product formed at the end of the reaction time. The basis of the proposed approach consists of mathematically approximating the continuous profile using discrete pulses. For the sake of computational convenience, the subintervals are each taken to be of equal duration. The pulse determined at the beginning of a particular subinterval approximately represents the
amount of substrate which would be added while using the continuous optimal profile over that subinterval. Model Formulation We now consider the operation of a biochemical reactor which contains biomass at some initial concentration X0. The substrate S is added continuously to the reactor during the period of operation. This feed contains neither the biomass nor the product. The evolution equations which described the reactor when the feed rate F is a continuous function of time are
X˙ ) µX -
FX V
F S˙ ) -σX + (SF - S) V P˙ ) πX -
FP V
V˙ ) F
(1a) (1b) (1c) (1d)
Here, X, S, and P represent the concentration of cell mass, substrate, and product, respectively, in the reactor, and V is the volume of the reactor. The conservation of volume eq 1d is more appropriate in biochemical reactions as it allows us to neglect the volume of the gases realeased on account of reaction i.e., CO2. Although the volume of CO2 in the gaseous phase may be very significant, the equivalent volume in the condensed phase is assumed to be negligible. µ, σ, π represent the rates at which biomass is produced, substrate is consumed, and the product is formed, respectively. The evolution of the states of the reactor are subject to the initial conditions (X0, S0, P0, and V0). The objective here is to determine the feed rate profile F(t) which maximizes the product concentration at the final time tf, i.e.,
Maximize P(tf) F(t) This problem was solved by Modak and Lim1 using Pontryagin’s maximum principle. The convergence of the numerical scheme while using this approach is very sensitive to the initial guess, as discussed earlier. This has been observed in the analysis of many nonsingular problems. In a singular problem, this is compounded because we have to integrate different systems of equations in different regions and the optimal profile can be discontinuous. With a view to overcoming these limitations, we use sequential quadratic programming to solve this problem in this paper. SQP can be used to determine a discrete set of variables. The continuous optimal policy is hence approximated mathematically as a discrete set of pulses, and the time of operation is divided into equal intervals. The choice of equal subintervals is only to make the problem algebraically tractable. For the case of equal subintervals and a given duration of operation, we have to find only the N pulses. This will determine the optimal policy for a given N. For the case of unequal subintervals, we have to find the N pulses as well as the duration of the N subintervals. For large values of N, both results appoximate the continuous feed rate policy and hence yield qualitatively identical results. In this work, we approximate the feed rate policy by discrete charges as explained earlier. The entire interval of operation is divided into N equal subintervals. The
2000 Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999
Figure 1. Schematic picture depicting the physical significance of the different variables as proposed in our approach.
optimization problem now reduces to determining the best possible approximation of the optimal profile using instantaneous pulses for a given number of pulses or subintervals. The reactor operation in each subinterval is in the batch mode. The evolution in the ith subinterval (where ti-1 < t < ti) is governed by the equations of the batch reactor
X˙ ) µX
(2a)
S˙ ) -σX
(2b)
P˙ ) πX
(2c)
These equations represent the reactor state during the batch period before the addition of the pulse at ti and after the addition of the pulse at ti-1. Figure 1 depicts the notation and the physical significance of the different variables which are being used and will be referred to later. Thus Vi represents the volume of the reactor in the interval (ti-1, ti) and ∆Vi is the pulse being added at ti-1. The evolution of the system is subject to initial conditions. These conditions specify the state at the beginning, i.e., t0- or before the addition of the pulse at t0.
X(to ) ) Xo
(3a)
S(t0 ) ) S0
(3b)
P(t0 ) ) P0
(3c)
V(to)
(3d)
) Vo
From the mass balance equations at the time instants ti-1 when the pulses ∆Vi are added, we have,
Vi-1 X(ti-1) Vi
(4a)
Vi-1 Vi - Vi-1 S(ti-1) + SF Vi Vi
(4b)
Vi-1 P(ti-1) Vi
(4c)
+ )) X(ti-1
+ S(ti-1 ))
+ P(ti-1 ))
Here the subscripts F and i stand for feed conditions + and for the ith time instant, respectively. Here ti-1 , ti-1 represent the instants immediately after and just before the addition of the pulse at ti-1. The pulse is assumed to contain only substrate and no biomass or product. Our objective is to seek the pulses ∆Vi’s (∆Vi ) Vi+1 Vi) to be added at ti’s for i ) 0, 1, ..., N - 1 which will maximize the product concentration P(tf). We assume throughout in this paper that the final volume Vf is fixed.
are of equal duration. The substrate feed addition, if any, is at the beginning of each subinterval. The optimization problem now is transformed to determining the optimal feed pulses to maximize final product concentration P(tf). Here the duration of each subinterval is known for a given N when the total cycle time is fixed. When solving the free time problem, we iterate on the time of reactor operation also maintaining the duration of the N subintervals equal. The differential equation 2a-c governs the system evolution. Constraints are imposed on the final reactor volume, i.e., at the final time (tf).
V(tf) ) Vf
(5)
In this approach, the different charges (∆Vi’s) to be found for the optimal performance have to satisfy the constraint in eq 5. The numerical method used to solve this optimal control problem is sequential quadratic programming. The convergence of this technique requires the accurate determination of derivatives of the objective function with respect to the manipulated variables in this case pulses (∆Vi’s). We have used the state trajectory sensitivity method proposed by Rosen and Luus,13 to calculate these derivatives which ensure convergence. We describe the application of this method to our system now. Calculation of Sensitivity with Respect to ∆Vi’s, tf’s. Differentiating eq 2a-c with respect to ∆Vj’s, we have for j ) 1, ..., N
∂ ∂ ∂ ∂X ∂S ∂P ∂X˙ ) (µX) + (µX) + (µX) (6a) ∂∆Vj ∂X ∂∆Vj ∂S ∂∆Vj ∂P ∂∆Vj ∂S˙ ∂X ∂S ∂P ∂ ∂ ∂ ) (-σX) + (-σX) + (-σX) ∂∆Vj ∂X ∂∆Vj ∂S ∂∆Vj ∂P ∂∆Vj (6b) ∂P˙ ∂ ∂ ∂ ∂X ∂S ∂P ) (πX) + (πX) + (πX) (6c) ∂∆Vj ∂X ∂∆Vj ∂S ∂∆Vj ∂P ∂∆Vj Equations 2a-c and 6a-c are to be solved simulta+ neously from tj-1 to tj where j ) 1, ..., N. This follows because the addition of the jth pulse ∆Vj influences only the reactor operation from tj-1 to tf. In particular, it cannot affect the reactor performance before tj-1. The initial conditions for the sensitivity with respect to the jth pulse ∆Vj at tj-1 is given by differentiating the updating conditions (in eq 4a-c)
Vj-1 ∂X + (t ) ) -X(tj-1 ) 2 ∂∆Vj j-1 V
(7a)
j
Vj-1 ∂S + (t ) ) (SF - S(tj-1 )) 2 ∂∆Vj j-1 V
(7b)
Vj-1 ∂P + (t ) ) -P(tj-1 ) 2 ∂∆Vj j-1 V
(7c)
j
j
Method of Solution We now discuss how the method can be applied to determine the optimal solution when the N subintervals
The updating conditions for the sensitivity at the subsequent junction points ti-1’s for i ) j + 1, ..., N are
Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999 2001
given by differentiating eq 4a-c and remembering that subsequent volumes are a function of ∆Vj’s
Table 1. Lysine Fermentation Case Study
Vi - Vi-1 Vi-1 ∂X ∂X + (ti-1) ) X(ti-1 ) + (t ) ∂∆Vj Vi ∂∆Vj i-1 V2
charges
(8a)
i
N
1 15
3
5
7
10
0 6.7378 8.2622
0 0.9613 4.4469 5.6417 3.9501
0 0 1.8789 3.3297 3.9884 4.6193 1.1837
0 0 0 1.4927 2.2145 2.5581 2.8924 3.2212 2.6211 0 34.654
[S(ti-1 ) - SF](Vi - Vi-1) Vi-1 ∂S ∂S + (ti-1) ) + (t ) ∂∆Vj Vi ∂∆Vj i-1 V2i (8b)
Vi - Vi-1 Vi-1 ∂P ∂P + (ti-1) ) P(ti-1 ) + (t ) ∂∆Vj Vi ∂∆Vj i-1 V2
(8c)
i
Integrating eqs 2a-c and 6a-c with these updating conditions from t+ j to tN yields the sensitivity ∂X(tf)/ ∂∆Vj. When solving the optimization problem as a free time problem, the sensitivity with respect to final time is also required. This is obtained directly from eq 2c when the right-hand side is evaluated at tf. The function “CONSTR” available in MATLAB is used to solve this problem using sequential quadratic programming. Results and Discussion Fixed Final Time and Volume Operation. The two examples chosen for the study in this work are lysine and alcohol fermentation. These problems have been studied in detail previously.1 For lysine fermentation, the rate expressions as given in Modak and Lim1 are µ ) 0.125S, π ) -384µ2 + 134µ, σ ) µ/0.135. We emphasize that the units of the substrate concentration while using the kinetic expressions for lysine fermentation should be in weight percent. For alcohol fermentaion, we again follow Modak and Lim1 and use the expressions given below
µ)
0.408S -0.028P e (0.22 + S)
π)
S e-0.015P (0.44 + S) σ ) µ/0.1
The substrate and product concentration in these expressions are given in g/L. We first consider the operation for which the final reactor volume and the duration of the entire cycle are both externally fixed at the values determined by the optimal profiles of Modak and Lim.1 In this method, all subintervals are considered to be of equal durations. The duration of each subinterval is hence determined once we fix N, the number of stages. The solution to this optimization problem hence consists of determining the N pulses (∆Vi’s) which maximize P(tf). The conditions of the operation were chosen to be the same as that by Modak and Lim.1 This facilitates a comparison of the numerical resuslts we obtain with those obtained by them. For the lysine fermentation problem, the initial volume of the reactor V0 was chosen as 5 L, the initial biomass concentration was 0.02 g/L, the substrate feed concentration SF and the initial substrate concentration S0 were 2.8 wt %, and the final volume was 20 L. The time of reactor operation was taken as 33.62 h. The optimal profile of the control variable depends on the initial concentration of the reactor. For the choice above,
obj fn P(tf) gm/1 25.2249 32.350
33.6959 34.216
Modak and Lim1 had observed that the optimal profile was characterized by a batch period followed by a singular arc trajectory and again followed by a batch period. The results of the optimization for the lysine fermentation problem using our scheme is depicted in Table 1. It can be clearly seen that as we increase the number of pulses, N the objective function (the product formed) increases. This is to be expected because increasing N results in a better approximation of the true optimal policy. The optimum policy for the chosen conditions consists of three segments and an initial batch period (because the charges in the first few intervals are 0) followed by a period during which the flow rate is nonzero (corresponding to nonzero charges in the middle subintervals) and then by another batch period.1 The discretized approximation that we have obtained clearly emulates the behavior obtained by Modak and Lim.1 In addition, the total product obtained using the optimal profile from our strategy for N ) 10 is 693.084 g. This compares favorably with the value of 687.2 g previously reported.1 Figure 2a,b shows the variations of the biomass and product concentrations as the reaction progresses in the reactor for the case of N (the number of pulses) being 7. It is clear that the introduction of the pulses at the beginning of each subinterval results in the observed discontinuities of the profiles. The problem formulated by Modak and Lim1 uses the total product formed with a fixed final volume minus a cost function as the objective function. The cost function was taken to be proportional to time of operation tf. We have maximized the product concentration in our approach using the final time obtained by Modak and Lim.1 These two are equivalent because in our operation the final volume of the reactor and the final time of operation are both fixed. For alcohol fermentation, the optimal policy as determined using the maximum principle1 also follows the same sequence, i.e., batch-singular arc-batch for some initial conditions. For this system, the initial biomass concentration was chosen as 0.2 g/L. The initial and the feed substrate concentrations were fixed at 100 g/L. The initial and final volumes of the reactor are again 5 and 20 L, respectively. The time of operation of the reactor was fixed at 20 h as read from the graph of Modak and Lim.1 The optimal pulse distribution for N ) 10 is able to capture the trend of the optimal profile, i.e., batchsingular arc-batch as can be seen from Table 2. Figure 3a,b depicts the variations of the biomass and product concentrations as the reaction progresses for N ) 7. The discontinuities in the profile are due to the addition of the pulse charges at the beginning of the subintervals. It is also clear from these results that the duration of the initial and final batch periods as found previously1 agree well with those found in this work.
2002 Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999
Figure 2. (a) Evolution of the biomass concentration with time for lysine fermentation problem (N ) 7) when optimal profile is maintained, (fixed final volume and time). (b) Evolution of the product concentration with time for lysine fermentation problem (N ) 7) when optimal profile is maintained (fixed final volume and time). Table 2. Alcohol Fermentation Case Study N charges
1 15
3
5
7
10
0 3.6821 11.317
0 0 1.3721 5.0659 8.5620
0 0 0 0 3.3691 5.0095 6.6214
0 0 0 0 0 0.2792 2.6670 3.8467 5.5300 2.6771 32.537
obj fn P(tf) gm/1 29.082 30.666
31.616
32.134
Effect of Initial Conditions. The reactor operation discussed in this paper is the nonrepeated fed-batch operation. The performance of this reactor depends on the choice of the initial conditions. Modak and Lim1 projected the singular arc trajectory in the twodimensional plane of alcohol and substrate concentration. They established that for a choice of initial conditions to the right of this curve the optimal profile consists of three segments, batch-singular arc-batch, and initial conditions to the left would not follow a singular arc. For the alcohol fermentation problem, we studied the effect of the initial conditions on the reactor and determined its effect on the optimal profile. The final
Figure 3. (a) Evolution of the biomass concentration with time for alcohol fermentation problem (N ) 7) when optimal profile is maintained (fixed final volume and time). (b) Evolution of the product concentration with time for alcohol fermentation problem (N ) 7) when optimal profile is maintained (fixed final volume and time).
volume was constrained to 20 L, and the objective function was taken to be P(tf). For S0 ) 30 g/L and SF ) 100 g/L, the optimal sequence of pulses is (0.465, 0, 0.858, 2.43, 3.15, 3.7, 4.4). This choice of initial conditions correspond to points to the left of the projected singular arc of Modak and Lim.1 The addition of the first pulse clearly indicates that the initial portion of the reactor is not a batch operation. This addition can be interpreted as forcing the operating point to the right of the singular arc projected on the two-dimensional plane. The reactor is then operated in a batch mode until the evolution of the trajectory intersects the singular arc. This singular arc is now approximated by the pulses obtained. Thus, it does not necessarily mean that for points corresponding to initial conditions to the left of the singular arc projection found by Modak the optimal trajectory will not contain a singular arc. This is only to be expected because we are studying the evolution of the system when it is subject to the addition of pulses. These change the location of the operating point in the phase plane. A second reason the projection cannot be used is that the biomass concentration also varies along the singular arc and this is not reflected in the projection. Free Volume, Fixed Time Problem. We also carried out the simulation of the optimization problem
Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999 2003
(maximizing product concentration) for the case where the final volume was left free i.e was not set to a predetermined value. For the alcohol fermentation problem with the same set of initial conditions as those used in obtaining Table 2 and for N ) 7, the resultant set of charges was (0, 0, 0, 0, 3.8467, 7.2876, 9.2584). This results in a final volume of 25.3927 and a final product concentration of 32.2395 g/L. The value of this product concentration is not very different from the value obtained when the final volume is fixed to 20 L. However as expected, the objective function is marginally more than that obtained when using a fixed final volume of 20 L. To verify the results of the free volume problem, the constrained optimization program was then simulated for the case when the final volume was constrained at 25.3927 L. This generated the same optimal profile and the same objective function as that obtained using the unconstrained problem, thereby validating the numerical code. When we maximized the total product formed [P(tf) V(tf)] for this case, the algorithm converged on values which resulted in an infinite volume of the reactor. It is hence necessary to constrain the final volume of the reactor when using the product formed as the objective function. This is equivalent to maximizing the product concentration. Fixed Volume, Free Time Problem. We have also obtained the optimal solution to our system for the case when the final volume of the reactor is fixed, but the final time of the reactor operation is not. This choice replicates the work of Modak and Lim.1 To compare our results with theirs, we used the same objective functions as they did. Thus, for the lysine fermentation, they chose the objective function as [P(tf)V(tf) - tf] and for the alcohol fermentation problem as [P(tf)]. The lysine fermentation problem was simulated for the case of free final time and fixed final volume. Two different choices of the objective function were used: (i) the maximum profit function P(t)fV(tf) - tf and (ii) product maximization P(tf). Here represents the cost of operation of the reactor. For the first choice with ) 10, the SQP algorithm converged to a unique solution for a wide range of initial guesses. This was characterized by the sequence of pulses for N ) 7 as (0, 0, 1.859, 3.256, 4.064, 4.769, 1.05) and a time of operation as 31.9341 h. The objective function obtained was 350.588 g and this compares favorably with that obtained by Modak and LIm1 (353.58 g). Clearly, our policy is suboptimal because N is finite. For the second choice of the objective function, i.e., with ) 0, the SQP algorithm converges to a solution of pulses (0, 1.891, 2.1376, 2.3817, 2.6236, 2.8637, 3.1023) and tf ) 407 h. The objective function was found to be P(tf) ) 810 g. The program was found to be insensitive to the value of the final time, beyond the tf determined above. Thus, initial guesses starting with a high value of tf converged to a high value of tf and gave rise to the same objective function. Initial guesses of tf beyond 407 h converged to an optimal solution which gave the same objective function as determined above. The insensitivity of the objective functions to tf can be explained as follows. Because we are constraining the final volume, the total amount of substrate added is fixed. The reaction does not occur after a certain time, because all substrate is consumed beyond this. The amount of product formed is a constant because this is
Table 3. Effect of µm on Optimal Trajectories and Objective Function for Seven Pulses µm 0.3
0.34
0.38
0.4
0.46
0.5
∆V1 15 10.8776 0.7959 0 0 0 ∆V2 0 0 0 0 0 0 ∆V3 0 0 0 0 0 0 ∆V4 0 0 0 0 0.8633 1.6765 ∆V5 0 0 1.667 3.1839 3.4979 3.3268 ∆V6 0 0 5.067 5.0231 4.6211 4.3542 ∆V7 0 4.1224 7.470 6.7929 6.0178 5.6425 obj. fn. 30.6680 37.6300 35.0432 32.9779 27.3912 24.6395
constrained by the amount of substrate added and stoichiometric considerations and not by the time of reaction. So allowing the reaction to proceed for a time greater than a critical value does not change the objective function. Hence, the time of reaction is not uniquely fixed for the latter choice of objective function. When we maximize a profit function with > 0 as in case i, the algorithm tries to minimize the time of operation, to reduce the cost. Hence in this case, the optimal solution is unique. We have also observed a similar trend in the behavior for the alcohol fermentation problem. We have also found the effect of parameters µm and km on optimal trajectories. The constant 0.408 occurring in the expression for µ in the alcohol fermentation problem is denoted as µm. It is varied from 0.3 to 0.5. Thus for the lower value of µm, i.e., 0.3, the optimal pulse sequence is a batch operation, with all of the charges added at the first stage (Table 3). As µm is increased, the optimal sequence changes to one for which the pulses have to be added at the later stages. This change in optimal feeding policy arises because of the sensitive dependence on µm of S*, the substrate concentration for which we have a maxima in the rate expression. As we vary µm, the optimal profile varies in such a way that on the average we maintain S near S*. Calculations have revealed that the optimal trajectory is not very sensitive to km (0.22 in the expression of µ). This occurs because the km value is significantly smaller than the range of values S attains [km is O(10-1) and S is O(102)]. Consequently the dependence of optimal trajectories on km is very weak. The objective function increases and decreases with µm as depicted in Table 3. This arises because the rate varies with time in a different manner as we change µ. This is also caused by the variation in the optimal trajectories and eventually results in a different value of the objective function. Conclusions The optimum feed policy we have determined consisting of discretised charges, discussed here, can be implemented practically in a laboratory scale as well as in the plant scale. It requires only the addition of the substrate at regular intervals of time in some predetermined fashion. This would enable us to dispense with having to use a stepper motor, etc., to vary the flow rate in a predetermined way. This would drastically reduce the capital cost of the reactor operation, thereby increasing the profit considerably. For a finite N, the policy obtained by our method is sub optimal, and it approaches optimal policy as N tends to infinity. The pulse profile obtained also provides insight into the sequence of profiles which have to be maintained
2004 Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999
for optimal performance when the flow rate into the reactor is a continuous variable. This information can now be used to obtain the solution to the optimal control problem by the use of the continuous maximum principle. In other words, sequential quadratic programming can be used to generate good initial guesses which could be used with Pontryagin’s maximum principle. It must be emphasized here that the convergence of SQP occurs to the same solution for a variety of initial guesses. This establishes that the method proposed here is very robust. However, this convergence depends on the accurate determination of the sensitivities of the objective function and the constraints with respect to the variables we are optimizing. We have studied the effect of different choices of objective functions and when constraints are imposed on the final volume and time. The optimal profiles obtained depend on both of these factors. The sequence of pulses obtained result in sudden and significant variations in substrate concentrations. This does not have any deleterious effect on the system performance because the optimal profile calculated was found using substrate-inhibited kinetics. Thus, the presence of large pulses can result in a more efficient operation arising from a larger rate of reaction on the average. Literature Cited (1) Modak, J. M.; Lim, H. C. Feed back optimization of Fedbatch Fermentation. Biotechnol. Bioeng. 1987, 30, 528. (2) Cazzador, L. On the Optimal Control of Fed-Batch Reactors with Substrate-Inhibited Kinetics. Biotechnol. Bioeng. 1988, 31, 670.
(3) Cuthrell, J. E.; Biegler, L. T. Simultaneous Optimization and Solution Method for Batch Reactor Control Profiles. Comput. Chem. Eng. 1989, 13, 49. (4) Jeawan, L.; Parulekar Satish, H. Enhanced Production of R-Amylase in Fed Batch Cultures of Bacillus Subtilis TN106. Biotechnol. Bioeng. 1993, 42, 1142. (5) Levin, K. L. Maximizing the Product Distribution in Batch Reactors: Reactions in Parallel. Chem. Eng. Sci. 1992, 47, 1751. (6) Weigand, W. A. Maximum Cell Productivity by Repeated Fed-Batch Culture Yield Case. Biotechnol. Bioeng. 1981, 23, 249. (7) Bonte, P.; Modak, J. M.; Tayeb, Y. J.; Lim, H. C. Computational Algorithms for Optimal Fed Rates for a Class of Fed Batch Fermentation: Numerical Results for Penicillin and Cell Mass Production. Biotechnol. Bioeng. 1986, 28, 1408. (8) Modak, J. M.; Lim, H. C.; Jayeb, Y. K. General Characteristics of Optimal Feed Rate Profiles for Various Fed Batch Fermentation Processes. Biotechnol. Bioeng. 1986, 28, 1396. (9) Biegler, L. T. Solution of Dynamic Optimization, Problems by Successive Quadratic Programming and Orthogonal Collocation. Comput. Chem. Eng. 1984, 8, 243. (10) Morrison, K. R.; Sargent, R. W. H. Optimization of Multistage Processes Described by Differential-Algebraic Equations. Lect. Notes Math. 1984, 12, 86. (11) Vassiliadis, V. S.; Sargent, R. W. H.; Pantelides, C. C. Solution of a Class of Multistage Dynamic Optimization Problem 1: Problems without Path Constraints. Ind. Eng. Chem. Res. 1984, 33, 2111. (12) Shukla, P. K.; Pushpavanam, S. Optimization of Biochemical Reactors: An Analysis of Different Approximations of Fed Batch Operation. Chem. Eng. Sci. 1998, 53, 341. (13) Rosen, O.; Luus, R. Evaluation of gradients for piecewise constant optimal control. Comput. Chem. Eng. 1991, 15, 273.
Received for review August 20, 1998 Revised manuscript received January 25, 1999 Accepted January 26, 1999 IE9805483
