Abstract

Based on the deterministic description of batch culture expressed in form of switched ordinary differential equations, we introduce a switched stochastic counterpart system with initial state difference together with uncertain switching instants and system parameters to model the process of glycerol biodissimilation to 1,3-propanediol (1,3-PD) induced by Klebsiella pneumoniae (K. pneumoniae). Important properties of the stochastic system are discussed. Our aim is to obtain the unified switched instants and system parameters under the condition of different initial states. To do this, we will formulate a system identification problem in which these uncertain switched instants and system parameters are regarded as decision variables to be chosen such that the relative error between experimental data and computational results is minimized. Such problem governed by the stochastic system is subject to continuous state inequality constraints and box constraints. By performing a time-scaling transformation as well as introducing the constraint transcription and local smoothing approximation techniques, we convert such problem into a sequence of approximation subproblems. Considering both the difficulty of finding analytical solutions and the complex nature of these subproblems, we develop a parallelized differential evolution (DE) algorithm to solve these approximation subproblems. From an extensive simulation, we show that the obtained optimal switched instants and system parameters are satisfactory with initial state difference.

1. Introduction

1,3-Propanediol (1,3-PD) is an important chemical product with numerous applications in cosmetics, adhesives, lubricants, and medicines. In particular, 1,3-PD has recently been used as a monomer to synthesize a new type of polyester called polyurethanes [1]. There are three common methods of the bioconversion of glycerol by K. pneumoniae to 1,3-PD: batch culture [2], continuous culture [35], and fed-batch culture [6]. The reason for the necessity of studying batch culture [7] is that 1,3-PD yield, which is defined as the ratio between the formation of 1,3-PD and the consumption of glycerol, is high in batch culture. The high concentration of glycerol leads to the high formation of 1,3-PD and the low formation of by-product in batch culture of glycerol [8]. In particular, note that the high target product yield in batch culture is only applicable to the bioconversion of glycerol. Not all the batch culture is of the high target product yield. Therefore, batch culture, in which the bacteria and substrate will be added to the bioreactor only at the beginning process, is now attracting significant interest in many research areas. Relevant literature includes [9], where identification and robustness analysis of nonlinear multistage enzyme-catalytic dynamical system are researched; the study [10], where sensitivity analysis and identification of kinetic parameters are investigated; the study [11], where strong stability of a nonlinear multistage dynamic system is considered; the study [12], where robust identification of enzymatic nonlinear dynamical systems is studied; the study [13], where modelling and parameter identification for a nonlinear time-delay system are researched; the study [14], where distributionally robust parameter identification of a time-delay dynamical system with stochastic measurements is carried out (many wireless network systems and hybrid dynamical systems [15, 16] have time delays and switchings); the study [17], where dynamic optimization for switched time-delay systems with state-dependent switching conditions is studied; the study [18], where pathway identification using parallel optimization for a nonlinear hybrid system is considered; the study [19], where optimality condition and optimal control for a two-stage nonlinear dynamical system are investigated; the study [20], where bi-objective dynamic optimization of a nonlinear time-delay system is carried out; the study [21], where robust biobjective optimal control is studied. To the best of the authors’ knowledge and by surveying mentioned literatures, stochastic influences are not taken into consideration in batch culture.

Biotechnical treatment of microorganisms is commonly described by deterministic systems in form of nonlinear ordinary differential equations [22] to avoid expensive experiments. This description includes an idealization of the technical system component and a qualitative characterization of the biological part. In fact, microbial fermentation cannot commonly avoid stochastic influences reflected on the uncertainty of certain parameters [23]. Since this random phenomenon reveals different patterns and new features, such randomnness will degrade the role of deterministic systems in the research [24]. Therefore, research interest on stochastic dynamic systems has been growing in recent years [25]. Biological phenomena have the dynamical behaviours that are intrinsically erratic, and they are concisely described by a stochastic model, rather than by a deterministic one. In [26, 27], the randomness is introduced by the parameter perturbation and it becomes a standard technique in stochastic population modelling. However, parameter uncertainty is still ignored in their model.

Parameter uncertainty is a key issue in practice because it is difficult (if not impossible) to determine the exact values of many parameters in the dynamic equations describing microbial conversion [2830]. In this paper, a switched stochastic counterpart with uncertain switched instants and system parameters will be taken to replace the commonly used deterministic description of batch culture with ordinary differential equations to describe the process of glycerol biodissimilation to 1,3-PD by K. pneumoniae. In order to estimate the uncertain switched instants and system parameters, we present a system identification problem governed by the stochastic system and subject to continuous state inequality constraints and box constraints to minimize the relative error between experimental data and computational results. In order to handle continuous state inequality constraints, such problem is converted into a sequence of approximation subproblems and one solves them by using time-scaling transformation, constraint transcription, and local smoothing approximation techniques. Since it is a very complicated task to solve these subproblems, we develop a parallelized differential evolution (DE) algorithm to solve these approximation subproblems. Our contribution of this paper includes (1) to propose a new stochastic switched models for batch fermentation to describe the process’s randomness characteristics, (2) to establish a optimal identification procedure to identify the parameters, and (3) to design a parallel DE algorithm to save the computational time efficiently. It is observed that the obtained optimal switched instants and system parameters are satisfactory with initial state difference via numerical simulations.

The rest of this paper is organized as follows. In Section 2, a nonlinear switched stochastic dynamical system is formulated to describe the process of glycerol biodissimilation to 1,3-PD by K. pneumoniae and some significant properties are discussed. In Section 3, a system identification problem is proposed. In Section 4, time-scaling transformation is performed and approximate subproblems are presented. In Section 5, an optimization algorithm is constructed to solve these subproblems. In Section 6, the numerical results are clarified. Then, conclusion remarks are presented in Section 7.

2. Nonlinear Switched Stochastic System and Its Properties

2.1. Deterministic System

To the best of our knowledge, in the actual batch culture process, there are two different switched instants and in different chemical reactions [31]. Accordingly, three different periods have been involved in the typical cell growth of the batch culture ( as follows:(i)The development period (or the first period, denoted by )(ii)The growth period (or the second period, denoted by )(iii)The stabilization period (or the third period, denoted by )

Nomenclature(i) denotes the transposition of the vector or matrix .(ii).(iii) denotes switched instants.(iv).(v)Based on the sensitivity analysis [10], we take as the system parameter vector to be optimized.(vi) denotes the set of experiment times in batch culture, where is the total experiment times.(vii) denotes the initial state of the first period of the th experiment.(viii), , denote both the end period of the th stage and the initial state of the th period for the th experiment.(ix) denotes the state trajectory vector whose components are, respectively, the concentrations of biomass, extracellular glycerol, extracellular 1,3-PD, acetate, ethanol of the th stage of the th experiment at time .

Based on [32], the nonlinear three-period dynamical system governing the batch culture can be described as follows:In system (1),where the specific cellular growth rate can be expressed as follows:Let [32]. The admissible range of system parameters and is defined asUnder anaerobic conditions at and pH=7.0, the concentrations of biomass, glycerol, and products are restricted in a certain range according to the practical production. and [32], where and are the lower and upper bound of , respectively. So the admissible range of , is

2.2. Stochastic System

Nomenclature(i) denotes a complete probability space, where is a -algebra of subsets of with a filtration satisfying the usual conditions (i.e., it is increasing and right continuous while contains all null subsets of ).(ii) denotes the expectation operator with respect to some probability measure (iii)For , which denotes a stochastic process whose components , denote the scalar stochastic process on biomass, glycerol, 1,3-PD, acetic acid, and ethanol of the stage in the experiment, respectively.

The stochastic counterpart of system (1) can be rewritten in the matrix formwhere . It is clear to see that

Each system parameter is stochastically perturbed as follows: where , , ; its components are Gaussian white noise of the stage and the given diffusion matrix satisfies the following conditions:

The process of batch culture with perturbations can be formulated as the following nonlinear stochastic dynamical system:where(i),(ii) ≔ [, …, ,(iii)

2.3. Properties of Solutions to Stochastic System

In this section, we will discuss some properties of the solutions to the stochastic system. According to the definition of the functions and , defined by Sections 2.1 and 2.2, we can easily carry out the proof of Property 1.

Property 1. The vector-valued functions and , are measurable for and

Similarly to [26], we have

Property 2. For the given vector-valued functions and , there exist positive constants and such that, for , the following conditions hold:(i)uniform Lipschitz condition(ii)growth continuous where is the Euclidean vector norm.

By the proof in Property 2, Theorem 5.2 in [33], and Theorem 5.4 in [34], we can get the following interesting properties.

Property 3 (existence and uniqueness). For , given the vector-valued functions and , system (11) has a unique solution, denoted by , satisfying the initial condition on . Furthermore, is continuous with respect to and satisfies the following integral equation:

Property 4 (Markov property and boundedness). For , the solution is a Markov process on the interval whose initial probability distribution at is the distribution of and has continuous paths. Moreover where the constant depends only on , , , , and

Property 5 (stochastic continuity). Almost all realizations of are continuous on

3. System Identification Problem

The system identification problem governed by a stochastic system is generally to adjust the values of switched instants and system parameters so that the discrepancy between predicted and observed system output is as small as possible. The purpose of this section is to establish a sysem identification problem governed by system (11) in batch culture.

In the process of batch culture, we have measured experimental data. Let be the vector function fitted by experiment data. It denotes the concentration of every extracellular ingredient at time point The system identification problem is to choose an optimal vector such that the expectation of distinction between stochastic process, denoted by , and , is minimized. Hence, the cost function can be defined byNow, for , we are in a position to propose a system identification (SI) problem as follows:The following theorem is to show the identifiability of the SI problem.

Theorem 6. The SI problem admits an optimal solution.

Proof. According to Property 3, we see that is continuous with respect to , so the cost function is continuous on . Moreover, is a closed bounded set, which indicates that the optimal solution of the SI problem exists. Then, the proof is complete.

4. Time-Scaling Transformation and Approximate Subproblems

The SI problem presents two major challenges for the existing numerical solution approaches:(i)In the SI problem, the switched instants and are decision variables and need to be optimized. It is cumbersome to integrate the state and costate systems numerically when the switched instants and are decision variables.(ii)The continuous state constraints , , , are not easy to be satisfied in practice.

For these reasons, these conventional numerical optimization algorithms struggle to handle these challenges. In the section, we use the time-scaling transformation to deal with the variable time points together with the constraint transformation and the local smoothing approximate technique to handle the continuous state constraints.

4.1. Time-Scaling Transformation

Firstly, to handle the first challenge, one way of circumventing the difficulties caused by variable switched instants is to apply time-scaling transformation [35]—originally called the control parameterization enhancing transform. The time-scaling transformation works by mapping the variable switched instants to fixed points in a new time horizon, thus yielding a new optimization problem with the fixed switched instants. To apply the time-scaling transformation, we first introduce a new time variables   as the duration of the subinterval in . That is,

Then, the time-scaling transformation will be implemented by introducing a new time variable and relating to through the equationwhere is the so-called time-scaling function and .

Note that the time-scaling function is nondecreasing, continuous, and piecewise-linear. Moreover,

The new decision parameters , , will satisfy

Let Any satisfying (21) is called an admissible duration vector. Now, since , for , Therefore, the time-scaling function maps the fixed integer to the switched instant

Let and If , then . Substituting (19) and (20) into (11) gives

LetIt denotes the solution of system (23) corresponding to the admissible pair Based on (16) and (21) it becomeswhere , is a given function. Then, for , the SI problem becomes the transformation system identification (TSI) problem as follows:

Note that the time-scaling transformation has replaced the variable switched instants   and in the original approximate problem with conventional decision parameter vector , in the equivalent problem. Since , in the equivalent problem, are fixed, this problem can be solved readily by any standard gradient-based optimization method or others.

The TSI problem is equivalent to the original problem. In fact, if , is a solution of the transformed problem, then the optimal switched instants for the SI problem are

4.2. Approximate Problem

In view of the definition of semi-infinite programming problem [17, 36] and the continuous state inequality constraint in the TSI problem, we know that the TSI problem is a semi-infinite programming problem. An efficient algorithm for solving optimization problem of this type, which involve the so-called constraint transcription technique, is discussed in [36]. As a matter of fact, the essential difficulty lies in the judgement of the constraint

= , , , To overcome these difficulties, letFor , the constraint is equivalently transcribed intoHowever, is nonsmooth in . Consequently, the standard optimization routines would have difficulties with this type of equality constraints. Next, we replace (31) withwhere and

Note that the equality constraints specified in (32) do not satisfy the constraint qualification. Thus it is not advisable to compute it as such numerically. For this reason, we replace (32) withwhere

For constructing the optimization algorithm, when and are not feasible, we can move and towards the feasible region in the direction of and , respectively. In this paper, we develop a scheme for computing the gradients of constraint with respect to the parameters and , respectively.

Theorem 7. Given and , for each , the gradients of the constraint functionals with respect to and areandwhere andis the solution of the costate system (39)with the continuous boundary

Proof. Let be an arbitrary but fixed vector and . Define where is an arbitrarily small real number such that , Therefore, can be expressed aswhere is yet arbitrary. Hence, it follows thatwhere is defined as in (37). Integrating (43) by parts and combining (37), (38), (39), (40), and (42), we have Since is arbitrary, conclusion (35) of the theorem follows.
Similarly to the proof of the above, we can prove The proof is complete.

Then, the TSI problem can be approximated by the following problem:

As is shown in [36], we can prove the following theorem.

Theorem 8. There exists a such that for all , , any feasible decision variables and to the problem are also feasible decision variables to the TSI problem.

Remark 9. Theorem 8 ensures that the corresponding for each in this sequence is finite.

5. Parallel Differential Evolution Algorithm

Since it is usually not possible to derive an analytical solution to system (23), numerical approaches are indispensable, especially for biological fermentation. To solve the TSI problem numerically, various optimization routes, such as gradient-based algorithms [37], can be used. Nonetheless, all those techniques are only designed to find local optima. Motivated by the mechanisms of natural selection, differential evolution (DE) was first proposed by Storn and Price [38] in 1997. DE is a recent optimization technique and an exceptionally simple and easy method used for the evolution strategy. It is a significantly fast and robust numerical optimization method and it is more likely to find the true global optimum.

Furthermore, one of main hurdles with solving the TSI problem numerically is that there exists a huge number of numerical computations of differential equations. This makes solving the TSI problem intolerable by a serial IDE algorithm. To improve computational efficiency, it is natural to construct a parallel optimization algorithm. In 2004, Tasoulis et al. [39] proposed a parallel DE algorithm using a ring-network topology, which can improve both the speed and the performance of the method. This point of view, which has been observed for a wide variety of experiments [4043], is being gradually accepted by experts in the field of parallel DE.

However, what we need to solve is an optimization problem with both box constraints and continuous state constraints, to which DE cannot be applied directly [41]. To handle such constraints, we introduce the gradients of the constraint functions into our algorithm (see Theorem 7). Definitions of variables are shown in Table 1.

In addition, we solve the proposed stochastic identification problem using the following Stochastic Euler-Maruyama (EM) scheme. Let for some positive integer Our numerical approximation to will be denoted . The EM method takes the formwhere is an independent random variable of the form and denotes a normally distributed random variable with zero mean and unit variance.

With this in mind, given , combining Theorem 8, we propose a parallel modified differential evolution (MDE) algorithm for solving the problem. The main steps of the parallel MDE are as follows. For convenience, the optimization vector is denoted by

Remark 10. The crossover factor is regulated by the following adaptive strategy: and the mutation factor is given by

Remark 11. In Algorithm 3, is a parameter controlling the accuracy of the smoothing approximation; is a parameter controlling the feasibility of constraint (5); and are two predefined positive parameters ensuring the termination of Algorithm 3; parameters and must be chosen less than 1.

6. Numerical Results

The below parameters are derived empirically after numerous experiments: ; system (23) is solved by using Euler method with a step size of (h). In Algorithms 13, we select some parameters listed in Table 2 based on [36, 44]. The problem is solved by using Algorithms 13.

Require:
The experimental data such as , , ;
Ensure:
The value of and ;
1: Randomly generate particle positions, denoted by . Based on (34), compute the
constraint functions, denoted by , , , ;
2: for  each   do
3: if  ,  then
4: compute  ,  and move according to ,
until , with the search direction and the step-size selected by Armijo line search;
5: end if
6: end for
7: Solve system (23) by (47). Based on (25), compute the cost function, denoted by ,
and let , where ;
8: Let and broadcast , and , , to all Processors;
9: Let ;
10: for  ; ;   do
11: if  ( and )  then
12: return   and ;
13: end if
14: Let ;
15: for  ; ;   do
16: if    then
17: Let ;
18: Finish the loop and start the loop;
19: end if
20: Receive () the updated information of Processor during iteration , including and ;
21: if    then
22: Let and Broadcast () and the updated to all slave processors;
23: else
24: Broadcast () to all slave processors;
25: end if
26: end for
27: end for
1: Let ;
2: for  ; ;   do
3: if    then
4: return Algorithm 1;
5: end if
6: Randomly choose and () from and generate the trial vector
  
7: if   violates boundary constraints  then
8: It is reflected back from the bound by
  
9: end if
10: Test the value of the constraint functions
, , , , For each ,
if , then, based on Theorem 7, compute  ,
and move according to , until ;
11: Solve system (23) by (47) and compute the cost function based on (25);
12: Send () and to Master processor;
13: end for
1: Choose initial values of ;
2: Solve the problem using Algorithm 1 to give ;
3: Check feasibility of , , for , ;
4: if   is feasible  then
5: Let ;
6: if    then
7: Go to 2;
8: else
9: return  ;
10: end if
11: else
12: Let ;
13: if    then
14: Go to 2;
15: else
16: return  ;
17: end if
18: end if

The corresponding optimal system parameters and time variables are shown in Table 3. For comparison, the relative errors between the computational values and the experimental data in this work and the ones in [12] are listed in Table 4, in which the relative errors are defined as

From Table 4, we can see that the relative errors are cut down greatly in comparison with the ones in [12]. The errors between experimental data and expectation of sample trajectory with different initial glycerol and biomass concentrations are displayed in Figures 14, where the horizontal axes represent time; the left vertical axes represent the concentrations of biomass, the center vertical axes represent the concentrations of glycerol, while the right vertical axes apply for 1,3-PD; the scattergrams denote experimental data, the blue solid lines represent sample trajectory, and the red solid lines display expectation of sample trajectory. The curves in Figures 14 are also confirmed that the nonlinear system with optimal kinetic parameters and time variables can describe the batch culture process reasonably.

7. Conclusion

In this paper, we introduce a switched stochastic counterpart with uncertain switched instants and system parameters to replace the commonly used deterministic description of glycerol biodissimilation to 1,3-PD by K. pneumoniae in form of ordinary differential equations. Some important properties of the stochastic system are discussed. Our aim is to identify the uncertain switched instants and system parameters under condition of different initial state. For this, taking the relative error between experimental data and computational results as the cost function, we propose a system identification problem governed by the stochastic system and subject to some constraints. In consideration of both the difficulty of finding analytical solutions and the complexity of solution procedure, a parallelized differential evolution algorithm is developed to solve the identification problem. An illustrative numerical example shows the appropriateness of the optimal switched instants and system parameters with initial state difference.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant nos. 11771008 and 61773086), the ARC discovery grant (Grant no. DP160102819), the National Science Foundation for the Youth of China (Grant nos. 11401073, 11501574, and 11701063), the National Science Foundation for the Tianyuan of China (Grant no. 11626053), the China Postdoctoral Science Foundation (Grant No. 11626053), the Natural Science Foundation of Shandong Province in China (Grant no. ZR2017MA005), the Fundamental Research Funds for the Centrol Universities (Grant nos. DUT16LK07 and DLMU3132018215), and Tianjin Thousand Talents Program and Xinghai Project of Dalian Maritime University.