Hybrid Differential Evolution for Problems of Kinetic ... - ACS Publications

2876

Ind. Eng. Chem. Res. 2001, 40, 2876-2885

Hybrid Differential Evolution for Problems of Kinetic Parameter Estimation and Dynamic Optimization of an Ethanol Fermentation Process Feng-Sheng Wang,* Tzu-Liang Su, and Horng-Jhy Jang Department of Chemical Engineering, National Chung Cheng University, Chia-Yi 621, Taiwan, R.O.C.

Hybrid differential evolution (HDE) is applied to estimate the kinetic model parameters of batch fermentation for ethanol and glycerol production using Saccharomyces diastaticus LORRE 316. In this study, we considered two kinetic models for describing the dynamic behaviors of S. diastaticus LORRE 316. In the parameter estimation problem, we used the worst observed error for all experiments as an objective function so that the parameter estimation problem becomes the min-max estimation problem. Several numerical methods have been employed to solve the min-max estimation problem for comparison. HDE could use a small population size to obtain a more satisfied solution as compared from these computations. To validate the two kinetic models, we respectively used the two kinetic models to determine the optimal feed rates for the fedbatch optimization problem. Because the fedbatch optimization problem was a constrained dynamic optimization problem, we introduced the HDE with a multiplier updating method including adaptive penalty parameters to obtain the feasible feed rates. Such feed rates are then applied for fedbatch experiments in a 5 L fermenter for model validation. From the comparison of batch and fedbatch experiments, we observed that the proposed kinetic model was more adequate than Monod’s model. 1. Introduction Mathematical models that accurately predict physical phenomena are essential in many fields of engineering and science. The model provides the basis for design, control, and optimization of process systems. The objective of mathematical modeling is to obtain an expression that quantitatively describes the behavior of the process under consideration. The generality of a model depends on several factors, which include the complexity and information available concerning the process. These models frequently contain adjustable parameters that need to be determined from available experimental data. Estimation of kinetic parameters plays an important role for the development of mathematical modeling. In bioprocess engineering, various kinds of kinetic models have been proposed to quantitatively describe the dynamic behaviors of biological reaction systems. Several conventional methods, such as the graphical method and the gradient-based nonlinear optimization methods, have been employed as estimation techniques to obtain the kinetic parameters.1-3 The graphical method can be applied only to those problems that can be converted to linear regression problems. Basically, a gradient-based nonlinear optimization method does not have such a restriction, but it requires gradient information about the error function with respect to the parameter estimates and often gets stuck at local minima. Hence, a new estimation algorithm requires that good parameter values can be found without using gradient information; i.e., it depends minimally on the form of the kinetic model. Evolutionary algorithms (EAs) are a class of stochastic search and optimization methods that includes * Corresponding author. E-mail: [email protected]. Phone: +886 5 2428153. Fax: +886 5 2721206.

genetic algorithms (GAs), evolutionary programming, evolution strategies, genetic programming, and their variants.4 The EAs have been successfully applied to parameter estimation and a wide range of applications including design, scheduling, control, and others.5-8 Compared to traditional search and optimization methods, such as gradient-based and enumerative strategies, the EAs are robust, global, and generally more straightforward to apply in situations where there is little or no a priori knowledge about the problem to be solved. Differential evolution (DE) was one of the most excellent EAs for solving the real-valued test function suite of the first International Contest on Evolutionary Computation.9,10 The convergent speed of DE is very fast. In general, the faster convergence leads to a higher probability toward obtaining a local optimum. This fact results from diversity of population descending faster during solution progress. This drawback could be overcome by using a larger population size. By doing so, much computational time should be expended to evaluate the fitness function. This fact is particularly serious by using DE to solve optimal control problems due to expending much CPU time to solve differential equations.11 For saving CPU time, a modified DE has been applied to solve dynamic optimization of a continuous polymer reactor.12 Chiou and Wang13 have developed a hybrid differential evolution (HDE) to overcome such drawbacks. The HDE has been successfully applied to solve problems of fedbatch optimization and fuzzy decision making of fermentation processes.14,15 In this paper, we apply HDE to estimate kinetic parameters of batch fermentation for a Saccharomyces diastaticus LORRE 316, which is high ethanol tolerance yeast. The yeast LORRE 316 utilizes glucose to simultaneously produce ethanol and glycerol. The kinetic model should consist of the dynamic behavior of biomass, glucose, ethanol, and glycerol. The HDE is applied

10.1021/ie000544+ CCC: $20.00 © 2001 American Chemical Society Published on Web 06/02/2001

Ind. Eng. Chem. Res., Vol. 40, No. 13, 2001 2877

Figure 1. Configuration of the experimental setup.

to obtain the optimal parameters for the kinetic model. Several methods are also employed to solve the estimation problem for comparison. To validate the kinetic model, we apply the model to a dynamic optimization problem for fedbatch fermentation. We introduce a Lagrangian-like method for HDE to handle the dynamic optimization problem with inequality constraints. The optimal solution is then compared with the fedbatch experiment. 2. Experimental Material and Methods In this work, six batch experiments were used to estimate the kinetic model parameters. Six fedbatch experiments were used to validate the kinetic model. The experimental material and methods used in this work are expressed as follows. S. diastaticus LORRE 316 was supplied by Prof. George T. Tsao, Laboratory of Renewable Resources Engineering. Purdue University. This culture is a superproductive and ethanoltolerant yeast and has been found to be capable of producing ethanol from corn starch, yielding a final concentration as high as 17.5% (v/v). The nutrient medium (NM) contained 6 g/L yeast extract, 6 g/L peptone, 2 g/L K2HPO4, 4 g/L KH2PO4, 2 g/L MgSO4, and 2.5 g/L (NH4)2SO4. Glucose of 100, 150, and 200 g/L with NM were respectively used for batch fermentation. For each glucose concentration, we perform two experiments. A glucose concentration of about 10 g/L with NM was used for the initial medium in fedbatch fermentation. The glucose concentration in the feeding solution was 200 g/L in the fedbatch fermentation. The fermentation in the fedbatch experiments was started in a 1.5 L batch mode until the glucose concentration dropped to a desired value. A homemade fermenter with a total working volume of 5 L was used. The temperature of 35 °C was maintained by controlling the circulation of the cooling

water. In this work, the cell mass was determined offline spectrophotometrically at 600 nm and the corresponding dry weight of cells was obtained from an established calibration. However, a personal computer based on-line monitoring and controlling system, developed by Liu et al.16 for the fermentation processes, was used to on-line measure the concentrations of glucose, ethanol, and glycerol by a high-performance liquid chromatograph (HPLC). Figure 1 shows a schematic diagram of the experimental setup. An autosampling device (BenchMate II) from Zymark Corp. was connected to the fermenter and a HPLC (Waters 600) as shown in Figure 1. The on-line HPLC system for the analysis of glucose, ethanol, and glycerol in the fermentation broth was connected to the fermenter via the autosampling equipment, which could perform the pipetting, filtration, and dilution of the sample and the final injection onto the HPLC through automation based on a programmed procedure. The broth in the fermenter was automatically sampled by the autosampling device. The cells in the sampled solution were retained by a 0.45 µm syringe filter. The filtrate was diluted and automatically injected into the HPLC to analyze the concentrations of glucose, ethanol, and glycerol. The concentrations of glucose, ethanol, and glycerol were analyzed by the Bio-Rad Aminex HPX-87H column and a refractive index detector (Waters 410). The column temperature was kept at 85 °C, and the mobile phase was the deionized water, which was applied at a flow rate of 0.6 mL/min. The temperature in the detector was 40 °C. 3. HDE The HDE is an extended DE method. DE is a very simple population-based stochastic function minimizer, introduced by Storn and Price.9,10 This method has turned out to be one of the best GAs for solving the real-

2878

Ind. Eng. Chem. Res., Vol. 40, No. 13, 2001

Table 1. Basic DE and HDE Operations DE

HDE

1. representation and initialization 2. mutation operation 3. crossover operation 4. selection and evaluation 5. repeating of steps 2-4

1. representation and initialization 2. mutation operation 3. crossover operation 4. selection and evaluation 5. acceleration if necessary 6. migration naturally or enforced if necessary 7. repeating steps 2-6

valued test function suite of the first International Contest on Evolutionary Computation. The basic DE operations are shown in Table 1, which are similar to the conventional evolutionary algorithm. The DE computational algorithm is very simple and easy to implement, with only a few parameters required to be set by a user. The mutation in DE uses the difference between two chosen random individuals as a search direction. This is different from the conventional evolutionary algorithms. A mutant individual is yielded from a parent individual and the perturbed mutation. In DE, a fixed mutation factor multiplies the difference vector to yield the perturbed mutation. In HDE, a random mutation factor is applied to obtain more diversified individuals. In DE and HDE, the fitness of an offspring competes one to one with that of its parent. This competition, which is also different from the conventional evolutionary algorithms, gives rise to a faster convergence rate. However, this faster convergence also leads to a higher probability of obtaining a local optimum because the population diversity descends faster during the solution progress. This drawback can be overcome by using a larger population size, although much computation time is expended to evaluate the fitness function. This fact is particularly relevant when using DE to solve optimal control problems.11 The HDE can avoid using a larger population size and has been applied to solve bioreaction process optimization problems.13-15 The basic operations of the HDE are also shown in Table 1. When using an evolutionary algorithm to optimize a function, an acceptable tradeoff between convergence and diversity must generally be determined. Convergence implies a fast convergence although it may be to a local optimum. On the other hand, diversity guarantees a high probability of obtaining the global optimum. When the population diversity is small, the candidate individuals can be closely clustered. Therefore, the mutation and crossover operations can no longer generate the next better individual because a premature solution is obtained. In the HDE, an accelerated operation and a migration are embedded in the original DE algorithm, and these two operations act as a tradeoff operator. The accelerated operation is used to speed up convergence. 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 fitness after several generations. Under this situation, the accelerated operation can be used to speed up convergence. When the best fitness in 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 solution. The rate of convergence can be improved by the accelerated operation. However, faster descent usually

results in finding a local minimum. In addition, performing this operation so frequently can make the candidate individuals gradually cluster around the best individual so that the diversity of the population decreases. Furthermore, the closely clustered individuals cannot reproduce better individuals by mutation and crossover operations. As a result, a migration operation is performed to escape this local cluster because the new candidate individuals are regenerated on the basis of the best individual in the current generation. Correspondingly, the diversity of the candidates can still be retained using such a regeneration procedure. The migration operation in HDE is performed only if a measure of the population diversity fails to satisfy the desired tolerance. Chiou and Wang14 proposed a measure called the degree of population diversity to check whether the migration operation should be performed. The degree of population diversity is between 0 and 1. A zero value implies that all of the genes are clustered around the best individual. On the other hand, the value of one indicates that the current candidate individuals are a completely diversified population. The desired tolerance for population diversity is accordingly assigned within this region. A zero tolerance implies that the migrating operation in HDE is switched off. A tolerance of 1 implies that the migration operation is performed in every generation. 4. Batch Model 4.1. Kinetic Formulation. The basic equations that describe the growth of microorganisms, the consumption of glucose, and the formation of products are given as follows:

dx ) rx dt

(1)

ds ) rs dt

(2)

dp1 ) rp 1 dt

(3)

dp2 ) rp 2 dt

(4)

where x is the concentration of cell mass, s is the concentration of glucose, p1 is the concentration of ethanol, and p2 is the concentration of glycerol. In this study, the reaction rates for cell, glucose, ethanol, and glycerol are expressed as follows:

rx ) µx rs ) -

qp1 Yp1/s

x-

(5) qp2 Yp2/s

x

(6)

rp1 ) qp1x

(7)

rp2 ) qp2x

(8)

where Yp1/s is the ethanol yield factor, Yp2/s is the glycerol yield factor, and the unstructured kinetic model for the


specific cell growth and product formation are respectively expressed in the forms

min J ) min max {Jk, k ) 1, ..., Nexp}

µ)

u

Kp1

µms

Kp2 2

Ks + s + s /KsI Kp1 + p1 + p1 /Kp1I Kp2 + p2 + p22/Kp2I 2

(9) νp1s

q p1 )

K′p1

(10)

K′s + s + s2/K′sI K′p1 + p1 + p12/K′p1I νp2s

q p2 )

K′p2

(11)

2

K′′s + s + s /K′′sI K′p2 + p2 + p22/K′p2I

There are 19 parameters that have to be estimated from experimental observations. Ethanol/glycerol fermentation is a typical example of a process in which utilization of the carbon source cannot be split into two processes, namely, production of biomass and formation of product. Ethanol/glycerol is the end product of energy metabolism associated directly with biomass growth. It is, however, an inhibitor capable of retarding metabolic activity. Most of the unstructured models are empirical and based on either Monod’s equation or its numerous modifications which take into account the inhibition of microbial growth by a high concentration of product and/ or substrate. In this study, we consider the glucose, ethanol, and glycerol inhibition on biomass growth as expressed in (9). However, the formation of ethanol and glycerol not only associates with microbial growth, as expressed in (7) and (8), but also inhibits by glucose and itself, as expressed in (10) and (11). The mathematical estimation of model parameters is based on minimization of some quantity that can be calculated and that is a function of the parameters to be estimated. If the model under consideration is linear, the estimation is generally an easy task. There exists, however, no general theory for nonlinear parameter estimation. The least-squares error is commonly employed as a criterion to inspect how close the computed profiles of the state variables come to the experimental observations. In this study, we consider the timeweighted least-squares error as a criterion for each experiment. The criterion is expressed in the form

Jk )

becomes a min-max problem.19 We therefore have

1

Ns

{

∑ tj

Nsj)1

[xe(tj) - x(tj)]2

+

[se(tj) - s(tj)]2

xemax2 [p1e(tj) - p1(tj)]2 p1emax2

semax2 +

+

}

[p2e(tj) - p2(tj)]2 p2emax2

(12)

where xe(tj), se(tj), p1e(tj). and p2e(tj) are the measured data at t ) tj; x(tj), s(tj), p1(tj), and p2(tj) are the concentrations calculated using the model, and xemax, semax, p1emax, and p2emax are the maximum measured concentrations. Here Ns is the number of the sampling data. The least-squares regression sums up every observed error in (12) to yield an objective function.17-19 In this study, we consider the worst observed error for all experiments as an objective function. This approach is a special case of multiobjective parameter estimation problems so that the parameter estimation problem

u

k

(13)

where Nexp is the number of experiments and u is a vector of the estimated parameters. Therefore, we can apply the HDE to solve the min-max estimation problem to obtain the kinetic parameters. 4.2. Estimated Results. All computations were performed on a Pentium II computer using Microsoft Windows NT. We use Compaq Visual Fortran to implement the HDE algorithm. The user has to provide four setting factors for the HDE. The setting factors used for all runs are listed as follows. The crossover factor is set to 0.5. Two tolerances used in the migration operation are set to 0.05. The population size of 5 is used in the computation. In HDE, the mutation factor is taken as a random number in [0, 1]. The HDE with various generations is used to solve the min-max estimation problem. The best objective function value of 3.495 235 × 10-1 as shown in Table 2 is obtained using the total function call of 334 247, which is equivalent to the CPU time of about 49 min on a Pentium II 400 computer. This computation uses a migration operation of 128 and an acceleration operation of 458. Table 2 also shows the solution progressing for the optimal objective function value and its corresponding total function call. The objective function value decreases progressively as the generation increases. Figures 2-4 show the optimal computed profiles of cell mass, glucose, ethanol, and glycerol for various initial concentrations of glucose. The two experimental data (data points) used for estimation are also shown in these figures for comparison. These computed profiles fit the experimental data satisfactorily. To compare the performance of the HDE, we use a direct search method, a subroutine BCPOL in IMSL Math/Library,20 to solve the min-max estimation problem of parameter estimation. Because the solution obtained by the BCPOL subroutine is strongly dependent on the provided initial guess and stopping tolerance, we therefore run the problem 10 times with randomly assigned initial guesses and a stopping tolerance of 10-7. The best objective function value of 10 runs is 7.680 123 × 10-1 using the total function call of 15 009. However, this value can be reached by using the HDE with the total function call of 2169. The total function call is equivalent to CPU time so that this situation indicates that the convergence rate of the HDE is quicker than that of the BCPOL solver. In this work, we also provide the solution shown in the second column of Table 2, which is obtained by HDE, as the initial guess for the BCPOL solver to find the more accurate solution. The refined solution and total function call used for BCPOL at various generations are shown in the third and fourth columns of Table 2 for comparison. Although BCPOL can refine a solution, it is a premature solution and unable to obtain a global solution. The original DE was also applied to solve this minmax estimation problem for comparison. In the original DE, the mutation factor is fixed and provided by a user. We tried various pairs of mutation and crossover factors for the original DE to solve the problem. However, we always obtained a premature solution using the fixed mutation factor. We then used the random mutation factor at every generation. The other setting factors are identical with those used in HDE except with the

2880


Table 2. Optimal Objective Function Value J(u*) Obtained by Various Approaches and Number of Total Function Calls (TFC) Used at Various Generationsa HDE(Np)5) generation 500 1000 2000 4000 8000 16000 32000 64000

J(u*)

HDE(Np)5) + BCPOL TFC

6.450 443 × 10-1 4 663 8 431 4.552 529 × 10-1 4.224 163 × 10-1 14 455 3.879 576 × 10-1 26 845 3.603 710 × 10-1 49 115 3.529 581 × 10-1 90 625 3.497 682 × 10-1 172 687 3.495 235 × 10-1 334 247

J(u*)

DE(Np)20)

TFC

J(u*)

4.628 923 × 10-1 4663 + 12 011 3.868 323 × 10-1 8431 + 12 695 4.151 873 × 10-1 14 455 + 9877 3.854 823 × 10-1 26 845 + 3178 3.587 645 × 10-1 49 115 + 9231 3.521 833 × 10-1 90 625 + 7170 3.497 508 × 10-1 172 687 + 2921 3.495 174 × 10-1 334 247 + 3031

simple GA(Np)50) TFC

4.823 250 × 10-1 10 020 4.038 545 × 10-1 20 020 3.633 052 × 10-1 40 020 3.560 734 × 10-1 80 020 3.524 761 × 10-1 160 020 3.518 701 × 10-1 320 020 3.511 426 × 10-1 640 020 3.506 484 × 10-1 1 280 020

J(u*)

TFC

5.443 097 × 10-1 25 050 5.116 637 × 10-1 50 050 4.970 216 × 10-1 100 050 3.888 264 × 10-1 200 050 3.865 303 × 10-1 400 050 3.732 605 × 10-1 800 050 3.673 204 × 10-1 1 600 050 3.620 793 × 10-1 3 200 050

a Optimal parameters obtained by HDE: µ ) 0.6397, K ) 4.895, K ) 90.35, K ) 10.0, ν ) 3.429, K′ ) 17.45, K′ ) 460.4, ν ) m s p1 p2 p1 p2 s p1 1.748, K′′s ) 499.4, K′p2 ) 32.26, KsI ) 1115.1, Kp1I ) 12.89, Kp2I ) 12.45, K′sI ) 110.3, K′p1I ) 3.71, K′′sI ) 27.78, K′p21I ) 0.3797, Yp1/s ) 0.505, and Yp2/s ) 0.1955.

Figure 2. Concentration of cell mass (x), glucose (s), ethanol (p1), and glycerol (p2) for the batch fermentation with an initial glucose of 100 g/L. Data points are for the experiments. Solid curves are for the optimal estimated results.


migrating and accelerated operations switched off. We still obtained the premature objective function value of 0.5855 using 320 005 total function calls. A satisfied solution was unable to be obtained by DE using the 5 population size because of the vanishing of the mutation

and crossover operations. A satisfied solution can be obtained by using the population size of 20. Table 2 shows the optimal objective function value and the used total function call at various generations. From this table, we observe that the original DE uses the total function call of 1 280 020 to obtain the optimal objective function value of 3.506 484 × 10-1. However, this solution can be reached by using the HDE with the total function call of 128 440. The DE requires use of about 10 times the function call to obtain a solution identical with that observed from this comparison. To compare the convergent rate, a Fortran GA was also applied to solve the min-max estimation problem. The Fortran GA solver is accessed from a public web site, http://www.staff.uiuc.edu/∼carroll/ga.html.21 The GA solver recommended some setting factors for the user. We use the recommended factors with the population size of 50 to solve the min-max estimation problem. Table 2 shows the progress of the best objective function value for the GA. At the generations of 64 000, the GA yields a premature solution of 3.620 793 × 10-1 by using the total function call of 3 200 050. This solution can be reached by using the HDE with the total function call of 44 405. The simple GA requires use of 72 times the function call to obtain a solution identical with that from this comparison. The kinetic model in (9)-(11) is not unique. Most of the unstructured models are empirical and based on Monod’s equation. In this study, we also used Monod’s model to fit the experimental results for comparison. The reaction rates for cell, glucose, ethanol, and glycerol are expressed as follows: rx ) k1xs/(k2 + s), rs ) -k3rx, rp1 ) k4rx, and rp2 ) k5rx. We then applied the HDE using the same setting factors as those in the previous run to obtain the optimal kinetic parameters as k1 ) 0.3857, k2 ) 0.1, k3 ) 12.44, k4 ) 4.9, and k5 ) 0.3852. The best objective function value of 8.034 949 × 10-1 was obtained using the total function call of 327 111. Monod’s model is less adequate than the proposed model in (5)(8) as compared with the objective function value.

5. Fedbatch Optimization 5.1. Fedbatch Model. An exact kinetic model provides the basis for design, control, and optimization of bioprocess systems. In this work, we use the estimated model to design an optimal feed rate for fedbatch fermentation. The dynamic mass balance equations for the fedbatch process are expressed as follows:


dx F ) rx - x dt V

(14)

ds F ) rs + (sF - s) dt V

(15)

F(t) ) F(j), j ) 1, ..., Nu

dp1 F ) rp1 - p1 dt V

(16)

dp2 F ) rp2 - p2 dt V

(17)

where Nu is the number of time partitions. The infinite dimensional problem (14)-(26) is then approximated as a parameter selection problem. To avoid a wide variation of feed control action, each action is bounded by

dV )F dt

(18)

where V is the working volume of the fermentor and F is the feed rate. The kinetic models for the specific cell growth and product formation are expressed in (5)-(8). Some physical constraints should be considered in the optimization problem. The total volume of the fermentor is bounded by

g1 ) V(t) - Vf e 0

g3 )

p1(t) V(t) - p1(0) V(0) [V(t) - V(0)]sF + V(0) s(0) - V(t) s(t) p2(t) V(t) - p2(0) V(0) [V(t) - V(0)]sF + V(0) s(0) - V(t) s(t)

-Y ˆ p1/s e 0 (20) -Y ˆ p2/s e 0 (21)

From a stoichiometric point of view, both theoretical ˆ p2/s, for the ethanol formation yield factors, Y ˆ p1/s and Y from glucose are 0.51. If the constraints (20) and (21) are not included in the optimization problem, the unrealistic predicted values of p1(tf) may be found.22 To reduce the separation cost, the residue glucose is restricted by

g4 ) s(t) - sr e 0

(22)

where sr is the desired residue glucose. The goal of this optimization problem is to determine the optimal feed rate F(t) and fermentation time tf to maximize the ethanol production rate. The objective function is therefore expressed as

max J ) [p1(tf) V(tf) - p1(0) V(0)]/tf F(t), tf

(23)

Here, the decision variables are bounded in the solution space as follows:

0 e F(t) e Fmax

(24)

tfmin e tf e tfmax

(25)

The feed rate F(t) in this dynamic optimization problem is in terms of time so that such a problem is an infinite

(26)

g4+l ) |F(l + 1) - F(l)| - 0.25Fmax e 0, l ) 1, ..., Nu - 1 (27) Penalty function methods are commonly applied to transform the constrained optimization problem into a unconstrained problem. We therefore have 4

Ja(u) ) J +

∫0 Rk(t) 〈gk(z(t),u)〉+2 dt + ∑ k)1 tf

Nu-1

(19)

In addition, the stoichiometry of the ethanol and glycerol formation from glucose for a cell must be obeyed in the fermentation process, so two constraints for each specific yield factor have been introduced to obtain a realistic solution in this optimization problem. According to the definition of the yield factor, these constraints are

g2 )

dimensional problem. To solve this problem efficiently, the feed rate is first represented by a finite set of control actions F(j) in the time interval tj-1 e t < tj by

R4+k〈g4+k〉+2 ∑ k)1

(28)

where the decision vector u consists of the control actions F(j) and the fermentation time tf and the state variable z(t) is expressed as z ) [x, s, p1, p2, V]T. Here, the penalty parameters Rk are positive constant, and the bracket operator in (28) is defined as 〈gk〉+ ) max {gk, 0}. Penalty terms associated with the constraints are added to the objective function. In such a way, penalty terms reflect the violation of the constraints and assign the high cost of the penalty function to candidate points far from the feasible region. While we use HDE to solve the penalty problem, any candidate individuals that violate the constraints would inherit a worse fitness and be hard to survive. The penalty function methods are easy to implement and are considered efficient. However, the main limitation of the penalty functions is the degree to which each constraint is penalized. Powell23 has noted that the classical optimization methods including penalty functions have certain weaknesses that become most serious when penalty parameters are large. More seriously, large penalty parameters make the penalty function illconditioned so that it is difficult to achieve a good solution. On the other hand, if the penalty parameters are too small, the constraint violation does not contribute a high cost to the penalty function. Thus, the optimal solution based on the penalty function may not be feasible. Therefore, how to choose appropriate penalty parameters is not trivial. Lagrange methods are a classical approach for solving real-valued constrained optimization problems. Such methods can significantly improve the drawbacks of penalty methods. In this paper, we introduce the HDE with multiplier updating to solve the dynamic optimization problem. The augmented Lagrange function for the dynamic optimization is defined as 4

Ja(u) ) J +

∑ ∫0 Rk(t) [〈gk(z(t),u) + σk〉+2 - σk2] dt + k)1 tf

Nu-1

∑ R4+k[〈g4+k + σk〉+2 - σk2]

k)1

(29)

2882


Table 3. HDE with a Multiplier Updating Method Including Adaptive Penalty Parameters Step 1. Set the initial iteration l ) 0, the initial multipliers σlk ) σ0k ) 0, and the initial penalty parameters Rk > 0. Set the tolerance of the maximum constraint violation, K (e.g., K ) 1032), and the scalar factors, ω1 > 1 and ω2 > 1. Step 2. Use HDE to minimize the augmented Lagrange function Ja(u,Rk,σlk). Let ulb be a minimum solution to the problem Ja(u,Rk,σlk). Step 3. Evaluate the maximum constraint violation as ˆ k ) maxk |max (gk, -σk)| and establish the following sets of inequality constraints whose violation does not improve by the factor ω1: II ) {k:||max (gk, -σk)| > k/ω1with k ) 1, ..., mi}. ) σlk/ω2 for all k ∈ II and go to step 7. Otherwise, go to step 5. Step 4. If ˆ K g K, let Rk ) ω2Rk and σl+1 k ) 〈gk(ulb) + σlk〉+. Step 5. Update the multipliers as follows: σl+1 k ) σl+1 Step 6. If ˆ K e K/ω1, let K ) ˆ K and go to step 7. Otherwise, let Rk ) ω2Rk and σl+1 k k /ω2 for all k ∈ II. Let K ) ˆ K and go to step 7. Step 7. If the maximum iteration is reached, stop. Otherwise, repeat steps 2-6.

where the corresponding multipliers are updated at each generation as

σG k (t)

)

σG-1 (t) k

+ max

{gG-1 (z(t),u), k

-σG-1 (t)}, k

k ) 1, ..., 4 (30)

and G G-1 G-1 G-1 ) σ4+k + max {g4+k , -σ4+k }, k ) 1, ..., Nu σ4+k (31)

where G is the generation index in the HDE algorithm. The penalty parameters in (29) are, in general, fixed for the evolutionary computation. However, when we use smaller penalty parameters for the augmented Lagrange function, an infeasible solution may still be obtained. To overcome such a drawback, we introduce a global convergence method of Powell23 into the HDE with a multiplier updating method to enforce global convergence for the constrained optimization problem. Powell’s method applied an adaptive penalty parameter strategy in a gradient-based optimization method to solve nonlinear programming problems with general constraints. In this study, we modify Powell’s algorithm to suit the evolutionary algorithm of HDE. Table 3 shows the algorithm for HDE with a multiplier updating method including adaptive penalty parameters (HDEMUM-APP). Steps 3, 4, and 6 in Table 3 are used to improve the constraint violation and update the penalty parameters. In step 4, if the constraint violations do not improve, i.e., ˆ K g K, we increase the penalty parameters by the factor ω2 (e.g., ω2 ) 10 in this study) and reduce the corresponding multipliers by the same factor, thus keeping the multipliers unchanged. In step 6, we use the factor ω1 (e.g., ω1 ) 4 in this study) to check whether the constraint violation has improved by the factor ω1. The penalty parameters and the corresponding multipliers are updated in this step only when the constraint violations do not improve by the factor ω1. 5.2. Computational Results. To solve this problem, the feed rate F(t) is first approximated by 10 time partitions. Following procedures similar to those discussed in the previous section, we respectively use HDE, DE, BCPOL, and GA to solve the dynamic optimization problem. We also used the same setting factors as those discussed in the previous section for HDE. Here, the maximum iterations of 5000 are used for HDE to obtain a minimum solution for the augmented Lagrange function as shown in step 2 in Table 3. The initial, fed, and residue concentrations of glucose are set as 10, 200, and 0.1 g/L in the computation. The maximum feed rate and total working volume are considered as Fmax ) 1.0 L/h and 5 L. We run HDEMUM-APP three times. Each run uses different initial penalty parameters but the same maximum iteration

Table 4. Maximum Ethanol Production Rate for Various Methods method

item

HDE-MUM-APP(Np)5) Rk ) 1 Rk ) 10 Rk ) 106 HDE-MUM-FPP(Np)5) Rk ) 1 Rk ) 10 Rk ) 106 HDE-PFM(Np)5) Rk ) 1 Rk ) 10 Rk ) 106 DE-PFM(Np)20) Rk ) 1 Rk ) 10 Rk ) 106 BCPOL-PFM Rk ) 1 Rk ) 10 Rk ) 106 simple GA-PFM(Np)50) Rk ) 1 Rk ) 10 Rk ) 106 a

max J (g/h)

SCVa

TFC

23.542 23.542 23.489 24.541 23.659 23.509 27.266 25.084 23.357 27.267 24.685 23.542 27.263 25.082 19.502 26.230 25.334 23.651

1.0617 × 10-8 3.8619 × 10-9 0.0 1.3156 × 10-1 1.0602 × 10-2 6.4923 × 10-6 3.1630 5.0366 × 10-1 0.0 3.1677 4.5598 × 10-1 8.8574 × 10-6 3.1636 5.0276 × 10-1 2.3587 × 10-5 3.3143 1.6961 5.5445 × 10-2

324 218 316 912 310 128 307 488 317 641 319 802 259 315 262 068 315 027 1 000 020 1 000 020 1 000 020 6193 8307 12 908 2 500 050 2 500 050 2 500 050

4 Nu-1 SCV ) ∑k)1 ∫t0f〈gk(z(t),u〉+ dt + ∑k)1 〈g4+k〉+.

of 10 in the outer loop for multiplier updating. In the computation, we define the sum of constraint violations (SCV) to inspect the constraint feasibility at the final solution. The SCV is zero for the initial penalty parameters of 106, as shown in Table 5. In this run, the penalty parameters are updated from 106 to 109. The SCV of zero indicates that the minimum solution is completely feasible. A maximum production rate of 23.489 g/h is obtained for this run. Such a computation uses an acceleration operation of 3164, a migrating operation of 41, and a total function call (TFC) of 319 802, which is equivalent to a CPU time of about 18 min on a Pentium II 400 computer. The optimal feed rate (F1*) and concentration of ethanol (p1) and glucose (s) are shown in Figures 5 and 6, respectively. We then run the problem with the penalty parameters of 1 and 10, respectively. For both runs, the penalty parameters are updated from their initial value to 104 at the final solution. Table 4 shows the maximum production rate, SCV, and TFC for each run. The maximum production rates for the three runs are nearly identical although the penalty parameters used for each run show a wide difference. We then apply the HDE-MUM including fixed penalty parameters (HDE-MUM-FPP) to solve the problem. The maximum production rate, SCV, and TFC for various penalty parameters are shown in Table 4. We are unable to obtain a complete feasible solution by HDE-MUM-FPP, as observed from Table 4. However, a nearly feasible solution of 23.509 g/h is obtained for the penalty parameters of 106. The HDE with the penalty function method (HDEPFM) is also applied to solve the problem. The maximum production rate, SCV, and TFC for various penalty parameters are also shown in Table 4 for comparison. The production rate obtained by HDE-PFM with

Ind. Eng. Chem. Res., Vol. 40, No. 13, 2001 2883 Table 5. Comparison of Fedbatch Experiments and Computational Results for Various Feed Rate and Feed Glucose Concentrationa F1* and sF ) 199 g/L

Fc and sF ) 200 g/L

F2* and sF ) 203 g/L

F2* and sF ) 217 g/L

expt

computation

expt

computation

expt

computation

expt

computation

item

mean

M0

M1

mean

M0

M1

mean

M0

M1

mean

M0

M1

production rate (g/h) final ethanol (g/L) residue glucose (g/L) variance, eq 14

24.02 63.43 0.0

23.49 61.56 0.1 0.3149

21.56 57.42 0.0 0.390

20.38 53.58 4.76

20.35 55.76 16.63 0.3991

21.68 56.33 0.02 0.7586

24.76 65.7 0.0

23.29 61.88 0.23 0.1702

23.26 57.93 0.0 0.2775

24.61 65.45 0.0

24.79 65.92 0.59 0.1471

23.30 62.06 0.0 0.4099

a F * indicates that the optimal feed rate obtained using the kinetic model (5)-(8), F * is for Monod’s model, and F is for the constant 1 2 c feed rate. M0 is for the kinetic model (5)-(8) and M1 for Monod’s model.


Figure 5. Optimal feed rates, F1* and F2*, and constant feed rate Fc used for the fedbatch fermentation. F1* is for the optimal feed rate obtained using the kinetic model (5)-(8), F2* is for Monod’s model, and Fc is for the constant feed rate.

penalty parameters of 1 and 10 is larger than that obtained by HDE-MUM-APP. However, the maximum production rate obtained by HDE-PFM with Rk ) 1 and 10 is infeasible. While the penalty parameter is 106, the penalty function method is able to obtain a completely feasible solution but smaller than that of HDE-MUMAPP. The original DE with the penalty function method (DE-PFM) was also applied to solve this dynamic optimization problem for comparison. In the original DE, the mutation factor is fixed and provided by a user. We tried various pairs of mutation and crossover factors

Figure 6. Concentration of glucose and ethanol for the fedbatch fermentation using the optimal feed rate F1* and constant feed rate Fc, respectively. Data points are for the experiments. Solid curves are for the optimal fedbatch fermentation. Dashed curves are for the constant fedbatch fermentation.

for the original DE to solve the problem. However, we always obtained a premature solution using the fixed mutation factor. We then used the random mutation factor in every generation. The other setting factors are identical with those used in HDE except with the migrating and acceleration operations switched off. We still obtained the premature infeasible solution of 17.49 g/h using 250 005 total function calls. A satisfied solution was unable to be obtained by DE using the five population size because of the vanishing of the mutation and crossover operations. A reasonable solution can be obtained by using the population size of 20. Table 4 shows the best solution for the three runs and the total function call used. From this table, we observe that the DE-PFM with Rk ) 1 and 10 is unable to obtain a feasible solution. However, the solution obtained by DE-PFM with Rk ) 106 is nearly feasible. We use BCPOL with the penalty function method (BCPOL-PFM) to solve the dynamic optimization problem. Because the solution obtained by the BCPOL subroutine is strongly dependent on the provided initial guess and stopping tolerance, we therefore run the problem 10 times with randomly assigned initial guesses and a stopping tolerance of 10-7. A maximum production rate of 19.502 g/h is obtained by BCPOL-PFM with penalty parameters of 106. The total function call of 12 908 is used in this run. However, the solution is smaller than that of the above-mentioned methods. From Table 4, we also observe that the BCPOL-PFM with Rk ) 1 and 10 is unable to obtain a feasible solution.

2884


Finally, we use the simple GA with a penalty function method (GA-PFM) to solve the dynamic optimization problem for comparison. A population size of 50 is used to solve the dynamic optimization problem. The solutions for various penalty parameters are expressed in Table 4. The situation is similar to that of the abovementioned methods with a penalty function method. We have to use larger penalty parameters for the penalty function method to yield a feasible solution. The optimal feed rate is in terms of an employed kinetic model because the fedbatch optimization problem depends on the kinetic model as observed from (14) to (17). In this study, we also use Monod’s model for the mass balance equations (14)-(17) to obtain the optimal feed rate F2* as shown in Figure 5. This optimal feed rate will be applied in the next subsection for fedbatch experiments to validate the process model. 5.3. Experimental Results. The dynamic optimization problem discussed in this work was validated by experiments. The materials and methods used for fedbatch experiments were discussed in section 2. In the experiments, we use various feed rates and feed glucose concentrations for fedbatch fermentation. Figure 5 shows two optimal feed rates and a constant feed rate used for the fedbatch experiments. The optimal feed rate F1* is obtained using the kinetic model (5)-(8) as discussed in the above-mentioned section. However, the optimal feed rate F2* is obtained using Monod’s model. A fedbatch experiment was started in a 1.5 L batch mode until the glucose concentration dropped to the desired value, which was the initial condition used in the fedbatch optimization problem. Figure 6 shows two experimental results (circle data points) using the optimal feed rate F1* and a feed glucose concentration sF of 199 g/L. Table 5 shows the mean value of the ethanol production rate, final ethanol concentration, and residue glucose for the experiments. These optimal fedbatch experiments yield an ethanol production rate of 24.02 g/h. We then use the feed rate F1* and sF ) 199 g/L to compute the concentration profiles using the kinetic model (5)-(8) and Monod’s model, respectively. The computational results are shown in Table 5. To compare the applicability of the two models, the variance in (13) was computed and shown in Table 5. The proposed model (5)-(8) is more adequate than Monod’s model as compared with that from the variance. Figure 6 shows the computational glucose and ethanol profiles (solid curves) for the proposed model. The experimental data can fit the computed results satisfactorily as observed from this figure. We proceed with two fedbatch experiments with a constant feed rate of 0.27 g/L. The experimental data and computed results using two different models are shown in Table 5. This constant fedbatch experiment yields the ethanol production rate of 20.38 g/h. Under the constant feed rate, S. diastaticus LORRE 316 is unable to completely consume glucose to produce ethanol during the fermentation time of 13 h. As a result, glucose at the final time remains. For the computational results, the residue glucose using the proposed kinetic model is overpredicted, as is observed from this table. However, the residue glucose using Monod’s model is underpredicted. The variance of 0.3991 for the proposed model is smaller than that for Monod’s model, as shown in Table 5. Figure 6 also shows two experimental results (triangle data points) and computed ethanol and glucose profiles using the proposed model.

Figure 7. Concentration of cell mass (x), glucose (s), ethanol (p1), and glycerol (p2) for the fedbatch fermentation using the optimal feed rate F2* and the feed glucose of 203 g/L. Data points are for the experiment. Solid curves are for the computed profiles using the kinetic model (5)-(8). Dashed curves are for the computed profiles using Monod’s model.

Figure 8. Concentration of cell mass (x), glucose (s), ethanol (p1), and glycerol (p2) for the fedbatch fermentation using the optimal feed rate F2* and the feed glucose of 203 g/L. Data points are for the experiment. Solid curves are for the computed profiles using the kinetic model (5)-(8). Dashed curves are for the computed profiles using Monod’s model.

To further validate the kinetic model, we also use Monod’s model for the mass balance equations (14)(17) to obtain the optimal feed rate F2* as shown in Figure 5. This optimal feed rate is then applied for two fedbatch experiments with various feed glucose concentrations. Table 5 summaries the production rate, final ethanol concentration, residue glucose, and variance for the two different models. Figures 7 and 8 show the experimental data and computational profiles for various feed glucose concentrations. The proposed model is more adequate than Monod’s model as compared from Table 5 and these figures. 6. Conclusions We have applied the HDE to estimate the kinetic parameters of S. diastaticus LORRE 316 and to determine the optimal feed rate for the fedbatch fermentation. Several methods have also been employed to solve both problems. From the comparison, the HDE can use a small population size to obtain a more satisfied solution. To validate the kinetic model, we have carried


out the fedbatch experiments with two optimal feed rates and a constant feed rate. From the compared experiments, the proposed kinetic model is more adequate than Monod’s model. Acknowledgment Financial support provided by the National Science Council of Republic of China under Grant NSC89-2214E194-002 is gratefully acknowledged. Nomenclature F(t) ) feed rate (L/h) Ks ) saturation coefficient for cell growth on glucose K′s ) saturation coefficient for ethanol production on glucose K′′s ) saturation coefficient for glycerol production on glucose KsI ) inhibition coefficient for cell growth on glucose K′sI ) inhibition coefficient for ethanol production on glucose K′′sI ) inhibition coefficient for glycerol production on glucose Kp1 ) saturation coefficient for cell growth on ethanol Kp2 ) saturation coefficient for cell growth on glycerol K′p1 ) saturation coefficient for ethanol production on ethanol K′p2 ) saturation coefficient for ethanol production on glycerol Kp1I ) inhibition coefficient for cell growth on ethanol Kp2I ) inhibition coefficient for cell growth on glycerol K′p1I ) inhibition coefficient for ethanol production on ethanol K′p2I ) inhibition coefficient for glycerol production on glycerol qp1 ) specific ethanol production rate qp2 ) specific glycerol production rate p1(t) ) concentration of ethanol (g/L) p2(t) ) concentration of glycerol (g/L) p1e(ti) ) measured concentration of ethanol at the sampling time ti p2e(ti) ) measured concentration of glycerol at the sampling time ti s(t) ) concentration of glucose (g/L) se(ti) ) measured concentration of glucose at the sampling time ti sF ) feed glucose concentration (g/L) tf ) fermentation time (h) V(t) ) working volume (L) Vf ) final working volume (L) x(t) ) concentration of cell mass (g/L) xe(ti) ) measured concentration of cell mass at the sampling time ti Yp1/s ) yield coefficient for ethanol from glucose Yp2/s ) yield coefficient for glycerol from glucose Greek Letters µ ) specific cell growth rate µm ) coefficient of maximum specific growth rate νp1 ) coefficient of maximum specific ethanol production rate νp2 ) coefficient of maximum specific glycerol production rate

Literature Cited (1) Baltes, M.; Schneider, R.; Sturm, C.; Reuss, M. Optimal Experimental Design for Parameter Estimation in Unstructured Growth Models. Biotechnol. Prog. 1994, 10, 480. (2) Cazzador, L.; Lubenova, V. Nonlinear Estimation of Specific Growth Rate for Aerobic Fermentation Processes. Biotechnol. Bioeng. 1995, 47, 626. (3) Thornkill, N. F.; Manela, M.; Campbell, J. A.; Stone, K. M. Two Methods of Selecting Smoothing Splines Applied to Fermentation Process Data. AIChE J. 1994, 40, 716. (4) Back, T.; Fogel, D.; Michalewicz, Z. Handbook of Evolutionary Computation; Oxford University Press: New York, 1997. (5) Androulakis, I. P.; Venkatasubramanian, V. A Genetic Algorithmic Framework for Process Design and Optimization. Comput. Chem. Eng. 1991, 15, 217. (6) McKay, B.; Willis, M.; Barton, G. Steady-State Modelling of Chemical Process Systems using Genetic Programming. Comput. Chem. Eng. 1997, 21, 981. (7) Edwards, K.; Edgar, T. F.; Maousiouthakis, V. I. Kinetic Model Reduction Using Genetic Algorithms. Comput. Chem. Eng. 1998, 22, 239. (8) Hugget, A.; Sebastian, P.; Nadeau, J. P. Global Optimization of A Dryer by Using Neural Networks and Genetic Algorithms. AIChE J. 1999, 45, 1227. (9) 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; IEEE: Piscataway, NJ, 1996; p 842. (10) Storn, R.; Price, K. V. Differential Evolution: A Simple and Efficient Heuristic for Global Optimization Over Continuous Spaces. J. Global Opt. 1997, 11, 341. (11) 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. (12) Lee, M. H.; Han, C.; Chang, K. S. Dynamic Optimization of a Continuous Polymer Reactor Using a Modified Differential Evolution Algorithm. Ind. Eng. Chem. Res. 1999, 38, 4825. (13) Chiou, J. P.; Wang, F. S. A Hybrid Method of Differential Evolution with Application to Optimal Control Problems of a Bioprocess System. IEEE Conference on Evolutionary Computation, Anchorage, AK, 1998; IEEE: Piscataway, NJ, 1998; p 627. (14) Chiou, J. P.; Wang, F. S. Hybrid Method of Evolution Algorithms for Static and Dynamic Optimization Problems with Application to a Fedbatch Fermentation Process. Comput. Chem. Eng. 1999, 23, 1277. (15) Wang, F. S.; Jing, C. H.; Tsao, G. T. Fuzzy Decision Making Problems of Fuel Ethanol Production Using A Genetically Engineered Yeast. Ind. Eng. Chem. Res. 1998, 37, 3434. (16) Liu, Y. C.; Wang, F. S.; Lee, W. C. On-line Monitoring and Controlling System for Fermentation Processes. Biochem. Eng. J. 2001, 7, 17. (17) Kim, I. W.; Liebman, M. J.; Edgar, T. F. Robust Error-inVariables Estimation Using Nonlinear Programming Techniques. AIChE J. 1990, 36, 985. (18) Ong, S. L. Least-Square Estimation of Batch Culture Kinetic Parameters. Biotechnol. Bioeng. 1983, 25, 2347. (19) Wang, F. S.; Sheu, J. W. Multiobjective Parameter Estimation Problems of Fermentation Processes Using A High Ethanol Tolerance Yeast. Chem. Eng. Sci. 2000, 55, 3685. (20) IMSL Math/Library User’s Manual; IMSL, Inc.: Houston, 1991. (21) Carroll, D. L. Chemical Laser Modeling with Genetic Algorithms. AIAA J. 1996, 34, 338. (22) Fu, C.; Barford, J. P. Non-Singular Optimal Control for Fed-Batch Fermentation Processes with A Differential-Algebraic System Model. J. Process. Control 1993, 3, 211. (23) Powell, M. J. D. Algorithms for Nonlinear Constraints that Use Lagrangian Functions. Math. Program. 1978, 14, 224.

Received for review June 1, 2000 Accepted April 4, 2001 IE000544+

Hybrid Differential Evolution for Problems of Kinetic ... - ACS Publications

Recommend Documents