# Optimal Harvest Control in a Singular Prey-Predator Fishery Model with Maturation Delay and Gestation Delay

Authors: Chao Liu; Longfei Yu; Qingling Zhang
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## Abstract

This paper presents a singular prey-predator fishery model, where maturation delay for prey and gestation delay for predator are considered. Fishing efforts are introduced to harvest prey and predator population, which are developed as control instruments to investigate optimal utilization of fishery resource. By analyzing associated characteristic equation, local stability analysis is studied due to combined variations of double time delays. Furthermore, Pontryagin’s maximum principle is utilized to characterize optimal harvest control, and the optimality system is numerically solved based on an iterative method.


---

## Body

## 1. Introduction

According to statistics from the Food and Agriculture Organization of United Nations [1], approximate fifty-three percent of fish stock under observation has experienced overexploitation or depletion, which reiterate the fact that fishery needs to be managed with an effective and carefully defined objective to prevent overexploitation and replenish depleted stock [1]. In order to ensure sustainable fishery, it requires precision in stock assessment and reliability in fishery modelling [2, 3]. In recent decades, there are growing research interests in formulating bioeconomic mathematical models to obtain more comprehensive indications of feedback effect between exploitation activity and fishery resources [4–13].In [5], authors establish a dynamic model where both prey and predator fishery resource population are exploitable, which is as follows:(1)x˙t=rxt1-xtK-βxtytα+xt-Etxt,y˙t=ytmβxtα+xt-d-Etyt,where x(t) and y(t) represent population density of prey and predator population with time t, respectively. r is the intrinsic growth rate of prey, K stands for maximum environmental capacity, β denotes predation rate of predator, m represents biomass transmission rate from predation activity, and d is death rate of predator.It is well known that commercial harvesting fluctuates with dynamical variation of economic interest [14, 15]. Although the harvested biological models have received great attention from both mathematical and theoretical biologists, little work has been done on dynamic effect of economic interest on harvested prey-predator fishery model [14, 15]. In 1954, Gordon proposes common-property resource economic theory [14], which investigates dynamic effect of commercial harvesting on ecosystem from an economic perspective [14, 15]. In [14], an algebraic equation is proposed to investigate economic interest of commercial harvesting:(2)NetEconomicRevenue=TotalRevenueTR-TotalCostTC.In order to show dynamical variation of unit price of harvested population due to change of amount of harvested population,a1/(b1+E(t)x(t)) is assumed to be unit price of commercially harvested prey, and a2/(b2+Etyt) is assumed to be unit price of commercially harvested predator. Hence, TR in (2) takes the following form TR=a1Etxt/(b1+Etxt)+a2Etyt/(b2+Etyt), where a1,a2,b1,b2 are all positive constants. Let c represent cost of commercial harvesting. Hence, TC in (2) takes the following form TC=cE(t).Recently, many theoreticians and experimentalists have investigated complex dynamics of prey-predator fishery system [16, 17]; it reveals that time delay may cause the loss of stability and other complicated dynamical behavior such as the periodic structure and bifurcation phenomenon [15, 18, 19].Keeping these aspects in view, we will extend work in [5] by incorporating maturation delay for prey and gestation delay for predator into system (1). In this paper, prey is considered to be delayed by maturation delay τ1≥0 due to population crowding [15, 16], and predator is assumed to be delayed by gestation delay τ2≥0 [16]. Based on (1) and (2), a singular prey-predator fishery model is established as follows:(3)x˙t=rxt1-xt-τ1K-βxtytα+xt-Etxt,y˙t=ytmβxt-τ2α+xt-τ2-d-Etyt,0=Etxta1b1+Etxt+yta2b2+Etyt-c-v,where c represents cost of commercial harvesting, v is assumed to be economic interest of commercial harvesting, E(t) represents commercial harvesting amount with respect to time t, and interpretations for other parameters and state variables share the same interpretations introduced in system (1) and (2). Furthermore, system (3) is investigated with initial conditions(4)xθ≥0,y0≥0,E0≥0,θ∈-τ,0,τ=max⁡τ1,τ2.System (3) is rewritten as follows:(5)Ξx˙ty˙t0=rxt1-xt-τ1K-βxtytα+xt-Etxtytmβxt-τ2α+xt-τ2-d-EtytEtxta1b1+Etxt+yta2b2+Etyt-c-v,where the third equation in system (3) does not contain any differentiated variables; hence the third row in leading matrix Ξ=100010000 has a zero row.The rest of the sections of this paper are organized as follows. Positivity of any solution of system (3) and uniform persistence of system (3) are investigated in the second section. Local stability analysis of system (3) around interior equilibrium is discussed in the third section. The optimal system is derived and solved numerically based on an iterative method with Runge-Kutta fourth order scheme. Optimal harvest control problems associated with maximizing total discounted net revenues from the fishery and minimizing harvest cost are discussed in the fourth section. Numerical simulations are carried out to support theoretical analysis in the fifth section. Finally, this paper ends with a conclusion.


## 2. Positivity and Uniform Persistence

Theorem 1.
Any solutions of system (3) with initial conditions (4) are positive.Proof.
Due to lemma in [20] and Theorem A.4 in [21], any solution of system (3) with initial conditions (4) exists uniquely and each component of solution remains within interval [0,X0) for some X0>0. Standard and simple arguments show that any solution of system (3) always exists and stays positive.Theorem 2.
Ifv>0, τ1 and τ2 are bounded, and the following inequalities (6) hold, and then system (3) with initial conditions (4) is uniformly persistent:(6)mβKerτ1>αd,a1b2βKerτ1+a2b1rα+Kerτ1>cb1b2β,cαb1>a1Kerτ1mβKerτ1-αd,cb2αβ>a2rα+Kerτ1mβKerτ1-αd.Proof.
By using Taylor series expansion [2], for x(t) and τ1>0,τ2>0 we have(7)xt-τ1=xt-τ1ddtxt-τ1x˙t+⋯,xt-τ2=xt-τ2ddtxt-τ2x˙t+⋯.
Consequently, it can be derived that(8)xt-τ1≤xt,xt-τ2≤xt.
By using (8) and the first equation of system (3), it gives that (9)xt≤xt-τ1erτ1,xt≤rxt1-xtKerτ1,which derives that(10)limsupt→∞⁡xt≤Kerτ1≔Q1.
Ifτ1 is bounded, then it is easy to show that Q1>0 and is bounded. According to (10), it can be derived that there exists T1>0 such that x(t)≤Q1 holds for t>T1+τ1.
Based on the first equation of system (3) and practical interpretations of prey survival, it can be derived that there exists T2>T1+τ1 such that y(t)≤r[α+x(t)]/β, and(11)limsupt→∞⁡yt≤rα+xtβ≤rα+Kerτ1β≔Q2.
Similarly, based on the second equation of system (3) and practical interpretations of predator survival, it is shown that there exists T3>T2 such that(12)Et≤mβxt-τ2α-d.
By using (8), if τ2 is bounded then it follows that E(t)≤(mβxt/α)-d, and(13)limsupt→∞⁡Et≤mβKerτ1-αdα≔Q3.
IfmβKerτ1>αd, then it can be derived that Q3>0 and is bounded. According to (13), it is shown that there exists T4>T3 such that E(t)≤Q3 holds for t>T4.
By virtue of the third equation of system (3), when v>0, it can be derived that there exists T5>T4 such that (14)vEt=xta1b1+Etxt+yta2b2+Etyt-c≤a1xtb1+a2ytb2-c,which implies that(15)liminft→∞⁡Et≥vb1b2βa1b2βKerτ1+a2b1rα+Kerτ1-cb1b2β≔W3.
Ifa1b2βKerτ1+a2b1r(α+Kerτ1)>cb1b2β, then it is easy to show that W3>0 and is bounded, and then there exists T6>T5 such that E(t)≥W3 holds for t>T6.
By virtue of the third equation of system (3), when v>0, it can be derived that there exists T7>T6 such that (16)Etxt>b1a1c-a2b2Etyt,Etyt>b2a2c-a1b1Etxt.
According to (10), (11), and (13), it can be obtained that (17)liminft→∞⁡xt≥cb2αβ-a2rα+Kerτ1mβKerτ1-αdαβ≔W1,liminft→∞⁡yt≥b2αcb1-a1Kerτ1mβKerτ1-αda2b1mβKerτ1-αd≔W2.
If the following two inequalities hold(18)cαb1>a1Kerτ1mβKerτ1-αd,cb2αβ>a2rα+Kerτ1mβKerτ1-αd,then W1>0 and W2>0 and is bounded, and there exists T8>T7 such that x(t)≥W1 and y(t)≥W2 hold for t>T8.
Based on the above analysis, it can be concluded that system (3) with initial conditions (4) is uniformly persistent.


## 3. Local Stability Analysis

Firstly, interior equilibriumP∗ of system (3) in the case of v>0 is computed as follows: P∗(x∗,y∗,E∗)=(x∗,α+x∗Kd+rK-x∗/Kβ-mx∗,mβx∗/(α+x∗)-d), and x∗ satisfy the following equation:(19)A0x∗6+A1x∗5+A2x∗4+A3x∗3+A4x∗2+A5x∗+A6=0,where Ak (k=0,1,…,6) are defined as follows:(20)A0=rd-mβ2a1+a2-v+d-mβ,A1=mβ-d23αdr+Kd+Kr-αr-mKβmβ-d+v-a1-a2-cb1r,A2=Kαd+rmβ-d3-3αdmβ-d2Kd+Kr-αr-mKβ-3α2d2rmβ-d+mβ-d2cb1Kd+Kr-αr-mKβ+Kβcb2-αKd+ra1+a2-v+αmβ-dmβ-3da1+a2-vKd+Kr-αr-mKβ-cb1r,A3=αKd+rcb1-3αdmβ-d2+rα3d3+Kcb1b2βmβ-d+vb1b2Kβ+αmβ-dKd+Kr-αr-mKβ3αd2+cb1mβ-3d+v-a1b2Kβmβ-d+rdα23d-2mβa1+a2-vKd+Kr-αr-mKβ-cb1r+αmβ-dmβ-3dKβcb2-αKd+ra1+a2-v,A4=α3d33Kd+rmβ-d-dKd+Kr-αr-mKβ+b2Kβα2mβ-dv-a1+cb1αkd+rmβ-dmβ-3dKβcb2+cb1Kr+Kd-aαr-mKβ+2a1-vb2Kαβmβ-d+2a2-vb1dα2Kd+r-mβ+b2Kαβv-a1α-2mβ-d,A5=cb1α3dKd+r3d-2mβ+b2Kdα3βa1-v+α3d3Kβcb2-αKd+ra1+a2-v+cb1Kr+Kd-αr-mKβ+b1v-a2α3Kd+rmβ-d-mβK-2α2dKd+r,A6=cb1α4Kd2d+r+a2-vb1α3Kd+r+Kβb1b2α3v-cd.Furthermore, it follows from the practical interpretations of interior equilibrium thatP∗ exists provided that y∗>0 and E∗>0, which derives that(21)rx∗2+αr+Kmβ-Kd-Krx∗-αKd+r<0,mβ-dx∗-αd>0.It follows from Routh-Hurwitz criterion [21], simple sufficient conditions for existence of P∗ of system (3) are as follows: (22)0<v<min⁡a1+a2+d-mβ,1-a2b2β-dαd+rb2β-d-r.Based on Jacobian matrix of system (3) evaluated around P∗ and Ξ defined in (5), characteristic equation of system (3) around P∗ is as follows:(23)λ2+B1λ+B2+B3λ+B4e-λτ1+B5e-λτ2=0,where Bi, i=1,2,…,5, are defined as follows:(24)B1=-βx∗y∗α+x∗2-a2b2E∗y∗b1+E∗x∗2+a1b1E∗x∗b2+E∗y∗2a1b1x∗b2+E∗y∗2+a2b2y∗b1+E∗x∗2-cb1+E∗x∗2b2+E∗y∗2,B2=a2b2βx∗y∗2E∗b1+E∗x∗2+a1b1βx∗y∗E∗α+x∗b2+E∗y∗2α+x∗2a1b1x∗b2+E∗y∗2+a2b2y∗b1+E∗x∗2-cb1+E∗x∗2b2+E∗y∗2,B3=rx∗K,B4=-ra2b2x∗y∗E∗b1+E∗x∗2Ka1b1x∗b2+E∗y∗2+a2b2y∗b1+E∗x∗2-cb1+E∗x∗2b2+E∗y∗2,B5=mαβx∗y∗α+x∗2βα+x∗-a2b2E∗b1+E∗x∗a1b1x∗b2+E∗y∗2+a2b2y∗b1+E∗x∗2-cb1+E∗x∗2b2+E∗y∗2.


### 3.1. Case I:τ1>0, τ2=0

Whenτ1>0 and τ2=0, (23) can be rewritten as follows:(25)λ2+B1λ+B2+B5+B3λ+B4e-λτ1=0.By substitutingλ=iω1 (ω1>0) into (25) and separating real and imaginary parts, it gives that(26)ω12-B2+B5=B4cos⁡ω1τ1+B3ω1sin⁡ω1τ1,B1ω1=-B3ω1cos⁡ω1τ1+B4sin⁡ω1τ1,which derives that(27)ω14+B12-2B2+B5-B32ω12+B2+B52-B42=0.According to Routh-Hurwitz criterion [21], a simple sufficient condition (B2+B5)2-B42<0 guarantees existence of positive root ω1∗ of (27). Hence, (25) has a pair of purely imaginary roots of the form ±iω1∗.By eliminatingsin⁡(ω1τ1) from (26), it can be calculated that τ1k∗ corresponding to ω1∗ is as follows:(28)τ1k∗=1ω1∗arccosB4-B1B3ω1∗2-B2+B5B4B42+B32ω1∗2+2kπω1∗,where k=0,1,….It follows from Butler’s lemma [22] that system (3) is locally stable around P∗ when 0<τ1<τ10∗,τ2=0.Based on the above analysis, local stability analysis of system (3) around P∗ is concluded in Theorem 3.Theorem 3.
If(B2+B5)2-B42<0, then system (3) is locally stable around P∗ when 0<τ1<τ10∗,τ2=0, and τ10∗ is defined in (28).


### 3.2. Case II:τ1=0, τ2>0

Whenτ1=0, τ2>0, by using similar arguments in Section 3.1, local stability analysis of system (3) around P∗ can be concluded in the following Theorem 4.Theorem 4.
If(B2+B4)2-B52<0, then system (3) is locally stable around P∗ when τ1=0, 0<τ2<τ20∗, and τ20∗ is defined as follows:(29)τ2k∗=1ω2∗arccosω2∗2-B2-B4B5+2kπω2∗,where k=0,1,… and ω2∗ is a positive root of the following equation:(30)ω24+B1+B32-2B2+B4ω22+B2+B42-B52=0.


### 3.3. Case III:τ1>0, τ2∈(0,τ20∗), τ1≠τ2

In this subsection,τ1 is considered to be bifurcation parameter and τ2 is regarded as fixed value τ~2∈(0,τ20∗), where τ20∗ is determined in Theorem 4.Letλ=iσ1 (σ1>0) be the root of (23), and two transcendental equations are obtained as follows:(31)σ12-B2+B5cos⁡σ1τ~2=B4cos⁡σ1τ1+B3σ1sin⁡σ1τ1,B1σ1-B5sin⁡σ1τ~2=B4sin⁡σ1τ1-B3σ1cos⁡σ1τ1,where Bk, k=1,2,…,5, are defined in (23).By eliminatingτ~2 from (31), it can be derived that(32)cos⁡σ1τ1=f10σ1,τ~2f12σ1,τ~2,sin⁡σ1τ1=f11σ1,τ~2f12σ1,τ~2,where f1k(σ1,τ~2), k=0,1,2, are defined as follows:(33)f10σ1,τ~2=B4-B1B3σ12-B2B4-B4B5cos⁡σ1τ~2+B3B5σ1sin⁡σ1τ~2,f11σ1,τ~2=B3σ13+B1B4-B2B3σ1-B3B5σ1cos⁡σ1τ~2-B4B5sin⁡σ1τ~2,f12σ1,τ~2=B42+B32σ12.Based onf1k(σ1,τ~2), k=0,1,2, defined in (32), we define(34)F1σ1,τ~2=f102σ1,τ~2+f112σ1,τ~2-f122σ1,τ~2=0.It is difficult to discuss root properties of transcendental equation (34) due to its complicated form. Without investigating the properties of roots of (34) in detail, it follows from simple computations that (23) has a pair of purely imaginary roots when (34) has at least two roots 0<σ10<σ11<⋯.By denotingσ1w=max⁡{σ1k,k=0,1,…} and τ1 which is considered to be bifurcation parameter, we can obtain the corresponding critical value τ1w:(35)τ1k=ψ1w+2kπσ1w,τ1w=min⁡τ1k,k=0,1,…,where k=0,1,…,ψ1c∈[0,2π) satisfies the following equations:(36)cos⁡ψ1w=f10σ1w,τ~2f12σ1w,τ~2,sin⁡ψ1w=f11σ1w,τ~2f12σ1w,τ~2.Based on the above analysis, local stability analysis aroundP∗ due to dynamical variation of bifurcation parameter τ1 can be summarized as follows.Theorem 5.
If (34) has at least two roots 0<σ10<σ11<⋯, and there exists τ1w defined in (35), then system (3) is locally stable around P∗ when 0<τ1<τ1w,τ2=τ~2.


### 3.4. Case IV:τ2>0, τ1∈(0,τ10∗), τ1≠τ2

In this subsection,τ2 is considered to bifurcation parameter and τ1 is regarded as fixed value τ~1∈(0,τ10∗), where τ10∗ is determined in Theorem 3. By using similar arguments in Theorem 5 of this paper, local stability analysis around P∗ can be concluded in Theorem 6.Firstly, we define following equation:(37)F2σ2,τ~1=f202σ2,τ~1+f212σ2,τ~1-f222σ2,τ~1=0,where f2k(σ2,τ~2), k=0,1,2, are defined as follows:(38)f20σ2,τ~1=σ22-B2-B4cos⁡σ2τ~1-B3σ2sin⁡σ2τ~1,f21σ2,τ~1=B1σ2+B3σ2cos⁡σ2τ~1-B4sin⁡σ2τ~1,f22σ2,τ~1=B5.Theorem 6.
If (37) has at least two roots 0<σ20<σ21<⋯, there exists a critical delay τ2w, and then system (3) is locally stable around P∗ when 0<τ2<τ2w,τ1=τ~1, where τ2w=min⁡{ψ2w+2kπ/σ2w,k=0,1,…},σ2w=max⁡{σ2k,k=0,1,…} and ψ2w∈[0,2π) satisfies the following equations:(39)cos⁡ψ2w=f20σ2w,τ~1f22σ2w,τ~1,sin⁡ψ2w=f21σ2w,τ~1f22σ2w,τ~1.


## 3.1. Case I:τ1>0, τ2=0

Whenτ1>0 and τ2=0, (23) can be rewritten as follows:(25)λ2+B1λ+B2+B5+B3λ+B4e-λτ1=0.By substitutingλ=iω1 (ω1>0) into (25) and separating real and imaginary parts, it gives that(26)ω12-B2+B5=B4cos⁡ω1τ1+B3ω1sin⁡ω1τ1,B1ω1=-B3ω1cos⁡ω1τ1+B4sin⁡ω1τ1,which derives that(27)ω14+B12-2B2+B5-B32ω12+B2+B52-B42=0.According to Routh-Hurwitz criterion [21], a simple sufficient condition (B2+B5)2-B42<0 guarantees existence of positive root ω1∗ of (27). Hence, (25) has a pair of purely imaginary roots of the form ±iω1∗.By eliminatingsin⁡(ω1τ1) from (26), it can be calculated that τ1k∗ corresponding to ω1∗ is as follows:(28)τ1k∗=1ω1∗arccosB4-B1B3ω1∗2-B2+B5B4B42+B32ω1∗2+2kπω1∗,where k=0,1,….It follows from Butler’s lemma [22] that system (3) is locally stable around P∗ when 0<τ1<τ10∗,τ2=0.Based on the above analysis, local stability analysis of system (3) around P∗ is concluded in Theorem 3.Theorem 3.
If(B2+B5)2-B42<0, then system (3) is locally stable around P∗ when 0<τ1<τ10∗,τ2=0, and τ10∗ is defined in (28).


## 3.2. Case II:τ1=0, τ2>0

Whenτ1=0, τ2>0, by using similar arguments in Section 3.1, local stability analysis of system (3) around P∗ can be concluded in the following Theorem 4.Theorem 4.
If(B2+B4)2-B52<0, then system (3) is locally stable around P∗ when τ1=0, 0<τ2<τ20∗, and τ20∗ is defined as follows:(29)τ2k∗=1ω2∗arccosω2∗2-B2-B4B5+2kπω2∗,where k=0,1,… and ω2∗ is a positive root of the following equation:(30)ω24+B1+B32-2B2+B4ω22+B2+B42-B52=0.


## 3.3. Case III:τ1>0, τ2∈(0,τ20∗), τ1≠τ2

In this subsection,τ1 is considered to be bifurcation parameter and τ2 is regarded as fixed value τ~2∈(0,τ20∗), where τ20∗ is determined in Theorem 4.Letλ=iσ1 (σ1>0) be the root of (23), and two transcendental equations are obtained as follows:(31)σ12-B2+B5cos⁡σ1τ~2=B4cos⁡σ1τ1+B3σ1sin⁡σ1τ1,B1σ1-B5sin⁡σ1τ~2=B4sin⁡σ1τ1-B3σ1cos⁡σ1τ1,where Bk, k=1,2,…,5, are defined in (23).By eliminatingτ~2 from (31), it can be derived that(32)cos⁡σ1τ1=f10σ1,τ~2f12σ1,τ~2,sin⁡σ1τ1=f11σ1,τ~2f12σ1,τ~2,where f1k(σ1,τ~2), k=0,1,2, are defined as follows:(33)f10σ1,τ~2=B4-B1B3σ12-B2B4-B4B5cos⁡σ1τ~2+B3B5σ1sin⁡σ1τ~2,f11σ1,τ~2=B3σ13+B1B4-B2B3σ1-B3B5σ1cos⁡σ1τ~2-B4B5sin⁡σ1τ~2,f12σ1,τ~2=B42+B32σ12.Based onf1k(σ1,τ~2), k=0,1,2, defined in (32), we define(34)F1σ1,τ~2=f102σ1,τ~2+f112σ1,τ~2-f122σ1,τ~2=0.It is difficult to discuss root properties of transcendental equation (34) due to its complicated form. Without investigating the properties of roots of (34) in detail, it follows from simple computations that (23) has a pair of purely imaginary roots when (34) has at least two roots 0<σ10<σ11<⋯.By denotingσ1w=max⁡{σ1k,k=0,1,…} and τ1 which is considered to be bifurcation parameter, we can obtain the corresponding critical value τ1w:(35)τ1k=ψ1w+2kπσ1w,τ1w=min⁡τ1k,k=0,1,…,where k=0,1,…,ψ1c∈[0,2π) satisfies the following equations:(36)cos⁡ψ1w=f10σ1w,τ~2f12σ1w,τ~2,sin⁡ψ1w=f11σ1w,τ~2f12σ1w,τ~2.Based on the above analysis, local stability analysis aroundP∗ due to dynamical variation of bifurcation parameter τ1 can be summarized as follows.Theorem 5.
If (34) has at least two roots 0<σ10<σ11<⋯, and there exists τ1w defined in (35), then system (3) is locally stable around P∗ when 0<τ1<τ1w,τ2=τ~2.


## 3.4. Case IV:τ2>0, τ1∈(0,τ10∗), τ1≠τ2

In this subsection,τ2 is considered to bifurcation parameter and τ1 is regarded as fixed value τ~1∈(0,τ10∗), where τ10∗ is determined in Theorem 3. By using similar arguments in Theorem 5 of this paper, local stability analysis around P∗ can be concluded in Theorem 6.Firstly, we define following equation:(37)F2σ2,τ~1=f202σ2,τ~1+f212σ2,τ~1-f222σ2,τ~1=0,where f2k(σ2,τ~2), k=0,1,2, are defined as follows:(38)f20σ2,τ~1=σ22-B2-B4cos⁡σ2τ~1-B3σ2sin⁡σ2τ~1,f21σ2,τ~1=B1σ2+B3σ2cos⁡σ2τ~1-B4sin⁡σ2τ~1,f22σ2,τ~1=B5.Theorem 6.
If (37) has at least two roots 0<σ20<σ21<⋯, there exists a critical delay τ2w, and then system (3) is locally stable around P∗ when 0<τ2<τ2w,τ1=τ~1, where τ2w=min⁡{ψ2w+2kπ/σ2w,k=0,1,…},σ2w=max⁡{σ2k,k=0,1,…} and ψ2w∈[0,2π) satisfies the following equations:(39)cos⁡ψ2w=f20σ2w,τ~1f22σ2w,τ~1,sin⁡ψ2w=f21σ2w,τ~1f22σ2w,τ~1.


## 4. Optimal Control Problem

In this section,τ1 and τ2 are considered to be a fixed value τ^1∈(0,τ1w) and τ^2∈(0,τ20∗), respectively; τ^1≠τ^2, τ1w, and τ20∗ have been determined in Theorems 5 and 4, respectively. It follows from Theorem 5 that system (3) with τ^1 and τ^2 is locally stable around P∗, and we will investigate optimal control E(t) with maximizing total discounted net revenues and minimizing commercial harvesting cost during [t0,tf]. The optimal harvest control problem is formulated as follows:(40)JE=∫t0tfe-δtEtxta1b1+Etxt+Etyta2b2+Etyt-cEtdt,where δ is instantaneous annual discount rate [7], and optimal control problem (40) is subject to the first and second delayed differential equations of system (3) with τ^1 and τ^2. Furthermore, convexity of J(E) with respect to E(t), linearity of the delayed differential equations with E(t), and range values compactness of state variables can be uniformly investigated when discussing existence of optimal harvest control E(t). By assuming Eδ as an optimal harvest control with corresponding optimal state variables xδ and yδ of system (3) with τ^1 and τ^2, and Eδ satisfies the following condition:(41)JEδ=max⁡JEt∣Et∈U,where U represents control set that is defined as follows: (42)U=E:t0,tf⟶0,Emax∣EisLebesguemeasurable..The Hamiltonian functionH(t) associated with system (3) with τ^1 and τ^2 is(43)Ht=e-δtEtxta1b1+Etxt+Etyta2b2+Etyt-cEt+λ1trxt1-xτ^1K-βxtytα+xt-Etxt+λ2tytmβxτ^2α+xτ^2-d-Etyt,where xτ^1≔x(t-τ^1), xτ^2≔x(t-τ^2), and λ1(t),λ2(t) are adjoint functions.Theorem 7.
There exists an optimal harvestEδ which is determined by (44); furthermore, there exist adjoint functions λ1(t),λ2(t) satisfying system (47) with λk(tf)=0,k=1,2.(44)a1b1xδb1+Eδxδ2+a2b2yδb2+Eδyδ2=c+xδλ1t+yδλ2teδt.Proof.
By differentiatingH(t), the adjoint system can be obtained as follows:(45)dλ1dt=-∂H∂xt-χ1t0,tf-τ^1t∂H∂xτ^1t+τ^1-χ2t0,tf-τ^2t∂H∂xτ^2t+τ^2,dλ2dt=-∂H∂yt,where(46)χ1t0,tf-τ^1t=1,t∈0,tf-τ^1,0,otherwise,χ2t0,tf-τ^2t=1,t∈0,tf-τ^2,0,otherwise.
By using the optimal conditions, (45) can be rewritten as follows:(47)dλ1dt=-e-δta1b1Etb1+Etxt2+λ1trxtK-1+αβytα+xt2+Et+χ1t0,tf-τ^1trxtλ2t+τ^1K-χ2t0,tf-τ^2tmαβytα+xt2λ2t+τ^2,dλ2dt=-e-δta2b2Etb2+Etyt2+βxtλ1tα+xt+λ2tEt+d-mβxtα+xt.
In optimal control problem (40), there is not terminal cost and final state is free. Consequently, transversality conditions (bounded conditions) for adjoint functions are as follows: λk(tf)=0,k=1,2.
It is easy to obtain characterization of optimal harvest controlEδ, which derives that ∂H/∂EE=Eδ=0. Consequently, it follows from simple computations that there exists an optimal harvest Eδ determined by(48)a1b1xδb1+Eδxδ2+a2b2yδb2+Eδyδ2=c+xδλ1t+yδλ2teδt,where λ1(t),λ2(t) satisfies system (47) with λk(tf)=0,k=1,2.Remark 8.
It follows from similar arguments in Theorem7 that symmetric optimal control problem on cases in Section 3.4 can be also obtained, which are omitted in this paper. τ1 and τ2 can be set to be a fixed value τ^1∈(0,τ10∗) and τ^2∈(0,τ2w), respectively; τ^1≠τ^2, τ10∗, and τ2w are determined in Theorems 3 and 6, respectively.In the following part, the optimal system is numerically solved based on an iterative method with Runge-Kutta fourth order scheme. Firstly, by assuming that there exist a step sizeρ>0 and integers i,j,n with τ1=iρ,τ2=jρ, and tf-t0=nρ. Secondly, we initialize optimal harvest control with a given initial value E0=0.09. By using forward difference approximation, we solve the first and second differential equation of system (3) with initial conditions x0=0.38, y0=0.43. Thirdly, by using backward difference approximation, we solve the adjoint functions (47) with transversality conditions. Finally, the optimal harvest control is updated by values of the state and adjoint variables, and the updated optimal harvest control is replaced with the initial value given in the second step. This process will be repeated until successive iterates of harvest control values become sufficiently close.By utilizing combinations of forward and backward difference approximations, it can be derived that fork=0,…,n-1(49)xk+1=xk+ρrxk1-xk-iK-βxkykα+xk-Ekxk,yk+1=yk+ρykmβxk-jα+xk-j-d-Ekyk,λ1n-k-1=λ1n-k-ρe-δn-ka1b1Ek+1b1+Ek+1xk+12+ρλ1n-krxk+1K-1+αβyk+1α+xk+12+Ek+1+ρχ1t0,tf-τ^1λ1n-k+irxk-jK-χ2t0,tf-τ^2λ2n-k+jmαβyk+1α+xk+12,λ2n-k-1=λ2n-k-ρe-δn-ka2b2Ek+1b2+Ek+1yk+12-βxk+1λ1n-kα+xk+1-λ2n-kEk+1+d-mβxk+1α+xk+1.In order to show effectiveness of optimal harvest control designed in Theorem7, parameters values of system (3) are partially from [5], r=2.5, K=1.4, m=0.6, β=1, α=1.4, d=0.09, a1=1, a2=0.5, b1=1.2, b2=1, c=0.05, δ=0.01, t0=0, and tf=120 with appropriate units. By using given values of parameters and simple computations, interior equilibrium P∗(0.2471,3.3911,0.7467) exists provided that 0<v<0.1491. In the following part, v=0.075 is utilized in the numerical simulations, which is arbitrarily selected within (0,0.1491) and is enough to merit theoretical analysis obtained in this paper. It follows from Theorem 3, Theorem 7, and simple computations that system (3) is locally stable around (0.2471,3.3911,0.7467) when 0<τ1<τ10∗=0.5632 and τ2=0. By using similar arguments and Theorems 4 and 7, it is revealed that system (3) is locally stable around (0.2471,3.3911,0.7467) when 0<τ2<τ20∗=0.3894 and τ1=0. As analyzed in Section 3.3, τ^2=0.15 is arbitrarily selected within (0,τ20∗), that is, (0,0.3894), and it follows from Theorems 5 and 7 that system (3) is locally stable around (0.2471,3.3911,0.7467) when 0<τ1<τ1w=0.4126 and τ^2=0.15. As analyzed in Section 3.4, τ^1=0.28 is arbitrarily selected within (0,τ10∗), that is, (0,0.5632), and it follows from Theorems 6 and 7 that system (3) is locally stable around (0.2471,3.3911,0.7467) when τ^1=0.28 and 0<τ2<τ2w=0.2219. Numerical simulations are made to investigate dynamic effects of double time delays on variation of optimal prey, predator biomass, and optimal harvest control, which can be found in Figures 1 and 2, respectively. The initial values of simulation work made in Figure 1 are as follows: x(0)=0.28,y(0)=0.31,E(0)=0.76, and initial values of simulation work made in Figure 2 are as follows: x(0)=0.42,y(0)=0.21,E(0)=0.81.Figure 1
Stable dynamical variation of optimal prey, predator biomass and optimal harvest control with the increasing time underτ1=0.15,0.25,0.35 and τ2=0.15, respectively.Figure 2
Stable dynamical variation of optimal prey, predator biomass and optimal harvest control with the increasing time underτ1=0.28 and τ2=0.1,0.15,0.2, respectively.


## 5. Conclusion

By incorporating maturation delay for prey population and gestation delay for predator population, we extend work done in [5]. In this paper, a singular prey-predator fishery model with double time delays is established. Fishing efforts are introduced to commercially harvest prey population and predator population, which are developed as control instruments to investigate optimal utilization of prey-predator fishery resource. Positivity of solutions and uniform persistence of system are discussed in Theorems 1 and 2. By analyzing associated characteristic transcendental equation, local stability around interior equilibrium is discussed due to combined variations of double time delays, which can be found in Theorem 3 to Theorem 6. Furthermore, an optimal harvest control is designed in Theorem 7. By using the proposed mathematical model and dynamical analysis, authors aim to obtain some results which are theoretically beneficial to discussing optimal harvesting strategies and sustainability mechanism of harvested prey-predator fishery system with double time delays. Furthermore, the theoretical results may be potentially constructive for administrative agencies to formulate regulatory policies to preserve economic optimality through ensuring global sustainability.



---

*Source: 4398527-2016-11-21.xml*