# Fractional Optimal Control with Fish Consumption to Prevent the Risk of Coronary Heart Disease

Authors: I. Ameen; M. Hidan; Z. Mostefaoui; H.M. Ali
Journal: Complexity
 (2020)

---

## Abstract

According to the World Health Organization (WHO), Chronic Heart Disease (CHD) is one of the greatest defies currently confronting humankind which is sweeping the whole globe, with an expanding trend in developing countries. In this paper, a mathematical model (MM) was proposed to study the connection between fish consumption and CHD mortality in Egypt, by considering a system of ordinary differential equations (ODEs) involving time-fractional derivative (FD). We considered here the study on Egypt for the ease of obtaining real data, but the method and approach adopted here is not limited to Egypt only and can be applied to any country in the world with the information of the real data related to the subject of the study. Additionally, the control function which represents the metabolic and the behavioural risk factors of CHD that help to reduce the number of mortality due to CHD is incorporated in the proposed MM. A fractional optimal control problem (FOCP) with a proposed control is formulated and studied theoretically using the Pontryagin maximum principle, to minimize the susceptible population and also to decrease the mortality rate of CHD. Moreover, firstly we discussed the positivity and boundedness of solutions; then, the model equilibria are determined and their local stability analysis was investigated; furthermore, we use the improved forward-backward sweep method (FBSM) based on the predictor-corrector method (PCM) in order to obtain the solution of proposed FOCP. In addition, some numerical simulations were performed to show the effect of the proposed optimal control (OC) besides the impact of fish consumption on the mortality of CHD.


---

## Body

## 1. Introduction

CHD is defined as lack of oxygen supply to the heart due to narrowing of the coronary arteries (so it is called coronary artery disease). Precisely, the CHD occurs as a result of deposition of fatty streaks “cholesterol” on the arterial wall which is called atherosclerosis. Sometimes rupture or erosion of fibrous cap of the coronary artery occurs resulting in stimulation of platelet adhesion and aggregation causing thrombus formation. This thrombus may cause partial obstruction of the arterial wall (unstable angina) or complete obstruction myocardial infarction (MI). CHD can lead to angina, which is symptomized by chest pain exaggerated by exertion and relieved by rest or nitrate, shortness of breath, sweating, nausea, and cramping. Generally, complications mainly due to MI are heart failure, myocardial rupture, rupture of free wall of left ventricle, arrhythmia, post-MI pericarditis, heart attack, and sudden death. In fact, there is no cure for CHD. In this case, treatment tends to involve making healthful lifestyle changes, such as quitting smoking and adopting a healthful diet. However, some people may need to take medications after undergoing medical procedures.CHD can be prevented by addressing risk factors, which is divided intouncontrollable risk: gender (in middle age, common in males than females due to estrogen protective effect, but after menopause, the percentage will be equal), age (incidence increases by age), congenital predisposition (positive family history), and controllable risk: metabolic risk factors such as hypertension, diabetes, obesity, and dyslipidemia (the most common forms of dyslipidemia involve high levels of low-density lipoproteins (LDL) or bad cholesterol and low levels of high-density lipoproteins (HDL) or good cholesterol), and behavioural risk factors such as smoking, alcohol, physical inactivity, and unhealthy diet (see for details [1]).CHD is the most dangerous cardiovascular disease (CVD), as it causes the most deaths of any heart disease in Egypt. Based on the data published in the WHO 2017, the mortality rate of CHD in Egypt reached 24.58% of total mortality. So, Egypt ranks 18th in the world for CHD deaths [2].It is known that eating fish has multiple benefits on our body as fish oils are rich in eicosapentaenoic acid (EPA) and docosahexaenoic acid (DHA) which are the main sources of polyunsaturated fatty acids (PUFAs) of the family omega-3 [3]. Omega-3 is an unsaturated fatty acid that has an anti-inflammatory effect, antithrombotic effect, and antiatherogenic effect. So, it has a protective effect on CHD.From the abovementioned view, the effect of fish consumption on cardiac diseases, several studies have been carried out, such as follows: in [4], Kris-Etherton et al. concluded that omega-3 acid supplements which are highly contained in fatty fish can reduce CVDs. In 2008, Virtanen et al. achieved that the consuming fish can reduce the risk of major chronic diseases (CVDs, cancer, etc) [5]. Raatz et al. in 2013 performed a statistical study on the USA population and found that CVD risk is decreased while eating fish frequently [6]. The effects of fish consumption through the relative risk measure on CHD is investigated in the paper of Hooper et al. [7]. Zheng et al., in [8], concluded that fish consumption has a significant protective effect on fatal CHD.In Egypt, the main fisheries are seas (Mediterranean Sea, 950 km (coastline length) [9] and Red Sea, 1500 km [9]), northern lakes (Mariout, 250 km2 (surface area) [10]; Edku, 70 km2 [11]; Burullus, 410 km2 [12]; and Manzala, 700 km2 [13], coastal lagoons (Port Fuad, 60 km2 [14] and Bardawil, 650 km2 [15]), and inland lakes (Bitter and Temsah (Suez Canal), 250 and 15 km2 [16]). With all that fisheries, there exist shortage in fish protein supply due to the high rate of local consumption of fish (for example, the average consumption in 2012 was 20.55 kg of fish products per capita [17]).Optimal control theory is another area of mathematics that arose after the Second World War with the formulation of the famous Pontryagin maximum principle [18], responding to practical needs of engineering, particularly in the field of aeronautics and flight dynamics [19]. During the last decades, MMs with OC have been largely used to epidemiological models and biomedicine, where it provides more insight into the dynamics of the diseases and provides preventive and appropriate strategies to combat these diseases and controlling the spread of infectious diseases. An example of some of the diseases that have been addressed using OC theory HIV/AIDS (see, e.g., [20–23]), HBV [24], Malaria and Cholera [25, 26], and many other diseases (see, e.g., [27–29]).Optimal control problems with fractional calculus are called fractional optimal control problems (FOCPs) and are considered the generalization of classical OC problems, in which the dynamics of the control system are described by fractional differential equations (FDEs), and the objective function may be given by a fractional integration operator [30]. Several research papers in the literature provided the theoretical basis and fundamentals of FOCPs. Moreover, these papers extensively studied how to formulate the FOCPs and derived the optimality conditions for several states and control variables using analytical and numerical methods (see, e.g., [31–41]). Recently, the FOCPs have been applied to epidemiological models for faster and more accurate behaviour to controlling of diseases, as the fractional-order depends on the memory. So, the FOCPs can be potential flexible tools for modeling epidemiological and biological systems related to memory. Ding et al. [42] applied fractional OC on the HIV-Immune system model and solve this problem by a forward-backward algorithm. Basir et al. [43]. presented a fractional OC of an enzyme kinetic model and solved it numerically. Kheiri et al. [44] formulated and discussed a fractional model for HIV/AIDS with treatment and included three control efforts (effective use of condoms, ART treatment, and behavioural change control) into the model aimed at controlling the spread of HIV/AIDS epidemic. Sweilam et al. [45] proposed OC of a fractional novel West Nile virus model and used two simple numerical methods to solve this problem. Ali and Ameen [46] formulated and studied the FOCP with three suggested controls to describe the transmission dynamics of pine wilt disease and controlling the spread of this disease.The main aim of this paper is to find the best strategy to significantly reduce CHD mortality by using the fractional OC technique with consumption of an appropriate amount of fish. This is carried out by formulating a fractional mathematical model (FMM) involving a control function, which depends on the relationship between fish consumption from the fish population living along the Egyptian coasts and the dynamics of a population at risk of CHD. To simplify the computational complexities, we used developed FBSM based on PCM to obtain the numerical solution of state and costate equations for suggested FOCP with left Caputo fractional derivative (CFD). Consequently, we have two trends to reduce the rate of individuals at risk of CHD: increasing fish consumption (depend on total production, imports, and exports of fishery products) and controlling metabolic and behavioural risk factors.The organization of this paper as follows. In Section2, we present a formulation of the FMM. In Section 3, we introduce some mathematical preliminaries of the fractional calculus which are needed to demonstrate our main results. Moreover, positivity and boundedness of solutions for the proposed model are presented. In Section 4, we discussed the equilibrium points (EPs) and stability analysis of the model. The formulation of the FOCP and the derivation of the necessary optimality conditions are given in Section 5. In Section 6, the numerical method is introduced for solving the proposed model with and without control. In Section 7, we present numerical simulations of our results to devise the best strategy for the reduction of the disease. In Section 8, the conclusions are given.


## 2. Model Formulation

The first available MM was proposed by Lamlili et al. [47] to study the impact of fish consumption on CHD mortality in Morocco. Extension of the results in reference [47] can be found in [48], where the authors presented an OC strategy which can be obtained by acting on the controllable risk factors. The MM proposed in [48] is based on ODEs. Here, we introduced a more generalized model that is governed by a system of ODEs involving time FD; moreover, the parameters in this model are modified parameters of fractional-order (for more details, see [49, 50]). The modified model is considered as the following:(1)D0CtαSt=Λα−γα+ηdα1−uSt+βαIt,D0CtαIt=−γα+δαRR+βαIt+ηdα1−uSt,D0CtαXt=bα1−XtKXt−cαEXt,subject to(2)St0=N1>0,It0=N2>0,Xt0=N3>0,where we assume the population of size N constituted of St, the number of persons without CHD risk, and It, the number of persons with CHD risk, which is given by the equation Pt≔Nt=St+It at any time t∈0,tf. Xt denoted the biomass of the fish population at any time t. In [47], the logistic function, which describes the likelihood of developing CHD in a group of individuals eating fish compared to individuals who have little or no consumption, is given by RRt=a1/a2+ea3μXt with μ=cαE/Ptaϒ (see Figure 1). We use the variable u to represent how to control risk factors (i.e., controllable risk). So, the factor 1−u is used to minimize the susceptible population and also to decrease the mortality rate of CHD and d∈0,1 be a real positive constant to measure the efficacy of u. The meaning of parameters for the fractional-order model (FOM) (1) is given in Table 1.Figure 1
The relative riskRR with a=1 for different values of α.Table 1
Meaning and values of the parameters in system (1).
Parameter
Description
Value
Ref.
Λ
α
The recruitment of persons without CHD risk
0.0268
α
[51]
γ
α
The natural mortality rate
0.0057
α
[51]
η5
The probability to have CHD
0.06
Assumed
β
α
The rate of patients with CHD who are cured
0.005
α
Assumed
δ
α
The mortality rate due to CHD
0.002
α
[51]
b
α
The biotic potential
1
α
[17]
K
The carrying capacity
112
×
0.4
×
10
9
[17]
c
α
The catchability coefficient
0.04
α
[17]
E
The total harvesting effort
9
[17]
a
1
,
a
2 and a3
The negative coefficients of nonlinear regression
a
1
=
−
0.73
,
a
2
=
−
1.71
,
[47]
a
3
=
−
0.0616
a
The rate of fish consumed from the total harvested
Varies
Assumed
ϒ
The coefficient that transforms the consumption of fish per capita and per year to a frequency per month
559.5
[52]


## 3. Preliminaries

Firstly, we introduce some basic concepts of fractional calculus which help us complete this research (see e.g., [53]). Then, in Section 3.1, we state the main theorem about the nonnegative solutions of the FOM (1).Definition 1.
For a given functionf:ℝ+⟶ℝ and α>0, the left and right Riemann–Liouville fractional integrals (RLFIs) are, respectively, defined by(3)Iatαft=1Γα∫att−τα−1fτdτ,Itbαft=1Γα∫tbτ−tα−1fτdτ,where Γ⋅ is the Euler gamma function defined by(4)Γz=∫0∞e−ttz−1dt,Rez>0.
Letf⋅ be absolutely continuous functions on a,b and n−1<α≤n, where n∈ℕ. Then, the left and right for Riemann–Liouville fractional derivatives (RLFDs) and CFDs are defined as follows.Definition 2.
The left and right RLFDs are, respectively,(5)Datαft=1Γn−αddtn∫att−τn−α−1fτdτ,Dtbαft=1Γn−α−ddtn∫tbτ−tn−α−1fτdτ.Definition 3.
The left and right CFDs are, respectively,(6)DaCtαft=1Γn−α∫att−τn−α−1fnτdτ,DtCbαft=−1nΓn−α∫tbτ−tn−α−1fnτdτ.
There is a relation between left and right for RLFDs and CFDs as described in the following theorem.Theorem 1.
Lett>0 and n−1<α≤n∈ℕ. Then, the following relation holds:(7)Datαft=DaCtαft+∑j=0n−1fjaΓj−α+1t−aj−α,Dtbαft=DtCbαft+∑j=0n−1fjbΓj−α+1b−tj−α.
Therefore,(8)iffa=f′a=⋯=fn−1a=0,thenDatαft=DaCtαft,iffb=f′b=⋯=fn−1b=0,thenDtbαft=DtCbαft.Now, we recall the definitions of Laplace transform of Caputo’s derivative and Mittag-Leffler function in two arguments.Definition 4.
LetFs be the Laplace transform of the function ft. Then, Laplace transform of the Caputo derivative is given by(9)ℒDαCft,s=sαFs−∑i=0n−1sα−i−1fi0,n−1<α≤n,n∈ℕ.Definition 5.
Forx∈ℝ, the Mittag-Leffler function El,mx is defined by(10)El,mx=∑n=0∞xnΓln+m,l>0,m>0.
Then, the Laplace transform of the functiontm−1El,m±λtl is defined as follows:(11)ℒtm−1El,m±λtl=sl−msl∓λ.


### 3.1. Positivity and Boundedness of Solutions

LetΩ=S,I,X∈ℝ+30≤S+I≤Λ/γα,0≤X≤K. Then, we can show that the biological feasible region Ω is positivity invariant and bounded by the following theorem.Theorem 2.
There is a unique solution for the FOM (1) at t≥0 (where t0=0) and it remains in the set Ω.Proof.
We know that from Theorem1 and Remark 3.2 in [54], the solution on 0,∞ is existent and unique. Now, it remains to prove that Ω is a positivity invariant. We can write the fractional-order differential equation representing the total population as follows:(12)D0CtαNt+γαNt=gt,where gt=Λα−δαRRIt is a constant function of time. Solving equation (12) using the Laplace transform method [53], then we have the following solution:(13)Nt=∫0tt−τα−1Eα,α−γαt−ταgτdτ≥0,where 0<α<1,γ>0, and Eα,α−γαt−τα is the two-parameter Mittag-Leffler function (see Definition 5). Since Mittag-Leffler function is an entire function [53], thus Eα,α−γαt−τα is bounded for all t>0. Therefore, as t⟶∞, we have N≤Λ/γα. For the last equation in the fractional system (1), we apply Laplace-Adomian decomposition method (LADM) (see, e.g., [55]) as follows:(14)sαℒX−sα−1X0=−bαKℒX2+bα−cαEℒX,ℒX=X0s−bαKsαℒX2+bα−cαEsαℒX.
The method assumes the solution as an infinite series, and the nonlinearityX2 is decomposed as follows:(15)Xt=∑n=0∞Xkt,X2t=∑n=0∞Ant,where An is the so-called Adomian Polynomial given as(16)An=1n!dndλn∑j=0nλjXj2λ=0.
Substituting from equations (15) and (16) into (14), the result is(17)ℒX0=X0s,(18)ℒX1=−bαKsαℒA0+bα−cαEsαℒX0,…ℒXn+1=−bαKsαℒAn+bα−cαEsαℒXn.
Applying the inverse fractional Laplace transform to equation (17), hence the value of X0 can be obtained. Substituting the value of X0 in equation (16), the first Adomian polynomial A0 is obtained, and then substituting X0 and A0 in equation (18) to obtain X1 and proceeding in a similar way, the other terms X2,X3,… can be computed recursively and we can write the solution as Xt=X0+X1+X2+…. This show that if Xt>0, then for the term in the series of solution (i.e., X1) to be a positive value; hence, it is necessary to have K>A0 this implies that Xt≤K for all t∈0,tf. Obviously, the subinterval for the biomass X (i.e., Xt≤K) is adaptive with the ecological assumption that the biomass X cannot be greater than the carrying capacity K.


## 3.1. Positivity and Boundedness of Solutions

LetΩ=S,I,X∈ℝ+30≤S+I≤Λ/γα,0≤X≤K. Then, we can show that the biological feasible region Ω is positivity invariant and bounded by the following theorem.Theorem 2.
There is a unique solution for the FOM (1) at t≥0 (where t0=0) and it remains in the set Ω.Proof.
We know that from Theorem1 and Remark 3.2 in [54], the solution on 0,∞ is existent and unique. Now, it remains to prove that Ω is a positivity invariant. We can write the fractional-order differential equation representing the total population as follows:(12)D0CtαNt+γαNt=gt,where gt=Λα−δαRRIt is a constant function of time. Solving equation (12) using the Laplace transform method [53], then we have the following solution:(13)Nt=∫0tt−τα−1Eα,α−γαt−ταgτdτ≥0,where 0<α<1,γ>0, and Eα,α−γαt−τα is the two-parameter Mittag-Leffler function (see Definition 5). Since Mittag-Leffler function is an entire function [53], thus Eα,α−γαt−τα is bounded for all t>0. Therefore, as t⟶∞, we have N≤Λ/γα. For the last equation in the fractional system (1), we apply Laplace-Adomian decomposition method (LADM) (see, e.g., [55]) as follows:(14)sαℒX−sα−1X0=−bαKℒX2+bα−cαEℒX,ℒX=X0s−bαKsαℒX2+bα−cαEsαℒX.
The method assumes the solution as an infinite series, and the nonlinearityX2 is decomposed as follows:(15)Xt=∑n=0∞Xkt,X2t=∑n=0∞Ant,where An is the so-called Adomian Polynomial given as(16)An=1n!dndλn∑j=0nλjXj2λ=0.
Substituting from equations (15) and (16) into (14), the result is(17)ℒX0=X0s,(18)ℒX1=−bαKsαℒA0+bα−cαEsαℒX0,…ℒXn+1=−bαKsαℒAn+bα−cαEsαℒXn.
Applying the inverse fractional Laplace transform to equation (17), hence the value of X0 can be obtained. Substituting the value of X0 in equation (16), the first Adomian polynomial A0 is obtained, and then substituting X0 and A0 in equation (18) to obtain X1 and proceeding in a similar way, the other terms X2,X3,… can be computed recursively and we can write the solution as Xt=X0+X1+X2+…. This show that if Xt>0, then for the term in the series of solution (i.e., X1) to be a positive value; hence, it is necessary to have K>A0 this implies that Xt≤K for all t∈0,tf. Obviously, the subinterval for the biomass X (i.e., Xt≤K) is adaptive with the ecological assumption that the biomass X cannot be greater than the carrying capacity K.


## 4. Equilibrium Points and Stability Analysis

For the FOM (1), EPs are defined such that the left CFDs of I,S, and X is equal to zero, namely,(19)D0CtαS=0,D0CtαI=0,D0CtαX=0,then our model (1) has two EPs.


### 4.1. First Equilibrium Point (EP) atXeq=X0=0

In the case ofX0=0, means that there is no biomass of fish population which implies that the mortality rate due to CHD is not affected by fish consumption:(20)E0=S0,I0,X0,where(21)S0=Λαβα+γα+δαΘ1βαγα+γ2α+γαηdα1−u+δαγα+ηdα1−uΘ1,I0=ηdα1−uΛαβαγα+γ2α+γαηdα1−u+δαγα+ηdα1−uΘ1,Θ1=RRtX0=0=a1a2+1.


### 4.2. Second EP atXeq=X∗≠0

We now consider the case where there is biomass of fish population; thus,(22)X∗=Kbαbα−cαE.By substituting in other equations in (1), then we have(23)E∗=S∗,I∗,X∗,where(24)S∗=Λαβα+γα+δαΘ2βαγα+γ2α+γαηdα1−u+δαγα+ηdα1−uΘ2,I∗=ηdα1−uΛαβαγα+γ2α+γαηdα1−u+δαγα+ηdα1−uΘ2,Θ2=RRtX0=X∗=a1a2+ea3μK/bαbα−cαE.It is clear that from the existence and positivity of equilibria forXt, the fishing effort (E) should be less than the Biotechnical Productivity b/cα. In other words, we exclude the case of E>b/cα because this mean that a problem of overfishing will take place, and also the exponential function in Θ2 will tend to infinity, and then finally this choice will contribute to the extinction of the fish population.We conclude the above results for the EPs as follows.Theorem 3.
The EPE0 of system (1) always exists without any constraints, but the EP E∗ represented in (22) and (24) exists if E<b/cα.In general, the stability region of the system contains FDs, which is greatest when compared with the stability region of the system with integer-order derivatives (see, e.g., [56]). Therefore, we investigate analytically the stability of the EPs E0 and E∗.Theorem 4.
The EPE0 of the FOM (1) is locally unstable.Proof.
Determining the Jacobian matrix of system (1) at E0, we have(25)JE0=−γα−ηdα1−uβα0ηdα1−u−γα−βα−δαΘ1δαa3μΘ12a1I000bα−cαE.
Obviously, one of the eigenvalues ofJE0 is λ3=bα−cαE>0 (from Theorem 3), then E0 will be linearly (locally) unstable (see, e.g., [57]).
For the second EPE∗, we have the following theorem.Theorem 5.
The EPE∗ of system (1) is locally stable.Proof.
The Jacobian matrix evaluated atE∗ gives(26)JE∗=−γα−ηdα1−uβα0ηdα1−u−γα−βα−δαΘ2δαa3μΘ22ea3μX∗a1I∗00−bα+cαE.
Then, we obtain the characteristic polynomial ofJE∗, by using this equation detJE∗−λI=0, where χj,j=1,2,3, are the eigenvalues of JE∗. Therefore, one of the eigenvalues is χ1=−bα+cαE, which is negative according to Theorem 3, and the last two eigenvalues are given by the quadratic equation(27)χ2+F1χ+F2=0,where(28)F1=2γα+ηdα1−u+βα+δαΘ2,F2=γ2α+γαηdα1−u+βα+δαΘ2ηdα1−u+γα.
Thus, the coefficientsF1 and F2 of polynomial (27) are positive (i.e., all coefficients of polynomial (27) have the same signal); then, the roots (eigenvalues) have negative real part (Routh’s criterion, see, e.g., [58]). We can conclude that the steady state E∗ of system (1) is locally stable.


## 4.1. First Equilibrium Point (EP) atXeq=X0=0

In the case ofX0=0, means that there is no biomass of fish population which implies that the mortality rate due to CHD is not affected by fish consumption:(20)E0=S0,I0,X0,where(21)S0=Λαβα+γα+δαΘ1βαγα+γ2α+γαηdα1−u+δαγα+ηdα1−uΘ1,I0=ηdα1−uΛαβαγα+γ2α+γαηdα1−u+δαγα+ηdα1−uΘ1,Θ1=RRtX0=0=a1a2+1.


## 4.2. Second EP atXeq=X∗≠0

We now consider the case where there is biomass of fish population; thus,(22)X∗=Kbαbα−cαE.By substituting in other equations in (1), then we have(23)E∗=S∗,I∗,X∗,where(24)S∗=Λαβα+γα+δαΘ2βαγα+γ2α+γαηdα1−u+δαγα+ηdα1−uΘ2,I∗=ηdα1−uΛαβαγα+γ2α+γαηdα1−u+δαγα+ηdα1−uΘ2,Θ2=RRtX0=X∗=a1a2+ea3μK/bαbα−cαE.It is clear that from the existence and positivity of equilibria forXt, the fishing effort (E) should be less than the Biotechnical Productivity b/cα. In other words, we exclude the case of E>b/cα because this mean that a problem of overfishing will take place, and also the exponential function in Θ2 will tend to infinity, and then finally this choice will contribute to the extinction of the fish population.We conclude the above results for the EPs as follows.Theorem 3.
The EPE0 of system (1) always exists without any constraints, but the EP E∗ represented in (22) and (24) exists if E<b/cα.In general, the stability region of the system contains FDs, which is greatest when compared with the stability region of the system with integer-order derivatives (see, e.g., [56]). Therefore, we investigate analytically the stability of the EPs E0 and E∗.Theorem 4.
The EPE0 of the FOM (1) is locally unstable.Proof.
Determining the Jacobian matrix of system (1) at E0, we have(25)JE0=−γα−ηdα1−uβα0ηdα1−u−γα−βα−δαΘ1δαa3μΘ12a1I000bα−cαE.
Obviously, one of the eigenvalues ofJE0 is λ3=bα−cαE>0 (from Theorem 3), then E0 will be linearly (locally) unstable (see, e.g., [57]).
For the second EPE∗, we have the following theorem.Theorem 5.
The EPE∗ of system (1) is locally stable.Proof.
The Jacobian matrix evaluated atE∗ gives(26)JE∗=−γα−ηdα1−uβα0ηdα1−u−γα−βα−δαΘ2δαa3μΘ22ea3μX∗a1I∗00−bα+cαE.
Then, we obtain the characteristic polynomial ofJE∗, by using this equation detJE∗−λI=0, where χj,j=1,2,3, are the eigenvalues of JE∗. Therefore, one of the eigenvalues is χ1=−bα+cαE, which is negative according to Theorem 3, and the last two eigenvalues are given by the quadratic equation(27)χ2+F1χ+F2=0,where(28)F1=2γα+ηdα1−u+βα+δαΘ2,F2=γ2α+γαηdα1−u+βα+δαΘ2ηdα1−u+γα.
Thus, the coefficientsF1 and F2 of polynomial (27) are positive (i.e., all coefficients of polynomial (27) have the same signal); then, the roots (eigenvalues) have negative real part (Routh’s criterion, see, e.g., [58]). We can conclude that the steady state E∗ of system (1) is locally stable.


## 5. Formulation of FOCP and Derivation of Optimality Conditions

In order to get necessary optimality conditions with the aim of reducing the rate of individuals at risk of CHD, we suggest FOCP as the following form:(29)PMinimizeWtfsubject toD0CtαSt=Λα−γα+ηdα1−utSt+βαIt,D0CtαIt=−γα+δαRR+βαIt+ηdα1−utSt,D0CtαXt=bα1−XtKXt−cαEXt,D0CtαWt=It+Bu2t,with convenient initial conditions(30)St0=N1>0,It0=N2>0,Xt0=N3>0,Wt0=0,where Wtf is the objective function given by Wtf=I0tfαIt+Bu2t, I0tfα is the left RLFI, D0Ctα is the left CFD of order 0<α≤1, and B represents the proportional weight constant related to control ut.We consider the setU is the set of admissible control functions and it is Lebesgue measurable, which is defined as(31)U=ut:0≤u≤umax≤1,t∈0,tf.Furthermore, we summarize the necessary optimality condition of the FOCPP as follows.Theorem 6.
IfS∗,I∗,X∗,W∗,u∗ be OC process which minimizes the objective function in the problem P, then there exist adjoint variables pi,i=1,2,3,4, satisfying(32)Dttfαp1=p2−p1ηdα1−u∗−p1γα,Dttfαp2=p1βα−p2γα+δRR∗+βα+p4,Dttfαp3=p3bα−2bX∗K−cαE+p2δαa2μea3μX∗a2+ea3μX∗RR∗I∗,Dttfαp4=0,with transversality conditions(33)pjtf=0,j=1,2,3,p4tf=1,where the operator Dttfα is right RLFD.
Moreover, we obtain the OCu∗ which minimizes the problem P over the region U as follows:(34)u∗=minmax0,u¯,1,where(35)u¯=p2−p1ηdαS∗2Bp4.Proof.
According to the Pontryagin maximum principle [18] and results announced in [30], we define the Hamiltonian H of our FOCP as the following:(36)HP⟶,G⟶,u,t=P⟶Tt.G⟶ς⟶,u,t,with(37)P⟶=p1p2p3p4,ς⟶=St,It,Xt,Wt,G⟶ς⟶,u,t=g1ς⟶,u,tg2ς⟶,u,tg3ς⟶,u,tg4ς⟶,u,t,where(38)g1ς⟶,u,t=Λα−γα+ηdα1−uSt+βαIt,g2ς⟶,u,t=−γα+δαRR+βαIt+ηdα1−uSt,ag3ς⟶,u,t=bα1−XtKXt−cαEXt,g4ς⟶,u,t=It+Bu2.
Let(39)G∗⟶ς∗⟶,u∗,t=g1ς∗⟶,u∗,tg2ς∗⟶,u∗,tg3ς∗⟶,u∗,tg4ς∗⟶,u∗,t,ς∗⟶=S∗t,I∗t,X∗t,W∗t.
Ifu∗∈U is optimal for problem (P) with nonnegative initial conditions (30), then there exists a nontrivial absolutely continuous mapping p:0,tf⟶ℝ4, pt=p1t,p2t,p3t,p4t, called the adjoint vector, such that(40)Dttfαp1=∂HP⟶,G∗⟶,u∗,t∂S∗=−p1γα+ηdα1−u∗+p2ηdα1−u∗,Dttfαp2=∂HP⟶,G∗⟶,u∗,t∂I∗=p1βα−p2γα+δRR∗+βα+p4,Dttfαp3=∂HP⟶,G∗⟶,u∗,t∂X∗=p3bα−2bX∗K−cαE+p2δαa1a2μea3μX∗a2+ea3μX∗2I∗,Dttfαp4=∂HP⟶,G∗⟶,u∗,t∂W∗=0,where H is the Hamiltonian defined in equation (36). Since the state variables Stf,Itf,andXtf are free at the terminal time tf, while the Wtf appear in the objective function, then the transversality conditions are given by (41)pjtf=0,j=1,2,3,p4tf=1.
For optimality conditions, we differentiate the Hamiltonian with respect to the control variablesu∗ in the interior of the set U. Then, we have(42)∂HP⟶,G∗⟶,u∗,t∂u∗=p1−p2ηdαS∗+2p4Bu∗=0.
Therefore,(43)u∗=p2−p1ηdαS∗2p4B,and impose the bounds 0≤u∗≤umax∗≤1 on the control to yield equation (34) as required.
By using equation (7), the system of adjoint variables (32) is equivalent to the following right CFD:(44)DtCtfαp1=p2−p1ηdα1−u∗−p1γα,DtCtfαp2=p1βα−p2γα+δRR∗+βα+p4,DtCtfαp3=p3bα−2bX∗K−cαE+p2δαa2μea3μX∗a2+ea3μX∗RR∗I∗,DtCtfαp4=−tf−t−αΓ1−α,with the same transversality conditions (33).
To rewrite the adjoint variables (44) as a system of equations with left CFD we use the following Lemma.Lemma 1.
For0<α≤1, we have(45)DtCtfαpt=D0Ctαptf−t.Proof.
From Definition3, we have(46)DtCtfαpt=−1Γ1−α∫ttfτ−t−αp′τdτ,by replacing t to tf−t, we have(47)Dtf−tCtfαptf−t=−1Γ1−α∫tf−ttfτ−tf−t−αp′τdτ.
Letξ=tf−τ, then(48)Dtf−tCtfαptf−t=−1Γ1−α∫t0t−ξ−αp′tf−ξ−dξ=1Γ1−α∫t0t−ξ−αptf−ξ′dξ=0CDtαptf−t,by replacing tf−t to t in the left-hand side, the lemma will be proved.
Therefore, we can rewrite the adjoint variables (44) as a system of equations with the left CFD as follows:(49)D0Ctαp1tf−t=p2tf−t−p1tf−tηdα1−u∗tf−t−p1tf−tγα,D0Ctαp2tf−t=p1tf−tβα−p2tf−tγα+δRR∗tf−t+βα+p4tf−t,D0Ctαp3tf−t=p3tf−tbα−2bX∗tf−tK−cαE+p2tf−tδαa2μea3μX∗tf−ta2+ea3μX∗tf−tRR∗tf−tI∗tf−t,D0Ctαp4tf−t=−t−αΓ1−α.
Then, the state and costate (adjoint variables) equations, which are presented in equations (29) and (49), are given in terms of the left CFD. This approach simplifies the use of fractional numerical methods to solve the state and costate equations. Therefore, to solve a FOCP, we can apply only one fractional numerical method in Forward-Backward algorithms.


## 6. Description of the Numerical Methods

This section deals with the numerical methods used in this paper without restrictions of the time interval. In the beginning, we intend to solve model (1) in the absence of control u by the following discrete method, which offers accurate solutions during a long time interval.


### 6.1. The Implementation of PCM

The implicit fractional Adams method of order 2 is a generalization of the classical trapezoidal rule (for more details, see [59–61]). For instance, consider the initial value problem:(50)D0Ctαyt=ft,yt,yt0=y0,0<α≤1,t0<t≤T.In a discrete numerical method, the time intervalt0,T is replaced by a discrete set of points tj=t0+jh, h=T−t0/N, and j=0,1,…,N so that the solution is approximated by a sequence yjj=0,1,…,N such that yj≈ytj. Thus, the general formula for the predictor-corrector implementation with use of corrector iterations is described as follows:(51)yk+1i=y0+hαΓα+2∑j=0kaj,k+1ftj,yj+ftk+1,yk+1i−1,(52)yk+10≔yk+1pre=y0+hαΓα+1∑j=0kbj,k+1ftj,yj,i=1,…,q,where(53)aj,k+1=kα+1−k−αk+1α,ifj=0,k−j+2α+1+k−jα+1−2k−j+1α+1,if1≤j≤k,1,ifj=k+1,(54)bj,k+1=k+1−jα−k−jα,where yk+1i denotes the approximation after i corrector steps, yk+10=yk+1pre is the predictor, and yk+1q is the final approximation after q corrector steps that we actually use. As a consequence of the assumptions of Theorem 4.2.6 in [62], the error can be estimated as follows.Theorem 7 [63].
The approximation computed by schemes (51) and (52) described above satisfies(55)max0≤j≤Nytj−yj=OhM,where M=min2,1+qα.We now present some numerical experiments for model (1) (in the absence of u), using different values for a, by applying formulas (51) and (52) with this value of initial conditions Pt=98.5×106,S0=0.72×Pt,I0=0.28×Pt,andX0=0.36×106 and the values of parameters in Table 1.In Figure2, we consider different scenarios depending on the parameter a in order to achieve one of the trends to reduce the rate of individuals at risk of CHD.Figure 2
Mortality due to CHD, wherea=0,0.3,1 with different values of α (in the absence of u).


### 6.2. FBSM Based on PCM

FBSM is an effective iterative approach to solve the optimality systems, and it has been developed using the RK4 scheme to solve the classic OC problems. Based on the predictor-corrector implementation described by (51) and (52), we improve FBSM to solve our FOCP which is illustrated in the following algorithm.Step 0: putting the initial conditions and the value of parameters.Step 1: divide the time interval intoN subintervals of equal length and set h=tf−t0/N, tj=t0+jh, and j=0,1,…,N. This means that we begin the main for loop in the algorithm.Step 2: computeuj according to the following formula:(56)uj=minmax0,p2j−p1jηdαSj2Bp4j,1,j=0,1,…,N,where the initial value of the control u is determined when j=0 with the above initial conditions (30) and transversality conditions (33). The remaining values of u (i.e., for j=1,2,…,N−1) can be obtained from the looping to solve the state and costate equations with the initial and transversality conditions as in the next steps.Step 3: solving state equation (29) by applying formulas (51) and (52) with the initial value of u and the same starting point S0,I0,X0,W0 to obtain the new starting point Sj,Ij,Xj,Wj.Step 4: solve costate equation (49) by formulas (51) and (52) subject to the terminal conditions, values of u, and vector filed ς as follows:(57)PN−k−1i=ΨN−k−1+hαΓα+2∂HtN−k−1,ςN−k−1,uN−k−1,PN−k−1i−1∂ς,ΨN−k−1=hαΓα+2∑j=0kaj,k+1∂HtN−j,ςN−j,uN−j,PN−j∂ς,PN−k−10≔PN−k−1pre=hαΓα+1∑j=0kbj,k+1∂HtN−j,ςN−j,uN−j,PN−j∂ς,i=1,…,q,where the coefficients of aj,k+1 and bj,k+1 are given by equations (53) and (54), respectively.Step 5: updating the value of the controlu for j=1,2,…,N by entering the new ς and P values into Step 2.Step 6: the stopping criterion is a convergence, that is, it means if the values of the variables in a current iteration and the last iteration (i=1,2,…,q) are negligibly close, this implies that the outputs of the algorithm are our solutions. If the values are not close, return to Step 3.


## 6.1. The Implementation of PCM

The implicit fractional Adams method of order 2 is a generalization of the classical trapezoidal rule (for more details, see [59–61]). For instance, consider the initial value problem:(50)D0Ctαyt=ft,yt,yt0=y0,0<α≤1,t0<t≤T.In a discrete numerical method, the time intervalt0,T is replaced by a discrete set of points tj=t0+jh, h=T−t0/N, and j=0,1,…,N so that the solution is approximated by a sequence yjj=0,1,…,N such that yj≈ytj. Thus, the general formula for the predictor-corrector implementation with use of corrector iterations is described as follows:(51)yk+1i=y0+hαΓα+2∑j=0kaj,k+1ftj,yj+ftk+1,yk+1i−1,(52)yk+10≔yk+1pre=y0+hαΓα+1∑j=0kbj,k+1ftj,yj,i=1,…,q,where(53)aj,k+1=kα+1−k−αk+1α,ifj=0,k−j+2α+1+k−jα+1−2k−j+1α+1,if1≤j≤k,1,ifj=k+1,(54)bj,k+1=k+1−jα−k−jα,where yk+1i denotes the approximation after i corrector steps, yk+10=yk+1pre is the predictor, and yk+1q is the final approximation after q corrector steps that we actually use. As a consequence of the assumptions of Theorem 4.2.6 in [62], the error can be estimated as follows.Theorem 7 [63].
The approximation computed by schemes (51) and (52) described above satisfies(55)max0≤j≤Nytj−yj=OhM,where M=min2,1+qα.We now present some numerical experiments for model (1) (in the absence of u), using different values for a, by applying formulas (51) and (52) with this value of initial conditions Pt=98.5×106,S0=0.72×Pt,I0=0.28×Pt,andX0=0.36×106 and the values of parameters in Table 1.In Figure2, we consider different scenarios depending on the parameter a in order to achieve one of the trends to reduce the rate of individuals at risk of CHD.Figure 2
Mortality due to CHD, wherea=0,0.3,1 with different values of α (in the absence of u).


## 6.2. FBSM Based on PCM

FBSM is an effective iterative approach to solve the optimality systems, and it has been developed using the RK4 scheme to solve the classic OC problems. Based on the predictor-corrector implementation described by (51) and (52), we improve FBSM to solve our FOCP which is illustrated in the following algorithm.Step 0: putting the initial conditions and the value of parameters.Step 1: divide the time interval intoN subintervals of equal length and set h=tf−t0/N, tj=t0+jh, and j=0,1,…,N. This means that we begin the main for loop in the algorithm.Step 2: computeuj according to the following formula:(56)uj=minmax0,p2j−p1jηdαSj2Bp4j,1,j=0,1,…,N,where the initial value of the control u is determined when j=0 with the above initial conditions (30) and transversality conditions (33). The remaining values of u (i.e., for j=1,2,…,N−1) can be obtained from the looping to solve the state and costate equations with the initial and transversality conditions as in the next steps.Step 3: solving state equation (29) by applying formulas (51) and (52) with the initial value of u and the same starting point S0,I0,X0,W0 to obtain the new starting point Sj,Ij,Xj,Wj.Step 4: solve costate equation (49) by formulas (51) and (52) subject to the terminal conditions, values of u, and vector filed ς as follows:(57)PN−k−1i=ΨN−k−1+hαΓα+2∂HtN−k−1,ςN−k−1,uN−k−1,PN−k−1i−1∂ς,ΨN−k−1=hαΓα+2∑j=0kaj,k+1∂HtN−j,ςN−j,uN−j,PN−j∂ς,PN−k−10≔PN−k−1pre=hαΓα+1∑j=0kbj,k+1∂HtN−j,ςN−j,uN−j,PN−j∂ς,i=1,…,q,where the coefficients of aj,k+1 and bj,k+1 are given by equations (53) and (54), respectively.Step 5: updating the value of the controlu for j=1,2,…,N by entering the new ς and P values into Step 2.Step 6: the stopping criterion is a convergence, that is, it means if the values of the variables in a current iteration and the last iteration (i=1,2,…,q) are negligibly close, this implies that the outputs of the algorithm are our solutions. If the values are not close, return to Step 3.


## 7. Simulation of the Fractional Optimal Control Model

The main purpose of this section is to discuss the numerical solutions of the FOCPP and the effect a proportion of the fish consumed from the total harvested on this problem in the classical (α = 1) and fractional (different values of α) cases. The numerical results are achieved using MATLAB with initial conditions, realistic hypothetical parameter values in Table 1, and the weight constant value B=3,500,000.In the sentences below, we constructed some strategies depending on the proportion of fish consumed from the total harvested (i.e., parametera) to decrease the number of total mortality due to CHD and compared it in two cases, the presence and absence of the control operator.Strategy 1 (a=0).
In this strategy, we assume there is no fish consumption at all and only consider the effect of the control as in Figure3. We observe that there is a significant decrease in mortality due to CHD when we use the OC other than no control. Moreover, we take different values of α to present the effect of the FD on the number of total mortality due to CHD.Figure 3
Mortality due to CHD, wherea=0 without control (black lines) and with control (blue lines) for α=1 (solid line), α=0.98 (dashed line), and α=0.95 (dotted line).Strategy 2 (a=0.3).
In this strategy, we consider the fish consumption around 30% of the total harvested fish, and this strategy is illustrated in Figure4. From this simulation, we find that there is a significant difference in total mortality due to CHD when compared with and without control. Furthermore, the effects of fractional (different values of α) have been studied in this case. Also, we find this strategy decreases much when compared to Strategy 1.Figure 4
Mortality due to CHD, wherea=0.3 without control (black lines) and with control (blue lines) for α=1 (solid line), α=0.98 (dashed line), and α=0.95 (dotted line).Strategy 3 (a=0.6).
In this strategy, we consider the fish consumption around 60% of the total harvested fish, and the simulation results of this strategy are shown in Figure5. This simulation shows that there is a big difference in total mortality due to CHD when compared with and without control. Also, the effects of fractional (different values of α) have been investigated in this case. We observe that these decrease much when compared to the other two cases.Figure 5
Mortality due to CHD, wherea=0.6 without control (black lines) and with control (blue lines) for α=1 (solid line), α=0.98 (dashed line), and α=0.95 (dotted line).Strategy 4 (a=1).
In this strategy, we consider the total harvested fish completely consumed, and the effect of this strategy can be seen in Figure6. We note that there is a big difference in total mortality due to CHD when compared with and without control. Also, the effects of fractional (different values of α) have been studied in this case. We observe that these decrease much when compared to the other three cases.
We conclude from these simulations that whenever the consumption of fish from the total harvested fish increased with avoiding the metabolic and behavioural risk factors of CHD, then the CHD mortality rate is reduced, and this is evident in Figure7. For fractional cases, when α is reduced from 1, the memory effect of the problem increases, and therefore the behaviour of mortality due to CHD grows slowly. Moreover, the impact of memory for FOM can be seen in Figures 4–6, where increasing fish consumption reduces the memory of CHD, and this implies that there is a disparity in the curves for fractional-order cases.Figure 6
Mortality due to CHD, wherea=1 without control (black lines) and with control (blue lines) for α=1 (solid line), α=0.98 (dashed line), and α=0.95 (dotted line).Figure 7
Mortality due to CHD with different value ofa for α=1.


## 8. Conclusion

In our work, we discussed that a mathematical model contains the time FD, which describes the relationship between fish consumption and CHD mortality in Egypt. Moreover, we introduced control function which represents the metabolic (obesity, diabetes, raised lipids, and high blood pressure) and behavioural (unhealthy diet, physical inactivity, tobacco use, and alcohol) risk factors of CHD that help to reduce the number of mortality due to CHD. We have investigated the equilibria and stability analysis of our model, where we found that in the absence of biomass of fish population (i.e.,X=0), the EP is unstable. Otherwise, in the presence of the biomass of fish population X∗ and under the condition that the fishing effort E is less than the Biotechnical Productivity b/cα, the EP is stable. We present a general formulation for FOCP and the fractional-order necessary optimality conditions are derived using Pontryagin maximum principle. The numerical method used for this purpose is FBSM based on PCM where it is considered a simplified method to solve FOCPs with left CFD. We proposed some strategies with different values of the rate of fish consumed from the total harvested with and without control in order to obtain the best strategy to reduce CHD mortality. It is found that, whenever the consumption of fish from the total harvested fish increased with avoiding the metabolic and behavioural risk factors of CHD, then the CHD mortality rate is significantly reduced in a specified period of time. The effect of FD α (where α be chosen close to 1) on the studied model is displayed in the above figures to reduce the mortality due to CHD besides fish consumption and suggested control.



---

*Source: 9823753-2020-02-03.xml*