Modelling and optimal control for a fed-batch fermentation process

Modelling and optimal control for a fed-batch fermentation process

Applied Mathematical Modelling 37 (2013) 695–706 Contents lists available at SciVerse ScienceDirect Applied Mathematical Modelling journal homepage:...

318KB Sizes 2 Downloads 37 Views

Applied Mathematical Modelling 37 (2013) 695–706

Contents lists available at SciVerse ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Modelling and optimal control for a fed-batch fermentation process Chongyang Liu a,⇑, Zhaohua Gong a, Bangyu Shen b, Enmin Feng c a b c

School of Mathematics and Information Science, Shandong Institute of Business and Technology, Yantai 264005, China School of Mathematical Science, Huaiyin Normal University, Huai’an 223300, China School of Mathematical Sciences, Dalian University of Technology, Dalian 116024, China

a r t i c l e

i n f o

Article history: Received 2 May 2011 Received in revised form 15 February 2012 Accepted 29 February 2012 Available online 13 March 2012 Keywords: Nonlinear dynamical system Optimal control Control parametrization enhancing transform Constraint transcription Fed-batch fermentation

a b s t r a c t The main control goal in fed-batch fermentation is to get a high concentration of production. In this paper, a new nonlinear dynamical system, in which the feed rate of glycerol is the control function and the switching instants between the batch and feed processes are the variables, is proposed to formulate the fed-batch fermentation process of glycerol bioconversion to 1,3-propanediol (1,3-PD). To maximize the concentration of 1,3-PD at the terminal time, an optimal control model involving the nonlinear system and subject to the continuous state constraints and the control constraint is presented. To seek the optimal solution of the constrained optimal control problem, the control parametrization enhancing transform together with the constraint transcription technique is first used to convert the constrained optimal control problem into a sequence of mathematical programming problems. An improved Particle Swarm Optimization is subsequently constructed to solve the resultant mathematical programming problem. Numerical results show that, by employing the obtained optimal strategy, 1,3-PD concentration at the terminal time can be increased considerably. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction Fed-batch cultivation technique is a mode of bioreactors operation that provides distinct advantages over the other operation modes and has a widespread industrial application [1]. Concentration of substrate in cultivation medium of fed-batch process can be externally manipulated by using the appropriate feed-rate profiles. Optimal control of fed-batch fermentation process is one of the relevant problems in biotechnology. Consequently, optimal control of fed-batch fermentation processes has been a topic of research for many years [2,3]. 1,3-Propanediol (1,3-PD) is a valuable chemical intermediate that is suitable as a monomer for polycondensations to produce polyesters, polyethers and polyurethanes [4]. The microbial conversion of glycerol to 1,3-PD is particularly attractive in that the process is relatively easy and does not generate toxic byproducts. Dissimilation of glycerol to 1,3-PD by Klebsiella pneumoniae (K. pneumoniae) has been widely investigated since the 1980s due to its high productivity [5]. The critical concentrations of glycerol, 1,3-PD and other byproducts in the fermentation are determined [6]. The accumulation of 3-hydroxypropionaldehyde is investigated [7]. High key enzymes activities are obtained on the basis of optimal fermentation conditions [8]. During the bioconversion of glycerol to 1,3-PD, the most efficient cultivation method is a fed-batch operation which corrects pH by alkali addition for glycerol supply [9]. The fed-batch of bioconversion glycerol to 1,3-PD begins with the cells being grown under the batch culture for some time, usually until close to the end of the exponential growth phase (i.e., a period in which the number of new bacteria appearing per unit time is proportional to the present population). At this ⇑ Corresponding author. E-mail addresses: [email protected] (C. Liu), [email protected] (Z. Gong). 0307-904X/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.apm.2012.02.044

696

C. Liu et al. / Applied Mathematical Modelling 37 (2013) 695–706

point glycerol and alkali are continuously fed into the reactor at some rates, without the removal of culture fluid. This helps to maintain a suitable environment for the cells’ growth. At the end of the feed process, another batch phase starts again. The above processes are repeated until the end of the final batch phase. Hence, it is decisive for improving the productivity of 1,3PD to optimize the feed rates and the switching instants between the batch and feed processes in fed-batch fermentation. Determining the mathematical model can often be a very vital task since using non-accurate model in the calculations of the optimal feed rate may lead to undesirable results. Initially, the batch model is extrapolated to fed-batch cultivation by incorporating the dilution factors. Unstructured and nonsegregational models with specific rates of cell growth, metabolite production, and substrate consumption have been used to model fed-batch fermentation [10–12]. The models have been used for optimal control studies by a number of researchers [13,14]. Recently, nonlinear dynamical systems have been extensively investigated to formulate the fermentation process. Based on an assumption that the feed of glycerol only occurs at impulsive instants, nonlinear impulsive systems have been extensively investigated to formulate the fermentation process [15–18]. Nonetheless, since the feed rate of glycerol is finite, it is not reasonable to describe the actual fed-batch fermentation by the impulsive dynamical systems. Moreover, these studies concentrate on deducing the optimality conditions of the impulsive optimal control problem. In contrast, taking the feed of glycerol as a time-continuous process, we have proposed a nonlinear multistage dynamical system to formulate the fed-batch process [19]. For this system, the parameter identification problem was then investigated. Numerical simulations indicated this multistage dynamical system could describe the microbial fed-batch fermentation better compared with the impulsive dynamical system [20]. Furthermore, taking the feed rate of glycerol as the control function, a multistage optimal control model was also presented. Some computational approaches were also developed to seek the optimal solution of the multistage optimal control problem [19,21]. Numerical results showed that, by employing obtained optimal strategies, the concentration of 1,3-PD at the terminal time can be increased considerably. However, in all of the above models, the switching instants between the batch and feed processes are decided a priori. In this paper, considering the feed process as a time-continuous process, we propose a new nonlinear dynamical system, in which the feed rate of glycerol is taken as the control function and the switching instants between the batch and feed processes as the variables, to formulate the fed-batch fermentation process. Especially, optimal control of such switched system have been an active research area over the past decade, see, for example [22–25]. Nevertheless, to our knowledge, the switched systems with continuous state constraints have rarely been considered. Optimal control in fed-batch fermentation process is to maximize the final 1,3-PD productivity. For comparison, the whole fermentation period is fixed. Therefore, the 1,3-PD concentration at the terminal time is taken as the performance index in this work. By the way, many studies have considered the same performance index in optimal control of the fed-batch fermentation process [21,26,27]. Then, an optimal control problem subject to the continuous state constraints and the control constraint is presented. Incidentally, constrained optimal control problems have been extensively studied in the literature. Many interesting theoretical results can be found in [28]. For numerical computation, several successful families of algorithms have been developed, see, for example [29–32]. In particular, the control parametrization enhancing transform [31] has been used extensively in [33–35]. In this paper, to seek the optimal feed rate as well as the optimal switching instants in the constrained optimal control problem, the control parametrization enhancing transform is first used to approximate the constrained optimal control problem. The constraint transcription and smoothing techniques [36] are then applied to dealing with the continuous state constraints. In the end, an improved Particle Swarm Optimization (PSO) algorithm is developed to solve the resultant mathematical programming problems. Numerical results show that, by employing the obtained optimal strategy, the concentration of 1,3-PD at the terminal time can be increased considerably compared with previous results. This paper is organized as follows. In the next section, Section 2, the mathematical model of the fed-batch fermentation process is presented. Section 3 gives the constrained optimal control problem. Section 4 develops a computational approach to solve the constrained optimal control problem, while Section 5 illustrates the numerical results. Finally, conclusions are provided in Section 6. 2. Mathematical models in fed-batch fermentation The fed-batch fermentation begins with a batch culture, then batch-fed glycerol and alkali are poured into the reactor in order to provide nutrition and maintain a suitable environment for the cells’ growth. According to the actual fermentation process, we assume that (H1) The concentrations of reactants are uniform in the reactor, while time delay and nonuniform space distribution are ignored. (H2) The feeding media includes only fixed concentrations of glycerol and alkali. Moreover, the velocity ratio r of adding alkali to glycerol is a constant. Under the above assumptions (H1) and (H2), mass balances of biomass, substrate and products in the batch process are written as follows:

C. Liu et al. / Applied Mathematical Modelling 37 (2013) 695–706

8 x_ 1 ðtÞ ¼ q1 ðxðtÞÞx1 ðtÞ; > > > > > x_ 2 ðtÞ ¼ q2 ðxðtÞÞx1 ðtÞ; > > > < x_ ðtÞ ¼ q ðxðtÞÞx ðtÞ; 3 1 3 > x_ 4 ðtÞ ¼ q4 ðxðtÞÞx1 ðtÞ; > > > > > x_ 5 ðtÞ ¼ q5 ðxðtÞÞx1 ðtÞ; > > : x_ 6 ðtÞ ¼ 0;

697

ð1Þ

where x1 ðtÞ; x2 ðtÞ; x3 ðtÞ; x4 ðtÞ; x5 ðtÞ and x6 ðtÞ are the concentrations of biomass, glycerol, 1,3-PD, acetate, ethanol and the volume of culture fluid at t in fermentor, respectively. xðtÞ :¼ ðx1 ðtÞ; x2 ðtÞ; x3 ðtÞ; x4 ðtÞ; x5 ðtÞ; x6 ðtÞÞT is the state vector. On the basis of the previous work [37], the specific growth rate of cells q1 ðxðtÞÞ is expressed by

q1 ðxðtÞÞ ¼ D1

x2 ðtÞ ; x2 ðtÞ þ k1

ð2Þ

where D1 is the maximum specific growth rate; k1 is the Monod saturation constant. The specific consumption rate of substrate q2 ðxðtÞÞ is

q2 ðxðtÞÞ ¼ m2 þ

q1 ðxðtÞÞ x2 ðtÞ þ D2 : Y2 x2 ðtÞ þ k2

ð3Þ

In (3), m2 is the maintenance term of substrate consumption under substrate-limited conditions. Y 2 is the maximum growth yield. D2 is the maximum increment of substrate consumption rate under substrate-sufficient conditions. k2 is the saturation constant for substrate. The specific formation rates q‘ ðxðtÞÞ; ‘ ¼ 3; 4, of 1,3-PD and acetate are defined as

q‘ ðxðtÞÞ ¼ m‘ þ q1 ðxðtÞÞY ‘ þ D‘

x2 ðtÞ ; x2 ðtÞ þ k‘

ð4Þ

where m‘ are the maintenance terms of product formations under substrate-limited conditions; Y ‘ are the maximum product yields; D‘ are the maximum increments of product formation rates under substrate-sufficient conditions; and k‘ are saturation constants for products. Moreover, the specific formation rate q5 ðxðtÞÞ of ethanol can be described by

 q5 ðxðtÞÞ ¼ q2 ðxðtÞÞ

 b1 b2 ; þ c1 þ q1 ðxðtÞÞx2 ðtÞ c2 þ q1 ðxðtÞÞx2 ðtÞ

ð5Þ

in which b1 ; b2 ; c1 and c2 are parameters for determination of yield of ethanol on glycerol. Due to the feed of glycerol and alkali in the fermentation process, there exist dilute effects on the concentrations of involving substances. Consequently, the batch model can be extrapolated to the feed process by incorporating the dilution factors. Namely, mass balances of biomass, substrate and products in the feed process are given below:

8 x_ 1 ðtÞ ¼ ðq1 ðxðtÞÞ  DðxðtÞ; uðtÞÞÞx1 ðtÞ; > > > cs0 > > x_ 2 ðtÞ ¼ DðxðtÞ; uðtÞÞð1þr  x2 ðtÞÞ  q2 ðxðtÞÞx1 ðtÞ; > > > < x_ ðtÞ ¼ q ðxðtÞÞx ðtÞ  DðxðtÞ; uðtÞÞx ðtÞ; 3 1 3 3 > _ 4 ðtÞ ¼ q4 ðxðtÞÞx1 ðtÞ  DðxðtÞ; uðtÞÞx4 ðtÞ; x > > > > > x_ 5 ðtÞ ¼ q5 ðxðtÞÞx1 ðtÞ  DðxðtÞ; uðtÞÞx5 ðtÞ; > > : x_ 6 ðtÞ ¼ ð1 þ rÞuðtÞ;

ð6Þ

where uðtÞ 2 R1 is the feed rate of glycerol in the feed process. cs0 > 0 denotes the concentration of initial feed of glycerol in the medium. r > 0 is the velocity ratio of adding alkali to glycerol. Since uðtÞ P 0, the volume x6 ðtÞ of solution is nondecreasing and x6 ðtÞ > 0 due to the positivity of the initial volume of culture fluid. Therefore, DðxðtÞ; uðtÞÞ is the dilution rate defined by

DðxðtÞ; uðtÞÞ ¼

ð1 þ rÞuðtÞ : x6 ðtÞ

ð7Þ

Now, let uðtÞ be the control function and suppose the switching instant si ; i 2 f1; 2; . . . ; 2Ng, between the batch and feed processes satisfies that 0 ¼ s0 < s1 < s2 <    < s2N < s2Nþ1 ¼ T. In particular, s2jþ1 ; j 2 K1 :¼ f0; 1; 2; . . . ; Ng is the moment of feed glycerol and alkali, s2j ; j 2 K2 :¼ f1; 2; . . . ; Ng is the moment of ending the feed, and N is a given constant in this paper. Furthermore, denote the right-hand item of the ith equation in the system (6) by fi ðxðtÞ; uðtÞÞ for each i 2 K3 :¼ f1; 2; . . . ; 6g and let

f ðxðtÞ; uðtÞÞ : ðf1 ðxðtÞ; uðtÞÞ; . . . ; f6 ðxðtÞ; uðtÞÞÞT : Then, the nonlinear dynamical system describing the whole process of fed-batch fermentation can be formulated as:

ð8Þ

698

C. Liu et al. / Applied Mathematical Modelling 37 (2013) 695–706

8 _ xðtÞ ¼ f ðxðtÞ; uðtÞÞ; > > > < uðtÞ 2 U ; t 2 ðs ; s ; i i1 i > xð s þÞ ¼ xð s Þ; i ¼ 1; 2; . . . ; 2N þ 1; i1 i1 > > : xð0Þ ¼ x0 ;

ð9Þ

where x0 is a given initial state, the notation + indicates the limit from the right, and

 Ui ¼

½ai ; bi ; if i even; f0g; if i odd:

ð10Þ

Here, a2j1 and b2j1 ; j 2 K2 , are positive constants which denote the minimal and the maximal rates of adding glycerol, respectively. Thus, we define the class of admissible control functions as

U :¼ fu 2 L2 ð½0; T; R1 ÞjuðtÞ 2 U i ;

i ¼ 1; 2; . . . ; 2N þ 1g;

ð11Þ

where L2 ð½0; T; R1 Þ is the space of square-integrable Lebesgue measurable functions from ½0; T into R1 . Since biological considerations limit the rate of switching, there are maximal and minimal time durations that are spent on each of the batch and feed processes. On this basis, define the set of admissible switching instants as

C :¼ fðs1 ; s2 ; . . . ; s2N ÞT 2 R2N jqi 6 si  si1 6 .i ; i 2 K :¼ f1; 2; . . . ; 2N þ 1gg;

ð12Þ

where qi and .i are the minimal and the maximal time durations, respectively. Accordingly, any s 2 C is regarded as an admissible vector of switching instants. There exist critical concentrations, outside which cells cease to grow, of biomass, glycerol, 1,3-PD, acetate and ethanol. Hence, it is biologically meaningful to restrict the concentrations of biomass, glycerol and products in a set W defined as

xT ðtÞ 2 W :¼

6 Q

½x‘ ; x‘ ;

8t 2 ½0; T:

ð13Þ

‘¼1

For the system (9), some important properties are discussed as follows. Theorem 1. The function f ð; Þ in the system (9) satisfies the following conditions: (a) f ð; Þ : R6þ 

2Nþ1 S i¼1

U i ! R6 , together with its partial derivatives with respect to x and u, are continuous on R6þ 

2Nþ1 S

Ui .

i¼1

(b) f ð; Þ is affine in control u, (c) There exists a constant K > 0 such that

kf ðxðtÞ; uðtÞÞk 6 KðkxðtÞk þ 1Þ;

8ðxðtÞ; uðtÞÞ 2 R6þ 

2Nþ1 S

Ui ;

ð14Þ

i¼1

where k  k is the Euclidean norm.

Proof (a) This conclusion can be obtained by the expressions of f in (1) and (6). (b) It is easy to verify that f is affine in control u by its definition. (c) We can complete the proof using a method similar to the proof of Property 1 in [38]. h

Theorem 2. For each u 2 U and s 2 C, the nonlinear dynamical system (9) has a unique continuous solution denoted by xðju; sÞ. Furthermore, xðju; sÞ satisfies the following integral equation

xðtju; sÞ ¼ x0 þ

Z

t

f ðxðsju; sÞ; uðsÞÞds;

8t 2 ½0; T

ð15Þ

0

and is continuous in u and s. Proof. This conclusion can be obtained from Theorem 1 and the theory of ordinary differential equations [39].

h

C. Liu et al. / Applied Mathematical Modelling 37 (2013) 695–706

699

Theorem 3. If xðju; sÞ is a solution of the system (9) with given initial condition x0 , then it is uniformly bounded. Proof. In view of Theorems 1 and 2, we obtain that, for each u 2 U and

kxðtju; sÞk 6 kx0 k þ

Z

t

kf ðxðsju; sÞ; uðsÞÞkds;

6 kx0 k þ K

Z

0

s 2 C,

t

kxðsju; sÞ þ 1kds;

8t 2 ½0; T:

0

By Gronwall inequality, it follows that

kxðtju; sÞk 6 M 01 ; where

M 01

8t 2 ½0; T;

:¼ ðkx0 k þ KTÞ expðKTÞ.

h

3. Constrained optimal control problems Basically, the control task of fed-batch fermentation lies in determination of the proper feed rate of glycerol and the switching instants between the batch and feed processes to obtain as much 1,3-PD as possible at the terminal time of the fermentation. Especially, physical limitations (11)–(13) are set during the fermentation process. Then, the constrained optimal control problem may now be stated formally as: Problem (P). Given the nonlinear dynamical system

8 _ xðtÞ ¼ f ðxðtÞ; uðtÞÞ; > > > < uðtÞ 2 U ; t 2 ðs ; s ; i i1 i > xð s þÞ ¼ xð s Þ; i ¼ 1; 2; . . . ; 2N þ 1; i1 i1 > > : xð0Þ ¼ x0 find a control u 2 U and a switching vector

s 2 C such that the constraint (13) is satisfied and the cost functional

Jðu; sÞ :¼ x3 ðTju; sÞ

ð16Þ

is minimized. By similar arguments as those given for Theorem 5 in [38], we confirm the existence of the optimal solution for Problem (P). Theorem 4. Problem (P) has at least one optimal solution.

4. Computational approaches In this section, we shall develop a numerical solution method to Problem (P). 4.1. Approximate problem (P(p)) For each pi P 1; i 2 f1; 2; . . . ; 2N þ 1g, let the time subinterval ½si1 ; si  be partitioned into npi subintervals with npi þ 1 partition points denoted by p i sp0i ; sp1i ; . . . ; spnipi ; sp0i ¼ si1 ; spnipi ¼ si ; and spk1 6 sk i :

Let npi be chosen such that npi þ1 P npi . The control is now approximated in the form of a piecewise constant function as follows:

up ðtÞ ¼

npi 2Nþ1 P P

rpi ;k vðspi

i¼1 k¼1

Here,

vðspi

k1

vðspi

k1

p

;sk i 

ðtÞ:

ð17Þ p

p

;ski 

p

p

i is the indicator function on the interval ðsk1 ; ski  defined by

;s i  k1 k

(

ðtÞ ¼

p

p

i 1; t 2 ðsk1 ; ski ;

0; otherwise:

P2Nþ1

Let i¼1 npi ¼ j. Then, rp ¼ ððrp1 ÞT ; . . . ; ðrp2Nþ1 ÞT ÞT 2 Rj , where mate control (17). From (11), it is clear that

rpi :¼ ðrpi ;1 ; . . . ; rpi ;npi ÞT defines the heights of the approxi-

rpi ;k 2 U i ; k ¼ 1; 2; . . . ; npi ; i ¼ 1; 2; . . . ; 2N þ 1: Let Np be the set of all those

rp which satisfy the constraints (18).

ð18Þ

700

C. Liu et al. / Applied Mathematical Modelling 37 (2013) 695–706 p

Note that ski ; k ¼ 1; 2; . . . ; npi ; i ¼ 1; 2; . . . ; 2N þ 1, taken as the decision variables will encounter numerical difficulties as mentioned in [31]. For this reason, a control parametrization enhance transform is introduced to map these variable time points into preassigned fixed knots in a new time scale. It is achieved by introducing a transform from t 2 ½0; T to s 2 ½0; 2N þ 1 as follows:

_ ¼ v p ðsÞ; tðsÞ where

v

p

tð0Þ ¼ 0;

ð19Þ

is given by np i 2Nþ1 P P

v p ðsÞ ¼

i¼1 k¼1

dki v

i ðsÞ:

p

k i1þk1 np ;i1þnp i

ð20Þ

i

p

In (20), dki P 0; k ¼ 1; 2; . . . ; npi ; i ¼ 1; 2; . . . ; 2N þ 1, are decision variables. Let dp be the vector whose components are p dki ; k ¼ 1; 2; . . . ; npi ; i ¼ 1; 2; . . . ; 2N þ 1, and Xp be the set of all such dp . Clearly,

tðsÞ ¼

Z

s

v p ðgÞdg ¼

0

np l i1 PP

p

dj l

l¼1j¼1 npi

þ

  k1 p p ; dj i þ dki s  i þ 1  npi j¼1

k1 P



8s 2 i  1 þ

k1 k ;i  1 þ npi npi

 ð21Þ

and

tð2N þ 1Þ ¼

np i 2Nþ1 P P i¼1 k¼1

p

dki ¼ T:

ð22Þ

Let

wp ðsÞ ¼ up ðtðsÞÞ: Then

wp ðsÞ ¼

npi 2Nþ1 P P

rpi ;k v

i¼1 k¼1

k i1þk1 np ;i1þnp i

i ðsÞ: i

Define

~xðsÞ :¼ ðxðsÞT ; tðsÞÞT and

~f ð~xðsÞ; rp ; dp Þ :¼ ððv p ðsÞf ðxðtðsÞÞ; wðsjrp ÞÞT ; v p ðsjdp ÞÞT : Let ~ xðjrp ; dp Þ be the solution of the following system corresponding to the control parameter vector ðrp ; dp Þ 2 Np  Xp :

8 p p ~ _ > < ~xðsÞ ¼ f ð~xðsÞ; r ; d Þ; ~xði  1þÞ ¼ ~xði  1Þ; > : ~xð0Þ ¼ ðxT0 ; 0ÞT :

s 2 ði  1; i; i ¼ 1; 2; . . . ; 2N þ 1;

Then, the constraint (13) can be rewritten as

ð~x1 ðsÞ; ~x2 ðsÞ; ~x3 ðsÞ; ~x4 ðsÞ; ~x5 ðsÞ; ~x6 ðsÞÞT 2 W;

8s 2 ½0; 2N þ 1:

ð23Þ

Now, we may specify the approximate problem (P(p)) as follows. Problem (P(p)). Subject to the system of differential equations

8 p p ~ _ > < ~xðtÞ ¼ f ð~xðsÞ; r ; d Þ; ~xði  1þÞ ¼ ~xði  1Þ; > : ~xð0Þ ¼ ðxT0 ; 0ÞT

s 2 ði  1; i; i ¼ 1; 2; . . . ; 2N þ 1;

find a combined vector ðrp ; dp Þ 2 Np  Xp such that the constraints (22) and (23) are satisfied and the cost functional

Jðrp ; dp Þ :¼ ~x3 ð2N þ 1jrp ; dp Þ

ð24Þ

is minimized. 4.2. Continuous state constraints Since constraint (23) in Problem (P(p)) is a continuous state constraint, Problem (P(p)) can be viewed as a semi-infinite programming problem. An efficient algorithm for solving optimization problems of this type is discussed in [36]. We will now brief discuss the application of this algorithm to Problem (P(p)).

C. Liu et al. / Applied Mathematical Modelling 37 (2013) 695–706

701

Let

g ‘ ð~xðsjrp ; dp ÞÞ :¼ ~x‘ ðsjrp ; dp Þ  ~x‘ ; g 6þ‘ ð~xðsjrp ; dp ÞÞ :¼ ~x‘  ~x‘ ðsjrp ; dp Þ;

‘ ¼ 1; 2; . . . ; 6:

The constraint (23) is equivalently transcribed into

Gðrp ; dp Þ ¼ 0;

ð25Þ P12 R T

p

p

where Gðr ; d Þ :¼ l¼1 0 maxf0; g l ð~ xðsjr ; d ÞÞgds. However, Gð; Þ is non-differentiable at the point g l ¼ 0; l 2 f1; 2; . . . ; 12g. We replace (25) with p

e e;c ðrp ; dp Þ :¼ c þ P G 12

p

Z 0

l¼1

where

2Nþ1

ue ðg l ð~xðsjrp ; dp ÞÞÞds P 0;

ð26Þ

e > 0; c > 0 and

ue ðgÞ ¼

8 > < g;

if g < e;

ðgþeÞ2 > 4e

:

ð27Þ

; if  e 6 g 6 e; if g > e:

0;

Thus, Problem (P(p)) is approximated by a sequence of Problems {ðP e;c ðpÞÞ} defined by replacing constraint (25) with (26). e e;c ð; Þ can be computed by the following theorem. Then, the gradient of the constraint G e e;c ðrp ; dp Þ given in (26), it holds that its gradients with respect to parameterized control rp and dp Theorem 5. For the constraint G are, respectively,

e e;c ðrp ; dp Þ @G ¼ @ rpi ;k

Z

i1þnkp

i

i1þk1 np

e ~xðsjrp ; dp Þ; rp ; dp ; ~kðsÞÞ @ Hð ds; @ rpi ;k

i

and

e e;c ðrp ; dp Þ @G pi ;k

@d

¼

Z

i1þnkp

i

e ~xðsjrp ; dp Þ; rp ; dp ; ~kðsÞÞ @ Hð @dpi ;k

i1þk1 np i

ds;

where

e ~xðsjrp ; dp Þ; rp ; dp ; ~kðsÞÞ ¼ Hð

12 P l¼1

ue ðg l ð~xðsjrp ; dp ÞÞÞ þ ~kT ðsÞ~f ð~xðsjrp ; dp Þ; rp ; dp Þ;

and

~kðsÞ ¼ ð~k1 ðsÞ; ~k2 ðsÞ; ~k3 ðsÞ; ~k4 ðsÞ; ~k5 ðsÞ; ~k6 ðsÞ; ~k7 ðsÞÞT is the solution of the costate system

!T e ~xðsjrp ; dp Þ; rp ; dp ; ~kðsÞÞ @ Hð _ ~kðsÞ ¼ ; @ ~x with the boundary conditions

~kð2N þ 1Þ ¼ ð0; 0; 0; 0; 0; 0; 0ÞT ; ~kðiÞ ¼ ~kðiþÞ; i ¼ 1; 2; . . . ; 2N: Proof. We can complete the proof using a method similar to the proof of Theorem 5.3. in [19]. h 4.3. Optimization algorithms Each of Problems {(P e;c ðpÞÞ} is a mathematical programming problem which can be solved by various optimization methods such as gradient-based techniques [29,30]. However, all those techniques are only designed to find local optimal solutions. In contrast, stochastic evolution methods generally lead to better results than deterministic ones. Among these stochastic evolution methods, PSO algorithm introducing by Kennedy and Eberhart in [40] exhibits distinct advantages. Presently, PSO has been applied to optimizing the fermentation processes [19,41–43]. In a typical PSO, each particle ‘flies’ over the search space to look for promising regions according to the experiences of both its own and those of the groups. Thus, the social sharing of information takes place and individuals profit from the discoveries and previous experiences of all other particles in a wide landscape during their search process around the better solutions. Traditionally, the original PSO method

702

C. Liu et al. / Applied Mathematical Modelling 37 (2013) 695–706

deals with unconstrained optimization problems. However, what we need to solve is an optimization problem with both control bounds and state constraints, to which the original PSO can not be applied directly. By the way, although there exist many constraint handling techniques in the evolutionary computation, see, for example [44], the treatment of continuous state constraints is rarely considered. We present a handling technique for this type of constraints in the sequel. As a result, the following strategies are added to the original PSO proposed in [40]. (I) (Handling the control constraint). If there is bound violation for a control parameter in the ith individual at the mth step, then that parameter is generated randomly between given lower and upper bounds using the following equation: p ;k i ;k i ;k  r i Þ; rpi i ;k ðmÞ ¼ rplow þ r pi i ;k ðrpupp k 2 f1; 2; . . . ; npi g; i 2 f1; 2; . . . ; 2N þ 1g; low p ;k

p ;k

i i where rlow and rupp are, respectively, the lower and the upper bounds of the control parameter which can be obtained by pi ;k (18). ri is a random parameter which is taken uniformly from ½0; 1, (II) (Dealing with the continuous state constraints). For the parameter of the ith individual at the mth step, test the value e e;c ðrpi ðmÞ; dp ðmÞÞ. If G e e;c ðrpi ðmÞ; dp ðmÞÞ ¼ 0, then the parameterized control is feasible. Otherwise, that is, of G

i

i

G e;c ðri ðmÞ;di ðmÞÞ e e;c ðrpi ðmÞ; dp ðmÞÞ > 0, move the parameter towards the feasible region in the direction of  e and G p i @ r ðmÞ p

p

i



p p e G e;c ðri ðmÞ;di ðmÞÞ p

@di ðmÞ

with Armijo line search.

(III) (Stopping criteria). The algorithm stops when the maximal iteration M is reached. Based on the above improved PSO algorithm, we can obtain an approximately optimal control and optimal switching instants for Problem (P) as shown in the following algorithm. Algorithm: 1 Step 1. Choose initial values of e and c. p; Step 2. Solve Problem ðP e;c ðpÞÞ using the improved PSO algorithm to give ðrp; e;c ; de;c Þ. p; xðsjrp; for s 2 ½0; 2N þ 1; l ¼ 1; 2; . . . ; 12. Step 3. Check feasibility of g l ð~ e;c ; de;c ÞÞ P 0 p; , where c  is a prespecified positive constant, go Step 4. If ðrp; e;c ; de;c Þ is feasible, go to Step 5. Otherwise, set c ¼ ac. If c < c

to Step 6. Otherwise go to Step 2. Step 5. Set e ¼ be. If e > e, where e is a prespecified positive constant, go to Step 3. Otherwise go to Step 6. npi P P, where P is a predefined positive constant, go to Step 7. Otherwise, go to Step 1 with npi Step 6. If min i2f1;2;...;2Nþ1g

increased to npi þ1 for each i. Step 7. Stop and construct up; and Then, ðu

p;

p;

;s

p; sp; from ðrp; e;c ; de;c Þ according to (17) and (21).

Þ obtained is an approximately optimal solution of Problem (P).

Remark 1. In the algorithm, e is a parameter controlling the accuracy of the smoothing approximation. c is a parameter controlling the feasibility of the constraint (23).

Table 1 The critical concentrations and the parameters in (2)–(5). ‘

m‘

Y‘

D‘

k‘

c‘

x‘

x‘

1 2 3 4 5 6

– 2.20 2.69 0.97 – –

– 113.6 67.69 33.07 – –

0.67 28.58 26.59 5.74 – –

0.28 11.43 15.50 85.71 – –

0.025 0.06 5.18 50.45 – –

0.01 15 0 0 0 4

6 2039 1036 1026 60.9 7

Table 2 The bounds of feed rates in Phs. I–IX [19]. Phases

I–II

III

IV–V

VI

VII

VIII–IX

Upper bounds [mL s1] Lower bounds [mL s1]

0.2524 0.1682

0.2390 0.1594

0.2524 0.1682

0.2657 0.1771

0.2924 0.1949

0.3058 0.2038

703

C. Liu et al. / Applied Mathematical Modelling 37 (2013) 695–706

. Especially, the Remark 2. It is important for the validity of the above algorithm to choose the parameters a; b; e and c  are two sufficient small values such that the algorithm is effective. parameters a and b must be chosen less than 1. e and c

5. Numerical results On the basis of the reactant composition, cultivation conditions, determination of biomass, substrate and metabolites were reported in [45]. Under anaerobic conditions at 37  C and pH 7.0, the critical concentrations x and x for the cells’ growth and the parameters in (2)–(5) are listed in Table 1. Moreover, to solve numerically the nonlinear dynamical system (9), the initial state, velocity ratio of adding alkali to glycerol, concentration of initial feed glycerol, fermentation time and the number of switchings are x0 ¼ ð0:1115 g L1 ; 495 mmol L1 ; 0; 0; 0; 5 LÞT ; r ¼ 0:75; cs0 ¼ 10762 mmol L1 ; T ¼ 24:16 hr, and 2N ¼ 1354, respectively. In order to save computational time, the fermentation process is partitioned into the first batch phase (Bat. Ph.) and phases I–IX (Phs. I–IX) according to the number of switchings. The same feed strategies are adopted in each one of Ph. I to Ph. IX. Furthermore, the time durations for two adjacent processes, i.e., a feed process and its succeeding batch process, s1 Þ in Ph. I to Ph. IX are equal and assumed to be 3600ðT s. It should be mentioned that this approach had been adopted to n obtain the experimental data in the actual fermentation process. Moreover, the bounds of feed rates in Phs. I–IX are listed in Table 2 [19]. The bounds of the time duration in each of Phs. I–IX are listed in Table 3 [27]. In the improved PSO algorithm, the number of initial particles swarm Np, the maximal iteration Mp, the inertia weight xp , and the acceleration constants c1 ; c2 are, respectively, 200, 100, 0.5, 2, 2. These parameters are derived empirically after numerous experiments. In Algorithm 1, the initial values of u and s are chosen as those in [19], the smoothing and feasible parameters were initially selected as e ¼ 0:1 and c ¼ 0:01, and then subsequently adjusted according to the guidelines in Algorithm 1. In particularly, the parameters a and b were chose as 0:1 and 0:01 until the solution obtained is feasible for  ¼ 1:0  107 . It is worth mentioned that in the original problem. The process was terminated when e ¼ 1:0  108 and c

Table 3 The bounds of time durations in the Bat. Ph. and Phs. I–IX [27]. Phases Bat. I ðj ¼ 1; . . . ; 28Þ II–V ðj ¼ 29; . . . ; 378Þ VI–VIII ðj ¼ 379; . . . ; 666Þ IX ðj ¼ 667; . . . ; 677Þ

Bounds

Values [s] 19080

q1 q2j q2jþ1 q2j q2jþ1 q2j q2jþ1 q2j q2jþ1

Bounds

.1 .2j .2jþ1 .2j .2jþ1 .2j .2jþ1 .2j .2jþ1

2 92 4 90 1 93 1 97

Table 4 The optimal switching instants in fed-batch fermentation process. Phases Bat.

s1

Switching instants

Optimal values [s] 19084.9

I (j ¼ 1; . . . ; 28) II (j ¼ 29; . . . ; 65) III (j ¼ 66; . . . ; 126) IV (j ¼ 127; . . . ; 245) V (j ¼ 246; . . . ; 378) VI (j ¼ 379; . . . ; 459) VII (j ¼ 460; . . . ; 522) VIII (j ¼ 523; . . . ; 666) IX (j ¼ 667; . . . ; 677)

s2j s2jþ1 s2j s2jþ1 s2j s2jþ1 s2j s2jþ1 s2j s2jþ1 s2j s2jþ1 s2j s2jþ1 s2j s2jþ1 s2j s2jþ1 ðj – 677Þ

19084.9 + 100.282 (j  1) 19092:543 þ 100:282j 21901.369 + 100.282 (j  29) 21892.8 + 100.282 (j  28) 25613.049 + 100.282 (j  66) 25603.2 + 100.282 (j  65) 31730.356 + 100.282 (j  127) 31720.4 + 100.282 (j  126) 43660.404 + 100.282 (j  246) 43654 + 100.282 (j  245) 56992.716 + 100.282 (j  379) 56991.6 + 100.282 (j  378) 65115.519 + 100.282j  460) 65114.5 + 100.282 (j  459) 71433.213 + 100.282 (j  523) 71432.2 + 100.282 (j  522) 85873.987 + 100.282 (j  667) 85872.9 + 100.282 (j  666)

Values [s] 19440 8 98 10 96 7 99 3 99

704

C. Liu et al. / Applied Mathematical Modelling 37 (2013) 695–706 0.3 Optimal Feeding Strategy (mL/s)

Optimal Feeding Strategy (L/h)

0.3 Bat. Ph.

0.25 0.2 0.15 0.1 0.05 0 0

1

2

3 Time (h)

4

Optimal Feeding Strategy (mL/s)

Optimal Feeding Strategy (mL/s)

0.25

0.3 0.2

0.2

0.1

0.15

0 0

0.1

2

4

6

8

10

0.05

40 60 Time (s)

80

4

6

8

10

20

40 60 Time (s)

80

100

Ph.III 0.25

0.3 0.2

0.2

0.1

0.15

0 0

0.1

2

4

6

8

10

0.05

20

40 60 Time (s)

80

100

0.3

0.25

Optimal Feeding Strategy (mL/s)

Optimal Feeding Strategy (mL/s)

0.3 0.2

0.2

0.1

0.15

0 0

0.1

2

4

6

8

10 12

0.05

Ph.V 0.25

0.3 0.2

0.2

0.1

0.15

0 0

0.1

2

4

6

8

10

0.05 0

20

40 60 Time (s)

80

100

0

20

40 60 Time (s)

80

100

0.3 Ph.VI

0.25

Optimal Feeding Strategy (mL/s)

Optimal Feeding Strategy (mL/s)

2

0.05

0

0.3

0.3 0.2

0.2

0.1

0.15

0 0

0.1

0.5

1

1.5

2

0.05

Ph.VII 0.25

0.3 0.2

0.2

0.1

0.15

0 0

0.1

0.5

1

1.5

2

0.05 0

20

40 60 Time (s)

80

100

0

0.3

20

40 60 Time (s)

80

100

0.3 Ph.VIII

0.25

Optimal Feeding Strategy (mL/s)

Optimal Feeding Strategy (mL/s)

0 0

0.1

100

Ph.IV

0.3 0.2

0.2

0.1

0.15

0 0

0.1

0.5

1

1.5

2

0.05 0 0

0.1

0.15

0 20

0.3

0 0

0.2

0.2

0.3 Ph.II

0 0

0.3

0 0

5

0.3

0 0

Ph.I 0.25

20

40 60 Time (s)

80

100

Ph.IX 0.25

0.3 0.2

0.2

0.1

0.15

0 0

0.1

0.5

1

1.5

2

0.05 0 0

20

40 60 Time (s)

80

100

Fig. 1. The optimal feeding strategy of glycerol in fed-batch fermentation process.

the first step, a small value of c was required to ensure feasibility. After that the c hardly changed as e was decreased. The specified constant P in Algorithm 1 is 2. Note also that only a small improvement (less than 0:01) was obtained by re-solving the problem with 5.

705

C. Liu et al. / Applied Mathematical Modelling 37 (2013) 695–706 1200 1,3−PD concentration in this paper. Experimental data [18]. 1,3−PD concentration in [19].

Concentration of 1,3−PD (mmolL−1)

1000

800

600

400

200

0

0

5

10 15 Fermentation Time (hr)

20

25

Fig. 2. The changes of 1,3-PD concentration with respect to time in fed-batch fermentation process.

Applying Algorithm 1 to the Problem (P), we obtain the optimal feed rates of glycerol in Phs. I–IX as shown in Fig. 1 and the optimal switching instants listed in Table 4. Here, all the computations are performed in Visual C++ 6.0 and numerical results are plotted by MATLAB 7.10.0. In particular, the ODEs in the computation process are numerically calculated by improved Euler method [46] with the relative error tolerance 104 . In detail, the blue1 line in the first subfigure of Fig. 1 indicates the feeding rate, which is identically equal to zero, of glycerol and the time duration in the Bat. Ph. Accordingly, the blue lines in the next 9 subfigures illustrate the feed rates of glycerol in conjunction with time durations of a feed process and its succeeding batch process in Ph. I to Ph. IX, respectively. To show the feed rates of glycerol in the feeding processes for Ph. I to Ph. IX better, 9 small subfigures are also incorporated in the corresponding 9 subfigures, respectively. Under the obtained optimal feed rates and the optimal switching instants, the computational concentration of 1,3-PD at the terminal time is 1025:3 mmol L1 which is increased by 28.64% in comparison with experimental result 797 mmol L1 appearing in [18] for numerical simulation. Furthermore, compared with the obtained 1,3-PD concentration 925:127 mmol L1 in [19], which is computed in case that the same number of phases is considered and the switching instants between the batch and feed processes are decided a priori, the concentration of 1,3-PD at the terminal time obtained in this paper is increased by 10.83%. Hence, it is decisive for enhancing the productivity of 1,3-PD to optimize the feed rate of glycerol and the switching instants between the batch and feed processes in fed-batch fermentation of glycerol to 1,3-PD. The concentration change of 1,3-PD obtained by the optimal strategy is shown in Fig. 2. For the purpose of comparison, the experimental data appearing in [18] and the 1,3-PD concentration obtained in [19] are also shown in Fig. 2. From Fig. 2, we conclude that 1,3-PD concentration at the terminal time in this paper is actually higher than the ones previously reported. 6. Conclusions In this paper, a nonlinear dynamical system was proposed to describe the bioconversion of glycerol to 1,3-PD in fed-batch fermentation. In order to obtain a high concentration of 1,3-PD at the terminal time, an optimal control mode subject to the continuous state constraints and the control constraint was presented. A computational approach was developed to seek the optimal solution of the constrained optimal control problem. Numerical results verified the validity of the mathematical model and the effectiveness of the computational method. Acknowledgements The supports of the Natural Science Foundation for the Youth of China (No. 11001153), the TianYuan Special Funds of the National Natural Science Foundation of China (No. 11126077), the Shandong Province Natural Science Foundation of China (No. ZR2010AQ016) and the National Natural Science Foundation of China (No. 11171050) are gratefully acknowledged. References [1] B.J. Minihane, D.E. Brown, Fed-batch culture technology, Biotechnol. Adv. 4 (1986) 207–218. [2] A. Mészáros, V. Bales, A contribution to optimal control of fed-batch biochemical processes, Bioproc. Eng. 7 (1992) 363–367. 1

For interpretation of colour in Fig. 1, the reader is referred to the web version of this article.

706

C. Liu et al. / Applied Mathematical Modelling 37 (2013) 695–706

[3] A. Korytowski, M. Szymkat, H. Maurer, G. Vossen, Optimal controlof a fedbatch fermentation process: numerical methods, sufficient conditions and sensitivity analysis, in: Proceedings of the 47th IEEE CDC, Cancun, Mexico, 2008, pp. 1551–1556. [4] H. Biebl, K. Menzel, A.P. Zeng, Microbial production of 1,3-propanediol, Appl. Microbiol. Biotechnol. 52 (1999) 289–297. [5] K. Menzel, A.P. Zeng, W.D. Deckwer, High concentration and productivity of 1,3-propanediol from continuous fermentation of glycerol by Klebsiella pneumoniae, Enzyme Microb. Technol. 20 (1997) 82–86. [6] H. Biebl, Glycerol fermentaion of 1,3-propanediol by Clostridium butyricum. Mesurement of product inhibition by use of a pH-auxostat, Appl. Microbiol. Biotechnol. 35 (1991) 701–705. [7] F. Barbirato, J.P. Grivet, P. Soucaille, A. Bories, 3-hydroxypropionaldehyde an inhibitory metabolite of glycerol fermentation to 1,3-propanediol by enterobacterial species, Appl. Environ. Microbiol. 62 (4) (1996) 1448–1451. [8] H.W. Chen, B.S. Fang, Z.D. Hu, Optimization of process parameters for key enzymes accumulation of 1,3-propanediol production from Klebisiella pneumoniae, Biochem. Eng. J. 25 (2005) 47–53. [9] A.P. Zeng, H. Biebl, Bulk-chemicals from biotechnology: the case of microbial production of 1,3-propanediol and the new trends, Adv. Biochem. Eng. Biotechnol. 74 (2002) 239–259. [10] P. Mhaskar, N.H. El-Farra, P.D. Christofides, Predictive control of switched nonlinear systems with scheduled mode transitions, IEEE Trans. Automat. Control 50 (2005) 1670–1680. [11] P. Mhaskar, N.H. El-Farra, P.D. Christofides, Stabilization of nonlinear systems with state and control constraints using Lyapunov-based predictive control, Syst. Control Lett. 55 (2006) 650–659. [12] P. Mhaskar, Robust model predictive control design for fault-tolearnt control of process systems, Ind. Eng. Chem. Res. 45 (2006) 8565–8574. [13] R. Luus, Application of dynamic programming to differential algebraic process systems, Comput. Chem. Eng. 17 (1993) 373–377. [14] H.S. Shin, H.C. Lim, Optimization of metabolite production in fed-batch cultures: use of sufficiency and characteristics of singular arc and properties of adjoint vector in numerical computation, Ind. Eng. Chem. Res. 46 (2007) 2526–2534. [15] C.X. Gao, E.M. Feng, Z.T. Wang, Z.L. Xiu, Nonlinear dynamical systems of bio-dissimilation of glycerol to 1,3-propanediol and their optimal controls, J. Ind. Manag. Optim. 3 (2005) 377–388. [16] C.X. Gao, K.Z. Li, E.M. Feng, Z.L. Xiu, Nonlinear impulsive system of fed-batch culture in fermentative production and its properties, Chaos Solitons Fract 28 (2006) 271–277. [17] H.Y. Wang, E.M. Feng, Z.L. Xiu, Optimality condition of the nonlinear impulsive system in fed-batch fermentation, Nonlinear Anal. R.W.A. 68 (2008) 12– 23. [18] G. Wang, E.M. Feng, Z.L. Xiu, Modelling and parameter identification of microbial biconversion in fed-batch cultures, J. Process Control 18 (2008) 458– 464. [19] C.Y. Liu, Z.H. Gong, E.M. Feng, H.C. Yin, Modeling and optimal control for nonlinear multistage dynamical system of microbial fed-batch culture, J. Ind. Manag. Optim. 5 (2009) 835–850. [20] Z.H. Gong, A multistage system of microbial fed-batch fermentation and its parameter identification, Math. Comput. Simulat. 80 (2010) 1903–1910. [21] C.Y. Liu, Optimal control for nonlinear dynamical system of microbial fed-batch culture, J. Comput. Appl. Math. 23 (2009) 252–261. [22] K.S. Narendra, J. Balakrishnan, Improving transient response of adaptive control systems using multiple models and switching, IEEE Trans. Automat. Control 39 (1994) 1861–1866. [23] A. Bemporad, A. Giua, C. Seatzu, A master-slave algorithm for the optimal control of continuous-time switched affine systems, in: Proceedings of the 41st IEEE CDC, Las Vegas, Nevada, 2002, pp. 1976–1981. [24] H.A.P. Blom, Y. Bar-Shalom, The interacting multiple model algorithm for systems with Markovian switching coefficients, IEEE Trans. Automat. Control 33 (8) (2002) 780–783. [25] X. Xu, P.J. Antsaklis, Optimal control of switched systems based on parameterization of the switching instants, IEEE Trans. Automat. Control 49 (2004) 2–16. [26] A. Ashoori, B. Moshiri, A. Jgaju-Sedigh, M.R. Bakhtiari, Optimal control of a nonlinear fed-batch fermentation process using model predictive approach, J. Process Control 19 (2009) 1162–1173. [27] Z.H. Gong, C.Y. Liu, E.M. Feng, L. Wang, Y.S. Yu, Modelling and optimization for a switched system in microbial fed-batch culture, Appl. Math. Model. 35 (2011) 3276–3284. [28] L. Cesari, Optimization-Theory and Applications, Springer-Verlag, New York, 1983. [29] A. Bryson, Y.C. Ho, Applied Optimal Control, Halsted Press, New York, 1975. [30] K.L. Teo, C.J. Goh, K.H. Wong, A Unified Computational Approach to Optimal Control Problems, Long Scientific Technical, Essex, 1991. [31] K.L. Teo, L.S. Jennings, H.W.J. Lee, V. Rehbock, The control parameterization enhancing transform for constrained optimal control problems, J. Aust. Math. Soc. Ser. B 40 (1999) 314–335. [32] V.K. Gorbunov, I.V. Lutoshkin, Development and experience of using the parameterization method in singular problems of dynamic optimization, J. Comput. Syst. Sci. Int. 43 (2004) 725–742. [33] R. Li, K.L. Teo, K.H. Wong, G.R. Duan, Control parameterization enhancing transform for optimal control of switched systems, Math. Comput. Model. 43 (2006) 1393–1403. [34] C.Z. Wu, K.L. Teo, Global impulsive optimal control computation, J. Ind. Manag. Optim. 2 (2007) 435–450. [35] R.C. Loxton, K.L. Teo, V. Rehbock, K.F.C. Yiu, Optimal control problems with a continous inequality constraint on the state and the control, Automatica 45 (10) (2009) 2250–2257. [36] K.L. Teo, V. Rehbock, L.S. Jennings, A new computational algorithm for functional inequality constrained optimization problems, Automatica 29 (1993) 789–792. [37] Z.L. Xiu, A.P. Zeng, L.J. An, Mathematical modelling of kinetics research on multiplicity of glyceroll bioconversion to 1,3-propanediol, J. Dalian Univ. Technol. 4 (2000) 428–433. [38] C.Y. Liu, E.M. Feng, H.C. Yin, Optimal switching control for microbial fed-batch cultures, Nonlinear Anal. H.S. 2 (2008) 1168–1174. [39] D.K. Arrowsmith, C.M. Place, Ordinary Differential Equations, Chapman and Hall, London, 1982. [40] J. Kennedy, R.C. Eberhart, Particle swarm optimization, in: Proceedings of IEEE International Conference on Neural Networks, Perth, Australia, 1995, pp. 1942–1948. [41] G. Wang, E.M. Feng, Z.L. Xiu, Vector measure as controls for explicit nonlinear impulsive system of fed-batch culture, J. Math. Anal. Appl. 351 (2009) 120–127. [42] F. Ferrari, S. Gutierrez, E.C. Biscaia Jr., Development of an optimal operation strategy in a sequential batch reactor (SBR) through mixed-integer particle swarm dynamic optimization (PSO), Comput. Chem. Eng. 34 (2010) 1994–1998. [43] J.X. Ye, Y.D. Zhang, E.M. Feng, Z.L. Xiu, H.C. Yin, Nonlinear hybrid system and parameter identification of microbial fed-batch culture with open loop glycerol input and pH logic control, Appl. Math. Model. 36 (2012) 357–369. [44] Z. Michalewicz, A survey of constraint handling techniques in evolutionary computation methods, in: Proceedings of the 4th Annual Conference on Evolutionary Programming, MIT Press, Cambridge, MA, 1995, pp. 135–155. [45] X. Chen, D.J. Zhang, W.T. Qi, S.J. Gao, Z.L. Xiu, P. Xu, Microbial fed-batch production of 1,3-propanediol by Klebsiella pneumoniae under microaerobic conditions, Appl. Microbiol Biotechnol. 63 (2003) 143–146. [46] E. Hairer, S.P. Nørsett, G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, Springer-Verlag, Berlin, Heidelberg, 2009.