Modelling and parameter identification for a two-stage fractional dynamical system in microbial batch process

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.


Introduction
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], viscoelastic 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.
2 Two-stage fractional dynamical system and its properties in batch process

Fractional calculus
Before establishing the fractional fermentation model, some of the most popular definitions in fractional calculus will be introduced, and we refer the reader to book [2] for detail.

Definition 1.
For a function f ∈ L 1 ([a, b], R), the Riemann-Liouville fractional integral of order α > 0 is defined by where L 1 ([a, b], R) is the set of Lebesgue measurable functions from [a, b] to R, and Γ(·) is the Gamma function.

Two-stage fractional dynamical system in batch culture
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].
In system (5) and (6), the fractional orders and kinetic parameters need to be identified. Thus, define the admissible set of fractional orders 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 W as where are the lower concentration thresholds of biomass, glycerol, 1,3-PD, ethanol, and acetate for cell growth, respectively; and x * 4 = 360.9 mmol L −1 , and x * 5 = 1026 mmol L −1 are the corresponding upper concentration thresholds [24].

Properties of the two-stage fractional dynamical system
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 t ∈ [0, T ], x(t) ∈ W , and p ∈ P, f , = 1, 2, defined by (5) and (6) satisfy: (i) f are continuously differentiable with respect to x(t) and p.
Proof. (i) According to the expressions of f , = 1, 2, in (3) and (4), this conclusion can be obtained. (ii) The boundedness of f can be verified by (1) and the compactness of W and P.

Parameter identification problem
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 twostage fractional dynamical system.
In batch culture, we have N s measured experimental data. Let y l = (y l 1 , y l 2 , . . . , y l 5 ) be the measured concentrations of biomass, glycerol, 1,3-PD, ethanol, and acetate at the measured moment t l , l ∈ {1, 2, . . . , N s }. Furthermore, considering that the measured data at the later stage are more practical than the data measured in the initial stage, we https://www.journals.vu.lt/nonlinear-analysis propose the following weighted least-squares error function [14]: Theorem 3. For problem (PIM), there exists at least one optimal pair (α * , p * ) ∈ F × P such that J(α * , p * ) J(α, p) ∀(α, p) ∈ F × P.
Then we define the feasible set of problem (PIM) as Clearly, V is a nonempty set. According to the compactness of sets F and P, V is a bounded set. Let {(α , p )} ∞ =1 ⊆ V denote any sequence, and let (α , p ) → (α, p) as → ∞. Based on Theorem 2 and the compactness of W , we have x(t|α, p) ∈ W , t ∈ [0, T ]. Thus, (α, p) ∈ V , which indicates that the feasible set V is closed. Furthermore, since the cost function J(α, p) is also continuous on F × P by (16), we conclude that problem (PIM) has at least one optimal pair such that (17) holds. The proof is complete.

Numerical solution approach
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.

Exact penalty method
In solving problem (PIM), it is computationally difficult since continuous state inequality constraints (7) must hold at infinite number of points in [0, T ]. 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).
Then we can transform problem (PIM) into the penalty problem as follows: Based on the derivation in [17], it is clear thatJ δ is an exact penalty function. As a result, the optimal solution of problem (PIM) can be obtained by solving a sequence of problem (PIM δ ).

Optimization algorithm
For each δ > 0, problem (PIM δ ) 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 (PIM δ ), the time cost of estimating candidate parameters is quite expensive. As a result, we will develop a PPSO algorithm to solve problem (PIM δ ) for each δ > 0.
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 (PIM δ ) is a penalty problem with penalty parameter δ. As a result, we propose the following PPSO algorithm to solve problem (PIM δ ) for each δ. Here σ = (α 1 , . . . , α 10 , c 1 , . . . , c 10 , d 1 , . . . , d 10 , ζ) ∈ R 31 denotes the decision vector of problem (PIM δ ).

Algorithm 1
Step 1. Allocate n p slave processors. Initialize the total number of particles p size . Compute the number of particles on each processor n s = p size /n p . Initialize data and variables on the master processor.
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 τ (τ = 1, 2, . . . , n p ).
(ii) If σ τ,s κ (q) 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).
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.

Numerical simulations
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.
The above parameters are determined on the basis of a large number of experiments. By running Algorithm 2 the obtained optimal values of fractional orders and kinetic parameters are listed in Table 1.
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 comparison, 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 the 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 twostage 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 (PIM δ ). We perform 30 test runs in solving problem (PIM δ ) 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). Biomass (gL -1 ) Simulation using two-stage fractional system. Simulation using one-stage fractional system in [24]. Experimental data in [35]. Simulation using two-stage fractional system. Simulation using one-stage fractional system in [24]. Experimental data in [35]. Simulation using two-stage fractional system. Simulation using one-stage fractional system in [24]. Experimental data in [35]. Simulation using two-stage fractional system. Experimental data in [35].

Conclusions
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.