Optimal harvesting of an abstract population model with interval biological parameters

We study the optimal harvesting policy for fishery in the marine protected and unreserved areas. In the literature, it is generally assumed that the fish population follows a concrete growth law. In contrast, we consider an abstract model with migration from the reserved area to the unreserved one. Then we examine and analyze the existence and stability of a nontrivial equilibrium point of the model. We also discuss the bionomic equilibrium. After that, we use the Pontryagin maximum principle to obtain the optimal harvest policy, where, instead of the well-known Hamiltonian function, we use the current Hamiltonian function to ease the calculation. Finally, we give some numerical examples to further illustrate our statements, where we also find that in practice the impreciseness of the parameters can influence the existence of the system positive equilibrium.


Introduction
A marine protected area (MPA) is essentially a space in the ocean where human activities are more strictly regulated than the surrounding waters, similarly to parks we have on land. These places are given special protections for natural or historic marine resources by local, state, territorial, native, regional, or national authorities.
The creation of the world biggest marine reserve in Antarctica, covering an area more than six times the UK, has been hailed as a "milestone for conservation. " Agreement was reached by the Commission for the Conservation of Antarctic Marine Living Resources (CCAMLR), made up of 24 countries including the UK and European Union, to protect 1.55 million square kilometers of the Ross Sea. However, zero harvesting of all species within a reserve is rare. Some 72 percent of the marine protected area will be a "no-take" zone, where all fishing is forbidden, whereas other areas will allow some harvesting of fish and krill for scientific research. Lindeboom [1] pointed out that marine reserves can provide an undisturbed system for scientific study. In fact, reserves or "no-take" areas [2] often form a part of larger MPAs that have less protection and may include areas that allow for some consumptive use [3]. In our model, we assume that there exists some harvest in the MPA.
Generally speaking, different species may follow different growth curves. Statistical analysis suggested that the Gompertz curve is the best equation to explain growth, in length, of mackerels during their first growth season [4,5]. A logistic curve is the best growth law to describe growth of micropterus salmoides [6]. Since establishing an MPA could change the species growth pattern [7,8], it is necessary to construct a model using an abstract growth function and analyze its dynamical properties. In Sect. 4, we give some examples with different kinds of growth functions, such as the Gompertz [4], logistic [9], and Smith growth laws [10]. The first two were studied by many researchers, whereas the last one was barely mentioned in the previous literature.
As we know, biological parameters are often obtained imprecisely. The impreciseness can result from various factors. They include internal factors of an ecosystem such as the change of climate, the destruction of habitats, and external factors such as data collection and measurement errors in experiments. There are different techniques to handle this problem, including the stochastic method [11], fuzzy method [12][13][14], and fuzzystochastic method [15]. In the stochastic method the uncertain parameters are replaced by stochastic variables with known probability distributions. In the fuzzy method the uncertain parameters are replaced by fuzzy sets with known membership functions. In the fuzzy-stochastic method a subset of uncertain parameters is replaced as in the fuzzy method, whereas the rest are taken as stochastic variables. Nevertheless, constructing a suitable membership function (or a suitable probability distribution) for the imprecise biological parameters is a difficult problem. Very recently, there are some excellent works about biological models with interval-value biological parameters [16][17][18][19][20][21]. Pal [19] and Wang [16] studied a prey-predator system with interval biological parameters, Sharma and Samanta [18] constructed a two-competing-species harvesting model with imprecise biological parameters, and later, You and Zhao [17] investigated the Gompertz population model with interval-value biological parameters. Indeed, the interval-valued functions give a way to contain all values in a given interval, which can be used to characterize the interval-value biological parameters. In this work, we replace the intrinsic growth rate of the population by an interval-valued function.
To determine the optimal policy of harvesting, we use the Pontryagin maximum principle [22]. In this process the Hamiltonian function is usually used by researchers. However, to make the calculation simpler, in this paper, we use, instead, the current Hamiltonian function [23].
It is usually assumed that there exists migration between both areas [24,25], from the marine reserve toward the zone of open access and vice versa, but in this study, we only consider the migration from the reserved area to the unreserved area, as in the previous literature [26].
This paper is organized as follows. In Sect. 2, we present some introductory material of our model. In Sect. 3, we discuss the existence of a nontrivial equilibrium point and its stability. Furthermore, at the end of this section, we prove the existence of bionomic equilibrium of the model and study the optimal harvesting policy. In Sect. 4, we give some numerical examples verifying our theoretical results. Finally, in Sect. 5. we give some brief conclusions.

Materials and methods
With objects mentioned before, we establish MPA and divide the region into two patches. We consider the two areas with the same kind species but different rates of growth and carrying capacities.
Let x(t) and y(t) represent the population density of the fish resource inside the reserved area and the unreserved area, respectively. Based on the model studied in [26], we make some assumptions: • Although different species may follow different growth curves, we will not consider them as concrete functions, but rather as abstract functions rxG(x) and syF(y) with some properties: (2) There exist positive numbers K , L such that G(K) = 0, F(L) = 0; K and L are the carrying capacity of the reserved area and the unreserved area, respectively. • We assume that the migration rate per unit time is a strictly increasing function m (x) such that m (x) > 0 and m(0) = 0 [27]. Here the model only accounts for migration from the protected area to the open zone. In the numerical section, we will take m(x) = mx 2 x 2 +a 2 [26], where m is the migration rate in (0, 1), a is the half saturation rate of migration, and when x = a, a proportion m 2 of the reserved population emigrates to the unreserved area. Due to its complexity, m(x) = mx 2 x 2 +a 2 is rarely used in other models. • As mentioned in [3], none harvesting of all species within the reserved area is rare, and it is reasonable to suppose that the capture in the reserved area also exists. Then our model can be described by the following 2-D nonlinear differential equations with interval coefficients: where r 0 ∈ [r l , r u ] and s 0 ∈ [s l , s u ]. We take catch functions h j (j = 1, 2) as the products of constant technological parameter q j , the efforts levels E j , and the stock levels in each area; thus h 1 = q 1 E 1 x(t) and h 2 = q 2 E 2 y(t) are harvestings in each area. As proved in [19], the parametric form of system (1) can be written as the following system: (2)

Equilibria
In this part, we prove the existence and uniqueness for the nontrivial steady state.

Theorem 1 System
(2) has a unique positive equilibrium point P( Then So there exists a unique y * in (0, +∞) such that f (y * ) = 0, and this completes the proof.
Theorem 1 shows that if q 1 E 1 is not very large, then model (2) has a positive steady state, which is significant for maintaining a sustainable ecosystem.

Analysis of steady state
Theorem 2 The steady state P(x * , y * ) is a node, and it is globally asymptotically stable.
, and the characteristic equation is Since by assumption G (x) < 0 and F (x) < 0, we conclude that According to [28], P(x * , y * ) is a node, and it is locally asymptotically stable. Let We choose the Dulac function Z(x, y) = 1 xy . Then m(x) y 2 < 0 and thus never changes sign and does not reach zero in the first octant. According to Bendixson-Dulac theorem [29], model (2) has no closed orbits contained entirely in the first octant. Since P(x * , y * ) is locally asymptotically stable, we can state that as t → ∞, every solution of model (2) tends to P(x * , y * ). Thus we complete the proof.

Bionomic equilibrium
In this part, we will show the existence of bionomic equilibrium of system (2).
Let c 1 , c 2 be the harvesting cost per unit efforts, and let p 1 , p 2 be the constant prices per unit biomass of the reserved and unreserved areas. The net revenue is given by Here we consider Theorem 3 If c 1 +θc 2 p 1 q 1 > K and c 1 +θc 2 θp 2 q 2 > L, then model (3) has a unique bionomic equilibrium.
Proof The economic profit is Bionomic equilibrium means both biologic and economic equilibria. Therefore we have π(x, y, E) = 0, dx(t) dt = 0, dy(t) dt = 0, that is, Using equations (5) and (6), we get which is equivalent to the system From (8) we get Since y > 0, we have x < c 1 +θc 2 p 1 q 1 .

Optimal harvesting policy
The aim of this part is to obtain an optimal harvesting policy, which maximizes the present value function. The problem can be described as subject to (3), where π(x, y, E) is given in (4), δ is the rate of annual discount, and the control variable satisfies Since (10) is dependent on time t merely through the discount term, this is an infinitehorizon autonomous problem. In [17,19] authors generally construct a Hamiltonian function. Here we consider the current value Hamiltonian function s given by where λ 1 (t) and λ 2 (t) are the current value multipliers associated with x(t) and y(t), respectively. Then we have Note that the current value Hamiltonian function is linear with respect to the control E, and we will consider the optimal harvesting as a combination of singular and bang-bang controls: When σ (t, l 1 , l 2 ) is not equal to zero, the control changes from one to another depending on the sign of σ (t, l 1 , l 2 ). The solution is called a bang-bang control, and σ (t, l 1 , l 2 ) is called the switching function. However, when σ (t, l 1 , l 2 ) equals zero, the optimal control cannot be obtained through the above procedure, and E δ is called a singular solution. In addition, according to the Pontryagin maximum principle [22], the adjoint system can be stated as follows: After a simple calculation, (11) becomes From (12b) and the interior equilibrium P(x * , y * ) we have For any C 1 , λ 2 (t) is a shadow price, which is bounded, so we can choose C 1 = 0. Then Using (12a)-(12b) and (13), we obtain Similarly, choosing C 2 = 0, we have From the interior equilibrium we get A > 0 and B > 0. Then These equations can be solved numerically to obtain x δ , y δ , E δ . When δ → ∞, equation (14) transforms to which indicates that π(x ∞ , y ∞ , E) = 0. In other words, when the discount rate tends to infinity, and when the net economic revenue tends to zero, then the fishery may close. This conclusion can be explained from the perspective of economics: a high interest rate may lead to a high inflation rate, and if the inflation rate increases rapidly, then the real value of the fishery resource may decrease. Hence the economic market plays an important role in the exploiting of fishery resources, and the owner of the fishery resources stock prefers to exploit it at a no-profit-no-lose basis [30]. From (14) we have and thus π is a decreasing function of δ (≥ 0). Therefore we can assert that when δ = 0, the net profit reaches its maximum.

Numerical examples
In this section, we give some examples to illustrate our statements and exhibit theoretical results. We choose the migration rate function m(x) = mx 2 x 2 +a 2 . We point out that the simulation here is considered as a qualitative, rather than a quantitative analysis. To ensure the rationality of parameters, we determine the value of parameters in our work based on [17,26]. We use Maple to help us finish all the calculations. In Appendix 2, we give a Maple code for Example 2. When 0 ≤ l 1 ≤ 0.3225138920, model (2) has no negative equilibrium point, and then the system tends to a "bad" steady state (see Figs. 1, 2), and the fish in the reserved area would die out. From Fig. 2 we can see that when l 1 = 0.3225138920, x = 0 and y = 7.100000000; and when 0.3225138920 < l 1 ≤ 1, the population of these two tend to a "good" steady state, and the system has a unique positive equilibrium point (see Figs. 3, 4). Example 2 In this example, we assume the two subpopulations follow the Smith growth law [10]: G(x) = K-x K+Dx and F(y) = L-y L+Dy , where D is a constant. We choose D = 1, K = 5, L = 8, m = 0.2, a = 0.8, r l = 3.6, r u = 5.2, s l = 5, s u = 6.3, q 1 E 1 = 3, q 2 E 2 = 4. Then, for any l 1 , l 2 ∈ [0, 1], the trivial equilibrium point (0, 0) obviously exists. The nontrivial equilibrium point P(x * , y * ), the eigenvalues of variational matrices at the corresponding equilibrium point, and the equilibrium point nature and its stability are given in Table 1, where for brevity, we take l 1 = l 2 .  Example 3 In this example, we present the optimal harvesting policy for two species with different growth laws, the logistic and Gompertz laws. To illustrate the difference coming from the two-species growth curve, we assume that they have the same parameters K = 5, L = 8, r l = 5, r u = 5.35, s l = 4, s u = 4.25, m = 0.2, a = 0.8, θ = 0.3, δ = 0.05, p 1 = 4, p 2 = 5,    q 1 = 0.2, q 2 = 0.1, c 1 = 2, c 2 = 3. Then the optimal equilibrium point (x δ , y δ ), the optimal harvesting effort E δ , and the net profit π of these two species are given in Tables 2 and 3. We can catch the fact that as l 1 , l 2 increase, the population in the unreserved area y δ decreases, in the reserved area x δ increases, and the optimal harvesting effort E δ and the net profit π increase as well.
Example 4 In this example, we illustrate that when the system reaches its equilibrium point, the population of the unreserved area is not always higher than the reserved area.

Discussion
In this work, we studied the optimal harvesting of fishery in the marine reserved area and its adjacent area. We supposed that the population is divided into two subpopulations and the growth curve of each population can be the same or different. Due to the fact that biological parameters are often obtained imprecisely, we considered the intrinsic growth rate of these two subpopulations in a certain range and replaced them with interval-valued functions, which contain all values in a given interval. To make the model less complicated, we assumed that the migration only exists from the marine protected area toward the unreserved area.
We proved that when r 1-l 1 l r l 1 u lim x→0 G(x)q 1 E 1 > 0, model (2) has a unique positive equilibrium point, which is a node and is globally asymptotically stable. In our model, the value of l 1 and l 2 has no influence on the nature of equilibrium, but in other models, such as the prey-predator model, the application of interval-valued functions may change the nature of equilibrium, and it is an interesting problem, which can be studied by the related researchers. Moreover, we analyzed the existence of a bionomic equilibrium of model (3); under the assumption c 1 +θc 2 p 1 q 1 > K , c 1 +θc 2 θp 2 q 2 > L, there exists a unique bionomic equilibrium. To determine the optimal harvesting, we also supposed that the different levels of effort are proportional, and in this case, it is significant to use the current Hamilton function.
According to the Pontryagin maximum principle, the optimal policy is a combination of bang-bang and singular controls. We obtained the level of effort E δ , under which the stock of both subpopulations can maintain at a equilibrium level and the current value can reach its maximum. Then we pointed out that if the discount rate increases to infinity, then the net economic revenue tends to zero, tending toward a bionomic equilibrium, and if the discount rate is zero, then it gains up to its maximum.
In our work, we considered an abstract growth function, which may be used in another biological model to derive some more general conclusions. We will study it in the future research.