# A Multistage Feedback Control Strategy for Producing 1,3-Propanediol in Microbial Continuous Fermentation

Authors: Honghan Bei; Lei Wang; Jing Sun; Liwei Zhang
Journal: Complexity
 (2019)

---

## Abstract

In this paper, we consider a multistage feedback control strategy for producing 1,3-propanediol in microbial continuous fermentation. Both the dilution rate and the concentration of glycerol in the input feed are used as control variables, and these variables are further assumed to be in the form of a linear combination of biomass and glycerol concentrations. Unlike the general form of linear feedback control, the coefficients of linear combination are continuous functions with respect to time. Inspired by the control parameterization method, we use the piecewise-constant functions to approximate the coefficient functions; then we get the multistage feedback control law by solving nonlinear mathematical programming problems. Numerical results indicate the flexibility and effectiveness of our strategy.


---

## Body

## 1. Introduction

Nowadays, the production of 1,3-propanediol(1,3-PD) by microbial fermentation is very popular. There are three main ways to produce 1,3-PD by microbial fermentation: continuous culture, batch culture, and fed-batch culture. The main method for producing 1,3-PD by continuous culture is widely concerned [1], since it has the advantages of high production strength, stable production, and high automation. Compared with chemical synthesis, microbial fermentation for producing 1,3-PD is more attractive because it is easy to implement and does not generate toxic byproducts. Unfortunately, the production of 1,3-PD by microbial fermentation is not up to the standard of industry. Therefore, more and more scholars focus on optimizing the yield of 1,3-PD via microbial fermentation [2–7].Previous studies on the optimal control of 1,3-PD only achieved the open-loop control, see, for example, [8–10]. This is a discount in practical production. Based on our previous study, we find that feedback control is more in line with actual production process or experimental process by realizing closed-loop control. The linear feedback optimal control [11] has been widely studied in theory and applications. Different from the traditional approach to determine an optimal feedback control such as solving the well-known Hamilton-Jacobi-Bellman partial differential equation, the sensitivity penalization approach for computing robust suboptimal controllers [12], and the neighboring extremal approach [13, 14], we consider the feedback control coefficient function with time dependence. In order to further increase the yield of 1,3-PD, we regard both dilution rate and the concentration of glycerol in the input feed as the controllers which are further assumed to be in the form of a linear combination of biomass and glycerol concentrations. Unlike the general form of linear feedback control, the coefficients of linear combination are continuous functions with respect to time. Inspired by the control parameterization method, we use the piecewise-constant functions to approximate the coefficient functions; then we get the multistage feedback control law by solving nonlinear mathematical programming problems. Finally, numerical results show the flexibility and effectiveness of our strategy.


## 2. Problem Statement

The dynamic model of the continuous culture process is based on the following assumptions.Assumption 1.
The material composition in the fermentation tank does not change with the position of space, and the solution in the reactor is sufficiently well mixed so that the concentrations of reactants are uniform.Assumption 2.
The continuously added medium contains glycerin only, and the substance in the reactor is exported at dilution rateD(t).Assumption 3.
The materials in the fermentation are fully mixed in which the concentrations are even, and the concentrations of the reactants change only with the change of reaction time.Under the above assumptions, the mass balance relationships for biomass, substrate, and products in the microbial continuous culture can be expressed as the following nonlinear dynamic system:(1)x1˙t=f1t,ut=μ-Dtx1t,x2˙t=f2t,ut=DtCs0t-x2t-q2x1t,x3˙t=f3t,ut=q3x1t-Dtx3t,x4˙t=f4t,ut=q4x1t-Dtx4t,x5˙t=f5t,ut=q5x1t-Dtx5t,and(2)xi0=x0ii=1,2,3,4,5,where x1(t), x2(t), x3(t), x4(t), and x5(t), respectively represent the concentrations of biomass, extracellular glycerol, extracellular 1,3-PD, acetate, and ethanol at time t; x0i are the initial concentrations of biomass, glycerol, 1,3-PD, acetate, and ethanol; t∈[0,tf], tf is the terminal time; μ is the specific growth rate of cells; q2 is the specific consumption rate of substrate; qi,i=3,4,5, are the specific formation rates of 1,3-PD, acetate, and ethanol, respectively; D(t) denotes the dilution rate; Cs0(t) denotes the concentration of glycerol in the input feed; and u(t)=(D(t),Cs0(t)). In particular,(3)μ=μmx2tx2t+ks∏i=251-xitxi∗,(4)q2=m2+μY2+Δq2x2tx2t+k2,(5)q3=m3+μY3+Δq3x2tx2t+k3,(6)q4=m4+μY4+Δq4x2tx2t+k4,(7)q5=q2b1c1+Dtx2t+b2c2+Dtx2t,where μm=0.67 is the maximum specific growth rate; ks=0.28 is the Monod saturation constant for substrate. Under anaerobic conditions at 37°C and pH=7.0, the values of other parameters used in (1) - (7) are listed in Table 1.Table 1
The values of parameters [15].
i m i Y i Δ q i k i bi c i 1 - - - - 0.025 0.06 2 2.20 0.0082 28.58 11.43 5.18 50.45 3 -2.69 67.69 26.59 15.50 - - 4 -0.97 33.07 5.74 85.71 - -Letx(t)=(x1(t),x2(t),x3(t),x4(t),x5(t))T, x0=(x01,x02,…,x05)T, and f(x(t),u(t))≔(f1(t,u(t)),…,f5(t,u(t)))T. Thus the nonlinear control system can be formulated as follows:(8)x˙t=fxt,ut,t∈0,tfx0=x0.For the actual bioprocess, it should be noted that there exist critical concentrations for the state vector. Therefore, it is natural to restrict the concentrations of biomass, glycerol, and products in a given setW defined as (9)xt∈W≔x∗,x∗=∏i=15xi∗,xi∗⊂R+5where xi∗ and xi∗, respectively, denote the upper and lower bounds of the corresponding state variables.Equation (9) can be equivalently transformed ascontinuous state inequality constraints by introducing the functions as follows:(10)hixt≤0,t∈0,tf,i=1,…,10,where hi(x(t))=xi(t)-xi∗ and hi+5(x(t))=x∗i-xi(t), i=1,…,5.In this paper, the dilution rateD(t) and the concentration of glycerol in the input feed Cs0(t) are chosen as the control variables u(t). It is obvious that the control variables are also constrained: (11)u∗≤ut≤u∗,t∈0,tf,where u∗ and u∗ are the lower and upper bounds of u(t).Letx(·∣u) be the solution of (8) on [0,tf]. We can describe the optimal control problem as follows.Problem P 0. Choose u to minimize the cost function(12)minJ0u=-x3tf∣us.t.x˙t=fxt,ut,xt∈W,u∗≤ut≤u∗,t∈0.tf.


## 3. Feedback Optimal Control Problem

In microbial fermentation, the most important factors to influence the final concentration of 1,3-PD are the concentrations of biomass and glycerol. And the linear state feedback control law is one of the most common feedback control structures [16]. Thus, the feedback controller is considered as a linear combination form of the biomass and the glycerol concentrations, i.e.,(13)ut=Dt,Cs0t)=φ1xt,ξt,φ2xt,ξt,t∈0,tf,where ξ(t)=(ξ1(t),ξ2(t),ξ3(t),ξ4(t))T; φ1(x(t),ξ(t))=ξ1(t)x1(t)+ξ2(t)x2(t); and φ2(x(t),ξ(t))=ξ3(t)x1(t)+ξ4(t)x2(t).The following bound constraints are imposed on thefeedback control coefficientsξ(t):(14)Uad=α1,β1×α2,β2×α3,β3×α4,β4,t∈0,tf.Substituting (13) into (8) gives(15)xi˙t=fi~xt,ξt,t∈0,tf,i=1,2,…,5,where(16)f~xt,ξt=f1~xt,ξt=μx1t-ξ1tx12t-ξ2tx1tx2t,f2~xt,ξt=ξ1tξ3tx12t+ξ1tξ4t-ξ1t+ξ2tξ3tx1tx2t+ξ2tξ4t-ξ2tx22t-q2x2t,fi~xt,ξt=qix1t-ξ1tx1txit-ξ2tx2txit,i=3,4,5.Consider system (16) with the initial condition of (8). Let x(·∣ξ) denote the solution of system (16) on [0,tf]. Then the constraint conditions (10) become(17)hi~xt≤0,t∈0,tf,i=1,…,10,Our goal is to present a state feedback control strategy to maximize the final concentration of 1,3-PD. We now consider the problem of choosing the feedback control coefficientsξk(t), k=1,2,3,4, to minimize the total system cost subject to constraints (17).Problem P. Choose ξ(t)∈Uad to minimize the cost function(18)min⁡J1ξ=-x3tf∣ξs.t.x˙t=f~xt,ξt,hi~xt≤0,i=1,…,10,x0=x0,ξt∈Uad,t∈0,tf.Problem P is a nonlinear optimization problem in which a finite number of decision variables (thefeedback control coefficients) needs to be optimized subject to a set of constraints. It is very difficult to solve Problem P, because each continuous inequality constraint in (17) actually constitutes an infinite number of constraints—one for each point in [0,tf]. Hence, Problem P can be viewed as asemi-infinite optimization problem. Then, we will use a penalty method to transform Problem P [17].The conditionx(t)∈W,t∈[0,tf] is equivalently transcribed into(19)Gξ=0,with(20)Gξ=∫0tf∑i=110max⁡h~ixt,0dt.Clearly, G(ξ)=0 if and only if x(t)∈W. However, the equality constraint (20) is nonsmooth at the points when h~i=0. Consequently, standard optimization routines would have difficulties in dealing with this type of equality constraints. Let(21)G~ξ=∫0tf∑i=110φϵ(h~ixtdt,where the smoothing parameter ϵ is a very small positive number, and φϵ:R→R is defined by(22)φϵη=ηifη>ϵ,η+ϵ24ϵif-ϵ≤η≤ϵ,0otherwise.Obviously, G~(ξ) is a smooth function in ξ. Then, the objective function of Problem P can be reformulated as(23)Jξ=J1ξ+ρG~ξ,where ρ>0 is the given penalty parameter. Problem P can be transformed into the following problem.Problem Q. Choose ξ to minimize the penalty function J(ξ).(24)min⁡Jξ=J1ξ+ρG~ξs.t.x˙t=f~xt,σ,t∈0,tf,x0=x0,ξt∈Uad.Similar to the work [18], we can get the following theorem.Theorem 1.
Letξϵ∗ be the optimal solution of Problem Q. Suppose that there exists an optimal solution ξ∗ of the original Problem P. Then(25)limϵ→0⁡Jξϵ∗=Jξ∗Theorem1 guaranteed that any local solution of the approximate problem can be used for generating a corresponding local solution of the original problem when the smoothing penalty parameter is sufficiently small.


## 4. Control Vector Parameterization Technique and Particle Swarm Adaptive Algorithm

To solve Problem Q numerically, the control vector parameterization approach is applied [19, 20], in which the feedback control variables ξ(t)=[ξ1(t),ξ2(t),ξ3(t),ξ4(t)] are discretized. Partition the time horizon [t0,tf] into p subintervals [tk-1,tk)(k=1,…,p) such that t0<t1<⋯<tp=tf. Using the piecewise-constant policy, the feedback control variable ξi(t) is approximated by (26)ξit≈ξ^it=∑k=1pσi,kχktwhere σi,k is the value of ξ^i(t) in the kth subinterval [tk-1,tk), and χk is defined as (27)χkt≔1,ift∈tk-1,tk,0,otherwise.Letσ=[σ1,σ2,σ3,σ4]T, where σi=[σi,1,…,σi,p].With theξ∈Uad, the differential equation (15) is of form(28)x˙t=f~xt,σ,where(29)f~xt,σ=f~xt,∑k=1pσi,kχi,kt.The initial condition remains:(30)x0=x0.Letx(·∣σ) be the solution of the system (28) corresponding to the control parameter vector σ. And in this way, Problem Q can be approximated by a sequence of nonlinear programming problems, in which σ is regarded as the decision vector.We may now specify the approximate Problem Q(p) as follows.Problem Q(p). Find a control parameter vector σ∈Uad to minimize the cost function.(31)min⁡Jσ=J1σ+ρG~σs.t.x˙t=f~xt,σ,t∈0,tf,x0=x0,σ∈Uad.Theorem 2.
Letξ^∗ be the optimal control of the approximate Problem Q(p). Suppose that the original Problem Q has an optimal control ξ∗. Then, (32)limp→∞⁡Jξ^∗=Jξ∗To solve the Problem P as mathematical programming problems, we require the gradient formulae for the functionJ(σ). We shall derive the required formulae as follows [21].Let the corresponding Hamiltonian function for the cost function be defined by(33)Ht,x,σ,λ=£xt,σ+λTfxt,σ,where (34)£xt,σ=ρ∑i=110φϵh~ixt,and(35)λt=λ1t,λ2t,λ3t,λ4t,λ5tTis the solution of the costate system(36)λ˙t=-∂Ht,x,σ,λT∂x,with the boundary condition(37)λtf=0,0,0,0,0T.Using the similar arguments in [19], we obtain the gradient of J is(38)∂Jσ∂σ=∫0tf∂Ht,xt∣σ,σ,λt∣σ∂σdt.During actual computation, the control parameterization is carried out on an partition of the interval[0,tf]; each component of (38) can be written in a more specific form:(39)∂Jσ∂σi,j=∫tj-1tj∂Ht,xt∣σ,σ,λt∣σ∂ξi^dt.Based on the above control vector parameterization approach, we adopt a gradient-based adaptive refinement method [22] for solving Problem Q(p). The algorithm is adaptive so as to obtain economic and effective discretization grids. In this way, a high-quality solution can be obtained with low computational cost.DefineJ∗l, ξ^i∗l=[σi,1∗l,…,σi,p∗l],(i=1,…,4), Δl=[t0l,…,tpl]T as the optimal objective function value, the optimal solution, and the corresponding discretization time grid found in iteration l. Δl′≔[t0l′,…,t2pl′]T is obtained by bisecting each subinterval in Δl with initial control variable ξil′=[σi,1∗l,σi,1∗l,…,σi,p∗l,σi,p∗l]T. Suppose J∗l′,ξi∗l′=[σi,1∗l′,…,σi,2p∗l′] are the optimal objective function value and the optimal solution in iteration l′, respectively. We hope to find a new discretization grid to make it better adapted to the solution.Define the sensitivity ofσi,jl′ as follows:(40)si,j=∂J∂σi,jl′,whereσi,jl′=σi,j+1/2∗l,where ⌊j+1/2⌋ denotes the maximum integer that does not exceed j+1/2. Suppose σi,K∗l-1 and σi,K∗l are the optimal values in time interval K≔[t2k-2l′,t2kl′] in iteration l-1 and iteration l, respectively.For a given valueε1>0, if(41)σi.K∗l-σi,K∗l-1<ε1,then let(42)si,2k-1=0andsi,2k=0.If the following conditions(43)si,2k-1>λ1si¯orsi,2k>λ1si¯,hold, in which(44)si¯=12p∑j=12psi,j,then the grid point t2k-1l′ in Δl′ is reserved; otherwise, eliminate it. When both t2k-1l′ and t2(k+1)-1l′ are removed, the grid point t2kl′ is also eliminated if(45)si,2k-1<λ2si¯,si,2k<λ2si¯,si,2k+1<λ2si¯,si,2k+1<λ2si¯,andσi,k+1∗l-σi,k∗l<ε2where λ1,λ2, and ε2 are given constants, and λ1>0,λ2∈(0,λ1],ε2>0.The main steps of this algorithm are as follows.Algorithm 3.
Step 0. Choose the initial time grids Δ0, the maximum number of iterations lmax≥1, error tolerance tolJ>0, constants ε1>0, ε2>0, λ1>0, λ2∈(0,λ1]. Initialize the parameters N,l,c1,c2,d1,d2,τ,ωmax,ωmin,Vmax,Vmin,Kmax.
Step 1. Set l=0.
Step 2. By using particle swarm optimization algorithm to obtain the optimal objective function value J∗l and the optimal solution ξ∗l, in which ξ∗l=[ξ1∗l,ξ2∗l,ξ3∗l,ξ4∗l]T, ξi∗l=[σi,1∗l,…,σi,p∗l], i=1,…,4, p is the interval number corresponding to the interval cross powder at this time.
Step 2.1. Randomly generate N particles with uniform distribution on Uad. Denote the position and velocity of particles by ξn=[ξ1n,ξ2n,ξ3n,ξ4n]T∈Uad, in which ξin=[σi,1n,…,σi,pn] and vn=[v1n,v2n,v3n,v4n], respectively, where vin=[θi,1n,…,θi,pn],θi,kn∈[Vmini,Vmaxi],k=1,…,p. Vmini and Vmaxi denote the ith components of Vmin and Vmax. Set Jpbestn is the best objective value found by the nth individual particle, ξn∗=[ξ1n∗,ξ2n∗,ξ3n∗,ξ4n∗] is the best position found by the nth individual particle. Let Jgbest denote the best objective value found by any member of the swarms, ξ∗=[ξ1∗,ξ2∗,ξ3∗,ξ4∗] denote the best position found by any member of the swarms.
Step 2.2. Let k=1, Jpbestn→+∞,Jgbest→+∞.
Step 2.3. For each n=1,2,…,N, use ξn to calculate the corresponding objective function values J(ξn).
Step 2.4. If J(ξn)<Jpbestn, then set Jpbestn=J(ξn) and ξn∗=ξn.
Step 2.5. If Jpbestn<Jgbest, then set Jgbest=Jpbestn and ξ∗=ξn∗.
Step 2.6. If k≤Kmax, then go to Step 2.7; otherwise, stop.
Step 2.7. Update the inertia term according to the following formula: (46)ω=ωmax-ωmin-d1exp⁡1Kmax+d2k-1.
Step 2.8. For each n=1,…,N, compute (47)θi,kn=ωθi,kn+c1r1σi,kn∗-σi,kn+c2r2σi,k∗-σi,kn,
wherer1,r2 obey the uniform distributions on [0,1].Step 2.9. For each n=1,2,…,N, compute (48)σi,kn=σi,kn+θi,knStep 2.10. Set k=k+1, and return go to Step 2.3.Step 3. Check the stopping criterion. If l=lmax or J∗l-J∗l-1/J∗l<tolJ(l>0), stop; otherwise, go to Step 4.Step 4. Refine time grids.Step 4.1. Bisecting each subinterval in Δl to obtain the temporary grids Δl′ and the corresponding control variables ξl′.Step 4.2. Compute the sensitivity according to (42), (43), and (44).Step 4.3. Eliminate unnecessary grid points according to (45), (46), and (47).Step 4.4. Let ξ^l+1=ξ^l′, Δl+1=Δl′.Step 5. Set l=l+1. If l=lmax, stop; otherwise, go to Step 2.Remark 4.
In the above algorithm,N denote the total number of particles in the swarm. c1 and c2 are the cognitive and social scaling parameters. ωmax and ωmin are the maximum and minimum inertia weights. Vmax and Vmin are vectors containing the maximum and minimum particle velocities. Kmax is the maximum number of iteration. d1 and d2 are control factors. k is the iteration index.


## 5. Numerical Results

In the microbial fermentation, we choose the boundary value of state vector asx∗=[0.001,100,0,0,0]T, x∗=[10,2039,939.5,1026,360.9]; the initial concentrations of biomass, glycerol, 1,3-PD, acetate, and ethanol are x10=0.404mmol/L, x20=440.8578mmol/L, x30=0.01mmol/L, x40=0.01mmol/L, and x50=0.01mmol/L, respectively. The control variable Cs0(t)∈[100,1800], D(t)∈[0.05,0.67]. For the parameters in the algorithm, we choose the following values: ε1=10-8, ε2=10-4, λ1=0.2, λ2=0.2, tolJ=10-4, c1=2, c2=2, d1=0.2, d2=0.7, ωmax=0.7, ωmin=0.4, N=8, and Kmax=20. The whole continuous fermentation was implemented with enough substrate. The total fermentation time is taken as 100h. As a matter of experience, for the range of ξ(t) we chose α1=0,β1=400,α2=0,β2=2,α3=0,β3=20,α4=0, and β4=0.0001. Using Algorithm 3, we get the concentration of 1,3-PD at the terminal time is 752.7951 mmol/L. The detailed evolution of time grids is illustrated in Figure 1, after six iterations, the optimal time grids is found, and in this case the results agree with experimental data. Figure 1 also shows the division of time grids in each iteration. The feedback control parameters are shown in Figures 2–5, respectively. The dilution rate and the glycerol concentration in feed are shown in Figures 6 and 7. The concentration changes of biomass, glycerol, 1,3-PD, acetate, and ethanol under the optimal feedback control are shown in Figure 8. The computational results verify the effectiveness of this method.Figure 1
Evolution of the time grids.Figure 2
The first feedback control coefficient.Figure 3
The second feedback control coefficient.Figure 4
The third feedback control coefficient.Figure 5
The fourth feedback control coefficient.Figure 6
The dilution rate.Figure 7
The glycerol concentration in feed.Figure 8
The concentration changes of biomass, glycerol, 1,3-PD, acetate, and ethanol.


## 6. Conclusions

In this paper, we have considered a feedback control strategy which is close-loop control for producing 1,3-PD in microbial continuous fermentation and developed a particle swarm adaptive algorithm to obtain the global solution. Numerical results show that the method is successful at producing high-quality control strategies.



---

*Source: 6252607-2019-03-24.xml*