Loading [MathJax]/jax/element/mml/optable/GeneralPunctuation.js
Research article

Extended DEA method for solving multi-objective transportation problem with Fermatean fuzzy sets

  • Received: 26 July 2022 Revised: 12 September 2022 Accepted: 15 September 2022 Published: 13 October 2022
  • MSC : 90C32, 90C70

  • Data envelopment analysis (DEA) is a linear programming approach used to determine the relative efficiencies of multiple decision-making units (DMUs). A transportation problem (TP) is a special type of linear programming problem (LPP) which is used to minimize the total transportation cost or maximize the total transportation profit of transporting a product from multiple sources to multiple destinations. Because of the connection between the multi-objective TP (MOTP) and DEA, DEA-based techniques are more often used to handle practical TPs. The objective of this work is to investigate the TP with Fermatean fuzzy costs in the presence of numerous conflicting objectives. In particular, a Fermatean fuzzy DEA (FFDEA) method is proposed to solve the Fermatean fuzzy MOTP (FFMOTP). In this regard, every arc in FFMOTP is considered a DMU. Additionally, those objective functions that should be maximized will be used to define the outputs of DMUs, while those that should be minimized will be used to define the inputs of DMUs. As a consequence, two different Fermatean fuzzy effciency scores (FFESs) will be obtained for every arc by solving the FFDEA models. Therefore, unique FFESs will be obtained for every arc by finding the mean of these FFESs. Finally, the FFMOTP will be transformed into a single objective Fermatean fuzzy TP (FFTP) that can be solved by applying standard algorithms. A numerical example is illustrated to support the proposed method, and the results obtained by using the proposed method are compared to those of existing techniques. Moreover, the advantages of the proposed method are also discussed.

    Citation: Muhammad Akram, Syed Muhammad Umer Shah, Mohammed M. Ali Al-Shamiri, S. A. Edalatpanah. Extended DEA method for solving multi-objective transportation problem with Fermatean fuzzy sets[J]. AIMS Mathematics, 2023, 8(1): 924-961. doi: 10.3934/math.2023045

    Related Papers:

    [1] Yichao Shao, Hengguo Yu, Chenglei Jin, Jingzhe Fang, Min Zhao . Dynamics analysis of a predator-prey model with Allee effect and harvesting effort. Electronic Research Archive, 2024, 32(10): 5682-5716. doi: 10.3934/era.2024263
    [2] Yujia Xiang, Yuqi Jiao, Xin Wang, Ruizhi Yang . Dynamics of a delayed diffusive predator-prey model with Allee effect and nonlocal competition in prey and hunting cooperation in predator. Electronic Research Archive, 2023, 31(4): 2120-2138. doi: 10.3934/era.2023109
    [3] Jianzhi Cao, Xuhong Ji, Pengmiao Hao, Peiguang Wang . Bifurcation analysis of an SIS model with a modified nonlinear incidence rate. Electronic Research Archive, 2025, 33(6): 3901-3930. doi: 10.3934/era.2025173
    [4] Zhuo Ba, Xianyi Li . Period-doubling bifurcation and Neimark-Sacker bifurcation of a discrete predator-prey model with Allee effect and cannibalism. Electronic Research Archive, 2023, 31(3): 1405-1438. doi: 10.3934/era.2023072
    [5] Xin Du, Quansheng Liu, Yuanhong Bi . Bifurcation analysis of a two–dimensional p53 gene regulatory network without and with time delay. Electronic Research Archive, 2024, 32(1): 293-316. doi: 10.3934/era.2024014
    [6] Fengrong Zhang, Ruining Chen . Spatiotemporal patterns of a delayed diffusive prey-predator model with prey-taxis. Electronic Research Archive, 2024, 32(7): 4723-4740. doi: 10.3934/era.2024215
    [7] Jiange Dong, Xianyi Li . Bifurcation of a discrete predator-prey model with increasing functional response and constant-yield prey harvesting. Electronic Research Archive, 2022, 30(10): 3930-3948. doi: 10.3934/era.2022200
    [8] Mengting Sui, Yanfei Du . Bifurcations, stability switches and chaos in a diffusive predator-prey model with fear response delay. Electronic Research Archive, 2023, 31(9): 5124-5150. doi: 10.3934/era.2023262
    [9] Mengxin He, Zhong Li . Dynamic behaviors of a Leslie-Gower predator-prey model with Smith growth and constant-yield harvesting. Electronic Research Archive, 2024, 32(11): 6424-6442. doi: 10.3934/era.2024299
    [10] Yanhe Qiao, Hui Cao, Guoming Xu . A double time-delay Holling Ⅱ predation model with weak Allee effect and age-structure. Electronic Research Archive, 2024, 32(3): 1749-1769. doi: 10.3934/era.2024080
  • Data envelopment analysis (DEA) is a linear programming approach used to determine the relative efficiencies of multiple decision-making units (DMUs). A transportation problem (TP) is a special type of linear programming problem (LPP) which is used to minimize the total transportation cost or maximize the total transportation profit of transporting a product from multiple sources to multiple destinations. Because of the connection between the multi-objective TP (MOTP) and DEA, DEA-based techniques are more often used to handle practical TPs. The objective of this work is to investigate the TP with Fermatean fuzzy costs in the presence of numerous conflicting objectives. In particular, a Fermatean fuzzy DEA (FFDEA) method is proposed to solve the Fermatean fuzzy MOTP (FFMOTP). In this regard, every arc in FFMOTP is considered a DMU. Additionally, those objective functions that should be maximized will be used to define the outputs of DMUs, while those that should be minimized will be used to define the inputs of DMUs. As a consequence, two different Fermatean fuzzy effciency scores (FFESs) will be obtained for every arc by solving the FFDEA models. Therefore, unique FFESs will be obtained for every arc by finding the mean of these FFESs. Finally, the FFMOTP will be transformed into a single objective Fermatean fuzzy TP (FFTP) that can be solved by applying standard algorithms. A numerical example is illustrated to support the proposed method, and the results obtained by using the proposed method are compared to those of existing techniques. Moreover, the advantages of the proposed method are also discussed.



    Predator-prey interactions have been an interesting and challenging issue that is frequently discussed in marine ecosystems, especially in fish populations. Predator-prey interactions are the most important component of ecology, determining various factors such as community composition, species behavior and dynamics. Mathematical modeling helps to provide insights into the dynamics of the system, which was investigated in the early pioneering work of Lotka[1] and Volterra[2]. In dynamical systems, continuity, equilibrium stability, bifurcation, and control problems are also often studied[3,4,5,6,7]. Traditionally, predators could be distinguished as specialists or generalists based on whether they ate only one or more types of prey. In general, a predator-prey model can be described by the following ordinary differential equation for both specialists and generalist predators[8]:

    ˙P=F(P)PG(P,N)N,˙N=γυ(G(P,N))N+H(N)N,

    where P and N denote the prey and predator's population densities at moment t, respectively. F() and H() represent the growth of the species in the absence of the other one. G() is known as the functional response, which characterizes the average individual of prey consumed by a predator; γ represents the conversion rate from prey consumption to predator, υ() is a monotonically increasing function. For specialist predators, there is H(N)=d<0; and for generalist predators, it is required that γυ(G(0,N))+H(N)>0. In literature, different type of functional responses were adopted to model different species, which can be prey dependent[9,10] or prey and predator dependent[11,12,13,14,15].

    The classical view of population dynamics claims that the higher the population density, the lower the overall growth rate due to the competition for resources. The lower the density, the higher the overall growth rate. However, when the population density is low, Allee[16] introduced the opposite view that the lower the population density, the lower the overall growth rate, that is, the Allee effect. The Allee effect is a common phenomenon in marine populations[17,18,19]. When the population density is low, it may affect population development due to pairing restrictions, dispersal, habitat changes, cooperative foraging, cooperative defense, and predator saturation. Therefore, the study of the Allee effect on ecosystems has attracted the attention of many scholars. In general, the Allee effect can be represented by a multiplier of the form P-L[20,21,22], where L is the threshold for the Allee effect. When L<0, it is a weak Allee effect, and the Allee effect is always positive no matter how much the prey growth rate decreases. When L>0, it is a strong Allee effect. The strong Allee effect indicates that in order for the population to grow, the population size or density must be higher than L; otherwise, the population will die out. Scholars have analyzed the dynamics of the system with the Allee effect of the form P-L, and discussed the existence of the equilibrium of the system and various bifurcation phenomena such as saddle-node bifurcation and B-T bifurcation[23,24,25].

    Fish is a kind of important ecological resources. In view of the fish resources development problems, scholars studied population behavior by adding harvest items on the basis of continuous systems. However, fishing activities are not continuous, so continuous dynamic systems cannot accurately describe the actual fishing process. In the process of fish harvesting, the periodic harvesting of fish is a kind of human-controlled behavior that can be described by an impulse differential equation, and it has been found that the impulse differential equation is more accurate in describing and portraying the dynamical behavior of the population[26]. The theory of semi-continuous dynamic systems has been widely used in modeling research on pest management modeling[27,28,29,30,31]. Based on the analysis of the literature, due to the fact that impulsive differential systems (semi-continuous dynamical systems) have the characteristics of both continuous and discrete dynamical systems, there are some studies applying the theory to the development of fish populations in deterministic environments[32,33,34,35,36] and uncertain environment[37,38,39]. In addition, most of the studies considered fish harvesting or fish protection only (unilateral control); in this study, a bilateral state control[40,41,42] is considered, and a predator-prey model for conservation and harvesting of two fish species is developed by constructing a semi-continuous dynamical system. When designing the state feedback control strategy, the number of objective fish was used as the state variable for feedback control. On the one hand, when the number of prey fish is small, the Allee effect will lead to their extinction, which will lead to the lack of enough food for prey fish and destroy the balance of the ecosystem. Therefore, it is necessary to release a certain amount of prey fish when the number of prey fish decreases to a certain threshold. When the number of prey fish is high, it is necessary to catch prey fish from the economic point of view. Since the fishing behavior will also lead to the harvesting of some predator fish, in order to maintain the ecological balance and avoid the extinction of predator fish caused by the fishing behavior, it is necessary to release a certain amount of predator fish larvae at the same time. Based on the above two aspects, we propose a bilateral control strategy to maintain the population size of the two species in a suitable range.

    This paper considers a predator-prey model in which the predator is a generalist and the growth of the prey is affected by the Allee effect. The structure is as follows: In Section Ⅱ, we describe the fishery model, non-negativity, persistent survivability and discuss the existence and stability of equilibria of the model. In Section Ⅲ, the bifurcation dynamics of the model are discussed using bifurcation theory. In Section Ⅳ, based on this model, we analyze the model of marine fish harvesting and conservation with bilateral controls. The existence and stability of the order-1 and order-2 periodic solutions of the system are analyzed by using the geometry theory of differential equations. In the fifth section, we performed numerical simulations using MATLAB to verify the correctness of the results.

    In this paper, we present a predator-prey model in which predators are generalists, prey growth rates are logical and subject to strong Allee effects, and the functional response is a Holling-Ⅰ type, while the conversion from prey consumption to predator species is saturated,

    {dPdT=rP(1PK)(PL)APN,dNdT=e(AP1+BP)N+(s1+fNd)N, (2.1)

    where p(T), N(T) denote the prey and predator's densities at the moment of T; K denotes the prey's environmental holding capacity; L is the threshold of the prey's Allee effect; r and s represent the intrinsic growth rates of prey and predator, respectively; A is the capture rate of prey by the predator; e is the efficiency with which prey is converted to predator; B is the half-saturation constant; f is the intensity of predator density dependence, and d is the predator mortality rate. Since the predators are generalist, it requires that eABd, s>d, and all parameters of model (2.1) are positive.

    To facilitate the analysis, let x=PK, y=N, t=rKT, α=ArK, β=LK, γ=eAr, δ=BK, s1=srK, d1=drK. Then system (2.1) is simplified to system (2.2):

    {dxdt=x(1x)(xβ)αxy,dydt=γxy1+δx+(s11+fyd1)y, (2.2)

    and from the biological point of consideration, the model (2.2) is limited in the region

    Ω={(x,y)R2+|0x1,y0}.

    On the other hand, as a renewable resource, fish species are closely related to human life. In order to maintain the balance of the ecosystem during the fishing process, we consider a fish stock control strategy with a combination of fishing and investment. First, to avoid the distinction of prey fish caused by the Allee effect, a quantity (η) of juvenile prey fish is released when the prey density declines to the level x=h1. But at a higher level x=h2, a proportion a of prey fish together with a proportion b of predator fish will also be caught for economic purposes, and simultaneously, a quantity τ of juvenile predator fish is released into the system to maintain the level of predator fish. Based on the control measures, the model can be described as follows:

    {dxdt=x(1x)(xβ)αxydydt=γxy1+δx+(s11+fyd1)y}h1<x<h2,Δx=ηΔy=0}x=h1,Δx=axΔy=by+τ}x=h2, (2.3)

    where η, a, b, τ are all positive, and a,b(0,1).

    For a given planar model

    {dxdt=f1(x,y),dydt=f2(x,y)ω(x,y)0,Δx=I1(x,y),Δy=I2(x,y)ω(x,y)=0, (2.4)

    Definition 1 (Order-k periodic solution[32,33,36]) The solution ˜z(t)=(˜x(t),˜y(t)) is called periodic if there exists m(1) satisfying ˜zm=˜z0. Furthermore, ˜z is an order-k T-periodic solution with kmin{j|1jm,˜zj=˜z0}.

    Lemma 1 (Analogue of Poincaré Criterion[32,33,36]). The order-k T-periodic solution z(t)=(ξ(t),η(t))T is orbitally asymptotically stable if |μq|<1, where

    μk=kj=1Δjexp(T0[f1x+f2y](ξ(t),η(t))dt),

    with

    Δj=f+1(I2yωxI2xωy+ωx)+f+2(I1xωyI1yωx+ωy)f1ωx+f2ωy,

    f+1=f1(ξ(θ+j),η(θ+j)), f+2=f2(ξ(θ+j),η(θ+j)) and f1, f2, I1x, I1y, I2x, I2y, ωx, ωy are calculated at (ξ(θj),η(θj)).

    In this section, the bounded-ness of the solution for Model (2.2) is discussed. Moreover, the existence, type, and local stability of the equilibrium as well as the bifurcation properties are verified.

    Define

    g1(x,y)=(1x)(xβ)αy,f1(x,y)=xg1(x,y);
    g2(x,y)=γx1+δx+(s11+fyd1),f2(x,y)=yg2(x,y).

    Theorem 1. The solution of Model (2.2) with non-negative initial values will remain non-negative for all time and is bounded on R2+.

    Proof. By Eq (2.2), it can be obtained that

    x(t)=x(0)exp[t0g1(x(s),y(s))ds],y(t)=y(0)exp[t0g2(x(s),y(s))ds],

    for all t0, as long as x(0) and y(0) are non-negative, then x(t) and y(t) are also non-negative.

    Next, we define a function u(t)=γαx(t)+y(t). Then

    ˙u=γα˙x+˙y=γα[x(1x)(xβ)αxy]+γxy1+δx+(s11+fyd1)yγα[x(1x)(xβ)αxy]+γxy+(s11+fyd1)yγ(1x)(xβ)αx+s1fd1y[γα((1β)24+d1)+s1f]d1u,

    which implies that

    u(t)1d1[γα((1β)24+d1)+s1f]+(u01d1[γα((1β)24+d1)+s1f])ed1t,

    so long as u0=γαx0+y0 is bounded, u(t) is bounded in Ω. To sum up, any solution of Model (2.2) starting with a non-negative bounded initial condition is non-negative and bounded in Ω.

    Obviously, four equilibria always exist

    O(0,0),E1(0,1f(s1d1d1)),E2(β,0),E3(1,0).

    Define

    y1(x)=1α(1x)(xβ), (3.1)
    y2(x)=1f[s1δ(1δ+x)d1(γd1δ)x1] (3.2)

    and denote ¯γ1d1δ. Due to the assumptions eABd and s>d, then γ<¯γ1, i.e., y2(x)<s1δ¯γ1γ. The positional relationship between y1(x) and y2(x) for different cases is shown in Figure 1.

    Figure 1.  Illustration of the positional relationship between y1(x) and y2(x) for different values of f.

    Let y1(x)=y2(x). Then it has

    a1x3+a2x2+a3x+a4=:g(x)=0, (3.3)

    where

    a1=f(δd1γ)>0,a2=f[(γδd1)(β+1)+d1],a3=α[δ(s1d1)+γ]+f[β(δd1γ)d1(β+1)],a4=α(s1d1)+βfd1>0.

    Define

    ¯αfβ,¯δ1d1(1+β)s1β,¯δ2fd1(1+β)fβd1+α(s1d1),¯δ3(1+β)fd1αs1,
    ¯γ2αδ(s1d1)+fd1(βδ(1+β))fβα,xda223a1a3a23a1,ρg(xd).

    Theorem 2. For any of the following cases: C1) α<¯α, δ<¯δ2, γ<¯γ1; C2) α<¯α, ¯δ2<δ<¯δ3, ¯γ2<γ<¯γ1; C3) α>¯α, δ<¯δ2, 0<γ<¯γ2, if ρ<0 holds, there exists two interior equilibria in Model (2.2); if ρ=0, there exists a unique interior equilibrium in Model (2.2); if ρ>0, there doesn't exists interior equilibrium in Model (2.2).

    Proof. Clearly, the existence of an interior equilibrium is equivalent to that of a positive root of Eq (3.3) in the interval (0, 1). Since

    g(0)=a4>0,g(1)=a1+a2+a3+a4=α[(δ+1)(s1d1)+γ]>0,

    and

    g(x)=3a1x2+2a2x+a3, (3.4)

    then for any of the cases C1)C3), there is a3<0. Moreover, g(0)<0, g(xd)=0, g and

    g'\left( 1 \right) = \left( {1 - \beta } \right)f\left( {\delta {d_1} - \gamma } \right) + \left( {1 - \beta } \right){d_1}f + \alpha \left[ {\delta \left( {{s_1} - {d_1}} \right) + \gamma } \right] > 0,

    so that for x \in (0, {x_d}) , g'(x) < 0 , for x \in({x_d}, 0) , g'(x) > 0 . If \rho < 0 , Eq (3.3) has two distinct positive roots x^*_i\in (0, 1) , i = 1, 2 . Denote y^*_i = y_1(x^*_i) , i = 1, 2 . Then two interior equilibria exist in Model (2.2), denoted by {E^*_{1}}\left({{x^*_{1}}, {y^*_{1}}} \right) and {E^*_{2}}\left({{x^*_{2}}, {y^*_{2}}} \right) . If \rho = 0 , Eq (3.3) has a unique positive root x^* = x_d\in (0, 1) , then a unique positive equilibrium exists in Model (2.2), denoted by {E^*}(x^*, y_1(x^*) ; when \rho > 0 , Eq (3.3) does not have positive root, and thus there does not exist interior equilibrium in Model (2.2).

    For any equilibrium \overline{E}(\overline{x}, \overline{y}) , there is

    J(\overline{E}) = \left[ {\begin{array}{*{20}{c}} {{ - 3{\overline{x}^2} + 2\left( {\beta + 1} \right){\overline{x}} - \beta - \alpha \overline{y}}}&{ - \alpha \overline{x}}\\ {\frac{{\gamma \overline{y}}}{{{{\left( {1 + \delta \overline{x}} \right)}^2}}}}&{\frac{{\gamma \overline{x}}}{{1 + \delta \overline{x}}} + \frac{{{s_1}}}{{{{\left( {1 + f\overline{y}} \right)}^2}}} - {d_1}} \end{array}} \right],

    its characteristic equation is

    {\lambda ^2} - \text{TR}(J(\overline{E}))\lambda + \text{DET}(J(\overline{E})) = 0,

    where

    \text{TR}(J(\overline{E})) = {{ - 3{\overline{x}^2} + 2\left( {\beta + 1} \right){\overline{x}} - \beta - \alpha \overline{y}}} + \frac{{\gamma \overline{x}}}{{1 + \delta \overline{x}}} + \frac{{{s_1}}}{{{{\left( {1 + f\overline{y}} \right)}^2}}} - {d_1},
    \text{DET}(J(\overline{E})) = \left[ {{ - 3{\overline{x}} + 2\left( {\beta + 1} \right){\overline{x}} - \beta - \alpha \overline{y}}}\right]\left[ {\frac{{\gamma \overline{x}}}{{1 + \delta \overline{x}}} + \frac{{{s_1}}}{{{{\left( {1 + f\overline{y}} \right)}^2}}} - {d_1}} \right] + \alpha \overline{x}\left( {\frac{{\gamma \overline{y}}}{{{{\left( {1 + \delta \overline{x}} \right)}^2}}}} \right).

    1) Boundary equilibria

    At O(0, 0) , there is

    J({O}) = \left[ {\begin{array}{*{20}{c}} {{-\beta}}&0\\\\ 0&{{s_1} - {d_1}} \end{array}} \right],

    since \lambda_1 = {{-\beta}} < 0 , \lambda_2 = s_1-d_1 > 0 , then {O} is an unstable higher-order singularity.

    At {E^1}\left(0, \frac{1}{f}\left({\frac{{{s_1} - {d_1}}}{{{d_1}}}} \right)\right) , there is

    J({E^1}) = \left[ {\begin{array}{*{20}{c}} { {{-\beta}}- \frac{\alpha }{f}\left( {\frac{{{s_1} - {d_1}}}{{{d_1}}}} \right)}&0\\ { \frac{\gamma }{f}\left( {\frac{{{s_1} - {d_1}}}{{{d_1}}}} \right)}&0 \end{array}} \right],

    since \lambda_1 = 0 , {\lambda _2} = - \frac{\alpha }{f}\left({\frac{{{s_1} - {d_1}}}{{{d_1}}}} \right) < 0 , so that {E^1} is locally stable.

    At {E^2}(\beta, 0) , there is

    J({E^2}) = \left[ {\begin{array}{*{20}{c}} {{\beta(1-\beta)}}&{ - \alpha {{\beta}}}\\ 0&{ \frac{{\gamma \beta }}{{1 + \delta \beta }} + {s_1} - {d_1}} \end{array}} \right],

    since \lambda_1 = {{\beta(1-\beta)}} > 0 , {\lambda _2} = \frac{{\gamma \beta }}{{1 + \delta \beta }} + {s_1} - {d_1} > 0 , then {E^2} is unstable.

    At {E^3}(1, 0) , there is

    J({E^3}) = \left[ {\begin{array}{*{20}{c}} {{\beta-1}}&{ - \alpha }\\ 0&{{ \frac{\gamma }{{1 + \delta }} + {s_1} - {d_1}} }\end{array}} \right],

    since \lambda_1 = {{\beta-1}} < 0 , \lambda _2 = \frac{\gamma }{{1 + \delta }} + {s_1} - {d_1} > 0 , then {E^3} is unstable.

    2) Interior equilibrium

    At E^*(x^*, y^*) , there are g_1(x^*, y^*) = 0 and g_2(x^*, y^*) = 0 , then

    J(E^*) = \left[ {\begin{array}{*{20}{c}} {(1 + \beta )x^* - 2{(x^*)^2}}&{ - \alpha x^*}\\ { \frac{{\gamma y^*}}{{{{\left( {1 + \delta x^*} \right)}^2}}}}&{ - \frac{{{s_1}fy^*}}{{{{\left( {1 + fy^*} \right)}^2}}}} \end{array}} \right],

    thus,

    \text{DET}(J(E)) = \left[ {x^*y^* \frac{{\partial {g_1}}}{{\partial y}}\frac{{\partial {g_2}}}{{\partial y}}\left( { \frac{{\text{d}{y_1}}}{{\text{d}x}} - \frac{{\text{d}{y_2}}}{{\text{d}x}}}\right)} \right]_{(x^*, y^*)},

    where x^*y^*\frac{{\partial {g_1}}}{{\partial y}}\frac{{\partial {g_2}}}{{\partial y}}|_{x^*, y^*} > 0 , then the sign of \text{DET}(J(E^*)) is identical to that of \frac{{\text{d}{y_1}}}{{\text{d}x}}|_{x = x^*} - \frac{{\text{d}{y_2}}}{{\text{d}x}}|_{x = x^*} . Next, it will discuss the sign of \frac{{\text{d}{y_1}}}{{\text{d}x}}|_{x = x^*} - \frac{{\text{d}{y_2}}}{{\text{d}x}}|_{x = x^*} for different cases in Theorem 2.

    ⅰ) When \rho < 0 holds, two interior equilibria {E^*_{1}}\left({{x^*_{1}}, {y^*_{1}}} \right) and {E^*_{2}}\left({{x^*_{2}}, {y^*_{2}}} \right) with 0 < {x^*_{1}} < {x^*_{2}} < 1 exists in Model (2.2), as illustrated in Figure 2(a). It can be easily checked that

    \text{SIGN}(J({E^*_{1}})) = \left[ {\begin{array}{*{20}{c}} + & - \\ + & - \end{array}} \right], \text{SIGN}(J({E^*_{2}})) = \left[ {\begin{array}{*{20}{c}} - & - \\ + & - \end{array}} \right].
    Figure 2.  Symbolic representation of Jacobian matrix elements for the case C_1 ).

    Besides, at {E^*_{{1}}} , there is \frac{{\text{d}{y_1}}}{{\text{d}x}}|_{x = x^*_1} > \frac{{\text{d}{y_2}}}{{\text{d}x}}|_{x = x^*_1} , thus

    \text{DET}(J({E^*_{1}})) = {\left[ {xy\frac{{\partial {f_1}}}{{\partial y}}\frac{{\partial {f_2}}}{{\partial y}}\left( {\frac{{\text{d}{y_1}}}{{\text{d}x}} - \frac{{\text{d}{y_2}}}{{\text{d}x}}} \right)} \right]_{({x^*_{1}}, {y^*_{1}})}} < 0,

    i.e., {E^*_{1}} is unstable. Similarly, at {E^*_{2}} , there is \frac{{\text{d}{y_1}}}{{\text{d}x}}|_{x = x^*_2} < \frac{{\text{d}{y_2}}}{{\text{d}x}}|_{x = x^*_2} , thus

    ({\lambda _1} + {\lambda _2})\left| {_{({x^*_{2}}, {y^*_{2}})}} \right. = \text{TR}(J({E^*_{2}})) < 0,
    ({\lambda _1}{\lambda _2})\left| {_{({x^*_{2}}, {y^*_{2}})}} \right. = \text{DET}(J({E^*_{2}})) = {\left[ {xy\frac{{\partial {f_1}}}{{\partial y}}\frac{{\partial {f_2}}}{{\partial y}}\left( {\frac{{\text{d}{y_1}}}{{\text{d}x}} - \frac{{\text{d}{y_2}}}{{\text{d}x}}} \right)} \right]_{({x^*_{2}}, {y^*_{2}})}} > 0,

    i.e., {E^*_{2}} is a locally asymptotically stable node.

    ⅱ) When \rho = 0 holds, system (2.2) has a unique interior equilibrium {E^*}\left({{x^*}, {y^*}} \right) , as shown in Figure 2(b). Let X = x - {x^*}, Y = y - {y^*} , then {E^*} is converted to the origin O(0, 0) , and the model is written as

    \begin{equation} \left\{ \begin{array}{l} \frac{{\text{d}X}}{{\text{d}t}} = {a_{11}}X + {a_{12}}Y + {A_1}{X^2} + {A_2}XY, \\ \frac{{\text{d}Y}}{{\text{d}t}} = {a_{21}}X + {a_{22}}Y + {B_1}{X^2} + {B_2}XY + {B_3}{Y^2} + {P_3}(X, Y), \end{array} \right. \end{equation} (3.5)

    where,

    \begin{array}{l} {a_{11}} = \left( {1 + \beta } \right){x^*} - 2(x^*)^2, {a_{12}} = - \alpha {x^*}, {a_{21}} = \frac{{\gamma {y^*}}}{{{{\left( {1 + \delta {x^*}} \right)}^2}}}, {a_{22}} = - \frac{{{s_1}f{y^*}}}{{{{\left( {1 + f{y^*}} \right)}^2}}}, \\ {A_1} = \beta + 1 - 3{x^*}, {A_2} = - \frac{\alpha }{2}, {B_1} = - \frac{{\gamma \delta {y^*}}}{{{{\left( {1 + \delta {x^*}} \right)}^3}}}, {B_2} = \frac{\gamma }{{{{\left( {1 + \delta {x^*}} \right)}^2}}}, {{{B_3} = - \frac{{{s_1}f}}{{{{\left( {1 + f{y^*}} \right)}^3}}}}} \end{array}

    and {P_3}(X, Y) is a function of (X, Y) with degree of three or higher. The Jacobian matrix at {E^*} is

    J({E^*}) = \left[ {\begin{array}{*{20}{c}} {{a_{11}}}&{{a_{12}}}\\ {{a_{21}}}&{{a_{22}}} \end{array}} \right],

    thus

    \text{DET}(J({E^*})) = {a_{11}}{a_{22}} - {a_{12}}{a_{21}} = 0, \text{TR}(J({E^*})) = {a_{11}} + {a_{22}}.

    a) If \text{TR}(J({E^*})) = {a_{11}} + {a_{22}} = 0 , then {\lambda _1} = {\lambda _2} = 0 . Make the transformation

    \left( {\begin{array}{*{20}{c}} X\\ Y \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{a_{11}}}&0\\ {{a_{21}}}&1 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{x_1}}\\ {{y_1}} \end{array}} \right),

    then Model (3.5) is converted into the following standard form:

    \left\{ \begin{array}{l} \frac{{\text{d}{x_1}}}{{\text{d}t}} = {{\bar a}_{12}}{y_1} + {{\bar A}_1}x_1^2 + {{\bar A}_2}{x_1}{y_1}, \\ \frac{{\text{d}{y_1}}}{{\text{d}t}} = {{\bar B}_1}x_1^2 + {{\bar B}_2}{x_1}{y_1} + {{\bar B}_3}y_1^2 + {{\bar P}_3}({x_1}, {y_1}), \end{array} \right.

    where

    {\bar a_{12}} = \frac{{{a_{12}}}}{{{a_{11}}}}, {\bar A_1} = {a_{11}}{A_1} + {a_{21}}{A_2}, {\bar A_2} = {A_2},
    {\bar B_1} = a_{11}^2{B_1} + {a_{11}}{a_{21}}\left( {{B_2} - {A_1}} \right) + a_{21}^2\left( {{B_3} - {A_2}} \right), {\bar B_2} = {a_{11}}{B_2} + {a_{21}}\left( {2{B_3} - {A_2}} \right), {\bar B_3} = {B_3},

    and {\bar P_3}(X, Y) is a function of three or more degrees about (X, Y) .

    Let \tau = {\bar a_{12}}t . For convenience, t is still used to represent \tau , then it has

    \begin{equation} \left\{ \begin{array}{l} \frac{{\text{d}{x_1}}}{{\text{d}t}} = {y_1} + {{\tilde A}_1}x_1^2 + {{\tilde A}_2}{x_1}{y_1}, \\ \frac{{\text{d}{y_1}}}{{\text{d}t}} = {{\tilde B}_1}x_1^2 + {{\tilde B}_2}{x_1}{y_1} + {{\tilde B}_3}y_1^2 + {{\tilde P}_3}({x_1}, {y_1}), \end{array} \right. \end{equation} (3.6)

    where

    {\tilde A_i} = \frac{{{{\bar A}_i}}}{{{{\bar a}_{12}}}}, i = 1, 2, {{\rm{\tilde B}}_{\rm{j}}} = \frac{{{{{\rm{\bar B}}}_{\rm{j}}}}}{{{{\bar a}_{12}}}}, j = 1, 2, 3.

    and {\tilde P_3}(X, Y) is a function of (X, Y) with three or higher degree. Model (3.6) can be transformed into the following form[43]:

    \left\{ \begin{array}{l} \frac{{\text{d}{x_1}}}{{\text{d}t}} = {y_1}, \\ \frac{{\text{d}{y_1}}}{{\text{d}t}} = {{\tilde B}_1}x_1^2 + \left( {{{\tilde B}_2} + 2{{\tilde A}_1}} \right){x_1}{y_1} + {{\tilde P'}_3}({x_1}, {y_1}), \end{array} \right.

    and if {\tilde B_1} \ne 0 , then {E^*} is a cusp point.

    Meanwhile, if {\tilde B_2} + 2{\tilde A_1} \ne 0 , {E^*} is a cusp of codimension two. If {\tilde B_2} + 2{\tilde A_1} = 0 , {E^*} is a cusp with at least codimension three.

    b) If \text{TR}(J({E^*})) = {a_{11}} + {a_{22}}\ne 0 , then {\lambda _1} = 0 , {\lambda _2} \ne 0 . Make the transformation

    \left( {\begin{array}{*{20}{c}} X\\ Y \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{a_{22}}}&{{a_{11}}}\\ { - {a_{21}}}&{{a_{21}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{x_1}}\\ {{y_1}} \end{array}} \right),

    system (3.5) is converted into the following standard form

    \left\{ \begin{array}{l} \frac{{\text{d}{x_1}}}{{\text{d}t}} = {{A'}_1}x_1^2 + {{A'}_2}{x_1}{y_1} + A'y_1^2, \\ \frac{{\text{d}{y_1}}}{{\text{d}t}} = {{a'}_{22}}{y_1} + {{B'}_1}x_1^2 + {{B'}_2}{x_1}{y_1} + {{B'}_3}y_1^2 + {{P'}_3}({x_1}, {y_1}), \end{array} \right.

    where

    \begin{array}{ll} {{A'}_1} & = \frac{{{a_{21}}a_{22}^2{A_1} - a_{21}^2{a_{22}}{A_2} - {a_{11}}a_{22}^2{B_1} + {a_{11}}{a_{21}}{a_{22}}{B_2} - {a_{11}}a_{21}^2{B_3}}}{{{a_{21}}\left( {{a_{11}} + {a_{22}}} \right)}}, \\ {{A'}_2} & = \frac{{2a_{11}^2{a_{22}} + a_{21}^2{a_{22}}{A_2} - {a_{11}}a_{21}^2{A_2} - 2a_{11}^2{a_{22}}{B_1} - {a_{11}}{a_{21}}{a_{22}}{B_2} + {a_{11}}a_{21}^2{B_2} + 2{a_{11}}a_{21}^2{B_3}}}{{{a_{21}}\left( {{a_{11}} + {a_{22}}} \right)}}, \\ {{A'}_3} & = \frac{{a_{11}^2{a_{21}}{A_1} + {a_{11}}a_{21}^2{A_2} - a_{11}^3{B_1} - a_{11}^2{a_{21}}{B_2} - {a_{11}}a_{21}^2{B_3}}}{{{a_{21}}\left( {{a_{11}} + {a_{22}}} \right)}}, \\ {{B'}_1} & = \frac{{{a_{21}}a_{22}^2{A_1} - a_{21}^2{a_{22}}{A_2} + a_{22}^3{B_1} - {a_{21}}a_{22}^2{B_2} + a_{21}^2{a_{22}}{B_3}}}{{{a_{21}}\left( {{a_{11}} + {a_{22}}} \right)}}, \\ {{B'}_2} & = \frac{{2{a_{11}}{a_{21}}{a_{22}}{A_2} + a_{21}^2{a_{22}}{A_2} - {a_{11}}a_{21}^2{A_2} + 2{a_{11}}a_{22}^2{B_1} + {a_{21}}a_{22}^2{B_2} - {a_{11}}{a_{21}}{a_{22}}{B_2} - 2a_{21}^2{a_{22}}{B_3}}}{{{a_{21}}\left( {{a_{11}} + {a_{22}}} \right)}}, \\ {{B'}_3} & = \frac{{a_{11}^2{a_{21}}{A_1} + {a_{11}}a_{21}^2{A_2} + {a_{11}}a_{22}^2{B_1} + {a_{11}}{a_{21}}{a_{22}}{B_2} + a_{21}^2{a_{22}}{B_3}}}{{{a_{21}}\left( {{a_{11}} + {a_{22}}} \right)}}. \end{array}

    Next, introduce a new variable \tau = {a'_{22}}t (for convenience, is represented by t ), then

    \left\{ \begin{array}{l} \frac{{\text{d}{x_1}}}{{\text{d}t}} = {{\hat A}_1}x_1^2 + {{\hat A}_2}{x_1}{y_1} + {{\hat A}_3}y_1^2, \\ \frac{{\text{d}{y_1}}}{{\text{d}t}} = {y_1} + {{\hat B}_1}x_1^2 + {{\hat B}_2}{x_1}{y_1} + {{\hat B}_3}y_1^2 + {{\hat P}_3}({x_1}, {y_1}), \end{array} \right.

    where, {\hat A_i} = \frac{{{{A'}_i}}}{{{{a'}_{22}}}}, {\hat B_i} = \frac{{{{B'}_i}}}{{{{a'}_{22}}}}, i = 1, 2, 3 .

    If {\hat A_1} \ne 0 then {E^*} is unstable. Meanwhile, {E^*} is a saddle node of attraction.

    To sum up, the following result hold.

    Theorem 3. For Model (2.2), 1) {O}(0, 0) is unstable, {E^1}(0, \frac{1}{f}\left({\frac{{{s_1} - {d_1}}}{{{d_1}}}} \right)) is locally stable, {E^2}(\beta, 0) , {E^3}(1, 0) is unstable; 2) For the interior equilibrium, when \rho > 0 , {E^*_{{1}}}\left({{x^*_{{1}}}, {y^*_{{1}}}} \right) is a saddle (unstable), and {E^*_{{2}}}\left({{x^*_{{2}}}, {y^*_{{2}}}} \right) is a stable node; when \rho = 0 , if \text{TR}(J({E^*})) = 0 and {\tilde B_1} \ne 0 , {E^*}\left({{x^*}, {y^*}} \right) is a cusp of codimension two in case of {\tilde B_2} + 2{\tilde A_1} \ne 0 , and a cusp with at least codimension three in case of {\tilde B_2} + 2{\tilde A_1} = 0 ; If \text{TR}(J({E^*})) \ne 0 and {\hat A_1} \ne 0 , {E^*}\left({{x^*}, {y^*}} \right) is an attractive saddle node.

    Let \alpha = {\alpha _0} satisfy

    \text{DET}(J({E^*})) = ({\lambda _1}{\lambda _2})\left| {_{({x^*}, {y^*})}} \right. = 0, \text{TR}(J({E^*})) = - ({\lambda _1} + {\lambda _2})\left| {_{({x^*}, {y^*})}} \right. \ne 0.

    Denote

    {\xi_1} = 2(x^*)^2 - (1 + \beta ){x^*}, {{{\xi_2} = - \frac{{\gamma y^*}}{{{{\left({1 + \delta x^*} \right)}^2}}}}}, {\omega _1} = \alpha {x^*}, {\omega_2} = \frac{{{s_1}fy^*}}{{{{\left( {1 + fy^*} \right)}^2}}}

    and define

    {\Phi_2} \buildrel \Delta \over = 2\omega _1^2{\xi_2}x^* - \frac{{{2\gamma \delta y^*}}}{{{{\left( {1 + \delta {x^*}} \right)}^3}}}\omega_1^2{\xi_1} - \frac{{2{s_1}{f^2}}}{{{{\left( {1 + f{y^*}} \right)}^3}}}\xi_1^3.

    Theorem 4. Let the parameters of Model (2.2) satisfy \text{TR}(J({E^*})) \ne 0 and {\hat A_1} \ne 0 . If {\Phi_2} \ne 0 , system (2.2) undergoes a saddle node bifurcation near {E^*}\left({{x^*}, {y^*}} \right) when \alpha = {\alpha _0} .

    Proof. At {E^*}\left({{x^*}, {y^*}} \right) , there is

    J\left( {{E^*}} \right) = \left( {\begin{array}{*{20}{c}} { - \xi {}_1}&{ - {\omega _1}}\\ { - \xi {}_2}&{ - {\omega _2}} \end{array}} \right).

    Let V = {\left({{V_1}, {V_2}} \right)^T} ( W = {\left({{W_1}, {W_2}} \right)^T} ) be the eigenvector corresponding to the zero eigenvalue of J\left({{E^*}} \right) ( {\left({J({E^*})} \right)^T} ). Then V = {\left({{\omega_1}, - {\xi_1}} \right)^T} , W = {\left(- {\xi_2}, {{\xi_1}} \right)^T} . Let g = {\left({{g_1}, {g_2}} \right)^T} . Then

    {g_\alpha }\left( {{E^*}, {\alpha_0}} \right) = \frac{{\partial g}}{{\partial \alpha }}\left( {{E^*}, {\alpha_0}} \right) = \left( {\begin{array}{*{20}{c}} {{-{y^*}}}\\ 0 \end{array}} \right).

    and

    \begin{align*} {D^2}g\left( {{E^*}, {\alpha _0}} \right)\left( {V, V} \right) & = \left( {\begin{array}{*{20}{c}} { \frac{{{\partial ^2}{g_1}}}{{\partial {x^2}}}V_1^2 + 2\frac{{{\partial ^2}{g_1}}}{{\partial x\partial y}}{V_1}{V_2} + \frac{{{\partial ^2}{g_1}}}{{\partial {y^2}}}V_2^2}\\ { \frac{{{\partial ^2}{g_2}}}{{\partial {x^2}}}V_1^2 + 2\frac{{{\partial ^2}{g_2}}}{{\partial x\partial y}}{V_1}{V_2} + \frac{{{\partial ^2}{g_2}}}{{\partial {y^2}}}V_2^2} \end{array}} \right) \\ & = {{\left( {\begin{array}{*{20}{c}} { - 2\omega_1^2}\\ { - \frac{{2\gamma \delta}}{{{{\left( {1 + \delta {x^*}} \right)}^3}}}\omega _1^2 - \frac{{2{s_1}{f^2}}}{{{{\left( {1 + f{y^*}} \right)}^3}}}\xi _1^2} \end{array}} \right).}} \end{align*}

    Since

    {\Phi_1} = {W^T}{g_\alpha }\left( {{E^*}, {\alpha _0}} \right) = {{-{{x^*}}{\xi_1} \ne 0}}

    and

    \begin{align*} {\Phi_2} & = {W^T}\left[ {{D^2}g\left( {{E^*}, {\alpha _0}} \right)\left( {V, V} \right)} \right] \\ & = 2\omega _1^2{\xi_2} - \frac{{2\gamma \delta}}{{{{\left( {1 + \delta {x^*}} \right)}^3}}}\omega _1^2{\xi_1} - \frac{{2{s_1}{f^2}}}{{{{\left( {1 + f{y^*}} \right)}^3}}}\xi _1^3 \ne 0, \end{align*}

    then according to Sotomayor's theorem[44], Model (2.2) undergoes a saddle-node bifurcation near {E^*}\left({{x^*}, {y^*}} \right) when \alpha = {\alpha_0} .

    Remark 1. For Model (2.2), it can be concluded that for a lager \alpha , the interior equilibrium does not exist. When \alpha decreases to \alpha = \alpha_0 , a unique equilibrium exists in the system. As \alpha decreases below \alpha_0 , the system undergoes a saddle-node bifurcation at the interior equilibrium {E^*} , giving rise to two interior equilibria {E^*_{{1}}} and {E^*_{{2}}} .

    According to Theorem 3, when \rho = 0 , \text{TR}(J({E^*})) = 0 , {\tilde B_1} \ne 0 and {\tilde B_2} + 2{\tilde A_1} \ne 0 , Model (2.2) to undergo a Bogdanov-Takens bifurcation of codimension two near {E^*} when \left({\alpha, \beta } \right) = \left({{\alpha _0}, {\beta _0}} \right) . Next, it will show the universal unfolding of the Bodmanov-Takens bifurcation of codimension two under parameter perturbation when \alpha and \beta are taken as bifurcation parameters.

    Define

    \begin{array}{l} {b_{01}} = {x^*}(1 - {x^*})\left( {{x^*} - {\epsilon_2}} \right) - {\epsilon_1}{x^*}{y^*}, {b_{11}} = \left[ {1 + \left( {\beta + {\epsilon_2}} \right)} \right]{x^*} - 2(x^*)^2, {b_{12}} = - \left( {\alpha + {\epsilon_1}} \right){x^*}, \\ {b_{10}} = 0, {b_{21}} = \frac{{\gamma {y^*}}}{{{{\left( {1 + \delta {x^*}} \right)}^2}}}, {b_{21}} = - \frac{{{s_1}f{y^*}}}{{{{\left( {1 + f{y^*}} \right)}^2}}}, {E_1} = \left( {\beta + {\epsilon_2}} \right) + 1 - 3{x^*}, {E_2} = - \frac{{\left( {\alpha + {\epsilon_1}} \right)}}{2}, \\ {F_1} = - \frac{{\gamma \delta {y^*}}}{{{{\left( {1 + \delta {x^*}} \right)}^3}}}, {F_2} = \frac{\gamma }{{{{\left( {1 + \delta {x^*}} \right)}^2}}}, {{{F_3} = - \frac{{{s_1}f}}{{{{\left( {1 + f{y^*}} \right)}^3}}}}}, \\ {H_1} = \frac{{ - b_{12}^2{b_{10}} + {b_{01}}{b_{12}}{b_{22}} - b_{01}^2{F_3}}}{{{b_{12}}}}, \\ {H_2} = \frac{{b_{01}^2E_2^2 + b_{12}^2{b_{10}}E_2^2 + b_{12}^3{b_{21}} - {b_{11}}b_{12}^2{b_{22}} + {b_{01}}b_{12}^2{F_2}}}{{b_{12}^2}}, \\ {H_3} = \frac{{{b_{11}}{b_{12}} - {b_{01}}{E_2} + {b_{12}}{b_{22}} - 2{b_{10}}{F_3}}}{{{b_{12}}}}, \\ {H_4} = - \frac{{{b_{01}}{b_{12}}{E_1}{E_2} - {b_{01}}{b_{11}}E_2^2 - b_{12}^2{b_{22}}{E_1} - b_{12}^2{F_1} + {b_{11}}b_{12}^2{F_2} - b_{11}^2{b_{12}}{F_3}}}{{b_{12}^2}}, \\ {H_5} = \frac{{2b_{12}^2E_1^2 - {b_{11}}{b_{12}}{E_2} - 2{b_{01}}E_2^2 + b_{12}^2{F_2} - 2{b_{11}}{b_{12}}{F_3}}}{{b_{12}^2}}, {H_6} = \frac{{{E_2} + {F_3}}}{{{b_{12}}}}, \\ {J_1} = {H_1}, {J_2} = {H_2} - 2{H_1}{H_6}, {J_3} = {H_3}\\ {J_4} = {H_1}H_6^2 - 2{H_2}{H_6} + {H_4}, {J_5} = {H_5} - {H_3}{H_6}, \\ \end{array}
    \begin{array}{l} {K_1} = - \frac{{{J_1}}}{{{J_4}}}, {K_2} = - \frac{{{J_2}}}{{{J_4}}}, {K_3} = - \frac{{{J_3}}}{{\sqrt { - {J_4}} }}, {K_4} = - \frac{{{J_5}}}{{\sqrt { - {J_4}} }}, \\ {L_1} = {K_1} + \frac{{K_2^2}}{4}, {L_2} = {K_3} + \frac{{{K_2}{K_4}}}{2}, {L_3} = {K_4}, \\ {O_1} = - {L_1}L_3^4, {O_2} = - {L_2}{L_3}, \\ \end{array}

    Let \left({{\epsilon_1}, {\epsilon_2}} \right) be a parameter vector in a small neighborhood of (0, 0) . Then

    Theorem 5. Let the parameters of Model (2.2) satisfy \rho = 0 , \text{TR}(J({E^*})) = 0 , {\tilde B_1} \ne 0 , {\tilde B_2} + 2{\tilde A_1} \ne 0 and \left| {\frac{{\partial \left({{O_1}, {O_2}} \right)}}{{\partial \left({{\epsilon_1}, {\epsilon_2}} \right)}}} \right| \ne 0 . When \left({\alpha, \beta } \right) varies in the neighborhood of \left({{\alpha _0}, {\beta _0}} \right) , Model (2.2) changes in the small neighborhood of E^*({x^*}, {y^*}) , and a codimensional 2 Bogdanov-Takens branching occurs.

    Proof. Consider the perturbation system

    \left\{ \begin{array}{l} \frac{{\text{d}x}}{{\text{d}t}} = x(1 - x)\left[ {x - \left( {\gamma + {\epsilon_2}} \right)} \right] - \left( {\alpha + {\epsilon_1}} \right)xy \triangleq F\left( {x, y} \right), \\ \frac{{\text{d}y}}{{\text{d}t}} = \frac{{\gamma xy}}{{1 + \delta x}} + \left( {\frac{{{s_1}}}{{1 + fy}} - {d_1}} \right)y \triangleq G\left( {x, y} \right), \end{array} \right.

    For \left({\alpha, \beta } \right) = \left({{\alpha _0}, {\beta _0}} \right) , there is \text{DET}(J({E^*})) = 0, \text{TR}(J({E^*})) = 0 . With the transformation {x_1} = x - {x^*} , {y_1} = y - {y^*} , we can obtain

    \begin{equation} \left\{ \begin{array}{l} \frac{{\text{d}{x_1}}}{{\text{d}t}} = {b_{01}} + {b_{11}}{x_1} + {b_{12}}{y_1} + {E_1}x_1^2 + {E_2}{x_1}{y_1}, \\ \frac{{\text{d}{y_1}}}{{\text{d}t}} = {b_{10}} + {b_{21}}{x_1} + {b_{22}}{y_1} + {F_1}x_1^2 + {F_2}{x_1}{y_1} + {F_3}y_1^2 + {N_3}({x_1}, {y_1}), \end{array} \right. \end{equation} (4.1)

    where {N_3}({x_1}, {y_1}, {\epsilon_1}, {\epsilon_2})\in{C^\infty } .

    Make the following transformation:

    \left\{ \begin{array}{l} {x_2} = {x_1}, \\ {y_2} = {b_{01}} + {b_{11}}{x_1} + {b_{12}}{y_1} + {E_1}x_1^2 + {E_2}{x_1}{y_1}, \end{array} \right.

    then Model (4.1) is converted to

    \left\{ \begin{array}{l} \frac{{\text{d}{x_2}}}{{\text{d}t}} = {y_2}, \\ \frac{{\text{d}{y_2}}}{{\text{d}t}} = {H_1} + {H_2}{x_2} + {H_3}{y_2} + {H_5}x_2^2 + {H_5}{x_2}{y_2} + {H_6}y_2^2 + {{N'_3}}({x_2}, {y_2}, {\epsilon_1}, {\epsilon_2}), \end{array} \right.

    where {N'_3}({x_2}, {y_2}, {\epsilon_1}, {\epsilon_2})\in{C^\infty } with coefficients depending smoothly on {\epsilon_1}, {\epsilon_2} .

    Next, introduce the variable \tau , denoted as \text{d}t = \left({1-{H_6}{x_2}} \right)\text{d}\tau (still denote \tau as t ), then

    \left\{ \begin{array}{l} \frac{{\text{d}{x_2}}}{{\text{d}t}} = \left( {1 - {H_6}{x_2}} \right){y_2}, \\ \frac{{\text{d}{y_2}}}{{\text{d}t}} = \left( {1 - {H_6}{x_2}} \right)\left( {{H_1} + {H_2}{x_2} + {H_3}{y_2} + {H_4}x_2^2 + {H_5}{x_2}{y_2} + {H_6}y_2^2 + {{N'}_3}({x_2}, {y_2}, {\lambda _1}, {\lambda _2})} \right), \end{array} \right.

    Let {x_3} = {x_2}, {y_3} = \left({1 - {H_6}{x_2}} \right){y_2} , then the above system of equations is transformed into

    \left\{ \begin{array}{l} \frac{{\text{d}{x_3}}}{{\text{d}t}} = {y_3}, \\ \frac{{\text{d}{y_3}}}{{\text{d}t}} = {J_1} + {J_2}{x_3} + {J_3}{y_3} + {J_4}x_3^2 + {J_5}{x_3}{y_3} + {{\tilde N}_3}({x_3}, {y_3}, {\epsilon_1}, {\epsilon_2}), \end{array} \right.

    where {\tilde N_3}({x_3}, {y_3}, {\epsilon_1}, {\epsilon_2})\in{C^\infty } with coefficients depending smoothly on {\epsilon_1}, {\epsilon_2} .

    (ⅰ) When {J_4} < 0 , the following transformations are applied to the variables:

    {x_4} = {x_3}, {y_4} = \frac{{{y_3}}}{{\sqrt { - {J_4}} }}, \tau = \sqrt { - {J_4}} t,

    still denote \tau as t , there is

    \left\{ \begin{array}{l} \frac{{\text{d}{x_4}}}{{\text{d}t}} = {y_4}, \\ \frac{{\text{d}{y_4}}}{{\text{d}t}} = {K_1} + {K_2}{x_4} + {K_3}{y_4} - x_4^2 + {K_4}{x_4}{y_4} + {M_3}({x_4}, {y_4}, {\epsilon_1}, {\epsilon_2}), \end{array} \right.

    where {M_3}({x_4}, {y_4}, {\epsilon_1}, {\epsilon_2})\in{C^\infty } with at least third order.

    Let {x_5} = {x_4} - \frac{{K_2}}{2} , {y_5} = {y_4} , and obtain

    \left\{ \begin{array}{l} \frac{{\text{d}{x_5}}}{{\text{d}t}} = {y_5}, \\ \frac{{\text{d}{y_5}}}{{\text{d}t}} = {L_1} + {L_2}{y_5} - x_5^2 + {L_3}{x_5}{y_5} + {{M'_3}}({x_5}, {y_5}, {\epsilon_1}, {\epsilon_2}), \end{array} \right.

    where {M'_3}({x_5}, {y_5}, {\epsilon_1}, {\epsilon_2})\in{C^\infty } with at least third order.

    If set {J_5} \ne 0 , then {L_3} \ne 0 . Define new variables: {x_6} = - L_3^2{x_5} , {y_6} = L_3^3{y_5} , \tau = - \frac{t}{{L{}_3}} , and denote {x_6} by x , {y_6} by y , and \tau by t , which yields that

    \begin{equation} \left\{ \begin{array}{l} \frac{{\text{d}x}}{{\text{d}t}} = y, \\ \frac{{\text{d}y}}{{\text{d}t}} = {O_1} + {O_2}y + {x^2} + xy + {{\tilde M}_3}(x, y, {\lambda _1}, {\lambda _2}), \end{array} \right. \end{equation} (4.2)

    where {\tilde N_3}({x_2}, {y_2}, {\epsilon_1}, {\epsilon_2})\in{C^\infty } with at least third order, and {O_1} , {O_2} can be represented by \epsilon_1 and \epsilon_2 .

    (ⅱ) When {J_4} > 0 , the following transformations are applied

    {x_4} = {x_3}, {y_4} = \frac{{{y_3}}}{{\sqrt { {J_4}} }}, \tau = \sqrt { {J_4}} t,

    still denote \tau by t , it has

    \left\{ \begin{array}{l} \frac{{\text{d}{x_4}}}{{\text{d}t}} = {y_4}, \\ \frac{{\text{d}{y_4}}}{{\text{d}t}} = {{{K'_1} + {{K'_2}}{x_4} + {{K'_3}}{y_4} + x_4^2 + {{K'_4}}{x_4}{y_4} + {N_3}({x_4}, {y_4}, {\epsilon_1}, {\epsilon_2})}}, \end{array} \right.

    where {N_3}({x_4}, {y_4}, {\epsilon_1}, {\epsilon_2})\in{C^\infty } with at least third order, and

    {K'_1} = \frac{{{J_1}}}{{{J_4}}}, {K'_2} = \frac{{{J_2}}}{{{J_4}}}, {K'_3} = \frac{{{J_3}}}{{\sqrt {{J_4}} }}, {K'_4} = \frac{{{J_5}}}{{\sqrt {{J_4}} }}.

    Let {x_5} = {x_4} + {{{{K'_2}}}}/{2}, {y_5} = {y_4} . Then it has

    \left\{ \begin{array}{l} \frac{{\text{d}{x_5}}}{{\text{d}t}} = {y_5}, \\ \frac{{\text{d}{y_5}}}{{\text{d}t}} = {{{L'_1} + {{L'_2}}{y_5} + x_5^2 + {{L'_3}}{x_5}{y_5} + {{N'_3}}({x_5}, {y_5}, {\epsilon_1}, {\epsilon_2})}}, \end{array} \right.

    where {N'_3}({x_5}, {y_5}, {\epsilon_1}, {\epsilon_2})\in{C^\infty } with at least third order, and

    {{{L'_1} = {K'_1} - \frac{{{{K'_2}}^2}}{4}, {L'_2} = {K'_3} - \frac{{{{K'_2}}{{K'_4}}}}{2}, {L'_3} = {K'_4}}}.

    If set {J_5} \ne 0 , then {L'_3} \ne 0 . Define new variables: {x_6} = {L'_3}^2{x_5} , {y_6} = {L'_3}^3{y_5} , \tau = {t}/{{{{L'_3}}}} , and still denote {x_6} by x , {y_6} by y , \tau by t , which yields that

    \begin{equation} \left\{ \begin{array}{l} \frac{{\text{d}x}}{{\text{d}t}} = y, \\ \frac{{\text{d}y}}{{\text{d}t}} = {{O'}_1} + {{O'}_2}y + {x^2} + xy + {{\tilde N}_3}({x_2}, {y_2}, {\epsilon_1}, {\epsilon_2}), \end{array} \right. \end{equation} (4.3)

    where {O'_1} = {L'_1}{L'_3}^4, {O'_2} = {L'_2}{L'_3} , and {\epsilon_1} , {\epsilon_2} can be represented by {O'_1} , {O'_2} .

    For convenience of discussion, {O'_1} , {O'_2} is still denoted by {O_1}, {O_2} . When \left| {\frac{{\partial \left({{O_1}, {O_2}} \right)}}{{\partial \left({{\epsilon_1}, {\epsilon_2}} \right)}}} \right| \ne 0 , Models (4.2) and (4.3) are the cardinal folds of the Bogdanov-Takens bifurcation[44], when \left({\alpha, \beta } \right) varies in the vicinity of \left({{\alpha_0}, {\beta_0}} \right) , Model (2.2) undergoes a codimension 2 Bogdanov-Takens bifurcation in a small neighborhood of {E^*}\left({{x^*}, {y^*}} \right) .

    It only focuses on the case of \rho < 0 in Theorem 3. For Model (2.3), there are

    {\mathcal{M}_1} = \left\{ {\left( {x, y} \right) \in {\rm{R}}_ + ^2:x = {h_1}, y \ge 0} \right\}, {\mathcal{M}_2} = \left\{ {\left( {x, y} \right) \in R_ + ^2:x = {h_2}, y \ge 0} \right\},
    {\mathcal{N}_1} = \left\{ {\left( {x, y} \right) \in R_ + ^2:x = {h_1} + \eta , y \ge 0} \right\}, {\mathcal{N}_2} = \left\{ {\left( {x, y} \right) \in R_ + ^2:x = \left( {1 - a} \right){h_2}, y \ge \tau } \right\}.

    Denote {l_1} = \frac{1}{\alpha }\left({1 - x} \right)\left({x - \beta } \right) as the prey isocline, {l_2} = \frac{{\left({{s_1} - {d_1}} \right)\left({1 + \delta x} \right) - \gamma {x}}}{{f\left[{{d_1} + \left({{d_1}\delta - \gamma } \right){x}} \right]}} as the predator's isocline, {l_3} , {l_4} , {l_5} , {l_6} as the saddle point separatrix of E^*_1 in different directions. The intersection point between l_1 and \mathcal{N}_2 is denoted by A_3 . The intersection point between l_4 and \mathcal{N}_2 is denoted by M . The trajectory from A_3 intersects \mathcal{M}_2 at the point {B_3} . Define \overline{\tau}_2\triangleq y_{M}-(1-b)y_{B_3} . Denote

    {\Omega _1} = \left\{ {\left( {x, y} \right)\left| {0 \le x \le {x^*_{1}}, y \ge 0} \right.} \right\}, {\Omega _2} = \left\{ {\left( {x, y} \right)\left| {{x^*_{1}} \le x \le 1, y \ge 0} \right.} \right\}.

    Definition 2 (Successor function). For a point S\in \mathcal{N}_1 , if the trajectory from S intersects \mathcal{M}_1 , then denote the intersection point by S^{-}\in \mathcal{M}_1 . Under the impulse effect, the point S^{-} is mapped to S^{+}\in \mathcal{N}_1 . In such a case, we can define f^{\text{I}}_{\text{SOR}_1} : \mathcal{N}_1\rightarrow R , S\rightarrow f^{\text{I}}_{\text{SOR}_1}(S)\triangleq y_{S^{+}}-y_{S} . If the trajectory from S intersects \mathcal{M}_2 , then denote the intersection point by S^{-} . Under the impulse effect, the point S^{-} is mapped to S^{+}\in \mathcal{N}_2 . If the trajectory from S^{+} intersects \mathcal{M}_1 , denote the intersection point by S^{+-}\in \mathcal{M}_1 . Under the impulse effect, the point S^{+-} is mapped to S^{++}\in \mathcal{N}_1 . In such cases, we can define f^{\text{II}}_{\text{SOR}_1} : \mathcal{N}_1\rightarrow R , S\rightarrow f^{\text{II}}_{\text{SOR}_1}(S)\triangleq y_{S^{++}}-y_{S} . Similarly, For a point S'\in \mathcal{N}_2 , define f^{\text{I}}_{\text{SOR}_2} : \mathcal{N}_2\rightarrow R , S'\rightarrow f^{\text{I}}_{\text{SOR}_2}(S')\triangleq y_{S'^{+}}-y_{S'} .

    Theorem 6. For system (2.3) with model parameters satisfying any one of C_1) C_3) and \rho > 0 , if D-1) {h_1} + \eta < {x^*_{1}} < \left({1 - a} \right){h_2} < {h_2} < {x^*_{2}} , y_2(h_1)\geq y_1(h_1+\eta) and \tau < \overline{\tau}_2 , an order-1 periodic solution exists in each \Omega_i , i = 1, 2 ; if D-2) {h_1} < \left({1 - a} \right){h_2} < {h_1} + \eta < {x^*_{1}} < {h_2} < {x^*_{2}} and y_2(h_1)\geq y_1(h_1+\eta) , an order-1 periodic solution exists in \Omega_1 ; if D-3) {h_1} < {x^*_{1}} < {h_1} + \eta < \left({1 - a} \right){h_2} < {h_2} < {x^*_{2}} and \tau < \overline{\tau}_2 , an order-1 periodic solution exists in \Omega_2 .

    Proof. For case D-1) {h_1} + \eta < {x^*_{1}} < \left({1 - a} \right){h_2} < {h_2} < {x^*_{2}} and y_2(h_1)\geq y_1(h_1+\eta) , denote the intersection point between l_2 and \mathcal{N}_1 by A_0 . Select a point A_1\in \mathcal{N}_1 above A_0 ; the trajectory from A_1 intersects \mathcal{M}_1 at B_1 . Under the impulse effect, the point B_1 is mapped to A^{+}_1\in \mathcal{N}_1 . Since g_1(x_{A_1}, y_{A_1}) < 0 , g_2(x_{A_1}, y_{A_1}) < 0 , then we have f^{\text{I}}_{\text{SOR}_1}(A_1) = y_{A^{+}_1}-y_{A_1} = y_{B_1}-y_{A_1} < 0 . Besides, denote the intersection between l_1 and \mathcal{N}_1 by A_2 . Since g_1(x_{A_2}, y_{A_2}) = 0 , g_2(x_{A_2}, y_{A_2}) > 0 and g_1(x_S, y_S) < 0 for S\in U(A_1, \epsilon) with x_S < h_1+\eta and y_S < y_{A_1} , then we have f^{\text{I}}_{\text{SOR}_1}(A_2) = y_{A^{+}_2}-y_{A_2} = y_{B_2}-y_{A_2} > 0 due to y_2(h_1)\geq y_1(h_1+\eta) . The continuity of f^{\text{I}}_{\text{SOR}_1} implies that a point S\in \overline{A_1A_2} exists so that f^{\text{I}}_{\text{SOR}_1}(S) = 0 , i.e., an order-1 periodic solution exists in \Omega_1 (Figure 3(a)). Similarly, it can be proved that an order-1 periodic solution exists in \Omega_1 for case D-2).

    Figure 3.  Schematic diagram of the trajectory trend of the system's (2.3) for the case of D-1).

    Under the impulse effect, {B_3} is mapped to {A_3}^{+}\in{\mathcal{N}_2} , where {y_{{A_3}^ + }} = \left({1 - b} \right){y_{{B_3}}} + \tau . Define {\overline{\tau}_1} \triangleq {y_{{A_3}}} - \left({1 - b} \right){y_{{B_3}}} .

    ⅰ) If \tau = \overline{\tau}_1 , then f^{\text{I}}_{\text{SOR}_2}({A_3}) = {y_{{A_3}^ + }} - {y_{{A_3}}} = 0 , i.e., an order-1 periodic solution exists in \Omega_1 (Figure 3(a)).

    ⅱ) If \tau < \overline{\tau}_1 , then f^{\text{I}}_{\text{SOR}_2}({A_3}) = {y_{{A_3}^ + }} - {y_{{A_3}}} < 0 . On the other hand, for {A_4}((1-a)h_2, \tau) , there is f^{\text{I}}_{\text{SOR}_2}({A_4}) = {y_{{A_4}^ + }} - {y_{{A_4}}} > 0 . The continuity of f^{\text{I}}_{\text{SOR}_2} implies that a point S'\in \overline{A_3A_4} exists so that f^{\text{I}}_{\text{SOR}_2}(S') = 0 , i.e., an order-1 periodic solution exists in \Omega_2 (Figure 3(b)).

    ⅲ) If \overline{\tau}_1 < \tau < \overline{\tau}_2 , then f^{\text{I}}_{\text{SOR}_2}({A_3}) = {y_{{A_3}^ + }} - {y_{{A_3}}} > 0 . According to the trend of the trajectory and the fact that any two trajectories cannot be intersected, select D\in \mathcal{N}_2 above and sufficiently close to the {A_3} , then f^{\text{I}}_{\text{SOR}_2}(D) = {y_{{D^ + }}} - {y_D} > 0 . On the other hand, for C = {A^{+}_3} , there is f^{\text{I}}_{\text{SOR}_2}(C) = {y_{{C^ + }}} - {y_C} < {y_{{A^{+}_{3}}}} - {y_C} = 0 . The continuity of f^{\text{I}}_{\text{SOR}_2} implies that a point S'\in \overline{CD} exists so that f^{\text{I}}_{\text{SOR}_2}(S') = 0 , i.e., an order-1 periodic solution exists in \Omega_2 (Figure 3(c)).

    Similarly, system (2.3) has an order-1 periodic solution in \Omega_2 for case D - 3 ).

    To sum up, an order-1 periodic solution exists in Model (2.3) for case D-1).

    Let \mathbf{z}_i(t) = (\phi_i(t), \varphi_i(t)) (k-1)T_i\leq t\leq kT_i , k\in \mathbb{N} be the order-1 periodic solution in \Omega_i , i = 1, 2 . For \mathbf{z}_1(t) = (\phi_1(t), \varphi_1(t)) (k-1)T_1\leq t\leq kT_1 , denote

    \phi_1\left( {T_1^ + } \right) = \phi_1\left( {{T_1}} \right) + \eta = {h_1} + \eta, \varphi_1\left( {T_1^ + } \right) = \varphi_1\left( {{T_1}} \right) = {\delta_1}.

    For \mathbf{z}_2(t) = (\phi_2(t), \varphi_2(t)) (k-1)T_2\leq t\leq kT_2 , denote

    \phi_2\left( {T_2^ + } \right) = \left( {1 - a} \right)\phi_2 \left( {{T_2}} \right) = \left( {1 - a} \right){h_2}, \varphi_2 \left( {T_2^ + } \right) = \left( {1 - b} \right)\varphi_2 \left( {{T_2}} \right) + \tau = {\delta_2}.

    Theorem 7. For the model parameters with any one of C_1) C_3) , \rho > 0 and D-1), \mathbf{z}_i(t) = (\phi_i(t), \varphi_i(t)) (k-1)T_i\leq t\leq kT_i is orbitally asymptotically stable if \mu_i < 1 , where

    \begin{array}{l} \mu_1 = \left| {\frac{{\left[ {1 - \left( {{h_1} + \eta } \right)} \right]\left[ {\left( {{h_1} + \eta } \right) - \beta } \right] - \alpha {\delta _1}}}{{\left( {1 - {h_1}} \right)\left( {{h_1} - \beta } \right) - \alpha {\delta _1}}}} \right| \exp \left\{ {\int_{{0^ + }}^{{T_1}} {\left[ {\phi_1(t)\left( {\beta + 1 - 2\phi_1(t)} \right) - \frac{{{s_1}f\varphi_1(t)}}{{{{\left( {1 + f\varphi_1(t)} \right)}^2}}}} \right]} dt} \right\}, \\ \mu_2 = \left| {\frac{{\left\{ {\left[ {1 - \left( {1 - a} \right){h_2}} \right]\left[ {\left( {1 - a} \right){h_2} - \beta } \right] - \alpha \delta_2} \right\}{(\delta _2-\tau)}}}{{\delta _2\left[ {\left( {1 - b} \right)\left( {1 - {h_2}} \right)\left( {{h_2} - \beta } \right) - \alpha \left( {{\delta _2} - \tau } \right)} \right]}}} \right|\Lambda\left( t \right) \end{array}

    with \Lambda\left(t \right) = \exp \left\{ {\int_{{0^ + }}^{{T_2}} {\left[{\phi_2(t)\left({\beta + 1 - 2\phi_2(t)} \right) -\frac{{{s_1}f\varphi_2(t)}}{{{{\left({1 + f\varphi_2(t)} \right)}^2}}}} \right]} dt} \right\} .

    Proof. For system (2.3), there are

    {f_1}(x, y) = x(1 - x)(x - \beta ) - \alpha xy, {f_2}(x, y) = \frac{{\gamma xy}}{{1 + \delta x}} + \left( {\frac{{{s_1}}}{{1 + fy}} - {d_1}} \right)y,

    and {\omega_1}\left({x, y} \right) = x - {h_1}, {I_{11}}\left({x, y} \right) = \eta, {I_{21}}\left({x, y} \right) = 0 , So it can be concluded that

    \frac{{\partial {f_1}\left( {x, y} \right)}}{{\partial x}} = \frac{{\dot x}}{x} + x\left( {\beta + 1 - 2x} \right), \frac{{\partial {f_2}\left( {x, y} \right)}}{{\partial x}} = \frac{{\dot y}}{y} - \frac{{{s_1}fy}}{{{{\left( {1 + fy} \right)}^2}}},
    \frac{{\partial {\omega_1}\left( {x, y} \right)}}{{\partial x}} = 1, \frac{{\partial {\omega_1}\left( {x, y} \right)}}{{\partial y}} = \frac{{\partial {I_{11}}\left( {x, y} \right)}}{{\partial x}} = \frac{{\partial {I_{11}}\left( {x, y} \right)}}{{\partial y}} = \frac{{\partial {I_{21}}\left( {x, y} \right)}}{{\partial x}} = \frac{{\partial {I_{21}}\left( {x, y} \right)}}{{\partial y}} = 0.

    Denote {f_{{1^ + }}}\triangleq {f_1}\left({\phi_1\left({T_1^ + } \right), \varphi_1\left({T_1^ + } \right)} \right), {f_{{2^ + }}}\triangleq {f_2}\left({\phi_1\left({T_1^ + } \right), \varphi_1 \left({T_1^ + } \right)} \right) . Then by Lemma 1, there are

    \begin{align*} {\kappa_1} & = \frac{{\left( { \frac{{\partial {I_{21}}}}{{\partial x}} \cdot \frac{{\partial {\omega_1}}}{{\partial x}} - \frac{{\partial {\beta_1}}}{{\partial x}} \cdot \frac{{\partial {\omega_1}}}{{\partial y}} + \frac{{\partial {\omega_1}}}{{\partial x}}} \right){f_{{1^ + }}} + \left( { \frac{{\partial {I_{11}}}}{{\partial x}} \cdot \frac{{\partial {\omega_1}}}{{\partial y}} - \frac{{\partial {I_{11}}}}{{\partial y}} \cdot \frac{{\partial {\omega_1}}}{{\partial x}} + \frac{{\partial {\omega_1}}}{{\partial y}}} \right){f_{{2^ + }}}}}{{ \frac{{\partial {\omega_1}}}{{\partial x}}{f_1} + \frac{{\partial {\omega_1}}}{{\partial y}}{f_2}}}\\ & = \frac{{\left( {{h_1} + \eta } \right)\left[ {1 - \left( {{h_1} + \eta } \right)} \right]\left[ {\left( {{h_1} + \eta } \right) - \beta } \right] - \alpha \left( {{h_1} + \eta } \right){\delta _1}}}{{{h_1}\left( {1 - {h_1}} \right)\left( {{h_1} - \beta } \right) - \alpha {h_1}{\delta _1}}} \end{align*}

    and

    \begin{align*} \mu_1& = |{\kappa_1}|\exp \left[ {\int_{{0^ + }}^{{T_1}} {\left( {\frac{{\partial {f_1}}}{{\partial x}} + \frac{{\partial {f_2}}}{{\partial y}}} \right)dt} } \right]\\ & = \left|\frac{{\left[ {1 - \left( {{h_1} + \eta } \right)} \right]\left[ {\left( {{h_1} + \eta } \right) - \beta } \right] - \alpha {\delta _1}}}{{\left( {1 - {h_1}} \right)\left( {{h_1} - \beta } \right) - \alpha {\delta _1}}}\right|\exp \left\{ {\int_{{0^ + }}^{{T_1}} {\left[ {\phi_1(t)\left( {\beta + 1 - 2\phi_1(t)} \right) - \frac{{{s_1}f\varphi_1(t)}}{{{{\left( {1 + f\varphi_1(t)} \right)}^2}}}} \right]dt} } \right\}. \end{align*}

    If \mu_1 < 1 , the order-1 periodic solution \mathbf{z}_1(t) = (\phi_1(t), \varphi_1(t)) (k-1)T_1\leq t\leq kT_1 is orbitally asymptotically stable.

    Similarly, for the order-1 periodic solution \mathbf{z}_2(t) = (\phi_2(t), \varphi_2(t)) (k-1)T_2\leq t\leq kT_2 , there is

    \mu_2 = \left| {\frac{{\left\{ {\left[ {1 - \left( {1 - a} \right){h_2}} \right]\left[ {\left( {1 - a} \right){h_2} - \beta } \right] - \alpha \delta_2} \right\}{(\delta _2-\tau)}}}{{\delta _2\left[ {\left( {1 - b} \right)\left( {1 - {h_2}} \right)\left( {{h_2} - \beta } \right) - \alpha \left( {{\delta _2} - \tau } \right)} \right]}}} \right|\Lambda\left( t \right),

    where \Lambda\left(t \right) \triangleq \exp \left\{ {\int_{{0^ + }}^{{T_2}} {\left[{\phi_2(t)\left({\beta + 1 - 2\phi_2(t)} \right) - \frac{{{s_1}f\varphi_2(t)}}{{{{\left({1 + f\varphi_2(t)} \right)}^2}}}} \right]} dt} \right\} . By Lemma 1, if \mu_2 < 1 , the order-1 periodic solution \mathbf{z}_2(t) = (\phi_2(t), \varphi_2(t)) (k-1)T_2\leq t\leq kT_2 is orbitally asymptotically stable. The proof is completed.

    Here only the case for {h_1} < \left({1 - a} \right){h_2} < {x^*_{1}} < {h_1} + \eta < {h_2} < {x^*_{2}} is considered. Let {G_1}\left({{h_1} + \eta, {y_{{G_1}}}} \right) be the intersection point between {l_1} and {\mathcal{N}_1} , {G_2}\left({{h_1}, {y_{{G_2}}}} \right) be the intersection point between {l_2} and {\mathcal{M}_1} . The trajectory that starts with {G_1} intersects {\mathcal{M}_2} at {B_1}\left({{h_2}, {y_{{B_1}}}} \right) . Define \overline{\tau}_3\triangleq y_{G_1}-(1-b)y_{B_1} .

    Theorem 8. For the model parameters with any one of C_1) - C_3) , \rho > 0 and D-4) {h_1} < \left({1 - a} \right){h_2} < {x^*_{1}} < {h_1} + \eta < {h_2} < {x^*_{2}} , if by_{G_2}\leq \tau\leq \min\{y_{G_2}, \overline{\tau}_3\} holds, an order-2 periodic solution exists in system (2.3).

    Proof. For {G_1}\in \mathcal{N}_1 , the trajectory starting from {G_1} intersects with {\mathcal{M}_2} at {B_1} . Then {B_1} is mapped to {B''_1}\in{\mathcal{N}_2} , and next intersects with \mathcal{M}_1 at \hat{B}_1 , and then \hat{B}_1 is mapped to G^{+}_1 under impulse effect. Since \tau < \overline{\tau}_3 , then {y_{B''_1}} = \left({1 - b} \right){y_{{B_1}}} + \tau < {y_{{G_1}}} . Since g_1(B''_1) < 0 and g_2(B''_1) < 0 , it has {y_{{{\hat B}_1}}} < {y_{{{B''}_1}}} , then {y_{G_1^ + }} = y_{\hat{B}_1} < y_{B''_1} = \left({1 - b} \right){y_{{B_1}}} + \tau < {y_{{G_1}}} , i.e., f^{\text{II}}_{\text{SOR}1}({G_1}) = {y_{G_1^ + }} - {y_{{G_1}}} < 0 .

    On the other hand, since \tau\leq y_{G_2} , then G_2\in \mathcal{N}_1 . The trajectory starting from A_1\in \mathcal{N}_1 with y_{A_1} = y_{G_2} intersects with \mathcal{M}_2 at A^{-}_1 , and then A^{-}_1 is mapped to A^{+}_1 . Then it intersects with \mathcal{M}_1 at A^{+-}_1 , and next A^{+-}_1 is mapped to A^{++}_1 . Since y_{A^{+}_1} = (1-b)y_{A^{-}_1}+\tau > (1-b)y_{G_2}+\tau > y_{G_2} , so f^{\text{II}}_{\text{SOR}1}({A_1}) = y_{A^{++}_1} - y_{A_1} = y_{A^{+-}_1}-y_{G_2} > 0 .

    The continuity of f^{\text{II}}_{\text{SOR}_1} implies that a point S\in \overline{A_1G_1}\subset \mathcal{N}_1 exists so that f^{\text{II}}_{\text{SOR}_1}(S) = 0 , i.e., an order-2 periodic solution exists (Figure 4(a)).

    Figure 4.  Schematic diagram of the trajectory trend of the system's (2.3) for case D-4).

    Similarly, a point S'\in \overline{A_2B2}\subset \mathcal{N}_2 exists so that f^{\text{II}}_{\text{SOR}_2}(S') = 0 , i.e., an order-2 periodic solution exists (Figure 4(b)).

    Let \mathbf{z}_{3}(t) = (\phi_3(t), \varphi_3(t)) (k-1)(T_1+T_2)\leq t\leq k(T_1+T_2) , k\in \mathbb{N} be the order-2 periodic solution. Denote

    \phi_3(0) = {h_1} + \eta, \phi_3(T_1) = {h_2}, \phi_3(T_1^+) = (1 - a)h_2, \phi_3(T_1 + T_2) = {h_1}, \phi_3((T_1 + T_2)^+) = {h_1} + \eta,
    \varphi_3(0) = {\delta_3}, \varphi_3(T_1) = \delta_4, \varphi_3(T_1^+) = (1-b)\delta_4 + \tau, \varphi_3(T_1 +T_2) = \delta_3, \varphi_3((T_1+T_2)^+) = \delta_3.

    Define

    \Theta_0\triangleq \frac{\delta_4}{(1-b)\delta_4+\tau}\frac{\gamma_1\gamma_3}{\gamma_2\gamma_4},

    where

    \begin{array}{l} \gamma_1\triangleq (1-h_1-\eta)(h_1+\eta-\beta)-\alpha\delta_3, \gamma_2\triangleq (1-h_2)(h_2-\beta)-\alpha\delta_4, \\ \gamma_3\triangleq (1-(1-a)h_2)((1-a)h_2-\beta)-\alpha((1-b)\delta_4+\tau), \gamma_4\triangleq (1-h_1)(h_1-\beta)-\alpha\delta_3. \end{array}

    Theorem 9. For the model parameters with any one of C_1) C_3) , \rho > 0 and D-4) {h_1} < \left({1 - a} \right){h_2} < {x^*_{1}} < {h_1} + \eta < {h_2} < {x^*_{2}} , and by_{G_2}\leq \tau\leq \min\{y_{G_2}, \overline{\tau}_3\} , if

    \mu_3 = \Theta_0 \exp \left\{ {\int_{{0^ + }}^{{T_1} + {T_2}} {\left[ {\phi_3(t)\left( {\beta + 1 - 2\phi_3(t)} \right) -\frac{{{s_1}f\varphi_3(t)}}{{{{\left( {1 + f\varphi_3(t)} \right)}^2}}}} \right]} dt} \right\} < 1,

    then \mathbf{z}_{3}(t) = (\phi_3(t), \varphi_3(t)) (k-1)(T_1+T_2)\leq t\leq k(T_1+T_2) is orbitally asymptotically stable.

    Proof. For convenience, denote the intersection point between \mathbf{z}_{3}(t) = (\phi_3(t), \varphi_3(t)) and \mathcal{N}_1 ( \mathcal{M}_2 , \mathcal{N}_2 , \mathcal{M}_1 ) by L_1(h_1+\eta, \delta_3) ( L_2(h_2, \delta_4) , L_3((1-a)h_2, (1-b)\delta_4+\tau) , L_4(h_1, \delta_3) ). Then, according to analogue of Poincaré Criterion, there are

    {\kappa_1} = \frac{f_1(L_3)}{f_1(L_2)} = \frac{(1-a)h_3[(1-(1-a)h_2)((1-a)h_2-\beta)-\alpha((1-b)\delta_4+\tau)]}{h_2[(1-h_2)(h_2-\beta)-\alpha\delta_4]},
    {\kappa_2} = \frac{f_1(L_1)}{f_1(L_4)} = \frac{(h_1+\eta)[(1-h_1-\eta)(h_1+\eta-\beta)-\alpha\delta_3]}{h_1[(1-h_1)(h_1-\beta)-\alpha\delta_3}

    and

    \begin{array}{ll} \int_{{0^ + }}^{{T_1} + {T_2}} {\left( {\frac{{\partial {f_1}}}{{\partial x}} + \frac{{\partial {f_2}}}{{\partial y}}} \right)dt} & = \ln \left( {\frac{{{h_2}}}{{{h_1} + \eta }}} \right) + \ln \left( {\frac{{{\delta_4}}}{{{\delta_{3}}}}} \right) + \ln \left( {\frac{{{h_1}}}{{\left( {1 - a} \right){h_2}}}} \right) + \ln\left(\frac{\delta_3}{(1-b)\delta_4+\tau}\right)\\ & + \int_{{0^ + }}^{{T_1} + {T_2}} {\left[ {\phi_3(t)\left( {\beta + 1 - 2\phi_3(t)} \right) - \frac{{{s_1}f\varphi_3(t)}}{{{{\left( {1 + f\varphi_3(t)} \right)}^2}}}} \right]dt}. \end{array}

    Then

    \begin{array}{ll} \mu_3& = |\kappa_1\kappa_2|\exp\left(\int_{{0^ + }}^{{T_1} + {T_2}} {\left( {\frac{{\partial {f_1}}}{{\partial x}} + \frac{{\partial {f_2}}}{{\partial y}}} \right)dt}\right)\\ & = \frac{\delta_4}{(1-b)\delta_4+\tau}\frac{\gamma_1\gamma_3}{\gamma_2\gamma_4}\exp\left(\int_{{0^ + }}^{{T_1} + {T_2}} {\left[ {\phi_3(t)\left( {\beta + 1 - 2\phi_3(t)} \right) - \frac{{{s_1}f\varphi_3(t)}}{{{{\left( {1 + f\varphi_3(t)} \right)}^2}}}} \right]dt}\right). \end{array}

    If \mu < 1 , then \mathbf{z}_{3}(t) = (\phi_3(t), \varphi_3(t)) (k-1)(T_1+T_2)\leq t\leq k(T_1+T_2) is orbitally asymptotically stable.

    For system (2.2) with model parameters \alpha = 0.09 , \beta = 0.15 , \gamma = 0.6 , \delta = 0.9 , {s_1} = 1.12 , {d_1} = 1.1 , f = 0.21 , there is \rho < 0 , then two interior equilibria exist in the system, and the phase diagram is presented in Figure 5(a); For the model parameters \alpha = 0.219 , \beta = 0.051 , \gamma = 0.6 , \delta = 0.9 , {s_1} = 1.12 , {d_1} = 1.1 , f = 0.21 , there is \rho = 0 , then a unique equilibrium exists in system (2.2), which is a sharp point. The phase diagram of the system (2.2) for such a case is presented in Figure 5(b). Moreover, system (2.2) undergoes a Bogdanov-Takens bifurcation of codimension two in a very small neighborhood of the unique interior equilibrium.

    Figure 5.  Illustration of the trajectory trend in system (2.2) for different parameters.

    While for system (2.2) with model parameters \alpha = 0.09 , \beta = 0.29 , \gamma = 0.6 , \delta = 0.9 , {s_1} = 1.12 , {d_1} = 1.1 , f = 0.21 , there is also \rho = 0 , in such case a unique equilibrium exists in system (2.2), which is a saddle node, and the phase diagram is presented in Figure 5(c). System (2.2) undergoes a saddle-node bifurcation of codimension one in a very small neighborhood of the interior equilibrium.

    For system (2.2) with model parameters \alpha = 0.09 , \beta = 0.35, \gamma = 0.6 , \delta = 0.9 , {s_1} = 1.12 , {d_1} = 1.1 , f = 0.21 , there is \rho > 0 . In such a case, system (2.2) does not have interior equilibrium, and the phase diagram of the system (2.2) is presented in Figure 5(d).

    For system (2.2) with parameters \beta = 0.051 , \gamma = 0.6 , \delta = 0.9 , {s_1} = 1.12 , {d_1} = 1.1 , f = 0.21 , the bifurcation diagrams of the residual dimension 1 for the system (2.2) are shown in Figure 6(a), (b) when \alpha is selected as the bifurcation parameter. The result shows that for larger \alpha , the interior equilibrium does not exist. When \alpha decreases to \alpha = \alpha_0 , a unique equilibrium exists in the system. As \alpha decreases below \alpha_0 , the system undergoes a saddle-node bifurcation at the interior equilibrium {E^*} , giving rise to two interior equilibria {E^*_{{1}}} and {E^*_{{2}}} .

    Figure 6.  Schematic diagram of system (2.3) bifurcation when \alpha is selected as a key parameter.

    For system (2.3) with model parameters \alpha = 0.1 , \beta = 0.15 , \gamma = 0.6 , \delta = 0.9 , {s_1} = 1.12 , {d_1} = 1.1 , f = 0.21 , the control parameters: {h_1} = 0.07 , \eta = 0.1 , {h_2} = 0.6 , a = 0.35 , b = 0.38 , \tau = 0.51 , an order-1 periodic solution can be formed in both \Omega_1 and \Omega_2 (Theorem 6), as presented in Figure 7(a); For model parameters \alpha = 0.1 , \beta = 0.15 , \gamma = 0.6 , \delta = 0.9 , {s_1} = 1.12 , {d_1} = 1.1 , f = 0.21 and control parameters {h_1} = 0.07 , \eta = 0.128 , {h_2} = 0.44 , a = 0.7 , b = 0.08 , \tau = 0.51 , an order-1 periodic solution can be formed in \Omega_1 (Theorem 6), as presented in Figure 7(b); while for control parameters {h_1} = 0.14 , \eta = 0.16 , {h_2} = 0.6 , a = 0.25 , b = 0.38 , \tau = 0.51 , an order-1 periodic solution can be formed (Theorem 6), as presented in Figure 7(c).

    Figure 7.  Illustration of the order 1 periodic solution. Schematic diagram of system (2.3) under different parameters.

    For system (2.3) with given model parameters \alpha = 0.09 , \beta = 0.215 , \gamma = 0.6 , \delta = 0.9 , {s_1} = 1.12 , {d_1} = 1.1 , f = 0.21 and control parameters {h_1} = 0.07 , \eta = 0.1 , {h_2} = 0.4 , a = 0.65 , b = 2 , \tau = 3 , an order-2 periodic solution can be formed in system (2.3) (Theorem 3.2), as presented in Figure 8(a); while for the control parameters {h_1} = 0.07 , \eta = 0.4 , {h_2} = 0.68 , a = 0.65 , b = 0.85 , \tau = 0.5 , a different order-2 periodic solution can be formed in system (2.3) (Theorem 8), as presented in Figure 8(b).

    Figure 8.  Illustration of the order 2 periodic solution. Schematic diagram of system (2.3) under different parameters.

    Considering that the Allee effect is an important mechanism in ecosystems and a realistic description of the interaction between species, we presented a model of prey-predator system with prey's Allee effect and generalist predator in the context of fish resources (Models (2.1) or (2.2)). We investigated the dynamic properties of Model (2.2) such as the type and stability of the boundary equilibria as well as the existence and stability of the interior equilibrium in detail (Theorems 1–3, Figure 5).

    To show the influence of the parameters on the dynamics of Model (2.2), we analyzed the bifurcations in the predation system by selecting the capture rate of prey by predator and Allee threshold as key parameters. We showed that Model (2.2) will undergo a saddle-node bifurcation as changing of the capture rate \alpha (Figure 6), and undergo a Bogdanov-Takens bifurcation of codimension at least 2 and 3 as changing of (\alpha, \beta) .

    To achieve sustainable and efficient exploitation of fish stocks, we adopted a bilateral intervention strategy, i.e., to avoid the distinction of prey fish caused by the Allee effect, releasing juvenile prey fish is adopted at a lower-level x = h_1 ; while for economic purposes, capturing both prey and predator fish is adopted at a higher-level x = h_2 . We obtained the conditions for the existence and stability of the order-1 periodic solution (Theorems 6, 7, Figure 7) and order-2 periodic solution (Theorems 8, 9, Figure 8) of the control system (2.3). The results showed that the extinction can be prevented by control even when the prey density is low, while in the case of the prey density increasing to a certain extent, fishing activities can be taken in a periodic way (periodic solution) to obtain the fish resources. Therefore, as long as the fish stocks are properly managed, the number of fish stocks can be controlled within an appropriate range, and the sustainable development and utilization of biological resources can be realized.

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    The work was supported by the National Natural Science Foundation of China (No. 11401068).

    The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.



    [1] A. Charnes, W. W. Cooper, E. Rhodes, Measuring the efficiency of decision making units, Eur. J. Oper. Res., 2 (1978), 429–444. https://doi.org/10.1016/0377-2217(78)90138-8 doi: 10.1016/0377-2217(78)90138-8
    [2] A. Charnes, W. W. Cooper, B. Golany, L. Seiford, J. Stutz, Foundations of data envelopment analysis for Pareto-Koopmans efficient empirical production functions, J. Econometrics, 30 (1985), 91–107. https://doi.org/10.1016/0304-4076(85)90133-2 doi: 10.1016/0304-4076(85)90133-2
    [3] L. A. Zadeh, Fuzzy sets, Inf. Control, 8 (1965), 338–353. https://doi.org/10.1016/S0019-9958(65)90241-X doi: 10.1016/S0019-9958(65)90241-X
    [4] L. Sahoo, An approach for solving fuzzy matrix games using signed distance method, J. Inf. Comput. Sci., 12 (2017), 73–80.
    [5] K. T. Atanassov, Intuitionistic fuzzy sets, Fuzzy Set. Syst., 20 (1986), 87–96. https://doi.org/10.1016/S0165-0114(86)80034-3 doi: 10.1016/S0165-0114(86)80034-3
    [6] R. R. Yager, Pythagorean membership grades in multi-criteria decision making, IEEE T. Fuzzy Syst., 22 (2014), 958–965. https://doi.org/10.1109/TFUZZ.2013.2278989 doi: 10.1109/TFUZZ.2013.2278989
    [7] R. R. Yager, Pythagorean fuzzy subsets, In: 2013 Joint IFSA world congress and NAFIPS annual meeting, 2013, 57–61. https://doi.org/10.1109/IFSA-NAFIPS.2013.6608375
    [8] T. Senapati, R. R. Yager, Fermatean fuzzy sets, J. Amb. Intel. Hum. Comp., 11 (2020), 663–674. https://doi.org/10.1007/s12652-019-01377-0 doi: 10.1007/s12652-019-01377-0
    [9] T. Senapati, R. R. Yager, Fermatean fuzzy weighted averaging/geometric operators and its application in multi-criteria decision making methods, Eng. Appl. Artif. Intel., 85 (2019), 112–121. https://doi.org/10.1016/j.engappai.2019.05.012 doi: 10.1016/j.engappai.2019.05.012
    [10] T. Senapati, R. R. Yager, Some new operations over Fermatean fuzzy numbers and application of Fermatean fuzzy WPM in multiple criteria decision making, Informatica, 30 (2019), 391–412.
    [11] L. Sahoo, Some score functions on Fermatean fuzzy sets and its application to bride selection based on TOPSIS method, Int. J. Fuzzy Syst. Appl., 10 (2021), 18–29. https://doi.org/10.4018/IJFSA.2021070102 doi: 10.4018/IJFSA.2021070102
    [12] L. Sahoo, Similarity measures for Fermatean fuzzy sets and its applications in group decision-making, Decis. Sci. Lett., 11 (2022), 167–180. https://doi.org/10.5267/j.dsl.2021.11.003 doi: 10.5267/j.dsl.2021.11.003
    [13] R. E. Bellman, L. A. Zadeh, Decision making in a fuzzy environment, Manage. Sci., 17 (1970), 141–164. https://doi.org/10.1287/mnsc.17.4.B141 doi: 10.1287/mnsc.17.4.B141
    [14] H. J. Zimmerman, Fuzzy programming and linear programming with several objective functions, Fuzzy Set. Syst., 1 (1978), 45–55. https://doi.org/10.1016/0165-0114(78)90031-3 doi: 10.1016/0165-0114(78)90031-3
    [15] T. Allahviranloo, F. H. Lotfi, M. L. Kiasary, N. A. Kiani, L. A. Zadeh, Solving fully fuzzy linear programming problem by the ranking function, Appl. Math. Sci., 2 (2008), 19–32.
    [16] M. Akram, I. Ullah, S. A. Edalatpanah, T. Allahviranloo, Fully Pythagorean fuzzy linear programming problems with equality constraints, Comput. Appl. Math., 40 (2021), 120. https://doi.org/ 10.1007/s40314-021-01503-9 doi: 10.1007/s40314-021-01503-9
    [17] M. Akram, I. Ullah, T. Allahviranloo, S. A. Edalatpanah, LR-type fully Pythagorean fuzzy linear programming problems with equality constraints, J. Inte. Fuzzy Syst., 41 (2021), 1975–1992. https://doi.org/ 10.3233/JIFS-210655 doi: 10.3233/JIFS-210655
    [18] M. Akram, G. Shahzadi, A. A. H. Ahmadini, Decision-making framework for an effective sanitizer to reduce COVID-19 under Fermatean fuzzy environment, J. Math., 2020 (2020), 3263407. https://doi.org/10.1155/2020/3263407 doi: 10.1155/2020/3263407
    [19] M. Akram, I. Ullah, M. G. Alharbi, Methods for solving LR-type Pythagorean fuzzy linear programming problems with mixed constraints, Math. Probl. Eng., 2021 (2021), 4306058. https://doi.org/10.1155/2021/4306058 doi: 10.1155/2021/4306058
    [20] M. Akram, S. M. U. Shah, M. A. Al-Shamiri, S. A. Edalatpanah, Fractional transportation problem under interval-valued Fermatean fuzzy sets, AIMS Mathematics, 7 (2022), 17327–17348. https://doi.org/ 10.3934/math.2022954 doi: 10.3934/math.2022954
    [21] M. A. Mehmood, M. Akram, M. G. Alharbi, S. Bashir, Solution of fully bipolar fuzzy linear programming models, Math. Probl. Eng., 2021 (2021), 9961891. https://doi.org/10.1155/2021/9961891 doi: 10.1155/2021/9961891
    [22] M. A. Mehmood, M. Akram, M. G. Alharbi, S. Bashir, Optimization of LR-type fully bipolar fuzzy linear programming problems, Math. Probl. Eng., 2021 (2021), 1199336. https://doi.org/10.1155/2021/1199336 doi: 10.1155/2021/1199336
    [23] J. Ahmed, M. G. Alharbi, M. Akram, S. Bashir, A new method to evaluate linear programming problem in bipolar single-valued neutrosophic environment, Comp. Model. Eng., 129 (2021), 881–906. https://doi.org/10.32604/cmes.2021.017222 doi: 10.32604/cmes.2021.017222
    [24] F. L. Hitchcock, The distribution of product from several resources to numerous localities, J. Math. Phys., 20 (1941), 224–230. https://doi.org/10.1002/sapm1941201224 doi: 10.1002/sapm1941201224
    [25] R. D. Banker, A. Charnes, W. W. Cooper, Some models for estimating technical and scale inefficiencies in data envelopment analysis, Manage. Sci., 30 (1984), 1078–1092. https://doi.org/10.1287/mnsc.30.9.1078 doi: 10.1287/mnsc.30.9.1078
    [26] T. Ahn, A. Charnes, W. W. Cooper, Some statistical and DEA evaluations of relative efficiencies of public and private institutions of higher learning, Socio-Econ. Plan. Sci., 22 (1988), 259–269. https://doi.org/10.1016/0038-0121(88)90008-0 doi: 10.1016/0038-0121(88)90008-0
    [27] Y. Roll, W. D. Cook, B. Golany, Controlling factor weights in data envelopment analysis, IIE Trans., 23 (1991), 2–9. https://doi.org/10.1080/07408179108963835 doi: 10.1080/07408179108963835
    [28] J. K. Sengupta, A fuzzy systems approach in data envelopment analysis, Comput. Math. Appl., 24 (1992), 259–266. https://doi.org/10.1016/0898-1221(92)90203-T doi: 10.1016/0898-1221(92)90203-T
    [29] C. Kao, S. T. Liu, Fuzzy efficiency measures in data envelopment analysis, Fuzzy Set. Syst., 113 (2000), 427–437. https://doi.org/10.1016/S0165-0114(98)00137-7 doi: 10.1016/S0165-0114(98)00137-7
    [30] S. Saati, M. Memariani, G. R. Jahanshahloo, Efficiency analysis and ranking of DMUs with fuzzy data, Fuzzy Optim. Decis. Ma., 1 (2002), 255–267. https://doi.org/10.1023/A:1019648512614 doi: 10.1023/A:1019648512614
    [31] S. Lertworasirikul, S. C. Fang, J. A. Joines, H. L. Nuttle, Fuzzy data envelopment analysis (DEA): A possibility approach, Fuzzy Set. Syst., 139 (2003), 379–394. https://doi.org/10.1016/S0165-0114(02)00484-0 doi: 10.1016/S0165-0114(02)00484-0
    [32] A. L. M. Zerafat, S. M. Saati, M. Mokhtaran, An alternative approach to assignment problem with nonhomogeneous costs using common set of weights in DEA, Far East J. Math. Sci., 10 (2003), 29–39.
    [33] W. W. Cooper, L. M. Seiford, K. Tone, Introduction to data envelopment analysis and its uses: With DEA-solver software and references, New York: Springer, 2006.
    [34] P. Zhou, B. W. Ang, K. L. Poh, A survey of data envelopment analysis in energy and environmental studies, Eur. J. Oper. Res., 189 (2008), 1–18. https://doi.org/10.1016/j.ejor.2007.04.042 doi: 10.1016/j.ejor.2007.04.042
    [35] P. Guo, Fuzzy data envelopment analysis and its applications to location problems, Inform. Sci., 179 (2009), 820–829. https://doi.org/10.1016/j.ins.2008.11.003 doi: 10.1016/j.ins.2008.11.003
    [36] F. H. Lotfi, G. R. Jahanshahloo, A. R. Vahidi, A. Dalirian, Efficiency and effectiveness in multi-activity network DEA model with fuzzy data, Appl. Math. Sci., 3 (2009), 2603–2618.
    [37] F. H. Lotfi, G. R. Jahanshahloo, M. Soltanifar, A. Ebrahimnejad, S. M. Mansourzadeh, Relationship between MOLP and DEA based on output-orientated CCR dual model, Expert Syst. Appl., 37 (2010), 4331–4336. https://doi.org/10.1016/j.eswa.2009.11.066 doi: 10.1016/j.eswa.2009.11.066
    [38] S. H. Mousavi-Avval, S. Rafiee, A. Mohammadi, Optimization of energy consumption and input costs for apple production in Iran using data envelopment analysis, Energy, 36 (2011), 909–916. https://doi.org/10.1016/j.energy.2010.12.020 doi: 10.1016/j.energy.2010.12.020
    [39] A. Amirteimoori, An extended transportation problem: A DEA-based approach, Cent. Eur. J. Oper. Res., 19 (2011), 513–521. https://doi.org/10.1007/s10100-010-0140-0 doi: 10.1007/s10100-010-0140-0
    [40] A. Amirteimoori, An extended shortest path problem: A data envelopment analysis approach, Appl. Math. Lett., 25 (2012), 1839–1843. https://doi.org/10.1016/j.aml.2012.02.042 doi: 10.1016/j.aml.2012.02.042
    [41] A. Nabavi-Pelesaraei, R. Abdi, S. Rafiee, H. G. Mobtaker, Optimization of energy required and greenhouse gas emissions analysis for orange producers using data envelopment analysis approach, J. Clean. Prod., 65 (2014), 311–317. https://doi.org/10.1016/j.jclepro.2013.08.019 doi: 10.1016/j.jclepro.2013.08.019
    [42] Z. Zhu, K. Wang, B. Zhang, Applying a network data envelopment analysis model to quantify the eco-efficiency of products: A case study of pesticides, J. Clean. Prod., 69 (2014), 67–73. https://doi.org/10.1016/j.jclepro.2014.01.064 doi: 10.1016/j.jclepro.2014.01.064
    [43] M. Azadi, M. Jafarian, S. R. Farzipoor, S. M. Mirhedayatian, A new fuzzy DEA model for evaluation of efficiency and effectiveness of suppliers in sustainable supply chain management context, Comput. Oper. Res., 54 (2015), 274–285. https://doi.org/10.1016/j.cor.2014.03.002 doi: 10.1016/j.cor.2014.03.002
    [44] G. H. Shirdel, A. Mortezaee, A DEA-based approach for the multi-criteria assignment problem, Croat. Oper. Res. Rev., 6 (2015), 145–154. https://doi.org/10.17535/crorr.2015.0012 doi: 10.17535/crorr.2015.0012
    [45] A. Azar, M. Z. Mahmoudabadi, A. Emrouznejad, A new fuzzy additive model for determining the common set of weights in data envelopment analysis, J. Inte. Fuzzy Syst., 30 (2016), 61–69. https://doi.org/10.3233/IFS-151710 doi: 10.3233/IFS-151710
    [46] A. Mardania, E. Kazimieras, Zavadskasb, Streimikienec, A. Jusoha, M. Khoshnoudia, A comprehensive review of data envelopment analysis (DEA) approach in energy efficiency, Renew. Sust. Energ. Rev., 70 (2017), 1298–1322. https://doi.org/10.1016/j.rser.2016.12.030 doi: 10.1016/j.rser.2016.12.030
    [47] A. Hatami-Marbini, A. Ebrahimnejad, S. Lozano, Fuzzy efficiency measures in data envelopment analysis using lexicographic multiobjective approach, Comput. Ind. Eng., 105 (2017), 362–376. https://doi.org/10.1016/j.cie.2017.01.009 doi: 10.1016/j.cie.2017.01.009
    [48] A. Hatami-Marbini, S. Saati, Efficiency evaluation in two-stage data envelopment analysis under a fuzzy environment: A common weights approach, Appl. Soft Comput., 72 (2018), 156–165. https://doi.org/10.1016/j.asoc.2018.07.057 doi: 10.1016/j.asoc.2018.07.057
    [49] R. M. Rizk-Allaha A. E. Hassanienb, M. Elhoseny, A multi-objective transportation model under neutrosophic environment, Comput. Electr. Eng., 69 (2018), 705–719. https://doi.org/10.1016/j.compeleceng.2018.02.024 doi: 10.1016/j.compeleceng.2018.02.024
    [50] M. Tavana, K. Khalili-Damghani, A new two-stage Stackelberg fuzzy data envelopment analysis model, Measurement, 53 (2014), 277–296. https://doi.org/10.1016/j.measurement.2014.03.030 doi: 10.1016/j.measurement.2014.03.030
    [51] S. A. Edalatpanah, F. Smarandache, Data envelopment analysis for simplified neutrosophic sets, Neutrosophic Sets Sy., 29 (2019), 215–226. https://doi.org/10.5281/zenodo.3514433 doi: 10.5281/zenodo.3514433
    [52] J. Liu, J. Song, Q. Xu, Z. Tao, H. Chen, Group decision making based on DEA cross-efficiency with intuitionistic fuzzy preference relations, Fuzzy Optim. Decis. Ma., 18 (2019), 345–370. https://doi.org/10.1007/s10700-018-9297-0 doi: 10.1007/s10700-018-9297-0
    [53] S. A. Edalatpanah, Data envelopment analysis based on triangular neutrosophic numbers, CAAI T. Intell. Techno., 5 (2020), 94–98. https://doi.org/10.1049/trit.2020.0016 doi: 10.1049/trit.2020.0016
    [54] M. Bagheri, A. Ebrahimnejad, S. Razavyan, F. H. Lotfi, N. Malekmohammadi, Solving the fully fuzzy multi-objective transportation problem based on the common set of weights in DEA, J. Inte. Fuzzy Syst., 39 (2020), 3099–3124. https://doi.org/10.3233/JIFS-191560 doi: 10.3233/JIFS-191560
    [55] M. R. Soltani, S. A. Edalatpanah, F. M. Sobhani, S. E. Najafi, A novel two-stage DEA model in fuzzy environment: Application to industrial workshops performance measurement, Int. J. Comput. Int. Sys., 13 (2020), 1134–1152. https://doi.org/10.2991/ijcis.d.200731.002 doi: 10.2991/ijcis.d.200731.002
    [56] L. Sahoo, A new score function based Fermatean fuzzy transportation problem, Results Control Optim., 1 (2021), 100040. https://doi.org/10.1016/j.rico.2021.100040 doi: 10.1016/j.rico.2021.100040
    [57] S. Ghosh, S. K. Roy, A. Ebrahimnejad, J. L. Verdegay, Multi-objective fully intuitionistic fuzzy fixed-charge solid transportation problem, Complex Intell. Syst., 7 (2021), 1009–1023. https://doi.org/10.1007/s40747-020-00251-3 doi: 10.1007/s40747-020-00251-3
    [58] A. Mondal, S. K. Roy, S. Midya, Intuitionistic fuzzy sustainable multi-objective multi-item multi-choice step fixed-charge solid transportation problem, J. Amb. Intel. Hum. Comp., 2021, 1–25. https://doi.org/10.1007/s12652-021-03554-6 doi: 10.1007/s12652-021-03554-6
    [59] B. K. Giri, S. K. Roy, Neutrosophic multi-objective green four-dimensional fixed-charge transportation problem, Int. J. Mach. Learn. Cyb., 13 (2022), 3089–3112. https://doi.org/10.1007/s13042-022-01582-y doi: 10.1007/s13042-022-01582-y
    [60] S. Ghosh, K-H. Kufer, S. K. Roy, G-W. Weber, Carbon mechanism on sustainable multi-objective solid transportation problem for waste management in Pythagorean hesitant fuzzy environment, Complex Intell. Syst., 8 (2022). https://doi.org/10.1007/s40747-022-00686-w doi: 10.1007/s40747-022-00686-w
    [61] M. Akram, S. M. U. Shah, T. Allahviranloo, A new method to determine the Fermatean fuzzy optimal solution of transportation problems, J. Intell. Fuzzy Syst., 2022. https://doi.org/10.3233/JIFS-221959 doi: 10.3233/JIFS-221959
    [62] M. Bagheri, A. Ebrahimnejad, S. Razavyan, F. H. Lotfi, N. Malekmohammadi, Fuzzy arithmetic DEA approach for fuzzy multi-objective transportation problem, Oper. Res., 22 (2022), 1479–1509. https://doi.org/10.1007/s12351-020-00592-4 doi: 10.1007/s12351-020-00592-4
    [63] Y. M. Wang, Y. Luo, L. Liang, Fuzzy data envelopment analysis based upon fuzzy arithmetic with an application to performance assessment of manufacturing enterprises, Expert Syst. Appl., 36 (2009), 5205–5211. https://doi.org/10.1016/j.eswa.2008.06.102 doi: 10.1016/j.eswa.2008.06.102
    [64] A. Mahmoodirad, T. Allahviranloo, S. Niroomand, A new effective solution method for fully fuzzy transportation problem, Soft Comput., 23 (2019), 4521–4530. https://doi.org/10.1007/s00500-018-3115-z doi: 10.1007/s00500-018-3115-z
    [65] M. Ehrgott, Multi-criteria optimization, Berlin, Heidelberg: Springer, 2005.
  • This article has been cited by:

    1. Xinrui Yan, Yuan Tian, Kaibiao Sun, Dynamic analysis of additional food provided non-smooth pest-natural enemy models based on nonlinear threshold control, 2025, 1598-5865, 10.1007/s12190-024-02318-7
    2. Yuan Tian, Xinlu Tian, Xinrui Yan, Jie Zheng, Kaibiao Sun, Complex dynamics of non-smooth pest-natural enemy Gomportz models with a variable searching rate based on threshold control, 2025, 33, 2688-1594, 26, 10.3934/era.2025002
    3. Xinrui Yan, Yuan Tian, Kaibiao Sun, Hybrid Effects of Cooperative Hunting and Inner Fear on the Dynamics of a Fishery Model With Additional Food Supplement, 2025, 0170-4214, 10.1002/mma.10805
  • Reader Comments
  • © 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(3021) PDF downloads(215) Cited by(41)

Figures and Tables

Figures(4)  /  Tables(8)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog