Microgrid (MG) has extensive properties to overcome common problems of local distribution system. Some of those problems are generation and demand difference, blackout and brownout, environmental concerns due to burning of natural resources in power stations (indirectly), and reliability issues. Research on microgrid is being conducted to enhance its features to mitigate power quality (PQ) problems associated with distribution system. Voltage sag and swell have been major power quality problems for decades, loads in distribution system are heavily affected due to these power quality problems. In the distribution system, microgrid and power quality compensation strategy should be existed in order to ensure reliability and voltage sag/swell mitigation. Dynamic voltage restorer (DVR) is comprehensive power electronics based Flexible Alternating Current Transmission System (FACTS) device, it is third-generation FACTS device as its control scheme selection flexibility and power line coupling approach make it advance when compare to first and second-generation FACTS devices. In this paper, an Enhanced Dynamic Voltage Restorer (EDVR) is presented to efficiently mitigate voltage sag/swell in grid connected microgrid. On the one side, the presence of microgrid structure ensures reliability of distribution system for local loads on the other side, EDVR ensures voltage sag/swell free power supply for loads. The control strategy of EDVR is based on enhanced synchronous reference frame (ESRF) approach and fuzzy technique system. ESRF is specially design for fast and precise operation of EDVR whereas; fuzzy technique system is responsible for standardized voltage supply for local loads. DC link voltage of EDVR is effectively regulated with the help of proposed control scheme at the time of voltage sag/swell compensation. Stability analysis of ESRF control has been done using modeling of VSC and eigenvalue analysis system. Simulation results on MATLAB/Simulink verified the performance of EDVR under presented control approach hence the specific loads in distribution system are more secure under proposed microgrid system with EDVR.
1.
Introduction
The interaction between two or more species has been an important and interesting issue in biology and ecology since the famous Lotka-Volterra model [1,2] was proposed. The interaction between different species will generate rich interesting dynamics of biological species, and exhibit the complexity and diversity [3,4]. Amensalism, as a typical type of interaction between the species, has been intensively considered in the last decades. Amensalism describes a basic biological interaction in nature, where one species inflicts harm on another not affected by the former, which means that it does not receive any costs or benefits to itself. The first pioneer work for the investigations of amensalism model is due to Sun [5], who, in 2003, proposed the following two-species amensalism model:
where x=x(t) and y=y(t) represent the population densities of two species at time t, respectively; r1 and k1 describe the intrinsic growth rate and the carrying capacity of the first species, respectively; r2 and k2 describe the intrinsic growth rate and the carrying capacity of the second species, respectively; α describes the effect of the second species on the first species. They explored the stability properties of all possible equilibria in system (1.1).
Setting a11=r1k1, a12=r1αk1 and a22=r2k2, we can see that model (1.1) can be rewritten as
where all parameters r1, r2, a11, a12, and a22 are positive real numbers.
Ever since the first amensalism model was presented, the complicated dynamics of amensalism models have been studied extensively (see [6,7,8,9,10,11,12,13,14,15] and the references cited therein). For example, in [6], Guan and Chen considered a two-species amensalism model with Beddington-DeAngelis functional response and gave some comprehensive bifurcation and global dynamics of this system. Recently, Luo et al. [12] proposed an amensalism model with Holling-Ⅱ functional response and weak Allee effect for the first species, and discussed its local dynamics and global structure.
In the real world, from the point of view of human needs, the management of renewable resources, the exploitation of biological resources, and the harvesting of populations are commonly practiced in forestry, fishery, and wildlife management [16,17,18]. Hence, it is necessary to introduce and to consider the harvesting of species in population models. During the last decade, population models with harvesting and the role of harvesting in the management of renewable resources have received significant attention from the researchers [19,20,21,22,23,24]. Generally speaking, there are three types of harvesting: 1) constant harvesting [25]; 2) linear harvesting [25]; 3) non-linear harvesting [26,27]. As is well known, non-linear harvesting is more realistic from biological and economical points of view [28]. In [26], Clark first proposed the non-linear harvesting term h(E,x)=qExcE+lx, which is called Michaelis-Menten type functional form of catch rate, here q is the catchability coefficient, E is the external effort devoted to harvesting, c and l are constants. After that, many scholars started to focus attention on the influence of the Michaelis-Menten type harvesting on the population systems [28,29,30,31,32,33,34]. For example, a reaction-diffusion predator-prey model with non-local delay and Michaelis-Menten type prey-harvesting was investigated by Zhang et al. [32], they obtained that the discrete and non-local delays are responsible for a stability switch in this system, and a Hopf bifurcation occurs as the delays pass through a critical value. In [33], Hu and Cao discussed a predator-prey system with the nonlinear Michaelis-Menten type predator harvesting and gave a detailed analysis of the stability and bifurcation for this system.
These studies have shown that harvesting activity has important influence to dynamics of population systems and the harvested models can exhibit richer dynamics compared to the model with no harvesting. However, seldom did scholars consider the dynamical behaviors of an amensalism model with a harvesting term. Accordingly, inspired by the previous works, based on the system (1.2), we now introduce the Michaelis-Menten type harvesting to the second species, then we propose the following amensalism model:
where the parameters r1, r2, a11, a12, and a22 are positive constants, and this term qEycE+ly represents Michaelis-Menten harvesting.
Setting ¯t=r1t, ¯x=a11r1x, ¯y=a12r1y, α=a12qElr21, δ=r2r1, β=a22a12, γ=a12cElr1, and dropping the bars, then system (1.3) is transformed into
The layout of this paper is arranged as follows. In the next section, the global dynamics including positivity and boundedness of solutions, number of equilibria, local asymptotical stability, codimension one bifurcations, dynamical behaviors near infinity, closed orbits analysis and global phase portraits are shown for system (1.4). Further, in Section 3, we also obtain the global dynamics of system (1.1) without harvesting term. Numerical simulations and discussions are displayed by Section 4 for demonstrating the theoretical results and the impact of harvesting term. Finally, a brief conclusion is presented in Section 5.
2.
Global dynamics of system (1.4)
In this section, we mainly investigate the global dynamics of system (1.4), which include the positivity and boundedness of solutions, the existence and local stability analysis of equilibria, all possible bifurcation behaviors and the global structure of system (1.4).
2.1. Positivity and boundedness of solutions
Here, we give the positivity and boundedness of the solutions of model (1.4) in the region R2+={(x,y):x≥0,y≥0}.
2.1.1. Positivity
Lemma 2.1. All solutions of model (1.4) with positive initial value are positive for all t≥0.
Proof. Solving model (1.4) with positive initial condition (x(0),y(0)) gives the result:
Therefore, we can easily see that each solution of model (1.4) starting from the positive initial condition (x(0),y(0)) still remains in the first quadrant.
2.1.2. Boundedness
Lemma 2.2. All solutions of system (1.4) with positive initial value are bounded in region Ω={(x(t),y(t)):0≤x(t)<1 and 0≤y(t)<δβ} for the large enough time t.
Proof. If the initial value is selected in R2+, each solution of model (1.4) remains positive from the above result. One can see
from the first equation of model (1.4). And a standard comparison theorem shows that
Also, the following inequality
is obtained from the second equation of system (1.4). It yields
Hence, there exists a sufficiently large time T, such that 0≤x(t)<1 and 0≤y(t)<δβ for t>T. In summary, all the solutions of system (1.4) are always bounded. Thus we complete the proof.
2.2. Existence and stability of equilibria
Theorem 2.1. System (1.4) always has two boundary equilibria E0(0,0) and E1(1,0) for all positive parameters. For the existence of other equilibria, we have
1) For the possible positive equilibria:
i) if α>α∗, then system (1.4) has no positive equilibria;
ii) if α=α∗ and 0<δ−βγ<2β, then there is a unique positive equilibrium E33(1−y∗3,y∗3);
iii) if α∗∗<α<α∗, when α<δ+δγ−β−βγ and δ−βγ>0, or α=δ+δγ−β−βγ and 0<δ−βγ<2β, there is one positive equilibrium E31(1−y∗1,y∗1); when α>δ+δγ−β−βγ and 0<δ−βγ<2β, there are two distinct positive equilibria E31(1−y∗1,y∗1) and E32(1−y∗2,y∗2);
iv) if α=α∗∗ and 0<δ−βγ<β, then E31 coincides with E1 and there is one positive equilibrium E32(1−y∗2,y∗2), where y∗2=δ−βγβ;
v) if 0<α<α∗∗ and α>δ+δγ−β−βγ, then there is one positive equilibrium E32(1−y∗2,y∗2).
2) For the other possible boundary equilibria:
i) if α>α∗, then system (1.4) has no other boundary equilibria; there is one boundary equilibrium E22(0,y∗2).
ii) if α=α∗ and δ−βγ>0, then there is a boundary equilibrium E23(0,y∗3);
iii) if α∗∗<α<α∗ and δ−βγ>0, then there are two distinct boundary equilibria E21(0,y∗1) and E22(0,y∗2);
iv) if α=α∗∗ and δ−βγ>0, then E21 coincides with E0 and there is one boundary equilibrium E22(0,y∗2), where y∗2=δ−βγβ;
v) if 0<α<α∗∗, then
Where α∗=(δ+βγ)24β, α∗∗=δγ, y∗1=δ−βγ−√(δ+βγ)2−4αβ2β, y∗2=δ−βγ+√(δ+βγ)2−4αβ2β, and y∗3=δ−βγ2β.
Proof. It is easy to see that system (1.4) always possesses two boundary equilibria given by E0(0,0) and E1(1,0) for all positive parameters.
1) If x≠0 and y≠0, system (1.4) has a positive equilibrium which satisfies
For the positive equilibria, y must satisfy 0<y<1. Let Δ denote the discriminant of the second equation of (2.1), namely,
and let
Obviously, when 0<α<α∗, the second equation of (2.1) has two roots
When α=α∗, the second equation of (2.1) has one root
Let
and denote by
Then we have
H1) if α>α∗∗, then k(0)>0;
H2) if α=α∗∗, then k(0)=0;
H3) if α<α∗∗, then k(0)<0.
In addition, 0<α∗∗≤α∗, and α∗∗=α∗ if and only if δ=βγ.
Based on the above analysis, combining with H1), H2) and H3), we can derive that:
ⅰ) If α>α∗, then Δ<0 which implies k(y)=0 has no real roots. So system (1.4) has no positive equilibria.
ⅱ) If α=α∗, then Δ=0, so k(y)=0 has a unique real roots y∗3=δ−βγ2β, combing with the condition 0<y∗3<1, we get 0<δ−βγ<2β. In this situation, there exists a unique positive equilibrium E33(1−y∗3,y∗3).
ⅲ) If α∗∗<α<α∗, then Δ>0 and k(0)>0, implying k(y)=0 two positive real roots y∗1 and y∗2 if δ−βγ>0. Moreover, when k(1)<0, namely, α<δ+δγ−β−βγ, we have 0<y∗1<1<y∗2, that is, there is one positive point E31(1−y∗1,y∗1); when k(1)=0 and the symmetry axis of k(y) is y=δ−βγ2β<1, namely, α=δ+δγ−β−βγ and 0<δ−βγ<2β, we have 0<y∗1<y∗2=1, which means there is one positive point E31(1−y∗1,y∗1); when k(1)>0 and the symmetry axis of k(y) is y=δ−βγ2β<1, namely, α>δ+δγ−β−βγ and 0<δ−βγ<2β, we have 0<y∗1<y∗2<1, that is there are two different positive equilibria E31(1−y∗1,y∗1) and E32(1−y∗2,y∗2).
ⅳ) If α=α∗∗, then Δ>0 and k(0)=0, hence k(y)=0 has two real roots y∗1=0 and 0<y∗2=δ−βγβ<1 if 0<δ−βγ<β, thus E31 coincides with E1 and there is one positive equilibrium E32(1−y∗2,y∗2).
ⅴ) If α<α∗∗, then Δ>0 and k(0)<0, so there are y∗1<0 and y∗2>0. When k(1)>0, i.e., α>δ+δγ−β−βγ, we get one positive root 0<y∗2<1, thus there exits one positive equilibrium E32(1−y∗2,y∗2).
The proof of the corresponding boundary equilibria in 2) is similar to that of 1), and hence omitted here. This completes the proof of Theorem 2.1.
In the following, we focus on the local stability of each equilibrium of system (1.4).
Theorem 2.2. For the boundary equilibria E0 and E1, the following statements are true.
1) If 0<α<α∗∗, then E0 is a hyperbolic unstable node and E1 is a hyperbolic saddle.
2) If α>α∗∗, then E0 is a hyperbolic saddle and E1 is a hyperbolic stable node.
3) If α=α∗∗, then
i) when δ≠βγ, E0 and E1 are both saddle-nodes;
ii) when δ=βγ, E0 is a non-hyperbolic saddle and E1 is a non-hyperbolic stable node.
Proof. The Jacobian matrix of system (1.4) evaluated at any equilibrium is
where H(x,y)=δ−2βy−αγ(γ+y)2.
For the equilibrium E0, the Jacobian matrix is
and the two eigenvalues of J(E0) are λ1(E0)=1>0 and λ2(E0)=δ−αγ. Obviously, if 0<α<α∗∗, then λ2(E0)>0, so E0 is a hyperbolic unstable node. If α>α∗∗, then λ2(E0)<0, thus E0 is a hyperbolic saddle. If α=α∗∗, then λ2(E0)=0, this means that the equilibrium E0 is non-hyperbolic, so it is hard to directly judge its type from their eigenvalues. We further discuss its stability properties by applying Theorem 7.1 in Chapter 2 in [35].
In order to change system (1.4) into a standard form, we expand system (1.4) in power series up to the fourth order around the origin
where Q0(y) represents a power series with the terms yi (i≥5). From dxdt=0, we obtain that there is a unique implicit function x=φ0(y)=0 such that φ0(y)+P(φ0(y),y)=0 and φ0(0)=φ′0(0)=0. Then substituting x=φ0(y)=0 into the second equation of (2.7), we get that
Because α=α∗∗=δγ, the coefficient at y2 is δγ−β. By Theorem 7.1 in [35], we have when δ≠βγ, the equilibrium E0 is a saddle-node.
When δ=βγ, (2.8) becomes
Employing the notations of Theorem 7.1 in Chapter 2 in [35], we have m=3 and am=−βγ<0, so the equilibrium E0 is a non-hyperbolic saddle.
For the equilibrium E1, the Jacobian matrix is
The two eigenvalues of the above matrix J(E1) are λ1(E1)=−1<0 and λ2(E1)=δ−αγ. Clearly, if α<α∗∗, then λ2(E1)>0, so E1 is a hyperbolic saddle. If α>α∗∗, then λ2(E1)<0, then E1 is a hyperbolic stable node. If α=α∗∗, then λ2(E1)=0, the equilibrium E1 is non-hyperbolic and its stability cannot be given directly from the eigenvalues. In order to determinate the stability of E1, we translate E1 to the origin by the translation (¯x,¯y)=(x−1,y) and expand system (1.4) in power series up to the fourth order around the origin, which makes system (1.4) to be the following form:
where Q1(¯y) represents a power series with terms ¯yi (i≥5).
To transform the Jacobian matrix into a standard form, we use the invertible translation
then system (2.11) becomes
By introducing a new time variable τ=−t, we get
Based on the implicit function theorem, from dudτ=0, we can deduce a unique function
which satisfies φ1(0)=φ′1(0)=0 and φ1(v)+P(φ1(v),v)=0. Then substituting it into the second equation of (2.14), we have
From α=δγ it follows that the coefficient at v2 is β−δγ. Thus, by using Theorem 7.1 in Chapter 2 in [35], we get that when δ≠βγ, the equilibrium E1 is a saddle node.
When δ=βγ, (2.15) can be written as
By Theorem 7.1 in [35] again, we obtain m=3 and am=βγ>0. Then E0 is a non-hyperbolic unstable node. Hence, E1 is a non-hyperbolic stable node due to that we have used the transformation τ=−t. Accordingly, Theorem 2.2 is proved.
Theorem 2.3. For the boundary equilibria E21, E22 and E23, the following statements are true.
1) Assume α=α∗ and δ−βγ>0, then there exists the boundary equilibrium E23(0,y∗3). Moreover,
i) If δ−βγ>2β, then E23 is a saddle-node, which includes a stable parabolic sector.
ii) If 0<δ−βγ<2β, then E23 is a saddle-node, which includes an unstable parabolic sector.
iii) If δ−βγ=2β, then E23 is non-hyperbolic.
2) Assume α∗∗<α<α∗ and δ−βγ>0, then there are two boundary equilibria E21(0,y∗1) and E22(0,y∗2). Moreover,
i) If α<δ+δγ−β−βγ, then E21 is a hyperbolic unstable node and E22 is a hyperbolic stable node.
ii) If α>δ+δγ−β−βγ and δ−βγ>2β, then E21 is a hyperbolic saddle and E22 is a hyperbolic stable node.
iii) If α>δ+δγ−β−βγ and 0<δ−βγ<2β, then E21 is a hyperbolic unstable node and E22 is a hyperbolic saddle.
iv) If α=δ+δγ−β−βγ and δ−βγ>2β, then E21 a saddle-node, which includes an unstable parabolic sector and E22 is a hyperbolic stable node.
v) If α=δ+δγ−β−βγ and 0<δ−βγ<2β, then E21 a hyperbolic unstable node and E22 is a saddle-node, which includes a stable parabolic sector.
3) Assume α=α∗∗ and δ−βγ>0, then there is one boundary equilibrium E22(0,y∗2), where y∗2=δ−βγβ. Moreover,
i) If 0<δ−βγ<β, then E22 is a hyperbolic saddle.
ii) If δ−βγ=β, then E22 is a saddle-node, which includes a stable parabolic sector.
iii) If δ−βγ>β, then E22 is a hyperbolic stable node.
4) Assume 0<α<α∗∗, then there is one boundary equilibrium E22(0,y∗2).
i) If α<δ+δγ−β−βγ, then E22 is a hyperbolic stable node.
ii) If α=δ+δγ−β−βγ, then E22 is a saddle-node, which includes a stable parabolic sector.
iii) If α>δ+δγ−β−βγ, then E22 is a hyperbolic saddle.
Proof. For the equilibrium E2i(i = 1, 2, 3), the Jacobian matrix is
where H(x∗i,y∗i)=δ−2βy∗i−αγ(γ+y∗i)2. The two eigenvalues of the above matrix J(E2i) are λ1(E2i)=1−y∗i and λ2(E2i)=H(x∗i,y∗i).
Because (x∗i,y∗i) satisfies δ−βy∗i−αγ+y∗i=0, we can derive
Then
and
In addition, for λ1(E2i)=1−y∗i, we discuss it in the following four cases.
Case 1. α=α∗ and δ−βγ>0.
If δ−βγ>2β, then y∗3>1, which implies λ1(E23)=1−y∗1<0, so E23 is non-hyperbolic. In order to determinate the stability of the equilibrium E23, we translate E23 to the origin by letting (¯x,¯y)=(x,y−y∗3), and expand the system in power series up to the third order around the origin, under which system (1.4) can be transformed into
where c0=δ−2βy∗3−αγ(γ+y∗3)2=0, c1=αγ(γ+y∗3)3−β, c2=−αγ(γ+y∗3)4, and Q2(y) represents for a power series with terms ¯yi (i≥4).
Then introducing a new time variable τ=(1−y∗3)t, we get
We see that the coefficient at ¯y2 is c11−y∗3>0. Hence, from Theorem 7.1 in Chapter 2 of [35], we have E23 is a saddle-node, which includes a stable parabolic sector and this parabolic sector is on the upper half-plane.
If δ−βγ<2β, then y∗3<1, which means λ1(E23)=1−y∗1>0. Same analysis as the above we can get E23 is also a saddle-node, which includes an unstable parabolic sector and the parabolic sector is on the lower half-plane.
If δ−βγ=2β, then y∗3=1, which means λ1(E23)=1−y∗1=0, so E23 is a non-hyperbolic critical point with two zero eigenvalues.
Case 2. α∗∗<α<α∗ and δ−βγ>0.
1) If α<δ+δγ−β−βγ, then y∗1<1<y∗2, that implies λ1(E21)=1−y∗1>0 and λ1(E22)=1−y∗2<0, so E21 is a hyperbolic unstable node and E22 is a hyperbolic stable node;
2) If α>δ+δγ−β−βγ and δ−βγ>2β, then 1<y∗1<y∗2, implying λ1(E21)=1−y∗1<0 and λ1(E22)=1−y∗2<0, thus E21 is hyperbolic saddle and E22 is hyperbolic stable node;
3) If α>δ+δγ−β−βγ and δ−βγ<2β, then y∗1<y∗2<1, that means λ1(E21)=1−y∗1>0, λ1(E22)=1−y∗2>0, hence E21 is a hyperbolic unstable node and E22 is a hyperbolic saddle;
4) If α=δ+δγ−β−βγ and δ−βγ>2β, then 1=y∗1<y∗2, which means λ1(E21)=1−y∗1=0 and λ1(E22)=1−y∗2<0, so E21 is non-hyperbolic and E22 is a hyperbolic stable node.
In order to determinate the stability of the equilibrium E21, we translate the equilibrium E21 to the origin by the translation (¯x,¯y)=(x,y−1), and expand system (1.4) in power series up to the third order around the origin, which makes system (1.4) to be the following form:
where d0=δ−2β−αγ(γ+1)2, d1=αγ(γ+1)3−β, d2=−αγ(γ+1)4, and Q3(¯y) represents for a power series with terms ¯yi satisfying i≥4.
Then introducing a new time variable τ=d0t, we get
From d¯ydτ=0, we can derive a unique implicit function ¯y=ϕ(¯x)=0, which satisfies ϕ(0)=ϕ′(0)=0 and ϕ(¯x)+P(¯x,ϕ(¯x))=0. Substituting ¯y=ϕ(¯x)=0 into the first equation of (2.20), we have that
The coefficient at ¯x2 is
By Theorem 7.1 in Chapter 2 of [35], we know E21 is a saddle-node, which includes an unstable parabolic sector and this parabolic sector is on the left half-plane.
5) If α=δ+δγ−β−βγ and δ−βγ<2β, then y∗1<y∗2=1, implying λ1(E21)=1−y∗1>0 and λ1(E22)=1−y∗2=0, therefore E21 a hyperbolic unstable node and E22 is non-hyperbolic. Similarly to the proof of E21 in 4), we can get that E22 is a saddle-node, which includes a stable parabolic sector.
Case 3. α=α∗∗ and δ−βγ>0.
If δ−βγ>β, then y∗2>1, that is λ1(E22)=1−y∗2<0, so E22 is a hyperbolic stable node. If δ−βγ<β, then y∗2<1, that is λ1(E22)=1−y∗2>0, thus E22 is a hyperbolic saddle. If δ−βγ=β, then y∗2=1, that is λ1(E22)=1−y∗2=0, thus E22 is a non-hyperbolic. Similar to the proof of E21 in 4), we can obtain E22 is a saddle-node, which includes a stable parabolic sector.
Case 4. α<α∗∗.
If α>δ+δγ−β−βγ, then y∗2<1, which implies λ1(E22)=1−y∗2>0, thus E22 is a hyperbolic saddle. If α<δ+δγ−β−βγ, then y∗2>1, which means λ1(E22)=1−y∗2<0, thus E22 is a hyperbolic stable node. If α=δ+δγ−β−βγ, then y∗2=1, which implies λ1(E22)=1−y∗2=0, thus E22 non-hyperbolic. Similar to the proof of 4) in Case 2, we can derive E22 is a saddle-node, which includes a stable parabolic sector.
This completes the proof of Theorem 2.3.
Theorem 2.4. For the positive equilibria E31, E32 and E33, the following statements are true.
1) If α∗∗<α<α∗, when α<δ+δγ−β−βγ and δ−βγ>0, or α=δ+δγ−β−βγ and 0<δ−βγ<2β, or α>δ+δγ−β−βγ and 0<δ−βγ<2β, then E31 is a hyperbolic saddle.
2) If 0<α<α∗∗ and α>δ+δγ−β−βγ, or α=α∗∗ and 0<δ−βγ<β, or α∗∗<α<α∗, α>δ+δγ−β−βγ and 0<δ−βγ<2β, then E32 is a hyperbolic stable node.
3) If α=α∗ and 0<δ−βγ<2β, then E33 is a saddle-node.
Proof. For the equilibrium E3i(i = 1, 2, 3), the Jacobian matrix is
where H(x∗i,y∗i)=δ−2βy∗i−αγ(γ+y∗i)2. The two eigenvalues of the above matrix J(E3i) are λ1(E3i)=y∗i−1<0 and λ2(E3i)=H(x∗i,y∗i). From the analysis of the proof of Theorem 2.3, it follows that λ2(E31)>0, λ2(E32)<0, and λ2(E33)=0.
Hence, when α∗∗<α<α∗, α<δ+δγ−β−βγ and δ−βγ>0, or α∗∗<α<α∗, α=δ+δγ−β−βγ and 0<δ−βγ<2β, or α∗∗<α<α∗, α>δ+δγ−β−βγ and 0<δ−βγ<2β, we have λ1(E31)<0 and λ2(E31)>0, which means E31 is a hyperbolic saddle.
When 0<α<α∗∗ and α>δ+δγ−β−βγ, or α=α∗∗ and 0<δ−βγ<β, or α∗∗<α<α∗, α>δ+δγ−β−βγ and 0<δ−βγ<2β, there exists the positive equilibrium E32 with two eigenvalues λ1(E32)<0 and λ2(E32)<0, hence E32 is a hyperbolic stable node.
When α=α∗ and 0<δ−βγ<2β, there exists the positive equilibrium E33 with two eigenvalues λ1(E33)<0 and λ2(E33)=0. Then we further to study the stability of E33 by transforming E33 to the origin by the translation (¯x,¯y)=(x−x∗3,y−y∗3) and expand system (1.4) in power series up to the third order around the origin, which makes system (1.4) to be the following form:
where c0=δy∗3−β(y∗3)2−α+αγγ+y∗3=0, c1=δ−2βy∗3−αγ(γ+y∗3)2=0, and Q4(¯y) stands for a power series with terms ¯yi (i≥4).
In order to transform the Jacobian matrix into a standard form, we use the invertible translation
then system (2.22) can be expressed as
Then introducing a new time variable τ=−x∗3t, we get
According to the implicit function theorem, from dudτ=0, we can derive a unique function
which satisfies φ4(0)=φ′4(0)=0 and φ4(v)+P(φ4(v),v)=0. Substituting it into the second equation of (2.25), we have that
The coefficient at the term v2 is
Therefore, on basis of Theorem 7.1 in [35], the equilibrium E33 of system (1.4) is a saddle-node.
The proof of Theorem 2.4 is finished.
For readers' convenience, the local dynamical properties of equilibria of system (1.4) are totally summarized in Table 1. It can be clearly seen that the number and stability of equilibria are different according to the value range of the parameters. Therefore, we divide (α,β,γ,δ)∈R4+ into the following twelve regions.
For R1=R11∪R12, system (1.4) possesses two equilibria E0 and E1. In R11, E0 is a saddle and E1 is a stable node. In R12, E0 and E1 represent saddle-nodes. And specifically, in R1, E1 is stable in the first quadrant. Therefore, the qualitative properties of these two equilibria can be seen in Figure 1(a).
For R2=R21∪R22, there are three equilibria for system (1.4). In R21=R2a∪R2b, system (1.4) admits E0 (unstable node or saddle-node), E1 (saddle or saddle-node) and E22 (stable node or saddle-node). In particular, in R21, E0 is unstable and E22 is stable in the first quadrant. In R22, there are E0 (saddle), E1 (stable node) and E23 (non-hyperbolic or saddle-node). The qualitative properties of these equilibria in R2 are displayed by Figure 1(b), (c).
For R3=R31∪R32∪R33, there exist four equilibria for system (1.4). Considering R31=R3a∪R3b, there are E0 (unstable node or saddle-node which is unstable in the first quadrant), E1 (saddle or saddle-node), E22 (saddle) and E32 (stable node). Considering R32, system (1.4) admits E0 (saddle), E1 (stable node), E21 (saddle or saddle-node) and E22 (stable node). Considering R33, there exist E0, E1, E23 and E33, which are saddle, stable node, saddle-node and saddle-node. The detailed image descriptions of the qualitative properties of the equilibria in R3 are illustrated by Figure 1(d)–(f).
For R4=R41∪R42, system (1.4) has five equilibria E0 (saddle), E1 (stable node), E21 (unstable node), E22 (stable node or saddle-node) and E31 (saddle). Particularly, E22 is stable in the first quadrant. As shown in Figure 1(g), the qualitative properties of these five equilibria are given.
For R5, system (1.4) possesses six equilibria E0, E1, E21, E22, E31 and E32, which represent saddle, stable node, unstable node, saddle, saddle and stable node, respectively. Hence, the qualitative properties of these six equilibria are displayed by Figure 1(h).
2.3. Bifurcation analysis
In this section, all possible bifurcations of system (1.4) will be found out. Moreover, the feasible conditions for the occurence of those bifurcations will be given.
Theorem 2.5. System (1.4) undergoes a saddle-node bifurcation at equilibrium E33 when the parameters satisfy the restriction α=αSN=(δ+βγ)24β along with the conditions 0<δ−βγ<2β and α>δ+δγ−β−βγ as given in Theorem 2.1.
Proof. Now we will verify the transversality condition for the occurrence of a saddle-node bifurcation at α=αSN by utilizing the Sotomayor's theorem [36]. From Section 2.2, we have the two eigenvalues of J(E33) are λ1(E33)<0 and λ2(E33)=0. Denote V and W the eigenvectors corresponding to the eigenvalue λ2(E33) for the matrices J(E33) and J(E33)T, respectively. Simple computations yield
Furthermore, we can obtain
and
Obviously, when 0<δ−βγ<2β, the vectors V and W satisfy the transversality conditions
and
Hence, by Sotomayor's theorem, when α=αSN, system (1.4) undergoes a saddle-node bifurcation at non-hyperbolic critical point E33. The number of positive equilibria of system (1.4) changes from zero to two as α passes from the right of α=αSN to the left. Thus, the proof of Theorem 2.5 is completed.
Theorem 2.6. System (1.4) undergoes a transcritical bifurcation at the boundary equilibrium E1 when the parameters satisfy the conditions α=αTC=δγ and 0<δ−βγ<2β.
Proof. Now we also apply Sotomayor's theorem [36] to prove the transversality condition for the occurrence of transcritical bifurcation at E1 as α=αTC. From Section 2.2, we know that λ1(E1)<0 and λ2(E1)=0. Let V and W represent the two eigenvectors corresponding to the zero eigenvalue λ2(E1) of the matrices J(E1) and J(E1)T, respectively, then they are given by
Furthermore, we can obtain
and
Clearly, the vectors V and W satisfy the transversality conditions
and
Therefore, by Sotomayor's theorem, if δ≠βγ, system (1.4) experiences a transcritical bifurcation at non-hyperbolic critical point E1 as the parameter α varies through the bifurcation value α=αTC.
Thus, the proof of Theorem 2.6 is completed.
2.4. Global structure
In this subsection, we pay some attention to the global structure of the plane nonlinear dynamic system (1.4).
2.4.1. The dynamics near infinity
For the sake of the research of the global dynamics for system (1.4), we have to investigate the qualitative properties of equilibria at infinity, which is important and useful to study the behavior of orbits when |x|+|y| → ∞.
Firstly, using Poincaré transformation of Chapter 5 in [35]:
system (1.4) can be transformed into
In the first quadrant of u-z plane, setting z=0, system (2.30) has one boundary equilibrium e(1β−1,0) if β>1 and no boundary equilibrium if β≤1 on the nonnegative u-axis.
Proceed to the next step, we apply the second type of Poincaré transformation:
system (1.4) becomes
In the nonnegative cone of v-z plane, there exist two boundary equilibria e1(0,0) and e2(β−1,0) if β>1 and only one e1(0,0) if β≤1 for system (2.31).
In this paper, of concern is the global dynamics of system (1.4) when β>1. The case β≤1 is similar, hence, it is omitted.
Theorem 2.7. In system (2.30), e(1β−1,0) is a saddle. In system (2.31), e1(0,0) is an unstable node and e2(β−1,0) is a saddle.
Proof. The Jacobian matrix at e(1β−1,0) is evaluated as follows:
obviously, which has two eigenvalues λ1(e)=−1<0 and λ2(e)=ββ−1>0. Hence, e is a saddle in system (2.30).
Similarly, the corresponding Jacobian matrices of system (2.31) evaluated at e1(0,0) and e2(β−1,0) are
Clearly, e1 is an unstable node and e2 is a saddle in system (2.31). This proof of Theorem 2.7 is completed.
From the point of view of Poincaré transformation, we get
Consequently, e and e2 are equivalent to an infinite singular point I for system (1.4), e1 is equivalent to an infinite singular point I1 for system (1.4). By combining the Poincaré transformation of Chapter 5 in [35] and Theorem 2.7, we illustrate the dynamical behaviors near infinity of model (1.4) with Figure 2.
2.4.2. Nonexistence of closed orbit
We, here, for further research of the global structure of system (1.4), need to discuss the existence of closed orbits in phase space.
Theorem 2.8. No closed orbit exists for model (1.4) in R2+.
Proof. On the contrary, assume that model (1.4) exhibits a closed orbit in R2+. From the discussion in Theorem 2.1, we can deduce that all positive equilibria will lie on invariant lines y=δ−βγ−√Δ2β or y=δ−βγ+√Δ2β. On the basis of Theorem 4.6 of [35], in the interior of the closed orbit, system (1.4) must admit a positive equilibrium. As such, the closed orbit has two points of intersection with y=δ−βγ−√Δ2β or y=δ−βγ+√Δ2β, which is contrary to the existence and uniqueness of solutions. Therefore, there is no closed orbit for model (1.4). We end the proof of Theorem 2.8.
2.4.3. Global phase portraits
System (1.4), according to the aforementioned Theorem 2.7, admits two equilibria I and I1 at infinity, which stand for a saddle and an unstable node, respectively. A step further, we derive that model (1.4) exists no closed orbit in R2+ from the above Theorem 2.8. It now follows that all possible cases of global phase portraits for model (1.4) can be given.
For the purpose of simplicity, ξI and ω(ξI) are utilized to refer the unstable manifold of I and the ω-limit set of ξI in the nonnegative cone of R2, respectively.
For (α,β,γ,δ) ∈ R1, it is easy to be checked that ω(ξI) is E1. See Figure 3(a) for the global phase portrait of model (1.4) related to the case R1.
For (α,β,γ,δ) ∈ R2, we can easily verify that ω(ξI) is E22 in R21 and ω(ξI) is E23 in R22. See Figure 3(b), (c) for the global phase portraits of model (1.4) related to the cases R21 and R22.
For (α,β,γ,δ) ∈ R3, we analyze the global phase portraits of model (1.4) from the three subregions R31, R32 and R33, respectively. When (α,β,γ,δ) ∈ R31, we can verify that ω(ξI) is E32. The global phase portraits of model (1.4) related to R31 are given by Figure 3(d), (e). When (α,β,γ,δ) ∈ R32, ω(ξI) is E22. The global phase portrait of model (1.4) related to R32 is shown by Figure 3(f). When (α,β,γ,δ) ∈ R33, it is clear that ω(ξI) is E33 and the global phase portraits of model (1.4) in this region are shown in Figure 3(g), (h).
For (α,β,γ,δ) ∈ R4, it is clear that ω(ξI) is E22. See Figure 3(i) for the global phase portrait of model (1.4) related to the case R4.
For (α,β,γ,δ) ∈ R5, it is also clear that ω(ξI) is E32. And the global phase portraits of model (1.4) in R5 are shown in Figure 3(j), (k).
3.
Global dynamics of system (1.1)
Now, we mainly consider the dynamic properties of all possible equilibria and bifurcation behaviors of model (1.1). More, in R2+, we analyze the global structure of this model and illustrate our analysis with the global phase portraits.
To facilitate this, system (1.1) takes the following form
On account of the dynamic properties of system (3.1) corresponding to system (1.1), we study the dynamics of system (3.1) in the following.
3.1. Existence and stability of equilibria
As it has been shown in [5], the existence and stability properties of equilibria of system (3.1) are clearly summarized in Table 2. Next, we divide (δ,β)∈R2+ into two regions as follows.
Consider (δ,β)∈G1, system (3.1) has three equilibria E0, E1 and E2. Specifically, E0 and E1 represent an unstable node and a saddle. E2 is a stable node in G11 and non-hyperbolic in G12, which is stable in the first quadrant. See Figure 4(a), the qualitative properties of these three equilibria are illustrated.
Consider (δ,β)∈G2, system (3.1) admits four equilibria E0, E1, E2 and E3, which are unstable node, saddle, saddle and stable node, respectively. See Figure 4(b), the qualitative properties of these four equilibria are illustrated.
3.2. Bifurcation analysis
Here, we will analyze some of all the possible bifurcation scenarios for model (3.1). According to the previous methods applied in Section 2.3, a heuristic discussion similar to the one used by Theorem 2.6 can easily get the conditions for a transcritical bifurcation. Therefore, the following theorem can be directly obtained.
Theorem 3.1 System (3.1) undergoes a transcritical bifurcation at E2 when the parameters satisfy the condition δ=β.
3.3. Global structure
In order to investigate the global structure of system (3.1), the methods proposed by Section 2.4 can be extended. By applying Poincaré transformation and analysis of closed orbit, Figure 5(a) describes the existence and stability of the singular points at infinity, as well as all the boundary equilibria. In addition, Figure 5(b)–(d) describes minutely the global phase portraits in R2+ similarly.
4.
Numerical simulations and discussions
Some numerical examples, in this section, are proposed so as to check the obtained parameter conditions and theoretical results of systems (1.4) and (3.1). Further, we are interested to show the role of harvesting on the local dynamical properties.
Example 4.1. We consider parameter values as below: α = 2, β = 1.2, γ = 0.8, δ = 2.5. Simple calculations lead to α=α∗∗, δ−βγ≥β and δ>β, which means both systems (1.4) and (3.1) have three boundary equilibria. The numerical phase portraits related to this case are shown in Figures 6(a) and 7(a). The impact of harvesting is shown by Figure 8.
Example 4.2. Some parameter values are taken α = 1.85, β = 0.3, γ = 2.5 and δ = 0.8. In this case, parameters satisfy the conditions 0<α<α∗∗ and α>δ+δγ−β−βγ. Consequently, system (1.4) has a positive equilibrium (global asymptotically stable) and three boundary equilibria. From δ>β, we can effortlessly observe there exists no positive equilibrium for model (3.1). See Figures 6(b) and 7(b), the numerical phase portraits corresponding to this situation are presented. The impact of harvesting is illustrated with curves in Figure 9.
Example 4.3. We select parameter values as follows: α = 0.8, β = 2.5, γ = 0.4, δ = 2. For these parameters, we derive α=α∗∗ and 0<δ−βγ<2β. Therefore, system (1.4) possesses the same types of equilibria as the above Example 4.2. Whereas, there is a positive equilibrium (global asymptotically stable) in system (3.1) for δ<β. See Figures 6(c) and 7(c), the corresponding numerical phase portraits are presented. The impact of harvesting is illustrated with curves in Figure 10.
It can be seen that the harvesting effect has a significant role on the population densities. From a biological point of view, the permanence and extinction of species are influenced to differential extent by the harvesting effect on the second species, for instance, greatly slowing the extinction of the first species (Figure 8) or guaranteeing the permanence of the system by preserving the first species from extinction (Figure 9). In addition, both species need much more time to reach their stable coexistence state due to the harvesting effect, which is shown in Figure 10.
5.
Conclusions
In this article, by combining local and global dynamical analysis, we have investigated the dynamics of the amensalism system with Michaelis-Menten type harvesting for the second species. In order to further explore the influence of harvesting effect, a complete qualitative analysis of the system with no harvesting is also performed. We could show some differences for the model (1.4) with harvesting on the second species and the system (3.1) with no harvesting.
1) When harvesting is present in the system, it has shown that model (1.4) has more different types of equilibria. System (3.1) has at most four equilibria including one interior point in R2+ and the unique interior point (if it exists) of system (3.1) is globally stable. However, when the Michaelis-Menten type harvesting is on the second species, we shown that the model (1.4) has at most six equilibria including two interior points in R2+ and the interior point E31 is a saddle.
2) The complex existence of equilibria leads to more bifurcation behaviors of model (1.4). Two kinds of bifurcation, under some parameter restrictions, saddle-node bifurcation and transcritical bifurcation will occur for system (1.4), while there exhibits only one saddle-node bifurcation for the system with no harvesting.
3) More types of equilibria and more bifurcation behaviors can induce potentially dramatic changes to the dynamics of the system. By the classifications of global phase portraits of model (1.4), we find that more complex dynamics occur for the appearance of harvesting on the second species. Such as, for model (1.4) with harvesting, the first species survives permanently rather than extinction or both species take a lot longer to approach the coexistence state.
These results reveal that the dynamical properties of system (1.4) get more complicated and richer compared to the system with no harvesting. Furthermore, it would also be interesting to study system (3.1) with both the first species and the second species with harvesting terms since we usually harvest, or would like to harvest, both populations.
Acknowledgments
The first author was supported by the National Natural Science Foundation of China (No.12001503) and the third author was supported by the Project of Beijing Municipal Commission of Education (KM 202110015001).
Conflict of interest
The authors declare there is no conflicts of interest.