# Optimal Harvesting Policies for a Stochastic Food-Chain System with Markovian Switching

Authors: Yanming Ge; Yifan Xu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

An optimal harvesting
problem for a stochastic food-chain system with Markovian switching is investigated in this paper. Firstly,
the existence, uniqueness, and positivity of the food-chain
system’s solution are proved. Secondly, persistent in mean of the
system is obtained. Then the optimal harvesting policy is
discussed. Finally, the main results are
illustrated by several examples.


---

## Body

## 1. Introduction

Optimal harvesting problem is an important and interesting topic from both biological and mathematical point of view. Since Clark’s works [1, 2], one of the most important area-optimal harvesting problems have received a lot of attention and been studied widely. Among these studies, a large number of literatures were focused on deterministic models [3–9], with some on the stochastic versions [10–17], but only a few on the food-chain systems. Furthermore, it is well known that the theory of food chains illustrated the balance of nature and that no animal or plant can exist independently. Motivated by these arguments presented above, we are interested in the optimal harvesting problems on the following stochastic food-chain system:(1)dx=r1ξt-e-a1x-a2ym1ξt+xxdt+σ1ξtxdB1t,dy=r2ξt-b1ym1ξt+x-b2zm2ξt+yydt+σ2ξtydB2t,dz=r3ξt-czm2ξt+yzdt+σ3(ξ(t))zdB3(t)with the initial value (x0,y0,z0). Where Bi(t), i=1,2,3, is a standard Brownian motion and e is the harvesting effort (control parameter), and x(t), y(t), z(t) represent the population densities of three species (resource, consumer, and predator) at time t, respectively. All parameters are positive constants and parametric functions are continuous and positive. ri(·), i=1,2,3, represent the intrinsic growth rate of species x, y, z, respectively; a1 measures the strength of competition among individuals of species x; a2 is the maximum value of the per capita reduction rate of x due to y; b1, b2, and c have similar meaning to a2; m1(·) measures the extent to which the environment provides protection to species x and y; m2(·) measures the extent to which the environment provides protection to species y and z; ξ(t) be a right continuous Markov chain; σi(·), i=1,2,3, represents the intensity of the white noise. This system is the extension of a predator-prey model with modified Leslie-Gower and Holling-type II schemes with stochastic perturbation which was discussed by Ji et al. [18], Song et al. [19], and Guo et al. [20], and there the factor of Markovian switching is not considered.In the most literatures [1, 2] the sustainable yield function is used as the harvesting function. Here the harvesting function associated with (1) is (2)F(e)=elimt→∞∫0t‍x(s)dst.This type of harvesting function is also used by some other papers, such as Wang [21, chapter 4] and Zou and Wang [22] and defined as the time averaging yield function. The optimal harvesting problem considered in this paper is then stated as follows. Find a harvesting effort e∗ such that (3)F(e∗)=supe≥0{F(e)}.Based on the aforementioned discussion, obviously, the first and most important duty is to discuss the existence oflimt→∞(∫0t‍xsds/t) and then the optimal harvesting problem. Therefore, the rest of the paper is organized as follows. In Section 2, we show that system has a global positive solution. In Section 3, we obtain some long time behavior of the solution, especially the property of persistent in mean, which ensures the existence of the time averaging yield function and its explicit expression is given. In Section 4, the optimal harvesting policies are investigated. In Section 5, we illustrate our main results through several numerical examples. Last but not least, conclusions are drawn in Section 6.On the other hand, for convenience, we give some notations and assumptions in the rest of this section.Throughout this paper, unless otherwise specified, let(Ω,F,{Ft}t≥0,P) be a complete probability space with a filtration {Ft}t≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while F0 contains all P-null sets). The standard Brownian motion Bi(t), i=1,2,3, is defined on this probability space.The right continuous Markov chainξ(t) on this probability space taking values in a finite-state space S={1,2,…,N} with the generator γ=(γkl)N×N is given by (4)Pξt+δ=l∣ξt=k=γklδ+o(δ)ifl≠k,1+γklδ+o(δ)ifl=k,where δ>0. Here γkl is the transition rate from k to l and γkl≥0 if l≠k, while (5)γkk=1-∑l≠k‍γkl.We assume that the Markov chain ξ(·) and the Brownian motion Bi(·) are independent of each other, i=1,2,3. As a standing hypothesis we also assume in this paper that the Markov chain ξ(t) is irreducible. This is very reasonable as it means that the system will switch from any regime to any other regime. This is equivalent to the condition that, for any k,l∈S, one can find finite numbers i1,…,in∈S such that γki1γi1i2⋯γinl>0. Under this condition, the Markov chain has a unique stationary distribution π=(π1,π2,…,πN)∈RN and πk>0 for any k∈S.In order to obtain some properties of the system, some assumptions are given in the following. These assumptions are conventional; they guarantee that the ecosystem is not collapsed as time lapses.Assumption 1.
Consider0<D1∶=minξ∈S{r1(ξ)-e-0.5σ12(ξ)}, 0<Di∶=minξ∈S{ri(ξ)-0.5σi2(ξ)}≤maxξ∈S{ri(ξ)-0.5σi2(ξ)}∶=Ai, i=2,3.Assumption 2.
D
2
-
b
2
(
P
3
+
ϵ
)
>
0
,
D
1
-
a
2
(
P
2
+
ϵ
)
>
0
,
P
1
>
0, where P3∶=(1/c)∑i=1N‍πi(r3(i)-0.5σ32(i)), P2∶=(1/b1)[∑i=1N‍πi(r2(i)-0.5σ22(i))-b2P3],  P1∶=(1/a1)[∑i=1N‍πi(r1(i)-e-0.5σ12(i))-a2P2], and ϵ is positive and sufficiently small.Letmˇi=maxξ∈S{m(ξ)}, m^i=minξ∈S{m(ξ)}, i=1,2, and K denotes a float constant in the rest of this paper, which expresses different constants in different positions.The key method used in this paper is the comparison theorem for stochastic equations. This theorem for stochastic differential equations was developed by Ikeda and Watanabe [23] and has been used by many authors [24–26].


## 2. Positive and Global Solutions

As the state of the system(x(t),y(t),z(t)) is the population density of species in the system at time t, it should be nonnegative. Moreover, in order for a stochastic differential equation to have a unique global (i.e., no explosion in a finite time) solution for any given initial data, the coefficients of the equation are generally required to satisfy the linear growth condition and local Lipschitz condition [25]. However, the coefficients of each equation in system obey neither the linear growth condition nor local Lipschitz continuous. In this section, we show existence and uniqueness of the positive solution.Lemma 3.
For any initial valuex0>0, y0>0, z0>0, system has a unique positive local solution (x(t),y(t),z(t)) for t∈[0,τe) almost surely (a.s.), where τe is the explosion time.Proof.
To begin with, consider the following equations:(6)dut=b1ξt-a1eut-a2evtm1ξt+eutdt+σ1ξtdB1t,dvt=b2ξt-b1evtm1ξt+eut-b2ewtm2ξt+evtdt+σ2ξtdB2t,dwt=b3ξt-cewtm2ξt+evtdt+σ3(ξ(t))dB3(t)on t≥0 with initial value u(0)=lnx0, v(0)=lny0, w(0)=lnz0, where b1(t)=r1(ξ(t))-e-0.5σ12(ξ(t)), bi(t)=ri(ξ(t))-0.5σi2(ξ(t)), i=2,3. Notice that the last equations’ coefficients satisfy the local Lipschitz condition; thus there is a unique solution (u(t),v(t),w(t)) on t∈[0,τe). Therefor, it follows from Itô’s formula that x(t)=eu(t), y(t)=ev(t), z(t)=ew(t) is the unique positive local solution of system with initial value x0>0, y0>0, z0>0.Lemma3 only tells us that there is a unique positive local solution to (1). Next, we show this solution is global, that is, τe=+∞. For convenience, we define six equations:(7)dxˇt=r1ξt-e-a1xˇxˇdt+σ1ξtxˇdB1t,xˇ0=x0,dyˇt=r2ξt-b1yˇm1ξt+xˇyˇdt+σ2ξtyˇdB2t,yˇ0=y0,dz^t=r3ξt-cz^m2ξtz^dt+σ3ξtz^dB3t,z^0=z0,dzˇt=r3ξt-czˇm2ξt+yˇzˇdt+σ3ξtzˇdB3t,zˇ0=z0,dy^t=r2ξt-b1y^m^1-b2zˇm^2y^dt+σ2ξty^dB2t,y^0=y0,dx^t=r1ξt-e-a1x^-a2yˇm^1x^dt+σ1ξtx^dB1t,x^0=x0.Obviously, whent∈[0,τe), by the comparison theorem for stochastic equations [27, Theorem  3.1], it yields (8)x^≤x≤xˇ,y^≤y≤yˇ,z^≤z≤zˇ.Furthermore, x^, xˇ, y^, yˇ, z^, zˇ are all existing on t≥0, and hence we have the following.Theorem 4.
There is a unique positive solutionx(t), y(t), z(t) on t≥0 to (1) a.s. for any initial value x0>0, y0>0, z0>0, and relations(9)x^≤x≤xˇ,y^≤y≤yˇ,z^≤z≤zˇare all satisfied on t≥0.


## 3. The Long Time Behavior

Theorem4 shows that the solution of the system (1) will remain in the positive cone R+3. This nice property provides us with a great opportunity to discuss how the solution varies in R+3 in detail. In this section we will give some long time behavior of the solution, especially the property of persistent in mean, which ensures the existence of limt→∞(∫0t‍xsds/t).Lemma 5 (see [28]).
If Assumption1 is satisfied, then one has (10)limt→∞lnxˇ(t)t=0,limt→∞lnz^(t)t=0a.s.Lemma 6.
If Assumption1 is satisfied, then one has (11)limt→∞lnyˇ(t)t=0,limt→∞lnzˇ(t)t=0a.s.Proof.
Firstly, we give an auxiliary equation(12)dyˇ¯t=r2ξt-b1yˇ¯mˇyˇ¯dt+σ2ξtyˇ¯dB2t,yˇ¯0=y0.Obviously, yˇ¯(t)≤yˇ(t), and using the similar method of Lemma 5, we have (13)limt→∞lnyˇ¯(t)t=0a.s.
Therefore,(14)0=liminft→∞lnyˇ¯(t)t≤liminft→∞lnyˇ(t)ta.s.
Next, we need only to provelimsupt→∞(lnyˇt/t)≤0a.s.
The quadratic variation of∫0t‍σ2(ξ(s))dB2(s) is ∫0t‍σ22(ξ(s))ds≤Kt, and by the strong law of large numbers for local martingales, we have (15)∫0t‍σ2(ξ(s))dB2(s)t⟶0a.s.,t⟶∞.
Therefore, for allϵ>0, ∃0<T<∞, we have (16)∫0t‍σ2ξsdB2s<ϵta.s.,t≥T.
From this, we have(17)∫st‍σ2ξsdB2s<ϵ(s+t)a.s.,t>s≥T.
On the other hand, from Lemma5, we have (18)-ϵt≤lnxˇt≤ϵta.s.,t≥T.
By the arguments as above, whent>s≥T, we can get(19)1yˇt=e-[∫Tt‍(r2(ξ(s))-0.5σ22(ξ(s)))ds+∫Tt‍σ2(ξ(s))dB2(s)]×1yˇT+∫Tt‍b1m1ξs+xˇ×e∫Ts‍(r2(ξ(τ))-0.5σ22(ξ(τ)))dτ+∫Ts‍σ2(ξ(τ))dB2(τ)ds1y˘T+∫Tt‍b1m1ξs+x˘≥∫Tt‍b1m1ξs+xˇe∫ts‍(r2(ξ(τ))-0.5σ22(ξ(τ)))dτ+∫ts‍σ2(ξ(τ))dB2(τ)ds≥∫Tt‍b1m1ξs+eϵse∫ts‍(r2(ξ(τ))-0.5σ22(ξ(τ)))dτ+∫ts‍σ2(ξ(τ))dB2(τ)ds≥∫Tt‍b1m1ξs+1e-ϵseA2s-te-ϵt+sds≥b1mˇ1+1e-A2+ϵt∫Tt‍eA2-2ϵsds=b1mˇ1+1A2-2ϵe-A2+ϵteA2-2ϵt-eA2-2ϵT.
Therefore, we obtain(20)A2+ϵ≥limsupt→∞lnyˇtt+A2-2ϵ;that is (21)limsupt→∞lnyˇ(t)t≤3ϵ.
For the arbitrary ofϵ>0, we must have (22)limsupt→∞lnyˇ(t)t≤0.
Hence,limt→∞(lnyˇt/t)=0a.s., and the second conclusion can be proved similarly.Theorem 7.
If Assumption1 is satisfied, then we have (23)limt→∞lnz(t)t=0a.s.,(24)limt→∞∫0t‍zs/m2ξs+ysdst=1c∑i=1N‍πir3i-0.5σ32i.Proof.
Following fromlimt→∞(lnz^t/t)=0a.s., limt→∞(lnzˇt/t)=0a.s., and z^≤z≤zˇ. Obviously, we have (25)limt→∞lnz(t)t=0a.s.
In the following, we prove the second conclusion.
LetV(z,ξ)=lnz; applying Itôs formula gives (26)dlnzt=r3ξt-0.5σ32ξt-cztm2ξt+ytdt+σ3ξtdB3t.
Hence(27)lnzt-lnz0=∫0t‍r3ξs-0.5σ32ξsds-∫0t‍czsm2ξs+ysds+∫0t‍σ3ξsdB3s.
Based on the first conclusion in this theorem, the strong law of large numbers for local martingales, and the ergodic property of Markov chain, the second conclusion is proved.Theorem 8.
If Assumptions1 and 2 are satisfied, then one has (28)limt→∞lny(t)t=0,limt→∞lnx(t)t=0,a.s.limt→∞∫0t‍(ys/m1ξs+xs)dst=1b1∑i=1N‍πir2i-0.5σ22i-b2P3,limt→∞∫0t‍xsdst=1a1∑i=1N‍πir1i-e-0.5σ12i-a2P2.Proof.
Based onlimt→∞(lnyˇt/t)=0 and y≤yˇ, we have (29)limsupt→∞lny(t)t≤limsupt→∞lnyˇ(t)t=0a.s.
Next, we need only to proveliminft→∞(lny(t)/t)≥0a.s.
The quadratic variation of∫0t‍σ2(ξ(s))dB2(s) is ∫0t‍σ22(ξ(s))ds≤Kt, and by the strong law of large numbers for local martingales, we have (30)∫0t‍σ2(ξ(s))dB2(s)t⟶0a.s.,t⟶∞.
Therefore, for allϵ>0, ∃0<T<∞, we have (31)∫0t‍σ2ξsdB2s<ϵta.s.,t≥T.
From this, we have(32)∫st‍σ2ξsdB2s<ϵ(s+t)a.s.,t>s≥T.
On the other hand, from Lemma5, we have (33)-ϵt≤lnxˇt≤ϵta.s.,t≥T.
By the arguments as above, whent>s≥T, we can get (34)1yt=1yTe-[∫Tt‍(r2(ξ(s))-0.5σ22(ξ(s)))ds-∫Tt‍(b2zs/(m2ξs+ys))ds+∫Tt‍σ2(ξ(s))dB2(s)]+∫Tt‍b1m1ξs+xse∫ts‍(r2(ξ(τ))-0.5σ22(ξ(τ)))dτ-∫ts‍(b2zτ/(m2ξτ+yτ))dτ+∫ts‍σ2(ξ(τ))dB2(τ)ds≤1yTe-[∫Tt‍(r2(ξ(s))-0.5σ22(ξ(s)))ds-∫Tt‍(b2zs/(m2ξs+ys))ds+∫Tt‍σ2(ξ(s))dB2(s)]+∫Tt‍b1m^1e∫ts‍(r2(ξ(τ))-0.5σ22(ξ(τ)))dτ-∫ts‍(b2zτ/(m2ξτ+yτ))dτ+∫ts‍σ2(ξ(τ))dB2(τ)ds=:I1+I2.
Based on the second conclusion of Theorem7, for all ϵ>0, ∃0<T<∞, when t>s>T, we have (35)∫Tt‍z(s)m2(ξ(s))+y(s)ds<(P3+ϵ)t-P3-ϵTa.s.,∫st‍z(τ)m2(ξ(τ))+y(τ)dτ<(P3+ϵ)t-(P3-ϵ)sa.s.
Thus,(36)I1≤1yTe-D2-b2P3+ϵt-T+t+2b2+1Tϵ≤Ket+2b2+1Tϵa.s.,I2≤∫Tt‍Ke[t+(2b2+1)s]ϵds≤Ke(2b2+2)tϵa.s.
Hence, we obtain1/y(t)≤2Ke(2b2+2)tϵa.s., and furthermore (37)limsupt→∞(-lny(t)t)≤(2b2+2)ϵ.In other words, (38)liminft→∞lny(t)t≥-(2b2+2)ϵ.
For the arbitrary ofϵ>0, we must have (39)liminft→∞lny(t)t≥0.
The first assertion is proved. The assertionlimt→∞(lnxt/t)=0,a.s. can be proved similarly.
Similarly to the proof of the second assertion of Theorem7, the last two assertions can be proved.Definition 9 (see [5]).
The system is said to be persistent in mean, if(40)liminft→∞∫0t‍x(s)dst>0,liminft→∞∫0t‍y(s)dst>0,liminft→∞∫0t‍z(s)dst>0a.s.Theorem 10.
If Assumptions1 and 2 are satisfied, then the system is persistent in mean.Proof.
l
i
m
i
n
f
t
→
∞
(
∫
0
t
‍
x
s
d
s
/
t
)
>
0 is already proved in Theorem 8.
From the third assertion of Theorem8, we have (41)liminft→∞∫0t‍ysdsm^1t≥limt→∞∫0t‍ys/m1ξs+xsdst=P2>0;
that is,liminft→∞(∫0t‍ysds/t)>0.
From the second assertion of Theorem7, we have (42)liminft→∞∫0t‍zsdsm^2t≥limt→∞∫0t‍zs/m2ξs+ysdst=P3>0;
that is,liminft→∞(∫0t‍zsds/t)>0.
This theorem is proved.


## 4. The Optimal Harvesting Policies

Based on the explicit expression of the time averaging yield function obtained in the last section, here we discuss the optimal harvesting problem mentioned in Section1.Theorem 11.
If Assumptions1 and 2 is satisfied, then the optimal harvesting effort is (43)e∗=min{D1¯-a2(P2+ϵ),e¯},where D1¯=minξ∈S{r1(ξ)-0.5σ12(ξ)}, e¯=0.5Σ-0.5a2P2, Σ=∑i=1N‍πi(r1(i)-0.5σ12(i)), and the optimal harvesting output is (44)F(e∗)=1a1Σ-a2P2e∗-e∗2.Proof.
Based on Theorem8, the optimization problem can be expressed as follows:(45)maxFe=1a1Σ-a2P2e-e2s.t.e≤D1¯Assumption1s.t.e≤D1¯-a2P2+ϵAssumption2s.t.e≤Σ-a2P2Assumption2s.t.e≥0.
From the definitions ofD1¯ and Σ, we get D1¯≤Σ. Therefore, the above optimization problem can be simplified as follows:(46)maxFe=1a1Σ-a2P2e-e2s.t.0≤e≤D1¯-a2P2+ϵ.
Because the objective function is concave, and we can obtain the unique maximum point easily as(47)e∗=min{D1¯-a2(P2+ϵ),e¯},here the e¯ is obtained by letting dF(e)/d(e)=0.
Substituting it into the harvesting function, we obtain the optimal harvesting output(48)F(e∗)=1a1[(Σ-a2P2)e∗-e∗2].
This theorem is proved.Remark 12.
(i) That the feasible zone of optimization problem (46) is nonempty is guaranteed by Assumptions 1 and 2.
(ii) From the explicit expression of the optimal harvesting effort, we can easily investigate how the parameters influence on it, such thate∗ is decreasing in a2, and this claim coincides with the fact that if the consumer’s (y(t)) consuming capacity is enhanced (a2 augments), the harvesting effort must reduce (e∗ go down), or the resource (x(t)) will be extinct and the whole ecosystem is crashed.


## 5. Numerical Results

We present numerical experiments in this section to show how the proposed model works in the constructive examples. The results enhance the readers to understand the theoretical conclusions from the practical applications.Here, we use the Milstein method [29] to construct the discretization equation of (1); that is,(49)xk+1=xk+xkr1ξkΔt-e-a1xk-a2ykm1ξkΔt+xkΔt+σ1ξkΔtxkΔtζk+0.5σ12ξkΔtxk2ζk2-1Δt,yk+1=yk+ykr2ξkΔt-b1ykm1ξkΔt+xk-b2zkm2ξkΔt+ykΔt+σ2(ξ(kΔt))ykΔtηk+0.5σ22(ξ(kΔt))yk2(ηk2-1)Δt,zk+1=zk+zkr3ξkΔt-czkm2ξkΔt+ykΔt+σ3(ξ(kΔt))zkΔtςk+0.5σ32(ξ(kΔt))zk2(ςk2-1)Δt,where ζk, ηk, and ςk, k=1,…,n, are the Gaussian random variables.For simplicity, assume that the random environments are modeled by a two-state Markov chain with state setS={1,2} and generator (50)Q=-7755.The other parameters are defined as follows:r1(1)=8, r1(2)=9, r2(1)=6, r2(2)=6, r3(1)=4, r3(2)=3, σ1(1)=0.8, σ1(2)=0.9, σ2(1)=0.6, σ2(2)=0.6, σ3(1)=0.4, σ3(2)=0.3, m1(1)=5, m1(2)=4, m2(1)=3, m2(2)=2, a1=a2=b2=1, b1=c=10, e=0. In this scenario, B1=7.68, B2=5.82, B3=2.955, P1=7.66, P2=0.55, P3=0.34,  B1-a2(P2+ɛ)=7.13-ɛ, and B2-b2(P3+ɛ)=5.48-ɛ. Obviously, Assumptions 1 and 2 are satisfied.Based on the aforementioned discussion, we obtain the following results.Figures1 and 2 show that the solutions of (1) are positive in the deterministic environment (without regime switching); that is, ξ(t)≡1 or 2. Figure 4 shows that the solutions of (1) are positive in the random environment (with regime swithcing); the random environment is described by Figure 3. They are all identical to Theorem 4. Figure 5 shows limt→∞(lnxt/t)=0, limt→∞(lnyt/t)=0, and limt→∞(lnzt/t)=0 in the random environment described by Figure 3, and this is consistent with Theorems 7 and 8.Figure 1
Sample path of (1) with initial condition (x0,y0,z0)=(1,1,1) and ξ(t)≡1.Figure 2
Sample path of (1) with initial condition (x0,y0,z0)=(1,1,1) and ξ(t)≡2.Figure 3
Sample path ofξ(t) with initial condition ξ(0)≡1.Figure 4
Sample path of (1) with initial condition (x0,y0,z0)=(1,1,1) and sample path of ξ(t) is Figure 3.Figure 5
The curves oflnx/t, lny/t, and lnz/t.


## 6. Conclusions

This paper studies an optimal harvesting problem for a food-chain system with markovian switching. Based on the properties, the food-chain system’s solution is existing, unique, and positive; the system is persistent in mean, and the rationality of the optimal harvesting problem is proved. Then the optimal harvesting policy is obtained.Nevertheless, there are rooms to continue work on this issue, such that more than one of control variables in the system are considered. The permanence and extinction of the system and the stability in distribution need to be investigated too.



---

*Source: 875159-2015-06-24.xml*