Time-delay optimal control of a fed-batch production involving multiple feeds

In this paper, we consider time-delay optimal control of 1, 3-propan-ediol (1, 3-PD) fed-batch production involving multiple feeds. First, we propose a nonlinear time-delay system involving feeds of glycerol and alkali to formulate the production process. Then, taking the feeding rates of glycerol and alkali as well as the terminal time of process as the controls, we present a time-delay optimal control model subject to control and state constraints to maximize 1, 3-PD productivity. By a time-scaling transformation, we convert the optimal control problem into an equivalent problem with fixed terminal time. Furthermore, by applying control parameterization and constraint transcription techniques, we approximate the equivalent problem by a sequence of finite-dimensional optimization problems. An improved particle swarm optimization algorithm is developed to solve the resulting optimization problems. Finally, numerical results show that 1, 3-PD productivity increases considerably using the obtained optimal control strategy.

1. Introduction. 1,3-Propanediol (1,3-PD) is an odorless viscous liquid for synthesis of polyesters, polyurethanes and heterocyclic compounds [9]. It has numerous applications in solvents, antifreeze, lubricants and other fields [24]. 1,3-PD can be produced from renewable resources such as glycerol using microorganisms with advantages of mild conditions, simple operation, and fewer by-products [2].
Glycerol can be converted to 1,3-PD via one of three microbial production modes: batch production, continuous production and fed-batch production. In particular, fed-batch production is a mixture of the batch and continuous productions. During the bioconversion of glycerol to 1,3-PD, the most efficient cultivation method appears to be a fed-batch culture which corrects pH by alkali addition with glycerol supply [34]. Fed-batch production is commonly used in industrial production due to its ability to overcome catabolite repression or glucose effect which usually occurs during production of fine chemicals. A proper feeding rate, with the right component constitution, is required in order to improve production during the process. Moreover, fed-batch production also gives the operator of freedom of manipulating the process via the feeding rates. Therefore, optimal control of feeding rates in fed-batch production has received extensive attention. Assuming the feeding of substrate only occurs at the impulsive instants, optimal impulsive control of fed-batch production is discussed in [1,4,27]. Nonetheless, since the feeding rate of substrate is finite, it is not reasonable to describe the fed-batch production by the impulsive process. In contrast, taking the feeding rate of substrate as a time-continuous process, optimal multistage control of fed-batch production is investigated in [17]. Furthermore, regarding the fed-batch production as switching between batch process and feeding process, optimal switched control is studied in [30]. However, time-delays are ignored in the above optimal control problems. In fact, like most real processes, fed-batch production is also influenced by timedelays since nutrient metabolization does not immediately lead to the production of new biomass [28]. Recently, time-delay optimal control of fed-batch production is discussed in [15,18,19]. More recently, a published book [16] summarizes some optimal control results arising in fed-batch production processes. Although the results obtained are interesting, the control action in the above optimal control problems only includes the feeding rate of substrate glycerol and the feeding rate of alkali is obtained by an empirical velocity ratio of alkali to glycerol. Obviously, this cannot guarantee the obtained feeding rate of alkali to be optimal.
In this paper, we propose a nonlinear time-delay system involving feeds of both glycerol and alkali to formulate 1,3-PD fed-batch production. Taking 1,3-PD productivity as the performance index and the feeding rates of glycerol and alkali as well as the terminal time of process as the controls, we present a time-delay optimal control model subject to control and state constraints. Then, by a time-scaling transformation, we convert the optimal control problem into an equivalent problem with fixed terminal time. Furthermore, by exploiting control parameterization and constraint transcription techniques [26], we approximate the equivalent problem by a sequence of finite-dimensional optimization problems. An improved particle swarm optimization (PSO) algorithm is developed to solve the resulting optimization problems. Finally, numerical results show that 1,3-PD productivity increases considerably by utilizing the obtained optimal control strategy.
The organization of this paper is as follows. In Section 2, a nonlinear time-delay system involving multiple feeds for 1,3-PD fed-batch production is formulated. Section 3 gives time-delay optimal control model and its equivalent form. Section 4 develops a numerical solution method to solve the equivalent problem, while Section 5 illustrates the numerical results. Finally, conclusions are provided in Section 6.

2.
Nonlinear time-delay systems. 1,3-PD fed-batch production switches between batch process and feeding process. In batch process, nothing is added to and removed from the reactor. In feeding process, glycerol and alkali are continuously added to provide nutrition and maintain a suitable environment for cells growth. In particular, time-delays exist in the production process [28]. According to the production process, we assume that  Under Assumptions 1 and 2, mass balances of biomass, substrate and products in fed-batch production can be formulated as the following nonlinear time-delay OPTIMAL CONTROL OF FED-BATCH PROCESS 1699 system: where x(t) := (x 1 (t), x 2 (t), x 3 (t), x 4 (t), x 5 (t)) ∈ R 5 is the state vector whose components are, respectively, the concentrations of biomass, glycerol, 1,3-PD, acetic acid and ethanol; u(t) := (u 1 (t), u 2 (t)) ∈ R 2 is, the feeding rate vector of glycerol and alkali, the control vector; h is a given delay argument; T is the free terminal time; φ : R → R 5 is a given history function; and In (2), c s0 > 0 denotes the concentration of initial feed of glycerol in the medium. D(t, u(t)) is the dilution rate defined as where V (t) := V 0 + t 0 (u 1 (s) + u 2 (s))ds, and V 0 is the initial volume of culture fluid in the reactor. Based on the previous work [14], the specific growth rate of cells µ(x(t)), the specific consumption rate of substrate q 2 (x(t)), and the specific formation rates of products q (x(t)), = 3, 4, 5, are expressed as the following equations: where , and c 4 are kinetic parameters; and x * are critical concentrations for cells growth. Let 2N + 1 be the total number of batch and feeding processes (N feeding processes and N + 1 batch processes, since the fed-batch process starts and finishes in batch process) in fed-batch process. Then, we denote the start moment of the batch process by t 2j , and the start moment of feeding process by t 2j+1 , j ∈ {0, 1, . . . , N }. Note that the number of 2N + 1 is determined by the terminal time T . Therefore, T ∈ R 1 is called an admissible terminal time if it satisfies the following bound constraint: where T min and T max are, respectively, the lower and upper bounds of terminal time.
Let T be the set of all such admissible terminal times. Furthermore, define where a 1 i > 0 and b 1 i > 0 are, respectively, the minimal and maximal feeding rates of glycerol; and a 2 i > 0 and b 2 i > 0 are, respectively, the minimal and maximal feeding rates of alkali. Thus, Any essentially bounded function u from [0, called an admissible control. Let U be the class of all such admissible controls. Any pair (u, T ) ∈ U × T is called an admissible pair for time-delay system (1). It should be noted that there exist critical concentrations of biomass, glycerol, 1,3-PD, acetate and ethanol, outside which cells cease to grow. Thus, it is biologically meaningful to restrict the concentrations of biomass, glycerol and products within a set W defined as where x * are the lower thresholds for cell growth of biomass, glycerol, 1,3-PD, acetic acid, and ethanol, respectively; and x * (x * 2 , x * 3 , x * 4 and x * 5 as used in the formula for µ(x(t))) are the corresponding upper thresholds.
3. Time-delay optimal control problems. Let x(·|u, T ) be the continuous solution of time-delay system (1) on [−h, T ] corresponding to each (u, T ) ∈ U × T . It is desired that 1,3-PD productivity is maximized at the terminal time in 1,3-PD fed-batch production process, where 1,3-PD productivity is defined as As a result, taking 1,3-PD productivity as the performance index and incorporating the constraints (8), (10) and (11), we can state the time-delay optimal control problem as follows: Problem (OCP) is a time-delay optimal control problem with free terminal time in system (1). It is difficult to solve numerically because time-delay system (1) must be integrated over a variable time horizon. To surmount this difficulty, we apply the following time-scaling transformation [20]: where s ∈ [−h, 1] is a new time variable withh = h/T . Then, time-delay system (1) can be transformed into an equivalent form as follows: ,ũ(s)); andφ(s) := φ(T s). In addition, s 2N +1 = 1, and the switching moments and the set of admissible controls becomes where L ∞ ([0, 1], R 2 ) denotes the Banach space of all essentially bounded functions Letx(·|ũ, T ) be the continuous solution of system (14) on [−h, 1] corresponding to each (ũ, T ) ∈Ũ × T . Then, constraint (11) turns intõ Thus, Problem (OCP) is converted into the following equivalent problem: Note that Problem (EOCP) is a time-delay optimal control problem with fixed terminal time in system (14). 4. Numerical solution methods. Problem (EOCP) is, in essence, an optimal control problem. It is well known that the control parameterization technique is very efficient in solving optimal control problems [6,5,7,13,26,31].
Let Ξ p denote the set of all σ p satisfying constraint (21). Substituting (19) into system (14) Letx(·|σ p , T ) be the continuous solution of system (22) on [−h, 1] corresponding to each (σ p , T ) ∈ Ξ p × T . Then, constraint (17) becomes Thus, Problem (EOCP) can be approximated by the following finite-dimensional optimization problem: Problem (EOCP(p)) is a semi-infinite programming involving continuous state inequality constraint (23). As is well known, it is difficult to numerically solve (EOCP(p)) directly. By the way, to numerically computation constraint (23), an integral penalty method [3] and an equivalent end-point constraints method [23] were introduced. Nevertheless, a common characteristic of all these techniques is that the penalty terms or end-point constraints introduced have zero gradients with respect to optimization variables at the optimal solution. This, in turn, can result in a reduced convergence rate near the optimal solution. This problem can be dealt with by two approaches: a constraint transcription technique [10,26,29]; and an exact penalty method [11,32]. In this paper, we choose to apply the constraint transcription technique [26] to approximate constraint (23) due to its effectiveness and simplicity. Thus, let However, constraint (24) is non-smooth since max{0, ·} is non-differentiable at the origin. Thus, we replace constraint (24) with the following smooth inequality constraint in canonical form: where > 0 and γ > 0 are two adjusting parameters; and For each > 0 and γ > 0, let Problem (EOCP ,γ (p)) be Problem (EOCP(p)) replacing constraint (23) with (25). Note that Problem (EOCP ,γ (p)) can be viewed as a standard mathematical programming. In addition, it can be shown, as in [26], that for each > 0, there exists a corresponding γ( ) > 0 such that whenever 0 < γ < γ( ), constraint (25) implies constraint (23).
In the numerical computation, the gradients ofG ,γ (σ p , T ) with respect to σ p and T are required. Definē Obviously, for almost all s ∈ [−h, 1], we haveẋ(s|σ p , T ) = ψ(s|σ p , T ). Note that, in the sequel, we will use ∂x to denote partial differentiation with respect tox(s −h).
The following theorem provides these required gradients.
Proof. The proof is similar to that given for Theorem 4 in [18].
Based on the above theorem, Problem (OCP) can be solved by a sequence of problems {(EOCP ,γ (p))}. Each of {(EOCP ,γ (p))} is a smooth mathematical programming problem which can be solved by gradient-based techniques [3,26]. However, these algorithms are designed to find local optimal solutions. To overcome this difficulty, we introduce an improved PSO algorithm to solve each of {(EOCP ,γ (p))}. PSO was developed by Kennedy and Eberhart in the study of artificial life [8,22]. At present, PSO has attracted wide attention in neural network, optimization and other fields [12,33]. In PSO, each solution of the considered optimization problem is known as a particle. Each particle flights at a certain velocity according to its own and the swarm's experience to dynamically adjust the velocity, flying to the global best position, making the optimization problem to achieve the optimal solution. Nevertheless, each of {(EOCP ,γ (p))} is an optimization problem with constraints (8), (21) and (25). As a result, traditional PSO cannot be exploited to solve each of {(EOCP ,γ (p))}. Based on Theorem 4.1, we propose some improved strategies to solve each of {(EOCP ,γ (p))}. Assume that N p particles are evolving, the position and velocity of the ıth particle can be expressed in terms of σ p ı = (σ p ı,1 , . . . ,σ p ı,κ ) and ν p ı = (ν p ı,1 , . . . , ν p ı,κ ) , where κ = 2n pi + 1. In particular, the lower bound and upper bound of the position are represented byσ low andσ upp , respectively. In detail, the implementation strategies in the improved algorithm are given as follows.
• Updating position and velocity For the th component of the ıth particle at the (k + 1)th iteration, update the position and velocity with the following strategy: , where pb p ı = (pb p ı,1 , . . . , pb p ı,κ ) is the best position for ı particle in history iteration; gb p ı = (gb p ı,1 , . . . , gb p ı,κ ) is the best position in the swarm; r 1 ı, , r 2 ı, and r 3 ı, are the random numbers within the range of [0, 1]; c 1 (k) and c 2 (k) are two coefficients; and ω(k) is an inertia weight. The coefficients and inertia weight are defined as where it max is the maximal number of iterations; ω max and ω min are, respectively, the maximal and minimal inertia weights.

• Handling the state constraints
For the position of the ıth particle at the kth iteration, test the value of G(σ p ı (k)). If G(σ p ı (k)) = 0, then the position is feasible. Otherwise, G(σ p ı (k)) > 0 and move the position towards the feasible region based on the gradient formulae in Theorem 4.1, where and γ are adjusted according to the − γ process in [20].

• Stopping criterion
The algorithm stops when the maximal iteration it max is reached.
5. Numerical results. Consider a 1,3-PD fed-batch production process by Klebsiella pneumoniae reported in [18]. In the numerical simulation, we use the same settings as those to obtain the experimental results to optimize the feeding rates and the terminal time. In particular, the maximal duration of fed-batch process is partitioned into the first batch phase (Ph. I) and phases II-X (Phs. II-X). Within each one of Phs. II-X, all batch processes have 100s minus the duration of the feeding processes, and the same feeding rates of glycerol and alkali are adopted. The characteristics of each phase are given in  (1), the history function, the initial volume of culture fluid, the concentration of initial feed of glycerol in the medium, and time-delay are, respectively, φ(t) = (0.1115gL −1 , 495mmolL −1 , 0, 0, 0) , V 0 = 5L, c s0 = 10762mmolL −1 , and h = 0.4652h. Under anaerobic conditions at 37 • C and pH 7.0, the kinetic parameters and the critical concentrations for cells growth in system (1) are listed in TABLE 2.
In the improved PSO algorithm, the number of initial particles swarm N p , the maximal iteration it max are 50 and 100, respectively. These parameters are derived empirically after numerous experiments. The lower and upper bounds of terminal time are T min = 11h and T max = 24.16h, respectively. In addition, the lower and upper bounds for the feeding rates of glycerol and alkali are listed in TABLE 3. In handling the state constraint, = 1.0 × 10 −2 and γ = 2.5 × 10 −3 are taken as the initial values. We reduce γ by a factor of 2 if the control satisfies (24), or reduce both and γ by a factor of 10 if the control does not satisfy (24). The − γ process is terminated when ≤ 10 −8 .
By applying the improved PSO algorithm with control parameterization method, we obtain the optimal terminal time T * = 13.497h and the obtained optimal feeding strategy of glycerol and alkali are plotted in FIGURE 1. Under these optimal feeding rates and terminal time, we obtain that 1,3-PD productivity at the optimal terminal time is 73.397mmolh −1 , which increases 50.793% compared with the experiment result 48.674mmolL −1 in [21]. We also plot the concentration changes of biomass, glycerol and 1,3-PD with respect to the fermentation time in FIGURE 2. Furthermore, the changes of 1,3-PD productivity with respect to the fermentation  time are shown in FIGURE 3. From FIGURE 3, we confirm that the obtained 1,3-PD productivity is better than the previous results in [21].   Figure 3. 1,3-PD productivity changes with respect to fermentation time. Stars represent the 1,3-PD productivity in experiment [21], and solid line denotes the 1,3-PD productivity in this work.
6. Conclusions. This paper has considered time-delay optimal control of 1,3-PD fed-batch production involving multiple feeds. Taking 1,3-PD productivity as the performance index and the feeding rates of glycerol and alkali together with the terminal time of process as the controls, we presented a time-delay optimal control model subject to control and state constraints. To the best of our knowledge, this is the first time-delay optimal control problem involving feeds of glycerol and alkali in the literature for optimizing the fed-batch production. By applying control parameterization and constraint transcription techniques, the optimal control problem is approximated by a sequence of finite-dimensional optimization problems. Furthermore, an improved PSO algorithm is developed to solve the resulting finitedimensional optimization problems. Numerical results show that, by employing the obtained optimal control strategy, 1,3-PD productivity increases considerably. In closing, we note that the time-delay in system (1) is assumed to be time-invariable. However, the time-delay in fed-batch process is generally time-variable. Thus, designing optimal feeding strategy in the presence of variable time-delay is an interesting area to pursue for future research.