On the optimal control of coronavirus (2019-nCov) mathematical model; a numerical approach

In this paper, a novel coronavirus (2019-nCov) mathematical model with modified parameters is presented. This model consists of six nonlinear fractional order differential equations. Optimal control of the suggested model is the main objective of this work. Two control variables are presented in this model to minimize the population number of infected and asymptotically infected people. Necessary optimality conditions are derived. The Grünwald–Letnikov nonstandard weighted average finite difference method is constructed for simulating the proposed optimal control system. The stability of the proposed method is proved. In order to validate the theoretical results, numerical simulations and comparative studies are given.


Introduction
The well-known coronavirus disease (COVID-19) pandemic can be consider as one of the serious pandemic diseases all over the world, for more details, see [1,2]. The spread of this disease has serious impact on human society and health. The modeling study of infectious diseases is very useful in making strategies to control this disease. Recently, many interesting papers on modeling the coronavirus have appeared, see for example [3][4][5][6][7].
In general, mathematical models involving the known ordinary differentiation could be used to capture dynamical systems of infectious disease, when only initial conditions are used to predict future behaviors of the spread. However, when the situation is unpredictable, which is the case of COVID-19, due to uncertainties associated with the pandemic, ordinal derivatives and their associated integral operators show deficiency. The fractional order differential equation (FODE) models seem more consistent with the real phenomena than the integer order models [8][9][10][11][12][13].
Moreover, one of the new topics in mathematics is the fractional optimal control (FOC). FOC can be defined using a variety of fractional derivative definitions. Riemann-Liouville and Caputo fractional derivatives [14][15][16][17][18] can be considered the most important fractional derivative definitions. Sweilam and Al-Mekhlafi introduced some numerical studies for FOC, for more details, see [19][20][21].
The main contribution of this work is to develop a numerical scheme to provide approximate solutions for the fractional optimal control problems (FOCPs). We consider the mathematical model in [4] with modified parameters. The fractional order derivatives are defined here in the Caputo sense. Moreover, we introduce two control variables, u p (t) and u ap (t), in order to minimize the number of the infected and the asymptotically infected. The Grünwald-Letnikov nonstandard weighted average finite difference method (GL-NWAFDM) is established to simulate the obtained fractional order system.
To the best of our knowledge, the fractional optimal control for coronavirus (2019-nCoV) mathematical model with GL-NWAFDM has never been explored.
This paper is organized as follows: The basic mathematical formulas are introduced in Sect. 2. The proposed model with FO and two controls is presented in Sect. 3. In Sect. 4, the formulation of the optimal control problem and the necessary optimality conditions are derived. In Sect. 5, the new numerical method GL-NWAFDM and the stability analysis are introduced. Numerical simulations are discussed in Sect. 6. Finally, the conclusions are presented in Sect. 7.

Preliminaries and notations
In this section, we recall some important definitions of fractional calculus used throughout the remaining sections of this paper. The fractional order derivative in the Caputo sense can be defined as follows [22]: where is the Euler gamma function and 0 < α < 1.

Fractional order model of coronavirus with control
Herein, we consider the new mathematical model of coronavirus given in [4] with modified parameters. Two controls, u p , u ap , are added to health care such as isolating patients in private health rooms and providing respirators and giving them treatments soothing regularly. This model consists of six nonlinear ordinary differential equations. Moreover, The total population N p = S p + E p + I p + R p Table 2 The parameter values for COVID model [4] Parameter Description Virus removing rate from M (0.01) α Table 1 describes the state variables and Table 2 describes the parameters. It is important to notice that the parameters depend on the fractional order α. To make the system consistent in the physical sense and more consistent with the reality, we must make sure that the right-hand sides of these equations have the same dimensions, for more details, see [15]. Let us assume that N p is a constant. The modified model is then represented by a system of fractional order differential equations: The existence and uniqueness of the solutions of (2) follow from the results given in [25]. The feasible region for model (2) is given by The basic reproduction number of the proposed model (2) is given as follows [4]: where , .
The endemic threshold is given at R 0 = 1, the disease will die out when R 0 < 1, and the endemic case when R 0 < 1, for more details, see [26]. In this work we consider R 0 > 1.

The FOCPs
Consider system (2) in R 6 , let be the admissible control set. We define the objective functional as follows: Now, the goal is to evaluate u p , u ap such that the following functional is minimum, subject to the constraints where and satisfying the initial conditions Consider a modified cost functional as follows [19]: The Hamiltonian is defined as follows: From (7) and (8), we have the FOPC necessary conditions: where Moreover, λ ι (T f ) = 0, ι = 1, 2, 3, . . . , 6.
Remark 1 Under some additional assumptions on the objective functional J and the righthand side of equation (6) must be added, for example, the convexity of J and the linearity of ξ in j and u p , u ap , then the necessary conditions of optimality are also sufficient, for more details, see [27].

GL-NWAFDM
In this section, we construct a novel numerical method called GL-NWAFDM as an extension to the method given in [24,28]. This method can be an explicit method (easy for coding) or an implicit method (more stable and efficient), depending on the weight factor ∈ [0, 1]. To approximate the solutions of system (19) using GL-NWAFDM, we first discretize the Caputo fractional operator (1) with replacing (t) by ϕ( t), where Then the discretization for system (19), where n = 0, 1, 2, 3, . . . , N , using GL-NWAFDM can be written as follows: Let f (t n ) = f n = ζ n be the approximate solution of this equation, then using GL-NWAFDM with (1) we rewrite equation (21) in the following form: Then we have So the proposed implicit scheme is stable.

Numerical simulations
In this section, numerical simulations of the proposed model (19) and (13) with and without optimal control are presented. Gl-NWAFDM is given to obtain the numerical results of the state equations (19) with the following initial conditions [4]: S p (0) = 8,065,518, E p (0) = 200,000, I p (0) = 282, A p (0) = 200, R p (0) = 0, M(0) = 50,000, and we consider N p = 8,266,000. Then, by using the implicit finite difference method [21], we solve (13) with (14) by using different values of 0 < α ≤ 1 with B 1 = B 2 = 500, 1 = 2 = 10, and 0 < ≤ 1. In the numerical simulations the time level is chosen in days, it is up to 10,000. The graphical results obtained through these figures demonstrate that in the case without control, the number of the infected and the asymptotically infected population is increasing, while the number of the population is decreasing in the controlled case as we can see in Fig. 1. This figure demonstrates the effectiveness of two control cases for the proposed model. Moreover, Table 3 reports the values of the objective functional obtained by the proposed method with and without controls and = 1, 0.5. We note that the results obtained in the case of fully implicit at ( = 1) are better than the results in the case = 0.5. Also, the best result of the control case is given at α = 0.8. Figure 2 shows how the behavior of the solutions in the control case is changing by using different values of α and T f = 150. Figure 3 shows the evolution of the approximate solutions for the control variables using different α. It is noted that the range of the solution remains between zero and one. The approximate solutions for S p and M with control case and different values of α are given in Fig. 4. Figure 5 illustrates the behavior of the approximate solutions of E p , I p , A p , R p at different values of α at big time, and it is demonstrated that the proposed method is unconditionally stable at = 1. The results are given by using MATLAB (R2015a).

Conclusions
In the present work, the optimal control of coronavirus model with fractional operator is presented. Also, the combination of fractional order derivative and optimal control in the model improves the dynamics and increases the complexity of the model. Two control variables, u p (t) and u ap (t), are added to health care such as isolating patients in private health rooms and providing respirators and giving them treatments soothing regularly. These have been implemented to minimize the number of the infected and the asymptotically infected as we can see in Fig. 1. Necessary optimality conditions are derived. GL-NWAFDM is constructed to study the behavior of the proposed problem. This method depends on the values of the factor . It can be explicit or implicit with large stability regions as we can see in our results. Moreover, the stability analysis of the proposed method is studied. It was shown that this method has good stability properties in the implicit case. Some simulations are presented to support our theoretical findings. It is concluded that GL-NWAFDM can be applied to solve such fractional optimality systems simply and effectively.