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,3propanediol (1,3PD) 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 timescaling 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,3Propanediol (1,3PD) is an important chemical product with numerous applications in cosmetics, adhesives, lubricants, and medicines. In particular, 1,3PD 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,3PD: batch culture [2], continuous culture [3–5], and fedbatch culture [6]. The reason for the necessity of studying batch culture [7] is that 1,3PD yield, which is defined as the ratio between the formation of 1,3PD and the consumption of glycerol, is high in batch culture. The high concentration of glycerol leads to the high formation of 1,3PD and the low formation of byproduct 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 enzymecatalytic 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 timedelay system are researched; the study [14], where distributionally robust parameter identification of a timedelay 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 timedelay systems with statedependent 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 twostage nonlinear dynamical system are investigated; the study [20], where biobjective dynamic optimization of a nonlinear timedelay 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 [28–30]. 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,3PD 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 timescaling 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,3PD by K. pneumoniae and some significant properties are discussed. In Section 3, a system identification problem is proposed. In Section 4, timescaling 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,3PD, acetate, ethanol of the th stage of the th experiment at time .
Based on [32], the nonlinear threeperiod 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,3PD, 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 vectorvalued functions and , are measurable for and
Similarly to [26], we have
Property 2. For the given vectorvalued 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 vectorvalued 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. TimeScaling 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 timescaling 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. TimeScaling Transformation
Firstly, to handle the first challenge, one way of circumventing the difficulties caused by variable switched instants is to apply timescaling transformation [35]—originally called the control parameterization enhancing transform. The timescaling 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 timescaling transformation, we first introduce a new time variables as the duration of the subinterval in . That is,
Then, the timescaling transformation will be implemented by introducing a new time variable and relating to through the equationwhere is the socalled timescaling function and .
Note that the timescaling function is nondecreasing, continuous, and piecewiselinear. Moreover,
The new decision parameters , , will satisfy
Let Any satisfying (21) is called an admissible duration vector. Now, since , for , Therefore, the timescaling 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 timescaling 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 gradientbased 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 semiinfinite programming problem [17, 36] and the continuous state inequality constraint in the TSI problem, we know that the TSI problem is a semiinfinite programming problem. An efficient algorithm for solving optimization problem of this type, which involve the socalled 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 gradientbased 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 ringnetwork 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 [40–43], 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 EulerMaruyama (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 1–3, we select some parameters listed in Table 2 based on [36, 44]. The problem is solved by using Algorithms 1–3.
