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

Authors: B. Dubey; Atasi Patra
Journal: Journal of Mathematics
 (2013)

---

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


---

## Body

## 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. Section2 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(1a)dBdt=a0B-a1B2=a0B(1-BM),
where B(t) is the density of resource biomass, a0 is its intrinsic growth rate, a1 is an interspecific interference coefficient, and M=a0/a1 is the carrying capacity.LetN(t) be the density of the population at any time t≥0 which utilizes the resource biomass B(t) 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 a0 and M in (1a) are functions of N. Thus, (1a) reduces to
(1b)dBdt=a0(N)B-f0B2M(N),
where f0 is a positive constant.We consider the following assumptions.(i)
The intrinsic growth ratea0(N) is a decreasing function of Nand it satisfies
(2)a0(0)=r>0,a0′(N)≤0forN≥0.
We take a particular form of a0(N) as
(3)a0(N)=r-α1N.(ii)
The carrying capacityM(N) is also a decreasing function of N and it satisfies
(4)M(0)=m0>0,M′(N)≤0forN≥0.
We take a particular form of M(N) as
(5)M(N)=m01+m1N.Let us denoteK=rm0/f0 and α2=m1f0/m0.Then (1b) can be rewritten as
(6)dBdt=rB(1-BK)-α1NB-α2NB2.Now, we consider a population of densityN(t) 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:
(7)dNdt=sN(1-NL)+β1NB+β2NB2,
where s is the intrinsic growth rate of population, L is its carrying capacity in the absence of resource biomass, β1 and β2 are the growth rates of population in the presence of resource biomass.We assume that the resource biomass is harvested with the harvesting rateh(t)=qEB, where q is a positive constant and in fishery resource it is known as catchability coefficient, and E 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:(8a)dBdt=rB(1-BK)-α1NB-α2NB2-qEB,(8b)dNdt=sN(1-NL)+β1NB+β2NB2,B(0)>0,N(0)>0.In Section3, 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(9)Ω={(B,N)∈R2+:0≤B≤K,0≤N≤L0}
is a region of attraction for all solutions initiating in the interior of the positive quadrant, where
(10)L0=Ls(s+β1K+β2K2).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,
(11)P0(0,0),P1(B~,0),P2(0,L),P*(B*,N*).The equilibrium pointsP0 and P2 always exist. For the existence of equilibrium point P1, we note that B~ is given by(12a)B~=Kr(r-qE),This shows thatP1 exists if
(12b)r>qE.For the fourth equilibrium pointP*, we note that B* and N* are the positive solutions of the following algebraic equations:(13a)r(1-BK)-α1N-α2NB-qE=0,(13b)s(1-NL)+β1B+β2B2=0.Substituting the value ofN, from (13b) into (13a), we get a cubic equation in B as follows:
(14)a1B3+a2B2+a3B+a4=0,
where
(15)a1=KLα2β2>0,a2=KL(α2β1+α1β2)>0,a3=sr+KLα2s+α1β1KL>0,a4=α1KLs+(qE-r)Ks.The above equation (14) has a unique positive solution B=B*, if the following inequality holds:
(16)(α1L+qE)<r.After knowing the value ofB*, the value of N* can then be calculated from the relation
(17)N*=Ls(s+β1B*+β2B*).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 Theorems2 to 5.Theorem 2.
(i)  If the equilibrium pointP1(B~,0) exists, then P0 is always unstable in the B-N plane.
(ii)  If the equilibrium pointP1(B~,0) does not exist, then P0 is a saddle point with stable manifold locally in the B-direction and with unstable manifold locally in the N-direction.Theorem 3.
The equilibrium pointP1(B~,0), whenever it exists, is a saddle point with stable manifold locally in the B-direction and with unstable manifold locally in the N-direction.Theorem 4.
(i)  Ifr>(α1L+qE), then P2(0,L) is a saddle point with unstable manifold locally in the B-direction and stable manifold locally in the N-direction.
(ii)  Ifr<(α1L+qE), then P2(0,L) is always stable in the B-N plane.Theorem 5.
The interior equilibriumP*, whenever it exists, is always locally asymptotically stable in the B-N plane.In the next theorem, we are able to find a sufficient condition forP* to be globally asymptotically stable.Theorem 6.
Let the following inequality holds inΩ:(18){α1+β1+(α2+β2)K+β2B*}2<4sL(rK+α2N*).
Then the interior equilibriumP* is globally asymptotically stable with respect to all solutions initiating in the interior of the region Ω defined in Lemma 1.Proof of Theorem6.
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.
LetH(B,N)=1/BN>0 is a continuously differential function in the interior of the positive quadrant of B-N plane.
Let(19)F1(B,N)=rB(1-BK)-α1NB-α2NB2-qEB,F2(B,N)=sN(1-NL)+β1NB+β2NB2.
Then,
(20)Δ(B,N)=∂∂B(HF1)+∂∂N(HF2)=-1N(rK+α2N)-sBL<0.This shows thatΔ(B,N) does not change sign and is not identically zero in the positive quadrant of the B-Nplane. 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 B-N 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 timet is given by
(21)π(B,E,t)=(pqB-c)E,
where p is price per unit biomass and c is utilized cost per unit resource biomass. The bionomical equilibrium is P∞(B∞,N∞,E∞), where B∞,N∞, and E∞ are the positive solutions of(22a)B•=N•=π=0.Solving (22a), we get
(22b)B∞=cpq,(22c)N∞=Ls[s+β1cpq+β2c2p2q2],(22d)E∞=rq(1-cKpq)-(α1+cα2pq)N∞q.It is clear that(23)E∞>0,ifrq(1-cKpq)>(α1+cα2pq)N∞q.Thus the bionomical equilibriumP∞(B∞,N∞,E∞) exists under condition (23).Remark 8.
From (22b) and (23), it may be noted that B∞=(c/(pq))<K.


## 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]
(24)hMSY0=rK4.If the resource biomass is subjected to the harvesting by a population, the sustainable yield is given by(25)h=qEB*=rB*(1-B*K)-α1N*B*-α2N*B*2.We note that(26)∂h∂B*=0givesB*=K(r-α1N*)2(r+α2KN*),∂2h∂B*2<0.Thus,hMSY=B*((rB*/K)+α2N*B*), when B*=K(r-α1N*)/2(r+α2KN*).From the above equations, it is interesting to note that whenN*=0, then B*=K/2 and hMSY=rK/4=hMSY0.This result matches to Clark [1].Ifh>hMSY, then it denotes the overexploitation of the resource and if h<hMSY, 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 effortu=E is considered as an input. The harvesting effort is applied to the stock and produces a yield Y=qB per unit effort. We assume that the total yield Y per unit effort is subject to the constraint Ymin<Y<Ymax. Then our objective is to construct an output feedback control E(Y)=E0+u(Y(t)) in such a way that the steady state P*(B*,N*) 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)X•=f(X)+ug(X),
where u=E(t),
(28)f(X)=[rB(1-BK)-α1NB-α2NB2sN(1-NL)+β1NB+β2NB2],g(X)=[-qB0].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 effortE0, there exists a w>0 such that P* is globally asymptotically stable through the output feedback control law u=w(Y-Y*), where Y*=qB*.


## 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 valueJ of a continuous time stream of revenues is given by
(29)J=∫0∞e-δt(pqB(t)-c)E(t)dt,
where δ is the instantaneous rate of annual discount. Thus our objective is to
(30)maxJ,subject tothe state equations(8a)-(8b)and to thecontrol constraints0≤E≤Emax.For this purpose, we use Pontryagin’s Maximum Principle. The associated Hamiltonian function is given by(31)H=e-δt(pqB-c)E+λ1(t)[rB(1-BK)-α1NB-α2NB2-qEB{(1-BK)}]+λ2(t)[sN(1-NL)+β1NB+β2NB2{sN(1-NL)}],
where λ1 and λ2 are adjoint variables and σ(t)=e-δt(pqB-c)-λ1qB is called a switching function.The optimal controlE(t) which maximizes H must satisfy the following conditions:(32a)E=Emax,whenσ(t)>0,thatis,λ1eδt<p-cqB,(32b)E=0,whenσ(t)<0,thatis,λ1eδt>p-cqB.The usual shadow price isλ1eδt and the net economic revenue on a unit harvest is p-(c/qB). This shows that if the shadow price is less than the net economic revenue on a unit harvest, then E=Emax and if the shadow price is greater than the net economic revenue on a unit harvest, then E=0. When the shadow price equals the net economic revenue on a unit harvest, that is, σ(t)=0, then the Hamiltonian H becomes independent of the control variable E(t), that is, ∂H/∂E=0.This is a necessary condition for the singular controlE*(t) to be optimal over control set  0<E*<Emax.Thus, the optimal harvesting policy is(33)E(t)={Emax,σ(t)>00,σ(t)<0E*,σ(t)=0,
when σ(t)=0, we have
(34)λ1=e-δt(p-cqB).Now, in order to find the path of singular control, we utilize the Pontrygin’s Maximum Principle. According to this principle, the adjoint variablesλ1 and λ2 must satisfy
(35)dλ1dt=-∂H∂B,dλ2dt=-∂H∂N.The above equations can be rewritten as(36)dλ1dt=-[e-δtpqE+λ1{r(1-2BK)-α1N-2α2NB-qE}+λ2{β1N+2β2NB}{{r(1-2BK)-α1N-2α2NB-qE}}],dλ2dt=λ1(α1B+α2B2)-λ2{s(1-2NL)+β1B+β2B2}.Using (13a)-(13b), the above equations can be rewritten as
(37)dλ1dt=-[e-δtpqE-λ1{rBK+α2NB}+λ2{β1N+2β2NB}{{rBK+α2NB}}],(38)dλ2dt=e-δt(p-cqB)(α1B+α2B2)+λ2sNL.Equation (38) can again be written as
(39)dλ2dt-A1λ2=-B1e-δt,
where
(40)A1=sNL,B1=-(p-cqB)(α1B+α2B2).Solving (38), we get
(41)λ2=B1A1+δe-δt+K0eA1t,We note that whent→∞, then the shadow price λ2eδt is bounded if K0=0.Thus, we have(42)λ2=B1A1+δe-δt.Now from (37), we get
(43)dλ1dt-A2λ1=-B2e-δt,
where
(44)A2=-rBK-α2NB,B2=pqE-B1A1+δ(β1N+2β2NB).Equation (43) yields
(45)λ1=B2A2+δe-δt+K1eA2t.Again for the shadow priceλ1eδt to be bounded, we must have K1=0.Thus,(46)λ1=B2A2+δe-δt.Equations (34) and (46) yield
(47)B2A2+δ=(p-cqB).Substituting the values ofA1,B1,A2, and B2into (47), we get
(48)E=1pq(p-cqB){{(α1B+α2B2)(β1N+2β2NB)((sN/L)+δ)}(δ+rBK+α2NB)+(α1B+α2B2)(β1N+2β2NB)((sN/L)+δ)}.Hence, solving (13a)-(13b) with the help of (48), we get an optimal solution (Bδ,Nδ) and the optimal harvesting effort E=Eδ.


## 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):
(49)r=1.6,s=1.2,K=100,L=100,E=100,p=0.1,q=0.01,α1=0.001α2=0.01,β1=0.01,β2=0.1,
with initial condition: B(0)=5 and N(0)=8.For this set of parameters, condition (16) for the existence of the interior equilibrium is satisfied. This shows that P*(B*,N*)  exists and it is given by
(50)B*=0.4790,N*=102.3108.Thus the interior equilibrium pointP*(B*,N*) 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 P* to be globally asymptotically stable, no conclusion can be drawn at this stage. The behavior of B and N with respect to time t 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 ofB and N versus t for the values of parameters given in (49).Now, we choose another set of values of parameters as follows:(51)r=1.6,s=1.2,K=100,L=100,E=100,p=0.1,q=0.01,α1=0.001α2=0.0001,β1=0.01,β2=0.0001
with different initial values.For the set of values of parameters given in (51), it may be noted that the positive equilibrium P*(B*,N*) exists and it is given by
(52)B*=17.4311,N*=117.0579.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 P*(B*,N*) is locally as well as globally asymptotically stable in the interior of the first quadrant. In Figure 2, we have plotted the behavior of B and N with different initial values. Figure 2 shows that all the trajectories starting from different initial points converge to the point P*(17.4311,117.0579). This shows that P*(17.4311,117.0579) is globally asymptotically stable.Figure 2
Global stability ofP*(B*,N*).It may be noted here thatα2 and β2 are important parameters governing the dynamics of the system. Therefore, we have plotted the behavior of B and N with respect to time for different values of α2 in Figure 3 and for different values of β2 in Figure 4.Plot ofB and N with respect to time t for different values of α2, and other values of parameters are the same as given in (51).
(a)
(b)Plot ofB and N 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.
(a)
(b)From Figure3(a), we note that B decreases as α2 increases. If α2 is very small, then B initially increases slightly and then decreases and converges to its equilibrium level. If α2 increases beyond a threshold value, then B always decreases and settles down at its equilibrium level. Figure 3(b) shows that N also decreases as α2 decreases. This is due to the fact that with the increase in α2, the equilibrium level of B decreases and since the population N is dependent on the resource biomass B, thus N also decreases. However, it is interesting to note here that N always increases with respect to time t and finally attains its equilibrium level.From Figure4(a), it is noted that B decreases as β2 increases. Figure 4(b) shows that N increases as β2 increases. It may also be noted that for all positive values of β2, N increases continuously with time and finally settles down at its steady state. However, the resource biomass B 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 c=0.001 and δ=5. Then we get the optimal equilibrium levels as given below:
(53)Bδ=1.40895,Nδ=101.191,Eδ=146.201.We observe thatE is a very important parameter which governs the dynamics of the system. The behavior of B with respect to time is shown in Figure 5 and the behavior of N with respect to time is shown in Figure 6 for different values of E. Figure 5 shows that if E is less than its optimal level (Eδ=146.201), then B increases and settles down at its equilibrium level. If E>Eδ, then B decreases and settles down at a lower equilibrium level. If E is large enough, then B tends to zero. Figure 6 shows that if E<Eδ, and then N increases with time and finally obtains its equilibrium level. If E>Eδ, then N again increases with time, but it settles down at a lower equilibrium level. This is due to the fact that when E>Eδ, then B decreases and hence N also decreases. This suggests that the harvesting effort E should be always kept less than Eδ (optimal harvesting level) so that both the resource and the population can be maintained at an optimal level.Figure 5
Plot ofB versus t for different values of E.Figure 6
Plot ofN versus t for different values of E.


## 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 equilibriumP* (whenever it exists) is always locally asymptotically stable. But P* is not always globally asymptotically stable. However, we have found a sufficient condition under which P* is globally asymptotically stable. This condition gives a threshold value of the specific growth rate (r) of resource biomass. If r is larger than this threshold value, then P* 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 (h) and maximum sustainable yield (hMSY) has been carried out. It has been shown that if h>hMSY, then the resource biomass will tend to zero and if h<hMSY, 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α2,β2 and E are important parameters governing the dynamics of the system. It has been found that if α2 increases, then the resource biomass and the population both decrease. However, if β2 increases, then the population density increases and the resource biomass density decreases. A threshold value for the optimal harvesting effort Eδ has been found theoretically as well as numerically. It has been shown that the harvesting effort E should be always kept less than Eδ to maintain the resource and the population at an optimal equilibrium level.



---

*Source: 613706-2013-01-17.xml*