# Optimal Control of a Nonlinear Time-Delay System in Batch Fermentation Process

Authors: Yongsheng Yu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The main control goal in batch process is to get a high yield of products. In this
paper, to maximize the yield of 1,3-propanediol (1,3-PD) in bioconversion of glycerol
to 1,3-PD, we consider an optimal control problem involving a nonlinear
time-delay system. The control variables in this problem include the initial concentrations
of biomass and glycerol and the terminal time of the batch process.
By a time-scaling transformation, we transcribe the optimal control problem into
a new one with fixed terminal time, which yields a new nonlinear system with
variable time-delay. The gradients of the cost and constraint functionals with
respect to the control variables are derived using the costate method. Then, a
gradient-based optimization method is developed to solve the optimal control
problem. Numerical results show that the yield of 1,3-PD at the terminal time
is increased considerably compared with the experimental data.


---

## Body

## 1. Introduction

Batch processing represents the natural way to scale-up processes from the laboratory to the production environment [1]. During the batch operation, no substrate is added to the initial charge and no product is removed until the end of the process [2]. Optimal control of batch processes has received attention recently because it is a choice for reducing production costs, improving product quality, and meeting safety requirements [3–5].1,3-Propanediol (1,3-PD) is one of the important products used in the chemical industry. Using glycerol for producing 1,3-PD is effective from both the economical and ecological point of view. The operations of glycerol bioconversion to 1,3-PD consist of batch, continuous, and fed-batch cultures. Compared with continuous and fed-batch cultures, glycerol fermentation in batch process can obtain the highest production concentration and molar yield 1,3-PD to glycerol [6]. For this process, many studies have been carried out including the quantitative description of the cell growth kinetics [7, 8], enzyme-catalytic kinetics [9], the parameter identification problem [10–13], and the pathway optimization problem [14]. In particular, considering the existence of time-delays in the fermentation process [15, 16], a nonlinear time-delay system was recently proposed in [13]. Numerical simulations showed that the nonlinear time-delay system can describe the batch process better than previous mathematical models. However, in batch process of glycerol bioconversion to 1,3-PD, the aim is to obtain as much 1,3-PD as possible at the terminal time. Mathematically, this is an optimal control problem. Although the achieved results are interesting, such optimal control problem is ignored in the above researches.In this paper, taking the initial concentrations of biomass and glycerol and the terminal time of the fermentation process as control variables and the yield of 1,3-PD as the performance index, we propose an optimal control model involving nonlinear time-delay system in [13] and subject to continuous state constraints. In fact, this problem is a free time delayed optimal control problem. By the way, there has been a mounting interest in optimal control of time-delay systems. A maximum principle for optimal control problem with a constant delay was provided in [17]. It should, however, be noted that many time-delay systems including the nonlinear time-delay system in this paper are highly nonlinear. Therefore, it is often impossible to obtain analytical solutions of the delayed optimal control problem and one has to resort to numerical solution methods. As a result, some successful algorithms, such as iterative dynamic programming method [18], control parameterization method [19], and measure theoretical approach [20], have been developed. However, terminal time in above optimal control problems is fixed. The presence of free terminal time makes the delayed optimal control problem much more complicated. Thus, to solve the optimal control problem in this work, these existing computational methods cannot be used directly and new computational techniques should be explored.In this paper, by a time-scaling transformation, we equivalently transcribe the free time delayed optimal control problem into one with fixed terminal time. It is significant to mention that this transformation yields a more complex dynamic system in which the time-delays are variable. By the costate method, we derived the gradients of the cost functional and constraints with respect to the control variables. Then, a gradient-based optimization technique is developed. Finally, numerical results show that the yield of 1,3-PD at the terminal time is increased considerably compared with the experiment data.


## 2. Nonlinear Time-Delay Systems

In batch process, a proper quantity of biomass and glycerol is added to the reactor only once and stirred uniformly under given conditions. Then, 1,3-PD is removed at the end of the process. Based on the work [13], the following nonlinear time-delay system can be used to describe the batch process:
(1)x˙(t)=f(x(t),x(t-h)),t∈(0,T],x(0)=ζ,x(t)=ϕ(t,ζ),t≤0,
where x(t):=(x1(t),x2(t),x3(t),x4(t),x5(t))⊤∈R5 is the state vector whose components are, respectively, the concentrations of biomass, glycerol, 1,3-PD, acetate, and ethanol in the reactor at time t; h is a delay argument; T>0 is the terminal time; ζ:=(ζ1,ζ2,ζ3,ζ4,ζ5)⊤ is the initial state; ϕ:R×R5→R5 is a given initial function; and the dynamics of the batch process is given by
(2)f(x(t),x(t-h)):=(μx1(t-h)-q2x1(t-h)q3x1(t-h)q4x1(t-h)q5x1(t-h)).
In (2), the specific growth rate of cells μ, the specific consumption rate of substrate q2, and the specific formation rates of products ql, l=3,4,5, can be expressed as
(3)μ=μmx2(t)k1+x2(t)∏l=25‍(1-xl(t)xl*),ql=ml+μYl,l=2,3,…,5,
where μm, k1, ml, and Yl are kinetic parameters and x* and x* are the critical concentrations for cells growth.The terminal time in system (1) is a control variable. Define
(4)T:={γ∈R:Tmin⁡≤γ≤Tmax⁡},
where Tmin⁡ and Tmax⁡ are the lower and upper bounds for the terminal time, respectively. Any T∈T is called anadmissible terminal time. Moreover, since no product comes into being at the initial point, ζl=0,l=3,4,5. Hence, let u=(ζ1,ζ2)⊤ and define
(5)U:={υ∈R2:ai≤υi≤bi},
where ai and bi are real numbers such that ai<bi. Any u∈U is called anadmissible initial vector of biomass and glycerol.For the nonlinear time-delay system (1), there exists a unique continuous solution x(·∣u) corresponding to each u∈U on [0,∞) [21]. Furthermore, there exist critical concentrations of biomass, glycerol, and products, outside which cells cease to grow. Hence, it is biologically meaningful to restrict the concentrations of biomass, glycerol, and products within a set W defined as
(6)x(t∣u)∈W:=∏l=15‍[x*l,xl*],t∈[0,T].


## 3. Optimal Control Problems

In the batch process, it is desired that the yield of 1,3-PD should be maximized at the end of the process. This is achieved by manipulating some control variablesu and T. Thus, the optimal control problem in batch process can be formulated as
(7)(OCP)min⁡J(u,T)=-x3(T∣u)Ts.t.x(t∣u)∈W,t∈[0,T],(u,T)∈U×T.Note that (OCP) is of nonstandard type because the terminal time and the initial values of biomass and glycerol are variables to be determined. Thus, the (OCP) is actually a free time optimal control problem. It is difficult to solve the (OCP) using existing numerical techniques [18–20]. The main difficulty is the implicit dependence of the system state on the terminal time. We now employ a time-scaling transformation from [0,T] to [0,1] as follows:
(8)t=t(s)=Ts,
where s∈[0,1] is a new time variable. Clearly, s=0 corresponds to t=0, and s=1 corresponds to t=T. Let x~(s)=x(t(s)). Then
(9)x~˙(s)=dds{x(t(s))}=dx(t(s))dtdt(s)ds=Tf(x(Ts),x(Ts-h))=Tf(x~(s),x~(s-hT-1)).
The initial conditions become
(10)x~(0)=ζ,x~(s)=ϕ(Ts,ζ),s≤0.
Thus, system (1) is transformed into the following form:
(11)x~˙(s)=Tf(x~(s),x~(s-hT-1)),s∈[0,1],x~(0)=ζ,x~(s)=ϕ(Ts,ζ),s≤0.
Let x~(·∣u,T) denote the solution of system (11) corresponding to each (u,T)∈U×T on [0,1]. Then, (OCP) can be transcribed into the following equivalently optimal control problem with fixed terminal time:
(12)(EOCP)min⁡J(u,T)=-x~3(1∣u,T)T(13)s.t.x~(s∣u,T)∈W,s∈[0,1],(u,T)∈U×T.Note that the time-delay in system (11) depends on the control variable T. Thus, when applied to time-delay systems, the time-scaling transformation (8) yields a more complex dynamic system in which the time-delay is variable.


## 4. Computational Approaches

(EOCP) is in essence an optimization problem with continuous state inequality constraint (13), which has an infinite number of constraints and can be viewed as a semi-infinite programming problem. An efficient algorithm transforming this type of problems to standard mathematical programming problems was discussed in [22]. We will now briefly discuss the application of this algorithm to (EOCP).Let(14)gl(x~(s∣u,T)):=xl*-x~l(s∣u,T),g5+l(x~(s∣u,T)):=x~l(s∣u,T)-x*l,l=1,2,…,5.
The continuous state constraint (13) becomes
(15)G(u,T)=0,
where G(u,T):=∑l=110∫01‍min⁡{0,gl(x~(s∣u,T))}ds. However, the equality constraint (15) is nondifferentiable at the points when gl=0. We replace (15) with the following inequality constraint:
(16)G~ϵ,τ(u,T):=τ+∑l=110‍∫01‍πϵ(gl(x~(s∣u,T)))ds≥0,
where ϵ>0, τ>0, and
(17)πϵ(η)={η,ifη<-ϵ,-(η-ϵ)24ϵ,if-ϵ≤η≤ϵ,0,ifη>ϵ.
Thus, (EOCP) is approximated by a sequence of standard mathematical programming problems {(EOCPϵ,τ)} defined by replacing constraint (15) with (16). Moreover, each of {(EOCPϵ,τ)} can be solved by a gradient-based optimization method (e.g., sequential quadratic programming (SQP) [23]). However, this optimization method requires the gradients of the cost functional and constraints. Now, we will derive these required gradients using the so-calledcostate method, which is a commonly used technique in the optimal control domain [22, 24, 25].Define(18)ψ(s):={∂ϕ(Ts,ζ)∂s,s≤0,Tf(x~(s),x~(s-hT-1)),s∈[0,1],χI(s):={1,s∈I,0,otherwise.
The gradients of the cost functional J(·,·) defined in (12) with respect to the control variables are given in the following theorem.Theorem 1.
Let(u,T)∈U×T. Then
(19)∂J(u,T)∂u=λ⊤(0)∂x~(0)∂u+∫-hT-10‍Tλ⊤(s+hT-1)×∂f(x~(s+hT-1),x~(s))∂x~(s)∂ϕ(Ts,ζ)∂uds,(20)∂J(u,T)∂T=x~3(1)T-2+∫01‍λ⊤(s)f(x~(s),x~(s-hT-1))ds+∫01‍hT-1λ⊤(s)×∂f(x~(s),x~(s-hT-1))∂x~(s-hT-1)ψ(s-hT-1)ds+∫-hT-10‍Tλ⊤(s+hT-1)×∂f(x~(s+hT-1),x~(s))∂x~(s)∂ϕ(Ts,ζ)∂Tds,
where λ(·) is the solution of the following costate system:(21)λ˙(s)=-T(∂f(x~(s),x~(s-hT-1))∂x~(s))⊤λ(s)-Tχ[0,1-hT-1](s)(∂f(x~(s+hT-1),x~(s))∂x~(s))⊤λ(s+hT-1),hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhs∈[0,1],with the initial conditions
(22)λ(1)=(0,0,-T-1,0,0)⊤,λ(s)=(0,0,0,0,0)⊤,s>1.Proof.
The derivations of the gradients of the cost functionalJ(u,T) with respect to u and T are similar. Thus, only the derivation of the gradient of the cost functional J(u,T) with respect to u is given below.
Letw:[0,∞]→R5 be an arbitrary function that is continuous and differentiable almost everywhere. Then, we may express the cost functional J as follows:
(23)J(u,T)=-T-1x~3(1)+∫01‍w⊤(s)(Tf(x~(s),x~(s-hT-1))-x~˙(s))ds=-T-1x~3(1)+∫01‍Tw⊤(s)f(x~(s),x~(s-hT-1))ds-∫01‍w⊤(s)x~˙(s)ds.
Applying integration by parts to the last integral term gives
(24)J(u,T)=-T-1x~3(1)+∫01‍(Tw⊤(s)f(x~(s),x~(s-hT-1))+w˙⊤(s)x~(s))ds-w⊤(1)x~(1)+w⊤(0)x~(0).
Differentiating (24) with respect to u yields(25)∂J(u,T)∂u=-T-1∂x~3(1)∂u-w⊤(1)∂x~(1)∂u+w⊤(0)∂x~(0)∂u+∫01‍(Tw⊤(s)∂f(x~(s),x~(s-hT-1))∂x~(s)+w˙⊤(s))∂x~(s)∂uds+∫01‍Tw⊤(s)∂f(x~(s),x~(s-hT-1))∂x~(s-hT-1)∂x~(s-hT-1)∂uds.Performing a change of variable in the last term on the right-hand side of (25) gives
(26)∫01‍Tw⊤(s)∂f(x~(s),x~(s-hT-1))∂x~(s-hT-1)∂x~(s-hT-1)∂uds=∫01‍Tχ[0,1-hT-1](s)w⊤(s+hT-1)×∂f(x~(s+hT-1),x~(s))∂x~(s)∂x~(s)∂uds+∫-hT-10‍Tw⊤(s+hT-1)×∂f(x~(s+hT-1),x~(s))∂x~(s)∂ϕ(Ts,ζ)∂uds.
Substituting (26) into (25) yields
(27)∂J(u,T)∂u=-T-1∂x~3(1)∂u-w⊤(1)∂x~(1)∂u+w⊤(0)∂x~(0)∂u+∫01‍(w˙⊤(s)+Tw⊤(s)∂f(x~(s),x~(s-hT-1))∂x~(s)+Tχ[0,1-hT-1](s)w⊤(s+hT-1)×∂f(x~(s+hT-1),x~(s))∂x~(s))∂x~(s)∂uds+∫-hT-10‍Tw⊤(s+hT-1)×∂f(x~(s+hT-1),x~(s))∂x~(s)∂ϕ(Ts,ζ)∂uds.
Choosing w(·)=λ(·) and substituting (21)-(22) into the above equation, we obtain the conclusion (19). The gradient formula (20) can be derived similarly. The proof is complete.The gradients of the constraintG~ϵ,τ(·,·) defined in (16) with respect to the control variables are given in the next theorem.Theorem 2.
Let(u,T)∈U×T. Then
(28)∂G~ϵ,τ(u,T)∂u=λ~⊤(0)∂x~(0)∂u+∫-hT-10‍Tλ~⊤(s+hT-1)×∂f(x~(s+hT-1),x~(s))∂x~(s)∂ϕ(Ts,ζ)∂uds,(29)∂G~ϵ,τ(u,T)∂T=∫01‍λ~⊤(s)f(x~(s),x~(s-hT-1))ds+∫01‍hT-1λ~⊤(s)∂f(x~(s),x~(s-hT-1))∂x~(s-hT-1)×ψ(s-hT-1)ds+∫-hT-10‍Tλ~⊤(s+hT-1)×∂f(x~(s+hT-1),x~(s))∂x~(s)∂ϕ(Ts,ζ)∂Tds,
where λ~(·) is the solution of the following costate system:
(30)λ~˙(s)=-∑l=110‍∂πϵ(gl(x~(s)))∂x~(s)-T(∂f(x~(s),x~(s-hT-1))∂x~(s))⊤λ~(s)-Tχ[0,1-hT-1](s)(∂f(x~(s+hT-1),x~(s))∂x~(s))⊤×λ~(s+hT-1),hhhhhs∈[0,1],
with the initial condition
(31)λ~(s)=(0,0,0,0,0)⊤,s≥1.Proof.
The derivations of the gradients of the constraintG~ϵ,τ(u,T) with respect to u and T are similar. Thus, only the derivation of the gradient of the constraint G~ϵ,τ(u,T) with respect to u is given below.
Letw:[0,∞]→R5 be an arbitrary function that is continuous and differentiable almost everywhere. Then, we may express the constraint G~ϵ,τ as follows:
(32)G~ϵ,τ(u,T)=τ+∑l=110‍∫01‍πϵ(gl(x~(s)))ds+∫01‍Tw⊤(s)f(x~(s),x~(s-hT-1))ds-∫01‍w⊤(s)x~˙(s)ds.
Applying integration by parts to the last integral term gives
(33)G~ϵ,τ(u,T)=τ-w⊤(1)x~(1)+w⊤(0)x~(0)+∑l=110‍∫01‍πϵ(gl(x~(s)))ds+∫01‍(Tw⊤(s)f(x~(s),x~(s-hT-1))+w˙⊤(s)x~(s)Tw⊤(s)f(x~(s),x~(s-hT-1)))ds.
Differentiating (33) with respect to u yields(34)∂G~ϵ,τ(u,T)∂u=-w⊤(1)∂x~(1)∂u+w⊤(0)∂x~(0)∂u+∫01‍(∑l=110‍∂πϵ(gl(x~(s)))∂x~(s)+Tw⊤(s)∂f(x~(s),x~(s-hT-1))∂x~(s)+w˙⊤(s)∑l=110‍∂πϵ(gl(x~(s)))∂x~(s))∂x~(s)∂uds+∫01‍Tw⊤(s)∂f(x~(s),x~(s-hT-1))∂x~(s-hT-1)∂x~(s-hT-1)∂uds.Performing a change of variable in the last term on the right-hand side of (34) gives
(35)∫01‍Tw⊤(s)∂f(x~(s),x~(s-hT-1))∂x~(s-hT-1)∂x~(s-hT-1)∂uds=∫01‍Tχ[0,1-hT-1](s)w⊤(s+hT-1)×∂f(x~(s+hT-1),x~(s))∂x~(s)∂x~(s)∂uds+∫-hT-10‍Tw⊤(s+hT-1)×∂f(x~(s+hT-1),x~(s))∂x~(s)∂ϕ(Ts,ζ)∂uds.
Substituting (35) into (34) yields
(36)∂G~ϵ,τ(u,T)∂u=-w⊤(1)∂x~(1)∂u+w⊤(0)∂x~(0)∂u+∫01‍(∑l=110‍∂πϵ(gl(x~(s)))∂x~(s)+w˙⊤(s)+Tw⊤(s)∂f(x~(s),x~(s-hT-1))∂x~(s)+Tχ[0,1-hT-1](s)w⊤(s+hT-1)×∂f(x~(s+hT-1),x~(s))∂x~(s))∂x~(s)∂uds+∫-hT-10‍Tw⊤(s+hT-1)×∂f(x~(s+hT-1),x~(s))∂x~(s)∂ϕ(Ts,ζ)∂uds.
Choosing w(·)=λ~(·) and substituting (30) and (31) into the above equation, we obtain the conclusion (28). The gradient formula (29) can be derived similarly. Thus, the proof is complete.On the basis of Theorems1 and 2, we can develop the following algorithm to solve the (OCP).Algorithm
1.Step 1.
Choose initial values ofu, T, ϵ, and τ.Step 2.
Solve (EOCPϵ,τ) using SQP [23] to give (uϵ,τ*,Tϵ,τ*).Step 3.
Check feasibility ofgl(x~(s∣uϵ*,Tϵ*))≥0, for l=1,2,…,10.Step 4.
If(uϵ,τ*,Tϵ,τ*) is feasible, then go to Step 5. Otherwise, set τ:=ατ, where α is a given positive constant. If τ<τ-, where τ- is a prespecified positive constant, then go to Step 6. Otherwise, go to Step 2.Step 5.
Setϵ:=βϵ, where β is a given positive constant. If ϵ>ϵ-, then go to Step 2. Otherwise, go to Step 6.Step 6.
Output(uϵ,τ*,Tϵ,τ*) and stop.At the conclusion of Steps1–6, (uϵ,τ*,Tϵ,τ*) is an approximate optimal solution of the (OCP).


## 5. Numerical Results

Algorithm  1 was applied to seek the optimal control variables in (OCP) and all computations were implemented in Fortran 6.5. We obtained thatu*=(0.201867,418.376)⊤ and T*=7.2578 h. Here, nonlinear time-delay system (1) was solved using the six-order Runge-Kutta method with Lagrange interpolation [26]. The time-delay h=0.26 h, the initial state ζ=(0.102,418.26087,0,0,0)⊤, the initial function ϕ(t)=ζ, the terminal time T=6.92 h, and the kinetic parameters and critical concentrations are listed in Table 1. In the optimization process, we assume that u1∈[0.02,0.202], u2∈[400,500], and T∈[5.0,8.0]. In addition, the smoothing and feasible parameters were initially selected as ϵ=0.1 and τ=0.01. The parameters α and β were chosen as 0.1 and 0.01 until the solution obtained was feasible for the (OCP). The process terminated when ϵ-=1.0×10-8 and τ-=1.0×10-7.Table 1
The kinetic parameters and critical concentrations in the system (1) [13].
μ
m
k
1
m
2
m
3
m
4
m
5
Y
2
Y
3
Y
4
Y
5
0.994
0.368
3.24
−3.679
−0.491
7.309
117.647
76
35.54
14.78
x
*
1
x
*
2
x
*
3
x
*
4
x
*
5
x
1
*
x
2
*
x
3
*
x
4
*
x
5
*
0.01
0
0
0
0
6
2039
939.5
1026
360.9Under the obtained optimal control strategy, we obtained the optimal yield of 1,3-PD at the optimal terminal time is78.1586 mmol h−1, which is increased by 90.04% compared with experiment data [13]. Furthermore, we plotted the optimal yield change of 1,3-PD with respect to the fermentation time in Figure 1. From Figure 1, we can see that the yield of 1,3-PD at the terminal time is actually increased.Figure 1
The yield change of 1,3-PD with respect to fermentation time.


## 6. Conclusion

In this paper, the optimal control problem in batch process was investigated. We presented the optimal control problem involving a nonlinear time-delay system and with free terminal time. We then transcribe the free time optimal control problem into a new one with fixed terminal time and variable time-delay. We developed a computational method based on the gradients of the cost and constraint functionals with respect to the control variables. Numerical results showed the effectiveness of the developed computational method.



---

*Source: 478081-2014-07-01.xml*