# Discretization Fractional-Order Biological Model with Optimal Harvesting

Authors: Barhoom Noor S. Sh.; Sadiq Al-Nassir
Journal: Abstract and Applied Analysis
 (2022)

---

## Abstract

In this paper, a discretization of a three-dimensional fractional-order prey-predator model has been investigated with Holling type III functional response. All its fixed points are determined; also, their local stability is investigated. We extend the discretized system to an optimal control problem to get the optimal harvesting amount. For this, the discrete-time Pontryagin’s maximum principle is used. Finally, numerical simulation results are given to confirm the theoretical outputs as well as to solve the optimality problem.


---

## Body

## 1. Introduction

Recently, population models have been given attention by researchers due to their importance and their relation to real life. In fact, there are different processes to describe population models, namely, ordinary and partial differential equations, difference equations, and fractional-order differential equations. Fractional-order derivative systems play a significant role in real-life problems due to their ability to provide accurate descriptions of different linear and nonlinear phenomena. We refer to the most recent works that are done by many researchers [1–5]. The mathematical behavior of fractional-order prey-predator models is also discussed and investigated in [6–10] and references therein. Discrete-time mathematical models hit many branches of pure and applied mathematics; also, they depict many real situations in life due to exhibiting and producing more plentiful dynamical behavior than those that have been seen in continuous time models of the same types. These types of models can be used to study and investigate problems in biology, economics, engineering, etc. For more details, see [11–17].In the last decades, there are many articles that deal with the subject of harvesting in different approaches for details (see [18–21] and references therein). We also refer to the excellent book that is written by Agrawal et al. [22]. In this book, recent important results of infection diseases and methods with new approximations are analyzed and developed. In [23], the authors investigated and discussed coronavirus disease in a mathematical model that deals with coronavirus, and an effective control policy for the COVID-19 outbreak is set.The study of biological models with harvesting, in particular, predator-prey models, is an important subject in the commercial harvesting industry and for many scientific communities including biology, ecology, and economics. The fundamental problem in commercial exploitation of renewable resources is to get the optimal gain between present and future harvests, so that employing optimal control theory in this field is a very useful tool that can help to make decisions in various problems in real world including biological situations [24, 25].This work is organized as follows. In Section2, a fractional-order system is discretized and all its equilibrium points are found. In Section 3, the local stability of all its equilibria is discussed. Section 4 deals with the employment of the discrete-time Pontryagin’s maximum principle to solve the optimal control solution and its corresponding optimal state solutions. In Section 5, numerical simulations are done to confirm the theoretical outcomes. Finally, conclusion is given in Section 6.


## 2. The Discrete-Time Model and Its Fixed Points

First, the authors considered and analyzed in [26] the following fractional-order three-dimensional prey-predator model with the functional response of Holling type III:
(1)Dαxt=x1−x−ax2y1+ax2,Dαyt=bx2y1+ax2−dy2z1+dy2−m1y−hy,Dαzt=cy2z1+dy2−m2z,whereX, Y, and Z are the density of the prey population, the middle predator population, and the top predator population. m1, m2, a, and b are the death rates for the middle predator population and top predator population, the capture rate of prey and middle predator, and the conservation rate of prey X, respectively. The parameters d and c are the conservation rates of middle predator Y to the top predatorZ. The parameter h represents the harvesting rate. For more details about the interpretations of the parameters, see [26]. Here, Dα is called the α-order Caputo differential operator [27, 28].In this work, we investigate the dynamic behavior of the discretization of the system (1); then, we extend the discretized system to the optimal control problem in the next section. The discretization process of the system (1) is given as follows:Assume thatx0=x0, y0=y0, and z0=z0 are the initial conditions of system (1), so that the discretization of the system (1) with piecewise constant arguments is given as
(2)Dαxt=xt/ss1−xt/ss−axt/ss2yt/ss1+axt/ss2,Dαyt=bxt/ss2yt/ss1+axt/ss2−dyt/ss2zt/ss1+dyt/ss2−m1yt/ss−hyt/ss,Dαzt=cyt/ss2zt/ss1+dyt/ss2−m2zt/ss.First, lett∈0,s, so t/s∈0,1. Thus, we obtain
(3)Dαxt=x01−x0−ax02y01+ax02,Dαyt=bx02y01+ax02−dy02z01+dy02−m1y0−hy0,Dαzt=cy02z01+dy02−m2z0.The solution of (3) is reduced to
(4)x1t=x0+Iαx01−x0−ax02y01+ax02=x0+tααΓαx01−x0−ax02y01+ax02,y1t=y0+Iαbx02y01+ax02−dy02z01+dy02−m1y0−hy0=y0+tααΓαbx02y01+ax02−dy02z01+dy02−m1y0−hy0,z1t=z0+Iαcy02z01+dy02−m2z0=z0+tααΓαcy02z01+dy02−m2z0.Second, lett∈s,2s, so t/s∈1,2. Hence, we get the following system:
(5)Dαxt=x11−x1−ax12y11+ax12,Dαyt=bx12y11+ax12−dy12z11+dy12−m1y1−hy1,Dαzt=cy12z11+dy12−m2z1,which have the following solution:(6)x2t=x1s+Isαx1s1−x1s−ax12sy1s1+ax12s=x1s+t−sααΓαx11−x1−ax12sy1s1+ax12,y2t=y1s+Isαbx12sy1s1+ax12s−dy12sz1s1+dy12s−m1y1s−hy1s=y1s+t−sααΓαbx12sy1s1+ax12s−dy12sz1s1+dy12s−m1y1s−hy1s,z2t=z1s+Isαcy12sz1s1+dy12s−m2z1s=z1s+t−sααΓαcy12sz1s1+dy12s−m2z1s,whereIsα≡1/Γα∫stt−τα−1dτ,α>0. Thus, after repeating the discretization process n times, we obtain
(7)xn+1t=xnns+t−nsααΓαxnns1−xnns−axn2nsynns1+axn2ns,yn+1t=ynns+t−nsααΓαbxn2nsynns1+axn2ns−dyn2nsz1s1+dyn2ns−m1ynns−hynns,zn+1t=z1ns+t−nsααΓαcyn2nsznns1+dyn2ns−m2znns,wheret∈ns,n+1s. For t⟶n+1s, then the previous system is reduced to
(8)xn+1=xn+sααΓαxn1−xn−axn2yn1+axn2,yn+1=yn+sααΓαbxn2yn1+axn2−dyn2zn1+dyn2−m1yn−hyn,zn+1=zn+sααΓαcyn2zn1+dyn2−m2zn.In order to get the fixed points of the considered system (8), we have to solve the following simultaneous equations:
(9)x=x+sααΓαx1−x−ax2y1+axn2,y=y+sααΓαbx2y1+ax2−dy2z1+dy2−m1y−hy,z=z+sααΓαcy2z1+dy2−m2z.Thus, all possible fixed points of system (8) are as follows:
(1)
The trivial fixed pointE0=0,0,0 and the fixed pointE1=1,0,0 exist without any restriction on the parameters(2)
The top predator free fixed pointE2=x¯,y¯,0 exists if b>am1+h and 0<x¯<1 where x¯=m1+h/b−am1+h and y¯=1/ax1−x¯1+ax¯2(3)
The interior fixed pointE3=x∗,y∗,z∗ exists if c>m3d and bx∗2/1+ax∗2>m1+h where y∗=m2/c−m2d, x∗ is the positive root of the equation ax3−ax2+ay∗+1x−1=0, and z∗=1+dy∗2/dy∗bx∗2/1+ax∗2−m1+hThen, the Jacobian matrixJ of the system (8) evaluated at any point x,y,z is given by
(10)x,y,z=1+m−2mx−2amxyk2−max2k02mbxyk21+mbx2k−2dyzk12−m1+h−mdy2k10cmyzk121+cmy2k1−mm2,wherek=1+ax2, k1=1+dy2, and m=sα/αΓα. Therefore, the characteristic equation of the Jacobian matrix is given by
(11)pλ=λ3+a2λ2+a1λ+a0,a2=−A1+A2+A3a1=A1A3+A2A3+A1A2+A5+A4,a0=−A1A2A3−A1A5−A3A4=−w−A3A4withw=A1A2A3+A1A5,A1=1+m−2mx−2maxyk2,A2=1+mbx2k−2dyzk12−m1+h,A3=1+cmy2k1−mm2,A4=max2k2mbxyk2,A5=2cdm2y3zk13.In order to study stability analysis of the fixed points of the model (8), we give Definition 1 and Lemma 2.Definition 1 (see [29]).
A fixed pointx∗ of system x⃑t+1=fx⃑t is said to be as follows:
(1)
Saddle point, if at least one eigenvalue of the Jacobian matrixJx∗ is in modulus>1(2)
Sink point, if all eigenvalues of the Jacobian matrixJx∗are in modulus<1(3)
Source point, if all eigenvalues of the Jacobian matrixJx∗are in modulus>1(4)
Nonhyperbolic point, if at least one eigenvalue of the Jacobian matrixJx∗ is in modulus=1Lemma 2 (see [29]).
LetFλ=λ2+pλ+q be a polynomial of degree two. Suppose that F1>0 and λ1,λ2 are the roots of F; then
(1)
λ1<1 and λ2<1, if and only if F−1>0 and q<1(2)
λ1>1 and λ2>1, if and only if F−1>0 and q>1(3)
λ1>1 and λ2<1 (or λ1<1 and λ2>1), if and only if F−1<0(4)
λ1=−1 and λ2≠1, if and only if F−1=0 and p≠0,2Theorem3 discusses the dynamical behavior of the fixed points E0, E1, E2, and E3 of the system (8), respectively.Theorem 3.
For the system (8), the fixed point E0 can be classified as follows:
(1)
The pointE0 will never be a sink point(2)
The pointE0 is a source, if and only if m>max2/m1+h,3/m2(3)
The pointE0 is a saddle point, if and only if m<min2/m1+h,3/m2(4)
The pointE0 is a nonhyperbolic point, if and only if m=2/m1+h or m=3/m2, where m=sα/αΓα, m1=c1k02R0/a1, and m2=R0/c2.Proof.
It is clear that the Jacobian matrix of system (8) at E0 can be written as follows:
(12)JE0=1+m0001−mm1+h0001−mm2.
The characteristic equation of the Jacobian matrixJE0 is given by
(13)pλ=1+m−λ1−mm1+h−λ1−mm2−λ.
So that the eigenvalues ofJE0 areλ1=1+m, λ2=1−mm1+h, andλ3=1−mm2. It is clear that λ1 is always greater than 1. Hence, the fixed point E0 will never be a sink point. It is clear that if 0<m<2/m1+h, then 1−mm1+h<1. Therefore, λ2<1. Now, if m<3/m2, then 1−mm2<1. Hence, λ3<1. Therefore, all results can be obtained.Theorem 4.
For the system (8), the fixed point E1 can be classified as follows:
(1)
E1 is a sink point if and only if m∈I1∩I2(2)
E1 is a source point if and only if m∈I3(3)
E1 is a saddle point if and only if one of the following holds
(14)m∈I4∩I2,m∈I1∩I5,m∈I1∩I6.:(4)
E1 is a nonhyperbolic point if and only if m=2 or m=3/m2 or m=3/m2, where I1=0,min2,3/m2, I2=1+am1+h−2/b,m1+h1+a/b, I3=max2,3/m2,m1+h1+a/b,∞, I4=max2,3/m2,∞, I5=m1+h1+a/b,∞, and I6=0,1+am1+h−2/b.Proof.
The Jacobian matrix of system (8) at E1 is
(15)JE1=1−m−ma1+a001+mb1+a−m1+h0001−mm2.
Therefore, the characteristic equation of the Jacobian matrixJE1 is given by
(16)pλ=1+m−λ1+mb1+a−m1+h−λ1−mm2−λ.
So that the eigenvalues ofJE1 are λ1=1−m, λ2=1+mb/1+a−m1+h, and λ3=1−mm2. Now, if <m<2, then 1−m<1. This implies that λ1<1 if 1+am1+h−2/b<m<m1+h1+a/b. This gives 1+mb/1+a−m1+h<1. Therefore, λ2<1. Assume that 0<m<2/m3. This implies that 1−mm3<1 and λ3<1. Therefore, all results are got.Theorem 5.
For the system (8), the fixed point E2 can be classified as follows:
(1)
E2 is a sink point, if and only if the following conditions hold:
(17)2mau+m2au−m2a2u21−x>m−2b,2au+2mabu−2ma2u21−x<b.(2)
E2 is a source if and only if 2au+2mabu−2ma2u21−x>b(3)
E2 is a saddle point if and only if 2mau+m2au−m2a2u21−x<m−2b(4)
E2 is a nonhyperbolic point if and only if 2mau+m2au−m2a2u21−x=m−2b or 2au+2mabu−2ma2u21−x=b, where =m1+h/b−am1+h, y=1/ax1−x1+ax2, k=1+ax2, and u=m1+h.Proof.
It is clear that the Jacobian matrix of system (8) at E2 can be written as follows:
(18)JE2=1−m+2mau1−xb−maub02m1−xb−au1−db1−x2a2u2+db21−x2001+mcy2k1−mm2.
Therefore, the characteristic polynomial of the Jacobian matrix atE2 is given by
(19)fλ=1−m+2mau1−xb−λFλ=0,
whereFλ=λ2+pλ+q such that p=2+m−2mau1−x/b and q=1−m+2mau1−x/b+2m2au1−x−2m2a2u2/b1−x.
Now, the first root of the characteristic polynomial isλ1=1+mcb21−x2/a2u2+db21−x2−mm2. If 0<m<2a2u2+db21−x2/c−dm2b21−x2−m2a2u2, then, we have 1+mcb21−x2/a2u2+db21−x2−mm2<1 and λ1<1.
In order to investigate the other roots, we use Lemma2 so that F1>0 if and only if 1+m−2−m+2mau1−x/b+1−m+2mau1−x/b+2m2au1−x−2m2a2u2/b1−x>0. Thus, F1>0 if and only if 2m2au1−x−2m2a2u2/b1−x>0. Since 0<x<1, we get that F1 is always greater than zero.
Now, we have to show thatF−1>0 and q<1. If 2mau+m2au−m2a2u21−x>m−2b, then 2mau/b+m2au−m2a2u2/b1−x>m−2. Hence, 2−m+2mau1−x/b+m2au1−x−m2a2u2/b1−x>0.
Therefore,4−2m+4mau1−x/b+2m2au1−x−2m2a2u2/b1−x>0.
This implies that1−m+2+2mau1−x/b+1−m+2mau1−x/b+2m2au1−x−2m2a2u2/b1−x>0 and F−1>0.
Now, if2au+2mabu−2ma2u21−x<b, then 2au/b+2mau−2ma2u2/b1−x<1. Hence, 2mau/b+2m2au−2m2a2u2/b1−x<m. Therefore, −m+2mau1−x/b+2m2au1−x−2m2a2u2/b1−x<0. This implies that 1−m+2mau1−x/b+2m2au1−x−2m2a2u2/b1−x<1 andq<1. Therefore, the results can be obtained.Before we study the behavior of the interior fixed point of the system (8), we need Lemma 6.Lemma 6 (see [29]).
Letpλ=λ3+a2λ2+a1λ+a0 be a cubic polynomial with real coefficients. Then, the necessary and sufficient conditions that all roots lie in the open disk λ<1 are as follows:
(20)a2+a0<1+a1,a2−3a0<3−a1,a02+a1−a0a2<1.Theorem 7.
For the system (8), the fixed point E3 is a stable point if the following conditions hold:
(21)c1<A3A4<c2,d1<A3A4<d2,A32A42+2w+a2A3A4+kk<0,
wherea2=−A1+A2+A3, kk=w2+a1+a2w−1, a1=A1A3+A2A3+A1A2+A5+A4, a0=−A1A2A3−A1A5−A3A4=−w−A3A4suchthatw=A1A2A3+A1A5, A1=1+m−2mx−2maxy/k2, A2=1+mbx2/k−2dyz/k12−m1+h, A3=1+cmy2/k1−mm2, A4=max2/k2mbxy/k2, A5=2cdm2y3z/k13, c1=−1−A1A3+A2A3+A1A2−A5−A4−A1+A2+A3−A1A2A3−A1A5, c2=1+A1A3+A2A3+A1A2+A5+A4−A1+A2+A3−A1A2A3−A1A5, d1=1/3−3+A1A3+A2A3+A1A2+A5+A4+A1+A2+A3−3A1A2A3−3A1A5, and d2=1/33−A1A3+A2A3+A1A2−A5−A4+A1+A2+A3−3A1A2A3−3A1A5.Proof.
The proof is based on Lemma6.
The Jacobian matrix of system (8) at E3 can be written as follows:
(22)JE3=1+m−2mx−2amxyk2−max2k02mbxyk21+mbx2k−2dyzk12−m1+h−mdy2k102cmyzk121+cmy2k1−mm2.
Therefore, the characteristic polynomial of the Jacobian matrixJE3 at E3 is
(23)pλ=λ3+a2λ2+a1λ+a0.
Now, letc1<A3A4<c2; this gives
(24)−1−A1A3+A2A3+A1A2−A5−A4−A1+A2+A3−A1A2A3−A1A5<A3A4<1+A1A3+A2A3+A1A2+A5+A4−A1+A2+A3−A1A2A3−A1A5.
Hence,−1−A1A3+A2A3+A1A2−A5−A4+A1+A2+A3+A1A2A3+A1A5<−A3A4<1+A1A3+A2A3+A1A2+A5+A4+A1+A2+A3+A1A2A3+A1A5, and −1−A1A3+A2A3+A1A2−A5−A4<−A1+A2+A3−A1A2A3−A1A5−A3A4<1+A1A3+A2A3+A1A2+A5+A4.
Therefore,−1−a1<a2+a0<1+a1 and condition (1) in Lemma 6 hold.
Now, assume thatd1<A3A4<d2. This gives
(25)13−3+A1A3+A2A3+A1A2+A5+A4+A1+A2+A3−3A1A2A3−3A1A5<A3A4<133−A1A3+A2A3+A1A2−A5−A4+A1+A2+A3−3A1A2A3−3A1A5.
Hence,−3+A1A3+A2A3+A1A2+A5+A4+A1+A2+A3−3A1A2A3−3A1A5<3A3A4<3−A1A3+A2A3+A1A2−A5−A4+A1+A2+A3−3A1A2A3−3A1A5, and −3+A1A3+A2A3+A1A2+A5+A4<−A1+A2+A3+3A1A2A3+3A1A5+3A3A4<3−A1A3+A2A3+A1A2−A5−A4.
Therefore,−3+a1<a2−3a0<3−a1 and condition (2) in Lemma 6 hold.
Finally, if we haveA32A42+2w+a2A3A4+kk<0, we get
(26)A32A42+2wA3A4+w2+a1+a2w+a2A3A4−1<0,and −w−A3A42+a1−a2−w−A3A4<1.
Hence, condition (3) in Lemma 6 holds. Therefore, all roots of the characterization polynomial are in the unit circle, and the point E3 is locally stable.


## 3. Optimal Harvesting Strategy

In this section, we study the optimal control strategy to get the optimal harvesting amount; this will be done through the discrete version of Pontryagin’s maximum principle, so that the form of the objective functional that we have to minimize is(27)Jht=∑t=0T−1c1htyt.The parameterht refers to the control variable, which represents the rate of harvesting at time t, with 0≤ht≤hMax, and hMax is the maximum harvesting. c1 is a positive constant. In this optimal control, the state equations are
(28)xt+1=xt+sααΓαxt1−xt−axt2yt1+axt2,yt+1=yt+sααΓαbxt2yt1+axt2−dyt2zt1+dyt2−m1yt−htyt,zt+1=zt+sααΓαcyt2zt1+dyt2−m2zt.The parameterss, a, b, d, m1, and m2 are defined as before, and they have the same interpretations.To solve the problem, we have to set the Hamiltonian function fort=0,1,⋯,T−1, which is given by
(29)Ht=c1htyt+λt+1xt+sααΓαxt1−xt−axt2yt1+axt2+μt+1yt+sααΓαbxt2yt1+axt2−dyt2zt1+dyt2−m1yt−htyt+ζt+1zt+sααΓαcyt2zt1+dyt2−m2zt,where λt, μt, and ζt are the adjoint functions.Now, according to the maximum principle of Pontryagin, the necessary conditions of the above problem are(30)λt=∂Ht∂xt=λt+11+sααΓα1−2xt−2axtyt1+axt22+μt+1sααΓα2bxtyt1+axt22,μt=∂Ht∂yt=c1ht+λt+1sααΓα−axt21+axt2+μt+11+sααΓαbxt21+axt2−2dytzt1+dyt22−m1−ht+ζt+1sααΓα2cytzt1+dyt22,ζt=∂Ht∂zt=μt+1sααΓα−dyt21+dyt2+ζt+11+sααΓαcyt21+dyt2−m2.The characterization of the optimal solutionht∗ will be
(31)ht∗=0,hMax,ifSt<0,ifSt>0,where St is the switching function. The optimal control solutions and the corresponding state solutions at that t will be determined numerically.


## 4. Numerical Simulation

This section concerns with numerical simulations that confirm the theoretical results. For the fixed pointE1, we use the following set of value parameters: a=0.7, c=.01, b=0.1, d=0.2, m1=0.2, m2=0.2, h=0.01, s=0.8,α=0.98, x1=1.9, y1=0.006, z1=0.005, and T=200, so that according to Theorem 4, the fixed point E1 is the local stability point. Figure 1 shows the local stability of E1.Figure 1
The plot shows the local stability ofE1 of the system (8).Table 1
This table compares the total optimal harvesting with other total harvesting strategies.
The control variableTotal net harvestingJht=h∗J=0.0531ht=0.2J=0.0515ht=0.3J=0.0520ht=0.35J=0.0518ht=0.4J=0.0516For the local stability of the fixed pointE2 of system (8), the following set of values of the parameter are used: a=0.5, c=0.2, b=0.3, d=0.4, m1=0.1, m2=0.2, h=0.01, s=0.5, α=0.98, x1=1.2, y1=0.7, z1=0.22, and T=2000, so that according to Theorem 5, the fixed point E2 is the local stability point. Figure 2 shows that local stability.Figure 2
The plot shows the local stability ofE2 of the system (8).For the local stability of the unique positive fixed pointE3 of system (8), we use the following values of parameters: a=0.7, c=0.1, b=0.2, d=0.2, m1=0.02, m2=0.2, h=0.01, s=0.5, α=0.98, x1=0.9, y1=1.7, z1=0.1, and T=2000, so that according to Theorem 5, the fixed point E3 is the local stability point. Figure 3 shows the local stability of E3.Figure 3
The plot shows the local stability ofE3 of the system (8).For the control problems (27) and (28), we follow an iterative method which is found in [3, 12, 19, 30]. We choose the values of parameters as follows: a=0.7, b=0.2,c=0.1, c1=0.01, d=0.2, m1=0.01, m2=0.05, α=0.9, and T=60. The initial conditions are (0.9, 1.7, 0.15); with these values of parameters, we get the total optimal harvesting of Jopt=0.0531. In Table 1, we compare the total optimal harvesting and other total harvesting strategies using the same values of the parameters. Figure 4 shows the optimal control variable as a function of time and the effect of optimal harvesting on the prey species, middle predator, and top predator species.Figure 4
(a) The optimal control variable as a function of time. (b) The effect of optimal harvesting on the predator species. (c) The effect of optimal harvesting on the prey species. (d) The effect of optimal harvesting on the top predator species.
(a)(b)(c)(d)


## 5. Conclusions

In this work, we have studied and investigated the first-order discretization of three-dimensional fractional-order differential biological systems. The existence and the stability of its equilibrium points are discussed. We have observed that some of them exist without any restriction on the parameters of the model. Optimal harvesting has been done to put the economic optimality at a certain level to get maximum profit gain. The analytical results that have been obtained in this paper are verified through numerical various numerical simulations.



---

*Source: 8171037-2022-05-28.xml*