Optimal Control for Multistage Nonlinear Dynamic System of Microbial Bioconversion in Batch Culture

Abstract

In batch culture of glycerol biodissimilation to 1,3-propanediol (1,3-PD), the aim of adding glycerol is to obtain as much 1,3-PD as possible. Taking the yield intensity of 1,3-PD as the performance index and the initial concentration of biomass, glycerol, and terminal time as the control vector, we propose an optimal control model subject to a multistage nonlinear dynamical system and constraints of continuous state. A computational approach is constructed to seek the solution of the above model. Firstly, we transform the optimal control problem into the one with fixed terminal time. Secondly, we transcribe the optimal control model into an unconstrained one based on the penalty functions and an extension of the state space. Finally, by approximating the control function with simple functions, we transform the unconstrained optimal control problem into a sequence of nonlinear programming problems, which can be solved using gradient-based optimization techniques. The convergence analysis and optimality function of the algorithm are also investigated. Numerical results show that, by employing the optimal control, the concentration of 1,3-PD at the terminal time can be increased, compared with the previous results.

1. Introduction

The bioconversion of glycerol to 1,3-propanediol (1,3-PD) has recently received more and more attention throughout the world due to its environmental safety, high region specificity, cheaply available feedstock, and relatively high theoretical molar yield [1]. Many researchs have been carried out including the quantitative description of the cell growth kinetics of multiple inhibitions, the metabolic overflow kinetics of substrate consumption and product formation [2–4], open-loop substrate input and pH logic control [5], enzyme-catalytic reductive pathway and transport of glycerol and 1,3-propanediol across cell membrane [6], parameter identification of biochemical systems [7] and feedback control and pulse feeding [8] for the models of the continuous cultures, feeding strategy of glycerol [9], and optimal control [10] and optimality condition [11] in fed-batch culture.

Compared with continuous and feed-batch cultures, glycerol fermentation in batch culture can obtain the highest production concentration and molar yield 1,3-PD to glycerol [12]. So nonlinear dynamical systems in this culture have been extensively considered in recent years [13–15]. In batch culture of glycerol biodissimilation to 1,3-propanediol (1,3-PD), the aim of adding glycerol is to obtain as much 1,3-PD as possible. In this paper, based on the previous model in [16], taking the yield intensity of 1,3-PD as the performance index and the initial concentration of biomass, glycerol and terminal time as the control vector, we propose an optimal control model subject to a multistage nonlinear dynamical system and constraints of continuous state. A computational approach is constructed to seek the solution of the above model in two aspects. On the one hand transform the optimal control problem into the one with fixed terminal time and transcribe it into an unconstrained one based on the penalty functions and an extension of the state space; on the other hand, by approximating the control function with simple functions, we transform the unconstrained optimal control problem into a sequence of nonlinear programming problems, which can be solved using gradient-based optimization techniques. The convergence analysis and optimality function of the algorithm are also investigated. Numerical results show that, by employing the optimal control, the concentration of 1,3-PD at the terminal time can be increased, compared with the previous results.

This paper is organized as follows. In Section 2, a nonlinear dynamical system of batch culture is proposed. In Section 3, we propose an optimal control model, develop a computational approach to solve the optimal control model, and prove the convergence of algorithm. Section 4 illustrates the numerical results. Finally, conclusions are provided in Section 5.

2. Nonlinear Dynamical System

On the basis of our previous literature(see [16]), mass balances of biomass, substrate, and products in batch culture can be formulated as the following nonlinear dynamical system:𝑥̇(𝑡)=𝑓(𝑡,𝑥(𝑡)),𝑡∈0,𝑡_𝑓(0)𝑓((𝑡=𝑓,𝑥=𝜉,𝑡,𝑥))₁(𝑡,𝑥(𝑡)),𝑓₂(𝑡,𝑥(𝑡)),𝑓₃(𝑡,𝑥(𝑡)),𝑓₄(𝑡,𝑥(𝑡)),𝑓₅(𝑡,𝑥(𝑡))^𝑇=𝜇𝑥₁(𝑡),−𝑞₂𝑥₁(𝑡),𝑞₃𝑥₁(𝑡),𝑞₄𝑥₁(𝑡),𝑞₅𝑥₁(𝑡)^𝑇, (2.1)(2.2) where 𝑥₁(𝑡), 𝑥₂(𝑡), 𝑥₃(𝑡), 𝑥₄(𝑡), and 𝑥₅(𝑡) are biomass, glycerol, 1,3-PD, acetate, and ethanol concentrations at time 𝑡 in the reactor, respectively. 𝜉=(𝜉₁,𝜉₂,𝜉₃,𝜉₄,𝜉₅)^𝑇 denotes the initial state, and 𝑡_𝑓 is the terminal time of the fermentation process. 𝑥≜(𝑥₁,𝑥₂,𝑥₃,𝑥₄,𝑥₅)^𝑇∈𝑅⁵₊ is as state vector. The specific growth rate of cells 𝜇, specific consumption rate of substrate, 𝑞₂ and specific formation rate of products 𝑞_𝑖, 𝑖=3,4,5, are expressed by the following equations on the basis of [13, 16]:𝜇=𝜇_𝑚−exp𝑡−𝑡_𝑚²2𝑡²_𝑙⁵_𝑖=2𝑥1−^𝑖𝑥^∗_𝑖,𝑞₂=𝑚₂+𝜇𝑌₂,𝑞_𝑖=𝑚_𝑖+𝜇𝑌_𝑖,𝑖=3,4,5. (2.3)

In batch culture, the initial concentrations of biomass, glycerol, and the terminal time can be chosen as control variables. Let 𝑢=(𝑢₁,𝑢₂,𝑢₃)^𝑇≜(𝜉₁,𝜉₂,𝑡_𝑓)^𝑇∈𝑅³₊ be the control vector. The solution of system (2.1) with respect to control vector is defined by 𝑥(⋅,𝑢).

Based on the factual fermentation, there exist critical concentrations, outside which cells cease to grow, of biomass, glycerol, 1,3-PD, acetate and ethanol. Hence, it is biologically meaningful to restrict the concentrations of biomass, glycerol, products, and the volume of culture fluid in a set 𝑊 and the control vector in a admissible control set 𝑈 defined respectively, as follows:𝑥(𝑡,𝑢)∈𝑊≜⁵_𝑖=1𝑥_∗𝑖,𝑥^∗_𝑖⊂𝑅⁵₊,∀𝑡∈𝐼=0,𝑡_𝑓.𝑢∈𝑈≜³_𝑖=1𝑢_∗𝑖,𝑢^∗_𝑖⊂𝑅³₊. (2.4)

Let 𝐶_𝑏([0,𝑇],𝑅⁵) denote the space of continuous bounded functions on [0,𝑇] with values in 𝑅⁵, equipped with the sup-norm topology, that is, for 𝑧∈𝐶_𝑏([0,𝑇],𝑅⁵), ‖𝑧‖_𝑐=sup{‖𝑧(𝑡)‖, 𝑡∈[0,𝑇]}, where ‖⋅‖ is the Euclidean norm.

3. Optimal Control Problem

The optimal control problem using the yield intensity of 1,3-PD at the terminal time as cost functional, based on the controlled multistage nonlinear dynamical system (2.1), can be formulated as follows:(𝑢)𝑥inf𝐽≜−³𝑢₃,𝑢𝑢₃𝑥̇(𝑡)((𝑥(0)𝑥(𝑡)s.t.=𝑓𝑡,𝑥𝑡,𝑢))=𝜉,∈𝑊,𝑡∈0,𝑢₃𝑢∈𝑈. (OCP)

From the theory on continuous dependence of solutions on parameters and our previous literature (see [16]), we know that 𝑥(⋅,𝑢) is continuous relative to 𝑢, so 𝐽(𝑢) is continuous on 𝑢∈𝑈. Moreover, 𝑈 is a closed bounded convex subset of 𝑅³₊. Hence we know the optimal control must exist, namely, ∃𝑢^∗∈𝑈 such that 𝐽(𝑢^∗)≤𝐽(𝑢), for all 𝑢∈𝑈.

3.1. Differentiability with respect to the Control Vector

In this subsection, our aim is to show the differentiability and the gradient information of solutions of the system (2.1) with respect to the control vector. To begin with, we discuss some properties of the function 𝑓(𝑡,𝑥(𝑡,𝑢)).

Proposition 3.1. For the system (2.1), 𝑓(𝑡,𝑥(𝑡,𝑢)) and (𝜕𝑓/𝜕𝑥_𝑖)(𝑡,𝑥(𝑡,𝑢))(𝑖=1,2,…,5) are continuous in (𝑡,𝑥) on an open set Δ in 𝑅₊×𝑊.

Proof. It follows by inspection that the function 𝑓∈𝐶²(𝑅₊×𝑊,𝑅⁵) by definition and (2.3).

Using Theorems I-1-4 and II-1-2 in [17], we can show that the system (2.1) has a unique solution 𝑥=𝑥(𝑡,𝑢), and the solution satisfies the integral equation𝑥(𝑡,𝑢)=𝜉+^𝑡₀𝑓(𝑠,𝑥(𝑠,𝑢))𝑑𝑠,𝑡∈𝐼. (3.1)𝑢₁=𝜉₁ and 𝑢₂=𝜉₂ are the initial value of the system (2.1), due to the differentiability of 𝑥(𝑡,𝑢) with respect to the initial vector, and we have𝜕𝑓𝜕𝜉_𝑗(𝑡,𝑥(𝑡,𝜉,𝑢))=𝑒_𝑗+^𝑡₀𝜕𝑓((𝜕𝑥𝑠,𝑥𝑡,𝜉,𝑢))𝜕𝑥𝜕𝜉_𝑗(𝑠,𝑥(𝑡,𝜉,𝑢))𝑑𝑠, (3.2) where 𝑒_𝑗 is the vector in 𝑅⁵ with entries 0 except for 1 at the 𝑗th entry, and 𝜕𝑓/𝜕𝑥 is the 5×5 matrix whose 𝑖th column is 𝜕𝑓/𝜕𝑥_𝑖, 𝑖,𝑗=1,2,…,5. From this speculation, we obtain the following result.

Proposition 3.2. Partial derivatives 𝜕𝑥/𝜕𝑢₁ and 𝜕𝑥/𝜕𝑢₂ exist and are continuous in (𝑡,𝑢). Furthermore, 𝜕𝑥/𝜕𝜉_𝑗 is the unique solution of the initial-value problem 𝑧̇=𝜕𝑓(((0)𝜕𝑥𝑡,𝑥𝑡,𝑢))𝑧,𝑧=−𝑒_𝑗. (3.3)

By virtue of the result of Proposition 3.2, we can obtain the value of 𝜕𝑥/𝜕𝑢₁ and 𝜕𝑥/𝜕𝑢₂. The following Proposition gives a formula to compute the value of 𝜕𝑥/𝜕𝑢₃.

Proposition 3.3. Partial derivatives 𝜕𝑥/𝜕𝑢₃ exist and are continuous in (𝑡,𝑢). Furthermore, 𝜕𝑥𝜕𝑢₃=𝑓𝑢₃,𝑥𝑢₃,𝑢. (3.4)

Proof. The existence and continuity of 𝜕𝑥/𝜕𝑢₃ can be directly obtained by the function 𝑓∈𝐶²(𝑅₊×𝑊,𝑅⁵) and implicit function theorem. Next, we derive the formula (3.4), for ∀𝑡∈𝐼, and we have: 𝑥𝑢₃,𝑢=𝜉+^𝑢3₀𝑓((𝑥𝑢𝑠,𝑥𝑠,𝑢))𝑑𝑠,₃+Δ𝑡,𝑢=𝜉+^𝑢3₀^+Δ𝑡𝑓(𝑠,𝑥(𝑠,𝑢))𝑑𝑠, (3.5) then, by the integral mean value theorem, there exists a constant 𝜃∈[0,1], such that 𝑥𝑢₃+Δ𝑡,𝑢−𝑥𝑢₃,𝑢=^𝑢3_𝑢^+Δ𝑡₃𝑓(𝑠,𝑥(𝑠,𝑢))𝑑𝑠=𝑓𝑢₃+𝜃Δ𝑡,𝑥𝑢₃+𝜃Δ𝑡,𝑢Δ𝑡. (3.6) Let Δ𝑡→0, then 𝜕𝑥𝜕𝑢₃=lim_Δ𝑡→0𝑥𝑢₃+Δ𝑡,𝑢−𝑥𝑢₃,𝑢Δ𝑡=lim_Δ𝑡→0𝑓𝑢₃+𝜃Δ𝑡,𝑥𝑢₃+𝜃Δ𝑡,𝑢=𝑓𝑢₃,𝑥𝑢₃,𝑢. (3.7) We obtain the desired result.

3.2. Model Transformation

The optimal control problem (3.5) is not a standard case because the terminal time 𝑡_𝑓 is free. Using the method in Section 6.8.1 of [18], the (3.5) can be transformed into the one with fixed terminal time. Treating 𝑡_𝑓 as an unknown parameter and using the transformation 𝑡=𝑡_𝑓𝜏, the (3.5) is converted to (3.8) as follows:(𝑢)𝑥inf𝐽≜−³(1,𝑢)𝑢₃s.t.𝑥̇(𝜏)=𝑡_𝑓𝑓𝑡_𝑓𝜏,𝑥𝑡_𝑓𝑥(0)𝑥(𝜏)[]𝜏,𝑢=𝜉,∈𝑊,𝜏∈0,1𝑢∈𝑈. (OCP^′)

3.3. Semi-Infinite Optimization with Inequality Constrained

For the optimal control problem 𝐽(𝑢), it is difficult to cope with the continuous state inequality constraints, that is, this is a semi-infinite optimization problem. To overcome the difficulty, let𝑔_𝑖(𝑥(𝜏,𝑢))∶=𝑥_𝑖(𝜏,𝑢)−𝑥^∗_𝑖,𝑔_5+𝑖(𝑥(𝜏,𝑢))∶=𝑥_∗𝑖−𝑥_𝑖(𝜏,𝑢),𝑖=1,2,…,5. (3.8) The condition, 𝑥(𝜏,𝑢)∈𝑊, for all 𝜏∈[0,1], is equivalently transcribed into𝐺(𝑢)=0, (3.9) where 𝐺(𝑢)=∑¹⁰_𝑙=1∫¹₀max{0,𝑔_𝑖(𝑥(𝜏,𝑢))}𝑑𝑡.

However, 𝐺(𝑢) is nonsmooth in 𝑔_𝑖(𝑥(𝜏,𝑢))=0. By the standard optimization routines [19, 20], the following smoothing technique is to replace 𝑔_𝑖,𝜖(𝑥(𝜏,𝑢)) with max{0,𝑔_𝑖(𝑥(𝜏,𝑢))}, where𝑔_𝑖,𝜖(𝑥(𝜏,𝑢))=0if𝑔_𝑖(𝑥(𝑔𝜏,𝑢))<−𝜖,_𝑖(𝑥(𝜏,𝑢))+𝜖²4𝜖if−𝜖≤𝑔_𝑖(𝑥(𝑔𝜏,𝑢))≤𝜖,_𝑖(𝑥(𝜏,𝑢))if𝑔_𝑖(𝑥(𝜏,𝑢))>𝜖. (3.10) Note that𝐺_𝜖(𝑢)=¹⁰_𝑙=1^𝑇₀𝑔_𝑙,𝜖(𝑥(𝜏,𝑢))𝑑𝑡 (3.11) is a smooth function in 𝑢. Let𝑊_𝜖≜𝑢∈𝑈∣𝐺_𝜖(𝑢)==0𝑢∈𝑈∣𝑔_𝑙(𝑥(𝜏,𝑢))≤−𝜖,𝑙=1,2,…,10,𝜏∈[0,1]. (3.12) Clearly, 𝑊_𝜖⊆𝑊∩𝑈 for each 𝜖>0.

We now define an approximate problem denoted by (3.13), where the smoothed state constraint is treated as the penalty function:inf𝐽_𝜖,𝛾(𝑢)∶=−𝑥₃(1,𝑢)+𝛾𝐺_𝜖(𝑢)s.t.𝑥̇(𝑡)=𝑡_𝑓𝑓𝑡_𝑓𝜏,𝑥𝑡_𝑓,𝑥(0)𝜏,𝑢=𝜉,𝑢∈𝑊_𝜖. (OCP_𝜖,𝛾)

By similar arguments as those given in [21], (3.13) is coincident with (3.5) as 𝜖→0. On this basis, (3.5) can be explored by solving a sequence of approximate (3.13). Each of these (3.13) is viewed as a smooth nonlinear mathematical programming problem.

3.4. Optimization Algorithm

In this subsection, similar with the approach based on gradient in [22, 23], we proposed an algorithm based on gradient of 𝐽_𝜖,𝛾(𝑢) to solve the (3.13). 𝜕𝐽_𝜖,𝛾(𝑢)/𝜕𝑢≜𝐽_𝑔(𝑢) can be derived by solving the ordinary differential equation (3.3) and computing (3.4). The admissible control set 𝑈 can be called a “box” because of its rectangular shape, we use the classic gradient projection method to cope with “box.”

The projection of an arbitrary u onto the feasible set 𝑈 is defined as follows. The 𝑖th component is given by𝑝𝑢,𝑢_∗,𝑢^∗_𝑖=𝑢_𝑖∗if𝑢_𝑖<𝑢_𝑖∗,𝑢_𝑖if𝑢_𝑖∈𝑢_𝑖∗,𝑢^∗_𝑖,𝑢^∗_𝑖if𝑢_𝑖>𝑢^∗_𝑖. (3.13)

Thus, the control variable 𝑢(𝑡) obtained by projecting the steepest descent direction at 𝑢 onto the feasible set 𝑈 is given by𝑢(𝜆)=𝑝𝑢−𝜆𝐽_𝑔,𝑢_∗,𝑢^∗, (3.14) where 𝜆 is optional step size.

On the basis of the above analysis, we can obtain an optimal control for (3.13) as shown in the following algorithm.

Algorithm 3.4. Step 1. Set constants 𝛼,𝛿∈(0,1), and 𝑟_max is positive constant. Set 𝑟=0, compute 𝐽_𝑔(𝑢_𝑟) by equations (3.3) and (3.4), if 𝐽_𝑔(𝑢_𝑟)<𝛿, stop. Else, then go to Step 2.Step 2. Compute the step-size 𝜆(𝑟)=𝛽^𝑘𝑟 with Armijo line search rules, where 𝑘_𝑟 is any integer such that 𝐽_𝜖,𝛾𝑢(𝑟)+𝛽^𝑘𝑟𝐽_𝑔𝑢_𝑟−𝐽_𝜖,𝛾(𝑢(𝑟))≤−𝛽^𝑘𝑟𝛼𝐽^𝑇_𝑔𝑢_𝑟𝐽_𝑔𝑢_𝑟,𝐽_𝜖,𝛾𝑢(𝑟)+𝛽^𝑘𝑟𝐽_𝑔𝑢_𝑟−𝐽_𝜖,𝛾(𝑢(𝑟))>−𝛽^𝑘𝑟−1𝛼𝐽^𝑇_𝑔𝑢_𝑟𝐽_𝑔𝑢_𝑟, (3.15)(3.16) go to Step 3.Step 3. If 𝑟>𝑟_max, stop. Otherwise, compute 𝑢(𝑟+1)=𝑝(𝑢(𝑟)−𝜆(𝑟)𝐽_𝑔(𝑢_𝑟),𝑢_∗,𝑢^∗) using the equation defined by (3.14), replace 𝑟 by 𝑟+1, and go to Step 1.
Note that, due to the boundedness of the function 𝐽_𝜖,𝛾(⋅), it is very easy to find a 𝑘_𝑟 satisfying (3.15) and (3.16), using the following subprocedure, which uses the last used step length 𝜆_𝑟−1=𝛽^𝑘𝑟−1, as the starting point for the computation of the next one.

Subprocedure of Algorithm 3.4
Step 1. If 𝑟=0, set 𝑘^′=0. Else, set 𝑘^′=𝑘_𝑟−1.Step 2. If 𝑘_𝑟=𝑘^′ satisfies (3.15) and (3.16). Else, set 𝑘^′=𝑘_𝑟−1, stop.Step 3. If 𝑘_𝑟=𝑘^′ satisfies (3.15) and not (3.16), replace 𝑘^′ by 𝑘^′−1, and go to Step 2. If 𝑘_𝑟=𝑘^′ satisfies (3.16) and not (3.15), replace 𝑘^′ by 𝑘^′+1, and go to Step 2.

For Algorithm 3.4, we see that −𝐽^𝑇_𝑔𝐽_𝑔(⋅) is continuous, that −𝐽^𝑇_𝑔𝐽_𝑔(𝑢)≤0 for all 𝑢∈𝑈, and that −𝐽²_𝑔(𝑢)=0 if and only if 𝐽_𝑔(𝑢)=0, that is, that −𝐽²_𝑔(⋅) is an optimality function for the problem (3.13). So, we have the following theorem to guarantee the convergence of the algorithm.

Theorem 3.5. If 𝑢_𝑟 is such that 𝐽_𝑔≠0, then 𝜆_𝑟 is computed by Algorithm 3.4 using a finite number of function evaluations and any accumulation point 𝑢 of this sequence satisfies 𝐽_𝑔(𝑢)=0.

Proof. We apply Theorem 1.2.24𝑎 in [24] with 𝜃(⋅)=−𝐽^𝑇_𝑔𝐽_𝑔(⋅), then the desired result can be obtained immediately.

Remark 3.6. Although the approach we are using here to deal with the constraints of continuous state is similar with the one mentioned in [19], There are still three main difference between them. First of all, our approach is applied for the batch culture, while the one in [19] is used for the fed-batch culture. Secondly, the control variables in the two approaches are different. Our variable controls the initial and terminal points, and their variable controls the switching time. Finally, we use the gradient-based algorithm to numerically solve the problem. Their algorithm is an improved Particle Swarm Optimization (PSO) algorithm, not gradient based.

4. Numerical Results

According to the model and algorithm mentioned above, we have programmed the software and applied it to the optimal control problem of microbial fermentation in batch culture. The system parameters are listed in Table 1 (see [9, 16]).

Table 1 

Parameters values in dynamical system (2.1).

The basic data are listed, respectively, as follows.

Boundary Value of Control Vector
𝑢_∗1=0.01 mmol/L, 𝑢^∗₁=1 mmol/L, 𝑢_∗2=200 mmol/L, 𝑢^∗₂=939.5 mmol/L, 𝑢_∗3=2ℎ, and 𝑢^∗₃=10 h.

Boundary Value of State Vector
𝑥_∗1=0.001 mmol/L, 𝑥^∗₁=2039 mmol/L, 𝑥_∗2=0.001 mmol/L, 𝑥^∗₂=939.5 mmol/L, 𝑥_∗3=0.01 mmol/L, 𝑥^∗₃=10 mmol/L. 𝑢_∗4=0.01 mmol/L, 𝑢^∗₄=1026 mmol/L, and 𝑢_∗5=200, 𝑢^∗₅=360.9 mmol/L.

We adopt 𝛼=0.4,𝛿=0.00001, and 𝑟_max=1000 in the procedure. Then, by Algorithm 3.4, the optimal control vector 𝑢 and objective function 𝐽_𝜖,𝛾(𝑢) are (0.973186,547.04,5.17509)^𝑇 and 54.5911, respectively. Numerical results show that, by employing the optimal control, the concentration of 1,3-PD at the terminal time can be increased, compared with the previous results.

5. Conclusions

In this paper, based on the previous model in [16], taking the yield intensity of 1,3-PD as the performance index and the initial concentration of biomass, glycerol, and terminal time as the control vector, we propose an optimal control model subject to a multistage nonlinear dynamical system and constraints of continuous state. A computational approach is constructed to seek the solution of the above model in two aspects. The convergence analysis and optimality function of the algorithm are also investigated. Numerical results show that, by employing the optimal control, yield intensity of 1,3-PD at the terminal time can be increased, compared with the previous results.

Our current tasks accommodate the modeling and simulation of the fermentation process. Moreover, the stability and reachability of the improved model need to be discussed.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (no. 10871033), the Natural Science Foundation for the Youth of China (no. 11001153) the Shandong Province Natural Science Foundation of China (no. ZR2010AQ016) and the Fundamental Research Funds for the Central Universities.