Design and Application of Intelligence Algorithms in Continuous Fermentation of Glycerol

The bioconversion of 1,3-propanediol from glycerol by Klebsiella pneumoniae can be described by a nonlinear dynamic system. Some work has been done on the identification and optimization of the system, in which the dilution rate of glycerol is considered as a constant. However, the demand of glycerol may vary at different fermentation stages, it is reasonable to view the glycerol metabolic system with dilution rate varying with time. In this paper, we model the glycerol metabolic process as a fourteen-dimensional nonlinear dynamical system, where the dilution rate is considered varying with time. Then an optimal discrete-valued control problem for maximizing the average concentration of 1,3-propanediol in the fermentation process is established. To solve the optimization problem, auxiliary control and an exact penalty function are used to convert this problem into a large-scale parameter optimization problem. For better balancing local and global search ability, a competitive particle swarm algorithm with time-varying control factors is proposed which is proved to be faster and more stable than the traditional competitive particle swarm algorithm. Numerical experiments are conducted to show the rationality, effectiveness and applicability of the method proposed.


I. INTRODUCTION
1,3-Propanediol (1,3-PD) is an important chemical material, in particular as a monomer for polyesters, polyethers and polyurethanes [1]. The production of 1,3-PD by chemical method has high production cost and easy to cause environmental pollution, so the microbial production is recently paid more attention for its low cost, high production and no pollution [2]. Glycerol, a by-product of the soap and detergent industry, can be converted to 1,3-PD by Klebsiella pneumoniae (K. pneumoniae) under anaerobic conditions [3]. This microbial production of 1,3-PD has been widely investigated in recent years because of its high productivity [4].
There are three fermentation methods for the production of 1,3-PD by microbial fermentation with glycerol as substrate, i.e., continuous fermentation, batch fermentation and fed-batch fermentation. In this paper, we consider a continuous fermentation process, in which glycerol is fed into the The associate editor coordinating the review of this manuscript and approving it for publication was Emre Koyuncu . fermentor continuously, the broth in fermentor pours out at some rates and the volume of the fermentation broth remains unchanged.
Glycerol dismutation by K . pneumoniae is a complex bioprocess, covering both extracellular and intracellular environments. In the intracellular environments, two parallel pathways which including oxidative pathway and reductive pathway are coupled. The goal product 1,3-PD is produced in the reductive pathway. In the coupled oxidative pathway, byproducts acetate and ethanol are generated. The microbial growth is subjected to multiple inhibitions of substrate and products. Sun et al. [5] presented a fourteen-dimensional nonlinear dynamic system to describe the continuous fermentation and multiplicity analysis, considering two regulated negative-feedback mechanisms of repression and enzyme inhibition. Ye et al. [6] studied the concentration robustness of this fourteen-dimensional glycerol metabolism system and proposed a robustness index to measure the concentration robustness of the considered system. By defining the timevarying metabolic network structure as an integer-valued VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ function, Wang et al. [7] modeled glycerol metabolism in continuous fermentation as a fourteen-dimensional nonlinear mixed-integer dynamic system and identified the dynamic network structure and kinetic parameters. In above works [5]- [7], dilution rate of the glycerol was both considered as a constant. However, the demand of glycerol may vary at different fermentation stages [8]- [10]. In this paper, based on the fourteen-dimensional nonlinear dynamical system, we will consider dilution rate varying with time and construct an optimal discrete-valued control problem for maximizing the average concentration of 1,3-PD in the fermentation process.
Classical optimal control methods are not applicable to the optimal discrete-valued control problem for its control being a discrete variable. In previous works [11], [12], methods including outer convexification, relaxation, rounding strategy and time transformation were used to solve the mixedinteger optimal control problem. In [13], Yu et al. transformed the optimal discrete-valued control problem into an equivalent continuous problem by introducing auxiliary controls and applying a time-scaling transformation, and the transformed problem was then solved by an exact penalty function approach. In this paper, we will use the methods which described in [13] to transform our problem into a large-scale parameter optimization problem. At present, it is popular to use algorithms under competition mechanism to solve largescale parameter optimization problems [14]- [16]. This work focuses on proposing a competitive particle swarm algorithm with time-varying control factors to solve the large-scale parameter optimization problem. Numerical results show that the method proposed is effective and applicable.
The rest of this paper is organized as follows. In Section II, a nonlinear hybrid dynamical system of glycerol continuous fermentation is established. Section III is devoted to proposing an optimal discrete-valued control problem with the objective function of maximizing the average concentration of 1,3-PD in the reaction process, and a large-scale parameter optimization problem is obtained. In Section IV, a competitive particle swarm optimization algorithm with time-varying control factors is constructed to solve the problem which we obtained. Numerical example and results are presented in Section V. Conclusion is made at the end of this paper.

II. NONLINEAR HYBRID DYNAMICAL SYSTEM
According to the factual experiments, we assume that (H1) Glycerol is the only substrate that added to the reactor during the process of continuous culture, and no medium is pumped inside or outside the reactor in the whole process.
(H2) The concentrations of reactants are uniform in the reactor, time delay and nonuniform space distribution are ignored.
Let x(t) = (x 1 (t), x 2 (t), · · · , x 14 (t)) T , the components of which represent the concentrations of biomass, extracellular glycerol, extracellular 1,3-PD, acetate, ethanol, intracellular glycerol, 3-HPA, intracellular 1,3-PD, m R (mRNA coding repressor), R (free repressor), m GDHt , GDHt, m PDOR , and PDOR at time t, respectively. According to the results of [6], [17], we assume that glycerol passes the membrane by passive diffusion and 1,3-PD by passive diffusion coupled with active transport, and the inhibition of 3-HPA on cell growth (GDHt activity, PDOR activity) exists at any concentration of 3-HPA. Under these assumptions, the velocity field of concentration changes during the process of glycerol continuous fermentation is formulated as follows according to the previous work [6]: Here, C s0 is the initial concentration of glycerol. p i (i = 1, 2, · · · , 26) are intracellular kinetic parameters, and the values of them are listed in Table 1 according to [5], [7]. According to [6], the specific cell growth rate µ(t), specific consumption rate of extracellular glycerol r 2 (t) and specific product formation rates r i (t)(i = 3, 4, 5) are expressed as follows: where µ m is the maximum specific growth rate, P i (i = 1, 2, · · · , 15) are extracellular kinetic parameters, and x * i (i = 2, 3, 4, 5) are critical concentrations. Their values are given in Table 2 according to the literature [4].
14, the right-hand side of the ith equation of (1) - (14). The nonlinear dynamical system of glycerol continuous fermentation can be described by Here, d(t) is dilution rate of glycerol in the fermentation process and the value of it will vary with time. But in experiments, d(t) is difficult to be continuously modified, so it is considered taking m different values d 1 , d 2 , · · · , d m in our model.
x 0 is the given initial state and [0, t f ] is the total fermentation time interval. The admissible set of x(·) is denoted Table 2 according to [4].
Based on the results given in [18], the following properties for system (20) can be easily proved.
Property 1: For a fixed x 0 ∈ X and a given d(·) ∈ D, the function f (x(·), d(·)) defined in (1) -(14) is Lipschitz continuous in x(·) on X , and there exist positive constants a, b such that the linear growth condition holds, i.e., Property 2: For a fixed x 0 ∈ X and a given d(·) ∈ D, there exists a unique solution to system (20), denoted by x(·, d(·)).

III. OPTIMAL CONTROL PROBLEM AND REFORMULATION A. OPTIMAL CONTROL PROBLEM
In this paper, for dynamic system (20), we consider maximization of the average concentration of 1,3-PD in fermentation process as the optimization objective and give the optimal control problem as follows: Here, d(t) is a discrete-valued control variable, so problem (P1) is an optimal discrete-valued control problem. Next, according to [13], we will transform problem (P1) to an equivalent continuous optimal control problem by introducing auxiliary controls.
m, are auxiliary control functions, and satisfy the following constraints: Obviously, v j (t) is a continuous map from [0, t f ] to [0, 1], but it only takes the value from {0, 1}.
Next, by letting d(t) = m j=1 v j (t)d j , problem (P1) then be transformed into the following problem (P2): Problem (P1) and (P2) are equivalent and the following theorem can be easily proved by the result of [13].
is an optimal control for problem (P2) if and only if d * (t) is an optimal control for problem (P1).

C. DISCRETIZATION
Suppose that system (20) has at most N − 1 switches, and the time is divided into N equidistance segments. Let τ k denote the kth switching time. Then where ξ jk is the value of v j (t) on [τ k−1 , τ k ), and χ I is the indicator function of I defined by The constraint condition of v j (t) become: Let ξ j = (ξ j1 , ξ j2 , · · · , ξ jN ) T ∈ R N , and ξ = [ξ 1 , ξ 2 , · · · , ξ m ] T ∈ R m×N . Problem (P2) can be written equivalently to problem (P3) as follows: Problem (P3) is a large-scale parameter optimization problem with m·N decision variables. For the quadratic inequality constraints in problem (P3) being difficult to be satisfied, problem (P3) is not easy to be solved directly. According to [13], we will use an exact penalty function to transform problem (P3) into an unconstrained optimization problem which is easy to be solved.

D. AN EXACT PENALTY FUNCTION
For problem (P3), according to [13], we construct the following exact penalty function: if ε = 0, and ξ is feasible for problem (P3), Here, ε is a new decision variable, (ξ, ε) is the constraint violation and be defined as follows: where α, β, γ are positive real numbers, and σ is a penalty parameter. To continue, we define It is clear that S 0 is the feasible solution set of problem (P3) for ε = 0. Now, we can consider the following optimization problem: Here, if ε = 0, and ξ is feasible for problem (P3), dt. According to [13], it is easy to get that under some appropriate assumptions, for a sufficiently small penalty parameter, a local optimal solution of problem (P4) is a local optimal solution of problem (P3). This solution can be used as a corresponding local solution of problem (P1).

IV. A COMPETITIVE PARTICLE SWARM ALGORITHM WITH TIME-VARYING CONTROL FACTORS
Considering that problem (P4) is a large-scale optimization problem in which the objective function is not differentiable, and the gradient-based methods cannot be used, we will use intelligent algorithms to solve it.
Competitive swarm optimizer (CSO), introduced by Cheng et al., is one of the most popular algorithms for solving large-scale optimization problems which uses the competitive mechanism between particles within a single swarm and updates half of the particles each time [14], [19], [20]. Mohapatra et al. [15] proposed an improved version of CSO (called MCSO), in which 2/3rd of the particles are updated each time. For increasing the convergence speed in early iterations, Wang et al. [7] proposed a new formula for inertia weight in the competitive particle swarm optimization algorithm. For balancing local search capability and global search capability, based on the research of Mohapatra et al. [15] and Wang et al. [7], we propose a competitive particle swarm algorithm with time-varying control factors.
Denote N size as the maximum iterations and N size as swarm size, which is a multiple of 3 that favors the tri-break-up mechanism. In each iteration, three particles are picked up arbitrarily from the swarm to undergo a tri-competition. The particle that has the highest fitness is named as the winner, and the rest are the losers (called loser 1 and loser 2). Let η w j (s),η l1 j (s),η l2 j (s) and v w j (s),v l1 j (s),v l2 j (s) be the position and velocity of the winner, loser 1 (l1) and loser 2 (l2) in the jth round of competition, here j ∈ I N size 3 and s is the generation number. The winner passes directly to the next generation. The update formulas of velocity and position of two losers in the basic version of MCSO are described as follows [15]: v a j (s + 1) = R a 1,j (s)v a j (s) + R a 2,j (s)(η w j (s) − η a j (s)) + R a 3,j (s)φ 1 (η(s) − η a j (s)), (22) η a j (s + 1) = η a j (s) + v a j (s + 1), (a = l1, l2).
Here R l1 i,j (s) and R l2 i,j (s)(i = 1, 2, 3) take different random values of (0,1) in the jth round of competition in iteration s. η(s) is the mean position value of the relevant particles. φ 1 is the parameter that controls the influence of η(s). In [7], Wang et al. changed the velocity update formula (22) to: where η best is the global best position value of the relevant particles and the inertia weight w(s) = 0.9 − 0.5 × s N ite . φ 2 is the parameter that controls the influence of η best . Now we expect that the local search ability is strong and the global search ability is weak at the initial stage of the algorithm, so that the particle search is more sophisticated. For another, we expect the local search ability to be weakened and the global search ability to be enhanced at the later stage of the algorithm. Based on this idea, the new velocity update formula (24) becomes: v a j (s + 1) = w(s) · v a j (s) + R a 1,j (s)φ 1 (η w j (s) − η a j (s)) + R a 2,j (s)φ 2 (η(s) − η a j (s)) + R a 3,j (s)φ 2 (η best − η a j (s)).
Here, φ 1 is constructed as a monotonically decreasing function that controls the influence of η w j (s), and φ 2 is constructed as a monotonically increasing function that controls the influence of η(s) as well as η best . The expressions of φ 1 and φ 2 are as follows: where φ 1start , φ 1end , and φ 2start , φ 2end are initial and final values of φ 1 and φ 2 respectively. Let (s) represents the set of particles that have not yet been picked up to compete in iteration s. The basic steps of the competitive particle swarm optimization algorithm with time-varying control factors (CPSOT) to solve problem (P4) are given in Algorithm 1.
For investigating the performance of CPSOT for problem (P4), we compare CPSOT with the algorithm MCSO proposed in [15] and the algorithm CPSO proposed in [7] by performing a set of numerical experiments. The experiments are implemented with MATLAB R2019b. Table 3 lists the mean time consumption (Time), mean objective value (Mean) and variance (Var) obtained by MCSO, CPSO and CPSOT through 10 independent runs for different N size (N size = 420, 600). As can be seen from Table 3, under the same maximum number of Fitness, CPSOT has the largest Fitness value, the least runtime and the smallest variance. The largest Fitness value and the least runtime mean that CPSOT has good solving accuracy and fast speed. The smallest variance indicates that CPSOT has good stability. Additionally, the result of N size = 600 is poorer than that of N size = 420 and this implies that the performance of the algorithms does not rely much on large swarm size [7], [14]. Fig. 1 shows the descending curves of the mean objective values of MCSO, CPSO and CPSOT. From Fig. 1, we can see that CPSOT converges faster than other two algorithms. Hence, CPSOT is effective and preferred for our proposed problem (P4). Algorithm 1 Pseudocode of the competitive particle swarm optimization algorithm with time-varying control factors (CPSOT).
Step 16. return η * = η best . According to the above experimental results, the CPSOT algorithm under N size = 420 is used to simulate the process of microbial fermentation to produce 1,3-PD, and the changes of concentrations of biomass, extracellular glycerol and 1,3-PD   with time are obtained (Fig. 2). The glycerol dilution rate switching sequence, i.e., the order of d(t)'s values along the time, is listed in Table 4. And the trajectories of d(t) are drawn in Fig. 3. Fig. 2 shows that the nonlinear dynamic system constructed by us can reasonably describe the real continuous fermentation process of glycerol. Additionally, the concentration curve shows periodic fluctuation. The concentration of 1,3-PD gradually increased to a peak as the glycerol concentration gradually decreased to a trough due to consumption. Such a periodic change process has a certain reference significance for the related experimental operation in the laboratory, and it has a reference value for improving the volume productivity and the substrate conversion rate.

VI. CONCLUSION
In this paper, the metabolic process of 1,3-PD produced by biological fermentation was studied, and an optimal discretevalued control problem was established and solved, in which the dilution rate was considered varying with time. To solve this problem, the method of using auxiliary control and an exact penalty function converted it into a large-scale parameter optimization problem, then we constructed a competitive particle swarm algorithm with time-varying control factors to solve the large-scale problem.
Numerical results show that it is feasible and effective to use exact penalty function method to solve the large-scale parameter optimization problems with quadratic inequality constraints, which provides a solution idea for optimal discrete-valued control problems. In addition, the competitive particle swarm algorithm with time-varying control factors is faster and more stable than the traditional competitive particle swarm algorithm in solving large-scale parameter optimization problems. According to the results of this paper, the proposed algorithm can be applied to other large-scale parameter optimization problems, which has certain practical significance.