Optimal Control of a Fed-Batch Fermentation Involving Multiple Feeds

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.

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 Section 2. 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𝑥̇(𝑡)=𝑓^𝑖((𝑡)(𝑡,𝑥𝑡𝑡,𝑥,𝑢))_𝑖−1+=𝑥𝑡_𝑖−1,𝑡∈𝑡_𝑖−1,𝑡_𝑖𝑥(0),𝑖=1,2,…,2𝑁+1,=𝑥₀, (2.1) where 𝑥(𝑡)∶=(𝑥₁(𝑡),𝑥₂(𝑡),𝑥₃(𝑡),𝑥₄(𝑡),𝑥₅(𝑡))^𝑇∈𝑅⁵₊,𝑡∈[0,𝑇], is the state vector whose components represent the extracellular concentrations of biomass, glycerol, 1,3-PD, acetic acid, and ethanol at time 𝑡 in the fermentor, respectively. 𝑇 is the terminal time of the fermentation, and 𝑥₀ is a given initial state. 𝑢(𝑡)∶=(𝑢₁(𝑡),𝑢₂(𝑡))^𝑇∈𝑅² is the control function whose components denote the feed rates of glycerol and alkali, respectively. 𝑡_𝑖,𝑖∈Λ∶={1,2,…,2𝑁+1} is the switching instant such that 0=𝑡₀,𝑡_𝑖−1<𝑡_𝑖,𝑖∈Λand𝑡_2𝑁+1=𝑇. In particular, 𝑡_2𝑗+1 is the moment of adding glycerol, at which the fermentation process switches to continuous culture from batch culture, and 𝑡_2𝑗+2 denotes the moment of ending the flow of glycerol, at which the fermentation process jumps into batch culture from continuous culture, 𝑗∈Λ₁∶={0,1,2,…,𝑁−1}. Furthermore, for 𝑡∈(𝑡_2𝑗,𝑡_2𝑗+1],𝑗∈Λ₂∶={0,1,…,𝑁},𝑓^2𝑗+1(𝑡,𝑥(𝑡),𝑢(𝑡))=𝜇𝑥₁(𝑡),−𝑞₂𝑥₁(𝑡),𝑞₃𝑥₁(𝑡),𝑞₄𝑥₁(𝑡),𝑞₅𝑥₁(𝑡)^𝑇, (2.2) for 𝑡∈(𝑡_2𝑗+1,𝑡_2𝑗+2],𝑗∈Λ₁,𝑓^2𝑗+2(𝑡,𝑥(𝑡),𝑢(𝑡))=𝜇𝑥₁(𝑡)−𝐷(𝑡)𝑥₁(𝑡)𝐷(𝑡)𝑐_𝑠0−𝑥₂(𝑡)−𝑞₂𝑥₁(𝑡)𝑞₃𝑥₁(𝑡)−𝐷(𝑡)𝑥₃(𝑡)𝑞₄𝑥₁(𝑡)−𝐷(𝑡)𝑥₄(𝑡)𝑞₅𝑥₁(𝑡)−𝐷(𝑡)𝑥₅(𝑡). (2.3) In (2.3), 𝑐_𝑠0 denotes the initial concentration of glycerol in feed. 𝐷(𝑡) is the dilution rate at time 𝑡 defined by𝐷(𝑡)=𝑢₁(𝑡)+𝑢₂(𝑡)𝑉(𝑡),𝑉(𝑡)=𝑉₀+^𝑡₀𝑢₁(𝑠)+𝑢₂(𝑠)𝑑𝑠. (2.4)(2.5) In (2.5), 𝑉₀ 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𝑥1−^ℓ(𝑡)𝑥^∗_ℓ^𝑛ℓ, (2.6) where Δ₁ is the maximum specific growth rate; 𝑘₁ is the Monod saturation constant; 𝑥^∗_ℓ are the maximal residual substrate and products concentrations; and 𝑛_ℓ are the exponents for the substrate and products. The specific consumption rate of substrate 𝑞₂ is𝑞₂=𝑚₂+𝜇𝑌₂+Δ₂𝑥₂(𝑡)𝑥₂(𝑡)+𝑘₂. (2.7) In (2.7), 𝑚₂ is the maintenance term of substrate consumption under substrate-limited conditions. 𝑌₂ is the maximum growth yield. Δ₂ is the maximum increment of substrate consumption rate under substrate-sufficient conditions. 𝑘₂ is the saturation constant for substrate. The specific formation rates 𝑞_ℓ,ℓ=3,4, of 1,3-PD and acetic acid are defined as𝑞_ℓ=𝑚_ℓ+𝜇𝑌_ℓ+Δ_ℓ𝑥₂(𝑡)𝑥₂(𝑡)+𝑘_ℓ, (2.8) where 𝑚_ℓ are the maintenance terms of product formations under substrate-limited conditions; 𝑌_ℓ are the maximum product yields; Δ_ℓ are the maximum increments of product formation rates under substrate-sufficient conditions; 𝑘_ℓ are saturation constants for products. Moreover, the specific formation rate 𝑞₅ of ethanol can be described by𝑞₅=𝑞₂𝑐₁𝑐₂+𝜇𝑥₂(𝑡)+𝑐₃𝑐₄+𝜇𝑥₂(𝑡), (2.9) in which 𝑐₁,𝑐₂,𝑐₃, and 𝑐₄ 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).

Define𝑈∶=𝑢(𝑡)∣𝑎_𝚤,𝑖≤𝑢_𝚤(𝑡)≤𝑏_𝚤,𝑖,𝚤=1,2,𝑡∈𝑡_𝑖−1,𝑡_𝑖,𝑖∈Λ, (2.10) where 𝑎_1,2𝑗+1,𝑎_2,2𝑗+1,𝑏_1,2𝑗+1 and 𝑏_2,2𝑗+1,𝑗∈Λ₂, are identically equal to zero. 𝑎_1,2𝑗+2,𝑎_2,2𝑗+2,𝑏_1,2𝑗+2, and 𝑏_2,2𝑗+2,𝑗∈Λ₁ 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 𝑢 from [0,𝑇] into 𝑅² with 𝑢(𝑡)∈𝑈.

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 set 𝑊 defined as𝑥^𝑇(𝑡)∈𝑊∶=⁵_ℓ=1𝑥_∗ℓ,𝑥^∗_ℓ,𝑡∈[0,𝑇]. (2.11)

For the system (2.1), we confirm that the system has a unique continuous solution on [0,𝑇], denoted by 𝑥(⋅;𝑢), for each 𝑢∈𝒰. Moreover, the solution 𝑥(⋅;𝑢) is continuous in 𝑢 and is uniformly bounded on [0,𝑇] [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) where 𝑥₃(⋅;𝑢) 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. (COCP)max𝐽(𝑢)s.t.𝑥^𝑇(𝑡;𝑢)∈𝑊,𝑡∈[0,𝑇],𝑢∈𝒰. (2.13)

3. Computational Methods

In solving the optimal control problem, a particularly attractive technique variously proposed [22, 23] is the parametrization of the control vector 𝑢 into a number of intervals so as to result in an expression as 𝑢^𝑝_𝑖𝑘(𝑡)=^𝑛𝑖_𝑗=1𝜎^𝑝_𝑖𝑗𝑘𝜙_𝑖𝑗(𝑡),𝑡∈𝑡_𝑖−1,𝑡_𝑖,𝑖∈{1,2,…,2𝑁+1},𝑘∈{1,2}, (3.1) that is, let the time subinterval [𝑡_𝑖−1,𝑡_𝑖] be partitioned into 𝑛_𝑖 subintervals with 𝑛_𝑖+1 partition points denoted by 𝜏^𝑖₀,𝜏^𝑖₁,…,𝜏^𝑖_𝑛𝑖,𝜏^𝑖₀=𝑡_𝑖−1,𝜏^𝑖_𝑛𝑖=𝑡_𝑖. (3.2) In the above equation, 𝜙_𝑖𝑗(𝑡) are known trial functions that need to be chosen. 𝑢^𝑝_𝑖𝑘(𝑡) is the value of the 𝑘th control variable in 𝑖th interval and 𝜎^𝑝_𝑖𝑗𝑘 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 𝜙_𝑖𝑗(𝑡) 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 𝑢(𝑡) is approximated by a piecewise constant value in a particular interval 𝑖, namely, 𝜙_𝑖𝑗(𝑡)=1,𝑡∈𝜏^𝑖_𝑗−1,𝜏^𝑖_𝑗,0,otherwise. (3.3)

Let 𝑥(⋅;𝜎^𝑝) be the solution of the system (2.1) replacing the control function 𝑢 with 𝑢^𝑝. Furthermore, denote the set of the admissible parameterized controls 𝜎^𝑝 by Ξ^𝑝, namely, Ξ^𝑝={𝜎^𝑝∣𝑢^𝑝∈𝒰}. Now, (COCP) can be approximated by the following parameter optimization problem COCP𝑝max𝐽̃𝜎^𝑝∶=𝑥₃𝑇;𝜎^𝑝s.t.𝑥^𝑇𝑡;𝜎^𝑝[],𝜎∈𝑊,𝑡∈0,𝑇_𝑝∈Ξ^𝑝. (3.4)

To solve (COCP) by the control vector parametrization method, we need to solve a sequence of problems {(COCP(𝑝))}^∞_𝑝=1. However, it is difficult to cope with the state constraints. To surmount these difficulties, let 𝑔_ℓ𝑥𝑡;𝜎^𝑝∶=𝑥_ℓ𝑡;𝜎^𝑝−𝑥^∗_ℓ,𝑔_5+ℓ𝑥𝑡;𝜎^𝑝∶=𝑥_∗ℓ−𝑥_ℓ𝑡;𝜎^𝑝,ℓ=1,2,…,5. (3.5) Then, the condition (2.11) is equivalently transcribed into 𝐺𝜎^𝑝=0, (3.6) where 𝐺(𝜎^𝑝)∶=∑¹⁰_𝑙=1∫^𝑇₀max{0,𝑔_𝑙(𝑥(𝑡;𝜎^𝑝))}𝑑𝑡. However, 𝐺(𝜎^𝑝) is nonsmooth in 𝜎^𝑝. Consequently, standard optimization routines would have difficulties in dealing with this type of equality constraints. The following smoothing technique is to replace max{0,𝑔_𝑙(𝑡;𝜎^𝑝)} with 𝑔_𝑙,𝜖(𝑡;𝜎^𝑝), where 𝑔_𝑙,𝜖𝑥𝑡;𝜎^𝑝∶=0,if𝑔_𝑙𝑥𝑡;𝜎^𝑝𝑔<−𝜖,_𝑙𝑥𝑡;𝜎^𝑝+𝜖²4𝜖,if−𝜖≤𝑔_𝑙𝑥𝑡;𝜎^𝑝𝑔≤𝜖,_𝑙𝑥𝑡;𝜎^𝑝,if𝑔_𝑙𝑥𝑡;𝜎^𝑝>𝜖. (3.7) In (3.7), 𝜖>0 is an adjustable parameter controlling the accuracy of the approximation. Note that 𝐺_𝜖𝜎^𝑝∶=¹⁰_𝑙=1^𝑇₀𝑔_𝑙,𝜖𝑥𝑡;𝜎^𝑝𝑑𝑡 (3.8) is a smooth function in 𝜃. The equality constraint (3.3) now can be approximated by 𝐺_𝜖𝜎^𝑝=0. (3.9)

As a result, (COCP(𝑝)) can be approximated by the following approximation problem as follows. COCP_𝜖𝑝max𝐽̂𝜎^𝑝∶=𝑥₃𝑇;𝜎^𝑝s.t.𝐺_𝜖𝜎^𝑝𝜎=0,_𝑝∈Ξ^𝑝. (3.10)

Then, the gradient of the constraint 𝐺_𝜖(⋅) can be computed by the following theorem.

Theorem 3.1. For the constraint 𝐺_𝜖(𝜎^𝑝) given in (3.9), it holds that its gradient with respect to parameterized control 𝜎^𝑝 is 𝜕𝐺_𝜖𝜎^𝑝𝜕𝜎^𝑝=^𝑇₀𝜕𝐻𝑡,𝑥(𝑡),𝜎^𝑝,𝜆(𝑡)𝜕𝜎^𝑝𝑑𝑡, (3.11) where 𝐻𝑡,𝑥(𝑡),𝜎^𝑝,𝜆(𝑡)=¹⁰_𝑙=1𝑔_𝑙,𝜖𝑥𝑡;𝜎^𝑝+𝜆^𝑇(𝑡)𝑓𝑡,𝑥(𝑡),𝜎^𝑝,𝜆(𝑡)=𝜆₁(𝑡),𝜆₂(𝑡),𝜆₃(𝑡),𝜆₄(𝑡),𝜆₅(𝑡)^𝑇 (3.12) is the solution of the costate system 𝜆̇(𝑡)=−𝜕𝐻𝑡,𝑥(𝑡),𝜎^𝑝,𝜆(𝑡)𝜕𝑥^𝑇, (3.13) with the boundary conditions 𝜆(𝑇)=(0,0,1,0,0)^𝑇,𝜆𝑡_𝜄+=𝜆𝑡_𝜄−,𝜄=1,2,…,2𝑁. (3.14)

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_𝜖(𝑝))}. Each of these {(COCP_𝜖(𝑝))} 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_𝜖(𝑝))}, 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_𝜖(𝑝))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: 𝜎^𝑝_𝜄𝑗(𝜅)=lower𝜎^𝑝_𝑗+rand^𝜄[0,1]×upper𝜎^𝑝_𝑗−lower𝜎^𝑝_𝑗,𝑗=1,…,𝐷, (3.15) where is 𝐷 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 𝐺_𝜖(𝜎^𝑝_𝜄(𝜅)). If 𝐺_𝜖(𝜎^𝑝_𝜄(𝜅))=0, then the parameter is feasible. Otherwise, that is, 𝐺_𝜖(𝜎^𝑝_𝜄(𝜅))>0, move the parameter towards the feasible region in the direction of −(𝜕𝐺_𝜖(𝜎^𝑝_𝜄(𝜅))/𝜕𝜎^𝑝_𝜄(𝜅)) with Armijo line search. (iii)(Stopping criteria). The algorithm stops when any of the following conditions holds: (a)the maximal iteration 𝑀 is reached; (b)the maximal deviation between the group’s best fitness values in the last 𝑀₁ 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_𝜖(𝜎^𝑝)) using the improved MDE algorithm to give 𝜎_𝜖^𝑝,∗.Step 3. Set 𝜖=𝛼𝜖. If 𝜖>𝜖, where 𝜖 is a prespecified positive constant, go to Step 2. Otherwise go to Step 4. Step 4. If min_{𝑖∈{1,2,…,2𝑁+1}}𝑛_𝑝𝑖≥𝑃, where 𝑃 is a predefined positive constant, go to Step 5. Otherwise, go to Step 1 with 𝑛_𝑝𝑖 increased to 𝑛_𝑝𝑖+1 for each 𝑖. Step 5. Construct 𝑢^𝑝,∗ from 𝜎_𝜖^𝑝,∗ by (3.1) and stop.

The piecewise constant control 𝑢^𝑝,∗ 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 are 𝑥₀=(0.1115gL⁻¹,495mmolL⁻¹,0,0,0)^𝑇, 𝑉₀=5L,𝑐_𝑠0=10762mmolL⁻¹, and 𝑇=24.16h, respectively. Fed-batch process begins at 𝑡₁=5.33h. The feeding moment 𝑡_2𝑗+1 and the feeding stopping moment 𝑡_2𝑗+2, 𝑗∈Λ₁={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.

In the improved MDE algorithm, the number of populations 𝑁𝑝, the maximal iteration 𝑀, and the parameters 𝐹, 𝐶𝑅, 𝑀₁, 𝜀, 𝜖, 𝜖 are, respectively, 86, 300, 0.5, 0.9, 10, 10⁻³, 0.1, and 1.0×10⁻⁸. The specified constant 𝑃 in Algorithm 1 is 1. These parameters are derived empirically after numerous experiments.

Applying Algorithm 1 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⁻⁴. 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⁻¹ which is increased by 32.17% in comparison with experimental result 797.23mmolL⁻¹. 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.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 1 

The optimal feed strategies of 1,3-PD and alkali in fed-batch fermentation process.

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.

Acknowledgments

This paper is supported by the Shandong Province Natural Science Foundation of China (Grant no. ZR2010AQ016 and Grant no. ZR2011AL003), the Natural Science Foundation for the Youth of China (Grant no. 11001153), and the TianYuan Special Funds of the National Natural Science Foundation of China (Grant no. 11126077).