A Mathematical Model for Optimal Management and Utilization of a Renewable Resource by Population

Abstract

A dynamical model is proposed and analyzed to study the effect of the population on the resource biomass by taking into account the crowding effect. Biological and bionomical equilibria of the system are discussed. The global stability behavior of the positive equilibrium is studied via the output feedback control. An appropriate Hamiltonian function is formed for the discussion of optimal harvesting of resource which is utilized by the population using Pontryagin's Maximum Principal. A numerical simulation is performed on the model to analyze the theoretical results.

1. Introduction

Renewable resources are a very important source of food and materials which are essential for the growth and survival of the biological population. The continuous and unplanned use of these resources may lead to the extinction of resources and thereby affecting the survival of resource-dependent species. There has been a considerable interest in the modeling of renewable resources such as fishery and forestry. Dynamic models for the commercial fishing have been studied extensively taking into account the economic and ecological factors [1, 2]. Based on the works of Clark [1, 2], several investigations have been conducted [3–12]. Leung and Wang [13] proposed a simple economic model and investigated the phenomena of nonexplosive fishing capital investment and nonextinctive fishery resources. Chaudhuri [14] proposed a model for two competing fish species each of which grows logistically. He examined the stability analysis and discussed the bionomical equilibrium and optimal harvesting policy. It was also shown [14] that there is no limit cycle in the positive quadrant. Ragozin and Brown [15] proposed a model in which the prey has no commercial value and the predator is selectively harvested. When the prey and predator both are harvested, then an optimal policy for maximizing the present value and an estimation to the true loss of resource value due to catastrophic fall in stock level has been discussed in detail by Mesterton-Gibbons [16]. In another paper, Mesterton-Gibbons [17] proposed a Lotka-Volterra model of two independent populations and studied an optimal harvesting policy. Fan and Wang [18] generalized the classical model of Clark [1, 2] by considering the time-dependent the Logistic equation with periodic coefficients and they showed that their model has a unique positive periodic solution, which is globally asymptotically stable for positive solutions. They also investigated the optimal harvesting policies for the constant harvest and the periodic harvest. The optimal harvesting policy of a stage structure problem was studied by Zhang et al. [19]. They found conditions for the coexistence and extinction of species. Song and Chen [20] proposed two-species competitive system and discussed the local and global stability analysis of the positive equilibrium point of the system. They have also discussed the optimal harvesting policy for the mature population. Dubey et al. [10] proposed a model where the fish population partially depends on a resource and is harvested. They examined stability analysis and the optimal harvesting policy with taxation as a control variable. Dubey et al. [11] discussed a model of a fishery resource system in an aquatic environment which was divided into two zones—the free fishing zone and the reserved zone. They discussed biological and bionomical equilibria and optimal harvesting policy. Dubey et al. [12] also proposed and analyzed an inshore-offshore fishery model where the fish population is being harvested in both areas. Then they investigated the stability analysis and optimal harvesting policy by taking taxation as a control instrument. Kar et al. [21] considered a two- prey one-predator model where both the preys grow logistically and harvested. Again Kar and Chottopadhayay [22] described a single species model, which has two stages: (i) a mature stage (ii) an immature stage; they discussed the existence of equilibrium points and their stability analysis. They proved that the optimal harvesting policy is much superior to the MSY policy and optimal paths always take less time than the suboptimal path to reach the optimal steady state. Kar et al. [23] proposed a prey-predator model with a nonmonotonic functional response and both the species are harvested. To obtain the strategies for the management of the system, they used the harvesting effort as a control variable. They found the stability condition of an interior equilibrium point in terms of harvesting effort and proved that there exists a super critical Hopf bifurcation.

From the above literature and to the best of our knowledge, it appears that harvesting of a renewable resource which is being utilized by a population for its own growth and development has not been considered by taking into account the effect of crowding. Hence, in this paper, we propose a mathematical model for the biological population which is being partially dependent on a renewable resource. This resource is further harvested for the development of the society.

The organization of the paper is as follows. Section 2 describes the development of the model and Section 3 gives a detailed outline of the stability analysis of the system. The bionomical equilibrium and the maximum sustainable yield are presented in Sections 4 and 5, respectively. The output feedback control is given in Section 6 and the optimal harvesting policy in Section 7. A numerical simulation experiment has been presented in Sections 8 and 9 is the concluding remarks followed by references.

2. Mathematical Model

Let us consider a renewable resource growing logistically in a habitat. Then the dynamics of this resource biomass is governed by𝑑𝐵𝑑𝑡=𝑎₀𝐵−𝑎₁𝐵²=𝑎₀𝐵𝐵1−𝑀, (1a) where 𝐵(𝑡) is the density of resource biomass, 𝑎₀ is its intrinsic growth rate, 𝑎₁ is an interspecific interference coefficient, and 𝑀=𝑎₀/𝑎₁ is the carrying capacity.

Let 𝑁(𝑡) be the density of the population at any time 𝑡≥0 which utilizes the resource biomass 𝐵(𝑡) for its own growth and development. Thus, the intrinsic growth rate and carrying capacity of the resource biomass will depend on the density of population. Hence, we assume that 𝑎₀ and 𝑀 in (1a) are functions of 𝑁. Thus, (1a) reduces to 𝑑𝐵𝑑𝑡=𝑎₀(𝑁)𝑓𝐵−⁰𝐵²𝑀(𝑁), (1b) where 𝑓₀ is a positive constant.

We consider the following assumptions.(i)The intrinsic growth rate 𝑎₀(𝑁) is a decreasing function of 𝑁and it satisfies 𝑎₀(0)=𝑟>0,𝑎^′₀(𝑁)≤0for𝑁≥0. (2)  We take a particular form of 𝑎₀(𝑁) as 𝑎₀(𝑁)=𝑟−𝛼₁𝑁. (3)(ii)The carrying capacity 𝑀(𝑁) is also a decreasing function of 𝑁 and it satisfies 𝑀(0)=𝑚₀>0,𝑀^′(𝑁)≤0for𝑁≥0. (4)  We take a particular form of 𝑀(𝑁) as 𝑀(𝑁)=𝑚₀1+𝑚₁𝑁. (5)

Let us denote 𝐾=𝑟𝑚₀/𝑓₀ and 𝛼₂=𝑚₁𝑓₀/𝑚₀.

 Then (1b) can be rewritten as 𝑑𝐵𝐵𝑑𝑡=𝑟𝐵1−𝐾−𝛼₁𝑁𝐵−𝛼₂𝑁𝐵². (6)

Now, we consider a population of density 𝑁(𝑡) which is also growing logistically. We assume that the growth rate and carrying capacity of the population depend on the resource biomass density. If the resource biomass density increases, the growth rate and carrying capacity of population also increase. Thus, the dynamics of population is governed by the following differential equation: 𝑑𝑁𝑁𝑑𝑡=𝑠𝑁1−𝐿+𝛽₁𝑁𝐵+𝛽₂𝑁𝐵², (7) where 𝑠 is the intrinsic growth rate of population, 𝐿 is its carrying capacity in the absence of resource biomass, 𝛽₁ and 𝛽₂ are the growth rates of population in the presence of resource biomass.

We assume that the resource biomass is harvested with the harvesting rate ℎ(𝑡)=𝑞𝐸𝐵, where 𝑞 is a positive constant and in fishery resource it is known as catchability coefficient, and 𝐸 is the harvesting effort, which is a control variable.

Keeping the above aspect in view, the dynamics of the system can be governed by the system of the following differential equations:𝑑𝐵𝐵𝑑𝑡=𝑟𝐵1−𝐾−𝛼₁𝑁𝐵−𝛼₂𝑁𝐵²−𝑞𝐸𝐵,𝑑𝑁𝑁𝑑𝑡=𝑠𝑁1−𝐿+𝛽₁𝑁𝐵+𝛽₂𝑁𝐵²,𝐵(0)>0,𝑁(0)>0. (8a)(8b)

In Section 3, we present the stability analysis of model (8a)-(8b).

3. Stability Analysis

First of all, we state the following lemma which is a region of attraction for the model system (8a)-(8b).

Lemma 1. The set Ω=(𝐵,𝑁)∈𝑅⁺₂∶0≤𝐵≤𝐾,0≤𝑁≤𝐿₀ (9) is a region of attraction for all solutions initiating in the interior of the positive quadrant, where 𝐿₀=𝐿𝑠𝑠+𝛽₁𝐾+𝛽₂𝐾². (10)

The above lemma shows that all the solutions of model (8a)-(8b) are nonnegative and bounded, and hence the model is biologically well behaved.

The proof of this lemma is similar to Freedman and So [24], Shukla and Dubey [25], hence omitted.

Now, we discuss the equilibrium analysis of the present model. It can easily be checked that model (8a)-(8b) has four non-negative equilibria, namely, 𝑃₀(0,0),𝑃₁𝐵,0,𝑃₂(0,𝐿),𝑃^∗𝐵^∗,𝑁^∗. (11)

The equilibrium points 𝑃₀ and 𝑃₂ always exist. For the existence of equilibrium point 𝑃₁, we note that 𝐵 is given by𝐾𝐵=𝑟𝑟−𝑞𝐸, (12a)

This shows that 𝑃₁ exists if 𝑟>𝑞𝐸. (12b)

For the fourth equilibrium point 𝑃^∗, we note that 𝐵^∗ and 𝑁^∗ are the positive solutions of the following algebraic equations:𝑟𝐵1−𝐾−𝛼₁𝑁−𝛼₂𝑠𝑁𝑁𝐵−𝑞𝐸=0,1−𝐿+𝛽₁𝐵+𝛽₂𝐵²=0. (13a)(13b)

Substituting the value of 𝑁, from (13b) into (13a), we get a cubic equation in 𝐵 as follows: 𝑎₁𝐵³+𝑎₂𝐵²+𝑎₃𝐵+𝑎₄=0, (14) where 𝑎₁=𝐾𝐿𝛼₂𝛽₂𝑎>0,₂=𝐾𝐿𝛼₂𝛽₁+𝛼₁𝛽₂𝑎>0,₃=𝑠𝑟+𝐾𝐿𝛼₂𝑠+𝛼₁𝛽₁𝑎𝐾𝐿>0,₄=𝛼₁𝐾𝐿𝑠+𝑞𝐸−𝑟𝐾𝑠. (15)

The above equation (14) has a unique positive solution 𝐵=𝐵^∗, if the following inequality holds: 𝛼₁𝐿+𝑞𝐸<𝑟. (16)

After knowing the value of 𝐵^∗, the value of 𝑁^∗ can then be calculated from the relation 𝑁^∗=𝐿𝑠𝑠+𝛽₁𝐵^∗+𝛽₂𝐵^∗. (17)

From (16), we note that for the coexistence of resource biomass and population together in a habitat, the intrinsic growth rate of the resource biomass must be larger than a threshold value. This threshold value depends upon the carrying capacity of the population and the harvesting effort.

Now we discuss the local and global stability behavior of these equilibrium points. For local stability analysis, first we find variational matrices with respect to each equilibrium point. Then by using eigenvalue method and the Routh-Hurwitz criteria, we can state Theorems 2 to 5.

Theorem 2. (i)  If the equilibrium point 𝑃₁(𝐵,0) exists, then 𝑃₀ is always unstable in the 𝐵-𝑁 plane.
(ii)  If the equilibrium point 𝑃₁(𝐵,0) does not exist, then 𝑃₀ is a saddle point with stable manifold locally in the 𝐵-direction and with unstable manifold locally in the 𝑁-direction.

Theorem 3. The equilibrium point 𝑃₁(𝐵,0), whenever it exists, is a saddle point with stable manifold locally in the 𝐵-direction and with unstable manifold locally in the 𝑁-direction.

Theorem 4. (i)  If 𝑟>(𝛼₁𝐿+𝑞𝐸), then 𝑃₂(0,𝐿) is a saddle point with unstable manifold locally in the 𝐵-direction and stable manifold locally in the 𝑁-direction.
(ii)  If 𝑟<(𝛼₁𝐿+𝑞𝐸), then 𝑃₂(0,𝐿) is always stable in the 𝐵-𝑁 plane.

Theorem 5. The interior equilibrium 𝑃^∗, whenever it exists, is always locally asymptotically stable in the 𝐵-𝑁 plane.

In the next theorem, we are able to find a sufficient condition for 𝑃^∗ to be globally asymptotically stable.

Theorem 6. Let the following inequality holds in Ω:𝛼₁+𝛽₁+𝛼₂+𝛽₂𝐾+𝛽₂𝐵^∗2<4𝑠𝐿𝑟𝐾+𝛼₂𝑁^∗. (18)
Then the interior equilibrium 𝑃^∗ is globally asymptotically stable with respect to all solutions initiating in the interior of the region Ω defined in Lemma 1.

Proof of Theorem 6. Proof of this theorem is given in the Appendix.

Our next result shows that the model under consideration cannot have any closed trajectories in the interior of the first quadrant.

Theorem 7. System (8a)-(8b) cannot have any limit cycle in the interior of the positive quadrant.

Proof. Let 𝐻(𝐵,𝑁)=1/𝐵𝑁>0 is a continuously differential function in the interior of the positive quadrant of 𝐵-𝑁 plane.
Let 𝐹₁()𝐵𝐵,𝑁=𝑟𝐵1−𝐾−𝛼₁𝑁𝐵−𝛼₂𝑁𝐵²𝐹−𝑞𝐸𝐵,₂()𝑁𝐵,𝑁=𝑠𝑁1−𝐿+𝛽₁𝑁𝐵+𝛽₂𝑁𝐵². (19) Then, Δ()=𝜕𝐵,𝑁𝜕𝐵𝐻𝐹₁+𝜕𝜕𝑁𝐻𝐹₂1=−𝑁𝑟𝐾+𝛼₂𝑁−𝑠𝐵𝐿<0. (20)

This shows that Δ(𝐵,𝑁) does not change sign and is not identically zero in the positive quadrant of the 𝐵-𝑁plane. By Bendixon-Dulac criteria, it follows that the system (8a)-(8b) has no closed trajectory, and hence no periodic solution in the interior of the positive quadrant of the 𝐵-𝑁 plane.

4. Bionomical Equilibrium

In this section, we study the bionomical equilibrium of the model system (8a)-(8b). Bionomical equilibrium is the level at which the total revenue (TR) obtained by selling the harvested biomass in an economic equilibrium case is equal to the total cost (TC) to the harvested biomass, that is, the economic rent is completely dissipated.

The net economic revenue at time 𝑡 is given by 𝜋(𝐵,𝐸,𝑡)=𝑝𝑞𝐵−𝑐𝐸, (21) where 𝑝 is price per unit biomass and 𝑐 is utilized cost per unit resource biomass. The bionomical equilibrium is 𝑃_∞(𝐵_∞,𝑁_∞,𝐸_∞), where 𝐵_∞,𝑁_∞, and 𝐸_∞ are the positive solutions of𝐵^•=𝑁^•=𝜋=0. (22a)

Solving (22a), we get 𝐵_∞=𝑐,𝑁𝑝𝑞_∞=𝐿𝑠𝛽𝑠+¹𝑐+𝛽𝑝𝑞²𝑐²𝑝²𝑞²𝐸_∞,=𝑟𝑞𝑐1−−𝛼𝐾𝑝𝑞¹+𝑐𝛼₂𝑁𝑝𝑞𝑞^∞. (22b)(22c)(22d)

It is clear that 𝐸_∞𝑟>0,if𝑞𝑐1−>𝛼𝐾𝑝𝑞¹+𝑐𝛼₂𝑁𝑝𝑞𝑞^∞. (23)

Thus the bionomical equilibrium 𝑃_∞(𝐵_∞,𝑁_∞,𝐸_∞) exists under condition (23).

Remark 8. From (22b) and (23), it may be noted that 𝐵_∞=(𝑐/(𝑝𝑞))<𝐾.

5. The Maximum Sustainable Yield

The maximum sustainable yield (MSY) of any biological resource biomass is the maximum rate at which it can be harvested and any larger harvest rate will lead to the depletion of resource eventually to zero. In the absence of any population, the value of MSY is given by [1] ℎ⁰_MSY=𝑟𝐾4. (24)

If the resource biomass is subjected to the harvesting by a population, the sustainable yield is given by ℎ=𝑞𝐸𝐵^∗=𝑟𝐵^∗𝐵1−^∗𝐾−𝛼₁𝑁^∗𝐵^∗−𝛼₂𝑁^∗𝐵^∗2. (25)

We note that 𝜕ℎ𝜕𝐵^∗=0gives𝐵^∗=𝐾𝑟−𝛼₁𝑁^∗2𝑟+𝛼₂𝐾𝑁^∗,𝜕²ℎ𝜕𝐵^∗2<0. (26)

Thus, ℎ_MSY=𝐵^∗((𝑟𝐵^∗/𝐾)+𝛼₂𝑁^∗𝐵^∗), when 𝐵^∗=𝐾(𝑟−𝛼₁𝑁^∗)/2(𝑟+𝛼₂𝐾𝑁^∗).

From the above equations, it is interesting to note that when 𝑁^∗=0, then 𝐵^∗=𝐾/2 and ℎ_MSY=𝑟𝐾/4=ℎ⁰_MSY.

This result matches to Clark [1].

If ℎ>ℎ_MSY, then it denotes the overexploitation of the resource and if ℎ<ℎ_MSY, then the resource biomass is under exploitation.

6. Output Feedback Control

The habitat under our consideration consists of a resource biomass which is utilized by a population. The resource biomass is being harvested and the harvesting effort 𝑢=𝐸 is considered as an input. The harvesting effort is applied to the stock and produces a yield 𝑌=𝑞𝐵 per unit effort. We assume that the total yield 𝑌 per unit effort is subject to the constraint 𝑌_min<𝑌<𝑌_max. Then our objective is to construct an output feedback control 𝐸(𝑌)=𝐸⁰+𝑢(𝑌(𝑡)) in such a way that the steady state 𝑃^∗(𝐵^∗,𝑁^∗) is globally asymptotically stable for the closed-loop system. Then model equations (8a)-(8b) can be written as in the vector matrix differential equation form 𝑋^•=𝑓(𝑋)+𝑢𝑔(𝑋), (27) where 𝑢=𝐸(𝑡), 𝑓(𝑋)=𝐵𝑟𝐵1−𝐾−𝛼₁𝑁𝐵−𝛼₂𝑁𝐵²𝑁𝑠𝑁1−𝐿+𝛽₁𝑁𝐵+𝛽₂𝑁𝐵²,𝑔(𝑋)=0.−𝑞𝐵 (28)

Now, under an analysis similar to Louartassi et al. [26] and Mazoudi et al. [27], one can prove the following theorem.

Theorem 9. For any constant fishing effort 𝐸⁰, there exists a 𝑤>0 such that 𝑃^∗ is globally asymptotically stable through the output feedback control law 𝑢=𝑤(𝑌−𝑌^∗), where 𝑌^∗=𝑞𝐵^∗.

7. Optimal Harvesting Policy

In this section, we discuss the optimal management of a renewable resource in the presence of population which is to be adopted by the regulatory agencies to protect the resource and to ensure the survival of the population with a sustainable development. The present value 𝐽 of a continuous time stream of revenues is given by 𝐽=^∞₀𝑒^−𝛿𝑡𝑝𝑞𝐵(𝑡)−𝑐𝐸(𝑡)𝑑𝑡, (29) where 𝛿 is the instantaneous rate of annual discount. Thus our objective is to ()-()max𝐽,subjecttothestateequations8a8bandtothecontrolconstraints0≤𝐸≤𝐸_max. (30)

For this purpose, we use Pontryagin’s Maximum Principle. The associated Hamiltonian function is given by 𝐻=𝑒^−𝛿𝑡𝑝𝑞𝐵−𝑐𝐸+𝜆₁(𝑡𝐵)𝑟𝐵1−𝐾−𝛼₁𝑁𝐵−𝛼₂𝑁𝐵²−𝑞𝐸𝐵+𝜆₂(𝑡𝑁)𝑠𝑁1−𝐿+𝛽₁𝑁𝐵+𝛽₂𝑁𝐵², (31) where 𝜆₁ and 𝜆₂ are adjoint variables and 𝜎(𝑡)=𝑒^−𝛿𝑡(𝑝𝑞𝐵−𝑐)−𝜆₁𝑞𝐵 is called a switching function.

The optimal control 𝐸(𝑡) which maximizes 𝐻 must satisfy the following conditions:𝐸=𝐸_max,when𝜎(𝑡)>0,thatis,𝜆₁𝑒^𝛿𝑡𝑐<𝑝−,(𝑡)𝑞𝐵𝐸=0,when𝜎<0,thatis,𝜆₁𝑒^𝛿𝑡𝑐>𝑝−.𝑞𝐵 (32a)(32b)

The usual shadow price is 𝜆₁𝑒^𝛿𝑡 and the net economic revenue on a unit harvest is 𝑝−(𝑐/𝑞𝐵). This shows that if the shadow price is less than the net economic revenue on a unit harvest, then 𝐸=𝐸_max and if the shadow price is greater than the net economic revenue on a unit harvest, then 𝐸=0. When the shadow price equals the net economic revenue on a unit harvest, that is, 𝜎(𝑡)=0, then the Hamiltonian 𝐻 becomes independent of the control variable 𝐸(𝑡), that is, 𝜕𝐻/𝜕𝐸=0.

This is a necessary condition for the singular control 𝐸^∗(𝑡) to be optimal over control set 0<𝐸^∗<𝐸_max.

Thus, the optimal harvesting policy is 𝐸(𝑡)=𝐸_max(𝑡)(𝑡)𝐸,𝜎>00,𝜎<0_∗,𝜎(𝑡)=0, (33) when 𝜎(𝑡)=0, we have 𝜆₁=𝑒^−𝛿𝑡𝑐𝑝−.𝑞𝐵 (34)

Now, in order to find the path of singular control, we utilize the Pontrygin’s Maximum Principle. According to this principle, the adjoint variables 𝜆₁ and 𝜆₂ must satisfy 𝑑𝜆₁𝑑𝑡=−𝜕𝐻,𝜕𝐵𝑑𝜆₂𝑑𝑡=−𝜕𝐻.𝜕𝑁 (35)

The above equations can be rewritten as 𝑑𝜆₁𝑒𝑑𝑡=−^−𝛿𝑡𝑝𝑞𝐸+𝜆₁𝑟1−2𝐵𝐾−𝛼₁𝑁−2𝛼₂𝑁𝐵−𝑞𝐸+𝜆₂𝛽₁𝑁+2𝛽₂𝑁𝐵,𝑑𝜆₂𝑑𝑡=𝜆₁𝛼₁𝐵+𝛼₂𝐵²−𝜆₂𝑠1−2𝑁𝐿+𝛽₁𝐵+𝛽₂𝐵². (36)

Using (13a)-(13b), the above equations can be rewritten as 𝑑𝜆₁𝑒𝑑𝑡=−^−𝛿𝑡𝑝𝑞𝐸−𝜆₁𝑟𝐵𝐾+𝛼₂𝑁𝐵+𝜆₂𝛽₁𝑁+2𝛽₂𝑁𝐵,𝑑𝜆₂𝑑𝑡=𝑒^−𝛿𝑡𝑐𝑝−𝛼𝑞𝐵¹𝐵+𝛼₂𝐵²+𝜆₂𝑠𝑁𝐿. (37)(38)

Equation (38) can again be written as 𝑑𝜆₂𝑑𝑡−𝐴₁𝜆₂=−𝐵₁𝑒^−𝛿𝑡, (39) where 𝐴₁=𝑠𝑁𝐿,𝐵₁𝑐=−𝑝−𝛼𝑞𝐵¹𝐵+𝛼₂𝐵². (40)

Solving (38), we get 𝜆₂=𝐵₁𝐴₁𝑒+𝛿^−𝛿𝑡+𝐾₀𝑒^𝐴1𝑡, (41)

We note that when 𝑡→∞, then the shadow price 𝜆₂𝑒^𝛿𝑡 is bounded if 𝐾₀=0.

Thus, we have 𝜆₂=𝐵₁𝐴₁𝑒+𝛿^−𝛿𝑡. (42)

Now from (37), we get 𝑑𝜆₁𝑑𝑡−𝐴₂𝜆₁=−𝐵₂𝑒^−𝛿𝑡, (43) where 𝐴₂=−𝑟𝐵𝐾−𝛼₂𝐵𝑁𝐵,₂𝐵=𝑝𝑞𝐸−¹𝐴₁𝛽+𝛿¹𝑁+2𝛽₂𝑁𝐵. (44)

Equation (43) yields 𝜆₁=𝐵₂𝐴₂𝑒+𝛿^−𝛿𝑡+𝐾₁𝑒^𝐴2𝑡. (45)

Again for the shadow price 𝜆₁𝑒^𝛿𝑡 to be bounded, we must have 𝐾₁=0.

Thus, 𝜆₁=𝐵₂𝐴₂𝑒+𝛿^−𝛿𝑡. (46)

Equations (34) and (46) yield 𝐵₂𝐴₂=𝑐+𝛿𝑝−.𝑞𝐵 (47)

Substituting the values of 𝐴₁,𝐵₁,𝐴₂, and 𝐵₂into (47), we get 1𝐸=𝑐𝑝𝑞𝑝−𝑞𝐵𝛿+𝑟𝐵𝐾+𝛼₂+𝛼𝑁𝐵₁𝐵+𝛼₂𝐵²𝛽₁𝑁+2𝛽₂𝑁𝐵)).((𝑠𝑁/𝐿+𝛿 (48)

Hence, solving (13a)-(13b) with the help of (48), we get an optimal solution (𝐵_𝛿,𝑁_𝛿) and the optimal harvesting effort 𝐸=𝐸_𝛿.

8. Numerical Simulations

In order to investigate the dynamics of the model system (8a)-(8b) with the help of computer simulations, we choose the following set of values of parameters (other set of parameters may also exist): 𝑟=1.6,𝑠=1.2,𝐾=100,𝐿=100,𝐸=100,𝑝=0.1,𝑞=0.01,𝛼₁𝛼=0.001₂=0.01,𝛽₁=0.01,𝛽₂=0.1, (49) with initial condition: 𝐵(0)=5 and 𝑁(0)=8.

For this set of parameters, condition (16) for the existence of the interior equilibrium is satisfied. This shows that 𝑃^∗(𝐵^∗,𝑁^∗) exists and it is given by 𝐵^∗=0.4790,𝑁^∗=102.3108. (50)

Thus the interior equilibrium point 𝑃^∗(𝐵^∗,𝑁^∗) is locally asymptotically stable. It is noted that condition (18) in Theorem 6 is not satisfied for the set of parameters chosen in (49). Since condition (18) is just a sufficient condition for 𝑃^∗ to be globally asymptotically stable, no conclusion can be drawn at this stage. The behavior of 𝐵 and 𝑁 with respect to time 𝑡 is plotted in Figure 1 for the set of values of parameter chosen in (49). From this figure, we see that the density of population increases, whereas the density of resource biomass decreases with respect to time, and both settle down at its equilibrium level.

Figure 1 

Plot of 𝐵 and 𝑁 versus 𝑡 for the values of parameters given in (49).

Now, we choose another set of values of parameters as follows: 𝑟=1.6,𝑠=1.2,𝐾=100,𝐿=100,𝐸=100,𝑝=0.1,𝑞=0.01,𝛼₁𝛼=0.001₂=0.0001,𝛽₁=0.01,𝛽₂=0.0001 (51) with different initial values.

For the set of values of parameters given in (51), it may be noted that the positive equilibrium 𝑃^∗(𝐵^∗,𝑁^∗) exists and it is given by 𝐵^∗=17.4311,𝑁^∗=117.0579. (52)

It may be cheeked that in this case (for the values of parameters given in (51)), condition (18) in Theorem 6 is satisfied. This shows that 𝑃^∗(𝐵^∗,𝑁^∗) is locally as well as globally asymptotically stable in the interior of the first quadrant. In Figure 2, we have plotted the behavior of 𝐵 and 𝑁 with different initial values. Figure 2 shows that all the trajectories starting from different initial points converge to the point 𝑃^∗(17.4311,117.0579). This shows that 𝑃^∗(17.4311,117.0579) is globally asymptotically stable.

Figure 2 

Global stability of 𝑃∗(𝐵∗,𝑁∗).

It may be noted here that 𝛼₂ and 𝛽₂ are important parameters governing the dynamics of the system. Therefore, we have plotted the behavior of 𝐵 and 𝑁 with respect to time for different values of 𝛼₂ in Figure 3 and for different values of 𝛽₂ in Figure 4.

(a)

(b)

Figure 3 

Plot of 𝐵 and 𝑁 with respect to time 𝑡 for different values of 𝛼2, and other values of parameters are the same as given in (51).

(a)

(b)

Figure 4 

Plot of 𝐵 and 𝑁 with respect to time for different values of 𝛽2 with 𝛼2=0.0001 and other values of parameters are the same as in Figure 2.

From Figure 3(a), we note that 𝐵 decreases as 𝛼₂ increases. If 𝛼₂ is very small, then 𝐵 initially increases slightly and then decreases and converges to its equilibrium level. If 𝛼₂ increases beyond a threshold value, then 𝐵 always decreases and settles down at its equilibrium level. Figure 3(b) shows that 𝑁 also decreases as 𝛼₂ decreases. This is due to the fact that with the increase in 𝛼₂, the equilibrium level of 𝐵 decreases and since the population 𝑁 is dependent on the resource biomass 𝐵, thus 𝑁 also decreases. However, it is interesting to note here that 𝑁 always increases with respect to time 𝑡 and finally attains its equilibrium level.

From Figure 4(a), it is noted that 𝐵 decreases as 𝛽₂ increases. Figure 4(b) shows that 𝑁 increases as 𝛽₂ increases. It may also be noted that for all positive values of 𝛽₂, 𝑁 increases continuously with time and finally settles down at its steady state. However, the resource biomass 𝐵 initially increases for some time, then decreases continuously and finally gets stabilized at its lower equilibrium level. This shows that if the population utilizes the resource without any control, then the resource biomass decreases continuously and it may be doomed to extinction.

In order to see the qualitative behavior of the optimal harvesting resource, we solve (13a), (13b), and (48) for the same set of values of parameters as given in (51) with the additional values as 𝑐=0.001 and 𝛿=5. Then we get the optimal equilibrium levels as given below: 𝐵_𝛿=1.40895,𝑁_𝛿=101.191,𝐸_𝛿=146.201. (53)

We observe that 𝐸 is a very important parameter which governs the dynamics of the system. The behavior of 𝐵 with respect to time is shown in Figure 5 and the behavior of 𝑁 with respect to time is shown in Figure 6 for different values of 𝐸. Figure 5 shows that if 𝐸 is less than its optimal level (𝐸_𝛿=146.201), then 𝐵 increases and settles down at its equilibrium level. If 𝐸>𝐸_𝛿, then 𝐵 decreases and settles down at a lower equilibrium level. If 𝐸 is large enough, then 𝐵 tends to zero. Figure 6 shows that if 𝐸<𝐸_𝛿, and then 𝑁 increases with time and finally obtains its equilibrium level. If 𝐸>𝐸_𝛿, then 𝑁 again increases with time, but it settles down at a lower equilibrium level. This is due to the fact that when 𝐸>𝐸_𝛿, then 𝐵 decreases and hence 𝑁 also decreases. This suggests that the harvesting effort 𝐸 should be always kept less than 𝐸_𝛿 (optimal harvesting level) so that both the resource and the population can be maintained at an optimal level.

Figure 5 

Plot of 𝐵 versus 𝑡 for different values of 𝐸.

Figure 6 

Plot of 𝑁 versus 𝑡 for different values of 𝐸.

9. Conclusions

In this paper, a mathematical model has been proposed and analyzed to study the effect of harvesting of a renewable resource. It has been assumed that the resource is being utilized by a population for its own growth and survival. The resource and the population both are growing logistically. The existences of equilibria and stability analysis have been discussed with the help of the stability theory of ordinary differential equations. It has been shown that the positive equilibrium 𝑃^∗ (whenever it exists) is always locally asymptotically stable. But 𝑃^∗ is not always globally asymptotically stable. However, we have found a sufficient condition under which 𝑃^∗ is globally asymptotically stable. This condition gives a threshold value of the specific growth rate (𝑟) of resource biomass. If 𝑟 is larger than this threshold value, then 𝑃^∗ is globally asymptotically stable. Using Bendixon-Dulac criteria, it has been observed that the model system has no limit cycle in the interior of the positive quadrant. We have also discussed the bionomical equilibrium of the model and it has been observed that the bionomical equilibrium of the resource biomass does not depend upon the growth rate and carrying capacity of the population. An analysis for sustainable yield (ℎ) and maximum sustainable yield (ℎ_MSY) has been carried out. It has been shown that if ℎ>ℎ_MSY, then the resource biomass will tend to zero and if ℎ<ℎ_MSY, then the resource biomass and population may be maintained at desired level. The global asymptotic stability behavior of the positive equilibrium has also been studied through the output feedback control method.

We have constructed a Hamiltonian function and then using Pontryagin’s Maximum Principle, the optimal harvesting policy has been discussed. An optimal equilibrium solution has been obtained. The results are validated based on the numerical simulation. Here we observed that 𝛼₂,𝛽₂ and 𝐸 are important parameters governing the dynamics of the system. It has been found that if 𝛼₂ increases, then the resource biomass and the population both decrease. However, if 𝛽₂ increases, then the population density increases and the resource biomass density decreases. A threshold value for the optimal harvesting effort 𝐸_𝛿 has been found theoretically as well as numerically. It has been shown that the harvesting effort 𝐸 should be always kept less than 𝐸_𝛿 to maintain the resource and the population at an optimal equilibrium level.

Appendix

Proof of Theorem 6. We consider the following positive definite function about 𝑃^∗: 𝑉=𝐵−𝐵^∗−𝐵^∗𝐵ln𝐵^∗+𝑁−𝑁^∗−𝑁^∗𝑁ln𝑁^∗. (A.1)
Differentiating 𝑉 with respect to time along the solutions of model (8a)-(8b), a little algebraic manipulation yields 𝑑𝑉1𝑑𝑡=−2𝑎₁₁𝐵−𝐵^∗2+𝑎₁₂𝐵−𝐵^∗𝑁−𝑁^∗−12𝑎₂₂𝑁−𝑁^∗2, (A.2) where 𝑎₁₁𝑟=2𝐾+𝛼₂𝑁^∗,𝑎₂₂=2𝑠𝐿,𝑎₁₂=−𝛼₁+𝛽₁−𝛼₂𝐵+𝛽₂𝐵+𝐵^∗. (A.3)
Sufficient condition for 𝑑𝑉/𝑑𝑡 to be negative definite is that the following inequality holds: 𝑎²₁₂<𝑎₁₁𝑎₂₂. (A.4)
We note that (18) implies (A.4), thus 𝑉 is Liapunov’s function for all solutions initiating in the interior of the positive quadrant whose domain contains the region of attraction Ω, proving the theorem.

Acknowledgments

The authors gratefully acknowledge the support received from UGC, New Delhi, Grant no. F.510/2/DRS/2010 (SAP-I). A. Patra also acknowledges the financial assistance of DST (India) under WOS-A scheme for carrying out this work.