3434
Ind. Eng. Chem. Res. 1998, 37, 3434-3443
Fuzzy-Decision-Making Problems of Fuel Ethanol Production Using a Genetically Engineered Yeast Feng-Sheng Wang* and Chang-Huei Jing Department of Chemical Engineering, National Chung Cheng University, Chia-Yi 621, Taiwan, Republic of China
George T. Tsao Laboratory of Renewable Resource Engineering, Purdue University, 1295 Potter Center, West Lafayette, Indiana 47907-1295
A fuzzy-decision-making procedure is applied to find the optimal feed policy of a fed-batch fermentation process for fuel ethanol production using a genetically engineered Saccharomyces yeast 1400 (pLNH33). The policy consisted of feed flow rate, feed concentration, and fermentation time. The recombinant yeast 1400 (pLNH33) can utilize glucose and xylose simultaneously to produce ethanol. However, the parent yeast utilizes glucose only. A partially selective model is used to describe the kinetic behavior of the process. In this study, this partially selective fermentation process is formulated as a general multiple-objective optimal control problem. By using an assigned membership function for each of the objectives, the general multiple-objective optimization problem can be converted into a maximizing decision problem. In order to obtain a global solution, a hybrid method of differential evolution is introduced to solve the maximizing decision problem. A simple guideline is introduced in the interactive programming procedures to find a satisfactory solution to the general multiple-objective optimization problem. 1. Introduction In process design, chemical engineers often need to make decisions when faced with competing objectives. However, the science of decision making with competing objectives has been studied mostly by economists, mathematicians, operation researchers, and others (Tamura and Yoshikawa, 1990; Sawaragi et al., 1985). Until very recently, almost all optimal control problems in chemical process engineering have been treated as problems of single-objective optimization. In addition to the economical efficiency, the performance of many process systems needs to be evaluated with a variety of other criteria (Sophos et al., 1980). For instance, the objective of time optimal control is to find a control policy that drives the given initial states to the desired states in a minimum time. This type of problem occurs often in batch polymerization processes (Thomas and Kiparissides, 1984; Butala, 1990). In such processes, reducing the batch reaction time and achieving the desired polymer properties are simultaneously required. They, therefore, become multiple-objective optimization problems. Even though multiple-objective optimization of polymerization processes has shown promising results (Choi and Butala, 1991; Farber, 1986; Wang and Shieh, 1997), little is known about applying such ideas to biochemical processes. Multiple-objective optimization provides a framework for understanding the relationships between the various objective functions and allows an engineer to make decisions on how to trade-off among the objectives to achieve performance considered “the best”. It is an inherently interactive algorithm, with the engineer constantly making decisions. Much of the decision * Corresponding author. Telephone: +886 5 2428153. Fax: +886 5 2721206. E-mail:
[email protected].
making in the real world takes place in an environment in which the goals, the constraints, and the consequences of possible actions are not known precisely. To deal quantitatively with imprecision, decision making in a fuzzy environment can be applied to handle these imprecise goals and constraints. In this paper we apply the ideas of multiple-objective optimization to the design of an optimal feed policy for a fed-batch fermentation process. Ethanol has been given considerable attention over the years as an octane booster, fuel extender, or neat liquid fuel. Today, there is heightened interest in ethanol as a transportation fuel. Enough ethanol could be made from renewable sources of cellulosic biomass to replace some gasoline consumption. Ethanol production from renewable resources can improve energy security, reduce accumulation of carbon dioxide, and decrease urban air pollution. When blended with gasoline, “neat” ethanol reduces the release of smogforming compounds. Thus, ethanol from lignocellulosic materials holds great promise as a new industry in the world and has the potential for making a significant contribution to the solution of major energy as well as environmental problems (Olsson and Hahn-Hagerdal, 1996; Lee, 1997; Krishnan et al., 1997). Although ethanol production using fermentation of sugar has been studied for many years, there are several bottlenecks for the economical production of fuel ethanol. One of them is ethanol inhibition, which is considered to be the principal factor restricting the fermentation rate and concentration of ethanol achievable in the production process. Lignocellulosic feedstocks like wood, waste paper, and fast-growing energy crops have been identified as economical starting materials for ethanol production. These raw materials contain glucose and xylose as the major fermentable sugars.
S0888-5885(97)00736-7 CCC: $15.00 © 1998 American Chemical Society Published on Web 07/10/1998
Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 3435
Although yeasts can ferment glucose to ethanol readily, the fermentation of xylose to ethanol represents the second bottleneck in the production process. For a process that converts biomass to ethanol to be economically competitive, it is mandatory to convert xylose to ethanol in addition to glucose. Another area which influences the conversion of biomass is the pretreatment of the lignocellulosic feedstock. In order to achieve high yields of fermentable sugars and high product yields from the starting material, an effective pretreatment is essential. A significant breakthrough has been achieved at the Laboratory of Renewable Resources Engineering with the development of a genetically engineered Saccharomyces yeast 1400 (pLNH33). Xylose and glucose can both be utilized as the substrate by the recombinant yeast 1400 (pLHN33). This new yeast has also a high ethanol tolerance. The objective of this paper is to apply fuzzy decision making to determine the optimal feed policy of a fedbatch fermentation process for ethanol production from lignocellulosic sugars using the genetically engineered yeast 1400 (pLNH33). The mathematical models of the yeast 1400 (pLNH33) are described in section 2. The objective functions under the crisp or fuzzy environment are discussed in section 3. In this section, the problem of general multiple-objective optimization is converted into a maximizing decision problem. A hybrid method of differential evolution is then introduced to solve the decision-making problem to obtain a global solution as discussed in section 4. The computational results for this fermentation process are discussed in section 5. The concluding remarks are given in section 6. 2. Problem Formulation 2.1. Kinetic Model. The recombinant Saccharomyces yeast 1400 (pLNH33) was developed by cloning the xylose reductase and xylitol dehydrogenase genes from Pichia stipitis and overexpressing the xylulokinase activity of the host yeast, which was the fusion product yeast 1400 (Ho and Tsao, 1995). The recombinant yeast 1400 (pLNH33) utilizes glucose and xylose simultaneously. The recombinant plasmids in 1400 (pLNH33) may be lost by natural processes during fermentation. However, the parent yeast 1400 utilizes glucose only. Accordingly, Krishnan (1996) has developed a partially selective model to describe the cell growth and product formation of 1400 (pLNH33) on glucose and xylose mixtures. The specific cell growth rates for 1400 (pLNH33) on glucose or xylose are respectively expressed by
µg+ )
µx+ )
µmg+sg
{1 - (pg/pmg)βg} sg2 Kg + sg + Kig µmx+sx
βx
{1 - (px/pmx) } sx2 Kx + sx + Kix
(1)
(2)
The specific cell growth rate for 1400 (pLNH33) on glucose and xylose mixtures is approximated as the one on glucose, i.e., µmix ≈ µg+. The specific production rates for 1400 (pLNH33) on glucose or xylose are respectively in the forms
νg+ )
νx+ )
νmg+sg
{1 - (pg/p′mg)γg} sg2 K′g + sg + K′ig νmx+sx K′x + sx +
sx2
{1 - (px/p′mx)γx}
(3)
(4)
K′ix
Since the parent yeast 1400 on xylose cannot survive, both the specific cell growth and specific production rate on xylose are zero. The specific cell growth and specific production rate for the parent yeast 1400 on glucose are respectively expressed in the forms
µg- )
νg- )
µmg-sg
{1 - (pg/pmg)βg} sg2 Kg + sg + Kig νmg-sg K′g + sg +
sg2
{1 - (pg/p′mg)γg}
(5)
(6)
K′ig
Here the notations in (1)-(6) and their corresponding data are expressed in the Nomenclature section of this paper. 2.2. Fed-Batch Model. Most fermentation processes were operated in batches. Fed-batch operation has been shown to yield improved productivity. Such an operation has been found to be particularly effective for processes where effects such as substrate inhibition, catabolite repression, production inhibition, and auxotrophic mutation are important (Cazzador, 1988; Hong, 1986). These phenomena lead to unimodal reaction rate expressions which exhibit a maximum with respect to a single reactant concentration or in terms of two or more reactant concentrations. In a fed-batch fermentation, the substrate concentration can be maintained at a fairly low level to avoid growth inhibition by high concentrations. Due to this advantage, fed-batch fermentation techniques have been developed for many bioprocesses (Lim et al., 1986; Lee and Ramirez, 1994; Hilaly et al., 1995). In such cases, depending on the cultivation process employed, the optimum feed strategy of the substrate should be investigated so as to obtain the maximum production. The dynamic mass balance equations of the yeast 1400 (pLNH33) on glucose and xylose mixtures in a fedbatch fermenter are expressed as
dx+ F ) µmixx+ - Rµmixx+ - x+ dt V
(7)
When glucose is vanished in the broth, this cell growth equation becomes
dx+ F ) µx+x+ - x+ dt V
(8)
Here the first term in (7) is the rate of cell growth in the substrate mixtures. The second term denotes the rate of plasmid loss. When a cell of 1400 (pLNH33) reverts to its parent yeast, it cannot survive on xylose so that the probability of plasmid loss R is reduced to
3436 Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998
zero. The cell balance equation for the parent yeast on glucose and xylose mixtures is therefore expressed as
dxF ) µg-x- + Rµmixx+ - xdt V
pg(t) V(t) - pg(0) V(0) (V(t) - V(0))sgF + V(0) sg(0) - V(t) sg(t)
- Yp/g e 0 (20)
(9)
The ethanol formation from glucose and xylose mixtures is expressed as
dp dpg dpx ) + dt dt dt
g4 )
(10)
where dpg/dt denotes the rate of ethanol formation from glucose and dpx/dt the rate of ethanol formation from xylose. These two rate equations are respectively expressed as
dpg F ) νg+x+ + νg-x- - pg dt V
(11)
dpx F ) νx+x+ - px dt V
(12)
g5 )
px(t) V(t) - px(0) V(0) (V(t) - V(0))sxF + V(0) sx(0) - V(t) sx(t)
- Yp/x e 0 (21)
The theoretical stoichiometric yield coefficient for ethanol formation from glucose or xylose is 0.51. However, the experimental coefficient is Yp/g ) 0.45 for ethanol from glucose and Yp/x ) 0.35 for ethanol from xylose. If the constraints in (20) and (21) are not included in the optimization problem, unrealistic predicted values of p(t) may be found. This fact has been observed in ethanol production processes (Fu and Barford, 1993; Wang and Shyu, 1997). The concentrations of glucose and xylose at the end of the fermentation process have to be limited to avoid possible adverse effects on downstream product separation. We have therefore
g6 ) sg(tf) - sgmin e 0
(22)
Rates of glucose and xylose consumption are respectively expressed as
g7 ) sx(tf) - sxmin e 0
(23)
dsg F 1 (ν x + νg-x-) + (sgF - sg) )dt Yp/g g+ + V
The lower bounded levels, sgmin and sxmin, will help reduce the separation cost in downstream processing.
(13)
3. Multiple-Objective Decision-Making Problem
dsx 1 F )ν x + (s - sx) dt Yp/x x+ + V xF )-
1 F ν x + (λs - sx) Yp/x x+ + V gF
(14)
where the constant λ is a ratio of feed concentration of xylose and glucose in the mixture. The working volume is governed by a total mass balance in the form
dV/dt ) F(t)
(15)
where F(t) is the inlet flow rate of glucose and xylose mixtures. 2.3. System Constraints. Nearly all engineering processes will have physical constraints. In this study, the flow rate is bounded and the volume of the bioreactor is constrained, i.e.
0 e F(t) e Fmax
(16)
g1 ) V(t) - Vf e 0
(17)
The concentration of glucose and xylose must be positive for all time; otherwise, an unrealistic solution in the optimization problem would be obtained. We have therefore
g2 ) -sg(t) e 0
(18)
g3 ) -sx(t) e 0
(19)
In addition, the stoichiometry of the ethanol formation from glucose or xylose must be obeyed, posing two constraints given as
3.1. Decision Making under a Crisp Environment. A production planning problem is considered in this study such that the decision maker (DM) designs a feed control policy for the fermentation of ethanol from glucose and xylose mixtures using the recombinant yeast 1400 (pLNH33). The objective of the problem is to find optimal feed flow rate, feed concentration, and fermentation time of the fed-batch process such that the ethanol production should be greater than or equal to some threshold value, the consumption of glucose or xylose should be less than or equal to some threshold value, and the fermentation time should also be less than or equal to some threshold value. According to these statements, the production planning problem becomes a multiple-objective decision-making problem. Two requirements must be satisfied in such a decisionmaking problem. The first requirement is to find the optimal feed flow rate, feed concentration, fermentation time, and the associated objective function values. Such an optimal solution can be obtained by multiple-objective optimization techniques. However, the second requirement is to check whether or not the optimal solution will satisfy the pre-assigned threshold values. If the optimal solution does not satisfy the threshold values, the DM has to trade-off some threshold values. The effort should be repeated to find another optimal solution. According to the above-mentioned procedures, the first requirement is, therefore, expressed as the following multiple-objective optimal control and optimal parameter selection problem. This problem is simply called the multiple-objective optimization problem (MOOP) and is expressed as
max J1 ) p(tf) V(tf)
F(t),sgF,tf
(24)
Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 3437
min J2 ) sgF(V(tf) - V(0)) + sg(0) V(0)
Ji(u) e Ji(u*)
i ) 1, ..., 4
for each of the objective functions to convert a multipleobjective optimization problem into a single-objective function problem. For example, the unit price of ethanol and substrates can be used as the weighting factor. This weighted single-objective function becomes the problem of maximum profit. Such a problem can be solved by a conventional optimization technique. The optimal solution of the weighted sum problem is one of the Pareto solutions to the multiple-objective optimization problems. However, such a solution may be a local solution due to a duality gap between the solutions of the weighted single-objective and multiple-objective optimization problems. Clearly, if the Pareto surface is nonconvex, the weighted sum method may yield poor designs no matter what weight or optimization method is used. Several methods have been proposed to overcome the drawback of such a nonconvex problem (Sawaragi et al., 1985; Tamura and Yoshikawa, 1990). After the optimal solution is obtained from a multipleobjective optimization technique, the second requirement in this decision-making problem is then performed to check whether or not the optimal solution satisfies the assigned threshold values. If the optimal solution does not satisfy the threshold values, the DM has to assign another threshold requirement. The problem should then be repeated to find another optimal solution. Interactive programming can be employed to solve the decision-making problem. In this study, the interactive fuzzy optimization is extended to solve the multiple-objective optimal control and optimal parameter selection problem. 3.2. Fuzzy-Decision-Making Problems. So far we have considered the decision-making problem under the crisp environment; that is, the optimal solution must absolutely satisfy the assigned threshold values. In this section, we will discuss the multiple-objective decisionmaking problem under the fuzzy environment. Assume that the DM has fuzzy goals for each of the objective functions in (24)-(27). The fuzzy goal means an interval of the assigned threshold instead of a point value in a crisp environment. As a result, the DM considers the fuzzy objective function such as J1 should be substantially greater than or equal to a threshold interval [JL1 , JU 1 ]. The second, third, and fourth goals should be substantially less than or equal to the assigned threshold interval [JLk , JU k ], k ) 2, 3, 4. The multiple-objective optimization problem (24)-(27) is now extended to the general multiple-objective optimization problem (GMOOP) given as
Jk(u) < Jk(u*)
for some k
fuzzy max J1 ) p(tf) V(tf)
F(t),sgF,tf
(25)
min J3 ) λsgF(V(tf) - V(0)) + sx(0) V(0) (26)
F(t),sgF,tf
min J4 ) tf
F(t),sgF,tf
(27)
The first objective function corresponds to the total price of ethanol production. The second and third objective functions are the cost of the substrates. The last objective function corresponds to the operation cost. Multiple-objective optimization is a natural extension of the traditional optimization of a single-objective function. If the multiple-objective functions are commensurate, minimizing with respect to one objective function will minimize with respect to all criteria and the problem can be solved using traditional optimization techniques. However, if the objective functions are incommensurate, or competing, then the minimization of one objective function requires a compromise in another objective function. The competition between multiple-objective functions gives rise to the distinction between multiple-objective optimization and traditional single-objective optimization. The problem is further complicated by the lack of a complete order for multiple objectives. The concept of Pareto optimality or noninferiority is, therefore, used to characterize a solution to the multiple-objective optimization problem. In order to concisely define the Pareto optimal solution, we introduce the following definitions. Definition 1. The feasible region in input space, Ω, is the set of all admissible control variables and the system parameters that satisfy the system constraints, i.e.
Ω ) {u ≡ [F(t), sgF, tf]T|z3 ) f(z,u), z(0) ) z0; Fmin e F(t) e Fmax; gk(z,u) e 0, k ) 1, ..., 7} Here the state equation, z3 ) f(z,u), consists of the fedbatch model in (7)-(15). The decision variable u consists of the feed flow rate, feed concentration, and fermentation time. We are now in a position to define Pareto optimal solutions to the combined optimal control and optimal parameter selection problem. Definition 2. A control action u* is a Pareto optimal policy if and only if there does not exist u ∈ Ω such that
The image of a Pareto optimal policy is a Pareto optimal solution. In general, there are an infinite number of Pareto policies for a given multiple-objective optimization problem. The collection of Pareto policies is the Pareto set. The image of this set is called the trade-off surface. The literature on multiple-objective optimization is huge, and we cannot hope to mention all the techniques that have been used to generate Pareto solutions. However, the weighted sum method in multiple-objective optimization textbooks is pervasive (Sawaragi et al., 1985; Tamura and Yoshikawa, 1990; Sakawa, 1993). This is one of the most commonly used techniques for solving problems of multiple-objective optimization of chemical processes. The DM assigns a weighting factor
F(t),sgF,tf
(28)
fuzzy min J2 ) sgF(V(tf) - V(0)) + sg(0) V(0) (29) F(t),sgF,tf
fuzzy min J3 ) λsgF(V(tf) - V(0)) + sx(0) V(0) (30) F(t),sgF,tf
fuzzy min J4 ) tf F(t),sgF,tf
(31)
The fuzzy requirement for each of all objective functions can be quantified by eliciting membership functions from the DM. In maximization, a fuzzy goal stated by the DM may be to achieve “substantially greater than or equal to some interval”, and the DM is asked to determine the subjective membership function which is a strictly monotonically decreasing function with respect
3438 Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998
to J1(u) in the following way:
{
1; JU 1 e J1 µ1(J1) ) d1; JL1 e J1 e JU 1 0; J1 e JL1
(32)
where JL1 or JU 1 represents the value of J1 such that the grade of the membership function µ1(J1) is 0 or 1 and the grades of the membership for the intermediate function values are expressed by a strictly monotonically decreasing function d1 with respect to J1. Following a similar procedure, the membership functions for minimizing goals of (29)-(31) are expressed as
{
1; JLk e Jk k ) 2, 3, 4 L U d ; µk(Jk) ) k Jk e Jk e Jk 0; JU k e Jk
(33)
where JLk or JU k represents the value of Jk such that the grade of the membership function µk(Jk) is 1 or 0 and the grades of the membership for the intermediate function values are expressed by a strictly monotonically decreasing function dk with respect to Jk. The membership function in (32) or (33) becomes the characteristic function if dk is set to 0 or 1. As a result, the fuzzy problems (28)-(31) become the crisp problems in (24)(27). Having elicited the membership functions from the DM for each of the objective functions, the GMOOP (28)-(31) can be converted into the fuzzy multipleobjective optimization problem (FMOOP) by
min [µ1(J1) µ2(J2) µ3(J3) µ4(J4)]T u∈Ω
(34)
By introducing a general aggregation function µD(µk(Jk)), a fuzzy multiple-objective decision-making problem (FMODMP) or maximizing decision problem can be defined by
max µD u∈Ω
(35)
Several aggregation functions have been used in a standard fuzzy nonlinear programming problem as discussed in the textbook of Sakawa (1993). In this study, the fuzzy decision or minimum operator of Bellman and Zadeh (1970) is selected as the aggregation function by
µD ) min {µk(Jk), k ) 1, ..., 4} k
Observe that the value of the aggregation function can be interpreted as representing an overall degree of satisfaction with the DM’s multiple fuzzy goals. Let us consider the fuzzy maximizing problem. While the objective function value is greater than the assigned upper bound, such a solution absolutely satisfies the DM. On the other hand, the objective function value is less than the lower bound. It must be rejected. While the objective function value is located between the threshold interval, the DM has some degree of satisfaction for the solution. Fundamental to the MOOP in (24)-(27) is the Pareto optimal concept, and thus the DM must select a compromise solution from many Pareto optimal solutions. The relationships between the optimal solutions
of the FMODMP and the Pareto optimal concept of the MOOP can be characterized by the following theorem. Theorem 1. If u* is a unique optimal solution to the FMODMP in (35), then u* is a Pareto optimal solution to the MOOP in (24)-(27). The proof of this theorem follows immediately from Definition 2 of the Pareto optimality by making use of contradiction arguments. This theorem is used to guarantee that the unique optimal solution of the FMODMP is a Pareto solution to the crisp multipleobjective optimal control problems (24)-(27). The statement of this theorem does not guarantee that the unique optimal solution to (35) is a Parateo solution to the GMOOP (28)-(31). Sakawa (1993) introduced the concept of fuzzy Pareto or M-Pareto optimal solutions for the general multiple-objective nonlinear programming problems. Such a definition can be extended to the combined optimal control and optimal parameter selection problem in this study. This is defined in terms of membership functions instead of the objective functions. Definition 3. u* ∈ Ω is said to be a M-Pareto optimal solution to GMOOP if and only if there does not exist another u ∈ Ω such that µk(Jk(u)) g µk(Jk(u*)) for all k and µj(Jj(u)) * µj(Jj(u*)) for at least one j. Note that the set of Pareto optimal solutions is a subset of the set of M-Pareto optimal solutions as observed from Definitions 2 and 3 and (32). Here M refers to membership. Using the concept of M-Pareto optimality, the fuzzy version of Theorem 1 can be obtained under slightly different conditions. Theorem 2. If u* is a unique optimal solution to the FMODMP (35), then u* is an M-Pareto optimal solution to the GMOOP (28)-(31). The proof of this theorem follows immediately from Definition 3 of the M-Pareto optimality by making use of contradiction arguments. Theorem 2 is used to guarantee that the unique optimal solution of the maximizing decision problem (35) is a M-Pareto optimal solution of the fuzzy problems (28)-(31). The key point for using this theorem is to find a unique optimal solution of the problem (35). A global optimization method must be employed to determine such a unique solution. A hybrid method of differential evolution, one of global optimization methods, will be introduced in the next section to solve the maximizing decision problem (35). Interactive programming techniques are tools for searching a satisfactory solution by interaction between the DM and a computer. It can be regarded as an interface between humans and computers. An interactive programming algorithm is introduced in this study and is listed in the following to find a satisfactory solution to the GMOOP. r 1. Assign the threshold intervals [JLk JU k] . 2. Elicit a membership function µk(Jk) from the DM for each of the objective functions. 3. Solve the maximizing decision problem (35). Let (ur, µrk(Jk)) be the M-Pareto optimal solution to the GMOOP. 4. If the DM is satisfied with the current levels of µrk(Jk), k ) 1, ..., 4, the current M-Pareto optimal solution (ur, µrk(Jk)) is the satisfactory solution for the DM. Otherwise, classify the objectives into three groups based on the DM’s preference, including (a) a class of the objectives that the DM wants to improve, (b) a class
Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 3439
of the objectives that the DM may possibly agree to relaxing, and (c) a class of the objectives that the DM accepts. The index set of each class is represented by Ir, Rr, and Ar, respectively. The new threshold intervals r+1 is reassigned in such a way that [JLk JU k] L U r r+1 for any k ∈ Ir, [JL JU]r+1 ⊂ [Jk Jk ] ⊂ [JLk JU k] k k L U r r+1 ) [JL JU]r for [Jk Jk ] for any k ∈ Rr, and [JLk JU k] k k any k ∈ Ar. Then repeat to step 2. Here, it should be stressed that any improvement for one of the objective functions can be achieved only at the expense of at least one of the other objective functions. 4. Hybrid Differential Evolution Since the feed flow rate F(t) as a function of time is to be chosen in this optimal control problem, such a problem is an infinite dimensional problem. To solve this problem efficiently, the feed flow rate is first represented by a finite set of control parameters F(j) in the time interval tj-1 e t < tj as
F(t) ) F(j),
j ) 1, ..., Nu
(36)
where Nu is the number of time partitions. The infinite dimensional problem (35) is then approximated as the parameter selection problem. Several optimization methods can be used to solve such a problem. In this study, a hybrid method of differential evolution (DE) is introduced to find a global solution. Differential evolution is a very simple populationbased stochastic function minimizer which is introduced by Storn and Price (1996). This method turned out to be the best genetic type of algorithm for solving the realvalued test function suite of the first International Contest on Evolutionary Computation which was held in Nagoya, Japan, in 1996. The original version of differential evolution is used to solve the unconstrained nonlinear programming problems. Wang and Chiou (1997) have extended DE to optimal control problems. The main structure of DE is given as follows: 1. Initialization. 2. Mutation operation. 3. Crossover operation. 4. Evaluation operation. 5. Repeat steps 2-4. This structure is a parallel direct search algorithm which utilizes Np vectors of the decision parameters u, which consists of the control parameters, the feed concentration, and the fermentation time, in the optimization problem, i.e., UG i , i ) 1, 2, ..., Np, as a population for each generation G. The initial population is chosen randomly and should try to cover the entire search space uniformly in the form
U0i
) umin + Fi(umax - umin),
i ) 1, ..., Np
(37)
where Fi is the uniformly distributed random number. Here, umin and umax are respectively expressed as the lower bound and upper bound of the decision parameters. The three main key operations in DE are briefly discussed in turn. The essential ingredient in the mutation operation is the difference vector. The mutation process at the (G - 1)th generation begins by , UG-1 , randomly selecting two or four individuals UG-1 j k
, and UG-1 for any j, k, l, and m. These four UG-1 l m individuals are then combined to form a difference vector Djklm as
Djklm ) Djk + Dlm ) (UG-1 - UG-1 )+ j k - UG-1 (UG-1 l m ) (38) A perturbed individual U ˆ G-1 is therefore generated i in the mutation based on the parent individual UG-1 p process by G-1 U ˆG + λmDjklm, i ) Up
i ) 1, ..., Np
(39)
where a scaling factor λm ∈ (0, 1.2] is introduced by Storn and Price (1996) in (39) to ensure the fastest possible convergence. In essence, the resulting indi. Here, the vidual in (39) is a noisy replica of UG-1 p is dependent on which strategy parent individual UG-1 p of the mutation operators is employed. Five strategies in the original version of DE have been introduced by Storn and Price (1996). In order to increase the diversity of the new individuals at the next generation, the perturbed individual U ˆG i G-1 and the current individual Ui are chosen by a binomial distribution to perform the crossover operation to generate the offspring. In this crossover operation, the jth gene of the ith individual at the next generation is ˆ 1, ..., produced from the perturbed individual U ˆG i ) [u G G-1 ) [u1, ..., un]G-1 uˆ n]i and the current individual Ui i as follows:
uG ji )
{
if a random number > CR uG-1 ji , G uˆ ji , otherwise, j ) 1, ..., n; i ) 1, ..., Np (40)
where the crossover factor CR ∈ [0, 1] is set by the user. The evaluation function of an offspring is one-to-one compared to that of its parent in DE. This competition means that the parent is replaced by its offspring if the fitness of the offspring is better than that of its parent. On the other hand, the parent is retained in the next generation if the fitness of the offspring is worse than that of its parent. Two selection steps are performed in this evaluation operation. The first selection step is a one-to-one competition. The next step is to select the best individual in the population. These steps are expressed as G-1 ), J(UG UG i ) arg min{J(Ui i )}, i ) 1, ..., Np G U ˆG b ) arg min{J(Ui ), i ) 1, ..., Np}
(41) (42)
The above four steps are repeated until the maximum iterations or the desired fitness is achieved. The extended computational algorithms to solve the optimal parameter selection problems have been shown by Wang and Chiou (1997) in detail. When using DE to optimize a function, one generally has to find a good trade-off between convergence and diversity. The convergence means fast convergence even to a local minimum. On the other hand, the higher diversity guarantees a higher probability of finding the global minimum. This trade-off can be achieved by slightly adjusting the scaling factor λm of the mutation operator and the population size Np. To increase Np and
3440 Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998
simultaneously reduce λm slightly, differential evolution is getting more robust. By doing so, a lot of computational time should be expended to evaluate the fitness function. This fact is particularly serious by using DE to solve optimal control problems or optimal parameter selection problems due to expending extensive CPU time to solve differential equations (Wang and Chiou, 1997). While the population diversity is small, the candidate individuals are gathered together so that the perturbed individual cannot be generated any longer by the mutation and crossover operations as observed from (39) and (40). This fact results in a premature convergence because the difference vector has diminished. Chiou and Wang (1997) have proposed a hybrid method of DE for parameter estimation of a bioreaction model to overcome such drawbacks. The main key operations of the hybrid differential evolution (HDE) are shown as follows: 1. Initialization. 2. Mutation operation. 3. Crossover operation. 4. Evaluation operation. 5. Acceleration phase if necessary. 6. Migrating phase if necessary. 7. Repeat steps 2-6. An accelerated phase and a migrating phase are embedded into the original version of DE. These two phases act as a trade-off operator. The accelerated phase is used to speed up convergence. However, faster convergence usually results in finding a local minimum. The migrating phase is used to escape this local point. A new population of candidate individuals are regenerated based on the best individual at the current generation. Accordingly, the diversity of the candidates can still be retained by such regenerated procedures. According to our experience by using DE to solve optimization problems, the best fitness does not descend continuously from generation to generation. It usually descends toward a better next one after several generations. In this situation, the accelerated operation can be used to overcome this drawback. When the best fitness at the present generation is not improved any longer by the mutation and crossover operations, a descent method is then applied to push the best individual toward obtaining a better point. The accelerated operation is, therefore, expressed as
UG b )
{
U ˆG if J(U ˆG ˆ G-1 ) b, b ) < J(U b G U ˆ b - R∇J, otherwise
(43)
where U ˆG b is the best individual obtained from (42). The gradient of the objective function, ∇J, can be approximately calculated with finite difference. The step size R ∈ (0, 1] in (43) is determined with the descent property. At first, R is set to 1 to obtain the new N individual UN b . The fitness function J(Ub ) is then G compared with J(Ub ). If the descent property is obeyed, i.e., G J(UN b ) < J(Ub )
(44)
then UN b , which becomes a candidate in the next generation, is added into this population to replace the worst individual. On the other hand, if the descent property fails, the step size is reduced, say, to 0.5 or 0.8, and the descent method is repeated to find UN b
until R∇J is small enough or has achieved some iterations. As a result, the best fitness should be at least ) as observed from (44). equal or smaller than J(UG-1 b The rate of convergence can be improved by the accelerated operation. However, faster descent usually results in finding a local minimum or getting premature convergence. In addition, performing this operation frequently can result in the fact that the candidate individuals are gradually gathered around the best individual so that the diversity of the population is decreased quickly. Furthermore, the gathered individuals cannot reproduce the next better individuals by the mutation and crossover operations. As a result of this matter, the migrating phase has to be performed to regenerate a new population. The new candidates are regenerated based on the best individual UG b . The jth gene of the ith individual is therefore regenerated by
{
ujb - ujmin ujmax - ujmin uji ) ujb + δji(ujmax - ujb), otherwise (45) ujb + δji(ujmin - ujb), if δˆ ji
2 0, otherwise
where 1 and 2 are respectively expressed as the desired tolerance for the population diversity and the gene diversity with respect to the best individual. Here, ηji is defined as an index of gene diversity. Its value of zero denotes that the jth gene of the ith individual is close to the jth gene of the best individual. 5. Results and Discussion All computations were performed on a Pentium Pro 200 computer using a Microsoft Windows NT version 4.0 operating system. The HDE algorithm was written by a Microsoft Fortran PowerStation version 4.0. The HDE source code is available free from the corresponding author. Since the physical constraints in (17)-(23) are included in the optimization problem, the penalty function method is used to handle the system constraints in HDE. The fitness function used in HDE is, therefore, defined as 5
max J ) µD -
F(j),sgF,tf
γi∫0 〈gi〉+2 dt - γ6〈g6〉+2 - γ7〈g7〉+2 ∑ i)1 tf
(47)
Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 3441
Figure 1. Assigned membership function for each of the objective functions. Table 1. Assigned Threshold Intervals and the Optimal Membership Function Value for Each of the Objective Functions in the First and Second Run First Run objective 1st 2nd 3rd 4th L U [80 kg 90 kg] [100 kg 140 kg] [60 kg 80 kg] [40 h 50 h] [Jk Jk ] 0.0 0.6997 1.0 0.8930 µ/k (Jk)
Figure 2. Optimal feed flow rate in the second run.
Second Run objective 1st 2nd 3rd 4th [70 kg 90 kg] [120 kg 140 kg] [60 kg 80 kg] [40 h 50 h] [JLk JU k] 0.5798 0.5761 0.5987 0.8331 µ/k (Jk)
The integration of the square penalty functions in (47) is used to cover the state variables on the whole time domain. The linear membership functions are used in this simulated example. The membership function for each of the objective functions is described in Figure 1. The DM is first assigned the threshold intervals as shown in the first run of Table 1. Now, the maximizing decision problem can be solved by HDE. By choosing 20 time partitions for the feed flow rate, 22 decision parameters, i.e., 20 control parameters, the feed concentration, and the fermentation time, have to be determined in the finite-dimensional optimization problem. In this work, population size Np ) 20, scaling factor λm ) 0.5, crossover factor CR ) 0.5, tolerance of population diversity 1 ) 0.3, and tolerance of gene diversity 2 ) 0.3 were used in HDE. The maximizing decision µD for these threshold requirements was zero after 1000 generations. The membership function values for each of the objectives are shown in Table 1. The ethanol product of 71.00 kg was less than the lower bound of the assigned threshold requirement so that its corresponding membership function value was zero. The minimum supplied amounts of glucose and xylose were 112.01 kg and 59.16 kg, respectively. The fermentation time was 41.07 h in this run. These three minimum values were within the desired requirements. Since the maximizing decision was zero, a trade-off procedure should be performed to find a satisfied solution. The DM wants to improve the degree of satisfaction for ethanol production so that the lower bound of the threshold intervals was reduced to 70.00 kg. According to the concept of multiple-objective optimization, any improvement of one objective function value could be achieved only at the expense of at least one of the other objective functions. Assume the third and fourth objectives were accepted by the DM so that the threshold intervals were unchanged. The lower bound of the desired glucose consumption was increased to 120 kg such that the degree of satisfaction for the glucose consumption should be decreased. According to these four requirements, the maximizing decision prob-
Figure 3. Concentration profiles of yeast 1400 (pLNH33) and its parent yeast in the second run.
lem (35) was then solved again. By using the same setting factors as in the aforementioned run, the maximizing decision µD for these threshold requirements was 0.5760 after 1000 generations. In total, 38 migrating operations and 235 accelerated operations were performed in those generations. The value of µD meant that 57.60% of satisfaction was achieved by each of the assigned requirements. The membership function values for each of the objectives are shown in Table 1. The maximum ethanol product of 81.60 kg was obtained. This amount was greater than that of the aforementioned run. The minimum supplied amounts for glucose and xylose were 128.48 and 68.03 kg, respectively. The fermentation time was 41.67 h in this run. The optimal feed flow rate for this case is shown in Figure 2. The optimal feed concentrations for glucose and xylose were 155.60 and 83.78 g/L, respectively. The cell concentration profiles for yeast 1400 (pLNH33) and its parent culture are shown in Figure 3. The concentration profiles for glucose and xylose are represented in Figure 4. The concentration of ethanol is shown in Figure 5. The ethanol production from glucose or from xylose is also illustrated in Figure 5. The concentration of glucose was nearly zero after t > 15 h, as observed in
3442 Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998
6. Conclusions
Figure 4. Concentration profiles of glucose and xylose in the second run.
The fuel ethanol production planning problem using a genetically engineered yeast has been discussed in this study. Such a production planning problem was formulated into a framework of the general multipleobjective optimization problem. Many of the multipleobjective optimization problems in the real world take place in an environment in which the goals, the constraints, and the consequences of possible actions are not known precisely. To deal quantitatively with imprecision, the problem in a fuzzy environment has been introduced in this study to handle these imprecise goals and constraints. Such fuzzy multiple-objective optimal control problems were converted into a maximizing decision problem through the subjective membership functions for each of the objective functions. The optimal solution for each of the membership functions was denoted as the degree of satisfaction with the assigned threshold requirements. In order to obtain a global optimal solution, a hybrid method of differential evolution was introduced to solve the maximizing decision problem. A simple guideline was introduced in the interactive programming procedures to find a satisfactory solution to the general multiple-objective optimization problem. Acknowledgment This work was supported by the National Science Council of Republic of China under Grant NSC86-2214E194-004. Nomenclature
Figure 5. Concentration profiles of ethanol in the second run.
Figure 4. However, the concentration of xylose approached a maximum value. Yeast 1400 (pLNH33) and its parent culture could utilize glucose to product ethanol so that both types of cells grew exponentially before 15 h, as observed in Figure 3. After 15 h, yeast 1400 (pLNH33) could still ferment xylose to ethanol, as observed in Figure 5. However, the parent yeast could not survive on xylose so its growth ceased. The original version of differential evolution was also used to solve the problem for comparison. By using the same setting factors and the threshold requirements as in the second case except switching off the accelerated and migrating phases, an infeasible solution was obtained after 1000 generations. Therefore, the population size was increase to 40. DE with this population size was performed again. The maximizing decision µD of 0.5537 was obtained after 1000 generations. The maximum ethanol product of 81.08 kg was obtained. This amount was little less than that obtained from HDE. The minimum supplied amounts of glucose and xylose were 127.87 and 67.70 kg, respectively. The fermentation time was 41.62 h in this case. The optimal feed concentrations of glucose and xylose were 154.84 and 83.37 g/L, respectively.
F(t) ) feed flow rate (L/h) Fmax ) maximum feed flow rate (200 L/h) Kg ) saturation coefficient for cell growth on glucose (0.565 g/L) K′g ) saturation coefficient for ethanol production on glucose (1.342 g/L) Kig ) inhibition coefficient for cell growth on glucose (283.7 g/L) K′ig ) inhibition coefficient for ethanol production on glucose (4890.0 g/L) Kix ) inhibition coefficient for cell growth on xylose (18.1 g/L) K′ix ) inhibition coefficient for ethanol production on xylose (22.8 g/L) Kx ) saturation coefficient for cell growth on xylose (3.4 g/L) K′x ) saturation coefficient for ethanol production on xylose (3.4 g/L) p(t) ) concentration of ethanol (g/L) pg(t) ) concentration of ethanol from glucose (g/L) pmg ) maximum ethanol concentration for cell growth on glucose (129.9 g/L) p′mg ) maximum ethanol concentration for ethanol production on glucose (136.4 g/L) pmx ) maximum ethanol concentration for cell growth on xylose (59.04 g/L) p′mx ) maximum ethanol concentration for ethanol production on xylose (60.2 g/L) px(t) ) concentration of ethanol from xylose (g/L) sg(t) ) concentration of glucose (g/L) sgF ) feed concentration of glucose (g/L) sgmin ) lower bound of residual glucose (0.5 g/L) sx(t) ) concentration of xylose (g/L) sxF ) feed concentration of xylose (g/L) sxmin ) lower bound of residual xylose (0.5 g/L)
Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 3443 V(t) ) working volume (L) Vf ) total working volume (1000 L) x+(t) ) concentration of the recombinant yeast 1400 (pLNH33) (g/L) x-(t) ) concentration of the parent yeast 1400 (g/L) Yp/g ) yield coefficient for ethanol from glucose (0.45) Yp/x ) yield coefficient for ethanol from xylose (0.35) Abbreviation DE ) differential evolution DM ) decision maker FMODMP ) fuzzy multiple-objective decision-making problem FMOOP ) fuzzy multiple-objective optimization problem GMOOP ) general multiple-objective optimization problem HDE ) hybrid differential evolution MOOP ) multiple-objective optimization problem Greek Letters R ) probability of plasmid loss (0.075) βg ) power of ethanol inhibition for cell growth on glucose (1.29) βx ) power of ethanol inhibition for cell growth on xylose (1.94) γg ) power of ethanol inhibition for ethanol production on glucose (1.42) γx ) power of ethanol inhibition for ethanol production on xylose (0.68) λ ) ratio for feed concentration of xylose and glucose in the feed mixture (0.35/0.65) µD ) aggregation function µg+ ) specific cell growth rate for yeast 1400 (pLNH33) on glucose µg- ) specific cell growth rate for yeast 1400 on glucose µi(Ji) ) membership function for each of the objective functions µmg+ ) maximum specific growth rate coefficient for yeast 1400 (pLNH33) on glucose (0.6) µmg- ) maximum specific growth rate coefficient for yeast 1400 on glucose (0.662) µmx+ ) maximum specific growth rate coefficient for yeast 1400 (pLNH33) on xylose (0.19) µx+ ) specific cell growth rate for yeast 1400 (pLNH33) on xylose νg+ ) specific ethanol production rate for yeast 1400 (pLNH33) on glucose νg- ) specific ethanol production rate for yeast 1400 on glucose νmg+ ) coefficient of maximum specific ethanol production rate for yeast 1400 (pLNH33) on glucose (2.0 h-1) νmg- ) coefficient of maximum specific ethanol production rate for yeast 1400 on glucose (2.005 h-1) νmx+ ) coefficient of maximum specific ethanol production rate for yeast 1400 (pLNH33) on xylose (0.25 h-1) νx+ ) specific ethanol production rate for yeast 1400 (pLNH33) on xylose
Literature Cited Bellman, R. E.; Zadeh, L. A. Decision Making in a Fuzzy Environment. Manage. Sci. 1970, 17, B141. Butala, D. Modeling and Optimal Control of Polymerization Reactors. Ph.D. Dissertation, University of Maryland, College Park, MD, 1990. Cazzador, L. On the Optimal Control of Fed-batch Reactors with Substrate-Inhibited Kinetics. Biotechnol. Bioeng. 1988, 31, 670. Chiou, J. P.; Wang, F. S. Hybrid Differential Evolution for Parameter Estimation of a Batch Bioprocess. IEEE Singapore International Symposium on Control Theory and Applications, Singapore, 1997; p 171.
Choi, K. Y.; Butala, D. N. An Experimental Study of Multiobjective Dynamic Optimization of a Semibatch Copolymerization Process. Polym. Eng. Sci. 1991, 31, 353. Farber, J. N. Steady State Multiobjective Optimization of Continuous Copolymerization Reactors. Polym. Eng. Sci. 1986, 26, 499. Fu, P. C.; Barford, J. P. Non-singular Optimal Control for Fedbatch Fermentation Processes with a Differential-Algebraic System Model. J. Proc. Control 1993, 3, 211. Hilaly, A. K.; Karim, M. N.; Linden, J. C. A Study on Real-Time Optimization of a Fed-batch Recombinant Escherichia coli Fermentation. Control Eng. Pract. 1995, 3, 485. Ho, N. M. Y.; Tsao, G. T. Recombinant yeasts for effective fermentation of glucose and xylose. U.S. PCT Patent No. W095/ 13362, 1995. Hong, J. Optimal Substrate Feed Policy for a Fed-batch Fermentation with Substrate and Product Inhibition Kinetics. Biotechnol. Bioeng. 1986, 28, 1421. Krishnan, M. S. Process Development of Fuel Ethanol Production from Lignocellulosic Sugars Using Genetically Engineered Yeasts. Ph.D. Dissertation, Purdue University, 1996. Krishnan, M. S.; Xia, Y.; Ho, N. W. Y.; Tsao, G. T. Fuel ethanol production from lignocellulosic sugars: studies using genetically engineered yeasts. Fuels and Chemicals from Biomass; ACS Symposium Series 666; American Chemical Society: Washington, DC, 1997; p 74. Lee, J. Biological Conversion of Lignocellulosic Biomass to Ethanol. J. Biotechnol. 1997, 56, 1. Lee, J.; Ramirez, W. F. Optimal Fed-batch Control of Induced Foreign Protein Production by Recombinant Bacteria. AIChE J. 1994, 40, 899. Lim, H. C.; Tayeb, Y. J.; Modak, J. M.; Bonte, P. Computational Algorithms for Optimal Feed Rates for a Class of Fed-batch Fermentation: Numerical Results for Penicillin and Cell Mass Production. Biotechnol. Bioeng. 1986, 28, 1408. Olsson, L.; Hahn-Hagerdal, B. Fermentation of Lignocellulosic Hydrolysates for Ethanol Production. Enzyme Microb. Technol. 1996, 18, 312. Sakawa, M. Fuzzy Sets and Interactive Multiobjective Optimization; Plenum Press: New York, 1993. Sawaragi, Y.; Nakayama, H.; Tanino, T. Theory of Multiobjective Optimization; Academic Press: London, 1985. Sophos, A.; Rotstein, E.; Stephanopoulos, G. Multiobjective Analysis in Modeling the Petrochemical Industry. Chem. Eng. Sci. 1980, 35, 2415. Storn, R.; Price, K. V. Minimizing the Real Function of the ICEC’96 Contest by Differential Evolution. IEEE Conference on Evolutionary Computation, Nagoya, Japan, 1996; p 842. Tamura, H.; Yoshikawa, T. Large-Scale Systems Control and Decision Making; Marcel Dekker: New York, 1990. Thomas, I. M.; Kiparissides, C. Computation of the Near-Optimal Temperature and Initiator Policies for a Batch Polymerization Reactor. Can. J. Chem. Eng. 1984, 62, 284. Wang, F. S.; Chiou, J. P. Optimal Control and Optimal Time Location Problems of Differential-Algebraic Systems by Differential Evolution. Ind. Eng. Chem. Res. 1997, 36, 5348. Wang, F. S.; Shieh, T. L. Extension of Iterative Dynamic Programming to Multiobjective Optimal Control Problems. Ind. Eng. Chem. Res. 1997, 36, 2279. Wang, F. S.; Shyu, C. H. Optimal Feed Policy for Fed-batch Fermentation of Ethanol Production by Zymomous mobilis. Bioprocess Eng. 1997, 17, 63.
Received for review October 22, 1997 Revised manuscript received May 8, 1998 Accepted May 11, 1998 IE970736D