Necessary optimality conditions of a reaction-diffusion SIR model with ABC fractional derivatives

Fractional derivatives give rise to theoretical models that allow a significant improvement in the fitting of real data when compared with analogous classical models [3]. For real data of Florida Department of Health from September 2011 to July 2014, some authors conclude that the absolute error between the solutions obtained statistically and that of fractional models are smaller than those obtained by models of integer derivatives [24]. In the fractional calculus literature, systems using fractional derivatives give a more realistic behavior [26, 25, 23]. There exists many definitions of fractional derivative [23]. Among the more well-known fractional derivatives, we can cite the Riemann–Liouville one. It is not always suitable for modeling physical systems, because the Riemann–Liouville derivative of a constant is not zero, and the initial conditions of associated Cauchy problems are expressed by fractional derivatives. Caputo fractional derivatives offers another alternative, where the derivative of a constant is null and initial conditions are expressed as in the classical case of integer order derivatives [25, 23, 13]. However, the kernel of this derivative has a singularity. Fractional derivatives that possess a non-singular kernel have aroused more interest from the scientific community. This is due to the non-singular memory of the Mittag–Leffler function and also to the non-obedience of the algebraic criteria of associativity and commutativity. The ABC fractional derivative is sometimes preferable for modeling physical dynamical systems, giving a good description of the phenomena of heterogeneity and diffusion at different scales [1, 5, 6].

Fractional calculus plays an important role in many areas of science and engineering. It also finds application in optimal control problems. The principle of mathematical theory of control is to determine a state and a control for a dynamic system during a specified period to optimize a given objective [27]. Fractional optimal control problems have been formulated and studied as fractional problems of the calculus of variations. Some authors have shown that fractional differential equations are more accurate than integer-order differential equations, and that fractional controllers work better than integer-order controllers [7, 20, 21, 28]. In [30], Yuan et al. have studied problems of fractional optimal control via left and right fractional derivatives of Caputo. A numerical technique for the solution of a class of fractional optimal control problems, in terms of both Riemann–Liouville and Caputo fractional derivatives, is presented in [8]. Authors in [9, 11] present a pseudo-state-space fractional optimal control problem formulation. Fixed and free final-time fractional optimal control problems are considered in [10, 12]. Guo [16] formulates a second-order necessary optimality condition for fractional optimal control problems in the sense of Caputo. Optimal control of a fractional-order HIV-immune system, in terms of Caputo fractional derivatives, is discussed in [14]. In [22], authors proposed a fractional-order optimal control model for malaria infection in terms of the Caputo fractional derivative. Optimal control of fractional diffusion equations has also been studied by several authors. For instance, in [2], Agrawal considers two problems, the simplest fractional variation problem and a fractional variational problem in Lagrange form. For both problems, the author developed Euler–Lagrange type necessary conditions, which must be satisfied for the given functional to have an extremum. In [26], authors prove necessary optimality conditions of a nonlocal thermistor problem with ABC fractional time derivatives.

Several infectious diseases confer permanent immunity against reinfection. This type of diseases can be modeled by the $SIR$ model. The total population $(N)$ is divided into three compartments with $N=S+I+R$, where $S$ is the number of susceptible (those able to contract the disease), $I$ is the number of infectious individuals (those capable of transmitting the disease), and $R$ is the number of individuals recovered (those who have recovered and become immune). Vaccines are extremely important and have been proved to be most effective and cost-efficient method of preventing infectious diseases, such as measles, polio, diphtheria, tetanus, pertussis, tuberculosis, etc. The study of fractional calculus with a non-singular kernel is gaining more and more attention. Compared with classical fractional calculus with a singular kernel, non-singular kernel models can describe reality more accurately, which has been shown recently in a variety of fields such as physics, chemistry, biology, economics, control, porous media, aerodynamics and so on. For example, extensive treatment and various applications of fractional calculus with non-singular kernel has been discussed in the works of Atangana and Baleanu [5], and Djida et al. [15]. It has been demonstrated that fractional order differential equations (FODEs) with non-singular kernels give rise to dynamic system models that are more accurately.

In this work, we consider an optimal control problem for the reaction-diffusion SIR system with Atangana–Baleanu fractional derivative in the Caputo sense (the ABC operator). Our aim is to study the effect of non-local memory and vaccination strategies on the cost, needed to control the spread of infectious diseases. Our results generalize to the ABC fractional setting previous studies of classical control theory presented in [17]. The considered model there does not explain the influence of a complete memory of the system. For that, we extend such nonlinear system of first order differential equations to a fractional-order one in the ABC sense. We have further improved the cost and effectiveness of proposed control strategy during a given period of time.

This paper is organized as follows. Some important definitions related to the ABC fractional derivative operator and its properties are presented in Section 2, while the underlying fractional reaction-diffusion SIR mathematical model is formulated in Section 3. This led to the necessity of proving existence and uniqueness of solution to the proposed fractional model as well as existence of an optimal control. These results are extensively discussed in Sections 4 and 5. Section 6 is devoted to necessary optimality conditions. Interesting numerical tests, showing the importance of choosing very well the fractional characterization $\alpha$, are given in Section 7. Finally, conclusions of the present study are widely discussed in Section 8.

We now recall some properties on the Mittag–Leffler function and the definition of ABC fractional time derivative. First, we define the two-parameter Mittag–Leffler function $E_{\alpha,\xi}(z)$, as the family of entire functions of $z$ given by

where $\Gamma(\cdot)$ denotes the Gamma function

Observe that the Mittag–Leffler function is a generalization of the exponential function: $E_{1,1}(z)=e^{z}$. For more information about the definition of fractional derivative in the sense of Atangana–Baleanu, the reader can see [5, 6].

The SIR model is one of the simplest compartmental models. It was first used by Kermack and McKendrick in 1927. It has subsequently been applied to a variety of diseases, especially airborne childhood diseases with lifelong immunity upon recovery, such as measles, mumps, rubella, and pertussis. We assume that the populations are in a spatially homogeneous environment and their densities depend on space, reflecting the spacial spread of the disease. Then, the model will be formulated as a system of reaction-diffusion equations. In this section, we formulate an optimal control of a nonlocal fractional SIR epidemic model with parabolic equations and boundary conditions. We consider that the movement of the population depends on the time $t$ and the space $x$. Furthermore, all susceptible vaccinates are transferred directly to the recovered class.

Let $Q_{T}=[0,T]\times\Omega$ and $\Sigma_{T}=[0,T]\times\partial\Omega$, where $\Omega$ is a fixed and bounded domain in $\mathbb{R}^{2}$ with smooth boundary $\partial\Omega$, and $[0,T]$ is a finite interval. The dynamic of the ABC fractional SIR system with control is given by

with the homogeneous Neumann boundary conditions

and the following initial conditions of the three populations, which are considered positive for biological reasons:

The positive constants $\mu$, $r$ and $d$ are respectively the birth rate, the recovery rate of the infective individuals and the natural death rate. Susceptible individuals acquire infection by the contact with individuals in the class $I$ at a rate $\beta SI$, where $\beta$ is the infection coefficient. Positive constants $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$ denote the diffusion coefficients for the susceptible, infected and recovered individuals. The control $u$ describes the effect of vaccination. It is assumed that vaccination transforms susceptible individuals to recovered ones and confers them immunity. The notation $\Delta=\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}$ represents the usual Laplacian operator in two-dimensional space; $\nu$ is the outward unit normal vector on the boundary with $\frac{\partial}{\partial\nu}=\nu$; and $\nabla$ is the normal derivative on $\partial\Omega$. The no-flux homogeneous Neumann boundary conditions imply that model (4) is self-contained and there is a dynamic across the boundary, but there is no emigration.

Since the vaccination is limited and represents an economic burden, one important issue and goal is to know how much we should spend in vaccination to reduce the number of infections and, at the same time, save the cost of vaccination program. This can be mathematically interpreted by optimizing the following objective functional:

where $\theta$ is a weight constant for the vaccination control $u$, which belongs to the set

of admissible controls. Let $y=(y_{1},y_{2},y_{3})=(S,I,R)$, $y^{0}=(y_{1}^{0},y_{2}^{0},y_{3}^{0})=(S_{0},I_{0},R_{0})$, $\lambda=(\lambda_{1},\lambda_{2},\lambda_{3})$, $\mathrm{L}^{2}(\Omega)=(L^{2}(\Omega))^{3}$ and $A$ be the linear diffusion operator defined by

where

We also set

with

The problem can be rewritten in a compact form as

where ${}^{abc}_{0}D_{t}^{\alpha}$ is the Atangana–Baleanu fractional derivative of order $\alpha\in(0,1)$ in the sense of Caputo with respect to time $t$. The symbol $\triangle$ denotes the Laplacian with respect to the spacial variables, defined on $H^{2}(\Omega)\bigcap H_{0}^{1}(\Omega)$.

Existence of solution is proved in the weak sense.

Integrating by parts, involving the ABC fractional-time derivative (see [15]), and using a straightforward calculation, one obtains the following result.

Using the boundary conditions (5), we immediately get the following Corollary.

We proceed similarly as in [15]. Let $V_{m}$ define a subspace of $H^{1}(\Omega)$ generated by the $w_{1}$, $w_{2}$, $\ldots$, $w_{m}$ space vectors of orthogonal eigenfunctions of the operator $\Delta$. We seek $u_{m}:t\in(0,T]\rightarrow u_{m}(t)\in V_{m}$, solution of the fractional differential equation

To continue the proof of existence, we recall the following auxiliary result.

Since $f(y_{m})\in L^{2}(Q_{T})$, Theorem 4.2 implies that $y_{m}$ is given in a explicit form. The existence of a solution is obtained by using the a priori estimate of Theorem 4.2 and the same arguments used to pass to the limit as those used by us below in the proof of Theorem 5.1.

We prove existence of an optimal control by using minimizing sequences.

In this section, our aim is to obtain optimality conditions. As we shall see, our necessary optimality conditions involve an adjoint system defined by means of the backward Atangana–Baleanu fractional-time derivative.

Let $y^{*}(u^{*})$ be an optimal solution and $u^{\varepsilon}=u^{*}+\varepsilon u\in U_{ad}$ be a control function such that $u\in U_{ad}$ and $\varepsilon>0$. Denote $y^{\varepsilon}=(y_{1}^{\varepsilon},y_{2}^{\varepsilon},y_{3}^{\varepsilon})=(y_{1},y_{2},y_{3})(u^{\varepsilon})$ and $y^{*}=(y_{1}^{*},y_{2}^{*},y_{3}^{*})=(y_{1},y_{2},y_{3})(u^{*})$ the solutions of (4)–(6) corresponding to $u^{\varepsilon}$ and $u^{*}$, respectively. Setting $y^{\varepsilon}=y^{*}+\varepsilon z^{\varepsilon}$ and subtracting the system corresponding to $y^{*}$ from the one corresponding to $y^{\varepsilon}$, we have

System (18) can be rewritten as

associated to Neumann boundary conditions

and initial condition

where $z^{\varepsilon}=(z_{1}^{\varepsilon},z_{2}^{\varepsilon},z_{3}^{\varepsilon})$ and

Set

and

Then, (18) can be reformulated in the following form:

Since the elements of the matrix $F^{\varepsilon}$ are uniformly bounded with respect to $\varepsilon$ and $(-uy_{1}^{*},0,uy_{1}^{*})$ is bounded in $L^{2}(Q_{T})$, it follows by Theorem 4.2 that $z^{\varepsilon}=\frac{y^{\varepsilon}-y^{*}}{\varepsilon}$ is bounded in $L^{\infty}(0,T,L^{2}(\Omega))\bigcap L^{2}(0,T,H^{1}(\Omega))$. Therefore, up to a subsequence of $z^{\varepsilon}$, there exists $z$ such that as $\varepsilon$ tends to zero we have

Put

Note that all the components of the matrix $F^{\varepsilon}$ tend to the corresponding ones of the matrix $F$ in $L^{2}(Q_{T})$ as $\varepsilon\rightarrow 0$. From equations satisfied by $y_{\varepsilon}$ and $y$, we have that