Optimization of A Fed-Batch Fermenter Producing Baker’s Yeast Using Simulated Annealing Method

INTRODUCTION

Biochemical processes produce high value end products like vitamins and antibiotics and bulk products such as ethanol and baker’s yeast. The traditional and large scale beer and ethanol manufacturing processes based on using Saccharomyces cerevisiae (also called Baker’s yeast). Novel applications of S. cerevisiae have been developed during recent years. An example of its new application is insulin production for diabetics. S. cerevisiae like most other yeasts multiply by budding. Its production is carried out in a fed-batch mode in which a feed stream containing substrate and/or nutrients is fed into the fermenter during operation to gain higher performance levels, as compared to the batch mode. Numerous of researches have been carried out by previous workers to develop suitable, fast and robust optimization techniques for these types of fermentation processes. Most of the optimization studies on fermentation use detailed theoretical techniques (Chen et al., 2004; Pushpavanam et al., 1999; Johnson, 1987), while a few investigations have reported the use of direct search techniques (Iyer et al., 1999). Often, these techniques require excellent initial predictions for either control variable trajectory or unknown boundary conditions to achieve convergence. Additionally, the switching structure (feeding sequence) has to be deduced correctly. Conventionally, dynamic optimization of bioprocesses (i.e., fermentation) is carried out by applying of deterministic models, using the gradient-based methods. As usual the deterministic model is described by a set of differential equations derived from mass balances. Unlike the stochastic process models, the deterministic models have no uncertainty in the model.
The gradient-based optimization algorithm, finds the optimum point with respect to the first differentiation and second differentiation. The Calculation of the differential equations for complex systems is very complicated and sometimes becomes impossible. Common techniques that are based on mathematical optimization include Greens theorem and Pontryagins’ maximum principle. These methods have, however, limitations concerning the number of control variables and constraint handling and model complexity. The gradient-based algorithm may also get trapped at a local optimum in nonlinear systems and this is a major disadvantage of such optimization algorithms. One of the main advantages of the presented optimization approach (SA) is its ability to optimize a non-differentiable complex model. (Ronen et al., 2002).
The feeding profile of a fed-batch process was determined using a non-gradient base optimization approach. This kind of optimization is very effective on complex systems, but need extensive computational calculation (Modak, 1988).
Due to the lake of an exact mathematical model for the fermentation processes, the use of the artificial intelligence based optimization techniques is needed. Genetic Algorithm (GA), neural network, fuzzy logic and evolutionary algorithms such as simulated annealing are some of intelligent methods which belong to this class of optimization methods. The recent researches on the optimization problems are performed by applying of the intelligent techniques. The results showed that there are good agreement between the experimental data and those obtained from the intelligent techniques.  (Cheng and Wang, 2008; Valencia et al., 2007; Renard et al., 2006; Faber et al., 2005; Ronen et al., 2002; Hirafuji and Hagan, 2000; RAO et al., 1996).
In the present study, a developed optimization method is proposed for prediction of dynamic behavior of bioreactor.

Numerical aspect
Simulated annealing algorithm: SA is a strong random searching method which was first applied in difficult optimization problems by Kirkpatrick and coworkers in 1982, but this idea was first inspired by Metropolis et al., in 1953 who involved in publishing industry. SA, which analogizes annealing processes of liquids and metals to the minimization of an objective function, is based on the first law of thermodynamics. Considering a liquid metal is cooled with low speed to crystallize regularly and reaches to the minimum level of energy; this method is called “annealing cooling”. However, if frozen melted solid is made by fast cooling then it has not the optimum energy level. (Corana et al., 1999). A Boltzmann distribution (Pr) for the relationship between temperature T and decides energy E was assumed as the following equation:

    (1)
 

where KB=1.380650×10-23 m2kgs-2k-1 is the Boltzmann constant, and Z(T) is normalization factor, respectively.
The value of Pr is high at the beginning of simulation and it is decreases while going on the simulation due to the effect of a positive parameter which is named annealing temperature. Hence, the theoretic point of view predicts that SA method is able to overcome the local optimum traps and find out the global optimum by this parameter changing.  The basic structure of the simulated annealing method is illustrated in Figure 1. According to this flowchart a computer program was designed in this research. The case that has been optimized is a continuous and one dimensional optimization problem. By having initial guesses for state xk=0 and initial annealing temperature TAm=0 the quantity of the objective function jk is evaluated. With the following equation amount of new state x-k+1 will be calculated stochastically by using a random number r in the neighborhood with radius of DSl (current step size) around:

x-k+1 = xk + r DSl         "r Î[-1, 1]                         (2)
 
The acceptance condition for the new state xk+1 as an optimal value is its application regarding to the previous amount of x-k+1 in the following inequality:

j (x-k+1) £ j (xk+1) ® xk+1 = x-k+1                                    (3)
 
Otherwise, the probability of accepting x-k+1 as a new state xk+1 is given by the Metropolis criterion:

    (4)
 

Where     is a random number between [0, 1] and TAm is current annealing temperature. Steps with an increasing objective function are called up-hill steps. The procedure proceeds within the inner loop until a given number of evaluation cycles Nk have been made.
The step length (i.e., search area) is considered as an adjustable parameter in the outer loop. During the progression of the algorithm, the step size is decreased in order to contain the search region around the current value of the objective (Faber et al., 2005).
The SA algorithm chooses the step size so that the q ratio near to 0.5 (q-ratio is defined as the accepted to evaluate steps). The high q-ratio is so high it means that having long unnecessary computations. On the other hand, small q-ratio could only result in local optimum (Feehery and Barton, 1998).
As temperature TA decreases, less deteriorated objective functions are accepted by leading to a smaller acceptance ratio q, resulting in the decreasing of the step length. There are many approaches for adjusting the step length; here is a simple but efficient method (Corana et al., 1987):

 (5)  

The parameter c which controls the step variation is usually set to 1. In the SA algorithm, the temperature  is decreased after a number of step length adjustments Nl completing a “Markov chain” in the outer loop m. It is supposed the optimal state which is found at the current temperature. Multiplying the actual temperature with a constant factor is the simplest way to decrease the temperature (Farber et al., 2005):

    (6) 

The computation is now restarted and preceded with the current optimal state xopt and the temperature which has been updated. In this project work, the stopping (convergence) criterion is that the difference between the last NÎ steps and the actual optimum xopt is less than a problematic specific factor ε:

         i= 0, ....., Ne-1         (7) 

It is possible that this stopping criterion is solved before the convergence achievement and then additional convergence information should be checked. The annealing temperature and the step length are also considered in this work (Farber et al., 2005):

  
    (8) 

Bioreactor model: Models are the key components of many advanced optimal control problems. Process operation, scheduling, optimization, and control all include a model in explicit or implicit form (Modak, 1988). The considered process in this study is a fed-batch bioreactor growing S. cerevisiae on glucose as a carbon source. The selected model unlike some models which have been previously used to describe fed-batch baker’s yeast fermentation has no detailed representation of the cellular level growth process (Modak and Lim, 1987; Dutta et al., 2005). The growth process in these models is described by lumped, unstructured representation in which the biological components are evaluated by a single variable, the specific growth rate. The representative equations, as reported in literatures are:

(9) 

Where X, G, and E are the concentrations of cell mass, glucose, and ethanol, respectively. F is the volumetric feed rate; V, is the fermenter volume and GF is the glucose feed concentration. The specific rates, mG, mE, d, p, and h are summarized in Table 1. And kinetic parameters which were used in culture simulation are summarized in Table 2.
Optimization of bioreactor: Computing optimal feeding targets for the semi-batch process is an important factor in having an economical operation. However, the amount of ethanol may increase as an undesired by-product during fermentation, depending on the initial and operating conditions. According to the results of open loop simulation an increase in the ethanol concentration may occur if the feed rate of the substrate increases to achieve the maximum yield for final product as a biomass (S. cerevisiae). So the quantity and quality of the product deteriorates due to the formation of ethanol. On the other hand a slow feeding of glucose results in a longer operating time, but a lower amount of ethanol is produced and the cellular yield is increased. Therefore a profit function which weighs the amount of yeast in the operating cost is chosen for the optimization study. Generally, the biomass quantity (i.e., yeast), is the desired variable which should be maximized at the end of the fermentation. So, the objective function consists of only the biomass concentration in the simplest form. The optimization problem is made to determine the optimal glucose feed rate, F (t). The objective function, IP (index of performance), is expressed as:

 (10) 

With subscript f represents the final culture conditions and tf, the final time (operating period) was chosen as 12 h (for comparison with previous investigations). In addition it is assumed that the system is constrained so that the reactor volume at the end of the batch is restricted to 5 l, which is an end-point constraint. Moreover, the feed-rate cannot be negative or more than 1 lit/min at any time, so that:
 

   (11) 

The optimal feeding policy representing the best performance was computed for this system in three sets of initial and operation conditions which are given in Table 3. These sets of initial conditions were selected randomly for the computation algorithm so that the results were comparable to the referential results.
RESULTS

After choosing appropriate model and different random initial conditions and comparable with previous literatures, a computer program was written. Then in summary the simulations were done for various operational sets of Table 3. The results are shown in Figures 2-5. Some parameters are necessary to initializing of SA method; these parameters are given in Table 4. In addition the number of computing loop and trails and errors to become converged is presented in this Table.
The optimal glucose feed rate profiles (Figs. 3A, 4A, 5A) are computed for various initial conditions to investigate the effect of initial conditions on the optimal operation of fed-batch backer yeast culture. For all cases the internal parameters of SA do not change, whereas initial condition for differential equations changed. The discrepancies between the graphs are natural because these graphs have obtained in different conditions.
The next issue to be addressed in this work is whether there is a scope for improving the optimization algorithm robustness. The proposed model input variables such as kinetic parameters may be changing, so need to refine the model data which were used in the optimization procedure. None of unwanted problems take place when a good optimization method was applied.
A sensitivity test was studied to discus about the effect some key variables identified as the decision variables in the optimization results; some of this variables are K1 and K2 which affect specific growth rate on glucose (mG) and K5 which affects specific growth rate on ethanol (mE). The values of each of these variables were varied within a reassigned domain (10% increase and decrease) such as presented in Table 5, while the others were kept constant, to note the effects of variation on some calculated optimum profiles of initial condition I. outcome of entire analysis is presented in Figure 6.

DISCUSSION

Evaluation of iteration of SA is shown in Figure 2. On the Figure 2. The objective functions which are accepted as the new points indicate that, metropolis acceptance decreases with rising number of evaluation, which definitely shows the impact of the annealing temperature on the metropolis value. As it is indicated from the Figure 2, it is likely that during the simulation process target equation falls in the local optimum amount. But the merit of the SA method is that it does not get stuck in these conditions and searching function continues until the inequalities conditions (equations 7, 8) are satisfied and the global optimum is found.
An important case which must be noticed is that the amount of the initial guess for next time step is the answer of previous time. Therefore, the significance of the choice of the first guess temperature and the initial values of the other internal parameters in the SA are taken into accounts in order that by the lower number of calculations we could arrive in the global optimum in conclusion according to Table 4.
     The optimal glucose feed strategies for initial conditions I (low cell mass and glucose concentration) and II (higher cell mass concentration) are shown in Figures 3 and 4, which indicate the effects of the cell mass concentration. According to these results it is necessary to start the yeast fermentation process with a high glucose concentration (2-3 g/l) when the cell mass concentration is low. The results of the simulation in Figures 3B, 4B, and 5B show that by reducing the initial amount of the cell mass that was injected into the bioreactor (the case of I, II instead of III), the scope of cell mass during the fermentation lessens, and the ultimate integral amount of yeast which was produced in the bioreactor in this case is lower than the other conditions.
With the initial changes in conditions from I to II and to III and according to the Figures  3C, 4C, and 5C, it is indicated that the  peak related to the glucose profiles is getting to removed by the increase both of the initial amount of glucose and cell mass concentration  in the simulations. Consequently, what remain are the steady trends of glucose profile reduction through the time.
Ethanol is affected by these variations, in parallel. And its time range is approximately halved. Yeast cells consume the accumulated ethanol so the ethanol concentration decreases steadily. Thus, in an optimal fed-batch operation yeast cells consume both glucose and ethanol for their growth simultaneously, unlike sequential utilization in batch cultures.
Figure 5A shows the optimal feed rate profile at the initial high levels of glucose (case III, Table 3). High glucose concentration eliminates the need for glucose addition; therefore, fermentation starts with such as a batch period. The glucose concentration decreases steadily until the singular hyper surface is reached. Therefore after this time increasing in the glucose concentration has no effect on the feed rate profile and the ethanol concentration.
Feed rate at the primary times of case II, III reaches its maximum and leads to activation of constraint by equation 11.a. after approximately 3 h, ethanol reaches its maximum concentration and glucose vanishes. It leads to decrease in the feed rate culminating in constraint 11.a becoming loose.
Figure 6 show that there is just a slight difference between the former profiles and the recent profiles obtained through consideration of the model mismatch. Of course, this little difference can not be so important because this occurs in the first or last growing stage of microorganism. And it should be remembered this delicate difference is a good reason to accept the robustness of the SA method.

CONCLUSION

The optimal glucose feeding target for aerobic fermentation of S. cerevisiae in fed-batch mode was computed by applying an intelligent approach named SA. In order to investigate the effect of initial conditions on the optimal operation of culture, three various initial conditions were tested. The feeding profiles obtained for initial conditions I, II, III are similar to those obtained by Modak et al. (1988) previously. So the algorithm is capable of reaching global targets with the required computing accuracy. In comparison with traditional methods the main advantage of the presented optimization approach is its ability to optimize a non-differentiable complicated model. The sensitivity analysis of the SA optimization method was tested and resulting trends compared with the first obtained formerly in this research. Findings show slight difference between the first basic trends and the others, therefore the SA algorithm is sufficiently fast, robust and acceptable.