Optimal Control Problem of a Non-integer Order Waterborne Pathogen Model in Case of Environmental Stressors

1. Introduction

Water is one of the most required natural resources for life and it is not only vital for drinking but also needed for other daily uses. Every individual deserves safe water; but, unfortunately it has been estimated that at least 2 billion people have to use contaminated drinking water [1]. According to “Water for health: taking charge” report of WHO, about 3.4 million people died due to water related illnesses [2]. Specifically, almost 6,000 children die every day due to lack of clean water and hygienic conditions [3]. On the other hand, 35% of health care facilities in low and middle-income countries do not have water and soap, and half of the world's population is predicted to live in water-stressed zones by 2025 [3].

Mathematical models offer new ways to investigate the population dynamics in case of pollution related diseases. With the use of such models, strategic decisions can be tested beforehand, and we can reduce cost and save time. For example, in 2003, influence of stress to spread of infectious diseases is investigated [4]. In 2010, a SIR model is extended by including the pathogen population and, person–person and person–water–person interactions are modeled. On the other hand, spread of cholera due to unclean water is discussed for a delayed differential equation (DE) together with stability analysis [5]. The effect of heterogeneity by neglecting person-to-person transmission is studied [6]. On the other hand, the spread of virus/worm in computer networks have been modeled with DEs and optimal attack strategies have been decided [7] and the result can be useful to set network security protocols up. Another model including the influence of socioeconomic classes is constructed for multi patches [8]. As a different way to model waterborne diseases, networks are incorporated to express both indirect environment-to-human and direct human-to-human transmission routes [9]. As a two-dimensional model, we can mention the study on a system of reaction-diffusion equations published in 2018 [10].

Fractional differentiation and integration operators can be defined as the generalizations of their integer counterparts, where we use the fractional order to replace with the integer order [11–13]. A review including a systematic classification and applications can be found in the studies [14, 15]. Indeed, the integer-order model is a special case of the generalized fractional order model where the solution of the fractional order model must converge to the solution of the integer order model as the order approaches one. On the other hand, fractional DEs are suitable to express the real world phenomena which are influenced by hereditary and memory properties that cannot be expressed with integer order models [16]. Some application areas can be noted as physics [17], viscoelasticity [18], control systems [19], economics [20], finance [21], biology [22, 23], and bioengineering [24]. Specifically, anomalous dispersion is modeled with Riemann–Liouville derivative in space [25]. A pn semiconductor diode is described under sinusoidal operational conditions with the use of fractional calculus and some experimental results are compared with the model [26]. Viscoelastic materials are modeled with a variable order DE [27] and Caputo-Fabrizio derivative [28]. Three different fractional operators are used to investigate the properties of blood ethanol concentration [29]. Coexistence of tuberculosis and diabetes mellitus has been examined [30], while a honeybee population growth model is presented [31]. On the other hand, the construction of diffusion equation is obtained with the Atangana–Baleanu fractional integral [32]. Existence and uniqueness of Caputo time fractional Navier-Stokes equation is proven [33], whereas flow characteristics of gas and oil are discussed for a heterogeneous media [34].

For some real world phenomena, data is available and it is quite straightforward to test contribution of fractional derivative. For example, Diethelm has compared the numerical solution of the integer order and fractional order models with the data and he has revealed that the fractional order model fits the data better [35]. Then, this study has been extended with the use of other fractional derivatives [36]. As another epidemic disease, namely Ebola, has been modeled with a Caputo fractional derivative and the numerical results agree well with the data [37]. On the other hand, economic growth of the Group of Seven (G7) in 1973–2016 period has been modeled with fractional order derivative and it has been underlined by comparison with the data that the fractional model gives better results than the integer order model due to memory effect [38]. In addition, the data of varicella disease outbreak in Shenzhen city of China in 2015 has been used to compare the numerical solution of the fractional order model and success of the fractional derivatives in modeling such epidemiological phenomena has been exemplified [29].

The topic of fractional optimal control problems (FOCPs) has started to gain interest. The number of studies on FOCPs is quite less when compared to the studies on the fractional DEs; but, it has been an emerging research field for a while. For example, we suggest the studies on optimization of biological reactive systems [39, 40], where construction of a model and derivation of its numerical solution have been discussed. Furthermore, we can mention some studies on tumor dynamics [41], trajectory tracking [42], fractional–order sliding mode controller [43], virus infection [44], an HIV/AIDS epidemic model [45], dynamics of Human African trypanosomiasis [46], the light tracking systems [47], and human respiratory syncytial virus infection [44].

In 2019, transmission of waterborne diseases is expressed for susceptible compartment, susceptible stressed compartment and infected compartment with the inclusion of pathogen population and stress related negative effects [48]. Motivated by this recent work, a new non-integer order model is constructed with the addition of recovered subgroup in this current study. For the new model, we prove that the solution is non-negative and bounded from above. Then, local stability of the disease free equilibrium point has been analyzed with the use of the reproduction number. On account of unfortunate facts reported by WHO [1, 3], it is inevitable to figure some strategies out to stop any harm of water pollution. Thus, we apply optimal control theory to decide the optimal intervention strategies to minimize the number of infected individuals and to reduce the pathogen population. We believe that this study can be a nice start to optimize stress related models to found healthier generations of all creatures. The rest of this study is organized as follows: In section 2, we construct a non-integer order model with constant controls. Then, we prove that the solution is bounded and we present stability analysis for the disease free equilibrium point. In section 3, we construct an OCP and we derive the optimality conditions. Section 4 is devoted to some illustrative examples for testing purposes. Then, we conclude the study.

2. Model formulation

In this section, we derive a non-integer order waterborne pathogen model which is an extended version of the one in the study [48] with the inclusion of recovered compartment. Environmental pollution is created in several years due to accumulated causes. Here, we use time fractional derivative to express temporal change both in the host and pathogen population.

We define the (left) generalized derivative with Caputo fractional derivative for t > a as Podlubny [12]

$\begin{array}{ll}{}_a{{D}}_t^{q}g(t)=\{\begin{array}{cc}\frac{1}{\Gamma (1-q)}{\int }_a^t\frac{g^{\prime }(s)}{{(t-s)}^q}ds,& 0<q<1,\\ g^{\prime }(t),\hfill & q=1,\end{array}& (2.1)\end{array}$

where the Gamma function Γ :(0, ∞) → ℝ is defined as $\Gamma (t)={\int }_0^{\infty }{e}^{-u}{u}^{t-1}du.$ The corresponding right differentiation operator is defined as

$\begin{array}{ll}{}_t{{D}}_b^{q}g(t)=\{\begin{array}{cc}-\frac{1}{\Gamma (1-q)}{\int }_t^b\frac{g^{\prime }(s)}{{(s-t)}^q}ds,& 0<q<1,\\ -g^{\prime }(t),\hfill & q=1.\end{array}& (2.2)\end{array}$

We mention the generalized right Riemann-Liouville (RL) derivative as

$\begin{array}{ll}{{\displaystyle }}_t^R{{D}}_{t_f}^{q}g(t)=\{\begin{array}{cc}-\frac{1}{\Gamma (1-q)}\frac{d}{dt}{\int }_t^b\frac{g(s)}{{(s-t)}^q}ds,& 0<q<1,\\ -g^{\prime }(t),\hfill & q=1.\end{array}& (2.3)\end{array}$

Before presenting the model, we note a relation between the generalized Caputo and RL derivatives as

$\begin{array}{ll}{{\displaystyle }}_t^R{{D}}_{t_f}^{q}g(t)={}_t{{D}}_{t_f}^{q}g(t)+\frac{g(t_f)}{\Gamma (1-q)}(t_f-t{)}^{-q}.& (2.4)\end{array}$

The model is constructed by dividing the total population N(t) into four compartments, namely, susceptible population S₁(t), susceptible population influenced by environmental stressors S₂(t), infected population I(t), and recovered population R(t). To discuss the contribution of pollution to infection, pathogen pollution (Cells/liter) is included and it is denoted by P(t). The birth rate of susceptible individuals is fixed as Λ. A proportion (1 − p) of newborns joins S₁ class, whereas the rest contributes to S₂ class. Due to stressors, susceptible individuals enter stressors group at a rate of θ. The natural mortality rate and disease induced death rate are given by μ and α, respectively. Infectious individuals join recovered subgroup at a rate of ϕ. People are recovered at a rate of γ and a proportion r of them enters into susceptible class, while the rest joins stressors compartment. On the other hand, the disease is transmitted at a rate of β and transmission is affected due to stress related parameters, namely level of pollution affecting transmission ω and effect of pollution ρ. Moreover, the pathogen population grows at a rate of ξ and it decays at a rate of η. In Table 1, we present the description, units and values of the parameters.

TABLE 1

Table 1. Parameter values.

Then, we reach the following model:

$\begin{array}{l}{}_0{{D}}_t^q{S}_1(t)=(1-p){\Lambda }^q-{\theta }^q{S}_1(t)-{\beta }^q{S}_1(t)P(t)\\ -{\mu }^q{S}_1(t)+r{\gamma }^{q}R(t),& (\mathrm{2.5a})\end{array}$

$\begin{array}{l}{}_0{{D}}_t^q{S}_2(t)=p{\Lambda }^q+{\theta }^q{S}_1(t)-{\beta }^q(1+\omega \rho (1-u_1))S_2(t)P(t)\\ -{\mu }^q{S}_2(t)+(1-r){\gamma }^{q}R(t),& (\mathrm{2.5b})\end{array}$

$\begin{array}{l}{}_0{{D}}_t^{q}I(t)={\beta }^q{S}_1(t)P(t)+{\beta }^q(1+\omega \rho (1-u_1))S_2(t)P(t)\\ -({\varphi }^q+{\alpha }^q+{\mu }^q)I(t),& (\mathrm{2.5c})\end{array}$

$\begin{array}{ll}{}_0{{D}}_t^{\alpha }R(t)={\varphi }^{q}I(t)-({\gamma }^q+{\mu }^q)R(t),& (\mathrm{2.5d})\end{array}$

$\begin{array}{ll}{}_0{{D}}_t^{\alpha }P(t)={\xi }^{q}I(t)-{\eta }^q(1+u_2)P(t),& (2.5e)\end{array}$

$\begin{array}{l}{S}_1(0)=S_{1,0},S_2(0)=S_{2,0},I(0)=I_0,R(0)=R_0,P(0)=P_0,\end{array}$

where u₁ and u₂ stand for the fixed control terms. We note that the model is expressed for a generalized derivative. Indeed, it is possible to investigate the model for both classical and fractional order derivative. The reason of using fractional derivative is its success in a good fit to available data as discussed in the literature.

Remark 2.1. We note that it is possible to obtain analytical solution of some epidemiological model (such as SIS models) via a Lie algebra methodology. Due to lack of symmetry in more complicated models, this method is not widely used in biological or social systems. For the readers interested in this technique, we can suggest some studies as Shang [49, 50].

2.1. Basic Properties

To prove the non-negative solution of the model, we need the generalized mean value theorem and the corollary.

Lemma 2.2. Odibat and Shawagfeh [51] Let g(x) ∈ C[a, b] and ${}_a^C{{D}}_t^{q}g(t)\in C(a,b]$ for 0 < q ≤ 1. Then, for a ≤ s ≤ b and ∀ x ∈ (a, b], the following estimate holds:

$\begin{array}{ll}g(x)=g(a)+\frac{1}{\Gamma (q)}{(}_a^C{{D}}_t^{q}g)(s)(x-a{)}^q.& (2.6)\end{array}$

Corollary 2.3. Let g(x) ∈ C[a, b] and ${}_a^C{{D}}_t^{q}g(t)\in C(a,b]$ for 0 < q ≤ 1. If ${}_a^C{{D}}_t^{q}g(t)$ is non-negative ∀ x ∈ (a, b), then g(x) is non-decreasing for each x ∈ [a, b]. If ${}_a^C{{D}}_t^{q}g(t)$ is non-positive ∀ x ∈ (a, b), then g(x) is non-increasing for each x ∈ [a, b].

To show that the solution is bounded from above, we use the Laplace transform. The Laplace transform of the (left) Caputo derivative is obtained as

$\begin{array}{ll}{L}{\{}_a^C{{D}}_t^{q}g(t)\}=s^{q}G(s)-g(0)s^{q-1}.& (2.7)\end{array}$

Moreover, the Laplace transform of the Mittag-Leffler function is given by

$\begin{array}{ll}{L}\left\{t^{p-1}{E}_{q,p}(-a{t}^q)\right\}=\frac{s^{q-p}}{s^q+a},& (2.8)\end{array}$

where $E_{q,p}(z)=\sum _{i=0}^{\infty }\frac{z^i}{\Gamma (qi+p)}$ and $E_{q,p}(z)=z{E}_{q,q+p}(z)+\frac{1}{\Gamma (p)}$.

Theorem 2.4. Let ${X}(t)=(S_1(t),S_2(t),I(t),R(t),P(t))$ be the unique solution of the model (2.5) for t ≥ 0. Then, the solution ${X}(t)$ remains in ${ℝ}_{+}^5$.

Proof: By the study [52, Thm. 3.1, Remark 3.2], the solution to the model (2.5) is unique for t > 0. On the other hand, we observe that

$\begin{array}{ll}{}_0{{D}}_t^q{S}_1{|}_{S_1=0}=(1-p){\Lambda }^q+r{\gamma }^{q}R(t)\ge 0& (\mathrm{2.9a})\end{array}$

$\begin{array}{ll}{}_0{{D}}_t^q{S}_2{|}_{S_2=0}=p{\Lambda }^q+{\theta }^q{S}_1(t)+(1-r){\gamma }^{q}R(t)\ge 0,& (\mathrm{2.9b})\end{array}$

$\begin{array}{ll}{}_0{{D}}_t^{q}I{|}_{I=0}={\beta }^q{S}_1(t)P(t)+{\beta }^q(1+\omega \rho (1-u_1))S_2(t)P(t)\ge 0,& (\mathrm{2.9c})\end{array}$

$\begin{array}{ll}{}_0{{D}}_t^{q}R{|}_{R=0}={\varphi }^{q}I(t)\ge 0,& (\mathrm{2.9d})\end{array}$

$\begin{array}{ll}{}_0{{D}}_t^{q}P{|}_{P=0}={\xi }^{q}I(t)\ge 0.& (2.9e)\end{array}$

It means that the vector field points to ${ℝ}_{+}^5$ and the solution lies in ${ℝ}_{+}^5$ for a non-negative initial condition [51].□

Theorem 2.5. Let ${X}(t)=(S_1(t),S_2(t),I(t),R(t),P(t))$ be the solution of the model (2.5) for t > 0. Then, the solution ${X}(t)$ is bounded from above.

Proof: We add the Equations (2.5a)-(2.5d) to reach

$\begin{array}{ll}{}_0{{D}}_t^{q}N(t)={\Lambda }^q-{\mu }^{q}N(t)-{\alpha }^{q}I(t)\le {\Lambda }^q-{\mu }^{q}N(t).& (2.10)\end{array}$

Using the Laplace transform (2.7), we obtain

$\begin{array}{ll}{s}^q{L}(N(t))-s^{q-1}N(0)\le \frac{{\Lambda }^q}{s}-{\mu }^q{L}(N(t)).& (2.11)\end{array}$

Arranging, we get

$\begin{array}{l}{L}(N(t))\le {\Lambda }^q\frac{s^{-1}}{s^q+{\mu }^q}+N(0)\frac{s^{q-1}}{s^q+{\mu }^q}={\Lambda }^q\frac{s^{q-(1+q)}}{s^q+{\mu }^q}\\ +N(0)\frac{s^{q-1}}{s^q+{\mu }^q}.& (2.12)\end{array}$

Applying inverse Laplace transform to (2.12), we obtain

$\begin{array}{l}N(t)\le {{L}}^{-1}\left\{{\Lambda }^q\frac{s^{q-(1+q)}}{s^q+{\mu }^q}\right\}+{{L}}^{-1}\left\{N(0)\frac{s^{q-1}}{s^q+{\mu }^q}\right\}\\ {\displaystyle \underset{(2.8)}{{\displaystyle \underbrace{\le }}}}{\Lambda }^q{t}^q{E}_{q,q+1}(-\mu t^q)+N(0)E_{q,1}(-\mu t^q)\\ \le \frac{{\Lambda }^q}{{\mu }^q}{\mu }^q{t}^q{E}_{q,q+1}(-\mu t^q)+N(0)E_{q,1}(-\mu t^q)\\ \le max\left\{\frac{{\Lambda }^q}{{\mu }^q},N(0)\right\}({\mu }^q{t}^q{E}_{q,q+1}(-\mu t^q)+E_{q,1}(-\mu t^q))\\ =\frac{C}{\Gamma (1)}=C,& (2.13)\end{array}$

where $C=max\{\frac{{\Lambda }^q}{{\mu }^q},N(0)\}$. Thus, N(t) is bounded from above which means that the solutions S₁(t), S₂(t), I(t) and R(t) are bounded from above.

We proceed with the equation for pathogen population (2.5e), which is written as

$\begin{array}{ll}{}_0{{D}}_t^{q}P(t)={\xi }^{q}I(t)-{\eta }^q(1+u_2)P(t)\underset{(2.13)}{\underbrace{\le }}{\xi }^{q}C-{\eta }^{q}P(t).& (2.14)\end{array}$

Then, the similar technique is applied to show that the pathogen pollution P(t) is bounded from above.□

2.2. Stability Analysis of the Disease Free Equilibrium Point

The disease free equilibrium (DFE) point ${{E}}^0=(S_1^0,S_2^0,0,0,0)$ is obtained by equating the right-hand side of the Equations (2.5a)–(2.5b) to zero. Then, we obtain

$\begin{array}{ll}{S}_1^0=\frac{(1-p)\Lambda }{\theta +\mu },S_2^0=\frac{\Lambda (p\mu +\theta )}{\mu (\theta +\mu )}.& (2.15)\end{array}$

To discuss the local stability of the DFE point ${{E}}^0$, we need to obtain the basic reproduction number, namely a parameter measuring the contribution of a single infection to secondary infections [53]. We define the matrices for the new infections and the other terms for the Equations (2.5c)–(2.5e) as

$\begin{array}{l}{F}=\left(\begin{array}{c}\beta S_{1}P+{\beta }^q(1+\omega \rho (1-u_1))S_{2}P\\ 0\end{array}\right),\\ {V}=\left(\begin{array}{c}({\gamma }^q+{\alpha }^q+{\mu }^q)I\\ -{\xi }^{q}I+{\eta }^q(1+u_2)P\end{array}\right).& (2.16)\end{array}$

We note that the contribution of the recovered compartment does not appear here, since recovered individuals do not directly join the infected subgroup. Then, partial derivatives of these matrices with respect to I and P are obtained and then they are evaluated at the DFE point to obtain

$\begin{array}{l}F=\left(\begin{array}{cc}0& \beta S_1^0+{\beta }^q(1+\omega \rho )S_2^0\\ 0& 0\end{array}\right),\\ V=\left(\begin{array}{cc}({\gamma }^q+{\alpha }^q+{\mu }^q)& 0\\ -{\xi }^q& {\eta }^q\end{array}\right).& (2.17)\end{array}$

We reach the next generation matrix FV⁻¹ and its spectral radius is given by Sharma and Kumari [48]

$\begin{array}{ll}{{R}}_0=\frac{{\xi }^q{\beta }^q{\Lambda }^q({\mu }^q+{\theta }^q(1+\omega \rho (1-u_1)))+(1-u_1)\omega \rho p{\mu }^q}{(1+u_2){\eta }^q{\mu }^q({\theta }^q+{\mu }^q)({\gamma }^q+{\alpha }^q+{\mu }^q)}.& (2.18)\end{array}$

We mention a result to derive stability analysis [54].

Lemma 2.6. The DFE point ${{E}}^0$ is locally asymptotically stable if all eigenvalues λ_i of the linearization matrix of the model (2.5) satisfy $|arg({\lambda }_i)|>q\frac{\pi }{2}$.

Now, we prove the following result for the DFE point.

Theorem 2.7. The DFE point ${{E}}^0$ of the model (2.5) is locally asymptotically stable if ${{R}}_0<1$ and unstable if ${{R}}_0>1$.

Proof: The Jacobian matrix associated with the model (2.5) is obtained as

$\begin{array}{ll}J({{E}}^0)=\left(\begin{array}{ccccc}-({\theta }^q+{\mu }^q)& 0& 0& r{\gamma }^q& -\beta S_1^0\\ {\theta }^q& -{\mu }^q& 0& (1-r){\gamma }^q& -{\beta }^q(1+\omega \rho (1-u_1))S_2^0\\ 0& 0& -({\varphi }^q+{\alpha }^q+{\mu }^q)& 0& {\beta }^q(1+\omega \rho (1-u_1))S_2^0\\ 0& 0& {\varphi }^q& -({\gamma }^q+{\mu }^q)& 0\\ 0& 0& {\xi }^q& 0& -{\eta }^q(1+u_2)\end{array}\right).& (2.19)\end{array}$

The characteristic equation of this matrix is given by p(λ) = (λ + μ^q)(λ + θ^q + μ^q)(λ + γ^q + d^q)f(λ) and three eigenvalues are found as −μ^q, −(θ^q + μ^q) and −(γ^q + μ^q), which satisfy Lemma 2.6. Moreover, the quadratic part is written as $f(\lambda )={\lambda }^2+c_1\lambda +c_0$, where $c_1={\eta }^q(1+u_2)+{\varphi }^q+{\alpha }^q+{\mu }^q$, $c_0={\eta }^q(1+u_2)({\varphi }^q+{\alpha }^q+{\mu }^q)(1-{{R}}_0)$. By Routh-Hurwitz criterion, we check that c₁ > 0 and c₀ > 0 with ${{R}}_0<1$, which shows that the roots of f(λ) are negative. Then, we prove the desired result.□

3. Optimal control problem

We investigate some suitable policies to reduce the reproduction number associated with the model. Specifically, we optimize the decay rate of pathogens and the stress related parameters with the aim of minimizing the number of infected individuals and reducing the pathogen population in size with the use of non-constant controls u₁(t) and u₂(t), respectively. Therefore, we determine their optimal values with the following cost functional:

$\begin{array}{l}{\displaystyle \underset{(u_1,u_2)\in U_{ad}}{min}}J(u_1,u_2)=I(t_f)+P(t_f)\\ +{\displaystyle {\int }_0^{t_f}}\left({\omega }_{1}I(t)+{\omega }_{2}P(t)+\frac{{\omega }_3}{2}{u}_1^2(t)+\frac{{\omega }_4}{2}{u}_2^2(t)\right)dt,& (3.1)\end{array}$

subject to the system

$\begin{array}{l}{}_0{{D}}_t^q{S}_1(t)=(1-p){\Lambda }^q-{\theta }^q{S}_1(t)-{\beta }^q{S}_1(t)P(t)-{\mu }^q{S}_1(t)\\ +r{\gamma }^{q}R(t),& (\mathrm{3.2a})\end{array}$

$\begin{array}{l}{}_0{{D}}_t^q{S}_2(t)=p{\Lambda }^q+{\theta }^q{S}_1(t)-{\beta }^q(1+\omega \rho (1-u_1(t)))S_2(t)P(t)\\ -{\mu }^q{S}_2(t)+(1-r){\gamma }^{q}R(t),& (\mathrm{3.2b})\end{array}$

$\begin{array}{l}{}_0{{D}}_t^{q}I(t)={\beta }^q{S}_1(t)P(t)+{\beta }^q(1+\omega \rho (1-u_1(t)))S_2(t)P(t)\\ -({\varphi }^q+{\alpha }^q+{\mu }^q)I(t),& (\mathrm{3.2c})\end{array}$

$\begin{array}{ll}{}_0{{D}}_t^{\alpha }R(t)={\varphi }^{q}I(t)-({\gamma }^q+{\mu }^q)R(t),& (\mathrm{3.2d})\end{array}$

$\begin{array}{ll}{}_0{{D}}_t^{\alpha }P(t)={\xi }^{q}I(t)-{\eta }^q(1+u_2(t))P(t),& (3.2e)\end{array}$

$\begin{array}{l}{S}_1(0)=S_{1,0},S_2(0)=S_{2,0},I(0)=I_0,R(0)=R_0,P(0)=P_0,\end{array}$

where the admissible space of controls is given by

$\begin{array}{l}{U}_{ad}=\{(u_1(t),u_2(t))|(u_1(t),u_2(t))\text{ is measurable with }\\ 0.05\le u_1(t)\le 0.95,0\le u_2(t)\le 1\text{ for all }t\in [0,t_f]\}.\end{array}$

Thus, the optimal control $u^{*}(t)=(u_1^{*}(t),u_2^{*}(t))\in U_{ad}$ is required so that $J(u^{*})={min}_{(u_1(t),u_2(t))\in U_{ad}}J(u_1(t),u_2(t))$ is reached. We observe that $\frac{\partial {{R}}_0}{\partial u_1}<0$ and $\frac{\partial {{R}}_0}{\partial u_2}<0$. Then, we can determine the optimal values of u₁ and u₂ to make ${{R}}_0<1$.

3.1. Optimality System

We use Pontryagin's maximum principle [55–57] to obtain the necessary optimality conditions where the adjoint (costate) functions attach the state equation to the cost functional J. Indeed, this technique allows us to express the control in terms of the state and the adjoint functions. Moreover, Pontryagin's maximum principle interprets the constraint minimization problem (3.1)-(3.2) as another minimization problem (Hamiltonian) with respect to the controls. Therefore, we construct the Hamiltonian as

$\begin{array}{l}{H}({X},{U},{P})=\left({\omega }_{1}I(t)+{\omega }_{2}P(t)+\frac{{\omega }_3}{2}{u}_1^2(t)+\frac{{\omega }_4}{2}{u}_2^2(t)\right)\\ +{\lambda }_1^T((1-p){\Lambda }^q-{\theta }^q{S}_1(t)-{\beta }^q{S}_1(t)P(t)\\ -{\mu }^q{S}_1(t)+r{\gamma }^{q}R(t))\\ +{\lambda }_2^T(p{\Lambda }^q+{\theta }^q{S}_1(t)-{\beta }^q(1+\omega \rho (1-u_1))S_2(t)P(t)\\ -{\mu }^q{S}_2(t)+(1-r){\gamma }^{q}R(t))\\ +{\lambda }_3^T({\beta }^q{S}_1(t)P(t)+{\beta }^q(1+\omega \rho (1-u_1))S_2(t)P(t)\\ -({\varphi }^q+{\alpha }^q+{\mu }^q)I(t))\\ +{\lambda }_4^T({\varphi }^{q}I(t)-({\gamma }^q+{\mu }^q)R(t))\\ +{\lambda }_5^T({\xi }^{q}I(t)-{\eta }^q(1+u_2)P(t)),& (3.3)\end{array}$

where ${\lambda }_i{(t)}^{\prime }$s denote the adjoint variables associated with their corresponding states. The first part of the Hamiltonian ${H}$ is related to the integrand in (3.1) and its second part is constructed by multiplying each adjoint function ${\lambda }_i{(t)}^{\prime }$s with the corresponding right-hand side of the model (3.2). Then, Pontryagin's maximum principle leads to the following optimality system.

Theorem 3.1. Given an optimal control $u^{*}(t)=(u_1^{*}(t),u_2^{*}(t))\in U_{ad}$ and the state solution ${{X}}^{*}(t)=(S_1^{*}(t),S_1^{*}(t),I^{*}(t),R^{*}(t),P^{*}(t))$ corresponding to (3.2) that minimize J(u₁, u₂) in (3.1), there exists adjoint variable ${P}(t)=({\lambda }_1(t),{\lambda }_2(t),{\lambda }_3(t),{\lambda }_4(t),{\lambda }_5(t))$ satisfying

$\begin{array}{ll}{{\displaystyle }}_t^R{{D}}_{t_f}^q{\lambda }_1(t)=-({\theta }^q+{\beta }^{q}P(t)+{\mu }^q){\lambda }_1(t)+r{\gamma }^q{\lambda }_4(t)-{\beta }^q{X}_1{\lambda }_5(t),& (\mathrm{3.4a})\end{array}$

$\begin{array}{l}{{\displaystyle }}_t^R{{D}}_{t_f}^{\alpha }{\lambda }_2(t)={\theta }^q{\lambda }_1(t)-({\beta }^q(1+\omega \rho (1-u_1(t))P(t)+{\mu }^q){\lambda }_2(t)\\ +(1-r){\gamma }^q{\lambda }_4(t)-{\beta }^q(1+\omega \rho (1-u_1(t)))X_2(t){\lambda }_5(t),& (\mathrm{3.4b})\end{array}$

$\begin{array}{l}{{\displaystyle }}_t^R{{D}}_{t_f}^{\alpha }{\lambda }_3(t)={\beta }^{q}P(t){\lambda }_1(t)+{\beta }^q(1+\omega \rho (1-u_1(t)))P(t){\lambda }_2(t)\\ -({\varphi }^q+{\alpha }^q+{\mu }^q){\lambda }_3(t),\\ +({\beta }^q{X}_1(t)+{\beta }^q(1+\omega \rho (1-u_1(t))X_2(t)){\lambda }_5(t)+{\omega }_1,& (\mathrm{3.4c})\end{array}$

$\begin{array}{ll}{{\displaystyle }}_t^R{{D}}_{t_f}^{\alpha }{\lambda }_4(t)={\varphi }^q{\lambda }_3(t)-({\gamma }^q+{\mu }^q){\lambda }_4,& (\mathrm{3.4d})\end{array}$

$\begin{array}{ll}{{\displaystyle }}_t^R{{D}}_{t_f}^{\alpha }{\lambda }_5(t)={\xi }^q{\lambda }_3(t)-{\eta }^q(1+u_2(t)){\lambda }_5+{\omega }_2,& (3.4e)\end{array}$

with the transversality conditions

$\begin{array}{ll}{\lambda }_1(t_f)={\lambda }_2(t_f)={\lambda }_4(t_f)=0,{\lambda }_3(t_f)={\lambda }_5(t_f)=1.& (3.5)\end{array}$

Furthermore, the pair of optimal controls is given by

$\begin{array}{ll}{u}_1^{*}=min\left(max\left(({\lambda }_3-{\lambda }_2)\frac{{\beta }^q(1+\omega \rho )X_2(t)P(t)}{{\omega }_3},0.05\right),0.95\right),& (\mathrm{3.6a})\end{array}$

$\begin{array}{ll}{u}_2^{*}=min\left(max\left({\lambda }_5\frac{\eta P(t)}{{\omega }_4},0\right),1\right).& (\mathrm{3.6b})\end{array}$

Proof: The existence of an optimal control $u^{*}(t)=(u_1^{*}(t),u_2^{*}(t))$ and the associated state solution ${{X}}^{*}(t)=(S_1^{*}(t),S_1^{*}(t),I^{*}(t),R^{*}(t),P^{*}(t))$ comes from the convexity of J(u₁, u₂) with respect to u_i's and the Lipschitz continuous constraint [58], Theorem 5. 1. On the other hand, differentiation of Hamiltonian ${H}$ with respect to the state variable gives the adjoint system (3.4) as

$\begin{array}{l}{{\displaystyle }}_t^R{{D}}_{t_f}^q{\lambda }_1(t)=\frac{\partial {H}}{\partial S_1},{}_t^R{{D}}_{t_f}^q{\lambda }_2(t)=\frac{\partial {H}}{\partial S_2},{}_t^R{{D}}_{t_f}^q{\lambda }_3(t)=\frac{\partial {H}}{\partial I},\\ {{\displaystyle }}_t^R{{D}}_{t_f}^q{\lambda }_4(t)=\frac{\partial {H}}{\partial R},{}_t^R{{D}}_{t_f}^q{\lambda }_5(t)=\frac{\partial {H}}{\partial P},\end{array}$

with λ₁(t_f) = λ₂(t_f) = λ₄(t_f) = 0, λ₃(t_f) = λ₅(t_f) = 1 due to the definition of the cost functional. The optimality conditions is obtained by differentiating Hamiltonian ${H}$ with respect to the control variables u₁ and u₂ and then we project the resulting equations onto U_ad to reach (3.6).□

4. Numerical results

In this section, we present some numerical results to test the contribution of optimal control theory. As the optimization algorithm, forward-backward sweep method is applied [59], Chap. 5 [60] and the OCP is discretized with L1-method with Δt = 0.02, which is a variant of implicit Euler method [61]. We solve the OCP for 100 days and we choose ω₁ = 100 and ω₂ = 1000, since it is more important to reduce the pathogen pollution, which is the main source of infection, than the number of infected individuals. Moreover, we balance the weight constants ω₃ and ω₄ according to the cost of the corresponding strategy. Here, we assume that the cost of controlling the decay rate of pathogens is higher than the cost of reducing the stress related parameters, since stress acts on a smaller population than pollution (including children). Thus, we fix ω₃ < ω₄. To guarantee the reduction in I(t) and P(t) at the final time, we add the term I(t_f) + P(t_f) to the cost functional. On the other hand, to optimize the values of the stress related parameters, we set the bounds for u₁(t) as [0.5, 0.95], since neither it is possible to eliminate stress totally nor it is a force acting fully all the time. We set the bounds for u₂(t) as [0, 1] to be more strict. We fix the initial subpopulations as S_{1, 0} = 100, S_{2, 0} = 250, I₀ = 50, R₀ = 5, and P₀ = 90, which are taken close to the equilibrium components in (Sec.10, [48]), section 10, except for R₀. Here, we set up a scenario where the number of susceptible stressed individuals is greater than the susceptible compartment and the number of recovered individuals is quite small, since our goal is to test the strength of optimal control theory. We discuss two different control strategies and both of them have different costs and contribution. In addition to applying them simultaneously, single control will be tested, too.

With the parameters in Table 1, the reproduction number is calculated as

${{R}}_0=\left\{1.5867,1.9709,2.4463,3.0341\right\},$

for q ∈ {0.7, 0.8, 0.9, 1}, respectively. As q increases, the disease becomes more endemic. Our goal is to find optimal values of the controls so that ${{R}}_0<1$ can be reached.

Before discussing the optimal control strategy, we present the numerical solutions of I(t) and P(t) with u₁ = u₂ = 0 in Figure 1, which can be called the uncontrolled problem. We observe that I(t) and P(t) increase for q = {0.9, 1}, while they decrease for q = {0.7, 0.8} as time passes. It is quite vital to apply optimal control for q = {0.9, 1} to eliminate the spread of the disease. Here, the order q may represent the severity of pollution. As q increases, the pollution becomes more dominant and it leads to a rapid increase in the number of infected people.

FIGURE 1

Figure 1. (A,B) I(t) and P(t) with u₁(t) = u₂(t) = 0.

We start with the single control u₁(t) and present the results in Figure 2. We immediately observe that there is not an obvious difference between the results associated with the uncontrolled problem and the current case. Therefore, trying to decrease the stress related parameters is not a successful action to reduce the number of infected individuals and pathogen population in size, although we apply the full control for almost 100 days.

FIGURE 2

Figure 2. (A,B) I(t) and P(t) with single optimal control u₁(t) (C).

We proceed with the control u₂(t) and the corresponding figures are depicted in Figure 3. The single control u₂(t) is more powerful than the control u₁(t) due to rapid decrease in both I(t) and P(t). Here, we can catch the impact of the fractional order to the choice of the control strategies. As we increase q, control must acts less. The decay rate of pathogen must be controlled strictly almost until the 60th day for q = 1. Then, as we decrease the order q, the control must be applied more strictly. Thus, q can be associated with awareness or reception toward some policies to eliminate pollution.

FIGURE 3

Figure 3. (A,B) I(t) and P(t) with single optimal control u₂(t) (C).

Now, we proceed with the last case, namely the combination of controls u₁(t) and u₂(t). From the Figure 4, we observe success of the optimal control. The contribution of u₁(t) is almost the same when we compare the numerical solutions in Figures 3, 4. Controlling the stress related parameters has a positive effect on the control of decay rate of pathogen. Due to this contribution, u₂(t) acts fully not until the 60th day. Indeed, it starts to decrease earlier with the impact of the control u₁(t). Moreover, for the model with a lower order of differentiation, u₁(t) and u₂(t) must act more strictly.

FIGURE 4

Figure 4. (A,B) I(t) and P(t) with optimal controls u₁(t), u₂(t) (C,D).

We revealed the success of the OCP in the figures above. Now, we examine the results in detail and show some values of the FOCP with respect to the order of differentiation in Table 2. We firstly observe that the reduction in the cost functional is the smallest one for u₁(t). For the control u₂(t) and the combined control, the reduction is almost the same; but, the combined control improves the result a little bit more. For the control u₂(t) and the pair of controls, as we increase the order, the value of the cost functional decreases and the reduction increases. We remember that the number of infected individuals and the pathogen population decrease in size for q∈{0.7, 0.8} as t approaches 100 days in case of no controls (see Figure 2). Despite of this decrease, the optimal control theory is more powerful for higher values of the fractional order.

TABLE 2

Table 2. Values of the cost functional J and reduction in J.

5. Summary and conclusion

We extend a non-integer order model for waterborne diseases motivated by the recent work [48]. The solution is proven to be bounded. On the other hand, we find the reproduction number ${{R}}_0$ and we show that the DFE point is locally asymptotically stable if ${{R}}_0<1$. Then, we construct an OCP to minimize the number of infected individuals and the pathogen pollution where the controls are used to optimize the stress related parameters and the decay rate of pathogen. We observe that optimizing only stress related parameters is not a successful strategy. Indeed, the value of the decay rate of pathogen is more critical to achieve the goal of this study. For the combined control, the best results are achieved; but, it is not more different than the ones obtained with the control u₂(t). Therefore, we deduce that more attention must be paid to control the decay rate of pathogen within a polluted habitat and optimizing the values of the stress related parameters can be a supportive strategy. In addition, we note that as the order q is increased, the disease becomes more endemic. Here, the order q may represent the severity of pollution. As q increases, the pollution affects the population worse and it leads to a rapid increase in the number of infected people. On the other hand, control acts fully for shorter time as we increase q. Thus, q can be associated with awareness or reception toward some policies to eliminate pollution. As a future work, comparison of our numerical results with the field data (if available) can be mentioned. Moreover, the use of other types of fractional derivatives can be investigated for a better data fit and for more accurate control policies.

Data Availability

All datasets generated for this study are included in the manuscript and/or the supplementary files.

Author Contributions

The author confirms being the sole contributor of this work and has approved it for publication.

Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The author is grateful to three referees for their valuable comments and helpful suggestions, which have helped to improve the presentation of this work significantly. Moreover, the author thanks the publisher for the full waiver.

References

1. WHO. Drinking-water Fact Sheet. Geneva: WHO (2018).

2. WHO. Water for Health : Taking Charge. Geneva: WHO (2001).

Google Scholar

3. WHO. Progress on Drinking Water, Sanitation and Hygiene. Geneva: WHO (2017).

Google Scholar

4. Lafferty KD, Holt RD. How should environmental stress affect the population dynamics of disease? Ecol Lett. (2003) 6:654–64. doi: 10.1046/j.1461-0248.2003.00480.x

CrossRef Full Text | Google Scholar

5. Misra A, Singh V. A delay mathematical model for the spread and control of water borne diseases. J Theor Biol. (2012) 301:49–56. doi: 10.1016/j.jtbi.2012.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Collins O, Govinder K. Incorporating heterogeneity into the transmission dynamics of a waterborne disease model. J Theor Biol. (2014) 356:133–43. doi: 10.1016/j.jtbi.2014.04.022

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Shang Y. Optimal control strategies for virus spreading in inhomogeneous epidemic dynamics. Can Math Bull. (2013) 56:621–9. doi: 10.4153/CMB-2012-007-2

CrossRef Full Text | Google Scholar

8. Collins O, Robertson SL, Govinder K. Analysis of a waterborne disease model with socioeconomic classes. Math Biosci. (2015) 269:86–93. doi: 10.1016/j.mbs.2015.08.016

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Wang Y, Cao J. Global dynamics of a network epidemic model for waterborne diseases spread. Appl Math Comput. (2014) 237:474–88. doi: 10.1016/j.amc.2014.03.148

CrossRef Full Text | Google Scholar

10. Zhou J, Yang Y, Zhang T. Global dynamics of a reaction–diffusion waterborne pathogen model with general incidence rate. J Math Anal Appl. (2018) 466:835–59. doi: 10.1016/j.jmaa.2018.06.029

CrossRef Full Text | Google Scholar

11. Oldham K, Spanier J. The Fractional Calculus Theory and Applications of Differentiation and Integration to Arbitrary Order. New York, NY; London: Academic Press (1974).

Google Scholar

12. Podlubny I. Fractional differential equations. In: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications. Vol. 198, Mathematics in Science and Engineering. San Diego, CA: Academic Press, Inc. (1999).

Google Scholar

13. Samko SG, Kilbas AA, Marichev OI. Fractional Integrals and Derivatives Vol. 1993 (1993).

Google Scholar

14. Teodoro GS, Machado JT, de Oliveira EC. A review of definitions of fractional derivatives and other operators. J Comput Phys. (2019) 388:195–208. doi: 10.1016/j.jcp.2019.03.008

CrossRef Full Text | Google Scholar

15. Sun H, Zhang Y, Baleanu D, Chen W, Chen Y. A new collection of real world applications of fractional calculus in science and engineering. Commun Nonlinear Sci Numer Simul. (2018) 64:213–31. doi: 10.1016/j.cnsns.2018.04.019

CrossRef Full Text

16. Du M, Wang Z, Hu H. Measuring memory with the order of fractional derivative. Sci Rep. (2013) 3:3431. doi: 10.1038/srep03431

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Hilfer R. Applications of Fractional Calculus in Physics. Singapore: World Scientific (2000).

Google Scholar

18. Mainardi F. Fractional Calculus and Waves in Linear Viscoelasticity: An Introduction to Mathematical Models. London: World Scientific (2010).

Google Scholar

19. Monje CA, Chen Y, Vinagre BM, Xue D, Feliu-Batlle V. Fractional-order Systems and Controls: Fundamentals and Applications. New York, NY: Springer (2010).

Google Scholar

20. Tarasova VV, Tarasov V. T. Concept of dynamic memory in economics. Commun Nonlinear Sci Numer Simul. (2018) 55:127–45. doi: 10.1016/j.cnsns.2017.06.032

CrossRef Full Text | Google Scholar

21. Hu Y, Øksendal B. Fractional white noise calculus and applications to finance. Infinite Dimens Anal Quant Probab Relat Top. (2003) 6:1–32. doi: 10.1142/S0219025703001110

CrossRef Full Text | Google Scholar

22. Ionescu C, Lopes A, Copot D, Machado J, Bates J. The role of fractional calculus in modelling biological phenomena: a review. Commun Nonlinear Sci Numer Simul. (2017) 51:141–59. doi: 10.1016/j.cnsns.2017.04.001

CrossRef Full Text | Google Scholar

23. Magin RL. Fractional calculus models of complex dynamics in biological tissues. Comput Math Appl. (2010) 59:1586–93. doi: 10.1016/j.camwa.2009.08.039

CrossRef Full Text | Google Scholar

24. Magin RL. Fractional Calculus in Bioengineering. Begell House Connecticut (2006).

Google Scholar

25. Kelly JF, Meerschaert MM. Ch. 6: The fractional advection–dispersion equation for contaminant transport. In: Tarasov VE, editor. Handbook of Fractional Calculus with Applications, Vol. 5. De Gruyter (2019).

Google Scholar

26. Machado JT, Lopes AM. Fractional–order modeling of a diode, Commun Nonlinear Sci Numer Simul. (2019) 70:343–53. doi: 10.1016/j.cnsns.2018.11.008

CrossRef Full Text | Google Scholar

27. Meng R, Yin D, Drapaca CS. Variable–order fractional description of compression deformation of amorphous glassy polymers. Comput Mech. (2019) 64:163–71. doi: 10.1007/s00466-018-1663-9

CrossRef Full Text | Google Scholar

28. Hristov JY. Linear viscoelastic responses: the Prony decomposition naturally leads into the Caputo–Fabrizio fractional operator. Front Phys. (2018) 6:135. doi: 10.3389/fphy.2018.00135

CrossRef Full Text | Google Scholar

29. Qureshi S, Yusuf A, Shaikh AA, Inc M, Baleanu D. Fractional modeling of blood ethanol concentration system with real data application. Chaos. (2019) 29:013143. doi: 10.1063/1.5082907

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Pinto CM, Carvalho AR. Diabetes mellitus and TB co-existence: Clinical implications from a fractional order modelling. Appl Math Model. (2019) 68:219–43. doi: 10.1016/j.apm.2018.11.029

CrossRef Full Text | Google Scholar

31. Akman Yıldız T. A fractional dynamical model for honeybee colony population. Int J Biomath. (2018) 11:29–48. doi: 10.1142/S1793524518500638

CrossRef Full Text | Google Scholar

32. Hristov J. On the Atangana–Baleanu derivative and its relation to the fading memory concept: the diffusion equation formulation In: Jos Francisco Gmez LT, Escobar R, editors. Fractional Derivatives with Mittag-Leffler Kernel, Vol 194, Studies in systems, Decision and Control, Springer (2019). p. 175–93.

Google Scholar

33. Peng L, Debbouche A, Zhou Y. Existence and approximations of solutions for time–fractional Navier–Stokes equations. Math Methods Appl Sci. (2018) 41:8973–84. doi: 10.1002/mma.4779

CrossRef Full Text | Google Scholar

34. Chang A, Sun H, Zhang Y, Zheng C, Min F. Spatial fractional Darcy's law to quantify fluid flow in natural reservoirs. Phys A Stat Mech Appl. (2019) 519:119–26. doi: 10.1016/j.physa.2018.11.040

CrossRef Full Text | Google Scholar

35. Diethelm K. A fractional calculus based model for the simulation of an outbreak of dengue fever. Nonlinear Dyn. (2013) 71:613–9. doi: 10.1007/s11071-012-0475-2

CrossRef Full Text | Google Scholar

36. Qureshi S, Atangana A. Mathematical analysis of dengue fever outbreak by novel fractional operators with field data. Phys A Stat Mech Appl. (2019) 526:121127. doi: 10.1016/j.physa.2019.121127

CrossRef Full Text | Google Scholar

37. Area I, Batarfi H, Losada J, Nieto JJ, Shammakh W, Torres Á. On a fractional order Ebola epidemic model. Adv Diff Equat. (2015) 2015:278. doi: 10.1186/s13662-015-0613-5

CrossRef Full Text | Google Scholar

38. Tejado I, Pérez E, Valério D. Fractional calculus in economic growth modelling of the group of seven. Fract Calculus Appl Anal. (2019) 22:139–57. doi: 10.1515/fca-2019-0009

CrossRef Full Text | Google Scholar

39. Toledo-Hernandez R, Rico-Ramirez V, Iglesias-Silva GA, Diwekar UM. A fractional calculus approach to the dynamic optimization of biological reactive systems. Part I: fractional models for biological reactions. Chem Eng Sci. (2014) 117:217–28. doi: 10.1016/j.ces.2014.06.034

CrossRef Full Text | Google Scholar

40. Toledo-Hernandez R, Rico-Ramirez V, Rico-Martinez R, Hernandez-Castro S, Diwekar UM. A fractional calculus approach to the dynamic optimization of biological reactive systems. Part II: numerical solution of fractional optimal control problems. Chem Eng Sci. (2014) 117:239–47. doi: 10.1016/j.ces.2014.06.033

CrossRef Full Text | Google Scholar

41. Akman Yıldız T, Arshad S, Baleanu D. New observations on optimal cancer treatments for a fractional tumor growth model with and without singular kernel. Chaos Solitons Fract.(2018) 117:226–39. doi: 10.1016/j.chaos.2018.10.029

CrossRef Full Text | Google Scholar

42. Razminia A, Asadizadehshiraz M, Shaker HR. Optimal trajectory tracking solution: Fractional order viewpoint. J Franklin Inst. (2019) 356:1590–603. doi: 10.1016/j.jfranklin.2018.11.024

CrossRef Full Text | Google Scholar

43. Rui X, Yin W, Dong Y, Lin L, Wu X. Fractional–order sliding mode control for hybrid drive wind power generation system with disturbances in the grid. Wind Energy. (2019) 22:49–64. doi: 10.1002/we.2269

CrossRef Full Text | Google Scholar

44. Rosa S, Torres DF. Optimal control of a fractional order epidemic model with application to human respiratory syncytial virus infection. Chaos Solitons Fract. (2018) 117:142–9. doi: 10.1016/j.chaos.2018.10.021

CrossRef Full Text | Google Scholar

45. Kheiri H, Jafari M. Optimal control of a fractional–order model for the HIV/AIDS epidemic. Int J Biomath. (2018) 11:1850086. doi: 10.1142/S1793524518500869

CrossRef Full Text | Google Scholar

46. Bonyah E, Gomez-Aguilar J, Adu A. Stability analysis and optimal control of a fractional human African trypanosomiasis model. Chaos Solitons Fract. (2018) 117:150–60. doi: 10.1016/j.chaos.2018.10.025

CrossRef Full Text | Google Scholar

47. Zhou X, Mao Y, Ren W. Fractional order robust control under variational working attitude for light tracking system. In: 2018 Chinese Automation Congress (CAC). IEEE (2018) p. 2117–22.

Google Scholar

48. Sharma S, Kumari N. Dynamics of a waterborne pathogen model under the influence of environmental pollution. Appl Math Comput. (2019) 346:219–43. doi: 10.1016/j.amc.2018.10.044

CrossRef Full Text | Google Scholar

49. Shang Y. A lie algebra approach to susceptible–infected–susceptible epidemics. Electron J Diff Equat. (2012) 2012:1–7. Available online at: http:/ejde.math.txstate.edu

Google Scholar

50. Shang Y. Lie algebraic discussion for affinity based information diffusion in social networks. Open Phys. (2017) 15:705–11. doi: 10.1515/phys-2017-0083

CrossRef Full Text | Google Scholar

51. Odibat Z M, Shawagfeh NT. Generalized Taylor's formula. Appl Math Comput. (2007) 186:286–93. doi: 10.1016/j.amc.2006.07.102

CrossRef Full Text | Google Scholar

52. Lin W. Global existence theory and chaos control of fractional differential equations. J Math Anal Appl. (2007) 332:709–26. doi: 10.1016/j.jmaa.2006.10.040

CrossRef Full Text | Google Scholar

53. Van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci. (2002) 180:29–48. doi: 10.1016/S0025-5564(02)00108-6

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Ahmed E, Elgazzar A. On fractional order differential equations model for nonlocal epidemics, Phys A Stat Mech Appl. (2007) 379:607–14. doi: 10.1016/j.physa.2007.01.010

CrossRef Full Text | Google Scholar

55. Ali HM, Pereira FL, Gama SM. A new approach to the pontryagin maximum principle for nonlinear fractional optimal control problems. Math Methods Appl Sci. (2016) 39:3640–9. doi: 10.1002/mma.3811

CrossRef Full Text | Google Scholar

56. Almeida R, Pooseh S, Torres DF. Computational Methods in the Fractional Calculus of Variations. London: Imperial College Press (2015).

Google Scholar

57. Kamocki R. Pontryagin maximum principle for fractional ordinary optimal control problems, Math Methods Appl Sci. (2014) 37:1668–86. doi: 10.1002/mma.2928

CrossRef Full Text | Google Scholar

58. Sweilam N, AL-Mekhlafi S. On the optimal control for fractional multi-strain TB model, Optim Control Appl Methods. (2016) 37:1355–74. doi: 10.1002/oca.2247

CrossRef Full Text | Google Scholar

59. Lenhart S, Workman JT. Optimal Control Applied to Biological Models. Boca Raton, FL: Chapman & Hall; CRC Press (2007).

Google Scholar

60. McAsey M, Mou L, Han W. Convergence of the forward-backward sweep method in optimal control. Comput Optim Appl. (2012) 53:207–26. doi: 10.1007/s10589-011-9454-7

CrossRef Full Text | Google Scholar

61. Lin Y, Xu C. Finite difference/spectral approximations for the time-fractional diffusion equation. J Comput Phys. (2007) 225:1533–52. doi: 10.1016/j.jcp.2007.02.001

CrossRef Full Text | Google Scholar