Stability and optimal harvesting of a predator-prey system combining prey refuge with fuzzy biological parameters

Abstract: In this manuscript, a novel predator-prey system combining prey refuge with fuzzy parameters is formulated. Sufficient conditions for the existence and stability of biological equilibria are derived. The existence of bionomic equilibria is discussed under fuzzy biological parameters. The optimal harvesting policy, by Pontryagin’s maximal principle, is also investigated under imprecise inflation and discount in fuzzy environment. Meticulous numerical simulations are performed to validate our theoretical analysis in detail.


Introduction
Predator-prey system is one of the most important systems in ecology to study the relationship between prey and predator. For the sake of avoiding predation and increasing the survival rate, prey takes different strategies, for example, changing colour, releasing smell, injecting poison, seeking refuges and so on. Since Gause et al. [1] and Smith [2] introduced a quantity involving prey refuge, the study of prey refuge changing the dynamics of predator-prey system has been widely concerned by researchers, see [3][4][5][6] in details. Olivares and Jiliberto [7] put forward that the quantity on prey refuge can be expressed as a constant number or depending on the density of prey. Kar [8], Huang et al. [9] and Tripathi [10] separately studied the impact of refuge on dynamical behavior of predator-prey systems with Holling type II, Holling type III and Beddington-DeAngelis type function responses. Han et al. [11] studied the dynamical behavior of a spatiotemporal predator-prey model with refuge and diffusion including boundedness, permanence, coexisting equilibrium and Turing instability parameter space. Qi and Meng [12] investigated the threshold behavior of a stochastic predator-prey system with refuge and fear effect, and obtained the thresholds determining the extinction and persistence in the mean of prey and predator and the existence of a unique ergodic stationary distribution.
Natural resources are the material basis for human survival and development, among which biological resources are the resources with the most unique development mechanism. They can be continuously renewed via their own reproduction, which is the material embodiment of earth's biological diversity. However, biological resources are also a kind of exhaustible natural resources, and the environmental deterioration or human unreasonable use or excessive development will not only make its quantity and quality decline, but may even lead to the extinction of species. With the demand for food and energy for the development of human beings being rapidly increasing, biological resources are being exploited more and more. Scientific management of the use of marketable biological resources such as fishery, forestry and wildlife is necessary to protect ecosystems. The exploitation of renewable resources used to base on maximum sustainable yield (MSY), but it is impractical to completely ignore the cost operation of resource exploitation. In the face of the insufficiency of MSY, people tried to replace MSY with the optimal sustainable yield (OSY). Lately the research on optimal management of renewable resources has aroused the interest of many scientists and researchers [13][14][15][16][17] and the references therein. Clark [18,19] laid the foundation in the area of work. Kar and Chaudhuri [20] researched a harvesting model on two competing preys and a predator. He and Zhou [21] studied optimal harvesting for a nonlinear hierarchical age-structured population model.
In the study of biological mathematical models, many works focus on the deterministic model, which is mainly based on the law of large numbers and the assumption that the number of biological individuals is sufficiently large, the behaviors of the system will present a relatively stable statistical regularity, so as to be approximately regarded as a deterministic model. However, any species in nature will inevitably be affected by the complexity of the ecosystem itself and the limitations of human cognition, for example, the uncertainty caused by nature environment or human society such as forest fire, flood disaster, debris flow, volcanic eruption, earthquake, climate warming and so on, the uncertainty caused by simplified hypothesis, the uncertainty caused by replacing infinite samples with limited experimental samples in statistical analysis, the uncertainty caused by test results restricted by objective experimental conditions, even the uncertainty caused by some unknown factors. To tackle these uncertainty factors, many researchers used to adopt interval approach [22][23][24][25][26] and stochastic approach [27][28][29][30][31]. In the interval approach the uncertain parameters are characterized by interval-valued functions. In the stochastic approach the uncertain parameters are replaced as stochastic variables with known probability distribution functions. Nevertheless, the interval approach is too simple to do not achieve ideal effect on complex problems when characterizing imprecise parameters. Also a question arising for stochastic approach is that for imprecise parameters, appropriate probability distribution functions will make the model more complex. To overcome these difficulties, we apply fuzzy set theory to depict imprecise parameters. Fuzzy set theory was first come up with by professor Zadeh [32] and he also considered the appliance of fuzzy differential equations is an inherent way to model dynamic systems under possibilistic uncertainty [33]. And Chang and Zadeh [34] first initiated the concept of the fuzzy derivative. And then Kaleva [35] first put forward the concept of differential equations in a fuzzy environment. All derivatives are considered as Hukuhara derivative or generalized derivative in FDE (fuzzy differential equation). Bede et al. [36] studied first order linear fuzzy differential equations under generalized differentiability, and Khastan and Nieto [37] analysed a boundary value problem for second order fuzzy differential equations.
Stefanini and Bede [38] proposed generalization of the Hukuhara difference for compact convex sets, to introduce and study a generalization of the Hukuhara differentiability to the case of interval-valued functions. Puri and Ralescu [39] generalized the concept of Hukuhara differentiability (H-derivative) for set-valued mappings to fuzzy mappings classes. Kaleva [35] developed the theory for FDE by using the H-derivative. Of course, different ways to solve FDE have a huge number of works in this field. Bede and Gal [40] studied the solutions of fuzzy differential equations based on generalized differentiability. Hullermeier [41] proposed an approach to modelling and simulation of uncertain dynamical systems, and showed that all (reasonable) fuzzy functions can be approximated to any degree of accuracy in this way. Allahviranloo and Ahmadi [42] applied Laplace transform method to solve linear FDE. Also the existence and uniqueness of solution for fuzzy initial value problems are investigated by Villamizar-Roa et al. [43]. Different forms of FDEs are come up with in dynamics and good results have been achieved [23,[44][45][46][47]. Bassanezi et al. [48] is the cornerstone of studying the stability of dynamical system by employing fuzzy differential equations. Mizukoshi et al. [49] studied the stability of fuzzy dynamical systems considering initial conditions under fuzziness. Guo et al. [50] established fuzzy impulsive functional differential equations and applied them in logistic model and Gompertz model. The existence and uniqueness of solution for fuzzy functional differential equations were obtained by Lupulescu [51]. Recently, Sadhukhan et al. [52] considered optimal harvesting of a food chain model in fuzzy environment, it is worth mentioning that based on the fuzzy instantaneous annual rate of discount, the optimal harvesting policy is discussed. Pal and Mahapatra [23] drew interval number and fuzzy number into prey and predator interaction's model. The existence and stability of biological and bionomic equilibria are investigated under interval parameters, and the optimal harvesting policy considering instantaneous annual rate of discount in fuzzy environment is studied. Besides, Pal et al. [44] discussed the stability of a fuzzy predator-prey harvesting model with toxicity. And Pal et al. [53] also put forward a predator-prey model with fuzzy optimal harvesting in the environment of toxicity. The former is a major consideration on fuzzy number applied in the parameters of model, the latter mainly considers the inflation and discount rates as fuzzy number, and obtains the fuzzy optimal harvesting. However, for all we know, no one tried to consider both fuzzy biological parameters and fuzzy optimal harvesting in a model. This paper first establishes a predator-prey model with refuge, simultaneously considers both fuzzy biological parameters and fuzzy optimal harvesting into model. The dynamical behaviors covering the existence and stability of biological equilibria are performed. Moreover, there exist bionomic equilibria under fuzzy biological parameters. The optimal harvesting policy is obtained under imprecise inflation and discount in fuzzy environment.
The remaining part of this article is emerged as follows. Firstly, the concepts on fuzzy set and weighted sum method are put forward in Section 2. Section 3 presents the suitable model under fuzzy biological parameters in imprecise environment. The positivity and boundedness of system, as well as the existence and stability of equilibria are investigated in Sections 4 and 5, respectively. In Section 6, we discuss all kinds of bionomic equilibria. The fuzzy optimal harvesting is presented in Section 7.
Some numerical examples verifying the theoretical results are displayed in Section 8. The last section gives a brief summary of the work of this paper.

Preliminaries
This section provides definitions, lemmas and methods that will be used in latter sections.

Concept of fuzzy set
Definition 2.1. [32] Fuzzy set: A fuzzy setÃ in a universe of discourse X is defined as the following set of pairsÃ = {(x, µÃ(x)) : x ∈ X}. The mapping µÃ : X → [0, 1] is called the membership function of the fuzzy setÃ and µÃ is called the membership value or degree of membership of x ∈ X in the fuzzy setÃ.
The arithmetic operations on fuzzy numbers are provided by the following lemmas. Firstly, let us define E as the set of all upper semi-continuous normal convex fuzzy numbers withÃ ∈ E, then the α-level set A α is a bounded closed interval which can be expressed as , whenever non-decreasing sequence α k converges to α ∈ (0, 1]; then, the family [A α L , A α R ], 0 < α ≤ 1, represents the α-level sets of a fuzzy numberÃ in E. Conversely, if [A α L , A α R ], 0 < α ≤ 1 are the α-level sets of a fuzzy numberÃ ∈ E, then the conditions (1) and (2) hold true.

Weighted sum method
In weighted sum method [59], a simple utility function is in the form of w i g i for ith objective, where w i stands for the weight of ith objective. And the total utility function U is expressed as follows where w i > 0 and l i w i = 1 are satisfied.

Model foundation
Motivated by the predator-prey system (3) in Wang et al. [5] as follows where r, K, c, m, d, s and e are positive constants. x(t) and y(t) depict the scale of prey and predator at t with initial densities x 1 (0) > 0 and x 2 (0) > 0. r is the growth rate of prey, K stands for the carrying capacity of prey, c is the predation rate for prey, m shows the proportion of quantity of prey using refuges, d expresses as the death rate of predator and s is the intraspecies competition rate of predator, c is the predation rate, e represents the product of predation rate and conversion rate, q 1 , q 2 and E 1 , E 2 respectively say the catchability coefficients and the harvesting efforts of prey and predator. On account of the concept of fuzzy set, we regard the imprecise parametersr,c,d,s andẽ as triangular fuzzy numbers. Throughout Hukuhara's concept of derivative [39], we present fuzzy differential equations to describe predator-prey interaction with harvesting items as follows Here, prey increases in the absence of predator with fuzzy growth rater while predator decreases with fuzzy mortalityd for lack of prey.s expresses as fuzzy decrease rate of predator attributed to intraspecific competition. The fuzzy quantity of the prey predated per one predator under unit time yields toc(1 − m)x, andẽ represents the fuzzy growth rate of predator attacking on prey.
Applying α-level to cut triangular fuzzy numbersr,c,d,s andẽ, system (3.2) can be expressed as follows By the method of weighted sum, the above system translates into where w 1 and w 2 express as weight satisfying w 1 + w 2 = 1, and w 1 , w 2 ≥ 0. The simplified form of system (3.4) is given by where Of course if we consider imprecise biological parameters with interval parameters rather than fuzzy parameters, then system (3.5) can be translated into system (5) in Wang [5]. Moreover, if we neglect the intraspecific competition rates of prey and predator, as well as refuges, system (3.5) can be simplified to Eqs (7) and (8) in Pal et al. [22].

Positivity and boundedness of system
Throughout the following theorem, in this section, we make sure the positivity and boundedness of system (3.5). Proof. The right hand of system (3.5) satisfies completely continuous and locally Lipschitzian on C, which guarantees that there exists a unique solution (x(t), y(t)) of system (3.5) under initial conditions x(0) > 0 and y(0) > 0 on [0, ξ), where 0 < ξ < +∞. It follows from system (3.5) with x(0) > 0 and y(0) > 0 that which ensures the solution of system (3.5) is positive. On the other side, the first differential equation of system (3.5) provides which can be written as It is easy to see lim sup t→∞ y(t) . Therefore, the solution of system (3.5) is bounded.

Existence and stability of equilibria
We firstly put forward the existence of equilibria of system (3.5). Besides the local stability and global stability of the above equilibria are investigated by the method of Jacobian matrix and Lyapunov function, respectively.

Local stability
On account of the following theorem, the local stability of equilibria S 0 , S 1 , S * for system (3.5) is discussed.
Theorem 5.1. The following conclusions hold: (1) Trivial equilibrium S 0 becomes stable node when A 1 < q 1 E 1 or saddle point when A 1 > q 1 E 1 .
(2) Axial equilibrium S 1 becomes stable node when or saddle point when Interior equilibrium S * is stable node or focus when Proof. The following Jacobian matrix of system (3.5) is expressed as (1) The Jacobian matrix M 0 at S 0 is represented as According to the characteristic equation about the above matrix (2) The Jacobian matrix M 1 at S 1 is depicted as In accordance with the characteristic equation of the above matrix (3) The Jacobian matrix M * at S * is written as The characteristic equation is expressed as It is easy to see the sum of the roots of Eq (5.6) is negative, i.e., λ 1 + λ 2 = − B 2 y * + A 2 x * K < 0 and the product of the roots of Eq (5.6) is positive, i.e., Thus the roots of quadratic Eq (5.6) are negative real number or complex conjugates with negative real parts. Therefore, (x * , y * ) is stable node or focus when the existence conditions

Global stability
We next discuss the global stability of interior equilibrium S * .
Theorem 5.2. If interior equilibrium S * exists, then it is globally asymptotically stable.
Proof. The Lyapunov function is established as follows where l > 0 is a constant. Differentiating both sides of the above equation in regard to time t, one has Since dV dt in some neighborhood of (x * , y * ) is negative semidefinite, the equilibrium (x * , y * ) is globally asymptotically stable.

Bionomic equilibria
As we have discussed in Section 5, the biological equilibrium yields to dx dt ≡ 0 and dy dt ≡ 0. And economic equilibrium is expressed by the total revenue equaling the total cost. Then bionomic equilibrium is a combination of biological equilibrium and economic equilibrium. In this section, we analyse the existence of all bionomic equilibria.
Firstly, c 1 , c 2 denote unit capture cost for prey x and predator y, respectively. And p 1 , p 2 show unit biomass price for prey x and predator y, respectively. Define the profit function as follows N 1 and N 2 separately represent the profit function of prey x and predator y. And the bionomic equilibrium ( The following we consider four cases to excavate the bionomic equilibria of system (3.5).
Case I Boundary bionomic equilibria with no catching species x: E 1∞ = 0. One has ) satisfies Eq (6.3), we have to omit the bionomic equilibrium owing to E 2∞ < 0.
Then solving the second equation for a nonzero x ∞ yields Substituting Eq (6.4) into the third equation of Eq (6.3) gives which is positive provided that holds. Therefore, there exists a boundary bionomic equilibrium (x ∞ , y ∞ , 0, E 2∞ ).
Case II Boundary bionomic equilibria with no catching species y: E 2∞ = 0. It follows from , 0) and we get a nonegative semi-trivial bionomic Then solving the third equation for a nonzero y ∞ yields satisfying B 3 c 1 (1 − m) > B 1 p 1 q 1 . Substituting Eq (6.8) into the second equation of Eq (6.7) gives which is positive provided that holds. Thus, there exists a boundary bionomic equilibrium (x ∞ , y ∞ , E 1∞ , 0). Case III Boundary bionomic equilibria with no catching species x and y: E 1∞ = 0 and E 2∞ = 0. In this case, the cost exceeds revenue for both species x and y, we have to stop harvesting both of them. The bionomic equilibrium is equal to the biological equilibrium with E 1∞ = E 2∞ = 0.
Case IV Positive bionomic equilibrium with effort fishing on both species x and y: E 1∞ > 0 and E 2∞ > 0.
Solving the third equation of Eq (6.2), we obtain x ∞ = c 1 p 1 q 1 , y ∞ = c 2 p 2 q 2 . Then, it follows from the second and third equations of Eq (6.2) that which exist if satisfying the following conditions Therefore, there exists the positive bionomic equilibrium (x ∞ , y ∞ , E 1∞ , E 2∞ ).
Summarizing the above, we put forward the theorem as follows.

Fuzzy optimal harvesting
This section studies the optimal harvesting strategy with fuzzy net discount rate of inflation. Previously, denotek andr the inflation and discount rates and they are considered as fuzzy number essentially. The net discount rate of inflationδ is the difference of the above two fuzzy numbers, that isδ =r −k, and it can be regarded as triangular fuzzy number, that isδ = (δ 1 , δ 2 , δ 3 ). Our target is to maximize the valueJ as follows Using Pontryagin's maximal principle [60] to solve the above fuzzy optimization problem, the harvesting strategies not only guarantee profit maximization, but also maintain an optimal level for species. The control variable E i (t) subject to 0 ≤ E i (t) ≤ E max i , i = 1, 2. According to Maity and Maiti [61], Sadhukhan et al. [52] and Pal and Mahapatra [23], using α to cut triangular fuzzy number δ = (δ 1 , δ 2 , δ 3 ), that is it regarded as an interval number [δ L , δ R ], the optimization problem Maximize Again, by the method of weighted sum, it yields to where w 1 , w 2 ≥ 0 are two weights such that w 1 + w 2 = 1. The Hamiltonian can be written as where λ i = λ i (t)(i = 1, 2) are the adjoint variables. Applying Pontryagin's maximum principle [60], we get the adjoint equations It follows from Eqs (7.4) and (7.5) that Substituting the interior equilibrium into Eq (7.6) yields Now, we neglect λ 1 in Eq (7.7) to get a reduced differential equation for λ 2 where The solution of Eq (7.8) as follows where C i (i = 1, 2) are arbitrary constants and m i (i = 1, 2) are the roots of the auxiliary equations a 0 m 2 + a 1 m + a 2 = 0 and From Eq (7.9), λ 2 is bounded if and only if m i < 0(i = 1, 2) or C i = 0. Since it is difficult to check whether m i < 0, we consider C i = 0. Eq (7.9) is simplified as Similarly, we get It is easy to see the shadow prices λ i e δ L t (i = 1, 2) of the two species satisfy the transversality condition at ∞ [18], i.e., they are bounded as t → ∞. The Hamiltonian in Eq (7.4) must be maximized for Suppose that the optimal equilibrium does not appear at either E i = 0 or E i = E max i , we therefore consider the singular controls that is, (w 1 e −δ R t + w 2 e −δ L t )(p 1 q 1 x − c 1 ) = λ 1 q 1 x, (w 1 e −δ R t + w 2 e −δ L t )(p 2 q 2 y − c 2 ) = λ 2 q 2 y. (7.14) We substitute λ 1 and λ 2 in Eq (7.11) and Eq (7.10) into Eq (7.14) and obtain where To eliminate e −(δ R −δ L )t in Eq (7.17), we consider the following two cases: Case I If P L c 1 w 1 and Q L c 2 w 1 , then Eq (7.17) is equal to Case II If P L = c 1 w 1 or Q L = c 2 w 1 , we have Motivated by the method in [25], firstly transforming Eq (7.14) into the following form Differentiating λ i (i = 1, 2) in Eq (7.19) with respect to t, it derives that (7.20) Substituting λ i (i = 1, 2) in Eq (7.19) into Eq (7.6) yields (7.21) Combined Eqs (7.20) and (7.21) with E 1 and E 2 omitted, one has (w 1 e −δ R t + w 2 e −δ L t ). (7.22) Consider N 1 = (p 1 q 1 x − c 1 )E 1 > 0 and N 2 = (p 2 q 2 y − c 2 )E 2 > 0, we divide both sides of the first equation by that of the second equation of Eq (7.22) Besides, the values of E 1 and E 2 at the interior equilibrium turn to x.

Numerical simulations
In this section, some hypothetical data are taken to illustrate analytical results of the previous sections by mainly using the software MATLAB. Due to the features characterized by simulations should be analysed from qualitative perspective rather than quantitative one, then the parameters in system are not in view of real observed values. We consider three numerical examples to account for the stability of equilibria, the existence of bionomic equilibria as well as fuzzy optimal harvesting. We discuss equilibria S 0 , S 1 and S * of the predator-prey model under different values of α, w 1 and w 2 , which are showed from Table 1 to Table 3. Table 1. The trivial equilibrium S 0 for q 1 = 0.2, q 2 = 0.2, E 1 = 15 and E 2 = 10.
Tables 1-3 demonstrate the trivial equilibrium is always fixed at (0, 0) with different combinations of w 1 , w 2 and α, the values of prey x and predator y stay at zero; for the axial equilibrium, the value of prey x decreases under the same α with increasing w 1 , predator y is invariant in zero; for the interior equilibrium, the value of prey x increases under the same α with increasing w 1 , and predator y decreases under the same α with increasing w 1 .  Table 3. The interior equilibrium S * for q 1 = 0.2, q 2 = 0.2, E 1 = 5 and E 2 = 2.     Tables 1-3 reflect. The initial fluctuate for both species gradually saturate to a steady state level over time.
The phase portrait of prey x and predator y corresponding to stabilities of different equilibria for all kinds of combinations of w 1 , w 2 and α are shown in Figures 4-6, respectively. Meanwhile, the equilibria are stable under different initial conditions.   Dynamical behavior of prey x and predator y under interior steady state level with respect to α for initial condition (x(0), y(0)) = (2.5, 6), and q 1 = 0.2, q 2 = 0.2, E 1 = 5, E 2 = 2 are presented in Figure 7.
The figures in Figure 7 reflect the interior equilibrium changes for different α which is in accordance with the results in Table 3. As α increases, prey x increases and predator y decreases for w 1 = 0, w 2 = 1 and w 1 = 0.2, w 2 = 0.8 and w 1 = 0.4, w 2 = 0.6. However, in the cases w 1 = 0.6, w 2 = 0.4 and w 1 = 0.8, w 2 = 0.2 and w 1 = 1, w 2 = 0, as α increases, prey x decreases and predator y increases. Therefore, the parameter α plays an major part in the stability of interior equilibrium.
In the mean time, the dynamical behavior of prey x and predator y about w 1 and w 2 under initial condition (x(0), y(0)) = (2.5, 6) for α = 0, 0.3, 0.6 and α = 0.9, and q 1 = 0.2, q 2 = 0.2, E 1 = 5, E 2 = 2 are presented in Figure 8. These figures depict the same results as that in Table 3, i.e., for the interior equilibrium, the value of prey x increases under fixed α with increasing w 1 , and predator y decreases under fixed α with increasing w 1 .   The bionomic equilibria for different combinations of w 1 , w 2 and α are emerged in Table 4. It is observed that x ∞ and y ∞ are invariable with respect to w 1 , w 2 and α, and E 1∞ and E 2∞ are decreasing with fixed α and increasing w 1 . As well as, together with the value of α increasing, E 1∞ and E 2∞ decreases for w 1 = 0, w 2 = 1 and w 1 = 0.2, w 2 = 0.8 and w 1 = 0.4, w 2 = 0.6. Nevertheless, in the cases w 1 = 0.6, w 2 = 0.4 and w 1 = 0.8, w 2 = 0.2 and w 1 = 1, w 2 = 0, as the value of α increases, E 1∞ and E 2∞ are increasing.  Table 5. Optimal equilibrium and optimal harvesting effort for w 1 = 0, w 2 = 1, α = 0.  Table 6. Optimal equilibrium and optimal harvesting effort for w 1 = 0, w 2 = 1, α = 0.5. The optimal equilibrium and harvesting effort for different combinations w 1 , w 2 , α and w 1 , w 2 , α are appeared in Tables 5-9, respectively. For different combinations w 1 , w 2 , α and w 1 , w 2 , α , we see the existence of optimal equilibrium.
In Table 5, when α = 0 it is observed that prey x˜δ and predator y˜δ are invariant, optimal harvesting effort of the prey E 1δ decreases and then increases, but the optimal harvesting effort of the predator E 2δ increases and then decreases with increasing w 1 ; when α = 0.3, x˜δ decreases and then increases, y˜δ is invariant, E 1δ decreases and then increases, and E 2δ increases and then decreases with increasing Table 7. Optimal equilibrium and optimal harvesting effort for w 1 = 0, w 2 = 1, α = 0.9.  Table 8. Optimal equilibrium and optimal harvesting effort for w 1 = 0.5, w 2 = 0.5, α = 0.9. w 1 ; when α = 0.6, x˜δ increases and then decreases, y˜δ is invariant, E 1δ fluctuates between 0.0907 and 0.0908, and E 2δ is invariant with increasing w 1 ; when α = 0.9, x˜δ and y˜δ are invariant, E 1δ decreases and then increases, and E 2δ is invariant with increasing w 1 . Besides, when fixed w 1 and w 2 , x˜δ increases, y˜δ decreases, E 1δ increases, and E 2δ decreases with the increase of α . In Table 6, when fixed α , x˜δ, y˜δ, E 1δ and E 2δ are invariant with the increase of w 1 ; when fixed w 1 and w 2 , x˜δ increases, y˜δ decreases, E 1δ increases, and E 2δ increases with the increase of α .
In Table 7, when fixed α , x˜δ, y˜δ, E 1δ and E 2δ are invariant with the increase of w 1 ; when fixed w 1 and w 2 , x˜δ increases, y˜δ is invariant, E 1δ decreases, and E 2δ increases with the increase of α .
In Table 8, when fixed α , x˜δ, y˜δ, E 1δ and E 2δ are invariant with the increase of w 1 ; when fixed w 1 and w 2 , x˜δ increases, y˜δ increases, E 1δ decreases, and E 2δ increases with the increase of α .
In Table 9, when α = 0, x˜δ increases and then decreases and y˜δ is invariant, E 1δ increases and then decreases, but E 2δ decreases and then increases with increasing w 1 ; when α = 0.3 or α = 0.6, we can see that x˜δ increases and then decreases, y˜δ, E 1δ and E 2δ are invariant with increasing w 1 ; when α = 0.9, x˜δ, y˜δ and E 1δ are invariant, E 2δ decreases and then increases with increasing w 1 . When fixed w 1 and w 2 , x˜δ decreases, y˜δ is invariant, E 1δ is almost invariant, and E 2δ increases with the increase of α .
From Tables 5-7, it is easy to see that when fixed w 1 , w 2 and w 1 , w 2 , α , the prey increases, the predator decreases, the optimal harvesting effort of the prey increases and the optimal harvesting effort of the predator decreases with increasing α. Table 9. Optimal equilibrium and optimal harvesting effort for w 1 = 1, w 2 = 0, α = 0.9.  Based on Tables 7-9, we observed that when fixed α and w 1 , w 2 , α , the prey increases, the predator decreases, the optimal harvesting effort of the prey decreases and the optimal harvesting effort of the predator decreases with increasing w 1 .

Conclusions
In the field of biomathematics most harvesting models are on account of the assumption with precise parameters, but any species in nature will inevitably be affected by the complexity of the ecosystem itself and the limitations of people's cognition leading to the uncertainty of the biological parameters. These uncertainties can be roughly classified into three categories: fuzziness, randomness, and unascertainty. In view of this point, we have modified a fuzzy parameters predator-prey model with refuge. Here, the fuzzy parameters are introduced into biotic potential of the prey (r), predator mortality (d), intraspecific competition rate of the predator (s), predation rate (c) and the product of predation rate and conversion rate (e).
We have discussed the positivity and boundedness of system (Theorem 4.1) and investigated the existence and stability of equilibria (Theorems 5.1 and 5.2). The tables and graphs for different equilibria are showed from Table 1 to Table 3 and from Figure 1 to Figure 3, respectively. Figures 4-6 reflect phase portrait of prey and predator about different w 1 , w 2 and α. Dynamical behaviors of prey and predator about α or w 1 and w 2 , separately, are presented in Figures 7 and 8. We have also analysed the existence of bionomic equilibria (Theorem 6.1) and the mathematical results for positive bionomic equilibrium are displayed in Table 4.
In this paper, we consider the fuzzification of biological parameters in model. However, the fuzzification of unit capture cost or unit biomass price for both species seems to be also feasible. Future work will focus on the imprecision of economic parameters to enhance our model. Finally, this paper uses α-cut of triangular fuzzy number to describe imprecise parameters in model, however other common forms of fuzzy number such as normal fuzzy number and trapezoidal fuzzy number are important in fuzzy set theory, and then we are going to apply other forms of fuzzy number depicting imprecise parameters in model.