- Email: [email protected]
Locating the equilibrium points of a predator–prey model by means of affine state feedback
ARTICLE IN PRESS
Journal of the Franklin Institute 345 (2008) 489–498 www.elsevier.com/locate/jfranklin
Locating the equilibrium points of a predator–prey model by means of affine state feedback Wieslaw Krajewskia,, Umberto Viarob a
Systems Research Institute, Polish Academy of Sciences, ul. Newelska 6, 01-447 Warsaw, Poland b Department of Electrical, Management and Mechanical Engineering, University of Udine, Via delle Scienze 208, 33100 Udine, Italy
Received 31 January 2006; received in revised form 3 November 2007; accepted 15 February 2008
Abstract A predator–prey model with prey-dependent functional response is considered. The set of all points in the positive quadrant of the state plane that can be made equilibrium points by means of an affine state-feedback control law is determined, and the values of the control parameters ensuring the desired equilibria are provided. It is shown how the asymptotic stability of the equilibrium points depends on simple geometric conditions. The problem of stabilizing unstable equilibrium points is also briefly discussed. r 2008 The Franklin Institute. Published by Elsevier Ltd. All rights reserved. Keywords: Predator–prey model; State feedback; Equilibrium point; Stability; Stabilization
1. Introduction It has been observed (cf., e.g. [1–3]) that the dynamics of many renewable resources, as well as the kinetics of certain chemical reactions, can be described by a nonlinear predator– prey model whose state equations are dx xðtÞyðtÞ ¼ axðtÞ bx2 ðtÞ c uðtÞ, dt 1 þ dxðtÞ
Corresponding author. Tel.: +48 22 8361990; fax: +48 22 8372772.
E-mail address: [email protected] (W. Krajewski). 0016-0032/$32.00 r 2008 The Franklin Institute. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.jfranklin.2008.02.001
(1)
ARTICLE IN PRESS 490
W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
dy xðtÞyðtÞ ¼ eyðtÞ fy2 ðtÞ þ g vðtÞ, dt 1 þ dxðtÞ
(2)
where the state components xðtÞ and yðtÞ denote, respectively, prey and predator concentrations, parameters a, b, c, d, e, f, g are positive constants, and uðtÞ, vðtÞ represent exogenous inputs. By their intrinsic nature, only positive values of x and y are of interest. Model (1)–(2) corresponds to the Rosenzweig–MacArthur model [4] with the addition of the quadratic term fy2 ðtÞ in Eq. (2) accounting for the fighting among predators [5], and can be regarded as an improved version of the classical Lotka–Volterra model [6–8] in which b, d and f were equal to zero. In particular, parameters a and e represent, respectively, the prey intrinsic growth rate and the predator intrinsic death rate, the ratio a=b is the prey carrying capacity measuring the steady-state prey density in the absence of predators, and parameter d accounts for the saturation effect according to the functional response theory [9]; precisely, the fractional terms at the right-hand sides of Eqs. (1) and (2) adopt a Holling type II prey-dependent functional response pðxÞ:¼x=ð1 þ dxÞ. The stability properties of predator–prey models either in the absence of control actions or in the presence of constant inputs have been studied rather extensively. The influence of linear state-feedback control has recently been considered in [10] with reference to the Holling type II functional response, and in [11,12] for different kinds of state equations. The main advantages of linear state-feedback control of predator–prey systems are the following: (i) a control effort proportional to the amount of population, that is, to the system state, can be implemented easily in most practical cases (in fisheries, the control effort is related to the number of fishing boats operating simultaneously) and does not require continual population monitoring (the catch itself provides a measure of the population), (ii) feedback control is fairly insensitive to disturbances and parameter variations. In Sections 2 and 3, we pursue the analysis of [10], where a list of analytic conditions that should ensure the existence of one or two positive equilibria as well as the presence of a limit cycle was provided. Precisely, using a graphic approach we determine the set of all points in the first quadrant of the state plane that can be made equilibrium points by means of an affine state-feedback control strategy of the form uðtÞ ¼ h þ kxðtÞ,
(3)
vðtÞ ¼ myðtÞ,
(4)
where the proportionality coefficients k and m represent the control efforts relating each state component to the derivative of the same state component, and h accounts for constant prey harvesting ðh40Þ or feeding ðho0Þ. As is often the case in practice, no constant predator harvesting is considered. Unless otherwise stated, the control parameters k and m are assumed to be nonnegative (even if, in some cases, it might be profitable to introduce preys or predators into the system). Once the admissible positive equilibria are found, it will be possible to choose the most convenient one on the basis of steady-state yield (which can be evaluated easily if the system is maintained at an equilibrium state). To facilitate the design task, simple formulas for finding the values of the control parameters that ensure a given equilibrium are provided in Section 2. Note that the analytic conditions presented in [10] are not exhaustive and do not allow us to easily determine the control parameters. Section 3 analyzes the stability of the equilibrium points and presents simple geometric conditions for their
ARTICLE IN PRESS W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
491
stability. It is also shown how the number and nature of the positive equilibrium points depend on the model parameters: in particular, an example of predator–prey model with three positive equilibrium points (a case that is not considered in [10]) is presented. The problem of stabilizing unstable equilibrium points is finally mentioned. 2. Equilibria Replacing uðtÞ and vðtÞ in Eqs. (1) and (2) by their respective expressions Eqs. (3) and (4), we get the homogeneous state equations dx xðtÞyðtÞ : ^ ¼ F ðx; yÞ¼axðtÞ bx2 ðtÞ c h, dt 1 þ dxðtÞ
(5)
dy xðtÞyðtÞ : ¼ Gðx; yÞ¼ e^yðtÞ fy2 ðtÞ þ g , dt 1 þ dxðtÞ
(6)
where : ^ a k, a¼
(7)
: e^¼e þ m.
(8)
The equilibrium points can be determined by setting dx=dt ¼ 0 and dy=dt ¼ 0. Therefore, they correspond to the intersections, if any, of the curves described analytically by F ðx; yÞ ¼ 0 and Gðx; yÞ ¼ 0. 2.1. Case of h ¼ 0 If the constant harvesting term h is equal to 0, F ðx; yÞ ¼ 0 both for x ¼ 0 (vertical axis) and for a^ bx cy=ð1 þ dxÞ ¼ 0, that is, : 1 y ¼ jðxÞ¼ ða^ bxÞð1 þ dxÞ, (9) c which is the equation of a parabola whose axis is vertical and whose vertex (maximum of y) is at ^ 2 ^ b ðb þ adÞ ad ; yv ¼ . 2bd 4bcd This parabola intersects the y-axis only at xv ¼
a^ c and the x-axis at yp ¼
(10)
(11)
1 a^ xp1 ¼ ; xp2 ¼ . (12) d b From Eqs. (7) and (8) it follows that, as the control parameter k increases, vertex (10) moves to the left and downwards. Note that: (i) xp1 is independent of k, whereas yp and xp2 decrease as k increases (and, thus, as a^ decreases), and (ii) yp ¼ 0 for k ¼ a so that for k4a the parabola does not cross the first quadrant and no positive equilibrium point exists.
ARTICLE IN PRESS W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
492
Concerning the predator rate, Gðx; yÞ ¼ 0 both for y ¼ 0 and for ^e fy þ gx= ð1 þ dxÞ ¼ 0, that is, gx : 1 e^ , (13) y ¼ gðxÞ¼ f 1 þ dx which is the equation of a hyperbola whose (mutually orthogonal) asymptotes are 1 1 g e^ . x¼ ; y¼ d f d
(14)
Note that: (i) the first of Eq. (14) is independent of control parameter m, whereas the horizontal asymptote moves downwards as m (and, thus, e^Þ increases, (ii) the left branch of hyperbola Eq. (13) does not belong to the first quadrant, (iii) its right branch may cross this quadrant only if g4d^e, that is, moðg=dÞ e, and (iv) the intersection of the right branch with the horizontal axis occurs at e^ xh ¼ (15) g d^e so that, as m increases, xh increases from the value e . x^ h ¼ g de
(16)
Parabola (9), hyperbola (13), the vertical axis x ¼ 0 (where dx=dt ¼ 0Þ and the horizontal axis y ¼ 0 (where dy=dt ¼ 0Þ are represented in Fig. 1 with solid lines with reference to a case in which strictly positive equilibria exist. Observe that, besides the origin, nonnegative 25
20
15
10
k=0 m=0
5 R0
0
−5 −10
−8
−6
−4
−2
0
2
4
6
8
10
Fig. 1. Determination of the equilibrium points and of the state-plane region R0 where positive equilibria are possible.
ARTICLE IN PRESS W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
493
equilibria may also occur at the intersection of parabola (9) with the positive horizontal semi-axis. Taking into account Eq. (7), we see that, 8x40, the value of y in (9) becomes smaller as control parameter k increases. Therefore, in the first quadrant the parabola corresponding to k ¼ 0 (that is, to a^ ¼ aÞ is above all the parabolas corresponding to positive values of k. Similarly, from Eqs. (8) and (13) it follows that the right branch of the hyperbola for m ¼ 0 (that is, for e^ ¼ eÞ, is above the right branches of all the hyperbolas corresponding to m40. A consequence of the previous considerations is the following. Proposition 2.1. The region R0 of the first quadrant where equilibria may occur is formed by the points below both the parabola (9) with a^ ¼ a and the hyperbola (13) with e^ ¼ e (shaded area in Fig. 1). To any point ðxe ; ye Þ 2 R0 there corresponds a unique pair of control parameters ðke ; me Þ according to cye ke ¼ a bxe , (17) 1 þ dxe me ¼
gxe fye e, 1 þ dxe
(18)
where the parameters appearing at the right-hand sides only depend on the controlled process. If ðxe ; ye Þ 2 R0 , then ke and me are necessarily nonnegative but the chosen equilibrium point may turn out to be unstable: this problem is dealt with in Section 3. Observe finally that the same pair ðke ; me Þ may give rise to two further equilibria. In fact, by equating the right-hand sides of Eqs. (9) and (13), we obtain a third-degree algebraic equation in x with either one or three real roots (cf. Section 3). 2.2. Case of ha0 If the constant prey-harvesting parameter h is different from 0, the locus where F ðx; yÞ ¼ 0 is no longer formed by parabola (9) and the y-axis, but it is described by 1 þ dx ^ bx2 hÞ, ðax (19) cx whereas the locus where the predator rate vanishes is still formed by hyperbola (13) and the x-axis. The function at thepright-hand ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi side of (19) has a pole at x ¼ 0 and three zeros at 1=d (like (9)) and at ða^ a^ 2 4bhÞ=ð2bÞ, which are real for hoa^ 2 =ð4bÞ. When h40, curve (19) consists of two branches. The left branch belongs to the halfplane xo0 above parabola (9): it tends to 1 as x ! 1 and to þ1 as x ! 0 . The right branch belongs to the half-plane x40 below parabola (9), has a unique maximum and tends to 1 both as x ! 0þ and as x ! þ1 (cf. the dashed line in Fig. 2 for h ¼ 4Þ. It is easily seen that the right branch of Eq. (19) moves downwards as h increases, and is entirely contained in the fourth quadrant for h4a^ 2 =ð4bÞ. As already mentioned, if preys are not too expensive and predators are highly remunerative, it might be advantageous to allow ho0 with the purpose of increasing predator harvest, that is, the product m y. In this case, the right branch of curve (19) takes the form of the dashed line corresponding to h ¼ 7 in Fig. 2. y ¼ jðxÞ ¼
ARTICLE IN PRESS W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
494
30
25
20
15
h = −7
10
h=0 h=4
5
0 0
1
2
3
4
5
6
7
8
9
10
Fig. 2. Curve y ¼ jðxÞ for a negative, a null and a positive value of h (dashed lines) and curve y ¼ gðxÞ (solid line).
An immediate consequence of the previous considerations is the following. Proposition 2.2. Given hoa^ 2 =ð4bÞ, the state-plane region Rh of the first quadrant where equilibrium can be achieved through a suitable choice of k and m consists of the points below both curves (13) and (19). Each equilibrium state ðxe ; ye Þ 2 Rh corresponds to the value (18) of m (which is independent of h) and to the following value of k khe ¼ a bxe
cye h , 1 þ dxe xe
(20)
which, obviously, reduces to Eq. (17) for h ¼ 0. Conversely, any point of R0 (see Fig. 1) can be made an equilibrium point by setting m ¼ me and choosing one of the infinite pairs ðh; khe Þ satisfying (20). 3. Stability The stability of the equilibrium points depends on the eigenvalues of the associated Jacobian matrix 2 3 qF qF 6 qx qy 7 6 7 (21) Jðx; yÞ ¼ 6 7. 4 qG qG 5 qx qy
ARTICLE IN PRESS W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
495
From Eqs. (5), (6), (13) and (19) we obtain F ðx; yÞ ¼
cx ½jðxÞ y, 1 þ dx
Gðx; yÞ ¼ fy½gðxÞ y.
(22) (23)
By taking the derivatives of Eqs. (22) and (23) and considering that at an equilibrium point ðxe ; ye Þ we have ye ¼ jðxe Þ ¼ gðxe Þ, for ðx; yÞ ¼ ðxe ; ye Þ matrix (21) particularizes to 2 3 cxe dj cxe 6 1 þ dxe dx 1 þ dxe 7 : 7, J e ¼Jðxe ; ye Þ ¼ 6 (24) 4 5 dg fye fye dx where the derivatives are computed at x ¼ xe . The eigenvalues of J e are the roots of l2 trðJ e Þl þ detðJ e Þ ¼ 0,
(25)
which, according to Descartes’ rule, have a negative real part if and only if trðJ e Þo0;
detðJ e Þ40.
(26)
Taking into account Eqs. (24) and (26), the following proposition can be stated. Proposition 3.1. A positive equilibrium point ðxe ; ye Þ is asymptotically stable if and only if dj fye o ð1 þ dxe Þ; dx cxe
dj dg g o ¼ . dx dx f ð1 þ dxe Þ2
(27)
From a geometric point of view, the second of inequalities (27) means that the slope of the tangent to jðxÞ at x ¼ xe must be less than the slope of the tangent to gðxÞ, so that, in the neighbourhood of xe , jðxÞ4gðxÞ for xoxe and jðxÞogðxÞ for x4xe . Assuming c, f, g and d positive, at a positive equilibrium point the right-hand sides of inequalities (27) are both positive. Therefore we can state the following sufficient condition for asymptotic stability. Proposition 3.2. A positive equilibrium point ðxe ; ye Þ is asymptotically stable if dj 1 h ^ b 2bdxe o0. ¼ þ ad dx c x2e
(28)
Example 3.1 (Jing and Yang [10]). Consider the system (5)–(6) with a^ ¼ 0:9, b ¼ 0:1, c ¼ 2, d ¼ 0:5, e^ ¼ 1, f ¼ 1, g ¼ 0:4, h ¼ 0:8. The (unique) positive equilibrium point ðxe ; ye Þ ¼ ð4:93861; 0:42352Þ is asymptotically stable because dj=dx ¼ 0:056 at x ¼ xe . When h ¼ 0, the condition of Proposition 3.2 is satisfied if xe lies to the right of the vertex of parabola (9). Therefore, taking into account (7) and (10), the following proposition holds. Proposition 3.3. For h ¼ 0, a positive equilibrium point ðxe ; ye Þ is asymptotically stable if xe 4
ða kÞd b . 2bd
(29)
ARTICLE IN PRESS 496
W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
The locus described by the vertices ðxv ; yv Þ of parabolas (9) as k varies from 0 to the positive value ða bÞ=d is represented by a dashed line in Fig. 1. If the equilibrium point is located to the right of this locus, then it will necessarily be asymptotically stable. In the case of Fig. 1 all of the equilibrium points internal to R0 turn out to be asymptotically stable. Concerning the number of strictly positive equilibrium points corresponding to a given set of model parameters in Eqs. (5) and (6), let us consider the equation obtained by equating Eqs. (13) and (19) : FðxÞ¼jðxÞ gðxÞ ¼ 0, (30) whose roots comprise all the strictly positive equilibrium points (as well as other equilibrium points, if any, whose coordinates are not strictly positive). Function FðxÞ is rational with numerator given by : N FðxÞ ¼f ð1 þ dxÞ2 ðax bx2 hÞ cx½ðg edÞx e (31) and denominator by : DF ðxÞ¼cfxð1 þ dxÞ,
(32)
which is a second-degree polynomial that vanishes for x ¼ 0 and x ¼ 1=do0. Since the degree of N FðxÞ is four, the number of its real zeros (not necessarily distinct or corresponding to positive values of both x and y) may be 4, 2 or 0. If, however, h ¼ 0, then x can be factored out of the numerator thus canceling the corresponding factor in the denominator and the number of positive zeros of FðxÞ may be either 3 or 1 (by counting multiplicities). In the latter case, if the horizontal asymptote of hyperbola (13) (see the second of Eq. (14)) is sufficiently low, the unique intersection occurs at a point where the slope of the tangent to parabola (9) is negative (stable equilibrium); if, instead, the horizontal asymptote is sufficiently high, the intersection occurs for negative values of x and no positive equilibrium exists; for intermediate values of the second of Eq. (14) the intersection occurs when x is positive but less than xv , and the stability of the equilibrium can be checked using the first of conditions (27) only. Three intersections of Eq. (9) with Eq. (13) are also possible but rare in practice. Example 3.2. Consider the system (5)–(6) with a^ ¼ 1:8535, b ¼ 0:25, c ¼ 2:2929, d ¼ 1:4142, e^ ¼ 1, f ¼ 2, g ¼ 10:5666, h ¼ 0. This system exhibits three positive equilibrium points, as shown in Fig. 3. According to Proposition 3.1, the equilibria at ð1; 1:6884Þ and ð3; 2:5233Þ are asymptotically stable, whereas the equilibrium at ð2; 2:2600Þ is unstable. If an equilibrium point, chosen according to convenience or profit criteria, turns out to be unstable, resort can be made to standard feedback-stabilization techniques, as shown in Fig. 4. To this purpose, the preliminary determination of a linearized model of the process may be necessary. Of course, it is important to evaluate the domain of attraction of the stabilized equilibrium point as well as the feasibility of the resulting stabilization policy: for instance, this policy might require the insertion of preys and/or predators, which is often forbidden. If more or all equilibrium points must be stabilized simultaneously, one can exploit the method suggested in [13] since the number of inputs of the considered system equals its state dimension: in this case, the structure of the (nonlinear) stabilizing controller does not change with the equilibrium point.
ARTICLE IN PRESS W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
497
3
2.5
(3,2.5233) (2,2.2600)
2 (1,1.6884) 1.5
1
0.5
0
0
0.5
1
1.5
2
2.5
3
3.5
4
Fig. 3. Predator–prey model with three positive equilibria: the dashed line represents y ¼ jðxÞ and the solid line represents y ¼ gðxÞ.
Fig. 4. Feedback stabilization scheme.
4. Conclusions The convex state-plane region where positive equilibria of predator–prey models of form (1) and (2) can be assigned through an affine state-feedback control law like Eqs. (3) and
ARTICLE IN PRESS 498
W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
(4), has been determined. The values of the control parameters corresponding to a given equilibrium state can easily be found (see Eqs. (17), (18) and (20)). Geometrical considerations have allowed us to derive simple conditions for the stability of a positive equilibrium state (see Eqs. (27) and (29)). The adopted approach integrates recent alternative results [10] regarding the impact of state-feedback control on predator–prey systems with functional response and, in these authors’ opinion, provides a more practical design tool. References [1] H.I. Freedman, Deterministic Mathematical Models in Population Ecology, Marcel Dekker, New York, 1980. [2] A. Cornish-Bowden, Fundamentals of Enzyme Kinetics, Portland, London, 1995. [3] G.F. Fussmann, S.P. Ellner, K.W. Shertzer, N.G. Hairston Jr., Crossing the Hopf bifurcation in a live predator–prey system, Science 290 (2000) 1358–1360. [4] M.L. Rosenzweig, R.H. MacArthur, Graphical representation and stability conditions of predator–prey interactions, Am. Naturalist 97 (1963) 209–223. [5] E.A. McGehee, E. Peacock–Lopez, Turing patterns in a modified Lotka–Volterra model, Phys. Lett. A 342 (2005) 90–98. [6] A.J. Lotka, Elements of Physical Biology, Williams and Wilkins, Baltimore, MD, 1925. [7] V. Volterra, Lecons sur la Theorie Methematique de la Lutte pour la Vie, Gauthier-Vilars, Paris, 1931. [8] F.M. Scudo, J.R. Ziegler, The Golden Age of Theoretical Ecology, Lecture Notes in Biomathematics, vol. 22, Springer, New York, 1978, pp. 1923–1940. [9] C.S. Holling, Some characteristics of simple types of predation and parasitism, Can. Entomol. 91 (1959) 385–398. [10] H. Jing, Z. Yang, The impact of state feedback control on a predator–prey model with functional response, Discrete Continuous Dynam. Syst.—Ser. B 4 (2004) 607–614. [11] M. Fan, Y. Kuang, Dynamics of a nonautonomous predator–prey system with the Beddington-De Angelis functional response, J. Math. Anal. Appl. 295 (2004) 15–39. [12] P. De Leenheer, D. Angeli, E.D. Sontag, On predator–prey systems and small gain theorems, Math. Biosci. Eng. 2 (2005) 25–42. [13] A. Ferrante, W. Krajewski, A. Lepschy, U. Viaro, Simultaneous stabilization of multiple equilibrium points, in: Proceedings of the 10th IEEE International Conference Methods and Models in Automation and Robotics (MMAR’04), Miedzyzdroje, Poland, August 2004, pp. 321–324.
Journal of the Franklin Institute 345 (2008) 489–498 www.elsevier.com/locate/jfranklin
Locating the equilibrium points of a predator–prey model by means of affine state feedback Wieslaw Krajewskia,, Umberto Viarob a
Systems Research Institute, Polish Academy of Sciences, ul. Newelska 6, 01-447 Warsaw, Poland b Department of Electrical, Management and Mechanical Engineering, University of Udine, Via delle Scienze 208, 33100 Udine, Italy
Received 31 January 2006; received in revised form 3 November 2007; accepted 15 February 2008
Abstract A predator–prey model with prey-dependent functional response is considered. The set of all points in the positive quadrant of the state plane that can be made equilibrium points by means of an affine state-feedback control law is determined, and the values of the control parameters ensuring the desired equilibria are provided. It is shown how the asymptotic stability of the equilibrium points depends on simple geometric conditions. The problem of stabilizing unstable equilibrium points is also briefly discussed. r 2008 The Franklin Institute. Published by Elsevier Ltd. All rights reserved. Keywords: Predator–prey model; State feedback; Equilibrium point; Stability; Stabilization
1. Introduction It has been observed (cf., e.g. [1–3]) that the dynamics of many renewable resources, as well as the kinetics of certain chemical reactions, can be described by a nonlinear predator– prey model whose state equations are dx xðtÞyðtÞ ¼ axðtÞ bx2 ðtÞ c uðtÞ, dt 1 þ dxðtÞ
Corresponding author. Tel.: +48 22 8361990; fax: +48 22 8372772.
E-mail address: [email protected] (W. Krajewski). 0016-0032/$32.00 r 2008 The Franklin Institute. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.jfranklin.2008.02.001
(1)
ARTICLE IN PRESS 490
W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
dy xðtÞyðtÞ ¼ eyðtÞ fy2 ðtÞ þ g vðtÞ, dt 1 þ dxðtÞ
(2)
where the state components xðtÞ and yðtÞ denote, respectively, prey and predator concentrations, parameters a, b, c, d, e, f, g are positive constants, and uðtÞ, vðtÞ represent exogenous inputs. By their intrinsic nature, only positive values of x and y are of interest. Model (1)–(2) corresponds to the Rosenzweig–MacArthur model [4] with the addition of the quadratic term fy2 ðtÞ in Eq. (2) accounting for the fighting among predators [5], and can be regarded as an improved version of the classical Lotka–Volterra model [6–8] in which b, d and f were equal to zero. In particular, parameters a and e represent, respectively, the prey intrinsic growth rate and the predator intrinsic death rate, the ratio a=b is the prey carrying capacity measuring the steady-state prey density in the absence of predators, and parameter d accounts for the saturation effect according to the functional response theory [9]; precisely, the fractional terms at the right-hand sides of Eqs. (1) and (2) adopt a Holling type II prey-dependent functional response pðxÞ:¼x=ð1 þ dxÞ. The stability properties of predator–prey models either in the absence of control actions or in the presence of constant inputs have been studied rather extensively. The influence of linear state-feedback control has recently been considered in [10] with reference to the Holling type II functional response, and in [11,12] for different kinds of state equations. The main advantages of linear state-feedback control of predator–prey systems are the following: (i) a control effort proportional to the amount of population, that is, to the system state, can be implemented easily in most practical cases (in fisheries, the control effort is related to the number of fishing boats operating simultaneously) and does not require continual population monitoring (the catch itself provides a measure of the population), (ii) feedback control is fairly insensitive to disturbances and parameter variations. In Sections 2 and 3, we pursue the analysis of [10], where a list of analytic conditions that should ensure the existence of one or two positive equilibria as well as the presence of a limit cycle was provided. Precisely, using a graphic approach we determine the set of all points in the first quadrant of the state plane that can be made equilibrium points by means of an affine state-feedback control strategy of the form uðtÞ ¼ h þ kxðtÞ,
(3)
vðtÞ ¼ myðtÞ,
(4)
where the proportionality coefficients k and m represent the control efforts relating each state component to the derivative of the same state component, and h accounts for constant prey harvesting ðh40Þ or feeding ðho0Þ. As is often the case in practice, no constant predator harvesting is considered. Unless otherwise stated, the control parameters k and m are assumed to be nonnegative (even if, in some cases, it might be profitable to introduce preys or predators into the system). Once the admissible positive equilibria are found, it will be possible to choose the most convenient one on the basis of steady-state yield (which can be evaluated easily if the system is maintained at an equilibrium state). To facilitate the design task, simple formulas for finding the values of the control parameters that ensure a given equilibrium are provided in Section 2. Note that the analytic conditions presented in [10] are not exhaustive and do not allow us to easily determine the control parameters. Section 3 analyzes the stability of the equilibrium points and presents simple geometric conditions for their
ARTICLE IN PRESS W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
491
stability. It is also shown how the number and nature of the positive equilibrium points depend on the model parameters: in particular, an example of predator–prey model with three positive equilibrium points (a case that is not considered in [10]) is presented. The problem of stabilizing unstable equilibrium points is finally mentioned. 2. Equilibria Replacing uðtÞ and vðtÞ in Eqs. (1) and (2) by their respective expressions Eqs. (3) and (4), we get the homogeneous state equations dx xðtÞyðtÞ : ^ ¼ F ðx; yÞ¼axðtÞ bx2 ðtÞ c h, dt 1 þ dxðtÞ
(5)
dy xðtÞyðtÞ : ¼ Gðx; yÞ¼ e^yðtÞ fy2 ðtÞ þ g , dt 1 þ dxðtÞ
(6)
where : ^ a k, a¼
(7)
: e^¼e þ m.
(8)
The equilibrium points can be determined by setting dx=dt ¼ 0 and dy=dt ¼ 0. Therefore, they correspond to the intersections, if any, of the curves described analytically by F ðx; yÞ ¼ 0 and Gðx; yÞ ¼ 0. 2.1. Case of h ¼ 0 If the constant harvesting term h is equal to 0, F ðx; yÞ ¼ 0 both for x ¼ 0 (vertical axis) and for a^ bx cy=ð1 þ dxÞ ¼ 0, that is, : 1 y ¼ jðxÞ¼ ða^ bxÞð1 þ dxÞ, (9) c which is the equation of a parabola whose axis is vertical and whose vertex (maximum of y) is at ^ 2 ^ b ðb þ adÞ ad ; yv ¼ . 2bd 4bcd This parabola intersects the y-axis only at xv ¼
a^ c and the x-axis at yp ¼
(10)
(11)
1 a^ xp1 ¼ ; xp2 ¼ . (12) d b From Eqs. (7) and (8) it follows that, as the control parameter k increases, vertex (10) moves to the left and downwards. Note that: (i) xp1 is independent of k, whereas yp and xp2 decrease as k increases (and, thus, as a^ decreases), and (ii) yp ¼ 0 for k ¼ a so that for k4a the parabola does not cross the first quadrant and no positive equilibrium point exists.
ARTICLE IN PRESS W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
492
Concerning the predator rate, Gðx; yÞ ¼ 0 both for y ¼ 0 and for ^e fy þ gx= ð1 þ dxÞ ¼ 0, that is, gx : 1 e^ , (13) y ¼ gðxÞ¼ f 1 þ dx which is the equation of a hyperbola whose (mutually orthogonal) asymptotes are 1 1 g e^ . x¼ ; y¼ d f d
(14)
Note that: (i) the first of Eq. (14) is independent of control parameter m, whereas the horizontal asymptote moves downwards as m (and, thus, e^Þ increases, (ii) the left branch of hyperbola Eq. (13) does not belong to the first quadrant, (iii) its right branch may cross this quadrant only if g4d^e, that is, moðg=dÞ e, and (iv) the intersection of the right branch with the horizontal axis occurs at e^ xh ¼ (15) g d^e so that, as m increases, xh increases from the value e . x^ h ¼ g de
(16)
Parabola (9), hyperbola (13), the vertical axis x ¼ 0 (where dx=dt ¼ 0Þ and the horizontal axis y ¼ 0 (where dy=dt ¼ 0Þ are represented in Fig. 1 with solid lines with reference to a case in which strictly positive equilibria exist. Observe that, besides the origin, nonnegative 25
20
15
10
k=0 m=0
5 R0
0
−5 −10
−8
−6
−4
−2
0
2
4
6
8
10
Fig. 1. Determination of the equilibrium points and of the state-plane region R0 where positive equilibria are possible.
ARTICLE IN PRESS W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
493
equilibria may also occur at the intersection of parabola (9) with the positive horizontal semi-axis. Taking into account Eq. (7), we see that, 8x40, the value of y in (9) becomes smaller as control parameter k increases. Therefore, in the first quadrant the parabola corresponding to k ¼ 0 (that is, to a^ ¼ aÞ is above all the parabolas corresponding to positive values of k. Similarly, from Eqs. (8) and (13) it follows that the right branch of the hyperbola for m ¼ 0 (that is, for e^ ¼ eÞ, is above the right branches of all the hyperbolas corresponding to m40. A consequence of the previous considerations is the following. Proposition 2.1. The region R0 of the first quadrant where equilibria may occur is formed by the points below both the parabola (9) with a^ ¼ a and the hyperbola (13) with e^ ¼ e (shaded area in Fig. 1). To any point ðxe ; ye Þ 2 R0 there corresponds a unique pair of control parameters ðke ; me Þ according to cye ke ¼ a bxe , (17) 1 þ dxe me ¼
gxe fye e, 1 þ dxe
(18)
where the parameters appearing at the right-hand sides only depend on the controlled process. If ðxe ; ye Þ 2 R0 , then ke and me are necessarily nonnegative but the chosen equilibrium point may turn out to be unstable: this problem is dealt with in Section 3. Observe finally that the same pair ðke ; me Þ may give rise to two further equilibria. In fact, by equating the right-hand sides of Eqs. (9) and (13), we obtain a third-degree algebraic equation in x with either one or three real roots (cf. Section 3). 2.2. Case of ha0 If the constant prey-harvesting parameter h is different from 0, the locus where F ðx; yÞ ¼ 0 is no longer formed by parabola (9) and the y-axis, but it is described by 1 þ dx ^ bx2 hÞ, ðax (19) cx whereas the locus where the predator rate vanishes is still formed by hyperbola (13) and the x-axis. The function at thepright-hand ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi side of (19) has a pole at x ¼ 0 and three zeros at 1=d (like (9)) and at ða^ a^ 2 4bhÞ=ð2bÞ, which are real for hoa^ 2 =ð4bÞ. When h40, curve (19) consists of two branches. The left branch belongs to the halfplane xo0 above parabola (9): it tends to 1 as x ! 1 and to þ1 as x ! 0 . The right branch belongs to the half-plane x40 below parabola (9), has a unique maximum and tends to 1 both as x ! 0þ and as x ! þ1 (cf. the dashed line in Fig. 2 for h ¼ 4Þ. It is easily seen that the right branch of Eq. (19) moves downwards as h increases, and is entirely contained in the fourth quadrant for h4a^ 2 =ð4bÞ. As already mentioned, if preys are not too expensive and predators are highly remunerative, it might be advantageous to allow ho0 with the purpose of increasing predator harvest, that is, the product m y. In this case, the right branch of curve (19) takes the form of the dashed line corresponding to h ¼ 7 in Fig. 2. y ¼ jðxÞ ¼
ARTICLE IN PRESS W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
494
30
25
20
15
h = −7
10
h=0 h=4
5
0 0
1
2
3
4
5
6
7
8
9
10
Fig. 2. Curve y ¼ jðxÞ for a negative, a null and a positive value of h (dashed lines) and curve y ¼ gðxÞ (solid line).
An immediate consequence of the previous considerations is the following. Proposition 2.2. Given hoa^ 2 =ð4bÞ, the state-plane region Rh of the first quadrant where equilibrium can be achieved through a suitable choice of k and m consists of the points below both curves (13) and (19). Each equilibrium state ðxe ; ye Þ 2 Rh corresponds to the value (18) of m (which is independent of h) and to the following value of k khe ¼ a bxe
cye h , 1 þ dxe xe
(20)
which, obviously, reduces to Eq. (17) for h ¼ 0. Conversely, any point of R0 (see Fig. 1) can be made an equilibrium point by setting m ¼ me and choosing one of the infinite pairs ðh; khe Þ satisfying (20). 3. Stability The stability of the equilibrium points depends on the eigenvalues of the associated Jacobian matrix 2 3 qF qF 6 qx qy 7 6 7 (21) Jðx; yÞ ¼ 6 7. 4 qG qG 5 qx qy
ARTICLE IN PRESS W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
495
From Eqs. (5), (6), (13) and (19) we obtain F ðx; yÞ ¼
cx ½jðxÞ y, 1 þ dx
Gðx; yÞ ¼ fy½gðxÞ y.
(22) (23)
By taking the derivatives of Eqs. (22) and (23) and considering that at an equilibrium point ðxe ; ye Þ we have ye ¼ jðxe Þ ¼ gðxe Þ, for ðx; yÞ ¼ ðxe ; ye Þ matrix (21) particularizes to 2 3 cxe dj cxe 6 1 þ dxe dx 1 þ dxe 7 : 7, J e ¼Jðxe ; ye Þ ¼ 6 (24) 4 5 dg fye fye dx where the derivatives are computed at x ¼ xe . The eigenvalues of J e are the roots of l2 trðJ e Þl þ detðJ e Þ ¼ 0,
(25)
which, according to Descartes’ rule, have a negative real part if and only if trðJ e Þo0;
detðJ e Þ40.
(26)
Taking into account Eqs. (24) and (26), the following proposition can be stated. Proposition 3.1. A positive equilibrium point ðxe ; ye Þ is asymptotically stable if and only if dj fye o ð1 þ dxe Þ; dx cxe
dj dg g o ¼ . dx dx f ð1 þ dxe Þ2
(27)
From a geometric point of view, the second of inequalities (27) means that the slope of the tangent to jðxÞ at x ¼ xe must be less than the slope of the tangent to gðxÞ, so that, in the neighbourhood of xe , jðxÞ4gðxÞ for xoxe and jðxÞogðxÞ for x4xe . Assuming c, f, g and d positive, at a positive equilibrium point the right-hand sides of inequalities (27) are both positive. Therefore we can state the following sufficient condition for asymptotic stability. Proposition 3.2. A positive equilibrium point ðxe ; ye Þ is asymptotically stable if dj 1 h ^ b 2bdxe o0. ¼ þ ad dx c x2e
(28)
Example 3.1 (Jing and Yang [10]). Consider the system (5)–(6) with a^ ¼ 0:9, b ¼ 0:1, c ¼ 2, d ¼ 0:5, e^ ¼ 1, f ¼ 1, g ¼ 0:4, h ¼ 0:8. The (unique) positive equilibrium point ðxe ; ye Þ ¼ ð4:93861; 0:42352Þ is asymptotically stable because dj=dx ¼ 0:056 at x ¼ xe . When h ¼ 0, the condition of Proposition 3.2 is satisfied if xe lies to the right of the vertex of parabola (9). Therefore, taking into account (7) and (10), the following proposition holds. Proposition 3.3. For h ¼ 0, a positive equilibrium point ðxe ; ye Þ is asymptotically stable if xe 4
ða kÞd b . 2bd
(29)
ARTICLE IN PRESS 496
W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
The locus described by the vertices ðxv ; yv Þ of parabolas (9) as k varies from 0 to the positive value ða bÞ=d is represented by a dashed line in Fig. 1. If the equilibrium point is located to the right of this locus, then it will necessarily be asymptotically stable. In the case of Fig. 1 all of the equilibrium points internal to R0 turn out to be asymptotically stable. Concerning the number of strictly positive equilibrium points corresponding to a given set of model parameters in Eqs. (5) and (6), let us consider the equation obtained by equating Eqs. (13) and (19) : FðxÞ¼jðxÞ gðxÞ ¼ 0, (30) whose roots comprise all the strictly positive equilibrium points (as well as other equilibrium points, if any, whose coordinates are not strictly positive). Function FðxÞ is rational with numerator given by : N FðxÞ ¼f ð1 þ dxÞ2 ðax bx2 hÞ cx½ðg edÞx e (31) and denominator by : DF ðxÞ¼cfxð1 þ dxÞ,
(32)
which is a second-degree polynomial that vanishes for x ¼ 0 and x ¼ 1=do0. Since the degree of N FðxÞ is four, the number of its real zeros (not necessarily distinct or corresponding to positive values of both x and y) may be 4, 2 or 0. If, however, h ¼ 0, then x can be factored out of the numerator thus canceling the corresponding factor in the denominator and the number of positive zeros of FðxÞ may be either 3 or 1 (by counting multiplicities). In the latter case, if the horizontal asymptote of hyperbola (13) (see the second of Eq. (14)) is sufficiently low, the unique intersection occurs at a point where the slope of the tangent to parabola (9) is negative (stable equilibrium); if, instead, the horizontal asymptote is sufficiently high, the intersection occurs for negative values of x and no positive equilibrium exists; for intermediate values of the second of Eq. (14) the intersection occurs when x is positive but less than xv , and the stability of the equilibrium can be checked using the first of conditions (27) only. Three intersections of Eq. (9) with Eq. (13) are also possible but rare in practice. Example 3.2. Consider the system (5)–(6) with a^ ¼ 1:8535, b ¼ 0:25, c ¼ 2:2929, d ¼ 1:4142, e^ ¼ 1, f ¼ 2, g ¼ 10:5666, h ¼ 0. This system exhibits three positive equilibrium points, as shown in Fig. 3. According to Proposition 3.1, the equilibria at ð1; 1:6884Þ and ð3; 2:5233Þ are asymptotically stable, whereas the equilibrium at ð2; 2:2600Þ is unstable. If an equilibrium point, chosen according to convenience or profit criteria, turns out to be unstable, resort can be made to standard feedback-stabilization techniques, as shown in Fig. 4. To this purpose, the preliminary determination of a linearized model of the process may be necessary. Of course, it is important to evaluate the domain of attraction of the stabilized equilibrium point as well as the feasibility of the resulting stabilization policy: for instance, this policy might require the insertion of preys and/or predators, which is often forbidden. If more or all equilibrium points must be stabilized simultaneously, one can exploit the method suggested in [13] since the number of inputs of the considered system equals its state dimension: in this case, the structure of the (nonlinear) stabilizing controller does not change with the equilibrium point.
ARTICLE IN PRESS W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
497
3
2.5
(3,2.5233) (2,2.2600)
2 (1,1.6884) 1.5
1
0.5
0
0
0.5
1
1.5
2
2.5
3
3.5
4
Fig. 3. Predator–prey model with three positive equilibria: the dashed line represents y ¼ jðxÞ and the solid line represents y ¼ gðxÞ.
Fig. 4. Feedback stabilization scheme.
4. Conclusions The convex state-plane region where positive equilibria of predator–prey models of form (1) and (2) can be assigned through an affine state-feedback control law like Eqs. (3) and
ARTICLE IN PRESS 498
W. Krajewski, U. Viaro / Journal of the Franklin Institute 345 (2008) 489–498
(4), has been determined. The values of the control parameters corresponding to a given equilibrium state can easily be found (see Eqs. (17), (18) and (20)). Geometrical considerations have allowed us to derive simple conditions for the stability of a positive equilibrium state (see Eqs. (27) and (29)). The adopted approach integrates recent alternative results [10] regarding the impact of state-feedback control on predator–prey systems with functional response and, in these authors’ opinion, provides a more practical design tool. References [1] H.I. Freedman, Deterministic Mathematical Models in Population Ecology, Marcel Dekker, New York, 1980. [2] A. Cornish-Bowden, Fundamentals of Enzyme Kinetics, Portland, London, 1995. [3] G.F. Fussmann, S.P. Ellner, K.W. Shertzer, N.G. Hairston Jr., Crossing the Hopf bifurcation in a live predator–prey system, Science 290 (2000) 1358–1360. [4] M.L. Rosenzweig, R.H. MacArthur, Graphical representation and stability conditions of predator–prey interactions, Am. Naturalist 97 (1963) 209–223. [5] E.A. McGehee, E. Peacock–Lopez, Turing patterns in a modified Lotka–Volterra model, Phys. Lett. A 342 (2005) 90–98. [6] A.J. Lotka, Elements of Physical Biology, Williams and Wilkins, Baltimore, MD, 1925. [7] V. Volterra, Lecons sur la Theorie Methematique de la Lutte pour la Vie, Gauthier-Vilars, Paris, 1931. [8] F.M. Scudo, J.R. Ziegler, The Golden Age of Theoretical Ecology, Lecture Notes in Biomathematics, vol. 22, Springer, New York, 1978, pp. 1923–1940. [9] C.S. Holling, Some characteristics of simple types of predation and parasitism, Can. Entomol. 91 (1959) 385–398. [10] H. Jing, Z. Yang, The impact of state feedback control on a predator–prey model with functional response, Discrete Continuous Dynam. Syst.—Ser. B 4 (2004) 607–614. [11] M. Fan, Y. Kuang, Dynamics of a nonautonomous predator–prey system with the Beddington-De Angelis functional response, J. Math. Anal. Appl. 295 (2004) 15–39. [12] P. De Leenheer, D. Angeli, E.D. Sontag, On predator–prey systems and small gain theorems, Math. Biosci. Eng. 2 (2005) 25–42. [13] A. Ferrante, W. Krajewski, A. Lepschy, U. Viaro, Simultaneous stabilization of multiple equilibrium points, in: Proceedings of the 10th IEEE International Conference Methods and Models in Automation and Robotics (MMAR’04), Miedzyzdroje, Poland, August 2004, pp. 321–324.