# Optimization for Nonlinear Uncertain Switched Stochastic Systems with Initial State Difference in Batch Culture Process

Authors: Jinlong Yuan; Jun Xie; Honglei Xu; Enmin Feng; Zhilong Xiu
Journal: Complexity
 (2019)

---

## 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 byKlebsiella 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.


---

## Body

## 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 byK. pneumoniae to 1,3-PD: batch culture [2], continuous culture [3–5], 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 [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,3-PD byK. 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 Section2, a nonlinear switched stochastic dynamical system is formulated to describe the process of glycerol biodissimilation to 1,3-PD byK. 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 instantstf1 and tf2 in different chemical reactions [31]. Accordingly, three different periods have been involved in the typical cell growth of the batch culture ([0,tf]⊂R+≔[0,+∞)) as follows:(i)
The development period (or the first period, denoted byD1≔[0,tf1]⊂R+)(ii)
The growth period (or the second period, denoted byD2≔[tf1,tf2]⊂R+)(iii)
The stabilization period (or the third period, denoted byD3≔[tf2,tf3=tf]⊂R+)Nomenclature (i)
A T denotes the transposition of the vector or matrix A. (ii)
I n ≔ { 1,2 , … , n }. (iii)
u t = ( t f 1 , t f 2 ) ∈ D ≔ [ a 1 , b 1 ] × [ a 2 , b 2 ] ⊂ R + 2 denotes switched instants. (iv)
0 = t f 0 < a 1 < b 1 < a 2 < b 2 < t f 3 = t f. (v)
Based on the sensitivity analysis [10], we take up≔((up1)T,(up2)T,(up3)T)T∈R15 as the system parameter vector to be optimized. (vi)
I N ≔ { 1,2 , … , N } denotes the set of experiment times in batch culture, where N is the total experiment times. (vii)
x 0 , l ∈ R 5 denotes the initial state of the first period of the lth experiment. (viii)
x i - 1 , l ( t f i - 1 ), i∈{2,3}, denote both the end period of the (i-1)th stage and the initial state of the ith period for the lth experiment. (ix)
x i , l ( t ) = ( x 1 i , l ( t ) , … , x 5 i , l ( t ) ) T ∈ R 5 denotes the state trajectory vector whose components are, respectively, the concentrations of biomass, extracellular glycerol, extracellular 1,3-PD, acetate, ethanol of the ith stage of the lth experiment at time t∈[0,tf].Based on [32], the nonlinear three-period dynamical system governing the batch culture can be described as follows:(1)x˙i,lt=hi,lxi,l,upi,t∈Dixi,ltfi-1=x0,l,i=1,xi-1,ltfi-1,i∈2,3.In system (1),(2)hi,lxi,l,upi=μi,lx1i,l,-upi2+μi,l0.0165x1i,l,upi3+41.2584μi,lx1i,l,upi4+4.541μi,lx1i,l,+upi53.046μi,lx1i,lT,where the specific cellular growth rate μi,l can be expressed as follows:(3)μi,l=upi1x2i,lx2i,l+0.28∏j=251-xji,lxj∗.Let up0≔[0.67,2.2,-2.69,-0.97,5.26]T [32]. The admissible range Uad of system parameters up1,up2 and up3 is defined as(4)Uad≔u_p,u¯p=∏j=15u_pj,u¯pj=∏j=15up0j-up0j2,up0j+up0j2⊂R5.Under anaerobic conditions at 37°C and pH=7.0, the concentrations of biomass, glycerol, and products are restricted in a certain range according to the practical production. x∗≔[0.0001,0.1,0,0,0]T and x∗≔[15,2039,939.5,1026,360.9]T [32], where xj∗≥0 and xj∗≥0,j∈I5 are the lower and upper bound of xji(t),i∈I3, respectively. So the admissible range S0 of xji(t),i∈I3, is(5)xit∈S0≔x∗,x∗=∏j=15xj∗,xj∗⊂R+5.


### 2.2. Stochastic System

Nomenclature (i)
( U a d , F , { F t } t ≥ 0 , P ) denotes a complete probability space, where F is a σ-algebra of subsets of Uad with a filtration {Ft}t≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while F0 contains all null subsets N of Uad). (ii)
E ( · ) denotes the expectation operator with respect to some probability measure P. (iii)
Fori∈I3,(6)Xi,lt≔X1i,lt,…,X5i,ltT∈S0,t∈tfi-1,tfi, which
denotes a stochastic process whose componentsXji,l(t),j∈I5, denote the scalar stochastic process on biomass, glycerol, 1,3-PD, acetic acid, and ethanol of the ith∈I3 stage in the lth∈IN experiment, respectively.The stochastic counterpart of system (1) can be rewritten in the matrix form(7)X˙i,lt=hi,lXi,lt,upi=Ai,lXti,l,t∈Di,Xi,ltfi-1=X0l,i=1,Xi-1,ltfi-1,i∈2,3,where Ai,l≔(aικi,l)5×5. It is clear to see that(8)a11i,l=μi,l,a21i,l=-upi2-μi,l0.0165,a31i,l=upi3+41.2584μi,l,a41i,l=upi4+4.541μi,l,a51i,l=upi5+3.046μi,l,aικi,l=0,ι,κ∈I5,κ≠1.Each system parameteraικi,l is stochastically perturbed as follows: (9)aικi,l←aικi,l+σικidWi,κt,where Wi(t)≔(Wi,1(t),…,Wi,5(t))T∈Ft, t∈[tfi-1,tfi], i∈I3; its components are Gaussian white noise of the ith∈I3 stage and the given diffusion matrix σi=(σικi)5×5 satisfies the following conditions:(10)σιιi>0,if1≤ι≤5,σικi≥0,ifι≠κ.The process of batch culture with perturbations can be formulated as the following nonlinear stochastic dynamical system:(11)dXi,lt=hi,lXi,lt,upidt+σiXi,ltdWit,t∈Di,i∈I3,Xi,ltfi-1=X0l,i=1,Xi-1,ltfi-1,i∈2,3,where(i)
σi(Xi,l(t))=σiXi,l(t),(ii)
hl(Xl(t),up) ≔ [h1,l(X1,l(t),up1), …, h3,l(X3,l(t),up3)]T,(iii)
E(X0l)=x0l.


### 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 functionshi,l(Xi,l(t),upi) and σi(Xi,l(t)),i∈I3,l∈IN, defined by Sections 2.1 and 2.2, we can easily carry out the proof of Property 1.Property 1.
The vector-valued functionshi,l(Xi,l(t),upi) and σi(Xi,l(t)),i∈I3,l∈IN, are measurable for t∈Di and Xi,l∈S0.Similarly to [26], we haveProperty 2.
For the given vector-valued functionshi,l(Xi,l(t),upi) and σi(Xi,l(t)),i∈I3,l∈IN, there exist positive constants K and K′ such that, for t∈Di, the following conditions hold:(i)
uniform Lipschitz condition(12)hi,lX^i,l,upi-hi,lXˇi,l,upi+σiX^i,l-σiXˇi,l≤KX^i,l-Xˇi,l,∀X^i,l,Xˇi,l∈S0;(ii)
growth continuous(13)hi,lX^i,l,upi+σiX^i,l≤K′1+X^i,l,∀X^i,l∈S0,where · is the Euclidean vector norm.By the proof in Property2, Theorem 5.2 in [33], and Theorem 5.4 in [34], we can get the following interesting properties.Property 3 (existence and uniqueness).
For∀l∈IN, given the vector-valued functions hi,l(Xi,l(t),upi) and σi(Xi,l(t)),i∈I3, system (11) has a unique solution, denoted by Xl(t;X0l,ut,up)≔(X1l(t;X0l,ut,up),…,X5l(t;X0l,ut,up))T, satisfying the initial condition X0l,l∈IN on [0,tf]. Furthermore, Xl(t;X0l,ut,up) is continuous with respect to (ut,up)∈D×Uad3 and satisfies the following integral equation:(14)Xlt;X0l,ut,up=X1,lt;X0l,up1=X0l+∫0th1,lX1,ls,up1ds+∫0tσ1X1,lsdW1s,t∈D1,X2,lt;X1,ltf1,up2=X1,ltf1+∫tf1th2,lX2,ls,up2ds+∫tf1tσ2X2,lsdW2s,t∈D2,X3,lt;X2,ltf2,up3=X2,ltf2+∫tf2th3,lX3,ls,up3ds+∫tf2tσ3X3,lsdW3s,t∈D3.Property 4 (Markov property and boundedness).
For∀l∈IN, the solution Xl(t;X0l,ut,up) is a Markov process on the interval [0,tf] whose initial probability distribution at t=0 is the distribution of X0l and Xl(t;X0l,ut,up) has continuous paths. Moreover (15)sup0≤t≤tfEXlt;X0l,ut,up2≤B1+EX0l2,where the constant B depends only on K, σ1, σ2, σ3, and tf.Property 5 (stochastic continuity).
Almost all realizations ofXl(t;X0l,ut,up) are continuous on [0,tf].


## 2.1. Deterministic System

To the best of our knowledge, in the actual batch culture process, there are two different switched instantstf1 and tf2 in different chemical reactions [31]. Accordingly, three different periods have been involved in the typical cell growth of the batch culture ([0,tf]⊂R+≔[0,+∞)) as follows:(i)
The development period (or the first period, denoted byD1≔[0,tf1]⊂R+)(ii)
The growth period (or the second period, denoted byD2≔[tf1,tf2]⊂R+)(iii)
The stabilization period (or the third period, denoted byD3≔[tf2,tf3=tf]⊂R+)Nomenclature (i)
A T denotes the transposition of the vector or matrix A. (ii)
I n ≔ { 1,2 , … , n }. (iii)
u t = ( t f 1 , t f 2 ) ∈ D ≔ [ a 1 , b 1 ] × [ a 2 , b 2 ] ⊂ R + 2 denotes switched instants. (iv)
0 = t f 0 < a 1 < b 1 < a 2 < b 2 < t f 3 = t f. (v)
Based on the sensitivity analysis [10], we take up≔((up1)T,(up2)T,(up3)T)T∈R15 as the system parameter vector to be optimized. (vi)
I N ≔ { 1,2 , … , N } denotes the set of experiment times in batch culture, where N is the total experiment times. (vii)
x 0 , l ∈ R 5 denotes the initial state of the first period of the lth experiment. (viii)
x i - 1 , l ( t f i - 1 ), i∈{2,3}, denote both the end period of the (i-1)th stage and the initial state of the ith period for the lth experiment. (ix)
x i , l ( t ) = ( x 1 i , l ( t ) , … , x 5 i , l ( t ) ) T ∈ R 5 denotes the state trajectory vector whose components are, respectively, the concentrations of biomass, extracellular glycerol, extracellular 1,3-PD, acetate, ethanol of the ith stage of the lth experiment at time t∈[0,tf].Based on [32], the nonlinear three-period dynamical system governing the batch culture can be described as follows:(1)x˙i,lt=hi,lxi,l,upi,t∈Dixi,ltfi-1=x0,l,i=1,xi-1,ltfi-1,i∈2,3.In system (1),(2)hi,lxi,l,upi=μi,lx1i,l,-upi2+μi,l0.0165x1i,l,upi3+41.2584μi,lx1i,l,upi4+4.541μi,lx1i,l,+upi53.046μi,lx1i,lT,where the specific cellular growth rate μi,l can be expressed as follows:(3)μi,l=upi1x2i,lx2i,l+0.28∏j=251-xji,lxj∗.Let up0≔[0.67,2.2,-2.69,-0.97,5.26]T [32]. The admissible range Uad of system parameters up1,up2 and up3 is defined as(4)Uad≔u_p,u¯p=∏j=15u_pj,u¯pj=∏j=15up0j-up0j2,up0j+up0j2⊂R5.Under anaerobic conditions at 37°C and pH=7.0, the concentrations of biomass, glycerol, and products are restricted in a certain range according to the practical production. x∗≔[0.0001,0.1,0,0,0]T and x∗≔[15,2039,939.5,1026,360.9]T [32], where xj∗≥0 and xj∗≥0,j∈I5 are the lower and upper bound of xji(t),i∈I3, respectively. So the admissible range S0 of xji(t),i∈I3, is(5)xit∈S0≔x∗,x∗=∏j=15xj∗,xj∗⊂R+5.


## 2.2. Stochastic System

Nomenclature (i)
( U a d , F , { F t } t ≥ 0 , P ) denotes a complete probability space, where F is a σ-algebra of subsets of Uad with a filtration {Ft}t≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while F0 contains all null subsets N of Uad). (ii)
E ( · ) denotes the expectation operator with respect to some probability measure P. (iii)
Fori∈I3,(6)Xi,lt≔X1i,lt,…,X5i,ltT∈S0,t∈tfi-1,tfi, which
denotes a stochastic process whose componentsXji,l(t),j∈I5, denote the scalar stochastic process on biomass, glycerol, 1,3-PD, acetic acid, and ethanol of the ith∈I3 stage in the lth∈IN experiment, respectively.The stochastic counterpart of system (1) can be rewritten in the matrix form(7)X˙i,lt=hi,lXi,lt,upi=Ai,lXti,l,t∈Di,Xi,ltfi-1=X0l,i=1,Xi-1,ltfi-1,i∈2,3,where Ai,l≔(aικi,l)5×5. It is clear to see that(8)a11i,l=μi,l,a21i,l=-upi2-μi,l0.0165,a31i,l=upi3+41.2584μi,l,a41i,l=upi4+4.541μi,l,a51i,l=upi5+3.046μi,l,aικi,l=0,ι,κ∈I5,κ≠1.Each system parameteraικi,l is stochastically perturbed as follows: (9)aικi,l←aικi,l+σικidWi,κt,where Wi(t)≔(Wi,1(t),…,Wi,5(t))T∈Ft, t∈[tfi-1,tfi], i∈I3; its components are Gaussian white noise of the ith∈I3 stage and the given diffusion matrix σi=(σικi)5×5 satisfies the following conditions:(10)σιιi>0,if1≤ι≤5,σικi≥0,ifι≠κ.The process of batch culture with perturbations can be formulated as the following nonlinear stochastic dynamical system:(11)dXi,lt=hi,lXi,lt,upidt+σiXi,ltdWit,t∈Di,i∈I3,Xi,ltfi-1=X0l,i=1,Xi-1,ltfi-1,i∈2,3,where(i)
σi(Xi,l(t))=σiXi,l(t),(ii)
hl(Xl(t),up) ≔ [h1,l(X1,l(t),up1), …, h3,l(X3,l(t),up3)]T,(iii)
E(X0l)=x0l.


## 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 functionshi,l(Xi,l(t),upi) and σi(Xi,l(t)),i∈I3,l∈IN, defined by Sections 2.1 and 2.2, we can easily carry out the proof of Property 1.Property 1.
The vector-valued functionshi,l(Xi,l(t),upi) and σi(Xi,l(t)),i∈I3,l∈IN, are measurable for t∈Di and Xi,l∈S0.Similarly to [26], we haveProperty 2.
For the given vector-valued functionshi,l(Xi,l(t),upi) and σi(Xi,l(t)),i∈I3,l∈IN, there exist positive constants K and K′ such that, for t∈Di, the following conditions hold:(i)
uniform Lipschitz condition(12)hi,lX^i,l,upi-hi,lXˇi,l,upi+σiX^i,l-σiXˇi,l≤KX^i,l-Xˇi,l,∀X^i,l,Xˇi,l∈S0;(ii)
growth continuous(13)hi,lX^i,l,upi+σiX^i,l≤K′1+X^i,l,∀X^i,l∈S0,where · is the Euclidean vector norm.By the proof in Property2, Theorem 5.2 in [33], and Theorem 5.4 in [34], we can get the following interesting properties.Property 3 (existence and uniqueness).
For∀l∈IN, given the vector-valued functions hi,l(Xi,l(t),upi) and σi(Xi,l(t)),i∈I3, system (11) has a unique solution, denoted by Xl(t;X0l,ut,up)≔(X1l(t;X0l,ut,up),…,X5l(t;X0l,ut,up))T, satisfying the initial condition X0l,l∈IN on [0,tf]. Furthermore, Xl(t;X0l,ut,up) is continuous with respect to (ut,up)∈D×Uad3 and satisfies the following integral equation:(14)Xlt;X0l,ut,up=X1,lt;X0l,up1=X0l+∫0th1,lX1,ls,up1ds+∫0tσ1X1,lsdW1s,t∈D1,X2,lt;X1,ltf1,up2=X1,ltf1+∫tf1th2,lX2,ls,up2ds+∫tf1tσ2X2,lsdW2s,t∈D2,X3,lt;X2,ltf2,up3=X2,ltf2+∫tf2th3,lX3,ls,up3ds+∫tf2tσ3X3,lsdW3s,t∈D3.Property 4 (Markov property and boundedness).
For∀l∈IN, the solution Xl(t;X0l,ut,up) is a Markov process on the interval [0,tf] whose initial probability distribution at t=0 is the distribution of X0l and Xl(t;X0l,ut,up) has continuous paths. Moreover (15)sup0≤t≤tfEXlt;X0l,ut,up2≤B1+EX0l2,where the constant B depends only on K, σ1, σ2, σ3, and tf.Property 5 (stochastic continuity).
Almost all realizations ofXl(t;X0l,ut,up) are continuous on [0,tf].


## 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 measuredN experimental data. Let yl(t) be the vector function fitted by experiment data. It denotes the concentration of every extracellular ingredient at time point t∈[0,tf]. The system identification problem is to choose an optimal vector (ut,up)∈D×Uad3 such that the expectation of distinction between stochastic process, denoted by Xl(t;X0l,ut,up), and yl(t),l∈IN, is minimized. Hence, the cost function can be defined by(16)Jut,up≔E∑l=1N∫0tfXlt;X0l,ut,up-ylt/yltdtN×tf.Now, for i∈I3,l∈IN, we are in a position to propose a system identification (SI) problem as follows:(17)minJut,ups.t.Xi,lt;Xi,ltfi-1,upi∈S0,∀t∈tfi-1,tfi,ut,up∈D×Uad3.The following theorem is to show the identifiability of theSI problem.Theorem 6.
TheSI problem admits an optimal solution.Proof.
According to Property3, we see that Xl(t;X0l,ut,up) is continuous with respect to (ut,up)∈D×Uad3, so the cost function J(ut,up) is continuous on (ut,up)∈D×Uad3. Moreover, D×Uad3 is a closed bounded set, which indicates that the optimal solution of theSI problem exists. Then, the proof is complete.


## 4. Time-Scaling Transformation and Approximate Subproblems

TheSI problem presents two major challenges for the existing numerical solution approaches:(i)
In theSI problem, the switched instants tf1 and tf2 are decision variables and need to be optimized. It is cumbersome to integrate the state and costate systems numerically when the switched instants tf1 and tf2 are decision variables.(ii)
The continuous state constraintsXi,l(t;Xi,l(tfi-1),upi)∈S0, ∀t∈[tfi-1,tfi], i∈I3, l∈IN 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 applytime-scaling transformation [35]—originally called the controlparameterization 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 newtime variables  θi as the duration of the ith subinterval in (tfi-1,tfi]. That is,(18)θi≔tfi-tfi-1,i∈I3.Then, the time-scaling transformation will be implemented by introducing a new time variables∈[0,3] and relating s to t through the equation(19)t=μ¯s=∑ȷ=0i-1θȷ+θis-i-1,s∈i-1,i,i∈I3,where μ is the so-calledtime-scaling function and θ0=0.Note that the time-scaling function is nondecreasing, continuous, and piecewise-linear. Moreover,(20)dμ¯sds=θi,s∈i-1,i,i∈I3.The new decision parametersθi, i∈I3, will satisfy(21)∑i=13θi=tf,0≤θi≤tf,i∈I3.Letθ≔[θ1,θ2,θ3]T. Any θ∈R3 satisfying (21) is called anadmissible duration vector. Now, since i-(i-1)=1, for i∈I3, (22)μ¯s~i=∑ȷ=1iθȷ=tfi.Therefore, the time-scaling function maps the fixed integer s=i to the switched instant t=tfi.LetX~i,l(s)≔Xi,l(μ¯(s)) and W~is≔Wiμ¯(s). If s∈[i-1,i], then μ¯(s)∈[tfi-1,tfi]. Substituting (19) and (20) into (11) gives(23)dX~i,ls≔h~i,lX~i,ls,upids+σ~iX~i,lsdW~is=dXi,lμ¯s=θihi,lX~i,l,upids+θiσiX~i,lsdW~is,s∈i-1,i,X~i,ls~i-1=X~0l,i=1,X~i-1,li-1,i∈2,3.Let(24)X~ls;X0l,θ,up=X~i,ls;X0l,θi,upi,i=1,X~i,ls;X~i-1,ls~i-1,θi,upi,i∈2,3.It denotes the solution of system (23) corresponding to the admissible pair (θ,up). Based on (16) and (21) it becomes(25)Jθ,up≔E∑l=1N∫03X~ls;X0l,θ,up-ylμ¯s/ylμ¯sds3N,(26)Φθ≔∑i=13θi-tf=0,0≤θi≤tf,i∈I3,where Φ:[0,tf]→R, is a given function. Then, for i∈I3,l∈IN, theSI problem becomes the transformation system identification (TSI) problem as follows:(27)min.Jθ,ups.t.X~i,ls;X~i-1,li-1,θi,upi∈S0,∀s∈i-1,i,θi,upi∈0,tf×Uad,Φθ=0.Note that the time-scaling transformation has replaced the variableswitched instants  tf1 and tf2 in the original approximate problem with conventional decision parameter vector θi,i∈I3, in the equivalent problem. Since i=0,1,2,3, in the equivalent problem, are fixed, this problem can be solved readily by any standard gradient-based optimization method or others.TheTSI problem is equivalent to the original problem. In fact, if (θi∗,upi∗),i∈I3, is a solution of the transformed problem, then the optimal switched instants for theSI problem are (28)tfi∗=∑ȷ=1iθȷ∗.


### 4.2. Approximate Problem

In view of the definition of semi-infinite programming problem [17, 36] and the continuous state inequality constraint in theTSI problem, we know that theTSI problem is a semi-infinite programming problem. An efficient algorithm for solving optimization problem of this type, which involve the so-calledconstraint transcription technique, is discussed in [36]. As a matter of fact, the essential difficulty lies in the judgement of the constraintX ~ i , l s ; X ~ i - 1 , l i - 1 , θ i , u p i = [X~1i,l(s;X~i-1,l(i-1),θi,upi),…,X~5i,l(s;X~i-1,l(i-1),θi,upi)]T∈S0, ∀s∈[i-1,i], i∈I3, l∈IN. To overcome these difficulties, let(29)gji,lX~i,ls;X~i-1,li-1,θi,upi≔X~ji,ls;X~i-1,li-1,θi,upi-xj∗,j∈I5,g5+ji,lX~i,ls;X~i-1,li-1,θi,upi≔xj∗-X~ji,ls;X~i-1,li-1,θi,upi,j∈I5.For i∈I3,l∈IN, the constraint (30)gji,lX~i,ls;X~i-1,li-1,θi,upi≤0,∀s∈i-1,i,is equivalently transcribed into(31)Gθi,upi,i,l≔∑j=110∫i-1imax⁡0,gji,lX~i,ls;X~i-1,li-1,θi,upids=0.However, G(θi,upi,i,l) is nonsmooth in (θi,upi). Consequently, the standard optimization routines would have difficulties with this type of equality constraints. Next, we replace (31) with(32)G~εθi,upi,i,l≔∑j=110∫i-1iφεgji,lX~i,ls;X~i-1,li-1,θi,upi,i,lds=0,where ε>0 and(33)φεgji,lX~i,ls;X~i-1,li-1,θi,upi,i,l≔0,ifgji,lX~i,ls;X~i-1,li-1,θi,upi<-ε,gji,lX~i,ls;X~i-1,li-1,θi,upi+ε24ε,if-ε≤gji,lX~i,ls;X~i-1,li-1,θi,upi≤ε,gji,lX~i,ls;X~i-1,li-1,θi,upi,ifgji,lX~i,ls;X~i-1,li-1,θi,upi>ε.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) with(34)G~ε,γθi,upi,i,l≔∑j=110∫i-1iφεgji,lX~i,ls;X~i-1,li-1,θi,upi,i,lds-γ≤0,where γ>0.For constructing the optimization algorithm, whenθi and upi are not feasible, we can move θi and upi towards the feasible region in the direction of -∂G~ε,γ(θi,upi,i,l)/∂θi and -∂G~ε,γ(θi,upi,i,l)/∂upi, respectively. In this paper, we develop a scheme for computing the gradients of constraint G~ε,γ(θi,upi,i,l) with respect to the parameters θi and upi, respectively.Theorem 7.
Givenl∈IN and i∈I3, for each ε>0,γ>0, the gradients of the constraint functionals G~ε,γ(θi,upi,i,l) with respect to upi and θi are(35)∂G~ε,γθi,upi,i,l∂upik≔∫i-1i∂HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l∂upikds,k∈I5,and(36)∂G~ε,γθi,upi,i,l∂θi≔∫i-1i∂HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l∂θids,where(37)HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l=∑j=110φεgji,lX~i,ls;X~i-1,li-1,θi,upi,i,l+λi,lsTh~i,lX~i,ls;X~i-1,li-1,θi,upi,upi+σ~iX~i,ls;X~i-1,li-1,θi,upiW~˙is, and(38)λi,ls=λ1i,ls,…,λ5i,lsT,is the solution of the costate system (39)(39)λ˙i,ls=-∂HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l∂X~i,ls;X~i-1,li-1,θi,upi,with the continuous boundary(40)λi,li=0,0,0,0,0T,i∈I3.Proof.
Letupi∈Uad be an arbitrary but fixed vector and ∀δk∈R,k∈I5. Define(41)upi,ϖk≔upi1,…,upik+ϖδk,…,upi5, where ϖ>0 is an arbitrarily small real number such that uk∗<upi(k)+ϖδj<uk∗, k∈I5. Therefore, G~ε,γ(θi,upi,ϖ,i,l) can be expressed as(42)G~ε,γθi,upi,ϖk,i,l≔γ+∑j=110∫i-1iφεgji,lX~i,ls;X~i-1,li-1,θi,upi,ϖ,i,lds+∫i-1iλi,lsT×hi,lX~i,ls;X~i-1,li-1,θi,upi,ϖ,upi,ϖ+σ~iX~i,ls;X~i-1,li-1,θi,upi,ϖW~˙is-X~˙i,ls;X~i-1,li-1,θi,upi,ϖds,where λi,l is yet arbitrary. Hence, it follows that(43)ΔG~ε,γθi,upi,i,l≔dG~ε,γθi,upi,ϖk,i,ldϖϖ=0=∂G~ε,γθi,upi,i,l∂upikδk=∫i-1i∂HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l∂X~i,ls;X~i-1,li-1,θi,upiΔX~i,ls;X~i-1,li-1,θi,upi+∂HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l∂upikδk-λi,lsTΔX~˙i,ls;X~i-1,li-1,θi,upids,where H(X~i,l(s;X~i-1,l(i-1),θi,upi),upi,λi,l(s),i,l) is defined as in (37). Integrating (43) by parts and combining (37), (38), (39), (40), and (42), we have(44)∂G~ε,γθi,upi,i,l∂upikδk≔∫i-1i∂HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l∂upikδkds. Since δk is arbitrary, conclusion (35) of the theorem follows.
Similarly to the proof of the above, we can prove(45)∂G~ε,γθi,upi,i,l∂θi≔∫i-1i∂HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l∂θids. The proof is complete.Then, theTSI problem can be approximated by the following problem:(46)TSIε,γ:min.Jθ,ups.t.G~ε,γθi,upi,i,l≤0,i∈I3,l∈IN,θi,upi∈t0,tf×Uad,i∈I3,Φθ-γ≤0.As is shown in [36], we can prove the following theorem.Theorem 8.
There exists aγ(ε)>0 such that for all γ, γ(ε)>γ>0, any feasible decision variables θε,γ∗ and upε,γ∗ to the TSIε,γ problem are also feasible decision variables to theTSI problem.Remark 9.
Theorem8 ensures that the corresponding γ(ε) for each ε in this sequence is finite.


## 4.1. Time-Scaling Transformation

Firstly, to handle the first challenge, one way of circumventing the difficulties caused by variable switched instants is to applytime-scaling transformation [35]—originally called the controlparameterization 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 newtime variables  θi as the duration of the ith subinterval in (tfi-1,tfi]. That is,(18)θi≔tfi-tfi-1,i∈I3.Then, the time-scaling transformation will be implemented by introducing a new time variables∈[0,3] and relating s to t through the equation(19)t=μ¯s=∑ȷ=0i-1θȷ+θis-i-1,s∈i-1,i,i∈I3,where μ is the so-calledtime-scaling function and θ0=0.Note that the time-scaling function is nondecreasing, continuous, and piecewise-linear. Moreover,(20)dμ¯sds=θi,s∈i-1,i,i∈I3.The new decision parametersθi, i∈I3, will satisfy(21)∑i=13θi=tf,0≤θi≤tf,i∈I3.Letθ≔[θ1,θ2,θ3]T. Any θ∈R3 satisfying (21) is called anadmissible duration vector. Now, since i-(i-1)=1, for i∈I3, (22)μ¯s~i=∑ȷ=1iθȷ=tfi.Therefore, the time-scaling function maps the fixed integer s=i to the switched instant t=tfi.LetX~i,l(s)≔Xi,l(μ¯(s)) and W~is≔Wiμ¯(s). If s∈[i-1,i], then μ¯(s)∈[tfi-1,tfi]. Substituting (19) and (20) into (11) gives(23)dX~i,ls≔h~i,lX~i,ls,upids+σ~iX~i,lsdW~is=dXi,lμ¯s=θihi,lX~i,l,upids+θiσiX~i,lsdW~is,s∈i-1,i,X~i,ls~i-1=X~0l,i=1,X~i-1,li-1,i∈2,3.Let(24)X~ls;X0l,θ,up=X~i,ls;X0l,θi,upi,i=1,X~i,ls;X~i-1,ls~i-1,θi,upi,i∈2,3.It denotes the solution of system (23) corresponding to the admissible pair (θ,up). Based on (16) and (21) it becomes(25)Jθ,up≔E∑l=1N∫03X~ls;X0l,θ,up-ylμ¯s/ylμ¯sds3N,(26)Φθ≔∑i=13θi-tf=0,0≤θi≤tf,i∈I3,where Φ:[0,tf]→R, is a given function. Then, for i∈I3,l∈IN, theSI problem becomes the transformation system identification (TSI) problem as follows:(27)min.Jθ,ups.t.X~i,ls;X~i-1,li-1,θi,upi∈S0,∀s∈i-1,i,θi,upi∈0,tf×Uad,Φθ=0.Note that the time-scaling transformation has replaced the variableswitched instants  tf1 and tf2 in the original approximate problem with conventional decision parameter vector θi,i∈I3, in the equivalent problem. Since i=0,1,2,3, in the equivalent problem, are fixed, this problem can be solved readily by any standard gradient-based optimization method or others.TheTSI problem is equivalent to the original problem. In fact, if (θi∗,upi∗),i∈I3, is a solution of the transformed problem, then the optimal switched instants for theSI problem are (28)tfi∗=∑ȷ=1iθȷ∗.


## 4.2. Approximate Problem

In view of the definition of semi-infinite programming problem [17, 36] and the continuous state inequality constraint in theTSI problem, we know that theTSI problem is a semi-infinite programming problem. An efficient algorithm for solving optimization problem of this type, which involve the so-calledconstraint transcription technique, is discussed in [36]. As a matter of fact, the essential difficulty lies in the judgement of the constraintX ~ i , l s ; X ~ i - 1 , l i - 1 , θ i , u p i = [X~1i,l(s;X~i-1,l(i-1),θi,upi),…,X~5i,l(s;X~i-1,l(i-1),θi,upi)]T∈S0, ∀s∈[i-1,i], i∈I3, l∈IN. To overcome these difficulties, let(29)gji,lX~i,ls;X~i-1,li-1,θi,upi≔X~ji,ls;X~i-1,li-1,θi,upi-xj∗,j∈I5,g5+ji,lX~i,ls;X~i-1,li-1,θi,upi≔xj∗-X~ji,ls;X~i-1,li-1,θi,upi,j∈I5.For i∈I3,l∈IN, the constraint (30)gji,lX~i,ls;X~i-1,li-1,θi,upi≤0,∀s∈i-1,i,is equivalently transcribed into(31)Gθi,upi,i,l≔∑j=110∫i-1imax⁡0,gji,lX~i,ls;X~i-1,li-1,θi,upids=0.However, G(θi,upi,i,l) is nonsmooth in (θi,upi). Consequently, the standard optimization routines would have difficulties with this type of equality constraints. Next, we replace (31) with(32)G~εθi,upi,i,l≔∑j=110∫i-1iφεgji,lX~i,ls;X~i-1,li-1,θi,upi,i,lds=0,where ε>0 and(33)φεgji,lX~i,ls;X~i-1,li-1,θi,upi,i,l≔0,ifgji,lX~i,ls;X~i-1,li-1,θi,upi<-ε,gji,lX~i,ls;X~i-1,li-1,θi,upi+ε24ε,if-ε≤gji,lX~i,ls;X~i-1,li-1,θi,upi≤ε,gji,lX~i,ls;X~i-1,li-1,θi,upi,ifgji,lX~i,ls;X~i-1,li-1,θi,upi>ε.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) with(34)G~ε,γθi,upi,i,l≔∑j=110∫i-1iφεgji,lX~i,ls;X~i-1,li-1,θi,upi,i,lds-γ≤0,where γ>0.For constructing the optimization algorithm, whenθi and upi are not feasible, we can move θi and upi towards the feasible region in the direction of -∂G~ε,γ(θi,upi,i,l)/∂θi and -∂G~ε,γ(θi,upi,i,l)/∂upi, respectively. In this paper, we develop a scheme for computing the gradients of constraint G~ε,γ(θi,upi,i,l) with respect to the parameters θi and upi, respectively.Theorem 7.
Givenl∈IN and i∈I3, for each ε>0,γ>0, the gradients of the constraint functionals G~ε,γ(θi,upi,i,l) with respect to upi and θi are(35)∂G~ε,γθi,upi,i,l∂upik≔∫i-1i∂HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l∂upikds,k∈I5,and(36)∂G~ε,γθi,upi,i,l∂θi≔∫i-1i∂HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l∂θids,where(37)HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l=∑j=110φεgji,lX~i,ls;X~i-1,li-1,θi,upi,i,l+λi,lsTh~i,lX~i,ls;X~i-1,li-1,θi,upi,upi+σ~iX~i,ls;X~i-1,li-1,θi,upiW~˙is, and(38)λi,ls=λ1i,ls,…,λ5i,lsT,is the solution of the costate system (39)(39)λ˙i,ls=-∂HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l∂X~i,ls;X~i-1,li-1,θi,upi,with the continuous boundary(40)λi,li=0,0,0,0,0T,i∈I3.Proof.
Letupi∈Uad be an arbitrary but fixed vector and ∀δk∈R,k∈I5. Define(41)upi,ϖk≔upi1,…,upik+ϖδk,…,upi5, where ϖ>0 is an arbitrarily small real number such that uk∗<upi(k)+ϖδj<uk∗, k∈I5. Therefore, G~ε,γ(θi,upi,ϖ,i,l) can be expressed as(42)G~ε,γθi,upi,ϖk,i,l≔γ+∑j=110∫i-1iφεgji,lX~i,ls;X~i-1,li-1,θi,upi,ϖ,i,lds+∫i-1iλi,lsT×hi,lX~i,ls;X~i-1,li-1,θi,upi,ϖ,upi,ϖ+σ~iX~i,ls;X~i-1,li-1,θi,upi,ϖW~˙is-X~˙i,ls;X~i-1,li-1,θi,upi,ϖds,where λi,l is yet arbitrary. Hence, it follows that(43)ΔG~ε,γθi,upi,i,l≔dG~ε,γθi,upi,ϖk,i,ldϖϖ=0=∂G~ε,γθi,upi,i,l∂upikδk=∫i-1i∂HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l∂X~i,ls;X~i-1,li-1,θi,upiΔX~i,ls;X~i-1,li-1,θi,upi+∂HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l∂upikδk-λi,lsTΔX~˙i,ls;X~i-1,li-1,θi,upids,where H(X~i,l(s;X~i-1,l(i-1),θi,upi),upi,λi,l(s),i,l) is defined as in (37). Integrating (43) by parts and combining (37), (38), (39), (40), and (42), we have(44)∂G~ε,γθi,upi,i,l∂upikδk≔∫i-1i∂HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l∂upikδkds. Since δk is arbitrary, conclusion (35) of the theorem follows.
Similarly to the proof of the above, we can prove(45)∂G~ε,γθi,upi,i,l∂θi≔∫i-1i∂HX~i,ls;X~i-1,li-1,θi,upi,upi,λi,ls,i,l∂θids. The proof is complete.Then, theTSI problem can be approximated by the following problem:(46)TSIε,γ:min.Jθ,ups.t.G~ε,γθi,upi,i,l≤0,i∈I3,l∈IN,θi,upi∈t0,tf×Uad,i∈I3,Φθ-γ≤0.As is shown in [36], we can prove the following theorem.Theorem 8.
There exists aγ(ε)>0 such that for all γ, γ(ε)>γ>0, any feasible decision variables θε,γ∗ and upε,γ∗ to the TSIε,γ problem are also feasible decision variables to theTSI problem.Remark 9.
Theorem8 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 theTSI 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 theTSI problem numerically is that there exists a huge number of numerical computations of differential equations. This makes solving theTSI 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 [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.Table 1
Definitions of variables in Section5.
Variable Representation Variable Representation CR(0) the maximal iteration taken from[0,1] F ( 0 ) the initial mutation factor taken from[0,1] N ~ the total number of genes in the population M T the maximal iteration r,rc,rl,rcr,rm random numbers R c , R m given constants h(kl,ȷ) the search direction ρ ( k l , ȷ ) the step-size selected by Armijo line search Fmax,Fmin the maximal and minimal mutation factors k∗ [ 0,0 , 0 , u _ p ] based on (5) k∗ [ t f , t f , t f , u ¯ p ]In addition, we solve the proposed stochastic identification problem using the following Stochastic Euler-Maruyama (EM) scheme. Let Δsi=s~i-s~i-1/Q for some positive integer Q,τιi=s~i-1+ιΔsi,ι∈IQ. Our numerical approximation to X~ji,l(τιi) will be denoted X~j,ιi,l. TheEM method takes the form(47)X~j,ιi,l=X~j,ι-1i,l+h~i,lX~j,ι-1i,l,upi+∑p=15σ~j,piX~j,ι-1i,ldWιi,p,i∈I3,j∈I5,l∈IN,ι∈IQ,where dWιi,p=Wi,p(τιi)-Wi,p(τι-1i) is an independent random variable of the form ΔsiN0,1 and N(0,1) 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 TSIε,γ problem. The main steps of the parallel MDE are as follows. For convenience, the optimization vector is denoted by (48)k≔θ1,θ2,θ3,upT≔θ1,θ2,θ3,up1T,up2T,up3TT.Remark 10.
The crossover factor is regulated by the following adaptive strategy:(49)CRl=r,ifrc≤RcCRl-1,otherwise,and the mutation factor F(l) is given by(50)Fl=Fmin+rlFmax,ifrm≤Rm,Fl-1,otherwise.Remark 11.
In Algorithm3, ε 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:tf≔7.75; system (23) is solved by using Euler method with a step size of 1/72000(h). In Algorithms 1–3, we select some parameters listed in Table 2 based on [36, 44]. The TSI problem is solved by using Algorithms 1–3.Table 2
The parameter value of Algorithms1–3.
Algorithms1-2 N ≔ 4, N~≔40; CR(0)≔0.5, F(0)≔0.4, Fmax≔0.9, Fmin≔0.1, Rc≔0.3, R m ≔ 0.3, MT≔600; Q=1000; σi,i=0.02, σi,p=0.0004, i≠p, i∈I5, p∈I5; x 0 1 = { 0.102,418.2609,0 , 0,0 , 0,0 , 0 }, x02={0.2025,441.337,0,0,0,0,0,0}, x 0 3 = { 0.173,402.9348,0 , 0,0 , 0,0 , 0 }, x04={0.2245,509.8913,0,0,0,0,0,0}  Algorithm3 ε ≔ 0.1, γ≔0.01, α≔0.1, β≔0.01, ε-≔1.0×10-7, γ-≔1.0×10-7Algorithm 1:Master processor (Processor 0).
Require: The experimental data such as x0l,yl(μ-(s)), s∈[0,3], l∈IN; Ensure: The value of kbest and Jbest; 1:Randomly generate N~ particle positions, denoted by k0,1,…,k0,N~. Based on (34), compute the constraint functions, denoted by G~ε,γ(k0,ȷ,l)≔(G~ε,γ(k0,ȷ,1,l), G~ε,γ(k0,ȷ,2,l),G~ε,γ(k0,ȷ,3,l))T, ȷ∈IN~, l∈IN; 2:for  each ȷ∈IN~  do 3:if  G~ε,γ(k0,ȷ,l)>0,l∈IN,  then 4:compute  h(k0,ȷ)≔∂G~ε,γk,l∂kk=k0,ȷ,  and move k0,ȷ according to k0,ȷ≔k0,ȷ-ρ(k0,ȷ)h(k0,ȷ), until G~ε,γ(k0,ȷ,l)≤0, with h(k0,ȷ) the search direction and ρ(k0,ȷ) 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 J(k0,ȷ), ȷ∈IN~ and let Jbest0≔J(kg), where kg≔argminȷ∈IN~{J(k0,ȷ)}; 8:Let kbest≔kg and broadcast (MPI_Broad), Jbest0,kbest and k0,ȷ, ȷ∈IN~, to all Processors; 9:Let l=1; 10:for  l=1; l<MT; l++  do 11:if  (l≥Mϱ and Jbestl-Mϱ-Jbestl<ϱ1)  then 12:return  kbest and Jbest; 13:end if 14:Let ȷ≔1; 15:for  ȷ=1; ȷ<N~; ȷ++  do 16:if  ȷ>N~  then 17:Let Jbestl≔J(kg),Jbest≔J(kg),kbest≔kg; 18:Finish the ȷth loop and start the (ȷ+1)th loop; 19:end if 20:Receive (MPI_Recv) the updated information of Processor ȷ during iteration l, including kl,ȷ and J(kl,ȷ); 21:if  J(kl,ȷ)<J(kg)  then 22:Let kg≔kl,ȷ and Broadcast (MPI_Broad) kl,ȷ and the updated kg to all slave processors; 23:else 24:Broadcast (MPI_Broad) kl,ȷ to all slave processors; 25:end if 26:end for 27:end forAlgorithm 2:Slave processors (Processor ȷ,ȷ=1,2,…,N~).
1:Let l=1; 2:for  l=1; l<MT; l++  do 3:if  l>MT  then 4:return Algorithm 1; 5:end if 6:Randomly choose p1 and p2 (ȷ≠p1≠p2) from {1,2,…,N~} and generate the trial vector kl,ȷ=kg+Flkl-1,p1-kl-1,p2,ifrcr≤CRlkl-1,ȷ,otherwise; 7:if  kl,ȷ violates boundary constraints  then 8:It is reflected back from the bound by kl,ȷ=2k∗-kl,ȷ,ifkl,ȷ∈-∞,k∗,2k∗-kl,ȷ,ifkl,ȷ∈k∗,+∞; 9:end if 10:Test the value of the constraint functions G~ε,γ(kl,ȷ,l) ≔ ( G ~ ε , γ ( k l , ȷ , 1 , l ), G~ε,γ(kl,ȷ,2,l), G~ε,γ(kl,ȷ,3,l))T, ȷ∈IN~, l∈IN. For each ȷ∈IN~, if G~ε,γ(kl,ȷ,l)>0,l∈IN, then, based on Theorem 7, compute  h(kl,ȷ)≔∂G~ε,γk,l∂kk=kl,ȷ, and move kl,ȷ according to kl,ȷ≔kl,ȷ-ρ(kl,ȷ)h(kl,ȷ), until G~ε,γ(kl,ȷ,l)≤0; 11:Solve system (23) by (47) and compute the cost function J(kl,ȷ),ȷ∈IN~ based on (25); 12:Send (MPI_Send) kl,ȷ and J(kl,ȷ),ȷ∈IN~ to Master processor; 13:end forAlgorithm 3:Combining Algorithms 1 and 2 with Theorem 8, the RIP problem can be solved by the algorithm.
1:Choose initial values of ε>0,γ>0; 2:Solve the TSIε,γ problem using Algorithm 1 to give kε,γ∗ ≔ (θ1∗,θ2∗,θ3∗,(up1∗)T,(up2∗)T,(up3∗)T)T; 3:Check feasibility of gji,l(X~i,l(s;X~i-1,l(s~i-1),θi∗,upi∗))≤0, j∈I10, l∈IN for ∀s∈[s~i-1,s~i], i∈I3; 4:if  kε,γ∗ is feasible  then 5:Let ε≔βε; 6:if  ε>ε-  then 7:Go to 2; 8:else 9:return  kε,γ∗; 10:end if 11:else 12:Let γ≔αγ; 13:if  γ>γ-  then 14:Go to 2; 15:else 16:return  kε,γ∗; 17:end if 18:end ifThe corresponding optimal system parametersup1∗,up2∗,up3∗ and time variables θ1∗,θ2∗,θ3∗ 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(51)eκl≔E∑j=13QX~κlsj;X0l,θ,up-yκlμsj∑j=13Qyκlsj,sj=j×1Q,j=1,2,…,3Q,κ=1,2,3.Table 3
The value ofupi∗,θi∗, i∈I3, tf1∗, and tf2∗.
i u p i ∗ 1 [0.620069, 2.03605, -1.86244, -0.67078, 7.50586] 2 [0.388043, 1.13123, -3.42823, -0.670187, 7.50937] 3 [0.352169, 1.65974, -3.28238, -0.669891, 7.51112] θ 1 ∗ = 2.6, θ2∗=2.4, θ3∗=2.75; By (18) and (21), tf1∗=2.6, tf2∗=5.Table 4
The relative errors under the optimal switched times and system parameters.
in the paper in [12] l e 1 l(%) e 2 l(%) e 3 l(%) e 1 l(%) e 2 l(%) e 3 l(%) l=1 4.20 4.21 16.93 41.87 10.21 19.11 l=2 5.73 6.99 4.89 38.14 13.31 27.87 l=3 4.09 4.40 5.73 43.25 12.02 15.94 l=4 4.80 3.24 4.43 23.47 16.26 24.16From Table4, 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 1–4, 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 1–4 are also confirmed that the nonlinear system with optimal kinetic parameters up1∗,up2∗,up3∗ and time variables θ1∗,θ2∗,θ3∗ can describe the batch culture process reasonably.Figure 1
The results of initial glycerol concentration of 418.2609mmol/L and biomass concentration of 0.102 g/L with 1000 sample trajectories.Figure 2
The results of initial glycerol concentration of 441.337mmol/L and biomass concentration of 0.2025g/L with 1000 sample trajectories.Figure 3
The results of initial glycerol concentration of 402.9348mmol/L and biomass concentration of 0.173g/L with 1000 sample trajectories.Figure 4
The results of initial glycerol concentration of 509.8913mmol/L and biomass concentration of 0.2245g/L with 1000 sample trajectories.


## 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 byK. 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.



---

*Source: 4979580-2019-02-03.xml*