Almost periodic solutions for a SVIR epidemic model with relapse

1.   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 Ⅱ, Holling type Ⅲ 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.

2.   Preliminaries

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

2.1.   Concept of fuzzy set

Definition 2.1. ^[32] Fuzzy set: A fuzzy set $\tilde{A}$ in a universe of discourse $X$ is defined as the following set of pairs $\tilde{A} = \{(x, \mu_{\tilde{A}}(x)):x\in X\}$. The mapping $\mu_{\tilde{A}}:X\rightarrow [0, 1]$ is called the membership function of the fuzzy set $\tilde{A}$ and $\mu_{\tilde{A}}$ is called the membership value or degree of membership of $x\in X$ in the fuzzy set $\tilde{A}$.

Definition 2.2. ^[54] $\alpha$-cut of fuzzy set: The $\alpha$-cut of a fuzzy set $\tilde{A}$ is a crisp set and it is defined by $A_\alpha = \{x:\mu_{\tilde{A}}(x)\geq\alpha \}, \alpha\in(0, 1]$. For $\alpha = 0$ the support of $\tilde{A}$ is defined as $A_0 = Supp(\tilde{A}) = \overline{\{x\in R, \mu_{\tilde{A}}(x) > 0\}}$.

Definition 2.3. ^[32] Convex fuzzy set: A convex fuzzy set $\tilde{A}$ is a fuzzy set on a continuous universe such that for all $\alpha$, $A_\alpha$ is a convex classical set.

Definition 2.4. ^[55] Fuzzy number: A fuzzy number is a convex fuzzy set with $X = R$.

Definition 2.5. ^[56] Triangular fuzzy number: A triangular fuzzy number (TFN) $\tilde{A}\equiv(a_1, a_2, a_3)$ is fuzzy set of the real line $R$ characterized by the membership function $\mu_{\tilde{A}}:R\rightarrow [0, 1]$ as follows

Therefore, the $\alpha$-cut of triangular fuzzy number is a bounded and closed internal $[A_L(\alpha), A_R(\alpha)]$, where $A_L(\alpha) = {\inf}\{x:\mu_{\tilde{A}}(x)\geq\alpha\} = a_1+\alpha(a_2-a_1)$ and $A_R(\alpha) = {\sup}\{x:\mu_{\tilde{A}}(x)\geq\alpha\} = a_3-\alpha(a_3-a_2)$.

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 $\tilde{A}\in E$, then the $\alpha$-level set $A_\alpha$ is a bounded closed interval which can be expressed as $A_\alpha = [A_L^\alpha, A_R^\alpha]$, where $A_L^\alpha = \inf\{x:\mu_{\tilde{A}}(x)\geq\alpha\}$ and $A_R^\alpha = \sup\{x:\mu_{\tilde{A}}(x)\geq\alpha\}$.

Lemma 2.1. ^[57] If $\tilde{A}, \tilde{B}\in E$, then for $\alpha\in(0, 1]$,

$\rm(1)$ $[\tilde{A}+\tilde{B}]_\alpha = [A_L^\alpha+B_L^\alpha, A_R^\alpha+B_R^\alpha]$,

$\rm(2)$ $[\tilde{A}\cdot\tilde{B}]_\alpha = [\min\{A_L^\alpha B_L^\alpha, A_L^\alpha B_R^\alpha, A_R^\alpha B_L^\alpha, A_R^\alpha B_R^\alpha\}, \max\{A_L^\alpha B_L^\alpha, A_L^\alpha B_R^\alpha, A_R^\alpha B_L^\alpha, A_R^\alpha B_R^\alpha\}]$,

$\rm(3)$ $[\tilde{A}-\tilde{B}]_\alpha = [A_L^\alpha-B_R^\alpha, A_R^\alpha-B_L^\alpha]$.

Lemma 2.2. ^[58] $[A_L^\alpha, A_R^\alpha], 0 < \alpha\leq1$ is a family of nonempty intervals. If

$(1)$ $[A_L^\alpha, A_R^\alpha]\supset[A_L^\beta, A_R^\beta]$ for $0 < \alpha\leq\beta$ and

$(2)$ $[\lim_{k\rightarrow \infty}A_L^{\alpha_k}, \lim_{k\rightarrow \infty}A_R^{\alpha_k}] = [A_L^\alpha, A_R^\alpha]$,

whenever non-decreasing sequence $\alpha_k$ converges to $\alpha\in(0, 1]$; then, the family $[A_L^\alpha, A_R^\alpha]$, $0 < \alpha\leq1$, represents the $\alpha$-level sets of a fuzzy number $\tilde{A}$ in $E$. Conversely, if $[A_L^\alpha, A_R^\alpha]$, $0 < \alpha\leq1$ are the $\alpha$-level sets of a fuzzy number $\tilde{A}\in E$, then the conditions ${\rm(1)}$ and ${\rm(2)}$ hold true.

2.2.   Weighted sum method

In weighted sum method ^[59], a simple utility function is in the form of $w_ig_i$ for $i$th objective, where $w_i$ stands for the weight of $i$th objective. And the total utility function $U$ is expressed as follows

where $w_i > 0$ and $\sum_i^l w_i = 1$ are satisfied.

3.   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 parameters $\tilde{r}, \tilde{c}, \tilde{d}, \tilde{s}$ and $\tilde{e}$ 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 rate $\tilde{r}$ while predator decreases with fuzzy mortality $\tilde{d}$ for lack of prey. $\tilde{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 to $\tilde{c}(1-m)x$, and $\tilde{e}$ represents the fuzzy growth rate of predator attacking on prey.

Applying $\alpha$-level to cut triangular fuzzy numbers $\tilde{r}, \tilde{c}, \tilde{d}, \tilde{s}$ and $\tilde{e}$, 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\geq0$. 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].

4.   Positivity and boundedness of system

Throughout the following theorem, in this section, we make sure the positivity and boundedness of system (3.5).

Theorem 4.1. Any solution of system ${\rm(3.5)}$ is positive and bounded satisfying initial conditions $x(0) > 0$ and $y(0) > 0$ for all $t > 0$.

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, \xi)$, where $0 < \xi < +\infty$. 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 implies $\limsup_{t\rightarrow \infty}x(t)\leq\frac{A_1K}{A_2}\equiv k_1$. From the second differential equation of system (3.5), one has

which can be written as

It is easy to see $\limsup_{t\rightarrow \infty}y(t)\leq\frac{k_1B_3(1-m)}{B_2}$. Therefore, the solution of system (3.5) is bounded.

5.   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.

By a simple calculation, the equilibria of system (3.5) are given by

(1) Trivial equilibrium: $S_0 = (0, 0)$.

(2) Axial equilibrium: $S_1 = (\frac{K(A_1-q_1E_1)}{A_2}, 0)$ exists if $A_1 > q_1E_1$.

(3) Axial equilibrium: $S_2 = (0, -\frac{q_2E_2+B_1}{B_2})$ does not exist since $-\frac{q_2E_2+B_1}{B_2} < 0$.

(4) Interior equilibrium: $S^* = (x^*, y^*)$, where

exists if $B_3(1-m)x^* > B_1+q_2E_2$ and $A_1B_2+A_3(1-m)(B_1+q_2E_2) > B_2q_1E_1$ are satisfied.

5.1.   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:

$\rm(1)$ Trivial equilibrium $S_0$ becomes stable node when $A_1 < q_1E_1$ or saddle point when $A_1 > q_1E_1$.

$\rm(2)$ Axial equilibrium $S_1$ becomes stable node when

or saddle point when

$\rm(3)$ 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

one has $\lambda_1 = A_1-q_1E_1, \ \ \lambda_2 = -B_1-q_2E_2$, hence $A_0$ is stable node when $A_1 < q_1E_1$ or saddle point when $A_1 > q_1E_1$.

(2) The Jacobian matrix $M_1$ at $S_1$ is depicted as

In accordance with the characteristic equation of the above matrix

we have $\lambda_1 = -(A_1-q_1E_1)$ and $\lambda_2 = -(B_1+q_2E_2)+\frac{KB_3(1-m)(A_1-q_1E_1)}{A_2}$. Obviously, $S_1$ exists if $A_1 > q_1E_1$ holds, hence $S_1$ is a stable node when $A_1 > q_1E_1$ and $A_2(B_1+q_2E_2) > KB_3(1-m)(A_1-q_1E_1)$ or saddle point when $A_1 > q_1E_1$ and $A_2(B_1+q_2E_2) < KB_3(1-m)(A_1-q_1E_1)$.

(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., $\lambda_1+\lambda_2 = -\Bigl(B_2y^*+\frac{A_2x^*}{K}\Bigr) < 0$ and the product of the roots of Eq (5.6) is positive, i.e., $\lambda_1\lambda_2 = A_3B_3(1-m)^2x^*y^* > 0$. 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 $B_3(1-m)x^* > B_1+q_2E_2$ and $A_1B_2+A_3(1-m)(B_1+q_2E_2) > B_2q_1E_1$ are satisfied.

5.2.   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

Choose $l = \frac{A_3}{B_3}$, it yields to

Since $\frac{dV}{dt}$ in some neighborhood of $(x^*, y^*)$ is negative semidefinite, the equilibrium $(x^*, y^*)$ is globally asymptotically stable.

6.   Bionomic equilibria

As we have discussed in Section 5, the biological equilibrium yields to $\frac{dx}{dt}\equiv0$ and $\frac{dy}{dt}\equiv0$. 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

where $N_1 = (p_1q_1x-c_1)E_1, N_2 = (p_2q_2y-c_2)E_2$. $N_1$ and $N_2$ separately represent the profit function of prey $x$ and predator $y$. And the bionomic equilibrium $(x_\infty, y_\infty, E_{1\infty}, E_{2\infty})$ is governed by

The following we consider four cases to excavate the bionomic equilibria of system (3.5).

Case Ⅰ Boundary bionomic equilibria with no catching species $x$: $E_{1\infty} = 0$. One has

Clearly, $(x_\infty, y_\infty, E_{1\infty}, E_{2\infty}) = (0, \frac{c_2}{p_2q_2}, 0, -\frac{p_2q_2B_1+B_2c_2}{p_2q_2^2})$ satisfies Eq (6.3), we have to omit the bionomic equilibrium owing to $E_{2\infty} < 0$.

Then solving the second equation for a nonzero $x_\infty$ yields

which exists if satisfying $A_1p_2q_2 > A_3(1-m)c_2$. 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_\infty, y_\infty, 0, E_{2\infty})$.

Case Ⅱ Boundary bionomic equilibria with no catching species $y$: $E_{2\infty} = 0$. It follows from

that $(x_\infty, y_\infty, E_{1\infty}, E_{2\infty}) = (\frac{c_1}{p_1q_1}, 0, \frac{A_1p_1q_1K-A_2c_1}{p_1q_1^2K}, 0)$ and we get a nonegative semi-trivial bionomic equilibrium $(x_\infty, 0, E_{1\infty}, 0)$ with condition $A_1p_1q_1K > A_2c_1$.

Then solving the third equation for a nonzero $y_\infty$ yields

satisfying $B_3c_1(1-m) > B_1p_1q_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_\infty, y_\infty, E_{1\infty}, 0)$.

Case Ⅲ Boundary bionomic equilibria with no catching species $x$ and $y$: $E_{1\infty} = 0$ and $E_{2\infty} = 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\infty} = E_{2\infty} = 0$.

Case Ⅳ Positive bionomic equilibrium with effort fishing on both species $x$ and $y$: $E_{1\infty} > 0$ and $E_{2\infty} > 0$.

Solving the third equation of Eq (6.2), we obtain $x_\infty = \frac{c_1}{p_1q_1}, y_\infty = \frac{c_2}{p_2q_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_\infty, y_\infty, E_{1\infty}, E_{2\infty})$.

Summarizing the above, we put forward the theorem as follows.

Theorem 6.1. The following conclusions hold:

$\rm(1)$ The boundary bionomic equilibrium $(x_\infty, y_\infty, 0, E_{2\infty})$ exists if $A_1p_2q_2 > A_3(1-m)c_2$ and Eq $(6.6)$ are satisfied, where $y_\infty = \frac{c_2}{p_2q_2}$ and $x_\infty$ and $E_{2\infty}$ are defined in ${Eqs}$ $\rm(6.4)$ and ${\rm(6.5)}$, respectively.

$\rm(2)$ The boundary bionomic equilibrium $(x_\infty, 0, E_{1\infty}, 0)$ exists if $A_1p_1q_1K>A_2c_1$ holds, where $x_\infty = \frac{c_1}{p_1q_1}$ and $E_{1\infty} = \frac{A_1p_1q_1K-A_2c_1}{p_1q_1^2K}$.

$\rm(3)$ The boundary bionomic equilibrium $(x_\infty, y_\infty, E_{1\infty}, 0)$ exists if $B_3c_1(1-m)>B_1p_1q_1$ and ${Eq}$ $(6.10)$ are satisfied, where $x_\infty = \frac{c_1}{p_1q_1}$ and $y_\infty$ and $E_{1\infty}$ are defined in ${Eqs}$ ${\rm(6.8)}$ and ${\rm(6.9)}$, respectively.

$\rm(4)$ The boundary bionomic equilibrium $(x_\infty, y_\infty, 0, 0)$ is equal to the biological equilibrium with $E_{1\infty} = E_{2\infty} = 0$ described in Section ${\rm5}$.

$\rm(5)$ The positive bionomic equilibrium $(x_\infty, y_\infty, E_{1\infty}, E_{2\infty})$ exists if ${Eq}$ ${\rm(6.12)}$ holds, where $x_\infty = \frac{c_1}{p_1q_1}$, $y_\infty = \frac{c_2}{p_2q_2}$ and $E_{1\infty}$ and $E_{2\infty}$ are defined in ${Eq}$ ${\rm(6.11)}$.

7.   Fuzzy optimal harvesting

This section studies the optimal harvesting strategy with fuzzy net discount rate of inflation. Previously, denote $\tilde{k}$ and $\tilde{r}$ the inflation and discount rates and they are considered as fuzzy number essentially. The net discount rate of inflation $\tilde{\delta}$ is the difference of the above two fuzzy numbers, that is $\tilde{\delta} = \tilde{r}-\tilde{k}$, and it can be regarded as triangular fuzzy number, that is $\tilde{\delta} = (\delta_1, \delta_2, \delta_3)$. Our target is to maximize the value $\tilde{J}$ 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\leq E_i(t)\leq E_i^{\max}, i = 1, 2$. According to Maity and Maiti ^[61], Sadhukhan et al. ^[52] and Pal and Mahapatra ^[23], using $\alpha'$ to cut triangular fuzzy number $\tilde{\delta} = (\delta_1, \delta_2, \delta_3)$, that is it regarded as an interval number $[\delta_L, \delta_R]$, the optimization problem Maximize $\tilde{J}(E_1, E_2)$ is equal to Maximize $[J_L(E_1, E_2), J_R(E_1, E_2)]$, where

Again, by the method of weighted sum, it yields to

where $w_1', w_2'\geq0$ are two weights such that $w_1'+w_2' = 1$. The Hamiltonian can be written as

where $\lambda_i = \lambda_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 $\lambda_1$ in Eq (7.7) to get a reduced differential equation for $\lambda_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_0m^2+a_1m+a_2 = 0$ and

From Eq (7.9), $\lambda_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

where

It is easy to see the shadow prices $\lambda_ie^{\delta_Lt}(i = 1, 2)$ of the two species satisfy the transversality condition at $\infty$ ^[18], i.e., they are bounded as $t\rightarrow \infty$. The Hamiltonian in Eq (7.4) must be maximized for $E_i\in[0, E_i^{\max}]$. Suppose that the optimal equilibrium does not appear at either $E_i = 0$ or $E_i = E_i^{\max}$, we therefore consider the singular controls

that is,

We substitute $\lambda_1$ and $\lambda_2$ in Eq (7.11) and Eq (7.10) into Eq (7.14) and obtain

where

Equation (7.15) is equivalent to the following equations

To eliminate $e^{-(\delta_R-\delta_L)t}$ in Eq (7.17), we consider the following two cases:

Case Ⅰ If $P_L\neq c_1w_1'$ and $Q_L\neq c_2w_1'$, then Eq (7.17) is equal to

Case Ⅱ If $P_L = c_1w_1'$ or $Q_L = c_2w_1'$, we have

Motivated by the method in ^[25], firstly transforming Eq (7.14) into the following form

Differentiating $\lambda_i(i = 1, 2)$ in Eq (7.19) with respect to $t$, it derives that

Substituting $\lambda_i(i = 1, 2)$ in Eq (7.19) into Eq (7.6) yields

Combined Eqs (7.20) and (7.21) with $E_1$ and $E_2$ omitted, one has

Consider $N_1 = (p_1q_1x-c_1)E_1 > 0$ and $N_2 = (p_2q_2y-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

Solving Eqs (7.18) (or $(7.18)'$), (7.23) and (7.24) to get the optimal equilibrium solution $x = x_{\tilde{\delta}}, y = y_{\tilde{\delta}}$ and the optimal harvesting efforts $E_1 = E_{1\tilde{\delta}}, E_2 = E_{2\tilde{\delta}}$.

8.   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.

Example 1 Consider the parameter values as follows: $\tilde{r} = (2.00, 2.20, 2.40)$, $\tilde{c} = (1.20, 1.35, 1.50)$, $\tilde{d} = (0.30, 0.40, 0.50)$, $\tilde{s} = (0.06, 0.07, 0.08)$, $\tilde{e} = (0.60, 0.70, 0.80)$, $K = 5$ and $m = 0.1$.

We discuss equilibria $S_0$, $S_1$ and $S^*$ of the predator-prey model under different values of $\alpha, w_1$ and $w_2$, which are showed from Table 1 to Table 3.

Tables 1–3 demonstrate the trivial equilibrium is always fixed at $(0, 0)$ with different combinations of $w_1$, $w_2$ and $\alpha$, the values of prey $x$ and predator $y$ stay at zero; for the axial equilibrium, the value of prey $x$ decreases under the same $\alpha$ with increasing $w_1$, predator $y$ is invariant in zero; for the interior equilibrium, the value of prey $x$ increases under the same $\alpha$ with increasing $w_1$, and predator $y$ decreases under the same $\alpha$ with increasing $w_1$.

Figures 1–3 are time series of prey $x$ and predator $y$ with initial values $(2, 6)$ for different $w_1$, $w_2$ and $\alpha$, which show the same results as what 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 $\alpha$ 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 $\alpha$ 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 $\alpha$ which is in accordance with the results in Table 3. As $\alpha$ 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 $\alpha$ increases, prey $x$ decreases and predator $y$ increases. Therefore, the parameter $\alpha$ 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 $\alpha = 0, 0.3, 0.6$ and $\alpha = 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 $\alpha$ with increasing $w_1$, and predator $y$ decreases under fixed $\alpha$ with increasing $w_1$.

Example 2 Consider the parameter values as follows: $\tilde{r} = (1.50, 1.55, 1.60)$, $\tilde{c} = (0.250, 0.275, 0.300)$, $\tilde{d} = (0.450, 0.475, 0.500)$, $\tilde{s} = (0.150, 0.175, 0.200)$, $\tilde{e} = (1.300, 1.325, 1.350)$, $K = 10, m = 0.1, q_1 = 0.92, q_2 = 0.95, p_1 = 20, p_2 = 25, c_1 = 30$ and $c_2 = 15$.

The bionomic equilibria for different combinations of $w_1$, $w_2$ and $\alpha$ are emerged in Table 4. It is observed that $x_\infty$ and $y_\infty$ are invariable with respect to $w_1$, $w_2$ and $\alpha$, and $E_{1\infty}$ and $E_{2\infty}$ are decreasing with fixed $\alpha$ and increasing $w_1$. As well as, together with the value of $\alpha$ increasing, $E_{1\infty}$ and $E_{2\infty}$ 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 $\alpha$ increases, $E_{1\infty}$ and $E_{2\infty}$ are increasing.

Example 3 Consider the parameter values as follows: $\tilde{r} = (2.80, 2.83, 2.85), \tilde{c} = (2.45, 2.55, 2.65), \tilde{d} = (0.012, 0.013, 0.015), \tilde{s} = (0.008, 0.009, 0.010), \tilde{e} = (0.20, 0.21, 0.22), K = 10, m = 0.1, q_1 = 0.45, q_2 = 0.55, p_1 = 30, p_2 = 25, c_1 = 25, c_2 = 12$ and $\tilde{\delta} = (0.8, 0.9, 1.0)$.

The optimal equilibrium and harvesting effort for different combinations $w_1, w_2, \alpha$ and $w_1', w_2', \alpha'$ are appeared in Tables 5–9, respectively. For different combinations $w_1, w_2, \alpha$ and $w_1', w_2', \alpha'$, we see the existence of optimal equilibrium.

In Table 5, when $\alpha' = 0$ it is observed that prey $x_{\tilde{\delta}}$ and predator $y_{\tilde{\delta}}$ are invariant, optimal harvesting effort of the prey $E_{1\tilde{\delta}}$ decreases and then increases, but the optimal harvesting effort of the predator $E_{2\tilde{\delta}}$ increases and then decreases with increasing $w_1'$; when $\alpha' = 0.3$, $x_{\tilde{\delta}}$ decreases and then increases, $y_{\tilde{\delta}}$ is invariant, $E_{1\tilde{\delta}}$ decreases and then increases, and $E_{2\tilde{\delta}}$ increases and then decreases with increasing $w_1'$; when $\alpha' = 0.6$, $x_{\tilde{\delta}}$ increases and then decreases, $y_{\tilde{\delta}}$ is invariant, $E_{1\tilde{\delta}}$ fluctuates between 0.0907 and 0.0908, and $E_{2\tilde{\delta}}$ is invariant with increasing $w_1'$; when $\alpha' = 0.9$, $x_{\tilde{\delta}}$ and $y_{\tilde{\delta}}$ are invariant, $E_{1\tilde{\delta}}$ decreases and then increases, and $E_{2\tilde{\delta}}$ is invariant with increasing $w_1'$. Besides, when fixed $w_1'$ and $w_2'$, $x_{\tilde{\delta}}$ increases, $y_{\tilde{\delta}}$ decreases, $E_{1\tilde{\delta}}$ increases, and $E_{2\tilde{\delta}}$ decreases with the increase of $\alpha'$.

In Table 6, when fixed $\alpha'$, $x_{\tilde{\delta}}$, $y_{\tilde{\delta}}$, $E_{1\tilde{\delta}}$ and $E_{2\tilde{\delta}}$ are invariant with the increase of $w_1'$; when fixed $w_1'$ and $w_2'$, $x_{\tilde{\delta}}$ increases, $y_{\tilde{\delta}}$ decreases, $E_{1\tilde{\delta}}$ increases, and $E_{2\tilde{\delta}}$ increases with the increase of $\alpha'$.

In Table 7, when fixed $\alpha'$, $x_{\tilde{\delta}}$, $y_{\tilde{\delta}}$, $E_{1\tilde{\delta}}$ and $E_{2\tilde{\delta}}$ are invariant with the increase of $w_1'$; when fixed $w_1'$ and $w_2'$, $x_{\tilde{\delta}}$ increases, $y_{\tilde{\delta}}$ is invariant, $E_{1\tilde{\delta}}$ decreases, and $E_{2\tilde{\delta}}$ increases with the increase of $\alpha'$.

In Table 8, when fixed $\alpha'$, $x_{\tilde{\delta}}$, $y_{\tilde{\delta}}$, $E_{1\tilde{\delta}}$ and $E_{2\tilde{\delta}}$ are invariant with the increase of $w_1'$; when fixed $w_1'$ and $w_2'$, $x_{\tilde{\delta}}$ increases, $y_{\tilde{\delta}}$ increases, $E_{1\tilde{\delta}}$ decreases, and $E_{2\tilde{\delta}}$ increases with the increase of $\alpha'$.

In Table 9, when $\alpha' = 0$, $x_{\tilde{\delta}}$ increases and then decreases and $y_{\tilde{\delta}}$ is invariant, $E_{1\tilde{\delta}}$ increases and then decreases, but $E_{2\tilde{\delta}}$ decreases and then increases with increasing $w_1'$; when $\alpha' = 0.3$ or $\alpha' = 0.6$, we can see that $x_{\tilde{\delta}}$ increases and then decreases, $y_{\tilde{\delta}}$, $E_{1\tilde{\delta}}$ and $E_{2\tilde{\delta}}$ are invariant with increasing $w_1'$; when $\alpha' = 0.9$, $x_{\tilde{\delta}}$, $y_{\tilde{\delta}}$ and $E_{1\tilde{\delta}}$ are invariant, $E_{2\tilde{\delta}}$ decreases and then increases with increasing $w_1'$. When fixed $w_1'$ and $w_2'$, $x_{\tilde{\delta}}$ decreases, $y_{\tilde{\delta}}$ is invariant, $E_{1\tilde{\delta}}$ is almost invariant, and $E_{2\tilde{\delta}}$ increases with the increase of $\alpha'$.

From Tables 5–7, it is easy to see that when fixed $w_1, w_2$ and $w_1', w_2', \alpha'$, 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 $\alpha$.

Based on Tables 7–9, we observed that when fixed $\alpha$ and $w_1', w_2', \alpha'$, 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$.

9.   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 $\alpha$. Dynamical behaviors of prey and predator about $\alpha$ 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.

The highlight of this article is to solve the optimal control problem with fuzzy inflation and discount, that is maximize the fuzzy objective functional $\tilde{J}$. The optimal steady state solution is displayed in Tables 5–9. We obtain optimal equilibrium and optimal harvesting effort as respect to $w_1', w_2'$ and $\alpha'$ under five cases including $w_1 = 0$, $w_2 = 1$, $\alpha = 0$ and $w_1 = 0$, $w_2 = 1$, $\alpha = 0.5$ and $w_1 = 0$, $w_2 = 1$, $\alpha = 0.9$ and $w_1 = 0.5$, $w_2 = 0.5$, $\alpha = 0.9$ and $w_1 = 1$, $w_2 = 0$, $\alpha = 0.9$.

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 $\alpha$-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.

Acknowledgements

The authors thank the editor and referees for their careful reading and valuable comments. The work is supported by the National Natural Science Foundation of China (No.12101211) and the Program for Innovative Research Team of the Higher Education Institution of Hubei Province (No.T201812) and the Scientific Research Project of Education Department of Hubei Province (No.B2020090).

Conflict of interest

The authors have no conflict of interest.