Theoretical Study and Optimization of the Biochemical Reaction Process by Means of Feedback Control Strategy

Abstract

The aim of this work is to present a theoretical analysis and optimization of a biochemical reaction process by means of feedback control strategy. To begin with, a mathematical model of the biochemical reaction process with feedback control is formulated. Then, based on the formulated model, the analysis of system's dynamics is presented. The optimization of the bioprocess is carried out, in order to achieve maximal biomass productivity. It is shown that during the optimization, the bioprocess with impulse effects loses the possibility of synchronization and strives for a simple continuous bioprocess. The analytical results presented in the work are validated by numerical simulations for the Tessier kinetics model.

1. Introduction

Chemical and bioprocess engineering play an important role in the production of many chemical products. In particular, bioreactor engineering as a branch of chemical engineering and biotechnology is an active area of research on bioprocesses, including among others development, control, and commercialization of new technology [1]. The reaching of optimal results and the obtainment of maximal profits require modern control strategies based on mathematical models or artificial intelligence methods [2]. There are many advantages with the use of mathematical models, and one of which is the possibility of testing process stability [3]. According to different reactions and differential control technologies, many dynamic biochemical models in a chemostat have been established [4–12]. However, there are many factors which affect the growth and reproduction of microorganisms in the bioprocess. The key variable is the biomass concentration, because, among other things, it provides information concerning biomass productivity. For that reason, almost all mathematical growth models contain the biomass as an important variable [13]. Moreover, for aerobic microbes, the dissolved oxygen concentration (DOC) is an important growth factor. During the growth of microorganisms, the dissolved oxygen concentration depends on several factors, among them the biomass concentration, the concentration and type of substrate used, and the bioprocess conditions. Because a low level of dissolved oxygen concentration decreases the biomass yield and the specific growth rate, it is necessary to monitor and control the DOC level so that it stays within the appropriate range [14].

Since many biological phenomena such as bursting rhythm models in, for example, medicine, biology, pharmacokinetics, and frequency modulated systems exhibit impulsive effects [15], impulsive differential equations, which appear as a natural description of observed evolution phenomena of several real world problems, have been introduced in different kind of biological systems, for example, in population dynamics [16–19] and in biochemical process [20, 21]. It should be pointed out that in these works, the authors were concerned about the fixed time impulse effects, which has the rationality in describing the biological phenomena. While in some cases, using the state-dependent impulse effects to describe the biological phenomena is more appropriate. As far as the state-dependent impulse effect is concerned, Tang and Chen [22] introduced a Lotka-Voterra model, which is constructed according to the practices of IPM. By using analytical method, it is shown that there exists an orbitally asymptotically stable periodic solution with a maximum value no larger than the given economic threshold. Further, the complete expression of period of the periodic solution is given. More researches on the applications of the state-dependent impulse effects can be found in [23–36]. In these researches, the authors analyzed the proposed system's dynamic behavior (e.g., the existence and stability of period-1 solution and the existence of period-2 solution) by applying the Poincaré principle and Poincaré-Bendixson theorem of the impulsive differential equation. Recently, in microbioprocess engineering, the dynamic properties of the kinetic models with impulse effects characterized by the universal microorganism growth rate and two different kinds of biomass yield are also analyzed theoretically by Sun et al. [37–39].

The objective of this work is to illuminate theoretical and practical aspects of the nonlinear analysis of a universal mathematical model of the biochemical reaction process. The paper is organized as follows. In Section 2, a universal mathematical model of the biochemical reaction process with any characteristics of growth kinetics is formulated under impulse effects. In Section 3, a qualitative analysis of the proposed model is presented. In Section 4, simulations with Tessier's kinetics are presented in order to verify the theoretical results and discuss the biological essence. Moreover, in Section 4, in order to optimize the biochemical reaction process, the conception of objective function is presented, and next the bioprocess optimization is put forward. Finally, in Section 5, we offer our conclusions.

2. Mathematical Model and Preliminaries

The general model of continuously culturing microorganism in a chemostat is given by the following form of differential equations [40]: 𝑑𝑆𝑆𝑑𝑡=𝐷_in−1−𝑆𝑌_𝑥/𝑆𝜇(𝑆)𝑥,𝑑𝑥=𝜇(𝑆)𝑆(0)𝑑𝑡−𝐷𝑥,=𝑆₀,𝑥(0)=𝑥₀, (1) where 𝑥=𝑥(𝑡) denotes the biomass concentration (g/L), and 𝑆=𝑆(𝑡) the substrate concentration (g/L) in the bioreactor medium at time 𝑡, 𝑥₀, and 𝑆₀ denotes the initial biomass concentration and substrate concentration, respectively; 𝐷 is the dilution rate (h⁻¹); 𝑆_in is the influent substrate concentration (g/L); 𝑌_𝑥/𝑆 is the biomass yield (g/g) defined as the ratio of the biomass 𝑥 produced to the amount of substrate 𝑆 consumed; 𝑌_𝑥/𝑆<1 for biological constraints; the function 𝜇(𝑆) describes the biochemical kinetics, which are characterized by the cell concentration (𝑥), depending on one limiting substrate with concentration 𝑆.

Crooke et al. [6] pointed out that the biochemical kinetics expression plays an important role in the intrinsic oscillation mechanisms. So many different assumptions for the kinetics models are given in the literature, for example, Monod-type kinetics [3], that is, 𝜇(𝑆)=𝜇_max𝑆/(𝐾_𝑆+𝑆), where 𝜇_max is the maximum specific growth rate and 𝐾_𝑆 is the substrate saturation constant; Tessier-type kinetics [41], that is, 𝜇(𝑆)=𝜇_max(1−𝑒^{−𝑆/𝐾𝑆}); Moser-type kinetics [41], that is, 𝜇(𝑆)=𝜇_max(1+𝐾_𝑆𝑆^−𝜗)⁻¹, where 𝜗 is a positive constant. It can be easily seen that these kinetic models satisfy the following properties: (i)(regularity) 𝜇: ℝ⁺→ℝ⁺ is continuously differential and 𝜇(0)=0, (ii)(monotonicity) 𝜇 is monotonically increasing, that is, 𝑑𝜇/𝑑𝑆≥0 for all 𝑆∈ℝ⁺; 𝜇^′(𝑆)=0 for some 𝑆∈ℝ⁺ means that 𝜇^′(𝑆)=0 for 𝑆≥𝑆, (iii)(convexity) 𝜇(𝜃𝑆₁+(1−𝜃)𝑆₂)≥𝜃𝜇(𝑆₁)+(1−𝜃)𝜇(𝑆₂) for all 𝜃∈[0,1] and 𝑆₁,𝑆₂∈ℝ⁺.

Thus, in the study that the kinetic models only satisfy the regularity, monotonicity, and convexity are considered.

According to the Herbert's and Pirt's models [41], if the substrate concentration is big enough (i.e., 𝑆≫𝐾_𝑆), the biomass yield is constant. Because in continuous bioprocesses the optimal biomass productivity is obtained for bigger substrate concentration, in this work the constant biomass yield is assumed, that is, 𝑌_𝑥/𝑆=𝑌_𝑥/𝑆max=𝜆<1. In addition, according to the design ideas of the bioreactor, the biomass concentration in the chemostat should be controlled to a set level 𝑥_set, where 0<𝑥_set≤𝑥_CRI, and 𝑥_CRI is the critical level of biomass concentration in the bioreactor medium. When the biomass concentration reaches the set level, part of the medium containing biomass and substrate is discharged from the bioreactor, followed by the input of the “clear” medium (i.e., the medium without substrate). Therefore, system (1) can be modified by introducing the impulsive state feedback control as follows: 𝑑𝑆𝑆𝑑𝑡=𝐷_in−1−𝑆𝜆𝜇(𝑆)𝑥𝑑𝑥=𝜇(𝑆)𝑥𝑑𝑡−𝐷𝑥<𝑥_set,Δ𝑆=−𝑊_𝑓𝑆Δ𝑥=−𝑊_𝑓𝑥𝑥=𝑥_set,𝑆0⁺=𝑆₀,𝑥0⁺=𝑥₀, (2) where 𝑊_𝑓 is the part of “clear” medium inputted into the bioreactor in each biomass oscillation cycle, which satisfies 0≤𝑊_𝑓≤𝑊_𝑓max<1, where 𝑊_𝑓max is the maximum part of “clear” medium inputted into the bioreactor in each biomass oscillation cycle.

Before presenting the main results, we recall the following definitions and lemmas first [31, 42, 43].

Definition 1. (𝜉(𝑡),𝜂(𝑡)) is said to be period-1 solution if in a minimum cycle time there is one impulse effect. Similarly, (𝜉(𝑡),𝜂(𝑡)) is said to be period-2 solution if in a minimum cycle time there are two impulse effects.

Definition 2. Γ is said to be orbitally stable, if for any 𝜀>0, there exists 𝛿>0, with the proviso that every solution (𝜉(𝑡),𝜂(𝑡)) of system (2) whose distance from Γ is less than 𝛿 at 𝑡=𝑡₀, will remain within a distance less than 𝜀 from Γ for all 𝑡≥𝑡₀. Such a Γ is said to be orbitally asymptotically stable if, in addition, the distance of (𝜉(𝑡),𝜂(𝑡)) from Γ tends to zero as 𝑡→∞. Moreover, if there exist positive constants 𝛼,𝛽, and a real constant 𝐶 such that 𝜌((𝜉(𝑡),𝜂(𝑡),Γ)<𝛼𝑒^−𝛽𝑡 for 𝑡>𝐶, then Γ is said to be orbitally asymptotically stable and enjoys the property of asymptotic phase.

Definition 3. Lambert 𝑊 function is defined to be a multivalued inverse of the function 𝑧→𝑧𝑒^𝑧 satisfying Lambert𝑊(𝑧)exp(Lambert𝑊(𝑧))=𝑧. (3)

In the followings, we simply denote 𝑊≜Lambert𝑊. Then, it follows from (3) that 𝑊^′(𝑧)=𝑊(𝑧)𝑧((𝑧,1+𝑊)) (4) when 𝑧≠0 and 𝑧≠−1/𝑒. First of all, the function 𝑧exp(𝑧) has the positive derivative (𝑧+1)exp(𝑧) if 𝑧>−1. Define the inverse function of 𝑧exp(𝑧) restricted on the interval [−1,∞) to be 𝑊(0,𝑧)≜𝑊(𝑧). Similarly, we define the inverse function of 𝑧exp(𝑧) restricted on the interval (−∞,−1] to be 𝑊(−1,𝑧)≜𝑊(𝑧). For this study, both 𝑊(0,𝑧) and 𝑊(−1,𝑧) will be employed only for 𝑧∈[−exp(−1),0) due to our practical problem (see Figure 1). For more details of the concepts and properties of the Lambert 𝑊 function, see Corless et al. [43].

Figure 1 

The property of Lambert 𝑊 function.

Lemma 4. Assume that there exists a bounded closed region 𝐺=𝐴𝐵𝐶𝐷𝐴, as shown in Figure 2. The boundaries 𝐴𝐷 and 𝐵𝐶 are the no cut-arcs of model (5), and on 𝐴𝐷 and 𝐵𝐶, the direction of the orientation field determined by model (5) is pointing to the interior of 𝐺. In addition, there is no singularity in the interior of 𝐺, and the boundary. One of the boundary 𝐶𝐷 is the impulse set of the model (5); the corresponding phase set satisfies 𝐼(𝐶𝐷)⊂𝐴𝐵; 𝐴𝐵 is also the no cut-arcs of model (5), and on 𝐴𝐵, the direction of the orientation field determined by model (5) is pointing to the interior of 𝐺. Then, there must exist at least one period-1 solution of model (5) in region 𝐺.

Figure 2 

Illustration of Bendixson field.

Lemma 5 (analogue of Poincaré’ criterion). The 𝑇-periodic solution 𝑆=𝜉(𝑡), 𝑥=𝜂(𝑡) of system 𝑑𝑆(),𝑑𝑡=𝑃𝑆,𝑥𝑑𝑥()()()()()𝑑𝑡=𝑅𝑆,𝑥,if𝜙𝑆,𝑥≠0,Δ𝑆=𝛼𝑆,𝑥,Δ𝑥=𝛽𝑆,𝑥,if𝜙𝑆,𝑥=0, (5) is orbitally asymptotically stable and enjoying the property of asymptotic phase if the multiplier 𝜇₂ satisfies the condition |𝜇₂|<1 and unstable if |𝜇₂|>1, where 𝜇₂=^𝑞_𝑘=1Δ_𝑘exp^𝑇₀𝜕𝑃𝜉(𝑡)(𝑡)+𝜕𝑆,𝜂𝜕𝑅𝜉(𝑡)(𝑡),Δ𝜕𝑥,𝜂𝑑𝑡_𝑘=𝑃₊𝜕𝛽/𝜕𝑋𝜕𝜙/𝜕𝑆−𝜕𝛽/𝜕𝑆𝜕𝜙/𝜕𝑥+𝜕𝜙/𝜕𝑆𝑃+𝑅𝜕𝜙/𝜕𝑆+𝑅𝜕𝜙/𝜕𝑥₊(𝜕𝛼/𝜕𝑆)𝜕𝜙/𝜕𝑥−(𝜕𝛼/𝜕𝑥)𝜕𝜙/𝜕𝑆+𝜕𝜙/𝜕𝑥𝑃,𝜕𝜙/𝜕𝑆+𝑅𝜕𝜙/𝜕𝑥 (6)𝑃₊=𝑃(𝜉(𝜏⁺_𝑘),𝜂(𝜏⁺_𝑘)), 𝑅₊=𝑅(𝜉(𝜏⁺_𝑘),𝜂(𝜏⁺_𝑘)), and 𝑃,𝑅,𝜕𝛼/𝜕𝑆,𝜕𝛼/𝜕𝑥,𝜕𝛽/𝜕𝑆,𝜕𝛽/𝜕𝑥,𝜕𝜙/𝜕𝑆,𝜕𝜙/𝜕𝑥 are calculated at the point (𝜉(𝜏_𝑘),𝜂(𝜏_𝑘)).

3. Properties Analysis of the Bioprocess

We focus our discussion on the case 𝑆_in>𝜇⁻¹(𝐷), 𝑥_set<𝜆(𝑆_in−𝜇⁻¹(𝐷)), and (𝑆₀,𝑥₀)∈Ω_𝑇={(𝑆,𝑥)∣0<𝑆≤𝑆_in,0<𝑥<𝑥_set}.

3.1. Existence of Period-1 Solution

The line 𝑥=𝑥_set intersects the isoclinal line 𝑑𝑥/𝑑𝑡=0 at the point 𝐴, intersects the curve 𝑑𝑆/𝑑𝑡=0 at the point 𝑀(𝑆_𝑀,𝑥_set), and intersects the line 𝑆=𝑆_in at the point 𝐵, where 𝑆_𝑀 satisfies that 𝜆𝐷𝑆_𝑀+𝜇𝑆_𝑀𝑥_set=𝜆𝐷𝑆_in. (7) The line 𝑥=(1−𝑊_𝑓)𝑥_set intersects the line 𝑑𝑥/𝑑𝑡=0 at the point 𝐸 and intersects the line 𝑆=𝑆_in at the point 𝐻. Denote 𝐶=𝐶(𝑆_𝐶,(1−𝑊_𝑓)𝑥_set) and 𝐷=𝐷(𝑆_𝐷,(1−𝑊_𝑓)𝑥_set), where 𝑆_𝐶=(1−𝑊_𝑓)𝜇⁻¹(𝐷) and 𝑆_𝐷=(1−𝑊_𝑓)𝑆_in. The impulsive set 𝑀_imp⊆𝐴𝐵={(𝑆,𝑥)∣𝜇⁻¹(𝐷)≤𝑆≤𝑆_in,𝑥=𝑥_set}, and the phase set 𝑁_pha=𝐈(𝑀_imp)⊆𝐶𝐷={(𝑆,𝑥)∣𝑆_𝐶≤𝑆≤𝑆_𝐷,𝑥=(1−𝑊_𝑓)𝑥_set}. Let the point 𝐾^′ be the primary image of 𝐸, and denote the trajectories of the system starting from the points 𝐶 and 𝐸 that intersect the segment 𝐴𝐵 at 𝐶^′ and 𝐸^′ respectively. Denote the trajectory of the system between 𝐶 and 𝐶^′ by 𝐶𝐶^′, the trajectory of the system between 𝐸 and 𝐸^′ by 𝐸𝐸^′.

Theorem 6. If one of the conditions holds (i)  𝑆_𝐷≤𝜇⁻¹(𝐷)<𝑆_𝑖𝑛 and 𝑥_𝑠𝑒𝑡<𝜆(𝑆_𝑖𝑛−𝜇⁻¹(𝐷)); (ii) If 𝜇⁻¹(𝐷)<𝑆_𝐷, 𝑆_𝐶′≤𝑆_𝐾′ and 𝑥_𝑠𝑒𝑡<𝜆(𝑆_𝑖𝑛−𝜇⁻¹(𝐷)); (iii) If 𝜇⁻¹(𝐷)<𝑆_𝐷, 𝑆_𝐸′≥𝑆_𝐾′ and 𝑥_𝑠𝑒𝑡<𝜆(𝑆_𝑖𝑛−𝜇⁻¹(𝐷)), then system (2) has a period-1 solution with (𝑆₀,𝑥₀)∈Ω_𝑇.

Proof. Firstly, it is obvious that all trajectories of system (2) starting from the region Ω_𝑇 must interest with the segment 𝐴𝐵 and then jump to the segment 𝐶𝐷 if 𝜇⁻¹(𝐷)<𝑆_in and 𝑥_set<𝜆(𝑆_in−𝜇⁻¹(𝐷)). Next, we construct the closed region 𝐺⊂Ω_𝑇 such that all solutions of system (2) starting from Ω_𝑇 enter into 𝐺 and retain there.
Case (i) (𝑆_𝐷≤𝜇⁻¹(𝐷)). It is obvious that the straight line 𝐴𝐶 is described by 𝑥=𝑥_set𝑆/𝜇⁻¹(𝐷), and 𝐵𝐷 is described by 𝑥=𝑥_set𝑆/𝑆_in, and both lines pass through the point 𝑂(0,0). The derivative of the straight line 𝐿(𝑆,𝑥) passing through the points 𝐴(𝜇⁻¹(𝐷),𝑥_set) and 𝐷((1−𝑊_𝑓)𝑆_in,(1−𝑊_𝑓)𝑥_set) along the trajectories of system (2) is 𝑑𝐿𝑑𝑡₍₂₎=𝑑𝑥𝑑𝑡₍₂₎−𝑊_𝑓𝑥_set𝜇⁻¹(𝐷)−1−𝑊_𝑓𝑆_in𝑑𝑆𝑑𝑡₍₂₎<0, (8) in terms of 𝑑𝑥/𝑑𝑡|₍₂₎<0, 𝑑𝑆/𝑑𝑡|₍₂₎>0 for the segment 𝐴𝐷. From the qualitative characteristic of system (2), we know that 𝑑𝐿/𝑑𝑡|₍₂₎>0 for 𝐴𝐷 and 𝑑𝑥/𝑑𝑡|₍₂₎<0 for 𝐶𝐷. Besides, 𝑑𝑆/𝑑𝑡|₍₂₎<0 for 𝑆=𝑆_in. On the other hand, it follows from system (2) that 𝑥(𝑡)=0 and 𝑆(𝑡)=𝑆_in−(𝑆_in−𝑆₀)𝑒^−𝐷𝑡 for 𝑡∈(0,+∞). Especially, when 𝑆₀=𝑆_𝐶, 𝐹^→𝐺 is the semitrivial solution of system (2). Therefore, we have found a closed region 𝐺, the boundary of which consists of 𝐴𝐷, 𝐷𝐶, 𝐶𝐹, 𝐹𝐺, 𝐺𝐵, and 𝐵𝐴; see Figure 3(a).
Case (ii) (𝑆_𝐷>𝜇⁻¹(𝐷) and 𝑆_𝐶′≤𝑆_𝐾′). Similar to the discussion of Case (i), the boundary of 𝐺 can be constructed by 𝐸𝐶, 𝐶𝐶^′, 𝐶^′𝐴, and 𝐴𝐸; see Figure 3(b).
Case (iii) (𝑆_𝐷>𝜇⁻¹(𝐷) and 𝑆_𝐸′≥𝑆_𝐾′). In this case, the boundary of 𝐺 can be constructed by 𝐸𝐷, 𝐷𝐻, 𝐻𝐵, 𝐵𝐸^′, and 𝐸𝐸^′; see Figure 3(c).
In addition, there is no singularity in the region 𝐺. Therefore, it follows from the qualitative characteristics of the system that all the trajectories satisfying the conditions of theorem enter the closed region 𝐺 and retain there. Then it follows from Lemma 4 that system (2) has a period-1 solution.

(a)

(b)

(c)

Figure 3 

Illustration of system (2) when 𝜇−1(𝐷)<𝑆in and 𝑥set<𝜆(𝑆in−𝜇−1(𝐷)): (a) 𝑆𝐷≤𝜇−1(𝐷); (b) 𝑆𝐷>𝜇−1(𝐷) and 𝑆𝐶′≤𝑆𝐾′; (c) 𝑆𝐷>𝜇−1(𝐷) and 𝑆𝐸′≥𝑆𝐾′.

3.2. Position of Period-1 Solution

Suppose the period of the period-1 solution is 𝑇. Let 𝑆(𝑡)=𝜉(𝑡),  𝑥(𝑡)=𝜂(𝑡) be such a 𝑇-periodic solution. Then 𝜁(𝑡)=𝜆𝜉(𝑡)+𝜂(𝑡) is also 𝑇-periodic. Denote 𝜉₀=𝜉(0⁺)=𝜉(𝑇⁺),   𝜉₁=𝜉(𝑇),  𝜂₀=𝜂(0⁺)=𝜂(𝑇⁺), and 𝜂₁=𝜂(𝑇). Then from the 𝑇-periodicity of the periodic solution and the third and fourth equations of system (2), we have 𝜉₀=(1−𝑊_𝑓)𝜉₁,  𝜂₀=(1−𝑊_𝑓)𝜂₁, and 𝜂₁=𝑥_set. Without loss of generality, for 𝑡∈[0,𝑇], 𝜁(𝑡) satisfies the relation 𝑑𝜁𝑑𝑡=𝐷𝜆𝑆_in−𝜁. (9) Then, we have 𝜁(𝑡)=𝜆𝑆_in+𝜁0⁺−𝜆𝑆_inexp(−𝐷𝑡). (10) In particular, for 𝑡=𝑇, we have 𝜁(𝑇)=𝜆𝑆_in(1−exp(−𝐷𝑇))+𝜁(0⁺)exp(−𝐷𝑇).

By (10), if 𝜁(0⁺)≥𝜆𝑆_in, then we have 𝜆𝑆_in≤𝜁(𝑡)<𝜆𝑆_in+(𝜁(0⁺)−𝜆𝑆_in)=𝜁(0⁺) for all 𝑡∈(0,𝑇], else we have 𝜆𝑆_in≥𝜁(𝑡)>𝜁(0) for all 𝑡∈(0,𝑇].

Let 𝑆_𝑅>𝜇⁻¹(𝐷) be the root of the following equation: 𝑆_𝑅−^𝑆𝑅_{𝜇−1(𝐷)}𝐷𝜇(𝜗)𝑑𝜗=𝜇⁻¹(𝐷)−𝑥_set𝜆𝑊_𝑓+ln1−𝑊_𝑓. (11) Then, we have the following results.

Theorem 7. Suppose that 𝑆_in>𝜇⁻¹(𝐷) and 𝑥_𝑠𝑒𝑡≤𝜆(𝑆_𝑖𝑛−𝜇⁻¹(𝐷)). Then the initial values of period-1 solution (𝜉(𝑡),𝜂(𝑡)) of system (2) satisfies (1−𝑊_𝑓)min{𝑆_𝑅,𝑆_𝑀}≤𝜉(0⁺)<(1−𝑊_𝑓)(𝑆_𝑖𝑛−𝑥_𝑠𝑒𝑡/𝜆) and 𝜂(0⁺)=(1−𝑊_𝑓)𝑥_𝑠𝑒𝑡, where 𝑆_𝑀 and 𝑆_𝑅 are determined by (7) and (11), respectively.

Proof. Since 𝜁(0⁺)=𝜆𝜉(0⁺)+𝜂(0⁺)=𝜆(1−𝑊_𝑓)𝜉(𝑇)+(1−𝑊_𝑓)𝜂(𝑇)=(1−𝑊_𝑓)(𝜆𝜉(𝑇)+𝜂(𝑇))=(1−𝑊_𝑓)𝜁(𝑇). Then, 𝜁(𝑇)=𝜆𝑆_in(1−exp(−𝐷𝑇))+𝜁(0⁺)exp(−𝐷𝑇)=𝜆𝑆_in(1−exp(−𝐷𝑇))+(1−𝑊_𝑓)𝜁(𝑇)exp(−𝐷𝑇), which implies that 𝜁(𝑇)=1−exp(−𝐷𝑇)1−1−𝑊_𝑓exp(−𝐷𝑇)𝜆𝑆_in<𝜆𝑆_in, (12) that is, 𝜉₁=𝜉(𝑇)<𝑆_in−𝑥_set/𝜆.
Next, we consider the following comparison system of system (2) without impulsive feedback control: 𝑑𝑆=1𝑑𝑡𝜆𝜇(𝑆)𝑥_set−𝑥,𝑑𝑥=𝜇(𝑆)𝑑𝑡−𝐷𝑥. (13) It is easy to verify that system (13) has one positive equilibrium (𝜇⁻¹(𝐷),𝑥_set), which is a center point. Then all of its solutions are closed trajectories which satisfy 𝑉₁(𝑆,𝑥)≜𝑥−𝑥_setln𝑥+𝜆𝑆−𝜆^𝑆_𝑆𝑀𝐷𝜇(𝜗)𝑑𝜗−𝑥_set+𝑥_setln𝑥_set−𝜆𝑆_𝑀=𝐶, (14) where 𝐶 is an arbitrary constant. The derivative of 𝑉₁(𝑆,𝑥) along the trajectories of system (2) is 𝑑𝑉₁𝑑𝑡₍₂₎=𝜕𝑉₁𝜕𝑥𝑑𝑥+𝑑𝑡𝜕𝑉₁𝜕𝑆𝑑𝑆=𝑥𝑑𝑡1−^set𝑥𝜇(𝑆)𝑥𝐷−𝐷+𝜆1−𝜇(𝑆)𝐷𝑆_in−1−𝑆𝜆𝜇(𝑆)𝑥=𝜆𝜇(𝑆)−𝐷𝜇(𝑆)𝐷𝑆_in−1−𝑆𝜆𝜇(𝑆)𝑥_set>0, (15) for 𝜇⁻¹(𝐷)<𝑆<𝑆_𝑀. This implies that the trajectories of system (2) intersecting with the trajectories of system (13) pass through the trajectories of system (13) from the left to the right. In addition, the trajectory of system (13) passing through 𝑀(𝑆_𝑀,𝑥_set) (i.e., 𝐶=0 in (14)) intersects the line 𝑆=𝜇⁻¹(𝐷) at two points 𝐿₁(𝜇⁻¹(𝐷),𝑥_𝐿1) and 𝐿₂(𝜇⁻¹(𝐷),𝑥_𝐿2), where 𝑥_𝐿𝑖(𝑖=1,2) satisfy that ln𝑥_𝐿𝑖−𝑥_𝐿𝑖𝑥_set=Λ₁≜𝜆𝜇⁻¹(𝐷)−𝜆∫^𝜇−1_𝑆^(𝐷)_𝑀𝐷/𝜇(𝜗)𝑑𝜗−𝑥_set+𝑥_setln𝑥_set−𝜆𝑆_𝑀𝑥_set, (16)that is, −𝑥_𝐿𝑖𝑥_set𝑒^{−𝑥𝐿𝑖/𝑥set}𝑒=−^Λ1𝑥_set. (17) Since Λ₁=ln(𝑥_𝐿𝑖)−𝑥_𝐿𝑖/𝑥_set<ln𝑥_set−1, then we have −𝑒^Λ1/𝑥_set>−1/𝑒. Thus, we have 𝑥_𝐿1=−𝑥_set𝑊𝑒−1,−^Λ1𝑥_set,𝑥_𝐿2=−𝑥_set𝑊𝑒0,−^Λ1𝑥_set, (18) and 𝑥_𝐿2<𝑥_set<𝑥_𝐿1, where 𝑊 denotes Lambert 𝑊 function defined in Definition 3.
If 𝑥_𝐿2≥(1−𝑊_𝑓)𝑥_set, then the closed region 𝐺 constructed in Theorem 6 can be reduced to 𝐺^′, the boundary of which consists of 𝑀𝐿₂, 𝐿₂𝐸, 𝐸𝐶, 𝐶𝐹, 𝐹𝐺, 𝐺𝐵, and 𝐵𝑀. Thus, the period-1 solution satisfies that 𝜉₁=𝜉(𝑇)≥𝑥_𝑀>𝜇⁻¹(𝐷). Else, we have 𝑥_𝐿2<(1−𝑊_𝑓)𝑥_set. Then, the trajectory of system (13) passing through 𝐸 intersects the segment 𝐴𝑀 at the point 𝑅(𝑆_𝑅,𝑥_set). In this case, the closed region 𝐺 constructed in Theorem 6 can be reduced to 𝐺^′, the boundary of which consists of 𝑅𝐿₂, 𝐿₂𝐸, 𝐸𝐶, 𝐶𝐹, 𝐹𝐺, 𝐺𝐵, and 𝐵𝑅. Thus, the period-1 solution satisfies that 𝜉₁=𝜉(𝑇)≥𝑆_𝑅>𝜇⁻¹(𝐷).

Corollary 8. Suppose that 𝜇⁻¹(𝐷)<𝑆_𝑖𝑛≤𝜇⁻¹(𝐷)/(1−𝑊_𝑓) and 𝑥_𝑠𝑒𝑡≤𝜆(𝑆_𝑖𝑛−𝜇⁻¹(𝐷)). Then, the initial values of the period-1 solution (𝜉(𝑡),𝜂(𝑡)) of system (2) satisfies that 𝜂(0⁺)=(1−𝑊_𝑓)𝑥_𝑠𝑒𝑡 and (1−𝑊_𝑓)min{𝑆_𝑅,𝑆_𝑀}≤𝜉(0⁺)<𝜇⁻¹(𝐷).

Corollary 9. Suppose that 𝑆_𝑖𝑛>𝜇⁻¹(𝐷)/(1−𝑊_𝑓) and max{𝜆[𝑆_𝑖𝑛−𝜇⁻¹(𝐷)/(1−𝑊_𝑓)],0}≤𝑥_𝑠𝑒𝑡≤𝜆(𝑆_𝐹−𝜇⁻¹(𝐷)). Then the initial values of the period-1 solution (𝜉(𝑡),𝜂(𝑡)) of system (2) satisfies that (1−𝑊_𝑓)min{𝑆_𝑅,𝑆_𝑀}≤𝜉(0⁺)<𝜇⁻¹(𝐷) and 𝜂(0⁺)=(1−𝑊_𝑓)𝑥_𝑠𝑒𝑡.

Corollary 10. Suppose that 𝑆_𝑖𝑛>𝜇⁻¹(𝐷), 𝑥_𝑠𝑒𝑡≤𝜆(𝑆_𝐹−𝜇⁻¹(𝐷)) and 𝑊_𝑓≥1−𝜇⁻¹(𝐷)/(𝑆_𝑖𝑛−𝑥_𝑠𝑒𝑡/𝜆). Then, the initial values of the period-1 solution (𝜉(𝑡),𝜂(𝑡)) of system (2) satisfies that (1−𝑊_𝑓)min{𝑆_𝑅,𝑆_𝑀}≤𝜉(0⁺)<𝜇⁻¹(𝐷) and 𝜂(0⁺)=(1−𝑊_𝑓)𝑥_𝑠𝑒𝑡.

3.3. Stability of Period-1 Solution

Let 𝑆^∗=𝜇⁻¹((2−𝑊_𝑓)𝐷/(1+(1−𝑊_𝑓)²)). Denote 𝑥^∗=𝜆𝐷𝑆_in−𝑆^∗𝜇(𝑆^∗),𝑥=−𝑥^∗𝑊𝑒0,−^Λ2𝑥^∗, (19) where Λ₂≜[𝜆𝜇⁻¹(𝐷)−𝜆∫^𝜇−1_𝑆^(𝐷)_∗(𝐷/𝜇(𝜗))𝑑𝜗−𝑥_set+𝑥^∗ln𝑥_set−𝜆𝑆^∗]/𝑥^∗, and 𝑊 denotes Lambert 𝑊 function defined in Definition 3. Then, we have the following result.

Theorem 11. Suppose that 𝑆_𝑖𝑛>𝜇⁻¹(𝐷), 𝑥_𝑠𝑒𝑡≤min{𝑥^∗,𝑥/(1−𝑊_𝑓)}. Then, the period-1 solution (𝜉(𝑡),𝜂(𝑡)) of system (2) is orbitally asymptotically stable and enjoys the property of asymptotic phase.

Proof. Let 𝑆^∗=𝜇⁻¹((2−𝑊_𝑓)𝐷/(1+(1−𝑊_𝑓)²)). Denote 𝑥^∗=𝜆𝐷𝑆_in−𝑆^∗𝜇(𝑆^∗),𝑥=−𝑥^∗𝑊𝑒0,−^Λ2𝑥^∗, (20) where Λ₂≜[𝜆𝜇⁻¹(𝐷)−𝜆∫^𝜇−1_𝑆^(𝐷)_∗(𝐷/𝜇(𝜗))𝑑𝜗−𝑥_set+𝑥^∗ln𝑥_set−𝜆𝑆^∗]/𝑥^∗, and 𝑊 denotes Lambert 𝑊 function defined in Definition 3. According to Lemma 5, we calculate the multiplier of system (2) in variations corresponding to the 𝑇-periodic solution (𝜉(𝑡),𝜂(𝑡)). Denote 𝑃₀(𝜉₀,𝜂₀), 𝑃₁(𝜉₁,𝜂₁), where 𝜉₀=(1−𝑊_𝑓)𝜉₁, 𝜂₀=(1−𝑊_𝑓)𝑥_set, 𝜇⁻¹(𝐷)<𝜉₁<𝑆_in−𝑥_set/𝜆, and 𝜂₁=𝑥_set. In system (2), since 𝑃(𝑆,𝑥)=𝐷(𝑆_in−𝑆)−𝜇(𝑆)𝑥/𝜆, 𝑅(𝑆,𝑥)=(𝜇(𝑆)−𝐷)𝑥, 𝛼(𝑆,𝑥)=−𝑊_𝑓𝑆, 𝛽(𝑆,𝑥)=−𝑊_𝑓𝑥, and 𝜙(𝑆,𝑥)=𝑥−𝑥_set, then 𝜕𝑃1𝜕𝑆=−𝐷−𝜆𝜇^′(𝑆)𝑥,𝜕𝑅(𝑆)𝜕𝑥=𝜇−𝐷,𝜕𝛼𝜕𝑆=−𝑊_𝑓,𝜕𝛼𝜕𝑥=0,𝜕𝛽𝜕𝑆=0,𝜕𝛽𝜕𝑥=−𝑊_𝑓,𝜕𝜙𝜕𝑆=0,𝜕𝜙𝜕𝑥=1. (21) Therefore, Λ₁=𝑃₊𝜕𝛽/𝜕𝑥𝜕𝜙/𝜕𝑆−𝜕𝛽/𝜕𝑆𝜕𝜙/𝜕𝑥+𝜕𝜙/𝜕𝑆𝑃+𝑅𝜕𝜙/𝜕𝑆+𝑅𝜕𝜙/𝜕𝑥₊(𝜕𝛼/𝜕𝑆)𝜕𝜙/𝜕𝑥−(𝜕𝛼/𝜕𝑥)𝜕𝜙/𝜕𝑆+𝜕𝜙/𝜕𝑥𝑃=𝑅𝜕𝜙/𝜕𝑆+𝑅𝜕𝜙/𝜕𝑥₊1−𝑊_𝑓𝑅=1−𝑊_𝑓²𝜇𝜉₀−𝐷𝜇𝜉₁,𝜇−𝐷₂=Λ₁exp^𝑇₀𝜕𝑃𝜉(𝑡)(𝑡)+𝜕𝑆,𝜂𝜕𝑅𝜉(𝑡)(𝑡)𝜕𝑥,𝜂𝑑𝑡=Λ₁exp^𝑇₀1−𝐷−𝜆𝜇^′(𝑆)𝑥𝑑𝑡+^𝑇₀𝜇(𝑆)−𝐷𝑑𝑡=Λ₁exp^𝜂1_𝜂01𝑥𝑑𝑥−^𝑇₀1𝐷+𝜆𝜇^′(𝑆)𝑥=𝑑𝑡1−𝑊_𝑓𝜇𝜉₀−𝐷𝜇𝜉₁−−𝐷exp^𝑇₀1𝐷+𝜆𝜇^′(𝑆)𝑥𝑑𝑡. (22) Note that 𝜉₁=𝜉(𝑇)≥min{𝑆_𝑀,𝑆_𝑅}>𝜇⁻¹(𝐷), that is, 𝜇(𝜉₁)>𝐷. If 𝜉₀=(1−𝑊_𝑓)𝜉₁≥𝜇⁻¹(𝐷), that is, the periodic solution lies in the region, where 𝑑𝑥/𝑑𝑡>0, then we have 0≤𝜇₂<1, and the periodic solution (𝜉(𝑡),𝜂(𝑡)) is orbitally asymptotically stable and enjoying the property of asymptotic phase. But for 𝜉₀<𝜇⁻¹(𝐷), that is, 𝜉₁<𝜇⁻¹(𝐷)/(1−𝑊_𝑓), we have 𝜇(𝜉₀)−𝐷<0. Since the expression of 𝜉₁ is hard to be given, we cannot determine whether the periodic solution is stable or not. However, we can give some sufficient conditions to guarantee −1<𝜇₂<0 when 𝜉₁<𝜇⁻¹(𝐷)/(1−𝑊_𝑓), that is, the periodic solution is orbitally asymptotically stable and enjoying the property of asymptotic phase.
Since 0<𝑊_𝑓<1, then we have 1<(2−𝑊_𝑓)/(1+(1−𝑊_𝑓)²)<1/(1−𝑊_𝑓), that is, 𝐷<(2−𝑊_𝑓)𝐷/(1+(1−𝑊_𝑓)²)=𝜇(𝑆^∗)<𝐷/(1−𝑊_𝑓), which implies that 𝜇⁻¹(𝐷)<𝑆^∗<𝜇⁻¹(𝐷)/(1−𝑊_𝑓). If 𝜉₁≥𝑆^∗, then we have 𝜇(𝜉₁)+(1−𝑊_𝑓)𝜇(𝜉₀)≥(1+(1−𝑊_𝑓)²)𝜇(𝜉₁)≥(2−𝑊_𝑓)𝐷, which is equivalent to (1−𝑊_𝑓)(𝐷−𝜇(𝜉₀))≤𝜇(𝜉₁)−𝐷 or (1−𝑊_𝑓)(𝐷−𝜇(𝜉₀))/(𝜇(𝜉₁)−𝐷)≤1. Thus, we get |𝜇₂|<1, that is, the periodic solution is orbitally asymptotically stable and enjoying the property of asymptotic phase. In the following, we will give the conditions under which 𝜉₁≥𝑆^∗.
From the expression of the isoclinal line 𝑑𝑆/𝑑𝑡=0, it follows that for 𝑆=𝑆^∗, we have 𝑥=𝜆𝐷(𝑆_in−𝑆^∗)/𝜇(𝑆^∗)∶=𝑥^∗.
Consider the following comparison system of system (2) without impulsive feedback control 𝑑𝑆=1𝑑𝑡𝜆𝜇(𝑆)𝑥^∗−𝑥,𝑑𝑥=𝜇(𝑆)𝑑𝑡−𝐷𝑥. (23) It is easy to verify that system (23) has one positive equilibrium (𝜇⁻¹(𝐷),𝑥^∗), which is a center point. Then, all of its solutions are closed trajectories which satisfy 𝑉₂(𝑆,𝑥)≜𝑥−𝑥^∗ln𝑥+𝜆𝑆−𝜆^𝑆_𝑆∗𝐷𝜇(𝜗)𝑑𝜗−𝑥_set+𝑥^∗ln𝑥_set−𝜆𝑆^∗=𝐶, (24) where 𝐶 is an arbitrary constant. The derivative of 𝑉₂(𝑆,𝑥) along the trajectories of system (2) is 𝑑𝑉₂𝑑𝑡₍₂₎=𝜕𝑉₂𝜕𝑥𝑑𝑥+𝑑𝑡𝜕𝑉₂𝜕𝑆𝑑𝑆=𝑥𝑑𝑡1−^∗𝑥𝜇(𝑆)𝑥𝐷−𝐷+𝜆1−𝜇(𝑆)𝐷𝑆_in−1−𝑆𝜆𝜇(𝑆)𝑥=𝜆𝜇(𝑆)−𝐷𝜇(𝑆)𝐷𝑆_in−1−𝑆𝜆𝜇(𝑆)𝑥^∗>0, (25) for 𝜇⁻¹(𝐷)<𝑆<𝑆^∗. This implies that the trajectories of system (2) intersecting with the trajectories of system (23) pass through the trajectories of system (23) from the left to the right. In addition, the trajectory of system (23) passing through 𝑁(𝑆^∗,𝑥_set) (i.e., 𝐶=0 in (24)) intersects the line 𝑆=𝜇⁻¹(𝐷) at two points 𝐾₁(𝜇⁻¹(𝐷),𝑥_𝐾1) and 𝐾₂(𝜇⁻¹(𝐷),𝑥_𝐾2), where 𝑥_𝐾𝑖(𝑖=1,2) satisfy thatln𝑥_𝐾𝑖−𝑥_𝐾𝑖𝑥^∗=Λ₂≜𝜆𝜇⁻¹(𝐷)−𝜆∫^𝜇−1_𝑆^(𝐷)_∗𝐷/𝜇(𝜗)𝑑𝜗−𝑥_set+𝑥^∗ln𝑥_set−𝜆𝑆^∗𝑥^∗, (26)that is, −𝑥_𝐾𝑖𝑥^∗𝑒^{−𝑥𝐾𝑖/𝑥∗}𝑒=−^Λ2𝑥^∗. (27) Since Λ₂=ln(𝑥_𝐾𝑖)−𝑥_𝐾𝑖/𝑥^∗<ln𝑥^∗−1, then we have −𝑒^Λ2/𝑥^∗>−1/𝑒. Thus, by Definition 3, we have 𝑥_𝐾1=−𝑥^∗𝑊𝑒−1,−^Λ2𝑥^∗,𝑥_𝐾2=−𝑥^∗𝑊𝑒0,−^Λ2𝑥^∗, (28) and 𝑥_𝐾2<𝑥^∗<𝑥_𝐾1.
Next, we will prove that 𝜉₁≥𝑆^∗ when the conditions of theorem are satisfied. Let 𝑁 be the intersection point of the line 𝑥=𝑥^∗ and 𝑑𝑆/𝑑𝑡=0, 𝐾₁, and let 𝐾₂ be the intersection points of the line 𝑆=𝜇⁻¹(𝐷) and the trajectory of system (23) passing through the point 𝑀. If 𝑥_set≤min{𝑥^∗,𝑥/(1−𝑊_𝑓)}, then we have (1−𝑊_𝑓)𝑥_set≤𝑥. Since 𝑑𝑉₂/𝑑𝑡|₍₂₎>0, 𝑑𝑆/𝑑𝑡|_𝐾1𝐸>0, 𝑑𝑥/𝑑𝑡|_𝐸𝐶<0, 𝑑𝑆/𝑑𝑡|_𝐶𝐹>0, 𝑑𝑆/𝑑𝑡|_𝐺𝐵>0, and 𝐹^→𝐺 is semi-trivial solution of system (2). Therefore, we get a closed region 𝐺, the boundary of which consists of 𝑁𝐾₁, 𝐾₁𝐸, 𝐸𝐶, 𝐶𝐹, 𝐹𝐺, 𝐺𝐵, and 𝐵𝑁. In addition, there is no singularity in the region 𝐺. Therefore, it follows from the qualitative characteristics of the system that all the trajectories satisfying the conditions of theorem enter the closed region 𝐺 and retain there. Then, it follows from Theorem 6 that system (2) has a period-1 solution in 𝐺. Furthermore, we have 𝑆^∗≤𝜉₁≤𝑆_in. If 𝜉₁≥𝜇⁻¹(𝐷), then 𝜉₀≥𝜇⁻¹(𝐷), in this case, we have 0≤𝜇₂<1. Else 𝑆^∗≤𝜉₁<𝜇⁻¹(𝐷)/(1−𝑊_𝑓), in this case, we have −1<𝜇₂<0. Therefore, if the conditions of the theorem are satisfied, then the period-1 solution (𝜉(𝑡),𝜂(𝑡)) of system (2) is orbitally asymptotically stable and enjoying the property of asymptotic phase.

3.4. Discussion on the Period-𝑘 (𝑘≥2) Solution

For the predator-prey model concerning IPM strategies, Tang and Cheke [23] proved that there is no periodic solution with order larger than or equal to three, except for one special case, by using the properties of the Lambert 𝑊 function and Poincaré map. Moreover, They showed that the existence of an order two periodic solution implies the existence of an order one periodic solution. Next, the existence of the period-𝑘 (𝑘≥3) solution is discussed.

Theorem 12. Suppose that 𝑆_𝑖𝑛>𝜇⁻¹(𝐷)/(1−𝑊_𝑓) and 𝜆(𝑆_𝑖𝑛−𝜇⁻¹(𝐷)/(1−𝑊_𝑓))≤𝑥_𝑠𝑒𝑡<𝜆(𝑆_𝑖𝑛−𝜇⁻¹(𝐷)). Then, system (2) does not have period-𝑘 (𝑘≥3) solution.

Proof. Let (𝑆,𝑥) be an arbitrary solution of system (2). Denote the first intersection point of the trajectory to the impulsive set 𝑀 by (𝑆₁,𝑥_set) (here without loss of generality, we assume that 𝑆₁>𝜉₁), and the corresponding consecutive points are (𝑆₂,𝑥_set), (𝑆₃,𝑥_set), …, respectively. Consequently, under the effect of impulsive function 𝐼, the corresponding points after pulse are (𝑆⁺₁,(1−𝑊_𝑓)𝑥_set),(𝑆⁺₂,(1−𝑊_𝑓)𝑥_set),(𝑆⁺₃,(1−𝑊_𝑓)𝑥_set),…. By the qualitative analysis of system (2), we know that 𝑥̇>0 for 𝑆>𝜇⁻¹(𝐷) and 𝑥̇<0 for 𝑆<𝜇⁻¹(𝐷). If there exists the 𝑗th impulse effect such that 𝑆⁺_𝑘≥𝜇⁻¹(𝐷) for 𝑘≥𝑗, then the sequence will be monotone, that is, 𝑆⁺_𝑗≥𝑆⁺_𝑗+1≥⋯≥𝜇⁻¹(𝐷) or 𝜇⁻¹(𝐷)≤𝑆⁺_𝑗≤𝑆⁺_𝑗+1≤⋯. In this case, system (2) has no period-𝑘(𝑘≥2) solution. Let 𝑙₁ be the line passing the point 𝐺(𝑆_in,0) with the slope −𝜆: 𝑙₁∶𝑥+𝜆𝑆−𝜆𝑆_in=0. (29) Since 𝑑𝑙₁/𝑑𝑡|₍₂₎=𝑑𝑥/𝑑𝑡+𝜆𝑑𝑆/𝑑𝑡=𝐷(𝜆𝑆_in−𝜆𝑆−𝑥)=0, then the trajectory starting from the point on the segment 𝐶𝐷 either does not intersect 𝑙₁ or retain on 𝑙₁ before reaching the segment 𝐴𝐵. Therefore, when 𝑥_set≥𝜆(𝑆_in−𝜇⁻¹(𝐷)/(1−𝑊_𝑓)), the period-1 solution determined in Theorem 6 satisfies that (1−𝑊_𝑓)min{𝑆_𝑅,𝑆_𝑀}≤𝜉₀<(1−𝑊_𝑓)(𝑆_in−𝑥_set/𝜆)≤𝜇⁻¹(𝐷), which also implies that all 𝑆⁺_𝑘≤𝜇⁻¹(𝐷).
The tr
ajectory of the solution (𝑆,𝑥) jumps to the point (𝑆⁺₁,(1−𝑊_𝑓)𝑥_set) from the point (𝑆₁,𝑥_set). By the dynamics of system (2), we have the following two sequences according to (𝑆₃,𝑥_set) lies on the right or left of (𝑆₁,𝑥_set): (a)  ⋯≤𝑆₄≤𝑆₂<𝜉₁<𝑆₁≤𝑆₃≤⋯; (b)  𝑆₂≤𝑆₄≤⋯≤𝜉₁≤⋯≤𝑆₃≤𝑆₁. Since the subsequences {𝑆_2𝑘−1} and {𝑆_2𝑘}   (𝑘=1,2,…) are monotone and bounded, then there exists a limit for them, respectively. Let 𝑆=lim_𝑘→∞𝑆_2𝑘 and 𝑆=lim_𝑘→∞𝑆_2𝑘−1. When sequence (a) occurs, it can be shown that 𝑆<𝜉₁<𝑆. Therefore, system (2) has a period-2 solution. When sequence (b) occurs, we have 𝑆≤𝜉₁≤𝑆. If 𝑆<𝜉₁<𝑆, then system (2) has a period-2 solution. Else we have 𝑆=𝜉₁=𝑆, in this case, the trajectory tends to the period-1 solution (𝜉,𝜂). But it can be obtained that there is no order 𝑘 solution (𝑘≥3) in system (2), thus the system is not chaos by Li-Yorke chaos Theorem [44].

Remark 13. From the proof of Theorem 12, it can be seen that there is a possibility to exist a period-2 solution. Since it is difficult to obtain the exact expression of the solution to (2), then it will be challenging to give the sufficient conditions under which system (2) has a period-2 solution or not.

4. Applied Instance

Here, we will take the Tessier kinetics model [41], that is, 𝜇(𝑆)=𝜇_max(1−exp(−𝑆/𝐾_𝑆)), as the example of specific growth rate to verify the theoretical results. In addition, 𝜇_max=0.5[1/h], 𝐾_𝑆=0.2[g/L], and 𝜆=𝑌_𝑥/𝑆max=0.5[g/g] are used for the demonstration of system behavior.

4.1. Numerical Simulations

The simulations are carried out by changing one main parameter and fixing all other parameters. The following numerical simulations are given by the programs of Maple and Matlab softwares. The time interval is set two days, that is, 48[h].

We assume in the following that 𝑆₀=𝑆(0)=0.2[g/L] and 𝑥₀=𝑥(0)=2[g/L]. The vector graph of system (1) for 𝐷=0.4[1/h] and 𝑆_in=20[g/L] are shown in Figure 4, from which it can be seen that the equilibrium (0.32,9.84) is a stable node.

Figure 4 

Illustration of vector graph of system (1) when 𝜇max=0.5[1/h], 𝐾𝑆=0.2[g/L], 𝜆=𝑌𝑥/𝑆max=0.5[g/g], 𝐷=0.4[1/h], and 𝑆in=20[g/L].

The time series and phase portrait for 𝑥_set=10[g/L] is shown in Figure 5, it can be seen that no impulse occurs when 𝑥_set>𝜆(𝑆_in+𝐾_𝑆ln(1−𝐷/𝜇_max))=9.84[g/L].

(a)

(b)

(c)

Figure 5 

The time series and phase portrait when 𝐷=0.4[1/ℎ], 𝑆in=20[g/L] and 𝑥set=10[g/L].

Figure 6 gives the time series and phase portrait when 𝑥_set=6[g/L] and 𝑊_𝑓=0.3. It can be seen that the trajectory tends to be periodic.

(a)

(b)

(c)

Figure 6 

The time series and phase portrait when 𝐷=0.4[1/h], 𝑊𝑓=0.3, 𝑆in=20[g/L], and 𝑥set=6[g/L].

Figure 7 shows the different positions of the period-1 solution for 𝐷=0.4[1/h] and 𝑊_𝑓=0.1,0.3,0.5.

(a)

(b)

(c)

Figure 7 

Dependence of the biomass and substrate concentration on 𝑊𝑓. The phase portraits when 𝐷=0.4[1/h], 𝑆in=20[g/L], 𝑥set=6[g/L], and 𝑊𝑓=0.1,0.3,0.5.

Figure 8 shows the different positions of the period-1 solution for 𝑊_𝑓=0.3, 𝑆_in=20[g/L], and 𝐷=0.3[1/h],0.4[1/h],0.45[1/ℎ].

(a)

(b)

(c)

Figure 8 

Dependence of the biomass and substrate concentration on 𝐷. The phase portraits when 𝑊𝑓=0.3, 𝑆in=20[g/L], 𝑥set=6[g/L], and 𝐷=0.3[1/h],0.4[1/h],0.45[1/h].

The numerical simulations are consistent with the theoretical results obtained and presented in Section 3.

4.2. Optimization of the Bioprocess

Next, we discuss aspects of the bioprocess optimization.

4.2.1. The Proposed Objective Function and the Constraints

(1)  The Objective Function𝑋_OUT𝐷,𝑊_𝑓,𝑆_in,𝑥_set=1𝑇^𝑇₀𝑊𝐷𝑥𝑑𝑡+^𝑓𝑥_set𝑇, (30) where 𝑋_OUT is the biomass productivity in the steady state.

(2)  The Constraints (a)0≤𝑊_𝑓≤𝑊_𝑓max<1 𝑊_𝑓—the part of “clear” medium inputted into the bioreactor in each biomass oscillation cycle,  𝑊_𝑓max—the maximum part of “clear” medium inputted into the bioreactor in each biomass oscillation cycle. (b)𝑆_min≤𝑆_in≤𝑆_CRI:  𝑆_in—the concentration of the inputted substrate,  𝑆_min—the minimum concentration of the inputted substrate,  𝑆_CRI—the critical level of the dosaged substrate concentration. (c)0<𝐷_min𝜑≤𝐷≤𝐷_max=𝜇_max(1−𝑒^{−𝑆in/𝐾𝑆}):  𝐷—the dilution rate of the chemostat,  𝐷_min𝜑—the minimum dilution rate. If the used microorganisms growth kinetics model satisfies regularity and monotonicity, then for 𝑆_in≫𝑆_min there is 𝜑=1,i.e.,for𝐷≥𝐷_min1,𝑌_𝑥/𝑆=𝑌_𝑥/𝑆max=𝜆,𝜑=2,thatis,for𝐷≥𝐷_min2,𝑌_𝑥/𝑆=𝑌_𝑥/𝑆max=𝜆,𝜇(𝑆)=𝜇_max, (31) where 𝐷_min2>𝐷_min1,  𝜇_max—the maximum specific growth rate,  𝐾_𝑆—the substrate saturation constant.(d)0<𝑥_set≤𝜆(𝑆_in+𝐾_𝑆ln(1−𝐷/𝜇_max))≤𝑥_CRI:  𝑥_set—the set level of the biomass concentration in the bioreactor medium,  𝑥_CRI—the critical level of biomass concentration in the bioreactor medium.

The aim of the optimization is to find the maximum of the objective function (30) under constraints (a)–(d), that is, 𝑋_OUT→max. (32)

4.2.2. The Analysis of the Optimization

(A) Case of a Chemostat with Impulse Effect. For 𝐷>𝐷_min2 and 𝑊_𝑓>0, a continuous bioprocess (a chemostat with impulse effect) is obtained, that is, 𝑑𝑆𝑆𝑑𝑡=𝐷_in−𝜇−𝑆^max𝜆𝑥𝑑𝑥=𝜇𝑑𝑡^max−𝐷𝑥𝑥<𝑥_set,Δ𝑆=−𝑊_𝑓𝑆Δ𝑥=−𝑊_𝑓𝑥𝑥=𝑥_set. (33) In the steady state, that is, for the period-1 solution (𝜉(𝑡),𝜂(𝑡)), by the second equation of system (33), we have 1𝑑𝑡=𝜇_max1−𝐷𝜂𝑑𝜂. (34) Then, travelling along 𝜂(𝜉) from the point 𝑃₀(𝜉₀,𝜂₀), with 𝑡=𝑡_𝑃0 to the point 𝑃₁(𝜉₁,𝜂₁), with 𝑡=𝑡_𝑃1 in the counterclockwise direction yields the period 𝑇, that is, 𝑇=𝑡_𝑃1−𝑡_𝑃0=^𝜂1_𝜂01𝜇_max1−𝐷𝜂𝑑𝜂=−ln1−𝑊_𝑓𝜇_max.−𝐷 (35) It is obviouse that 𝑇→0 when 𝑊_𝑓→0. In addition, by (34), for 𝑡∈(0,𝑇], we have 𝜂(𝑡)=(1−𝑊_𝑓)𝑥_set𝑒^{(𝜇max−𝐷)𝑡}, thus, ^𝑇₀(𝑡)𝐷𝐷𝜂𝑑𝑡=𝜇_max𝑊−𝐷^𝑓𝑥_set. (36) Therefore, the objective function can be formulated as 𝑋⁽²⁾_OUT𝐷,𝑊_𝑓,𝑆_in,𝑥_set=𝐷𝜇_max𝑊−𝐷^𝑓𝑥_set𝑇+𝑊_𝑓𝑥_set𝑇. (37) For given 𝑊_𝑓=𝑊_𝑓set>0 and 𝑥_set, we have 𝑋⁽²⁾_OUToptimal=lim_{𝐷→𝜇max}𝑋⁽²⁾_OUT𝐷,𝑊_𝑓set,𝑆_in,𝑥_set=lim_{𝐷→𝜇max}𝑥_set𝑊_𝑓set𝑇1−𝐷/𝜇_max𝜇→−^max𝑥_set𝑊_𝑓setln1−𝑊_𝑓set, (38) where 𝑥_set<𝜆𝑆_CRI and 𝑊_𝑓set∈(0,𝑊_𝑓max].

(B) Case of Synchronization Loss. In this case the optimization of the function (37) for any value of 𝑊_𝑓 is presented. We should find the maximum of the objective function (37) under the constraints 0≤𝑊_𝑓≤𝑊_𝑓max, 𝐷_min≤𝐷≤𝐷_max, 𝑆_min≤𝑆_in≤𝑆_CRI, and 0<𝑥_set<𝜆𝑆_in.

Let 𝜑(𝑊_𝑓)=−𝑊_𝑓/ln(1−𝑊_𝑓). Take the first derivative of 𝜑(𝑊_𝑓) with respect to 𝑊_𝑓, we have 𝜑^′𝑊_𝑓=−ln1−𝑊_𝑓+𝑊_𝑓/1−𝑊_𝑓ln²1−𝑊_𝑓. (39) Let 𝜓(𝑊_𝑓)=ln(1−𝑊_𝑓)−𝑊_𝑓/(1−𝑊_𝑓). Then, we have 𝜓(0)=0. Take the first derivative of 𝜓(𝑊_𝑓) with respect to 𝑊_𝑓, we have 𝜓𝑊_𝑓=𝑊_𝑓1−𝑊_𝑓²>0, (40) then 𝜓(𝑊_𝑓)≥𝜓(0)=0, which implies that 𝜑^′(𝑊_𝑓)≤0, that is, 𝜑(𝑊_𝑓) is decreasing on 𝑊_𝑓.

Therefore, for given 𝑥_set, we have 𝑋⁽³⁾_OUToptimal=lim_𝑊𝑓→0𝑋⁽²⁾_OUToptimal=lim_{𝐷→𝜇max,𝑊𝑓→0}𝑋⁽²⁾_OUT→𝜇_max𝑥_set, (41) where 𝑥_set<𝜆𝑆_CRI.

5. Conclusions

In the paper the dynamic behavior of a universal biochemical reaction process with feedback control was analyzed, and it was shown that the stability of the bioprocess (i.e., the existence of the positive period-1 solution), depended on both the biomass yield and the microorganism growth rate. Furthermore, the conditions for the existence and stability of the system's period-1 solution were obtained (i.e., Theorems 6 and 11). It also pointed out that the system (2) was not chaotic according to an analysis of period-2 solution existence. According to the theoretical results, the production of the microorganisms tended to be periodic. The microorganisms in the considered chemostat always kept the suitable growth rate. The biomass concentration could have been controlled to a given level between the critical and the optimal dissolved oxygen concentration. Moderation of the main bioprocess parameters in the selected biomass oscillation cycles for the Tessier kinetics model was presented.

After the analysis of the system stability, bioprocess optimization was covered in the work. In the first case, the optimization of a simple continuous bioprocess (i.e., bioprocess in which 𝐷>0 and 𝑊_𝑓=0) was presented. In the second case of optimization, it was demonstrated that theoretically it was possible to obtain almost optimal biomass productivity without lost of synchronization of simultaneously performed bioprocesses (see (38) for small values of 𝑊_𝑓set). This was possible, when the substrate was dosed by both continuous (i.e., 𝐷>0) and pulse (i.e., 𝑊_𝑓>0) method. In this case, the impulsive chemostat was obtained, where the biomass productivity mainly depended on 𝐷, and the synchronization of simultaneously performed bioprocesses depended on 𝑊_𝑓. The third case of optimization showed that during the optimization, the impulsive chemostat lost possibility of synchronization and strived for a simple continuous bioprocess (see (41)). Moreover, it was shown that improving the biomass productivity and synchronization of simultaneously performed bioprocesses were conflict criteria.

Acknowledgment

This research is supported in part by National Natural Science Foundation of China (Grant no. 11101066) and the Fundamental Research Funds for the Central Universities.