Abstract: In this paper, we consider mathematical modelling and parameter identification problem in bioconversion of glycerol to 1,3-propanediol by Klebsiella pneumoniae. In view of the dynamic behavior with memory and heredity and experimental results in batch culture, a two-stage fractional dynamical system with unknown fractional orders and unknown kinetic parameters is proposed to describe the fermentation process. For this system, some important properties of the solution are discussed. Then, taking the weighted least-squares error between the computational values and the experimental data as the performance index, a parameter identification model subject to continuous state inequality constraints is presented. An exact penalty method is introduced to transform the parameter identification problem into the one only with box constraints. On this basis, we develop a parallel Particle Swarm Optimization algorithm to find the optimal fractional orders and kinetic parameters. Finally, numerical results show that the model can reasonably describe the batch fermentation process, as well as the effectiveness of the developed algorithm.
Keywords: fractional dynamical system, parameter identification, parallel optimization, batch fermentation.
Articles
Modelling and parameter identification for a two-stage fractional dynamical system in microbial batch process*
Chongyang Liu chongyangliu@aliyun.com
Recepción: 10 Febrero 2021
Revisado: 17 Julio 2021
Publicación: 02 Febrero 2022
Fractional calculus, also known as noninteger-order calculus in the literature, is a generalization of the ordinary calculus. From the perspective of mathematical classification it is a branch of mathematical analysis [5]. Fractional calculus has been widely used in numerous areas including control systems [1], bioengineering [21], viscoelasticmechanics [23], and image encryption [32]. However, its development was very slow before 1970. In 1970s, researchers discovered that a series of phenomena or processes, such as fractal geometry, power law phenomena, and memory processes, is closely related to fractional calculus. In particular, fractional differential equations (FDEs) are naturally related to system with memory and heredity, and this phenomenon is often reflected in the field of bioengineering [21, 30]. As a result, the applications of FDEs in complex biological system with memory and heredity have received particular attention during the past decades [11, 21, 30, 31]. It is worth noting that there are two most popular fractional derivatives, i.e., the Caputo fractional derivative and the Riemann–Liouville fractional derivative. Compared with the initial value problem of FDEs in the Riemann– Liouville sense, differential equations in the Caputo sense are initialized by the integer- order derivatives, which enhances its practicability in mathematical modelling of real- world processes. Thus, FDEs in the Caputo sense will be considered in this paper.
1,3-propanediol (1,3-PD), a colorless, hygroscopic viscous liquid, is widely used in the production of polyester, polyurethane, and heterocyclic compounds [15]. There are two main routes to produce 1,3-PD: chemical synthesis and microbial fermentation. The latter one is becoming increasingly attractive in industrial 1,3-PD production since it has the advantages of mild conditions, few by-products, and environment friendly. 1,3-PD microbial production can be used one of three fermentation modes: batch culture, contin- uous culture, and fed-batch culture. In a batch culture, a quantity of microorganisms and substrate is poured into the reactor only once at the beginning of the process and nothing is added to, or removed from the reactor during the culture process. In a continuous culture, the substrate is continuously injected into the reactor at a certain rate. At the same time, the culture fluid is taken out at the same rate, and thus the volume of the fermentation broth in the reactor remains unchanged. The fed-batch culture process is implemented by switching between batch and feeding processes. Compared with the continuous and fed- bach cultures, the batch culture is rather simple to operate, and it can be regarded as a part of fed-batch culture. In particular, glycerol fermentation to 1,3-PD in batch culture can obtain the highest productivity and yield [9]. Therefore, in this paper, we will focus on the 1,3-PD batch production. For 1,3-PD batch production, extensive studies have been carried out in the literature. Mathematical models of glycerol fermentation by Klebsiella pneumoniae were proposed in [27,35]. On the basis of these models, various identification problems and optimization algorithms have been investigated in [18–20,39,40]. However, these parameter identification problems involve only one-stage systems. Recently, nonlin- ear multistage systems have been proposed to formulate the batch fermentation process in [6, 12, 37]. Nevertheless, all the dynamic systems mentioned above are based on the classical ordinary differential equations. As is well known, the physical and chemical properties of the fermentation process will lead to the existence of memory and hereditary behavior, and the fractional calculus can reasonably represent the behavioral systems with memory and heredity in [11, 30]. By the way, various forms of fractional differential equations and their practical applications have been discussed in [3,10,33]. Also, there are many successful numerical methods for solving the fractional optimal control problems; see, for example, [7, 16, 25]. More recently, a fractional dynamic system in batch culture and its parameter identification problem have been investigated in [24]. Although the fractional system in [24] can describe the first two stages (developmental phase and growth phase) of the batch culture well, the simulation results in the stationary phase fail to meet the expectation. Hence, fractional dynamic system with only one stage cannot describe the whole batch fermentation process well.
In this paper, a two-stage fractional differential dynamical system in the sense of Caputo is proposed to formulate the batch fermentation process. For this system, some important properties of the solution are discussed. By taking the weighted least-squares error function as the cost function, we then present a parameter identification model subject to continuous state inequality constraints. To solve this parameter identification problem, an exact penalty method is used to transform the problem into the one only with box constraints. Furthermore, considering the fact that the number of fractional orders and kinetic parameters to be identified is very large, we develop a parallel Particle Swarm Optimization (PPSO) algorithm to solve the transformed problem. Numerical simulations show that the two-stage fractional dynamical system can reasonably describe the fermentation process, and the developed algorithm is applicable and effective.
The remaining of the paper is organized as follows. A two-stage fractional dynamical system of batch culture is established and some properties of the solution are discussed in Section 2. A parameter identification model is given in Section 3. A numerical solution approach based on the exact penalty method and PPSO algorithm is proposed in Section 4. Numerical simulations are given in Section 5. Conclusions are finally drawn in Section 6.
Before establishing the fractional fermentation model, some of the most popular defini- tions in fractional calculus will be introduced, and we refer the reader to book [2] for detail.
Definition 1. For a function 
, the Riemann–Liouville fractional integral of order 
is defined by
where 
is the set of Lebesgue measurable functions from 
, and 
is the Gamma fuction
Definition 2. For a function 
, if 
, then the Caputo fractional derivative of order 
is defined by
where 
is the smallest integer greater than or equal to 
, and . 
denotes the standard 
-order derivative of function 
.
Note that the Caputo fractional derivative of order 
can be converted to the one of order within range 
[42]. Thus, for simplicity, we assume that 
throughout the paper.
At the beginning of the batch fermentation process, an appropriate number of substrates, i.e., biomass and glycerol, is added to the reactor, and no input or output is performed during the process. In particular, the memory and hereditary effects in enzyme reactions and microbial growth can be reasonably described by the characteristics of fractional calculus [11, 21, 30].
Based on the work [6, 24], we formulate the batch fermentation process as the following two-stage fractional differential equations:
where 
, are, respectively, the concentrations of biomass, glycerol, 1,3-PD, ethanol, and acetate at time 
, are fractional orders; 
is a given switching time; 
is a given terminal time of the fermentation process; 
, denote the kinetic parameters of biomass growth, glycerol consumption, 1,3-PD formation, ethanol formation, and acetate formation, respectively; and 
, denote the inhibitory effects of cell death on cell growth, glycerol consumption, 1,3-PD formation, ethanol formation, and acetate formation, respectively. For brevity, we define 
. Moreover, denote the right-hand side of equations (1) and (2) by
and
Then equations (1) and (2) can be rewritten as the following two-stage fractional dynamical system:
where 
is a given initial state; and 
is the right-hand limit of 
. Note that the two-stage dynamical system (5) and (6) is essentially similar to the short-memory fractional differential equations mentioned in [34], both of which make full use of the memory characteristics of the Caputo derivative in a given time horizon.
In system (5) and (6), the fractional orders and kinetic parameters need to be identified.
Thus, define the admissible set of fractional orders as
where 
is a given constant. Furthermore, define the admissible set of kinetic parameters as
According to the biological significance, there are critical concentrations of biomass, glycerol, 1,3-PD, ethanol, and acetate in the batch fermentation process. Outside these critical concentrations, cells will cease to grow or even die. Therefore, define the admissible set 
as
where 
, are the lower concen- tration thresholds of biomass, glycerol, 1,3-PD, ethanol, and acetate for cell growth, respectively; and
are the corresponding upper concentration thresholds [24].
For system (5) and (6), some important properties will be discussed, e.g., the existence and uniqueness of the solution and the continuity of the solution with respect to fractional orders and kinetic parameters.
Property 1. For any 
, defined by (5)and (6) satisfy
are continuously differentiable with respect to 
and 
.
such that, for any 
for all 
were 
is the Euclidean norm.
Proof. (i) According to the expressions of 
, in (3) and (4), this conclusion can be obtained. (ii) The boundedness of 
can be verified by (1) and the compactness of 
and 
.
Property 2. The functions 
, defined by (3)and (4)are Lipschitz continuous, that is, for any 
, there exists a constant 
such that
Proof. By Property 1 and the differential mean value inequality [26] we have, for any 
where 
. Furthermore, by the compactness of 
and 
there exists a constant 
such that
Combining (9), (10) with (8) gives
for any 
Theorem 1. For each 
, the two-stage fractional system defined by (5)and (6)has a unique continuous solution 
. Moreover, 
, satisfy the following integral equations:
Proof. Using a similar proof as given for Theorem 3.4 in [42] with Properties 1 and 2, we can complete the proof.
Theorem 2.The solution 
of system (5)and (6)is continuous with respect to 
Proof. Case 1: 
For any 
, there exists two real numbers 
and 
such that 
. Furthermore,
According to Properties 1 and 2, the first part of inequality (11) can be rewritten as
where
By generalized Gronwall inequality in [36] we have
Furthermore, the second part of inequality (11) can be expressed as
where
Similarly, we have
Note that 
. Thus, from (12) and (13) we have 
as,
Case 2: 
.
From Case 1, for 
, we have
Then, using the similar proof as given for Case 1, we conclude
where
Similarly,
where
Note that 
. Thus, from (14) and (15) we have 
wich imply 
Based on Cases 1 and 2, we obtain that 
as 
, wich completes the proof
Parameter identification is the problem of adjusting the values of parameter to make the predicted values of the system consistent with the experimental data as much as possible. In this section, we will discuss the parameter identification problem involving the two- stage fractional dynamical system.
In batch culture, we have 
measured experimental data. Let 
be the measured concentrations of biomass, glycerol, 1,3-PD, ethanol, and acetate at the measured moment 
.Furthermore, considering that the measured data at the later stage are more practical than the data measured in the initial stage, we propose the following weighted least-squares error function [14]:
where 
is the calculated value for the 
component at time 
; and 
is the maximum measured concentration of the 
component. Then the parameter identification model (PIM) in the batch culture can be stated as
(PIM) min 
such that 
Theorem 3. For problem (PIM), there exists at least one optimal pair 
such that
Proof. By Theorem 2 the solution 
of system (5) and (6) is continuous on 
. Then we define the feasible set of problem (PIM) as
Clearly, 
is a nonempty set. According to the compactness of sets 
and , 
is a bounded set. Let
denote any sequence, and let 
as 
. Based on Theorem 2 and the compactness of 
, we have 
. Thus, 
, which indicates that the feasible set 
is closed. Furthermore, since the cost function 
is also continuous on 
by (16), we conclude that problem (PIM) has at least oe optimal pair such that (17) holds. The proof is complete.
Problem (PIM) is a parameter optimization problem subject to continuous state inequality constraints (7). In this section, a numerical solution approach will be developed based on an exact penalty method and a PPSO algorithm.
In solving problem (PIM), it is computationally difficult since continuous state inequality constraints (7) must hold at infinite number of points in 
. Fortunately, constraints transcription techniques [8, 29] and exact penalty methods [17, 38] have been proposed to overcome these difficulties. In the exact penalty method [17], it only requires that the penalty parameter is large enough but finite, and its adjustable parameters are fewer. Thus, the exact penalty method [17] is applied to problem (PIM) to deal with continuous state inequality constraints (7).
Let
Then the continuous state inequality constrains (7) can be equivalently converted to the following form:
where 
Let 
be a given constant, and let 
be a new decision variable. The penalty fuction is defined as
where 
is a penalty parameter; and 
and 
are two positive constants satisfying 
Then we can transform problem (PIM) into the penalty problem as follows:
Based on the derivation in [17], it is clear that 
is an exact penalty function. As a result, the optimal solution of problem (PIM) can be obtained by solving a sequence of problem 
For each 
, problem 
is a parameter optimization problem. Various optimization methods can be selected to find the optimal fractional orders and kinetic parameters such as gradient-based algorithm [28]. Nevertheless, the gradient-based optimization method is easy to trap into the local optimum, which is obviously not desired. Furthermore, because there are a large number of parameters to be identified when solving problem 
, the time cost of estimating candidate parameters is quite expensive. As a result, we will develop a PPSO algorithm to solve problem 
for each 
.
The Particle Swarm Optimization (PSO) algorithm, proposed by Professor Eberhart and Dr. Kennedy in 1995 [13], is a swarm cooperation-based random search algorithm developed by simulating the foraging behavior of birds in a group. In a standard PSO algorithm, each particle is regarded as a feasible solution to the optimization problem. Each particle has a current velocity and flies at that speed in a given space. In order to get the optimal solution of the optimization problem, each particle should fly to the globally optimized position through dynamically adjusting the flight speed by its own and the group’s experiences. With the rapid growth of computing size, the serial PSO algorithm will produce high computational cost. Thus, some parallel PSO algorithms were developed in [22, 39, 41]. However, problem 
is a penalty problem with penalty parameter 
As a result, we propose the following PPSO algorithm to solve problem
for each 
. Here 
denotes the decision vector of problem 
Algorithm 1
Step 1. Allocate 
slave processors. Initialize the total number of particles 
. Compute the number of particles on each processor 
. Initialize data and variables on the master processor.
and 
, the tolerance parameter 
, the cognitive and social parameters 
and 
, a sufficient large number Λ, and the maximal iteration 
.
, the global optimal fitness value 
Step 2. Broadcast data and variables on master processor to all slave processors.
Step 3. Execute the following steps on each slave processor, and denote the identifier of each slave processor by 
.
decision vectors from the horizon 
, randomly initialize 
velocities from the uniform distribution 
, denote the decision vector and velocity of particles by 
.
and 
by
(iv) if 
and go step 3 (iv). Otherwise, go to sptep 4.
Step 4. Gather 
, from slave processor 
into the master processor.
Step 5. Assign the master processor to execute the following operations.
and record 
as follows:
2 If 
and stop. Otherwise, set 
, brodcast 
and 
to all slave processors
Step 6. Execute the following procedure on each slave processor.
(i) Update particle velocities and positions by the following iterative formulas:
were 
random numbers in [0, 1]; and 
If 
violates the bound constraint, then execute the following operations:
and go to step 3 (iii).
On the basis of Algorithm 1, we propose the following algorithm to solve problem (PIM).
Algorithm 2
Step 1. Choose initial values of 
Step 2. Solve problem
by Algorithm 1 to give 
Step 3. If 
, then set 
and go to Step 4. Otherwise, take 
as an optimal solution of problem (PIM) and stop.
Step 4. If 
, then stop and output abnormal exit. Otherwise, set 
.
Remark 1. In Algorithm 2, 
is a penalty parameter; 
is an increment factor; 
is a maximum penalty parameter; 
is an error tolerance. If Algorithm 2 is abnormally terminated in Step 4, the parameters 
and 
can be modified to resume Algorithm 2.
In the numerical simulation, Algorithm 2 is used to solve problem (PIM) to find the optimal fractional orders and kinetic parameters, and all computations are implemented in Matlab 2018a environment on a Intel Core i5-7400 (64-bit, 8GB RAM, 3.4GHz) machine. Here the two-stage fractional dynamical system is solved by the implicit trapezoidal product-integration rule [4], and the step-size for integration is taken as 
. The initial state, the switching moment, and the terminal time are taken as 
, respectively [24]. In Algorithm 1,
. In Algorithm 2, the initial values of 
, and 
are, respectively, 
. The above parameters are determined on the basis of a large number of experiments.
Based on the obtained optimal fractional orders and kinetic parameters, we plot the concentrations of biomass, glycerol, 1,3-PD, ethanol, and acetate in Fig. 1. For compari- son, the one-stage fractional dynamic system with fractional orders and kinetic parameters in [24] is also solved. The obtained concentrations of biomass, glycerol, and 1,3-PD are also plotted in Fig. 1. From Fig. 1 we can see that, compared with the results in [24], our computed concentrations of biomass, glycerol, 1,3-PD, ethanol, and acetate can better describe the experimental data in [35]. In addition, Table 2 shows the relative errors between the calculated values and experimental data in this work, as well as relative errors in [24], where the relative errors are defined as
From Table 2 it can be seen that the relative errors of the two-stage fractional model in this work are significantly smaller than those in [24]. This also confirms that our two- stage fractional dynamical system can reasonably describe the batch fermentation process. To further test the performance of our proposed PPSO algorithm, a serial PSO (SPSO) algorithm with the same parameters as in PPSO algorithm is also developed to solve problem 
. We perform 30 test runs in solving problem 
with SPSO and PPSO algorithms. The obtained optimal costs, worst costs, average costs, and average iteration time are listed in Table 3. From Table 3 it can be seen that the optimal cost, the worst cost, the average cost, and the average iteration time obtained by our PPSO algorithm are all superior to those obtained by the SPSO algorithm. The convergence curves of the cost function by the PPSO and SPSO algorithms are also depicted in Fig. 2. From Fig. 2 we can see that the convergence of PPSO algorithm is faster than SPSO algorithm. From the above results we confirm that the developed PPSO algorithm is highly effective and efficient in solving problem (PIM).
In this paper, we have considered the parameter identification problem in batch process. A two-stage fractional dynamical system in the sense of Caputo is proposed to describe the batch process. Taking the error between the calculated values and the experimental data as the performance index, we present a parameter identification model subject to continuous state inequality constraints. By applying an exact penalty method the parameter identification problem is transformed into the one only with box constraints. A parallel Particle Swarm Optimization algorithm is developed to solve the resulting problem. Numerical results show that our proposed two-stage fractional dynamical system is reasonable and the proposed parallel algorithm is efficient. For further research, it will be of interest to investigate the optimal switching control problem involving the proposed two-stage fractional dynamical system.
