Analysis and constrained optimal impulsive control of a Holling-II type trophic system with two sources

Available online at www.sciencedirect.com

Journal of the Franklin Institute 352 (2015) 2728–2749 www.elsevier.com/locate/jfranklin

Analysis and constrained optimal impulsive control of a Holling-II type trophic system with two sources Luca Galbuseraa,n, Sara Pasqualib a

European Commission, DG Joint Research Centre (JRC), Institute for the Protection and Security of the Citizen (IPSC), Via E. Fermi 2749, 21027 Ispra (VA), Italy b CNR, Istituto di Matematica Applicata e Tecnologie Informatiche “Enrico Magenes”, Via E. Bassini 15, 20133 Milano, Italy Received 10 May 2013; received in revised form 20 March 2015; accepted 23 March 2015 Available online 13 April 2015

Abstract In this paper, we study a multi-trophic prey–predator system composed by top- and intermediate-level consumers and two food sources. Biomass fluxes across trophic levels are modeled by means of Holling-II type functional responses and one of the sources is subject to periodic impulsive events consisting in biomass injections. Referring to this system, we analyze the local and global stability properties of periodic eradicated solutions and establish permanence properties for the populations in terms of relations between the model parameters and the intensity of the impulses. Secondly, we formulate an optimal control problem for this system as well as an Iterative Dynamic Programming (IDP) scheme for its solution. By means of the proposed algorithm we specify the intensity of the impulsive controls to be applied periodically in order to regulate the biomass levels along time. In our formulation, the intensity of the periodic impulsive events depends on the biomass of the upper-level consumer population, as we can often see in an agricultural framework wherein the quantity of human workforce can significantly affect the size of crop periodically collected. & 2015 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.

1. Introduction In the recent literature, significant effort has been devoted to optimal natural resource exploitation strategies for human welfare and economic growth [1,2], which motivated a number of studies wherein n

Corresponding author. E-mail addresses: [email protected] (L. Galbusera), [email protected] (S. Pasquali).

http://dx.doi.org/10.1016/j.jfranklin.2015.03.031 0016-0032/& 2015 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

2729

the role of human presence is explicitly taken into account in modeling ecosystems, see for instance [3–5]. In this framework, control theory has been employed to face many natural resource management problems, such as harvesting [6,7], fishery management [8–10] and pest control [11–14]. Many authors dealt with the application of control theory to various classes of trophic chains, see for instance: [15] for optimal harvesting in single growing species; [16] for stabilization and synchronization problems in Lotka–Volterra models; [17] for optimal foraging in predator–prey systems with adaptive behaviors; [18,19], discussing the optimal control of a Rosenzweig– MacArthur tritrophic system in order to achieve sustainable management strategies with either top-down and bottom-up control layouts; [20], which accounts for various types of functional responses used to describe predator–prey interactions; [21], discussing stability and optimal control for different classes of Lotka–Volterra tritrophic systems. Furthermore, in recent studies on trophic systems impulses are often used to represent fast biomass injections at some levels of a trophic chain, or sudden events such as the birth of new individuals within a population or external perturbations. Models displaying this feature have been proposed, for instance, in the following papers: [22], dealing with a tritrophic Lotka– Volterra chain with impulses on the mid-level predator; [23], based on a tritrophic system with Holling-IV type functional responses and impulses acting on the top predator biomass; [24], studying a system with impulses on both prey and predator populations, local and global stability of periodic solutions and permanence conditions for the species; [25], dealing with a Lotka– Volterra system with competing species; see also references [26–30] for recent developments. From a broader point of view, the study of systems subject to impulses has been attracting much interest in the control community during the last years, and lead to a number of results concerning impulsive control [31,32], with applications in a number of fields including aerospace [33], multi-agent systems [34], finance [35], biology and medicine [36–39]. Optimal impulsive control, in particular, was successfully applied in order to define both the impulsive instants and intensities in a number of different contexts, including ecosystem management problems [40,41]. The main focus of this paper consists in the analysis and optimal impulsive control of a trophic system which differs from those most frequently considered in the previous literature in terms of both its structure and the features of the considered impulsive events. It comprises two food sources, an intermediate- and an upper-level consumer. While the first food source is accessed by the intermediate consumer, the second is used by the top-level consumer. A real-world example of such kind of system is represented by the collection of humans, herbivores and two plant sources representing the pasture used for herbivore feeding and the crop stocked for human use, respectively. In this framework, humans are assumed to feed on both the crop stock and the herbivores. Trophic interactions are represented by means of Holling-II models, while the growth of the pasture will be modeled as a continuous process along time. Here we opt for a Holling-II functional response, too, because we assume saturation for high levels of source density. In case saturation is not needed, a nonmonotonic Holling-IV functional response could be considered; in this case an inhibitory effect decreases the specific growth rate of the consumer, after it has reached a maximum (see [42] and related references). Furthermore, we assume that the crop stock undergoes a periodic filling-depletion cycle, determined by the agricultural cycle making its products available according to a specific seasonality. Considering this system, we will first discuss its structural features, including the boundedness of its trajectories, and the stability properties of specific periodic solutions of the system which influence the permanence of all the populations. Then, we will formulate an optimal control problem for this system, wherein the control action consists in the choice of the intensity of the periodic impulsive biomass injections. We will show how the presence of the crop stock, subject

2730

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

to impulses, can be exploited in order to regulate the biomass levels along time. We will also discuss an algorithm for the solution of such a control problem, based on an Iterative Dynamic Programming (IDP) approach, see [43]. IDP is a numeric method based on an iterated application of the Dynamic Programming approach, which involves: the choice of grid points for the control variables for each time instant inside the control horizon and the forward integration of the model; the application of Bellman's optimality principle proceeding backwards in time in order to choose the best control policy among the considered ones; the iteration of both these operations by refining the size of the search grid and focusing it around the best control policy found at each iteration. IDP has been applied for instance to chemical applications and optimization of biotechnological processes [44–47], vehicle control [48] and motion planning [49]. It has been applied to optimal control problems involving various classes of nonlinearities and high-dimensional systems while accounting for issues like the low sensitivity of the objective function with respect to the optimization variables and the splitting of the optimization horizon into stages of time-varying length. Therefore, it can be regarded as an interesting technique also in our application and in its possible extensions to different classes of trophic interaction models and aperiodic impulsive events. The paper is organized as follows: in Section 2 we introduce the trophic model under analysis and some preliminary results are summarized in Section 3; in Section 4 we discuss some fundamental structural properties of the system under analysis and introduce the associated periodic eradicated solutions assuming periodic impulses of constant intensity; their (local and global) stability properties are analyzed in Section 5, while permanence is the focus of Section 6; in Section 7 we introduce our optimal control problem formulation and the proposed IDP numerical algorithm for its solution; in Section 8 we propose a numerical example and finally there are some concluding remarks. Notation: N is the set of naturals and N0 ¼ N [ f0g; Rþ ¼ ½0; 1Þ and Rnþ ¼ ð0; 1Þ; n nn Rþ ¼ fx A Rn : x Z 0g and Rþ ¼ fxA Rn : x40g; PCðRþ ; RÞ (PC 1 ðRþ ; RÞ) is the class of real piecewise continuous (continuously differentiable) functions defined on Rþ ; M n ðRÞ is the set of all n  n real matrices; symbol (T) denotes the transpose, x_ ðtÞ the time derivative of x and modðÞ the modulo function. 2. Model formulation Define a countable set T  Rþ of time instants t k ; k A N0 such that t 0 ¼ 0, t kþ1 4t k ; 8k and limk-1 t k ¼ 1. Based on this definition, we can introduce the following continuous-time, impulsive dynamical model wherein the state variables represent the biomasses of mutually interacting food sources/consumers:
 8 x1 ðtÞ > > x_ 1 ðt Þ ¼ r 1 x1 ðt Þg  D13 f 31 ðx1 ðt ÞÞx3 ðt Þ > > K1 > > > > < x_ 2 ðtÞ ¼  q2 x2 ðtÞ  D24 f 42 ðx2 ðtÞÞx4 ðtÞ; 8 t A Rþ \T ð1Þ x2 ðt þ 8 k A N0 > k Þ ¼ x2 ðt k Þ þ hðt k Þ; >   > > > x_ 3 ðtÞ ¼ D31 f 31 ðx1 ðtÞÞ  q3 x3  D34 f 43 ðx3 ðtÞÞx4 ðtÞ > > > : x_ ðtÞ ¼ D f ðx ðtÞÞ þ D f ðx ðtÞÞ  q x ðtÞ 4 42 42 2 43 43 3 4 4 Defining x ¼ ½x1 ; x2 ; x3 ; x4 T , we denote initial condition of the system as xð0Þ ¼ x0 A R4þ . In our model, states x1 and x2 refer to two food sources, x3 to the intermediate-level consumer

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

2731

and x4 to the top-level consumer. The second food source represents a stock of biomass which is accessed exclusively by the top-level consumer and is subject to externally forced impulses hðt k ÞZ 0, representing biomass injections taking place at instants t k A T . In model (1), parameter r 1 40 is the specific growth rate of the first source, K 1 40 its carrying capacity and
 x1 x1 g ð2Þ ¼ 1 K1 K1   For all pairs ði; jÞA ð1; 3Þ; ð2; 4Þ; ð3; 4Þ , trophic interactions are modeled by means of a HollingII type functional response, i.e. f ji ðxi Þ ¼

xi aji þ xi

ð3Þ

where aji 40 is the half saturation constant, corresponding to the biomass of population i necessary to produce f ji ðxi Þ ¼ 12. The biomass flow across two different trophic levels depends on the functional response, the biomass of the upper trophic level (consumer) and the positive constants Dij representing the specific food demand of the consumer per time instant. For the same pairs as above, we also define Dji ¼ θji Dij where θji A ð0; 1 describes the efficiency of the biomass conversion associated to the trophic process. Quantities q2 ; q3 ; q4 40 are specific loss rates and, in particular, q2 describes the deterioration rate of the stocked biomass. In view of the properties of Eq. (3), we assume D31  q3 40 and D42 þ D43  q4 40: these are necessary conditions to avoid trivial dynamics in the ODEs associated to x3 and x4, respectively. A schematic representation of the system under analysis is provided in Fig. 1. Observe that, as detailed below in this section, in our optimal control formulation the magnitude of the impulses will depend on x4. From now on in this paper we assume the periodicity of the impulsive events, i.e. t kþ1  t k ¼ P;

8 k A N0

ð4Þ

where P40 is a constant. Periodic impulsive events correspond to many real-world situations, as the one in which crop becomes periodically available due to the natural cycle and is then stocked. Furthermore, for analysis purposes, in most of Sections 4 and 5 we will also assume hðt k Þ ¼ h;

8 k A N0

ð5Þ

Here h40 is a constant depending on the production efficiency. Instead, in the last part of this paper we will allow for time-varying impulse intensities, for control purposes. In particular, we

Fig. 1. System (1): trophic interactions and impulsive component.

2732

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

will assume that hðt k Þ is a bounded control variable 0 r hðt k Þr h

x4 ðt k Þ ; d þ x4 ðt k Þ

8 k A N0

ð6Þ

where d40 is the half saturation constant. In the real-world situation exemplified above, such an upper bound can be used to represent the dependency of the crop being produced on the availability of workforce. Observe that, for simplicity, in Eq. (6) the impulse acting on x2 at time tk just depends on the value of x4 at the same time instant. 3. Preliminary results This section is organized in two parts: Section 3.1 reports literature results related to comparison systems to be used for studying the solutions of ODEs, while Section 3.2 summarizes some mathematical tools relevant for studying the local stability properties of periodic solutions of impulsive systems. 3.1. Comparison system The following lemma provides a bound on the solutions of a system of impulsive ODEs [50,30]: Lemma 1. Consider a function u A PC 1 ðRþ ; RÞ satisfying the inequalities 8 _ ð Z ÞpðtÞuðtÞ þ f ðtÞ; 8 t A Rþ \T > < uðtÞr þ uðt k Þ r ðZ Þdk uðt k Þ þ hk ; 8 k A N0 > : uð0þ Þ r ðZ Þu 0 where p; f APCðRþ ; RÞ and dk Z 0, hk and u0 are constants. Then, for t40, ! R ! R Z t t t pðτÞ dτ pðτÞ dτ uðtÞr ð Z Þu0 ∏ dk e 0 þ ∏ dk e s f ðsÞ ds 0ot k ot

þ

X 0ot k ot

! R t pðτÞ dτ ∏ d j e tk hk

0

s r t k ot

t k ot j ot

3.2. Stability result for impulsive linear differential systems Here we briefly present a stability result for impulsive systems and periodic systems of ODEs that will be employed to perform local stability analysis of particular solutions of system (1). Now introduce the following impulsive linear differential system: ( X_ ðtÞ ¼ AðtÞXðtÞ; 8 t A Rþ \T ð7Þ Xðt þ 8k A N0 k Þ ¼ Xðt k Þ þ Bk XðtÞ; where (a) AðÞ A PCðR; M n ðRÞÞ and there is T40 such that Aðt þ TÞ ¼ AðtÞ, 8 t Z 0; (b) Bk A M n ðRÞ, detðI n þ Bk Þa 0 for k A N0 ; (c) there exists qA N such that Bkþq ¼ Bk , t kþq ¼ t k þ T for k A N0 .

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

2733

Let Ψ ðtÞ be a fundamental matrix associated to system (7). Then there is a unique nonsingular monodromy matrix Φ A M n ðRÞ such that Ψ ðt þ TÞ ¼ Ψ ðtÞΦ; 8 t A R. All monodromy matrices of Eq. (7) are independent of t and share the same eigenvalues λ1 ; …; λn , which are called Floquet multipliers of Eq. (7). These determine the stability properties according to the following lemma [50]: Lemma 2. If assumptions (a)–(c) hold, system (7) is

  

stable iff all Floquet multipliers λi ; i ¼ 1; …; n, satisfy jλi j r 1 and, if jλi j ¼ 1, then to λi there corresponds a simple elementary divisor; asymptotically stable iff all Floquet multipliers λi ; i ¼ 1; …; n, satisfy jλi jo1; unstable if there is a Floquet multiplier λi such that jλi j41.

4. Invariant sets, solution boundedness and periodic eradicated solutions system In view of Lemma 1, we can prove the following statement about system (1): n4 Lemma 3. Rþ is invariant for system (1). n4 . We can observe that Proof. Consider a saturated solution xðtÞ : ½0; t 0 Þ↦R4 with xð0ÞA Rþ 8 x_ 1 ðtÞZ p1 ðtÞx1 ðtÞ > > > > < x_ 2 ðtÞZ p2 ðtÞx2 ðtÞ x_ 3 ðtÞZ p3 ðtÞx3 ðtÞ > > > > : x_ 4 ðtÞ ¼ p ðtÞx4 ðtÞ 4  1 as long as the solutions remain positive and t 2 = T , where p1 ðt Þ ¼ r 1 g xK1 ðtÞ x3 ðt Þ,  D13 a31 1 1 1 p2 ðtÞ ¼  q2  D24 a42 x4 ðtÞ, p3 ðtÞ ¼ D31 f 31 ðx1 ðtÞÞ  q3  D34 a43 x4 ðtÞ and p4 ðtÞ ¼ D42 f 42 ðx2 ðtÞÞ þ D43 f 43 ðx3 ðtÞÞ  q4 . Consequently, taking into account the non-negativity of hðt k Þ; 8 k A N0 , by Lemma 1 we have 8 Rt p ðτÞ dτ > > 0 1 x ðtÞZ x ð0Þe > 1 1 > R > t > > < x2 ðtÞZ x2 ð0Þe 0 p2 ðτÞ dτ Rt p ðτÞ dτ > > 0 3 x ðtÞZ x ð0Þe > 3 3 > R > t > > : x ðtÞZ x ð0Þe 0 p4 ðτÞ dτ

4

4

n4 n4 Therefore, as xð0ÞA Rþ , by the latter inequalities it results xðtÞA Rþ , 8 t A ½0; t 0 Þ.

□

Furthermore, defining qm ¼ min qi i A f2;3;4g

we can prove the following result on trajectory boundedness: n4 Lemma 4. All solutions xðÞ of Eq. (1) subject to Eq. (5) and with x0 A Rþ are bounded.

ð8Þ

2734

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

n4 Proof. Consider a solution x(t) of system (1) with x0 A Rþ and define

VðxÞ ¼ θ43 θ31 x1 þ θ42 x2 þ θ43 x3 þ x4

ð9Þ

For t A Rþ \T it results


 x1 ðtÞ V_ ðxðt ÞÞ ¼ θ43 θ31 r 1 x1 ðt Þg  θ42 q2 x2 ðt Þ θ43 q3 x3 ðt Þ q4 x4 ðt Þ K1

It is easy to verify that, 8 t A Rþ \T ,


  x1 ðtÞ V_ ðxðt ÞÞ þ qm V ðxðt ÞÞr θ43 θ31 r 1 g þ qm x1 ðt Þ K1 Consequently, the following upper bound holds: V_ ðxðtÞÞ þ qm VðxðtÞÞ r C with

ð10Þ





  x1 ðr 1 þ qm Þ2 K 1 þ qm x1 ðt Þ ¼ θ43 θ31 C ¼ max θ43 θ31 r 1 g x1 K1 4r 1

ð11Þ

On the other side, 8 k A N0 , Vðxðt þ k ÞÞ ¼ Vðxðt k ÞÞ þ θ 42 h

ð12Þ

Finally, formulas (10) and (12) imply X

1  e  qm t V ðxðt ÞÞr V x 0þ e  qm t þ C þ e  qm ðt  tk Þ θ42 h qm 0ot ot

ð13Þ

k

In view of the latter inequality, it results that V is bounded and x is also bounded as a consequence. □ Now define C θ42 h G h ¼ þ qm 1 e  qm P

ð14Þ

The next corollary about the ultimate boundness of the solutions of Eq. (1) follows immediately from Eqs. (9) and (13): Corollary 1. The solutions of Eq. (1) subject to Eqs. (4) and (5) fulfil lim supt-1 1 1 x1 ðtÞr ðθ43 θ31 Þ  1 GðhÞ, lim supt-1 x2 ðtÞr θ42 GðhÞ, lim supt-1 x3 ðtÞr θ43 GðhÞ and lim supt-1 x4 ðtÞr GðhÞ. To complete this section, we address the existence of periodic solutions for the eradicated forms of system (1) subject to Eqs. (4) and (5), here defined as the ones in which at least one among populations 1, 3 and 4 is extinct. Theorem 1. System (1) subject to Eqs. (4) and (5) admits the following periodic eradicated solutions: 0

S : x0 ðtÞ ¼ ð0; x 2 ðtÞ; 0; 0Þ

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

2735

1

S : x1 ðtÞ ¼ ðK 1 ; x 2 ðtÞ; 0; 0Þ and, assuming D31 f 31 ðK 1 Þ q3 40 (equivalently 1  ðD31q3 aq31 ÞK 1 40), 3

S

13

13 : x13 ðtÞ ¼ ðx 13 1 ; x 2 ðtÞ; x 3 ; 0Þ

where x 2 ðt Þ ¼ x 02 eð  q2 Þt mod P ;

x 02 ¼

h 1  e  q2 P

ð15Þ

q3 a31 q3 a31 θ31 r1 a31 13 D31  q3 , x 3 ¼ D31  q3 ð1  ðD31  q3 ÞK 1 Þ. Furthermore, x2 ðtÞZ x 2 ðtÞ if x02 Z x 02 and x2 ðtÞox 2 ðtÞ if x02 ox 02 .

and x 13 1 ¼

0

1

limt-1 jx2 ðtÞ x 2 ðtÞj ¼ 0. Finally,

13

Proof. Referring to either S , S or S , the time evolution of the state variable x2 in Eq. (1) is governed by equations ( x_ 2 ðtÞ ¼  q2 x2 ðtÞ; 8 t A Rþ \T ð16Þ x2 ðt þ 8 k A N0 k Þ ¼ x2 ðt k Þ þ h; Thus, 8 t; t 0 A ðt k ; t kþ1  and 8k A N0 it results x2 ðtÞ ¼ e  q2 ðt  t 0 Þ x2 ðt 0 Þ and, by assumption (4), þ þ  q2 P x2 ðt þ x2 ðt þ kþ1 Þ ¼ x2 ðt kþ1 Þ þ h ¼ e k Þ þ h. Imposing x2 ðt kþ1 Þ ¼ x2 ðt k Þ; 8 k A N0 , we get the þ periodicity condition x2 0 ¼ 1  eh q2 P , from which Eq. (15) follows naturally. The solutions for x1 and x3 correspond to the different equilibria obtained with x4 ¼ 0. As for the second part of the theorem's statement, observe that the solution of Eq. (16) is given by h  q2 ðt  t k Þ x2 ðt Þ ¼ ðx2 t þ þ x 2 ðt Þ for t A ðt k ; t kþ1 , 8 k A N0 . Hence k  1  e  q2 P Þe 0 0 limt-1 jx2 ðtÞ x 2 ðtÞj ¼ 0 and x2 ðtÞZ x 2 ðtÞ if x2 ð0ÞZ x 2 , while x2 ðtÞox 2 ðtÞ if x2 ð0Þox 2 . □ Observe that system (1) with x2 ¼ x4 ¼ 0 is the two-species Rosenzweig–MacArthur (RMA) model which can admit limit cycles, see for instance [51]; the detailed discussion of this topic is beyond the scope of our analysis. 5. Stability of the periodic eradicated solutions This section is devoted to the discussion of the local and global stability properties of the periodic eradicated solutions of system (1) reported in Theorem 1. 5.1. Local stability In order to study the local stability properties of the periodic solutions reported in Theorem 1, we introduce the following variables: ξ1 ðtÞ ¼ x1 ðtÞ x 1 ;

ξ2 ðtÞ ¼ x2 ðtÞ x 2 ðtÞ;

ξ3 ðtÞ ¼ x3 ðtÞ x 3 ;

ξ4 ðtÞ ¼ x4 ðtÞ

Here the periodic function x 2 ðtÞ is defined as in Eq. (15) and the constants x 1 ; x 3 may be associated with any of the periodic eradicated solutions discussed above. Substituting the latter relationships in system (1) and linearizing it around ðx 1 ; x 2 ðtÞ; x 3 ; 0Þ leads to the following

2736

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

system: _ ¼ MðtÞξðtÞ ξðtÞ where

3 ξ1 6 7 6 ξ2 7 7 ξ¼6 6 ξ3 7; 4 5 ξ4 2

ð17Þ 2

m11

6 0 6 MðtÞ ¼ 6 4 m31 0

0

m13

m22

0

0 m42

m33 m43

i with mij ¼  Dij ajixþx for ioj, mij ¼ Dij i

0

3

m24 7 7 7 m34 5 m44 aij

 1 x i for i4j, m11 ¼ r 1 1  2 Kx 11  θ31 m31 ,

ðaij þx j Þ 1 1 m42 , m33 ¼  θ31 m13  q3  θ43 m43 and m44 ¼  θ42 m24  θ43 m34  q4 . m22 ¼  q2  θ42 Observe that m24, m42 and m44 vary with time, since they depend on x 2 ðtÞ. Also notice that model (17) is not impulsive, thanks to the properties of the P-periodic function x 2 ðtÞ. Defining ν1 ¼

q4 D42

and ν13 ¼

q4 D42

D43 x 13 3 13 42 ða43 þx 3 Þ

D

2

, we can introduce the following theorem: 0

1

Theorem 2. The periodic eradicated solutions S , S and S and (5) fulfil the following local stability properties:

13

of system (1) subject to Eqs. (4)

0

 S is unstable; 1  S is LAS (locally asymptotically stable) iff D31 f 31 ðK 1 Þ q3 o0 ð18Þ νqP



e 1 2  1 1  e  q2 P hoa42 if ν1 o1 ð19Þ 1  eðν1  1Þq2 P 13  S , which exists when D31 f 31 ðK 1 Þ q3 40, is LAS iff ( m11 þ m33 o0 ð20Þ m11 m33  m13 m31 40 and ν qP



e 13 2  1 1  e  q2 P hoa42 if ν13 o1 ð21Þ 1  eðν13  1Þq2 P and

Proof. In view of Lemma 2, the local stability properties of the periodic eradicated solutions reported in the theorem statement can be assessed by means of Floquet multipliers, consisting in the eigenvalues of the monodromy matrix ΦðPÞ such that ξðt þ PÞ ¼ ΦðPÞξðtÞ in system (17). We address the different cases separately: 0

 S : M(t) is upper triangular and consequently the corresponding monodromy matrix Φ0 ðtÞ has the same structure, with   Rt D42 f 42 ðx 2 ðτÞÞ  q4 dτ 0 r1 t  q2 t  q3 t diagðΦ ðtÞÞ ¼ e ; e ; e ; e 0 Denote by λ0;i the eigenvalues of Φ0 ðPÞ, for i ¼ 1; …; 4. It is easy to observe that λ0;1 41, which implies the instability of the periodic solution of S 0 .

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

2737

1

 S : M(t) is upper triangular and the corresponding monodromy matrix Φ1 ðtÞ has   Rt D42 f 42 ðx 2 ðτÞÞ  q4 dτ 1  r1 t  q2 t ðD31 f 31 ðK 1 Þ  q3 Þt diagðΦ ðtÞÞ ¼ e ; e ; e ; e 0 We proceed similar to the previous case, evaluating the associated eigenvalues λ1;i ; i ¼ 1; …; 4 of Φ1 ðPÞ. In particular, it results λ1;1 ; λ1;2 o1. Furthermore, λ1;3 o1 iff condition (18) holds. As for term λ1;4 , we can proceed as follows: Z P RP D42 f 42 ðx 2 ðτÞÞ  q4 dτ 0 λ1;4 ¼ e o1 3 D42 f 42 ðx 2 ðτÞÞ  q4 dτo0 0



3x 02 1  eðν1  1Þq2 P oa42 eν1 q2 P  1 from which Eq. (19) follows by studying the feasibility of the last inequality. 13  S : in this case, M(t) has m22 ¼  q2 ; m33 ¼  θ31 m13  q3 and m42 ; m43 ¼ 0. Thus, the characteristic polynomial associated to the monodromy matrix Φ13 ðPÞ is

PðΦ13 ðPÞÞ ¼ ð  λ þ m22 Þð  λ þ m44 Þ λ2  ðm11 þ m33 Þλ þ ðm11 m33  m13 m31 Þ Observe that λ13;1 ¼ m22 ¼  q2 o0, while term λ2  ðm11 þ m33 Þλ þ ðm11 m33  m13 m31 Þ has negative roots iff the Routh–Hurwitz conditions (20) hold. Finally, since term  λ þ m44 depends on the periodic signal x 2 ðtÞ, it is necessary to evaluate the corresponding Rt D42 f 42 ðx 2 ðτÞÞþD43 f 43 ðx 13 3 Þ  q4 dτ 0 . Provided term in the monodromy matrix Φ13 ðtÞ, i.e., ϕ13 44 ¼ e that the conditions (20) hold, the following additional one is required to ensure that the considered periodic solution is LAS: RP D f ðx ðτÞÞþD43 f 43 ðx 3 Þ  q4 dτ λ13;4 ¼ e 0 42 42 2 o1 Z P D42 f 42 ðx 2 ðτÞÞ þ D43 f 43 ðx 13 3 3 Þ q4 dτo0 0



3x 02 1  eðν13  1Þq2 P oa42 eν13 q2 P  1 Similar to the previous case, Eq. (21) follows by studying the feasibility of the last inequality. □

It is now possible to provide an interpretation of the LAS conditions provided in Theorem 2. Condition (18) describes the inability of source 1, even at its carrying capacity, to allow the survival of population 3. Instead Eq. (19) describes the inability of source 2 to sustain population 4 at the considered impulsive levels and periodicity. Condition (20) is the stability condition of the two13 species RMA model, while Eq. (21) is the analogue of Eq. (19) referred to the case of S . 5.2. Global stability This section is devoted to the analysis of the global stability properties of the periodic eradicated solutions specified in Theorem 1. The next theorem clearly shows the relationships between these properties and the local stability of the considered solutions, studied in Section 5.

2738

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749 1

Theorem 3. The periodic eradicated solutions S and S and (5) fulfil the following stability properties:

13

of system (1) subject to Eqs. (4)

1

n3  in the subspace ðx1 ; x2 ; x4 Þ ARþ (with x3 ¼ 0), S is GAS (globally asymptotically stable) iff condition (19) holds. 13 K1 n4  S is GAS in ðx1 ; x2 ; x3 ; x4 ÞA Rþ if condition (21) holds and there exist xGAS 1 4 2 and t Z 0 is sufficiently large so that x1 ðtÞ Z xGAS 8tZt ð22Þ 1 ;

Proof. Also in this case, we discuss the different cases separately. 1

 S : (Sufficiency) First observe that, since  we are assuming x3 ¼ 0, K1 is a GAS equilibrium point of the ODE x_ 1 ¼ r 1 x1 1  Kx11 on Rnþ . Now introduce the following system in order to study the dynamics of the second population: ( x~_ 2 ðtÞ ¼  q2 x~ 2 ðtÞ; 8 t A Rþ \T ~ 2 ðt k Þ þ h; x~ 2 ðt þ k Þ¼x

8 k A N0

Assuming x2 ð0Þ ¼ x~ 2 ð0Þ A Rnþ , by a standard comparison argument it results x2 ðtÞr x~ 2 ðtÞ. Furthermore, by the same rationale as in the proof of Lemma 1, we can prove that limt-1 j~x 2 ðtÞ x 2 ðtÞj ¼ 0, where x 2 ðtÞ is defined as in Eq. (15). Therefore, ( t 2 40 such that x2 ðtÞr x 2 ðtÞ þ ϵ; 8t Z t 2 . As for state variable x4, the assumption x3 ¼ 0 and the relationships studied above imply     x_ 4 ðtÞr D42 f 42 ðx2 ðtÞÞ  q4 x4 ðtÞr D42 f 42 ðx 2 ðtÞ þ ϵÞ q4 x4 ðtÞ where the latter inequality holds for t sufficiently large and ϵ40 arbitrarily small, as discussed above. Consequently, 8 n A N0 we have R ðnþ1ÞP D42 f 42 ðx 2 ðτÞþϵÞ  q4 dτ x4 ððn þ 1ÞPÞr x4 ðnPÞe nP Now R P observe that, since ϵ40 is chosen arbitrarily small, Eq. (19) is equivalent to 0 D42 f 42 ðx 2 ðτÞ þ ϵÞ q4 dτo0. Consequently, limn-1 x4 ðnPÞ ¼ 0. Furthermore, 0ox4 ðtÞox4 ðnPÞeD42 P ;

8 t A ½nP; ðn þ 1ÞP;

8 nA N0

thus also limt-1 x4 ðtÞ ¼ 0. As a consequence of the latter property, we also have limt-1 x2 ðtÞ ¼ limt-1 x~ 2 ðtÞ, thus limt-1 jx2 ðtÞ x 2 ðtÞj ¼ 0. 1 (Necessity) When Eq. (19) does not hold, some components of S are unstable in view of Theorem 2. 13 1  S : considering the state variable x2 in system (1), we can proceed as in the case of S , concluding that (t 2 40 such that x2 ðtÞr x 2 ðtÞ þ ϵ; 8t Z t 2 for ϵ40 arbitrarily small and limt-1 j~x 2 ðtÞ x 2 ðtÞj ¼ 0, where x 2 ðtÞ is defined as in Eq. (15). Furthermore, consider the subsystem of Eq. (1) composed by populations 1 and 3 only, assuming 13 x4 ¼ 0. The global stability properties of the associated coexistence equilibrium ðx 13 1 ; x3 Þ can be assessed by means of a Lyapunov approach, through a procedure similar to the one presented in [52].  To this end, let usintroduce the following Lyapunov function on the domain G13 ¼ ðx1 ; x3 Þ : x1 ; x3 A Rnþ :

 D31 x1 D13 x3 13 13 13 13 V 13 ¼ x  x  x ln x  x  x ln þ 1 3 1 1 3 3 a31 a31 þ x 13 x 13 x 13 1 1 3

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

2739

The Lyapunov function's derivative along the system's positive solutions fulfils 13 V_ ¼



 D31 x 13 D13 x 13 1 3 1 1 x_ 1 ðt Þ þ x_ 3 ðt Þ x1 a31 x3 a31 þ x 13 1

Furthermore, replacing x_ 1 ðtÞ and x_ 3 ðtÞ by model (1) and exploiting the equilibrium  D x 13 x 13 D x 13 relationships 0 ¼ r 1 1 K11  a 13þx313 and 0 ¼ a 31þx113  q3 , it is possible to demonstrate 31

1

31

1

that the following property holds:


 2

13 D31 ðx1  x 13 x1 þ x 13 D13 D31 x 13 x 13 1 Þ 1 1 x3 1 r1 1  1 þ 2 x1 K1 x1 a31 þ x 13 ða31 þ x 13 1 1 Þ
 13 13 D13 D31 x 13 x D D x x D x 13 31 1 3 13 1 3 þ  q3 x3 1 3  a ða þ x Þ a x3 a31 ða31 þ x 13 Þ 31 31 1 31 1  2

13 13 13 D31 ðx1  x 13 Þ x þ x D D x 1 13 31 1 x 3 1 1 ¼ r 1 1 þ 13 x1 K1 a31 þ x 1 ða31 þ x 13 1 Þa31

 13 13 x 1 ða31 þ x1 Þ x1 ða31 þ x 1 Þ 2  13 x 1 ða31 þ x1 Þ x1 ða31 þ x 13 1 Þ

13 V_ ¼

Since the arithmetic mean is greater or equal to the geometric mean, it results 2

x 13 x1 ða31 þ x 13 1 ða31 þ x1 Þ 1 Þ  r0 13 13 x1 ða31 þ x 1 Þ x 1 ða31 þ x1 Þ

13 where the equality only holds for x1 ¼ x 13 1 (and x3 ¼ x 3 , consequently). Furthermore, by our assumptions x1 ðtÞZ xGAS for t Z t with t sufficiently large, which obviously also 1 13 13 GAS Z x . Therefore, condition (22) implies V_ r 0 and in particular V_ ¼ 0 implies x 13 1 1

13 13 13 only for ðx1 ; x3 Þ ¼ ðx 13 1 ; x 3 Þ. Thus, by the LaSalle principle, we conclude that ðx 1 ; x 3 Þ is GAS for the considered subsystem under the specified assumptions. As it easy to verify, this property implies that also inequalities (20) hold. Furthermore, if we perturb the dynamics of this subsystem by means of x4 40, by comparison arguments it is possible to conclude that, still under the assumption (22), 13 lim inf t-1 x1 ðtÞZ x 13 1 and lim supt-1 x3 ðtÞ r x 3 , and consequently ( t 1 ; t 3 40 such that 13 13 x1 ðtÞZ x 1  ϵ; 8t Z t 1 and x3 ðtÞ rx 3 þ ϵ; 8 t Z t 3 and ϵ40 arbitrarily small. A consequence of the latter argument is that    x_ 4 ðtÞr D42 f 42 ðx2 ðtÞÞ þ D43 f 43 ðx3 ðtÞÞ  q4 x4 ðtÞr D42 f 42 ðx 2 ðtÞ þ ϵÞ  þD43 f 43 ðx 13 3 þ ϵÞ q4 x4 ðtÞ

where the latter inequality holds for t sufficiently large and ϵ40 arbitrarily small. Therefore, 8 nA N0 we have R ðnþ1ÞP D42 f 42 ðx 2 ðτÞþϵÞþD43 f 43 ðx 13 3 þϵÞ  q4 dτ x4 ððn þ 1ÞPÞr x4 ðnPÞe nP Now observe that, since ϵ40 can be chosen arbitrarily small, Eq. (21) is equivalent to RP 13 0 D42 f 42 ðx 2 ðτÞ þ ϵÞ þ D43 f 43 ðx 3 þ ϵÞ  q4 dτo0. Then, we conclude limn-1 x4 ðnPÞ ¼ 0. Furthermore, 0ox4 ðtÞox4 ðnPÞeðD42 þD43 ÞP ; 8 t A ½nP; ðn þ 1ÞP, thus also limt-1

2740

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

x4 ðtÞ ¼ 0. The latter limit implies that limt-1 x1 ðtÞ ¼ x~ 1 ðtÞ ¼ x 13 limt-1 1 , 13 x2 ðtÞ ¼ x 2 ðtÞ ¼ 0 and limt-1 x3 ðtÞ ¼ x~ 3 ðtÞ ¼ x 3 . Therefore the proof is completed. □ Remark 1. It is possible to provide a sufficient condition ensuring the fulfilment of the GAS 13

condition (22) referred to the case of S . To this end, recall that Corollary 1 proves 1 GðhÞ, with GðhÞ defined as in Eq. (14). Consequently, in view of lim supt-1 x3 ðtÞ r θ43 h  i D13 Lemma 1 we have lim inf t-1 x_ 1 ðt ÞZ x1 ðt Þ r 1 1  xK1 ðtÞ G h , which yields  a31 θ43 1  ¼ K 1 1 Da3113θGðhÞ Z 0. In particular, imposing , provided that 1 Da3113θGðhÞ lim inf t-1 x1 ðt ÞZ xGAS 1 43 r 1 43 r 1 1 Da3113θGðhÞ 4 12, equivalently 43 r 1

 ðr 1 þ qm Þ2 K 1 h r 1 a31 1 D13 θ43 θ43 θ31 þ o  q P m 1 e 4r 1 qm 2 condition (22) is verified. 6. Permanence In this section, we study the coexistence of all the populations composing system (1), that is a situation in which no population becomes extinct or grows indefinitely. In particular, we refer to the following notion of permanence for system (1) (see also [24] and related references): n4 , there exist Definition 1. System (1) is said to be permanent if, for any solution with xð0Þ A Rþ n4 min max min max constants t40 and x ; x A Rþ such that x r xðtÞr x , 8 t Z t.

Observe that permanence of system (1) excludes the asymptotic stability of the periodic eradicated solutions reported above. Next, we introduce a sufficient condition for permanence in system (1). Theorem 4. System (1) is permanent if GðhÞ 40 a31


    r1 x1 M 1 ða31 þ x1 Þ x3 x1  xm 1 r 0; D13 K1 Z P D42 f 42 ðx 2 ðτÞÞ þ D43 f 43 ðxm 3 Þ q4 dτ40 D31 f 31 ðK 1 Þ  q3  D34

ð23aÞ 8 x1 A Rnþ

ð23bÞ ð23cÞ

0

M where, denoting by xm 1 and x1 the unique values of x1 defining the isoclines D31 f 31 ðx1 Þ  q3 ¼ 0 GðhÞ 13 13 M M m and D31 f 31 ðx1 Þ  q3  D34 a43 ¼ 0, ðxm 1 ; x3 Þ ¼ ðx 1 ; x 3 Þ and ðx1 ; x3 Þ are their intersections with x3 x1 the isocline r 1 g K 1  D13 a31 þx1 ¼ 0.

Proof. The existence of a finite xmax fulfilling Definition 1 is proven by Lemma 4. Consequently, to n4 complete the proof of the theorem we need to demonstrate the existence of xmin A Rþ fulfilling the same definition. To this end, first observe that Eq. (23a) excludes the stability of the equilibrium with n2 ðx1 ; x3 Þ ¼ ðK 1 ; 0Þ. Notice that lim supt-1 x1 ðtÞr K 1 and define ðxn1 ðx4 Þ; xn3 ðx4 ÞÞ A Rþ as the x3 x1 x4 intersection between the isoclines r 1 gðK 1 Þ D13 a31 þx1 ¼ 0 and D31 f 31 ðx1 Þ q3  D34 a43 þx3 ¼ 0 for a

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

2741

specified value of x4. Then, introduce the Lyapunov function x4
 Z x1 D31 f 31 ðx1 Þ q3  D34 x3 a43 þ x3 13 dx1 þ x3  xn3 ðx4 Þ xn3 ðx4 Þ ln n V ¼ D13 f 31 ðx1 Þ x3 ðx4 Þ xn1 ðx4 Þ   on the domain G13 ¼ ðx1 ; x3 Þ : x1 ; x3 A Rnþ . The Lyapunov function's derivative along the system's trajectories is  0 1
 x r 1  x1 1 1 K1 x 13 4 @  xn3 ðx4 ÞA ð24Þ V_ ¼ D31 f 31 ðx1 Þ q3  D34 D13 f 31 ðx1 Þ a43 þ x3 It results that, for x4 ¼ 0, Eq. (23b) ensures V_ 13 r 0 and in particular V_ 13 ¼ 0 only in fðx1 ; x3 Þ : x1 ¼ xn1 ð0Þ; x3 40g. As ðxn1 ð0Þ; xn3 ð0ÞÞ is the largest invariant set in it, from the LaSalle n2 principle ðxn1 ð0Þ; xn3 ð0ÞÞ is GAS in Rþ . Now observe that, by Corollary 1 and Eq. (23a), it follows that n2 there exists an isocline intersection ðxn1 ðx4 Þ; xn3 ðx4 ÞÞ A Rþ for any admissible asymptotic value of x4. 3 Furthermore, as Eq. (23b) requires that the isocline r 1 gðKx11 Þ D13 a31xþx ¼ 0 is decreasing in 1 n n n ðx1 ð0Þ; x3 ð0ÞÞ and the same property holds for all x1 Z x1 ð0Þ, it also ensures the GAS property of ðxn1 ðx4 Þ; xn3 ðx4 ÞÞ for any constant x4 asymptotically admissible. Therefore, we can construct a set m n n Gn13 ¼ fðx1 ; x3 Þ : x1 Z xm 1 ; x3 Z x3 g, which encloses all points ðx1 ðx4 Þ; x3 ðx4 ÞÞ and results to be attractive. We conclude that for an arbitrarily small ϵ40 and t sufficiently large, it results m min m x1 ðtÞ4xmin 1 ¼ x1  ϵ40 and x3 ðtÞ4x3 ¼ x3  ϵ40. The proof proceeds by exploiting again Corollary 1 and the properties of Eq. (3) to observe that lim inf x_ 2 ðtÞZ η2 x2 ; t A Rþ \T t-1 η P with η2 ¼  q2  Da4224 G h . Consequently, we have x2 ðtÞ4~x 2 ðtÞ ϵ for ϵA ð0; 1 e 2eη2 P hÞ arbitrarily eðη2 Þt mod P small and t large enough, where x~ 2 ðt Þ ¼ 1  eη2 P h, by the same argument as in Theorem 1. eη2 P Therefore we can fix xmin 2 ¼ 1  eη2 P h  ϵ40. In order to complete the proof, it remains to demonstrate the existence of a value xmin fulfilling Definition 1. To this end, let us first fix a 4 positive value x min such that 4 R ðnþ1ÞP ϱ¼e with x^ 20 ¼

nP

D42

τ modP ð  q2  ðD24 =a42 Þx min 4 Þ  q2  ðD24 =a42 Þx min Þτ modP ð 4 a42 þ^x 20 e x^ 20 e

þD43 a

xm 3 m 31 þx3

 q4 dτ

41

ð25Þ

h

. Observe that this is always possible in view of assumption (23c). Þ and x4 ðtÞox min for t4t 1 and consider system Now assume that ( t 1 40 such that x3 ðtÞZ xmin 3 4
 8 D24 min > > x x^ 2 ðt Þ; 8 t A Rþ \T x^_ 2 ðt Þ ¼  q2  > > a42 4 > < 1  eð

 q2  D24 =a42 x min P 4

x^ 2 ðt þ 8 k A N0 k Þ ¼ x2 ðt k Þ þ h; >
 > > ^ xmin x > 2 3 > : x^_ 4 ðt Þ ¼ D42 ^ 4 ðt Þ þ D43  q 4 x a42 þ x^ 2 a43 þ xmin 3

ð26Þ

By the same rationale used in Theorem 1, we conclude that this system admits the following periodic solution for x^ 2 : min x^ 2 ðtÞ ¼ x^ 20 eð  q2  D24 =a42 x 4 Þt mod P

2742

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

Furthermore, limt-1 j^x 2 ðtÞ  x^ 2 ðtÞj ¼ 0. As for x^ 4 ðtÞ, we have R ðnþ1ÞP min D42 x^ 2 ðτÞ=ða42 þ^x 2 ðτÞÞþD43 xmin 3 =ða43 þx3 Þ  q4 dτ x^ 4 ððn þ 1ÞPÞ ¼ x^ 4 ðnPÞe nP so that x^ 4 ððn þ kÞPÞ-ϱk x^ 4 ðnPÞ-1 for n-1 and k finite, in view of the properties of x^ 2 ðtÞ and x3 ðtÞ stated above and Eq. (25). Comparing Eq. (26) with system (1), we have x2 ðtÞZ x^ 2 ðtÞ and x4 ðtÞZ x^ 4 ðtÞ, 8 t Z t 1 , assuming equal initial conditions at t 1 . Thus, we conclude limn-1 x4 ðnPÞ ¼ 1, which is a contradiction with the ultimate boundedness of trajectories claimed by Corollary 1. Consequently, we exclude that (t 1 40 such that x4 ðtÞox min 4 for t4t 1 . In order to complete the proof, also consider the possibility that, choosing any t 1 40 such that min min x4 ðt 1 Þ ¼ x min 4 , ( t 2 4t 1 such that x4 ðt 2 Þ ¼ x 4 and, for any t 1 otot 2 , x4 ðtÞox 4 . First notice that, min if x4 ðtÞox 4 for any t4t 2 , by the same argument as in the last paragraph, we obtain a contradiction with the boundedness property of the solution x4 ðtÞ. Therefore, fix t 3 ¼ infft : x4 ðtÞZ x min 4 ; t4t 2 g, so that t 3 Z t 2 . If we iterate this interval construction procedure forward in time and the process is forced to stop in finite time, the proof is completed by the same rationale as above. Differently, define an interval sequence ½t 2k  1 ; t 2k ; k A N with min t 2k  1 ot 2k r t 2kþ1 such that x4 ðtÞox min 4 ; t A ðt 2k  1 ; t 2k  and x4 ðt n Þ ¼ x 4 ; nA N. Then, it is possible to show that ( T 0 40 finite such that supft 2k  t 2k  1 ; k A Ng ¼ T 0 . In fact, assume by contradiction that there is a subsequence t 2ki  t 2ki  1 -1, k i -1. Then, by an argument similar to the one used in studying x^ 4 ðnPÞ above, it follows that x4 ðtÞ grows boundlessly, which leads to _ 4 Z x4 A, where a contradiction since x4 ðt 2k Þ ¼ x min 4 . Finally, observe that lim inf t-1 x xmin

xmin

min AT 0 min min 3 2 ^ min ^ 4 g it A ¼ D42 a42 þx . Putting xmin min þ D43 a þxmin  q4 . Then, fix x 4 ¼ minfx 4 ; x 4 ¼ x4 e 43 2

3

for t sufficiently large. Thus, the proof is completed. follows x4 ðtÞZ xmin 4

□

Observe that, in the statement of Theorem 4, conditions (23a) and (23b) ensure the existence of a xmin 1 40 for populations 1 and 3. They represent an extension of the GAS conditions expressed in Theorem 3.2 of [53]. Observe, in particular, that Eq. (23a) implies the instability of 1 the eradicated solution S and involves an adequate compensation between the biomass incoming from the first source and the intensity of the impulses, which support the growth of the top-level consumer also accessing biomass 3. Finally, condition (23c) implies the instability of the eradicated equilibrium which would exclude the presence of population 4 (see Theorem 3).

7. Optimal control In this section we introduce an optimal control problem associated to system (1) together with a numerical algorithm for its solution. Our objective will be to exploit the presence of the second biomass source in order to improve the levels of the biomasses x1, x3 and x4 along time by means of an impulsive control action. Therefore, from now on we consider periodic impulsive biomass injections in the system to be assigned along the optimization horizon. Defining t1 as the first impulsive instant, we define the following optimal control problem: min

hðt k Þ; 8 k A ½1;n  1

s:t:

J ðxðt 2 Þ; …; xðt n Þ; hðt 1 Þ; …; hðt n  1 ÞÞ

ð1Þ and ð6Þ

ð27Þ

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

where J ðxðt 2 Þ; …; xðt n Þ; hðt 1 Þ; …; hðt n  1 ÞÞ ¼

nX 1 k¼1

"

X



n 2

αi xi ðt kþ1 Þ xi

2743

#

! þ βh2 ðt k Þ

i ¼ 1;3;4

In this formulation, αi ; β are weights and xni Z 0 are target values for the respective state variables, for i ¼ 1; 3; 4, while n40 is the number of periods considered in the optimization problem. The algorithm proposed here to solve the optimal control problem (27) is based on a partition of the time horizon ½0; nP in stages separated by the impulsive events and exploits the Iterative Dynamic Programming (IDP) approach [43], to find the optimal solution. The algorithm refers to Bellman's principle of optimality and is defined by the following steps, for a given initial condition xð0þ Þ40. 1. Given the initial condition xð0þ Þ, integrate system (1) and calculate xðt 1 Þ accordingly; fix  positive constants N and M and an ordered set H 1 ¼ hi1 ; iA ½1; N  : 0 ¼ h11 o…o x4 ðt 1 Þ g of feasible values for hðt 1 Þ; correspondingly, compute N values xðiÞ ðt þ hN1 ¼ h dþx 1Þ 4 ðt 1 Þ assumed by xðt þ 1 Þ in the different cases, for i ¼ 1; …; N. 2. For all k A ½0; n 2 and for all iA ½1; N, integrate system (1) along ðk þ 1Þ-th time stage starting from xðiÞ ðt þ k Þ, choose hðt kþ1 Þ as the i-th  value in the ordered set  xðiÞ ðt Þ ðiÞ ij iN 4 kþ1 H kþ1 ¼ hkþ1 ; jA ½1; N  : 0 ¼ hi1 and calculate xðiÞ ðt þ ðiÞ kþ1 o⋯ohkþ1 ¼ h kþ1 Þ accorddþx4 ðt kþ1 Þ

xðiÞ ðt þ k Þ; ðiÞ þ

8 k A ½0; n 1 and 8 iA ½1; N, associated to each ingly. Finally, store the states impulsive input sequence, where x ð0 Þ ¼ xð0þ Þ; 8 iA ½1; N. ðiÞ þ 3. Referring to stage n, integrate system (1) starting from P each state x ðtnn2 1 Þ, i A½1; N computed at step 2 and calculate the cost-to-go J n ðjÞ ¼ l ¼ 1;3;4 αl ðxl ðt n Þ  xl Þ . 4. Step backward to stage n  1 and, for all iA ½1; N, integrate Eq. (1) starting from xðiÞ ðt þ n  2Þ

Fig. 2. System (1) exhibiting an asymptotic oscillatory regime in correspondence to the prescribed sufficient conditions for persistence.

2744

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

computed at step 2 along the stage; furthermore, for each j A ½1; M, choose hðt n  1 Þ as the   j-th ij i1 iM ðiÞ value in the ordered set H~ n  1 ¼ h~ k ; jA ½1; M  : 0 ¼ h~ k o…oh~ k ¼ h

x4ðiÞ ðt n  1 Þ

dþx4ðiÞ ðt n  1 Þ

and

denote it as hij ðt n  1 Þ. Select the value uj A ½0; N such that the point xðuj Þ ðt þ n  1 Þ defined at step as an approximate trajectory continuation. 2 is closest to the resulting state value at t þ n1 Finally, for each i A ½1; N select the hvalue j A½1; M minimizing thei cumulative cost-to-go

2 P P J n  1 ðiÞ ¼ J n ðuj Þ þ nk ¼1n  2 l ¼ 1;3;4 αl xl ðt kþ1 Þ  xnl þ βhij ðt n  1 Þ ; iA ½1; M. 5. Iterate backward in time according to the previous step, until t 0 ¼ 0. Store the values of hn ðt k Þ such that hðt k Þ ¼ hn ðt k Þ; k A ½1; n; produces the trajectory with the minimum cumulative cost-to-go. 6. Until the maximum number of iterations allowed is reached, iterate the procedure from step 2 on, while reducing the range of the sets of impulse values according to the optimal control input determined at the current iteration. A numerical application of the IDP algorithm is presented in the next section. 8. Numerical examples 8.1. Permanence We test the permanent behavior of the system under the conditions provided in Theorem 4 by using the following parameters: r1 ¼ 0.05, K1 ¼ 1.5, q2 ¼ 0.01, q3 ¼ 0.01, q4 ¼ 0.005, D13 ¼ 0:25, D24 ¼ 1, D34 ¼ 0:01, θ31 ¼ 0:1, θ42 ¼ 0:5, θ43 ¼ 0:05, a31 ; a42 ; a43 ¼ 1, h ¼ 0:5 and P ¼ 365. As 13 13 m shown in Fig. 2, the intersection point ðxM 1 ; x3 Þ ¼ ðx 1 ; x 3 Þ in the region where the population 1's isocline has negative slope, as prescribed by the condition (23b). All the remaining conditions of

Fig. 3. Time evolution of system (1) displaying permanence.

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

Fig. 4. The optimal states and impulsive sequence according to the IDP algorithm.

Fig. 5. Regime maxima of the optimal state trajectories according to the IDP algorithm as functions of h.

2745

2746

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

Theorem 4 as well, and the system trends towards an oscillatory regime which is ultimately lower min bounded by ðxmin 1 ; x3 Þ  ð0:66; 0:11Þ in the phase plane ðx1 ; x3 Þ, as shown in the figure for a specific trajectory starting from initial conditions xi ð0þ Þ ¼ x 0i ¼ 0:1; 8 i and first impulsive event at t¼ P. The corresponding behavior in time is represented in Fig. 3.

8.2. Optimal control As a second case, we consider an optimal control problem for the specified system with parameters r1 ¼ 0.1, K1 ¼ 1.1, q2 ¼ 0.004, q3 ¼ 0.01, q4 ¼ 0.002, D13 ¼ D24 ¼ 0:20, D34 ¼ 0:04, θ31 ¼ θ42 ¼ θ43 ¼ 0:20, a31 ¼ a42 ¼ a43 ¼ 1 and P¼ 365. System (1) admits the following coexistence equilibrium for the populations 1, 3 and 4 alone: ð~x 1 ; x~ 3 ; x~ 4 Þ  ð0:66; 0:33; 0:20Þ. We assume that our objective is to increment the biomass level of population 4 while maintaining the biomasses of levels 1 and 3 suitably close to the equilibrium values specified above. In particular, in Eq. (27) we fix xni ¼ x~ i ; i ¼ 1; 3 and xn4 ¼ 1:2~x 4 . As for the other parameters appearing in Eq. (27), we fix n¼ 10, h ¼ 5, d¼ 1 and the weights α1 ¼ α3 ¼ 0:1, α4 ¼ 10 and β ¼ 0:1, so that we prioritize the enhancement of the x4 level with respect to the initial equilibrium computed without the impulsive source. By applying the IDP method described above we obtain the results depicted in Fig. 4, where the values of hðt k Þ fulfil constraint (6). It can be observed that the chosen control strategy proves able to enhance the biomass level of population 4, leading x4 to oscillate around the value xn4 . Also notice that, under the considered parameter selection, in Eq. (19) we have ν1 ¼ 0:05, where ν1 is defined as in Theorem 2. Therefore, the local asymptotic stability of the periodic eradicated 1 solution S would require ho0:0775, assuming a constant impulsive intensity. Instead, the regime values for hðt k Þ found by the algorithm exceed this value, excluding the asymptotic 1 13 stability of both the periodic eradicated solutions S and S and ensuring the permanence of all the system's components. Furthermore, in Fig. 5 we report the outcomes of a parametric study, based on choosing different values for h appearing in Eq. (6) and computing the corresponding regime solutions to the proposed optimal control problem over a suitable optimization horizon, obtained by means of the proposed IDP algorithm.

9. Conclusions In this paper, we studied a multi-trophic system composed of two sources, an intermediateand an upper-level consumer, subject to periodic impulses on the second source's biomass. We analyzed some key features of the proposed system, particularly the local and global stability properties of the periodic eradicated solutions of the system and the conditions for the permanence of all the populations. Furthermore, we formulated an optimal impulsive control problem for this system. In our setup, the impulsive biomass injections are assumed to be periodic and the control action consists in defining the intensity of each impulse along time. Our formulation includes the presence of a constraint expressing the dependence of the new biomass available at the end of each period on the workforce contributing to its production. We also proposed a numerical algorithm for the solution of the proposed optimal control problem, which is based on the concept of Iterative Dynamic Programming. Its effectiveness was shown through a numerical example.

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

2747

Acknowledgements The authors would like to thank Ivano Azzini from JRC-Ispra for fruitful discussion. The research was partially funded by the Milano City Council under the program “Milano per la difesa, incremento e valorizzazione della biodiversità. Contributi a favore della solidarietà e della cooperazione internazionale - anno 2009”.

References [1] C. Becker, E. Ostrom, Human ecology and resource sustainability: the importance of institutional diversity, Annu. Rev. Ecol. Syst. 26 (1995) 113–133. [2] J. Baumgärtner, G. Gilioli, G. Tikubet, A. Gutierrez, Eco-social analysis of an East African agro-pastoral system: management of tsetse and bovine trypanosomiasis, Ecol. Econ. 65 (2008) 125–135. [3] B. Fath, Distributed control in ecological networks, Ecol. Modell. 179 (2004) 235–245. [4] J. Castilla, Coastal marine communities: trends and perspectives from human-exclusion experiments, Trends Ecol. Evol. 14 (1999) 280–283. [5] E. Neumayer, The human development index and sustainability: a constructive proposal, Ecol. Econ. 39 (2001) 101–114. [6] X. Song, L. Chen, Optimal harvesting and stability for a two-species competitive system with stage structure, Math. Biosci. 170 (2001) 173–186. [7] R. Lande, S. Engen, B.-E. Saether, Optimal harvesting, economic discounting and extinction risk in fluctuating populations, Nature 11 (1994) 88–90. [8] K. Chaudhuri, A bioeconomic model of harvesting a multispecies fishery, Ecol. Modell. 32 (1986) 267–279. [9] T. Pradhan, K. Chaudhuri, A dynamic reaction model of a two-species fishery with taxation as a control instrument: a capital theoretic analysis, Ecol. Modell. 121 (1999) 1–16. [10] T. Das, R. Mukherjee, K. Chaudhuri, Harvesting of a prey–predator fishery in the presence of toxicity, Appl. Math. Modell. 33 (2009) 2282–2292. [11] T. Christiaans, T. Eichner, R. Pethig, Optimal pest control in agriculture, J. Econ. Dyn. Control 31 (2007) 3965–3985. [12] M. Rafikov, J. Balthazar, H. von Bremen, Mathematical modeling and control of population systems: applications in biological pest control, Appl. Math. Comput. 200 (2008) 557–573. [13] L. Sun, H. Xiao, D. Yang, S. Li, Fisheries ecosystem management based on optimization algorithm, in: International Conference on Environmental Science and Information Application Technology, vol. 3, 2009, pp. 19–22. [14] M. Rafikov, T. Angelelli, Optimization of biological pest control of sugarcane borer, in: 2009 IEEE International Conference on Control Applications and International Symposium on Intelligent Control, 2009, pp. 1254–1258. [15] A. Leung, S. Stojanovic, Optimal control for elliptic Volterra–Lotka type equations, J. Math. Anal. Appl. 173 (1993) 603–619. [16] A. El-Gohary, M. Yassen, Optimal control and synchronization of Lotka–Volterra model, Chaos Solitons Fractals 12 (2001) 2087–2093. [17] V. Křivan, J. Eisner, Optimal foraging and predator–prey dynamics III, Theor. Popul. Biol. 63 (2003) 269–279. [18] Y. Shastri, U. Diwekar, Sustainable ecosystem management using optimal control theory: part 1 (deterministic systems), J. Theor. Biol. 241 (2006) 506–521. [19] Y. Shastri, U. Diwekar, Sustainable ecosystem management using optimal control theory: part 2 (stochastic systems), J. Theor. Biol. 241 (2006) 522–532. [20] N. Apreutesei, An optimal control problem for a prey–predator system with a general functional response, Appl. Math. Lett. 22 (2009) 1062–1065. [21] L. Galbusera, S. Pasquali, G. Gilioli, Stability and optimal control for some classes of tritrophic systems, Math. Biosci. Eng. 11 (2014) 257–283. [22] S. Zhang, L. Chen, Chaos in three species food chain system with impulsive perturbations, Chaos Solitons Fractals 24 (2005) 73–83. [23] S. Zhang, F. Wang, L. Chen, A food chain model with impulsive perturbations and Holling IV functional response, Chaos Solitons Fractals 26 (2005) 855–866.

2748

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

[24] W. Wang, J. Shen, J.J. Nieto, Permanence and Periodic Solution of Predator-Prey System with Holling Type Functional Response and Impulses, Discrete Dynamics in Nature and Society 2007 (2007) 15 Article ID 81756. [25] B. Liu, Z. Teng, W. Liu, Dynamic behaviors of the periodic Lotka–Volterra competing system with impulsive perturbations, Chaos Solitons Fractals 31 (2007) 356–370. [26] P. Georgescu, G. Moroşanu, Impulsive perturbations of a three-trophic prey-dependent food chain system, Math. Comput. Modell. 48 (2008) 975–997. [27] E. Braverman, R. Mamdani, Continuous versus pulse harvesting for population models in constant and variable environment, J. Math. Biol. 57 (2008) 413–434. [28] S. Zhang, D. Tan, Permanence in a food chain system with impulsive perturbations, Chaos Solitons Fractals 40 (2009) 392–400. [29] S. Nundloll, L. Mailleret, F. Grognard, Two models of interfering predators in impulsive biological control, J. Biol. Dyn. 4 (2010) 102–114. [30] P. Georgescu, H. Zhang, An impulsively controlled pest management model with n predator species and a common prey, Biosystems 110 (2012) 162–170. [31] A. Bensoussan, J.-L. Lions, Impulse Control and Quasivariational Inequalities, Gauthier-Villars, Paris, France, 1984. [32] T. Yang, Impulsive Control Theory, Lecture Notes in Control and Information Sciences, vol. 272, Springer, Berlin, Germany, 2001. [33] Z. Li, X. Yang, H. Gao, Autonomous impulsive rendezvous for spacecraft under orbital uncertainty and thruster faults, Journal of the Franklin Institute 350 (9) (2013) 2455–2473. [34] H. Jiang, J. Yu, C. Zhou, Consensus of multi-agent linear dynamic systems via impulsive control protocols, Int. J. Syst. Sci. 42 (2011) 967–976. [35] S. Baccarin, Optimal impulse control for a multidimensional cash management system with generalized cost functions, Eur. J. Oper. Res. 196 (2009) 198–206. [36] L. Zou, Z. Xiong, Z. Shu, The dynamics of an eco-epidemic model with distributed time delay and impulsive control strategy, J. Frankl. Inst. 348 (2011) 2332–2349. [37] H. Yu, S. Zhong, R.P. Agarwal, S.K. Sen, Effect of seasonality on the dynamical behavior of an ecological system with impulsive control strategy, J. Frankl. Inst. 348 (2011) 652–670. [38] S. Hou, K. Wong, Optimal impulsive control problem with application to human immunodeficiency virus treatment, J. Optim. Theory Appl. 151 (2011) 385–401. [39] L. Wang, Y. Xie, J. Fu, The dynamics of natural mortality for pest control model with impulsive effect, J. Frankl. Inst. 350 (2013) 1443–1461. [40] Y. Xiao, D. Cheng, H. Qin, Optimal impulsive control in periodic ecosystem, Syst. Cont. Lett. 55 (2006) 558–565. [41] L. Mailleret, F. Grognard, Global stability and optimisation of a general impulsive biological control model, Math. Biosci. 221 (2009) 91–100. [42] D. Xiao, S. Ruan, Global analysis in a predator–prey system with nonmonotonic functional response, SIAM J. Appl. Math. 61 (2001) 1445–1472. [43] R. Luus, Iterative Dynamic Programming, 1st edition, CRC Press, Inc., Boca Raton, FL, USA, 2000. [44] R. Luus, Piecewise linear continuous optimal control by iterative dynamic programming, Ind. Eng. Chem. Res. 32 (1993) 859–865. [45] S. Dadebo, K. McAuley, Dynamic optimization of constrained chemical engineering problems using dynamic programming, Comput. Chem. Eng. 19 (1995) 513–525. [46] B. Bojkov, R. Luus, Optimal control of nonlinear systems with unspecified final times, Chem. Eng. Sci. 51 (1996) 905–919. [47] Y. Lei, S. Li, X. Zhang, Q. Zhang, L. Guo, Optimal control of polymer flooding based on mixed-integer iterative dynamic programming, Int. J. Control 84 (2011) 1903–1914. [48] H.-G. Wahl, F. Gauterin, An iterative dynamic programming approach for the global optimal control of hybrid electric vehicles under real-time constraints, in: 2013 IEEE Intelligent Vehicles Symposium, 2013, pp. 592–597. [49] P. Ziemniak, Mobile sensor motion planning for identification of a contamination source using iterative dynamic programming, Solid State Phenom. 198 (2013) 102–107. [50] D.D. Baı̆nov and P.S. Simeonov, Impulsive Differential Equations: Periodic Solutions and Applications, Pitman Monographs and Surveys in Pure and Applied Mathematics, vol. 66, Longman Scientific & Technical, John Wiley & Sons, Harlow, UK, New York, NY, USA, 1993. [51] M. Myerscough, B. Gray, W. Hogarth, J. Norbury, An analysis of an ordinary differential equation model for a twospecies predator–prey system with harvesting and stocking, J. Math. Biol. 30 (1992) 389–411.

L. Galbusera, S. Pasquali / Journal of the Franklin Institute 352 (2015) 2728–2749

2749

[52] X. Tian, R. Xu, Global dynamics of a predator–prey system with Holling type II functional response, Nonlinear Anal.: Modell. Control 16 (2011) 242–253. [53] S.-B. Hsu, On global stability of a predator–prey system, Math. Biosci. 39 (1978) 1–10.