# Optimal Control of a Fed-Batch Fermentation Involving Multiple Feeds

Authors: Chongyang Liu; Zhaohua Gong; Zhaoyi Huo; Bangyu Shen
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

A nonlinear dynamical system, in which the feed rates of glycerol and alkali are taken as the control functions, is first proposed to formulate the fed-batch culture of 1,3-propanediol (1,3-PD) production. To maximize the 1,3-PD concentration at the terminal time, a constrained optimal control model is then presented. A solution approach is developed to seek the optimal feed rates based on control vector parametrization method and improved differential evolution algorithm. The proposed methodology yielded an increase by 32.17% of 1,3-PD concentration at the terminal time.


---

## Body

## 1. Introduction

1,3-propanediol (1,3-PD) possesses potential applications on a large commercial scale, especially as a monomer of polyesters or polyurethanes, its microbial production is recently paid attention to for its low cost, high production, and no pollution, and so forth. It is considered to be one of the bulk chemicals, which is likely to be produced by bioprocesses on large scales [1]. During the bioconversion of glycerol to 1,3-PD, the most efficient cultivation method appears to be a fed-batch culture which corrects pH by alkali addition for glycerol supply [2]. In the fed-batch cultivation process of glycerol bioconversion 1,3-PD, several factors influence the process state and govern the process behavior [3]. Along with the concentration of substrate, the process is influenced by the products of metabolism that accumulate in cultivation broth and cause an inhibition of cells’ growth. To maintain a suitable environment for cells’ growth, alkali is also intermittently added to the fermentor. Consequently, two control actions should be considered for modelling and optimal control the fed-batch cultivation process: the feeding with glycerol and the adding with alkali.Optimal control of fed-batch culture by feed rates has received extensive attention [4–17]. Initially, the batch model is extrapolated to fed-batch cultivation by incorporating the dilution factors. Unstructured and nonsegregational models with specific rates of cell growth rate, metabolite production rate, and substrate consumption rate have been used to model fed-batch fermentation. The models have been used for optimal control studies by a number of researchers [7, 9–11]. Recently, assuming the feed of substrate only occurs at the impulsive instants, nonlinear impulsive systems [4] have been extensively investigated to formulate the fermentation process. Subsequently, the properties [5], parameter identification problem [12], and optimal control problem [4, 13] have been investigated. Nonetheless, since the feed rate of substrate is finite, it is not reasonable to describe the actual fed-batch fermentation process by the impulsive system. In contrast, taking the feed rate of substrate as a time-continuous process, a nonlinear multistage was proposed in [8]. The parameter identification problem [6] and optimal control problem [8] for the system are investigated. Moreover, regarding the fed-batch fermentation as switching between the batch process and the feed process, the switched systems and their optimal control are discussed in [14–17]. Although the achieved results are interesting, the control action in the above nonlinear dynamical systems and optimal control problems only include the feed rate of the substrate.Typically, a fed-batch process could have more than one control variable that need to be optimized. This leads to a fairly complex optimization problem when the controlled nonlinear system involving more than one control variables and physical constraints on state variables. In this paper, we propose a nonlinear dynamical system, in which the feed rates of glycerol and alkali are taken as the control functions, to describe 1,3-PD production in fed-batch cultivation process. To maximize the 1,3-PD concentration at the terminal time, a constrained optimal control problem is then presented. Incidently, there exist many methods to solve the feed rates optimization problem involving more than one feed rates, such as maximum principle of Pontryagin [18], Luus-Jaakola search method [19], and genetic algorithm [10]. However, these methods are all applied to the fed-batch process in which the substrates are fed to the fermentor continuously. In the actual fermentation, glycerol and alkali are intermittently fed into fermentor. As a result, the computation is more complex and it is necessary to develop a fast and robust algorithm to solve the complex-constrained optimal control problem. A solution approach is developed to seek the optimal feed rates of glycerol and alkali based on the control vector parametrization method and an improved differential evolution algorithm. Numerical results show that the concentration of 1,3-PD at the terminal time can be increased considerably compared with previous results.The rest of the paper is organized as follows. The nonlinear dynamical system and constrained optimal control model in fed-batch fermentation process are formulated in Section2. Section 3 develops a computational approach to solve the constrained optimal control model, while Section 4 illustrates the numerical results. Finally, conclusions are provided in Section 5.


## 2. Problem Formulation

Based on the nonlinear dynamical system [6] and taking the feed rates of glycerol and alkali as the control functions, mass balances of biomass, substrate, and products in fed-batch culture can be formulated as(2.1)ẋ(t)=fi(t,x(t),u(t)),x(ti-1+)=x(ti-1),t∈(ti-1,ti],i=1,2,…,2N+1,x(0)=x0,
where x(t):=(x1(t),x2(t),x3(t),x4(t),x5(t))T∈R+5,t∈[0,T], is the state vector whose components represent the extracellular concentrations of biomass, glycerol, 1,3-PD, acetic acid, and ethanol at time t in the fermentor, respectively. T is the terminal time of the fermentation, and x0 is a given initial state. u(t):=(u1(t),u2(t))T∈R2 is the control function whose components denote the feed rates of glycerol and alkali, respectively. ti,i∈Λ:={1,2,…,2N+1} is the switching instant such that 0=t0,ti-1<ti,i∈Λandt2N+1=T. In particular, t2j+1 is the moment of adding glycerol, at which the fermentation process switches to continuous culture from batch culture, and t2j+2 denotes the moment of ending the flow of glycerol, at which the fermentation process jumps into batch culture from continuous culture, j∈Λ̅1:={0,1,2,…,N-1}. Furthermore, for t∈(t2j,t2j+1],j∈Λ̅2:={0,1,…,N},(2.2)f2j+1(t,x(t),u(t))=(μx1(t),-q2x1(t),q3x1(t),q4x1(t),q5x1(t))T,
for t∈(t2j+1,t2j+2],j∈Λ̅1,(2.3)f2j+2(t,x(t),u(t))=(μx1(t)-D(t)x1(t)D(t)(cs0-x2(t))-q2x1(t)q3x1(t)-D(t)x3(t)q4x1(t)-D(t)x4(t)q5x1(t)-D(t)x5(t)).
In (2.3), cs0 denotes the initial concentration of glycerol in feed. D(t) is the dilution rate at time t defined by(2.4)D(t)=u1(t)+u2(t)V(t),(2.5)V(t)=V0+∫0t(u1(s)+u2(s))ds.
In (2.5), V0 is the initial volume of solution in the fermentor. On the basis of the previous work [20], the specific growth rate of cells μ is expressed by(2.6)μ=Δ1x2(t)x2(t)+k1∏l=25(1-xl(t)xl*)nl,
where Δ1 is the maximum specific growth rate; k1 is the Monod saturation constant; xℓ* are the maximal residual substrate and products concentrations; and nℓ are the exponents for the substrate and products. The specific consumption rate of substrate q2 is(2.7)q2=m2+μY2+Δ2x2(t)x2(t)+k2.
In (2.7), m2 is the maintenance term of substrate consumption under substrate-limited conditions. Y2 is the maximum growth yield. Δ2 is the maximum increment of substrate consumption rate under substrate-sufficient conditions. k2 is the saturation constant for substrate. The specific formation rates qℓ,ℓ=3,4, of 1,3-PD and acetic acid are defined as(2.8)ql=ml+μYl+Δlx2(t)x2(t)+kl,
where mℓ are the maintenance terms of product formations under substrate-limited conditions; Yℓ are the maximum product yields; Δℓ are the maximum increments of product formation rates under substrate-sufficient conditions; kℓ are saturation constants for products. Moreover, the specific formation rate q5 of ethanol can be described by(2.9)q5=q2(c1c2+μx2(t)+c3c4+μx2(t)),
in which c1,c2,c3, and c4 are parameters for determination of yield of ethanol on glycerol.Under anaerobic conditions at 37°C and pH 7.0, the critical concentrations for cells growth and the parameters in (2.6)–(2.9) are listed in Table 1.Table 1
The parameters and the values of critical concentrations in the controlled nonlinear system (2.1).
ℓmℓYℓΔℓkℓnℓcℓx*ℓxℓ*1——0.670.28—0.0250.01622.200.008828.5811.4310.061520393−2.6967.6926.5915.5035.18010364−0.9733.075.7485.71350.45010265————3—060.9Define(2.10)U:={u(t)∣aı,i≤uı(t)≤bı,i,ı=1,2,t∈(ti-1,ti],i∈Λ},
where a1,2j+1,a2,2j+1,b1,2j+1 and b2,2j+1,j∈Λ̅2, are identically equal to zero. a1,2j+2,a2,2j+2,b1,2j+2, and b2,2j+2,j∈Λ̅1 are positive constants which denote the minimal and maximal rates of adding glycerol and alkali, respectively. Let 𝒰 be the class of all the measurable functions u from [0,T] into R2 with u(t)∈U.There exist critical concentrations, outside which cells cease to growth, of biomass, glycerol, 1,3-PD, acetic acid, and ethanol. Hence, it is biologically meaningful to restrict the concentrations of biomass, glycerol, products, and the volume of culture fluid in a setW defined as(2.11)xT(t)∈W:=∏l=15[x*l,xl*],t∈[0,T].For the system (2.1), we confirm that the system has a unique continuous solution on [0,T], denoted by x(·;u), for each u∈𝒰. Moreover, the solution x(·;u) is continuous in u and is uniformly bounded on [0,T] [21].In fed-batch fermentation of glycerol bioconversion to 1,3-PD, controlling the feed rates of glycerol and alkali is to maximize the final 1,3-PD productivity. Therefore, the 1,3-PD concentration at the terminal time is taken as the cost functional, that is,(2.12)J(u):=x3(T;u),
where x3(·;u) is the third component of the solution to the system (2.1). Furthermore, the control constraint (2.10) and the state constraint (2.11) should be satisfied.Thus, a constrained optimal control model in the fed-batch fermentation can be formulated as follows.(2.13)(COCP)max⁡J(u)s.t.xT(t;u)∈W,t∈[0,T],u∈U.


## 3. Computational Methods

In solving the optimal control problem, a particularly attractive technique variously proposed [22, 23] is the parametrization of the control vector u into a number of intervals so as to result in an expression as (3.1)uikp(t)=∑j=1niσijkpϕij(t),t∈(ti-1,ti],i∈{1,2,…,2N+1},k∈{1,2},
that is, let the time subinterval [ti-1,ti] be partitioned into ni subintervals with ni+1 partition points denoted by (3.2)τ0i,τ1i,…,τnii,τ0i=ti-1,τnii=ti.
In the above equation, ϕij(t) are known trial functions that need to be chosen. uikp(t) is the value of the kth control variable in ith interval and σijkp is the coefficient that need to be determined from optimization. A number of trial functions based on the use of different basis functions have been used in various work. However, the choice of trial functions ϕij(t) affects the evaluation of the exact optimization and needs to be carefully chosen. In this paper, we assume that each control variable in the control vector u(t) is approximated by a piecewise constant value in a particular interval i, namely, (3.3)ϕij(t)={1,t∈(τj-1i,τji],0,otherwise.Letx(·;σp) be the solution of the system (2.1) replacing the control function u with up. Furthermore, denote the set of the admissible parameterized controls σp by Ξp, namely, Ξp={σp∣up∈𝒰}. Now, (COCP) can be approximated by the following parameter optimization problem (3.4)(COCP(p))max⁡J̃(σp):=x3(T;σp)s.t.xT(t;σp)∈W,t∈[0,T],σp∈Ξp.To solve (COCP) by the control vector parametrization method, we need to solve a sequence of problems{(COCP(p))}p=1∞. However, it is difficult to cope with the state constraints. To surmount these difficulties, let (3.5)gl(x(t;σp)):=xl(t;σp)-xl*,g5+l(x(t;σp)):=x*l-xl(t;σp),l=1,2,…,5.
Then, the condition (2.11) is equivalently transcribed into (3.6)G(σp)=0,
where G(σp):=∑l=110∫0Tmax{0,gl(x(t;σp))}dt. However, G(σp) is nonsmooth in σp. Consequently, standard optimization routines would have difficulties in dealing with this type of equality constraints. The following smoothing technique is to replace max{0,gl(t;σp)} with ĝl,ϵ(t;σp), where (3.7)ĝl,ϵ(x(t;σp)):={0,ifgl(x(t;σp))<-ϵ,(gl(x(t;σp))+ϵ)24ϵ,if-ϵ≤gl(x(t;σp))≤ϵ,gl(x(t;σp)),ifgl(x(t;σp))>ϵ.
In (3.7), ϵ>0 is an adjustable parameter controlling the accuracy of the approximation. Note that (3.8)Ĝϵ(σp):=∑l=110∫0Tĝl,ϵ(x(t;σp))dt
is a smooth function in θ. The equality constraint (3.3) now can be approximated by (3.9)Ĝϵ(σp)=0.As a result, (COCP(p)) can be approximated by the following approximation problem as follows. (3.10)(COCPϵ(p))max⁡⁡Ĵ(σp):=x3(T;σp)s.t.Ĝϵ(σp)=0,σp∈Ξp.Then, the gradient of the constraintĜϵ(·) can be computed by the following theorem.Theorem 3.1.
For the constraintĜϵ(σp) given in (3.9), it holds that its gradient with respect to parameterized control σp is
(3.11)∂Ĝϵ(σp)∂σp=∫0T∂H(t,x(t),σp,λ(t))∂σpdt,
where
(3.12)H(t,x(t),σp,λ(t))=∑l=110ĝl,ϵ(x(t;σp))+λT(t)f̃(t,x(t),σp),λ(t)=(λ1(t),λ2(t),λ3(t),λ4(t),λ5(t))T
is the solution of the costate system
(3.13)λ̇(t)=-(∂H(t,x(t),σp,λ(t))∂x)T,
with the boundary conditions
(3.14)λ(T)=(0,0,1,0,0)T,λ(tι+)=λ(tι-),ι=1,2,…,2N.Proof.
The proof can be completed using the method of Chapter  3 in [24].On this basis, (COCP) can be solved by a sequence of approximation problems{(COCPϵ(p))}. Each of these {(COCPϵ(p))} is viewed as a smooth nonlinear mathematical programming problem solved by various optimization methods such as gradient-based techniques [23]. However, all those techniques are only designed to find local optimal solutions. Furthermore, in solving these {(COCPϵ(p))}, the evaluation of candidate feed rates is a computationally expensive operation because of solving the controlled nonlinear system (2.1). Consequently, finding the global optimum or a good suboptimal solution with traditional algorithms based on gradient method or nontraditional search and optimization techniques based on natural phenomenon such as genetic algorithm [25], evolution strategies [26], and simulation annealing [27] is too consuming, or even impossible within the time available.Differential evolution (DE), a recent optimization technique, is an exceptionally simple and easy to use evolution strategy, which is significantly faster and robust at numerical optimization and is more likely to find a function’s true global optimum [28]. Differential evolution algorithm has been used in the recent past to solve many engineering optimization problem, see, for example, [29, 30]. When using any population-based search algorithm in general and DE in particular to optimize a function, an acceptable trade-off between convergence rate (with reference to locating optimum) and robustness (with reference to not missing the global optima) must generally be determined. To increase the convergence speed of DE without compromising with the robustness, a modified differential evolution (MDE) is developed to solve nonlinear unconstrained optimization problems encountered in chemical engineering [31]. Nevertheless, the (COCPϵ(p))is a nonlinear optimization problem with constraints in state and control parameters, which MDE can not be applied directly to solving it. Hence, the following strategies are added to the MDE [31]. (i)
(Handling the control constraint). If there is bound violation for a parameter in theιth individual at the κth step, then that parameter is generated randomly between given lower and upper bound using the following equation:
(3.15)σιjp(κ)=lower(σjp)+randι[0,1]×(upper(σjp)-lower(σjp)),j=1,…,D,
where is D is the number of parameters.(ii)
(Dealing with the continuous state constraint). For the parameter of theιth individual at the κth step, test the value of Ĝϵ(σιp(κ)). If Ĝϵ(σιp(κ))=0, then the parameter is feasible. Otherwise, that is, Ĝϵ(σιp(κ))>0, move the parameter towards the feasible region in the direction of -(∂Ĝϵ(σιp(κ))/∂σιp(κ)) with Armijo line search.(iii)
(Stopping criteria). The algorithm stops when any of the following conditions holds:
(a)
the maximal iterationM is reached;(b)
the maximal deviation between the group’s best fitness values in the lastM1 iterations is less than ε, where ε is a predefined constant.Now, we can obtain an approximately optimal control for (COCP) as shown in the following algorithm.Algorithm 1.
Step 1.
Choose initial value ofϵ.Step 2.
Solve approximate problem(COCPϵ(σp)) using the improved MDE algorithm to give σϵp,*.Step 3.
Setϵ=αϵ. If ϵ>ϵ̅, where ϵ̅ is a prespecified positive constant, go to Step 2. Otherwise go to Step 4.Step 4.
Ifmini∈{1,2,…,2N+1}npi≥P, where P is a predefined positive constant, go to Step 5. Otherwise, go to Step 1 with npi increased to npi+1 for each i.Step 5.
Constructup,* from σϵp,* by (3.1) and stop.The piecewise constant controlup,* obtained is an approximately optimal solution of (COCP).


## 4. Numerical Results

The initial state, the initial volume of fermentor, the initial concentration of glycerol in feed, and the fermentation time arex0=(0.1115gL-1,495mmolL-1,0,0,0)T, V0=5L,cs0=10762mmolL-1, and T=24.16h, respectively. Fed-batch process begins at t1=5.33h. The feeding moment t2j+1 and the feeding stopping moment t2j+2, j∈Λ̅1={0,1,…,676} are determined by the experiment.In order to save computational time, the fermentation process is partitioned into the first batch phase (Bat.Ph.) and phases I-IX (Phs.I-IX) according to the actual fermentation process. In each of Phs.I-IX, the same feeding strategy is adopted. Moreover, the durations of the feeding processes in Phs. I-IX are 5, 7, 8, 7, 6, 4, 3, 2, 1 seconds in each 100 seconds, leaving 95, 93, 92, 93, 96, 97, 98, and 99 seconds for batch cultures, respectively. It should be mentioned that this approach had been adopted to obtain the experimental data [12]. In addition, the bounds of feed rates of glycerol and alkali in Phs.I-IX are listed in Table 2.Table 2
The bounds of feed rates in Phs.I-IX.
PhasesUpper (u1)Lower (u1)Upper (u2)Lower (u2)Phs.I-II0.25240.16820.19050.12615Ph.III0.23900.15940.179250.11955Phs.IV-V0.25240.16820.19050.12615Ph.VI0.26570.17710.1992750.132825Ph.VII0.29240.19490.21930.146175Phs.VIII-IX0.30580.20380.229350.15285In the improved MDE algorithm, the number of populationsNp, the maximal iteration M, and the parameters F, CR, M1, ε, ϵ, ϵ̅ are, respectively, 86, 300, 0.5, 0.9, 10, 10-3, 0.1, and 1.0×10-8. The specified constant P in Algorithm 1 is 1. These parameters are derived empirically after numerous experiments.Applying Algorithm1 to the optimal control model, we obtain the optimal feed strategies of glycerol and alkali. Here, the computational process is coded using Microsoft Visual C++ on PC with Intel Core 2 Duo CPU, 2.93 GHz RAM. Especially, the ODEs are numerically calculated by the improved Euler method with the relative error tolerance 10-4. The optimal feed rates of glycerol and alkali are plotted using MATLAB R2010a in Figure 1. To show the feed rates of glycerol and alkali in the feeding processes for Ph. I to Ph. IX better, 9 small subfigures are also incorporated in Figure 1, respectively. Under the obtained optimal feed rates, the computational concentration of 1,3-PD at the terminal time is 1053.67mmolL-1 which is increased by 32.17% in comparison with experimental result 797.23mmolL-1. In particular, the concentration change of 1,3-PD with respect to fermentation time under the optimal feed rates of glycerol and alkali is shown in Figure 2. From these curves in Figure 2, we conclude that 1,3-PD concentration at the terminal time in this paper is actually higher than the one previously reported.The optimal feed strategies of 1,3-PD and alkali in fed-batch fermentation process.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 2
The concentration changes of 1,3-PD with respect to fermentation time.


## 5. Conclusions

Taking the feed rates of glycerol and alkali as the control functions, a nonlinear dynamical system was proposed to describe a fed-batch fermentation. A constrained optimal control model involving the nonlinear dynamical system was then presented to maximize the concentration of 1,3-PD at the terminal time. Based on the control vector parametrization method and the improved MDE algorithm, a computational method was developed to seek the optimal solution of the optimal control model. Numerical results verified the validity of the mathematical model and the effectiveness of the computational method.



---

*Source: 245315-2012-03-08.xml*