Dynamics and Optimal Taxation Control in a Bioeconomic Model with Stage Structure and Gestation Delay

Abstract

A prey-predator model with gestation delay, stage structure for predator, and selective harvesting effort on mature predator is proposed, where taxation is considered as a control instrument to protect the population resource in prey-predator biosystem from overexploitation. It shows that interior equilibrium is locally asymptotically stable when the gestation delay is zero, and there is no periodic orbit within the interior of the first quadrant of state space around the interior equilibrium. An optimal harvesting policy can be obtained by virtue of Pontryagin's Maximum Principle without considering gestation delay; on the other hand, the interior equilibrium of model system loses as gestation delay increases through critical certain threshold, a phenomenon of Hopf bifurcation occurs, and a stable limit cycle corresponding to the periodic solution of model system is also observed. Finally, numerical simulations are carried out to show consistency with theoretical analysis.

1. Introduction

Recently, the dynamics of a class of stage-structured prey-predator models with gestation delay have been studied by several authors [1–9]. Especially, there is a well-developed theory of stage structured models which incorporates time delay into maturity of population [1]. The prey-predator models with stage structure and gestation delay play an important role in the modelling of multispecies population dynamics. Generally, individuals in each stage are identical in biological characteristics, and the reproduction of mature predator population after predating the prey is not instantaneous but mediates by some discrete time lag required for gestation of mature predator. Xu and Ma [9] proposed the following model, 𝑥̇(𝑡)=𝑥(𝑡)𝑟−𝛼𝑥(𝑡)−𝑎₁𝑦₂(𝑡),𝑦₁̇(𝑡)=𝑎₂𝑥(𝑡−𝜏)𝑦₂(𝑡−𝜏)−𝑟₁𝑦₁(𝑡)−𝑑𝑦₁(𝑡),𝑦₂̇(𝑡)=𝑑𝑦₁(𝑡)−𝑟₂𝑦₂(𝑡)−𝛽𝑦²₂(𝑡), (1.1) where 𝑥(𝑡) represents the density of prey population at time 𝑡, 𝑦₁(𝑡) and 𝑦₂(𝑡) represent the density of immature and mature predator population at time 𝑡, respectively; 𝑟 is the intrinsic growth rate of prey population, 𝑟₁ is death rate of immature predator and 𝑟₂ is the death rate of mature predator, 𝛼 and 𝛽 represents the intraspecific competition rate of prey and mature predator population, respectively; 𝑎₂/𝑎₁ is the rate of converting prey population into new immature predator population. The constant 𝜏≥0 denotes the gestation delay of mature predator, and 𝜏≥0 is based on the assumption that the reproduction of predator population after predating the prey population is not instantaneous but mediates by some discrete time lag required for gestation of mature predator population. 𝑑>0 denotes the proportional transforming rate from immature predator population to mature predator population. All the parameters mentioned above are positive constants.

It is well known that exploitation of several biological resources has been increased by the growing human needs for more food and energy, which attracts a global concern to protect the limited biological resources. Consequently, regulation of exploitation of biological resources has become a problem of major concern in view of dwindling resource stocks and deteriorating environment. It should be noted that some techniques and issues associated with bioeconomic exploitation have been discussed in details by Clark [10]. Due to its economic flexibility, taxation is usually considered as possible governing instruments in regulation for harvesting to keep the damage to the ecosystem minimal.

Recently, there has been considerable interest in the modeling of harvesting of biological resources. In these models, the harvesting effort is considered to be a dynamic variable, several kinds of harvesting policies are utilized to study the dynamical behavior of the model system. Furthermore, optimal harvesting policies with taxation are also discussed. Ganguly and Chaudhuri [11], Krishna et al. [12], and Dubey et al. [13, 14] investigated the optimal harvesting of a class of models of a single fishery species with taxation as a control. Chaudhuri et al. [15–17], Pradhan and Chaudhuri [18], and Kar et al. [19–25] studied the optimal taxation policies for harvesting of the prey-predator system. However, from the above literature survey, it may be pointed out that no attempt has been made to study the optimal taxation policy of a stage-structured prey-predator system. Furthermore, taxation instrument is discussed to control overharvesting from prey-predator system with gestation delay in [26]. However, stage structure of predator population is not considered, and the periodic orbit within the interior of the first quadrant of state space around interior equilibrium is also not investigated in [26].

The stability analysis of interior equilibrium is performed in the third section. It reveals that when gestation delay is zero, the interior equilibrium is locally asymptotically stable. It is also found that equilibrium switch occurs due to variation of gestation delay. Furthermore, an optimal harvesting policy for mature predator is also discussed in the absence of gestation delay. We aim to find an optimal harvesting policy which guarantees an ever-lasting exploitation of the biological resource and maximizes the benefits resulting from the harvesting. Numerical simulations are provided to support the analytical findings in this paper. Finally, this paper ends with a conclusion.

2. Model Formulation

Based on the above analysis, the work done by Xu and Ma in [9] is extended by incorporating harvest effort on mature predator, and taxation is chosen to control the conservation of biological resource. In this paper, a prey-predator model with gestation delay and stage structure for predator is established. It is assumed that mature predator is subject to a dynamic harvesting. To conserve the population in the prey-predator ecosystem, the regulatory agency imposes a taxation 𝜎>0 per unit biomass of mature predator (𝜎<0 denotes the subsidies given to the harvesting effort). Based on the above aspects, the model can be governed by the following differential equations: 𝑥̇(𝑡)=𝑥(𝑡)𝑟−𝛼𝑥(𝑡)−𝑎₁𝑦₂(𝑡),𝑦₁̇(𝑡)=𝑎₂𝑥(𝑡−𝜏)𝑦₂(𝑡−𝜏)−𝑟₁𝑦₁(𝑡)−𝑑𝑦₁(𝑡),𝑦₂̇(𝑡)=𝑑𝑦₁(𝑡)−𝑟₂𝑦₂(𝑡)−𝛽𝑦²₂(𝑡)−𝑞𝐸(𝑡)𝑦₂(𝑡),𝐸̇(𝑡)=𝛼₀𝐸(𝑡)𝑝−𝜎𝑞𝑦₂(𝑡)−𝑐, (2.1) where initial conditions are as follows: 𝑥(𝑡)=𝜓₁(𝑡)>0,𝑦₂(𝑡)=𝜓₂(𝑡)[),𝑦>0,𝑡∈−𝜏,0₁(0)=𝜓₃(0)>0,𝐸(0)=𝜓₄(0)>0. (2.2)

The harvesting term 𝐸(𝑡) is assumed to be proportional to both stock level and effort, which follows the catch per unit effort hypothesis [10]. The constant 𝑞 is the catchability coefficient, 𝑝 is the fixed price per unit of predator species, 𝑐 is the fixed cost of harvesting per unit of effort, and 𝛼₀ is called stiffness parameter measuring the strength of reaction of harvesting effort. The parameters mentioned above are all positive constants.

3. Qualitative Analysis of Model System

From the view of ecological management, we only concentrate on the interior equilibrium of the model system in this paper, since the biological meaning of the interior equilibrium implies that juvenile preys, mature preys, predators, and harvesting effort on predators all exist, which are relevant to our study.

It can be obtained that the only interior equilibrium of the model system (2.1) is 𝑃^∗(𝑥^∗,𝑦^∗₁,𝑦^∗₂,𝐸^∗), where 𝑥^∗=(𝑟−𝑎₁𝑦^∗₂)/𝛼, 𝑦^∗₁=𝑎₂𝑦^∗₂(𝑟−𝑎₁𝑦^∗₂)/𝛼(𝑟₁+𝑑), 𝐸^∗=(𝑎₂𝑑(𝑟−𝑎₁𝑦^∗₂)−𝛼(𝑟₂+𝛽𝑦^∗₂)(𝑟₁+𝑑))/𝛼𝑞(𝑟₁+𝑑), and 𝑦^∗₂=𝑐/(𝑝−𝜎)𝑞. It is easy to show that interior equilibrium exists, provided the following conditions are satisfied: 0<𝜎<min𝑝−𝑐𝑎₁𝑐𝑎𝑟𝑞,𝑝−¹𝑎₂𝑑+𝛼𝛽𝑟₁+𝑑𝑞𝑎₂𝑑𝑟−𝛼𝑟₂𝑟₁+𝑑, (3.1) which provides the range of taxation for the existence of interior equilibrium. This range of taxation may be utilized when the regulatory agency establishes relevant agencies for harvesting.

The model system (2.1) can be interpreted as the matrix form: 𝑋̇(𝑡)=𝐻(𝑋(𝑡)), (3.2) where 𝑋(𝑡)=(𝑥(𝑡),𝑦₁(𝑡),𝑦₂(𝑡),𝐸(𝑡))^𝑇∈ℝ⁴, and 𝐻(𝑋(𝑡)) is given as follows, 𝐻(𝑋(𝑡=𝐻))¹(𝑋(𝑡𝐻))₂(𝑋(𝑡𝐻))₃(𝑋(𝑡𝐻))₄(𝑋(𝑡=𝑥(𝑡)(𝑡)))𝑟−𝛼𝑥−𝑎₁𝑦₂(𝑡)𝑎₂𝑥(𝑡−𝜏)𝑦₂(𝑡−𝜏)−𝑟₁𝑦₁(𝑡)−𝑑𝑦₁(𝑡)𝑑𝑦₁(𝑡)−𝑟₂𝑦₂(𝑡)−𝛽𝑦²₂(𝑡)−𝑞𝐸(𝑡)𝑦₂(𝑡)𝛼₀𝐸(𝑡)𝑝−𝜎𝑞y₂(𝑡)−𝑐. (3.3) Let ℝ⁴₊=[0,∞)⁴ be the nonnegative octant in ℝ⁴, then 𝐺∶ℝ₊⁴⁺¹→ℝ⁴ is locally Lipschitz and satisfies the condition 𝐻_𝑖(𝑋(𝑡))|_𝑋∈ℝ4≥0.

Due to lemma in [27] and Theorem A.4 in [28], any solution of the model system (2.1) with positive initial conditions exist uniquely, and each component of the solution remains within the interval [0,𝑏) for some 𝑏>0. Furthermore, if 𝑏<∞, then limsup[𝑥(𝑡)+𝑦₁(𝑡)+𝑦₂(𝑡)+𝐸(𝑡)]=∞. Hence, this completes the positivity for the solutions of model system (2.1).

Now we consider boundedness of positive solutions (𝑥(𝑡),𝑦₁(𝑡),𝑦₂(𝑡),𝐸(𝑡)), and firstly choose the function 𝑊₁(𝑡)=𝑥(𝑡−𝜏)+𝑦₁(𝑡).

For 𝑡>𝑇₁+𝜏,𝑎₁−𝑎₂>0 (𝑇₁ is some fixed positive time), by calculating the time derivative of 𝑊₁(𝑡) along the solutions of model system (2.1), we get 𝑊̇₁(𝑡)=𝑟𝑥(𝑡−𝜏)−𝛼𝑥²(𝑡−𝜏)−𝑎₁−𝑎₂𝑥(𝑡−𝜏)𝑦₂(𝑡−𝜏)−𝑟₁+𝑑𝑦₁(𝑡). (3.4)

By virtue of positiveness of solution 𝑥(𝑡−𝜏), it is easy to show that 𝑟𝑥(𝑡−𝜏)−𝛼𝑥²()≤𝑟𝑡−𝜏².4𝛼 (3.5)

Based on the positiveness of solution 𝑥(𝑡−𝜏), 𝑦₂(𝑡−𝜏) and assumption 𝑎₁−𝑎₂>0, it is easy to show that 𝑎₁−𝑎₂𝑥(𝑡−𝜏)𝑦₂(𝑡−𝜏)>0. (3.6)

Hence, it is easy to show that 𝑊̇₁(𝑡)<𝑟²−𝑟4𝛼¹+𝑑𝑦₁(𝑡)<𝑟²−𝑟4𝛼¹+𝑑𝑥(𝑡−𝜏)+𝑦₁(𝑡), (3.7) which follows that there exists a positive quantity 𝑀₁ such that 0<𝑊₁(𝑡)<𝑀₁ for all large 𝑡>𝑇₁+𝜏. It proves the boundedness of positive solution 𝑥(𝑡), 𝑦₁(𝑡).

Let 𝑊₂(𝑡)=𝑦₂(𝑡), by calculating the time derivative of 𝑊₂(𝑡) along the solutions of model system (2.1), we have 𝑊̇₂(𝑡)≤𝑑𝑦₁(𝑡)−𝑟₂𝑦₂(𝑡). (3.8)

By virtue of the positivity of the solutions of model system (2.1) and the boundedness of 𝑦₁(𝑡) mentioned above, it follows that there exists a positive quantity 𝑀₂ such that 𝑊̇₂(𝑡)<𝑀₂−𝑟₂𝑦₂(𝑡) for all large time 𝑡>𝑇₂ (𝑇₂ is some fixed positive time). From the above differential inequality it follows that, there exists a positive quantity 𝑀₃ such that 0<𝑊₂(𝑡)<𝑀₃ for all large 𝑡>𝑇₂, which proves the boundedness of positive solution 𝑦₂(𝑡).

Let 𝑊₃(𝑡)=𝑦₂(𝑡)+𝐸(𝑡), by calculating the time derivative of 𝑊₃(𝑡) along the solutions of model system (2.1), we have 𝑊̇₃(𝑡)=𝑑𝑦₁(𝑡)−𝑟₂𝑦₂(𝑡)−𝛽𝑦²₂(𝑡)−1−𝛼₀𝑝−𝜎𝑞𝐸(𝑡)𝑦₂(𝑡)−𝑐𝛼₀𝐸(𝑡). (3.9)

By virtue of positiveness and boundedness of solution 𝑦₁(𝑡) and 𝑦₂(𝑡), it follows that there exists a positive quantity 𝑀₄ such that 𝑑𝑦₁(𝑡)−𝑟₂𝑦₂(𝑡)−𝛽𝑦²₂(𝑡)≤𝑀₄. Furthermore, under the following assumption: 1−𝛼₀𝑝−𝜎>0 (3.10) it is easy to show that 𝑊̇₃(𝑡)≤𝑀₄−𝑐𝛼₀𝐸(𝑡) for for all large 𝑡>𝑇₃ (𝑇₃ is some fixed positive time), which derives that there exists a positive quantity 𝑀₅ such that 0<𝑊₃(𝑡)<𝑀₅ for all large 𝑡>𝑇₃, which proves boundedness of positive solution 𝐸(𝑡).

Remark 3.1. Since the components (𝑥(𝑡),𝑦₁(𝑡),𝑦₂(𝑡)) of solution of model system (2.1) represent the population in the prey-predator system, the positivity implies that the population survives, and the boundedness reveals a natural restriction to growth as a consequence of limited resources. Furthermore, with the purpose of maintaining the sustainable development of prey-predator system, the harvesting cannot increase without any restriction. As analyzed above, the assumption (3.10) provides the range of taxation for the boundedness of harvesting effort. It is an inspiration for people to regulate the harvesting effort by means of economic instrument.
The Jacobian of model system (2.1) evaluated at the only interior equilibrium 𝑃^∗ leads to the following characteristic equations: 𝜆+𝛼𝑥^∗0𝑎₁𝑥^∗0−𝑎₂𝑦^∗₂𝑒^−𝜆𝜏𝜆+𝑟₁+𝑑−𝑎₂𝑥^∗𝑒^−𝜆𝜏00−𝑑𝜆+𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂𝑞𝑦^∗₂00−𝛼₀𝐸^∗𝑀(𝜆)(𝜆)𝑒𝑝−𝜎𝑞𝜆=0,+𝑁^−𝜆𝜏=0, (3.11)(3.12) where 𝑀(𝜆)=𝜆⁴+𝑚₁𝜆³+𝑚₂𝜆²+𝑚₃𝜆+𝑚₄,𝑁(𝜆)=𝑛₂𝜆²+𝑛₃𝑚𝜆,₁=𝑟₁+𝑑+𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝛼𝑥^∗,𝑚₂=𝑟₁+𝑑𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝛼𝑥^∗𝑟₁+𝑑𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝑐𝛼₀𝑞𝐸^∗,𝑚₃=𝛼𝑥^∗𝑟₁+𝑑𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝑐𝛼₀𝑞𝐸^∗𝑟₁+𝑑+𝛼𝑥^∗,𝑚₄=𝑐𝛼₀𝑞𝐸^∗𝛼𝑥^∗𝑟₁,𝑛+𝑑₂=−𝑎₂𝑑𝑥^∗,𝑛₃=𝑎₂𝑑𝑥^∗𝑎₁𝑦^∗₂−𝛼𝑥^∗. (3.13)

3.1. Case I: Gestation Delay 𝜏=0

In absence of gestation delay, stability of interior equilibrium 𝑃^∗ is investigated, and an optimal harvesting policy with taxation control is also investigated.

3.1.1. Local Stability Analysis

In the absence of gestation delay (𝜏=0), model system (2.1) is written as follows: 𝑥̇(𝑡)=𝑥(𝑡)𝑟−𝛼𝑥(𝑡)−𝑎₁𝑦₂(𝑡),𝑦₁̇(𝑡)=𝑎₂𝑥(𝑡)𝑦₂(𝑡)−𝑟₁𝑦₁(𝑡)−𝑑𝑦₁(𝑡),𝑦₂̇(𝑡)=𝑑𝑦₁(𝑡)−𝑟₂𝑦₂(𝑡)−𝛽𝑦²₂(𝑡)−𝑞𝐸(𝑡)𝑦₂(𝑡),𝐸̇(𝑡)=𝛼₀𝐸(𝑡)𝑝−𝜎𝑞𝑦₂(𝑡)−𝑐, (3.14) and (3.12) can be written as follows: 𝜆⁴+𝑚₁𝜆³+𝑚₂+𝑛₂𝜆²+𝑚₃+𝑛₃𝜆+𝑚₄=0. (3.15)

It can be shown that 𝑚₁>0,𝑚₄𝑚>0,₁𝑚₂+𝑛₂−𝑚₃+𝑛₃=𝛼𝑥^∗𝑟₁+𝑑²+𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂²+𝑟₁+𝑑𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝛼𝑥^∗2𝑟₁+𝑑+𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝑐𝛼₀𝑞𝐸^∗𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝛼𝛽𝑟₁+𝑑𝑦₂^∗2𝑟₁+𝑑+𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂𝑚>0,₁+𝑛₁𝑚₂+𝑛₂𝑚₃+𝑛₃−𝑚₃+𝑛₃−𝑚²₁𝑚₄=𝛼𝑥^∗𝑟₁+𝑑𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂×𝑟₁+𝑑+𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂𝑟₁+𝑑𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝛼𝑥^∗𝑟₁+𝑑+𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝛼𝑥^∗+𝑎₂𝑑𝑥^∗𝑟₁+𝑑+𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝛼𝑥^∗+𝑎₁𝑦^∗₂+𝑐𝛼₀𝑞𝐸^∗𝑟₁+𝑑+𝛼𝑥^∗×𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂𝑟₁+𝑑𝑟₁+𝑑+𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝑐𝛼₀𝑞𝐸^∗+𝑎₂𝑑𝑥^∗𝑟₁+𝑑+𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝛼𝑥^∗+𝑎₁𝑦^∗₂+𝛼𝑥^∗𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂²+𝑎₂𝑑𝑥^∗𝑎₁𝑦^∗₂−𝛼𝑥^∗𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂𝑟₁+𝑑𝑟₁+𝑑+𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝑐𝛼₀𝑞𝐸^∗+𝛼𝑥^∗𝑟₁+𝑑+𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝛼𝑥^∗𝑟₁+𝑑+𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝑎₂𝑑𝑥^∗2𝑟₁+𝑑+𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂+𝛼𝑥^∗>0. (3.16)

Based on the above analysis, it can be concluded that the roots of (3.15) have negative real parts by using the Routh-Hurwitz criteria [10]. Consequently, the interior equilibrium 𝑃^∗ is locally asymptotically stable in absence of gestation delay.

Furthermore, let 𝐽^∗ represent the variational matrix of the model system (3.14) at 𝑃^∗, then tr𝐽^𝑇₀𝑥^∗,𝑦^∗₁,𝑦^∗₂,𝐸^∗d𝑡=^𝑇₀tr−𝛼𝑥^∗0−𝑎₁𝑥^∗0𝑎₂𝑦^∗₂−𝑟₁+𝑑𝑎₂𝑥^∗00𝑑−𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂−𝑞𝑦^∗₂00𝛼₀𝐸^∗𝑝−𝜎𝑞0d𝑡=−^𝑇₀𝛼𝑥^∗+𝑟₁+𝑑+𝑑𝑦^∗₁𝑦^∗₂+𝛽𝑦^∗₂d𝑡<0. (3.17)

It can be easily verified that −(𝛼𝑥^∗+𝑟₁+𝑑+(𝑑𝑦^∗₁/𝑦^∗₂)+𝛽𝑦^∗₂)<0 based on the positivity of the solutions of model system.

Hence, ∫^𝑇₀tr(𝐽(𝑥^∗,𝑦^∗₁,𝑦^∗₂,𝐸^∗))d𝑡<0, which eliminates the existence of Hopf bifurcating periodic solution in the vicinity of 𝑃^∗.

Subsequently, we will show the nonexistence of periodic orbit encircling 𝑃^∗. Let ℎ(𝑥(𝑡),𝑦₁(𝑡),𝑦₂(𝑡),𝐸(𝑡))=1/𝑥(𝑡)𝑦₁(𝑡)𝑦₂(𝑡)𝐸(𝑡). According to the positivity of solutions of the model system (3.14), it is obvious that ℎ(𝑥(𝑡),𝑦₁(𝑡),𝑦₂(𝑡),𝐸(𝑡))>0.

Define Δ(𝑥(𝑡),𝑦₁(𝑡),𝑦₂(𝑡),𝐸(𝑡))=(𝜕/𝜕𝑥)(𝐻₁ℎ)+(𝜕/𝜕𝑦₁)(𝐻₂ℎ)+(𝜕/𝜕𝑦₂)(𝐻₃ℎ)+(𝜕/𝜕𝐸)(𝐻₄ℎ), where 𝐻_𝑖, 𝑖=1,2,3,4 have been defined before, then we have Δ𝑥(𝑡),𝑦₁(𝑡),𝑦₂(𝑡)(𝑡)𝛼,𝐸=−𝑦₁𝑦₂𝐸−𝑎₂𝑦²₁𝐸−𝑑𝑥𝑦²₂𝐸−𝛽𝑥𝑦₁𝐸<0 (3.18) for 𝑥(𝑡), 𝑦₁(t), 𝑦₂(𝑡), 𝐸(𝑡)>0, since all other parameters are strictly positive.

Therefore, there will be no periodic orbit within the interior of the first quadrant of state space around 𝑃^∗ based on Benedixon-Dulac criterion [29].

3.1.2. Optimal Harvesting Policy

With the purpose of planning harvesting and keeping sustainable development of ecosystem, we design an optimal harvesting policy to maximize the total discounted net revenue from the harvesting using taxation as a control instrument. The path traced out by (𝑥(𝑡),𝑦₁(𝑡),𝑦₂(𝑡),𝐸(𝑡)) with optimal taxation 𝜎(𝑡) is also investigated.

Net economic revenue to the society 𝜋(𝑥(𝑡),𝑦₁(𝑡),𝑦₂(𝑡),𝐸(𝑡),𝜎,𝑡) = Net economic revenue of harvesting + Net economic revenue to the regulatory agency = (𝑝−𝜎(𝑡))𝑞𝑦₂(𝑡)𝐸(𝑡)−𝑐𝐸(𝑡)+𝜎𝑞𝑦₂(𝑡)𝐸(𝑡) = (𝑝𝑞𝑦₂(𝑡)−𝑐)𝐸(𝑡).

Our objective is to maximize the following optimization problem: max^∞₀𝑒^−𝛿𝑡𝑝𝑞𝑦₂(𝑡)−𝑐d𝑡, (3.19) where 𝛿 is the instantaneous annual rate of discount, and the optimization problem is subject to the model system (3.14).

By using the Pontryagin's Maximum Principle [10], the associated Hamiltonian function is constructed by 𝐻𝑥(𝑡),𝑦₁(𝑡),𝑦₂(𝑡),𝐸(𝑡),𝜎(𝑡),𝑡=𝑒^−𝛿𝑡𝑝𝑞𝑦₂(𝑡)−𝑐𝐸(𝑡)+𝜆₁(𝑡)𝑥(𝑡)𝑟−𝛼𝑥(𝑡)−𝑎₁𝑦₂(𝑡)+𝜆₂(𝑡)𝑎₂𝑥(𝑡)𝑦₂(𝑡)−𝑟₁𝑦₁(𝑡)−𝑑𝑦₁(𝑡)+𝜆₃(𝑡)𝑑𝑦₁(𝑡)−𝑟₂𝑦₂(𝑡)−𝛽𝑦²₂(𝑡)−𝑞𝐸(𝑡)𝑦₂(𝑡)+𝜆₄(𝑡)𝛼₀𝐸(𝑡)𝑝−𝜎(𝑡)𝑞𝑦₂(𝑡)−𝑐, (3.20) where 𝜆₁, 𝜆₂, 𝜆₃, 𝜆₄ are adjoint variables. 𝜎 is the control variable satisfying the constraints 𝜎_min≤𝜎≤𝜎_max. 𝜎_max and 𝜎_min represent a feasible upper and lower limit of the taxation for harvesting effort, respectively. Specially, 𝜎_min<0 implies that subsidies have the effect of increasing the rate of expansion of the harvesting. According to [29], the condition for a singular control to be optimal can be obtained, that is, 𝜕𝐻/𝜕𝜎=0, from which we get 𝜆₄(𝑡)=0. (3.21)

For adjoint variables 𝜆_𝑖(𝑡), 𝑖=1,2,3,4, we have d𝜆₁d𝑡=−𝜕𝐻,𝜕𝑥d𝜆₂d𝑡=−𝜕𝐻𝜕𝑦₁,d𝜆₃d𝑡=−𝜕𝐻𝜕𝑦₂,d𝜆₄d𝑡=−𝜕𝐻,𝜕𝐸d𝜆₁=d𝑡2𝛼𝑥+𝑎₁𝑦₂−𝑟𝜆₁−𝑎₂𝑦₂𝜆₂,d𝜆₂=𝑟d𝑡¹+𝑑𝜆₂−𝑑𝜆₃,d𝜆₃d𝑡=𝑎₁𝑥𝜆₁−𝑎₂𝑥𝜆₂+2𝛽𝑦₂+𝑟₂+𝑞𝐸𝜆₃−𝑝𝑞𝐸𝑒^−𝛿𝑡,d𝜆₄d𝑡=𝑞𝑦₂𝜆₃−𝑒^−𝛿𝑡𝑝𝑞𝑦₂−𝑐. (3.22)(3.23)(3.24)(3.25)(3.26)

Based on (3.21), it follows from (3.26) that 𝜆₃(𝑡)=𝑒^−𝛿𝑡𝑐𝑝−𝑞𝑦₂. (3.27)

In order to obtain an optimal equilibrium solution, by considering the interior equilibrium 𝑃^∗ and solving (3.25), 𝜆₁(𝑡)=𝑎₂𝜆₂−𝐴₁𝑒^−𝛿𝑡𝑎₁, (3.28) where 𝐴₁=(𝑝(2𝛽𝑦^∗₂+𝑟₂)−(𝑝−(𝑐/𝑞𝑦^∗₂))(2𝛽𝑦^∗₂+𝑟₂+𝑞𝐸^∗+𝛿))/𝑥^∗. By virtue of (3.28), (3.23) can be rewritten as follows: d𝜆₂=d𝑡2𝛼𝑥^∗−𝑟𝜆₂−𝐴₂𝑒^−𝛿𝑡, (3.29) where 𝐴₂=(2𝛼𝑥^∗+𝑎₁𝑦^∗₂−𝑟+𝛿)𝐴₁/𝑎₂.

It is easy to obtain the solution of the above linear differential equation 𝜆₂(𝑡)=𝐴₃𝑒^−𝛿𝑡, (3.30) where 𝐴₃=(2𝛼𝑥^∗+𝑎₁𝑦^∗₂−𝑟+𝛿)𝐴₁/𝑎₂(𝛿𝑥^∗+2𝛼𝑥^∗+𝑟), based on (3.30), by solving (3.24) it follows that 𝜆₃(𝑡)=𝑟₁+𝑑+𝛿𝑑𝐴₃𝑒^−𝛿𝑡. (3.31) Substituting (3.31) into (3.27), we have 𝑐𝑝−𝑞𝑦^∗₂=𝑟1+¹+𝛿𝑑𝐴₃, (3.32) which provides an equation to the singular path and gives the optimal equilibrium levels of population 𝑥^∗=𝑥_𝛿,𝑦^∗₁=𝑦_1𝛿,𝑦^∗₂=𝑦_2𝛿. Then the optimal equilibrium levels of harvesting effort and taxation can be obtained as follows: 𝐸_𝛿=𝑎₂𝑑𝑟−𝑎₁𝑦_2𝛿−𝛼𝑟₂+𝛽𝑦_2𝛿𝑟₁+𝑑𝛼𝑞𝑟₁,𝜎+𝑑_𝛿𝑐=𝑝−𝑞𝑦_2𝛿. (3.33)

Remark 3.2. According to [30], 𝜆_𝑖(𝑡)𝑒^𝛿t(𝑖=1,2,3,4) represent unusual shadow prices along the singular path. From (3.30), (3.32), and (3.36), it may be concluded that these shadow prices remain constant over time interval in an optimum equilibrium when they strictly satisfy the transversality condition at ∞ [31]. Furthermore, they remain bounded as 𝑡→∞.
Considering the interior equilibrium, (3.27) can be written as 𝜆₃(𝑡)𝑞𝑦^∗₂=𝑒^−𝛿𝑡𝑝𝑞𝑦^∗₂−𝑐=𝑒^−𝛿𝑡𝜕𝜋,𝜕𝐸 (3.34) which implies that the user's total cost of harvesting per unit effort is equal to the discounted values of the future price at the steady state effort level.

3.2. Case II: Gestation Delay 𝜏>0

In this section, a stability switch in model system (2.1) due to gestation delay is investigated. Furthermore, a phenomenon of Hopf bifurcation occurs, and a stable limit cycle corresponding to the periodic solution of model system (2.1) is observed.

3.2.1. Local Stability Analysis

Let 𝜆=𝑖𝜔 be a root of (3.12), where 𝜔 is positive. Substitute 𝜆=𝑖𝜔 into (3.12), and separate the real and imaginary parts, then two transcendental equations can be obtained as follows: 𝜔⁴−𝑚₂𝜔²+𝑚₄=−𝑛₃𝜔sin(𝜔𝜏)+𝑛₂𝜔²(),𝑚cos𝜔𝜏₁𝜔³−m₃𝜔=𝑛₃𝜔cos(𝜔𝜏)+𝑛₂𝜔²sin(𝜔𝜏). (3.35)(3.36)

By squaring and adding these two equations, it can be obtained that, 𝜔⁸+𝐵₁𝜔⁶+𝐵₂𝜔⁴+𝐵₃𝜔²+𝐵₄=0, (3.37) where 𝐵₁=𝑚²₁−2𝑚₂, 𝐵₂=𝑚²₂+2𝑚₄−2𝑚₁𝑚₃−𝑛²₂, 𝐵₃=𝑚²₃−2𝑚₂𝑚₄−𝑛²₃, 𝐵₄=𝑚²₄ and 𝑚_𝑖, 𝑛_𝑗 (𝑖=1,2,3,4;𝑗=2,3) have been defined in (3.12).

According to the values of 𝐵_𝑖 (𝑖=1,2,3,4) and the Routh-Hurwitz criteria [10], a simple assumption of the existence of a positive root for (3.37) is 𝐵₃<0. If 𝐵₃<0 holds, then (3.37) has a positive root 𝜔₀, and (3.12) has a pair of purely imaginary roots of the form ±𝑖𝜔₀. Consequently, it can be obtained by eliminating sin(𝜔𝜏) from (3.35) and (3.36) that ()=𝜔cos𝜔𝜏⁴+𝑚₁𝜔³−𝑚₂𝜔²−𝑚₃𝜔+𝑚₄𝑛₂𝜔²+𝑛₃𝜔, (3.38) where the 𝜏_𝑘 corresponding to 𝜔₀ is as follows, 𝜏_𝑘=1𝜔₀𝜔arccos⁴+𝑚₁𝜔³−𝑚₂𝜔²−𝑚₃𝜔+𝑚₄𝑛₂𝜔²+𝑛₃𝜔+2𝑘𝜋𝜔₀,𝑘=0,1,2,…. (3.39) By virtue of Butler's lemma [32], it can be concluded that the interior equilibrium remains stable for 𝜏<𝜏₀, as 𝑘=0.

3.2.2. Hopf Bifurcation

In this section, the condition for Hopf bifurcation in [29] is utilized to investigate whether there is a phenomenon of Hopf bifurcation as 𝜏 increases through 𝜏₀. As stated above, 𝜆=𝑖𝜔₀ represents a purely imaginary root of (3.12), and it follows from the above analysis that |𝑀(𝜆)|=|𝑁(𝜆)|, which determines a set of possible values of 𝜔₀. We will determine the direction of motion of 𝜆=𝑖𝜔₀ as 𝜏 is varied, namely, d()Θ=signRe𝜆d𝜏_{𝜆=𝑖𝜔0}=signRed𝜆d𝜏⁻¹_{𝜆=𝑖𝜔0}. (3.40)

Theorem 3.3. Model system (2.1) undergoes Hopf bifurcation at the interior equilibrium 𝑃^∗ when tau=𝜏_𝑘, 𝑘=0,1,2,…. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium 𝑃^∗ when 𝜏>𝜏₀ and ‖𝜏−𝜏₀‖≪1.

Proof. Differentiating (3.12) with respect to 𝜏, we get d𝜆d𝜏⁻¹=4𝜆³+3𝑚₁𝜆²+2𝑚₂𝜆+𝑚₃𝜆⁴+𝑚₁𝜆³+𝑚₂𝜆²+𝑚₃𝜆+𝑚₄+2𝑛₂𝜆+𝑛₃𝑛₂𝜆²+𝑛₃𝜆−𝜏𝜆,Θ=signRed𝜆d𝜏⁻¹_{𝜆=𝑖𝜔0}=signRe4𝜆³+3𝑚₁𝜆²+2𝑚₂𝜆+𝑚₃𝜆⁴+𝑚₁𝜆³+𝑚₂𝜆²+𝑚₃𝜆+𝑚₄+2𝑛₂𝜆+𝑛₃𝑛₂𝜆²+𝑛₃𝜆−𝜏𝜆_{𝜆=𝑖𝜔0}=1𝜔²₀𝑚sign¹𝜔⁶₀+𝑚₁𝑚₂−2𝑚₃𝜔⁴₀+𝑚₂𝑚₃−3𝑚₁𝑚₄𝜔²₀+𝑚₃𝑚₄𝜔⁴₀−𝑚₂𝜔²₀+𝑚₄²+𝑚₃𝜔₀−𝑚₁𝜔³₀². (3.41) According to the values of 𝑚₁,𝑚₂,and𝑛₁ defined in (3.12), it is easy to show that 𝑚₁𝑚₂−2𝑚₃>0 and 𝑚₂𝑚₃−3𝑚₁𝑚₄>0. Consequently, it follows that sign[d(Re𝜆)/d𝜏]_{𝜆=𝑖𝜔0}>0, which implies there exists at least one eigenvalue with positive real part for 𝜏>𝜏₀, and the condition for Hopf bifurcation in the reference [29] is also satisfied yielding the required periodic solution.
Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium 𝑃^∗ when 𝜏>𝜏₀ and ‖𝜏−𝜏₀‖≪1.

4. Numerical Simulation

With the help of MATLAB, numerical simulations are provided to understand the theoretical results, which have been established in the previous sections of this paper.

4.1. Numerical Simulation for Optimal Harvesting Policy

In this subsection, values of parameters are taken from [9] which are used in Example 1 of [9] and set in appropriate units, 𝑟=5, 𝛼=4, 𝛽=1, 𝑎₁=3, 𝑎₂=2, 𝑟₁=0.1, 𝑟₂=0.1, 𝑑=2, 𝑝=13, 𝑞=0.18, 𝛼₀=0.08, 𝑐=1. For the model system (3.14), the range of the taxation 𝜎∈(0.5,7.0848) can be obtained based on (3.1) and (3.10). According to [10], we take the instantaneous annual rate of discount 𝛿=0.05 in appropriate units. According to the given values of parameters, (3.32) can be numerically computed and three roots can be obtained as follows: 𝑦^∗₂=0.5212,1.6323,2.2935. Based on (𝑝−𝜎)𝑞=𝑦^∗₂, the corresponding 𝜎_𝑖(𝑖=1,2,3) can be calculated, 𝜎₁=2.3405,𝜎₂=9.5965, and 𝜎₃=10.5777, respectively. It is obvious that only 𝜎₁=2.3405 satisfies the range (0.5,7.0848). Consequently, the optimal taxation is 𝜎_𝛿=𝜎₁=2.3405, then the optimal equilibrium levels of population and harvest effort can be also obtained (𝑥_𝛿,𝑦_1𝛿,𝑦_2𝛿,𝐸_𝛿)=(0.8591,0.4245,0.5212,5.6399), which are indicated in Figure 1.

Figure 1 

Dynamical responses of model system (3.14) with the optimal taxation.

4.2. Numerical Simulation for the Hopf Bifurcation

In this subsection, values of parameters are taken from [9], which are used in Example 3 of [9] and set in appropriate units, 𝑟=2, 𝛼=0.5, 𝛽=0.5, 𝑎₁=3, 𝑎₂=2, 𝑟₁=1, 𝑟₂=0.1, 𝑑=1, 𝑝=13, 𝑞=0.18, 𝛼₀=0.08, and𝑐=1. It follows from (3.1) and (3.10) that the taxation range is (0.5, 3.7407), and it can be obtained that population densities in model system (2.1) is (𝑥^∗,𝑦^∗₁,𝑦^∗₂)=(0.4,0.24,0.6) with 𝜎=2.3. It should be noted that 𝜎=2.3 is arbitrarily selected from the interval (0.5,3.7407), which can guarantee the existence of interior equilibrium of model system (2.1). Furthermore, it can be also calculated that 𝐵₃<0, which satisfies the assumption of the existence of a positive root for (3.37), and then 𝜏₀=0.8734 is calculated based on (3.39). By virtue of Butler's lemma [32], it can be concluded that the interior equilibrium remains stable for 𝜏<𝜏₀, which can be seen in Figure 2. It should be noted that 𝜏=0.3 is randomly selected in the interval (0,0.8734), which is enough to merit the above mathematical study.

Figure 2 

Dynamical responses of model system (2.1) with discrete time delay 𝜏=0.3.

According to Theorem 3.3 in this paper, a periodic solution caused by the phenomenon of Hopf bifurcation and a limit cycle corresponding to this periodic solution occurs as 𝜏 increases through 𝜏₀, which are shown in Figures 3 and 4, respectively.

Figure 3 

Dynamical responses of model system (2.1) with discrete time delay 𝜏=0.8734.

Figure 4 

A limit cycle for the model system (2.1) corresponding to the periodic solution shown in Figure 3.

5. Conclusion

In this paper, a bioeconomic model is proposed to investigate dynamics of the effects of a stage-structured prey-predator system with harvesting effort and gestation delay. Theoretical analysis shows that the interior equilibrium is locally asymptotically stable around interior equilibrium when the model system is in absence of discrete time delay. By using Pantryagin's Maximum Principle, an optimal harvesting policy with taxation is derived to ensure the sustainable development of biological resource and prosperous commercial harvesting. It reveals that the user's total cost of harvest per unit effort must be equal to the discounted value of the future price at the steady state level. In the case of gestation delay, the stability analysis reveals that gestation delay is responsible for the stability switch of model system. A phenomenon of Hopf bifurcation occurs as the discrete time delay increases through a certain threshold.

It should be noted that taxation instrument is discussed to control overharvesting from prey-predator system with gestation delay in [26]. Compared with work done in [26], stage structure of predator population is considered, and the periodic orbit within the interior of the first quadrant of state space around interior equilibrium is also investigated in this paper. The work done in [9] is extended by incorporating the harvesting effort into the prey-predator system, and taxation is adopted as a controlling instrument to regulate harvesting of predator. From the qualitative analysis of the model, the effect of harvesting effort is extensively investigated with and without discrete time delay. Compared with model system investigated in this paper, the harvesting effort is not considered in [9], the interior equilibrium becomes unstable as 𝜏=0.7. However, the interior equilibrium of the model system (2.1) remains stable as 𝜏=0.7 based on the analysis in this paper. It implies that the harvesting effort has an effect of stabilizing the interior equilibrium, and the cyclic behavior can be prevented by applying the harvesting effort into the model system.

Acknowledgments

This work is supported by National Natural Science Foundation of China, Grant no. 60974004 and by Liaoning Provincial Foundation of Science and Technology, Grant no. 20082023.