Numerical treatments for the optimal control of two types variable-order COVID-19 model.

In this paper, a novel variable-order COVID-19 model with modified parameters is presented. The variable-order fractional derivatives are defined in the Caputo sense. Two types of variable order Caputo definitions are presented here. The basic reproduction number of the model is derived. Properties of the proposed model are studied analytically and numerically. The suggested optimal control model is studied using two numerical methods. These methods are non-standard generalized fourth-order Runge-Kutta method and the non-standard generalized fifth-order Runge-Kutta technique. Furthermore, the stability of the proposed methods are studied. To demonstrate the methodologies' simplicity and effectiveness, numerical test examples and comparisons with real data for Egypt and Italy are shown.

Introduction

Coronaviruses are a large family of viruses that are known to cause illness ranging from the common cold to more severe diseases such as severe acute respiratory syndrome (SARS). The World Health Organization (WHO) has described this variant as a variant of serious concern. The pandemic started in late December 2019 with patients admitted to hospitals with an initial diagnosis of pneumonia in Wuhan, the capital of Hubei, China. In the months that followed, the World Health Organization recorded 1,133,758 total cases and 62,784 fatalities globally in a report dated 5 April 2020. Italy has been severely affected. The first Italy’s COVID-19 patient was detected on 21 February 2020 in a little town near Milan, in the northern region of Lombardy, for more details see [1].

Recently, scientific researches focus on the study the coronavirus infection in different scenarios, see for example [2], [3], [4].

Also, optimal control theory has successful applications in biological and medical problems. Furthermore, fractional optimal control problems are a subclass of optimal control whose dynamics are described by fractional differential equations. We can minimize the impact of COVID-19 pandemics by minimizing the number of detected asymptomatic infected people who acquire life threatening symptoms and the threatened percentage populations that become extinct. Fractional derivatives are non-local in nature but the integer order derivative can be used to characterize the short memory in time-dependent systems, whereas the variable-order fractional derivative (VOFD) is known for explaining the impacts of long variable memory [5], [6], [7], [8], [9], [10], [11], [12], [13] in time-dependent systems. Sweilam et al. recently published several research papers in filed of variable-order optimal control problems (VOCPs), for more details see [14], [15], [16].

In this article, we will study numerically the optimal control for two types of variable-order COVID-19 epidemic models. We consider the available data of daily confirmed cases in Italy, from of February to of April and in Egypt from March to June,  [17]. The variable-order fractional derivative is defined in the Caputo sense. Two methods will be used to study the suggested model. Non-standard generalized Runge–Kutta of fourth-order method (NGRK4M) and the non-standard generalized Runge–Kutta of fifth-order method (NGRK5M) are constructed. Comparative studies between these methods are presented. Comparisons with real data for Egypt and Italy are implemented.

The remainder of this paper is organized as follows: Two types of variable-order fractional and some definitions are given in Section “Basic Notations”. In Section “COVID-19 Variable Order Mathematical Model”, the COVID-19 model of variable-order with required analysis of equilibrium points is given. Formulation of variable-order optimal control problem is introduced in Section “Formulation of the Control Problem”. Two numerical approaches are introduced for solving the suggested model in Section “Numerical Methods”. We simulated numerically the proposed model in Section “Numerical Simulations”. In Section “Conclusions”, the conclusions are presented .

Basic notations

In the following problem:

We introduce the Caputo variable-order fractional derivative, where , , as follows [8]:

The right-left side Caputo’s derivatives of order type one (CVOT1) where is continuous function are defined respectively as follows [18]:

The right-left side Caputo’s derivatives of order type two (CVOT2), where is continuous function are defined respectively as follows [18]:

COVID-19 variable order mathematical model

In the following, the COVID-19 model described in [2] will be developed here to variable-order fractional model. The model’s variables and the parameters are listed in Table 1, Table 2 respectively. In the following, the updated nonlinear variable-order mathematical model: Let, and is the total population. The initial conditions are given as follows:

Table 1

Variables of the proposed model [2].

Variables	Description
S(t)	The susceptible population.
I(t)	The infected population with asymptomatic infected and undetected.
I_D(t)	The diagnosed population with asymptomatic infected and detected.
I_A(t)	The ailing population with symptomatic infected and undetected.
I_R(t)	The recognized population with symptomatic infected and detected.
I_T(t)	The threatened population which infected with life-threatening symptoms and detected.
I_H(t)	The recovered population.
I_E(t)	The extinct population.

Table 2

All of the parameters in the proposed model, as well as their meanings [2].

Symbol	Definition	Values
t	Time t≥0.
β^α(t)	The transmission rate between a susceptible and an infected population	(0.57)^α(t)
βDα(t)	The transmission rate between a susceptible and diagnosed population.	(0.0114)^α(t)
γ^α(t)	The transmission rate between a susceptible and ailing population.	(0.456)^α(t)
δ^α(t)	The transmission rate between a susceptible and a recognized population.	(0.011)^α(t)
ɛ^α(t)	The detection rate, with respect to (w. r. t) relative to asymptomatic cases.	(0.171)^α(t)
ɛAα(t)	The detection rate w. r. t relative to symptomatic cases.	(0.3705)^α(t)
ζ^α(t)	The rate which infected is not mindful of being infected creates clinically pertinent side effects and is comparable within the nonappearance of a particular treatment.	(0.1254)^α(t)
η^α(t)	The rate which infected mindful of being infected, creates clinically pertinent indications, and is comparable within the nonappearance of a particular treatment.	(0.1254)^α(t)
μ^α(t)	If there is no known specific medicine that is effective against the illness, the rate at which infected persons are unaware creates life-threatening symptoms.	(0.0171)^α(t)
ν^α(t)	the rate at which infected aware develop life-threatening symptoms; they are comparable if there is no known specific treatment that is effective against the disease.	(0.0274)^α(t)
ρIα(t)	The recovered infected population rate.	(0.0342)^α(t)
κ^α(t)	The recovered ailing population rate.	(0.0171)^α(t)
σ^α(t)	The recovered threatened population rate.	(0.017)^α(t)
ξ^α(t)	The recovered recognized population rate.	(0.017)^α(t)
ρDα(t)	The recovered diagnosed population rate.	(0.0342)^α(t)
d^α(t)	The death rate.	(0.01)^α(t)

The basic reproduction number

To describe the fundamental of reproductive number (), for more information, see [19]. Consider the matrices and where, the matrix representing the new infection terms and representing the remaining transfer terms associated with the baseline model are given, respectively, by: where is the individuals transfer rate into or out infected class  from subsystem . The order of infected variables and and is the appearance rate of new infected individuals in class ,

Then , where is the spectral radius of the matrix . Also, and

At , where is the disease-free equilibrium point. The basic reproduction number for the system (4)–(11) is [20]: where , , and .

The effective reproductive number () is the average number of secondary cases of infected population produced in completely susceptible population. The effective reproduction number can be estimated by the product of the basic reproductive number and the susceptible population . So: We have that stability of the system occurs if: , and the system becomes unstable if  .

Equilibrium points

The aim now is to determine the equilibrium points for (4)–(11) and explore their local asymptotic behavior:

. Then, from (1), we have where, denotes any equilibrium point.

Disease-free equilibrium

If , and the right-hand side of the system (4)–(11) are both equal to zero, we get the disease-free point as follows:

Endemic equilibrium

If the infected variables are not zero, the subsystem (6)–(10) has an endemic equilibrium, [19]. If the right-hand side of the subsystem (6)– (11) is equal to zero, then the endemic equilibrium point is given as follows:

From Eq. (6), from equations , we have By substitution in (12) From Eq. (10) and .

Formulation of the control problem

We will extend the necessary conditions of Pontryagin’s maximum principle to variable-order fractional differentiation equations. The COVID-19 infection can be controlled and minimized in a community by maximizing the number of healed people and minimizing the number infected individuals. For this purpose, we use three suitable control variables depending on time for the reduction of diagnosed asymptomatic infected; for the reduction of developing life-threatening symptoms and for the percentage of threatened population being extinct. Let, be the cost (objective) functional which will be minimize and and are the weighting constants of recognized symptomatic infected and threatened acutely symptomatic infected individuals respectively, whereas and are the positive weight constants used for the treatment of diagnosed asymptomatic infected, recognized symptomatic infected and threatened acutely symptomatic infected individuals, respectively. subjected to the constraint (4)–(11) as follows: and satisfying the initial conditions: The Hamiltonian is given as following: where is the co-state variables and right hand side of system (19). From (18), (20), we can derive the necessary conditions as follows: and it is also required that: , where ,; are the Lagrange multipliers. Eqs. (22), (23) describe the necessary conditions in terms of a Hamiltonian for the optimal control problem defined above.

If and are optimal controls with corresponding state ; then there are adjoint variables , , :

From (22) , the control functions and : and

Now, by substituting , and in (19), we have the state equations as follows:

Numerical methods

The mathematical models (25) and (24) will be studied in this section, moreover, two non-standard methods are constructed in the next sections.

A numerical method is considered non-standard if at least one of the following conditions is satisfied [21]:

1-It is decided to use the non-local approximation.

2- Discretization of a derivative is not traditional and use a nonnegative function, We will give a quick rundown of these methods:

NGRK4M

Let a collection of mesh points , such that , ,   is the final time and .

Now, we can write the approximation solution of Eq. (1) [22]: For more details on the fractional NGRK4M stability analysis, see [10]. We use the proposed method to solve (25) with initial condition. Also, we discretize (24) with transitivity condition using the same technique with backward in time discretization, we get system of algebraic equations which can be solved easily.

NGRK5M

Let a collection of mesh points , such that , and , .

In the proposed GRK5M [22], to approximate solution for Eq. (1), we will substitute the step size in GRK5M, by a function in , and this function is continuous:

The general formula for NGRK5M type 1 is given as follows: We use the proposed method to solve (25) with initial condition. Also, we discretize (24) with transitivity condition using the same technique with backward in time discretization, we get system of algebraic equations which can be solved easily.

Stability analysis of NGRK5M type 1

To study the stability of NSRK5M. Consider , then, we can approximate (1) using NGRK5M type 1 as: In order to study the stability of proposed method (28), the test problem is used as follows: The solution is asymptotically stable if is a negative number and becomes unstable if is positive and increasing, [9]. NGRK5M is applied to (29), the equation becomes: Therefore the stability condition:

The general formula for NGRK5M type 2 is given as follows: and We use the proposed method to solve (25) with initial condition. Also, we discretize (24) with transitivity condition using the same technique with backward in time discretization, we get system of algebraic equations which can be solved easily.

Stability of NGRK5M type 2

In the following, we can approximate (1) using NGRK5M type 2 as: Using (29), (32), we have: Therefore, the stability condition:

The algorithm of NGRK5M

This algorithm consists of the following steps:

Consider and represent the vector approximations for the state and adjoint systems respectively.

Step 1. Make an initial guess for and

Step 2. Solve forward NGRK5M in time, using the initial condition and the values for and .

Step 3. Solve backward NGRK5M in time using the transversality condition and the values for and and .

Step 4. By entering the new and values into the characterization of the optimal control, update and .

Step 5. Output the current values as solutions, if values of the variables in this iteration and the last iteration are negligibly close. Return to Step 2, if values are not close.

Numerical simulations

The solutions of the variable-order fractional optimality systems (24) and (25) are simulated using two numerical methods in this section. These methods are NGRK4M and NGRK5M forward with time to solve the state system (25) with initial condition and backward with time to solve the co-state systems (24) with transversality conditions: . Also, we will use the modified parameters from Table 2. We consider the initial states of the population at day is set as [2]: , and , where Italian population and Egyptian population . Fig. 1, Fig. 2, Fig. 3, Fig. 4, show the comparison between the behavior of approximation solutions and real data [2], using GRK5M and different . The number of diagnosed recovered cases represent as: Also, the number of cumulative diagnosed cases are the sum of the following variables , , , , .

Fig. 1

Real data compared to the approximate solutions using GRK5M, when .

Fig. 2

Approximate solutions behavior compared to real data in Egypt, when .

Fig. 3

Behavior of the approximate solutions compared to real data, when type one using GRK5M.

Fig. 4

Behavior of the approximate solutions compared to real data, when type two using GRK5M.

Table 3, shows the values of the objective functional for different with and without controls using NGRK5M when . Fig. 5, shows the behavior of the solutions the threatened population , the recognized population and the growth rate of recovered population using different methods in case and . We noted that the results obtained by NGRK5M are better than the result obtained by NGRK4M. Fig. 6, illustrates the solution for using NGRK5M with linear and nonlinear when and at three different cases of i.e.,  is constant order and is defined using type 1 of NGRK5M and type 2 of NGRK5M. We noted that the results obtained by type 2 of NGRK5M are convergent to the solution obtained by constant fractional order better than the results obtained by type 1 of NGRK5M. Fig. 7, Fig. 8 show the impact of and on behavior of at different values of . We noted that the model of variable-order is general model than the model of fractional order and integer order, new behavior of the solution is appeared by using different values of .

Table 3

The values of objective functional and using NGRK5M.

α(t)	J values without controls	J values with controls
α(t)=1	20.3132	10.4946
α(t)=1−0.009t	17.5503	9.6116
α(t)=1−0.001(cos^2(t))	20.3113	10.4971
α(t)=1−0.05cos(10t)	19.4945	10.4954
α(t)=0.98	20.0478	10.4552
α(t)=0.98−0.001t	15.4357	8.7975
α(t)=0.85−0.0001t	10.582	6.6575
α(t)=0.99−0.001sin(t)	20.2403	10.4981
α(t)=0.90−0.001sin^2(t)	15.3228	8.7522
α(t)=0.95−0.05cos^2(t)	18.5981	9.9883

Fig. 5

Behavior of the approximate solutions of , and the growth rate of at different methods when .

Fig. 6

Optimization of the approximation solutions of at different types of .

Fig. 7

The impact of and on behavior of at . Case(a) when and case (b) when .

Fig. 8

The impact of and on behavior of at different value of .

Conclusions

The current paper analyzed a novel optimal control variable-order COVID-19 pandemic with modified parameters. Moreover, two variable order definitions in the Caputo’s sense are presented to extended the COVID-19 pandemic model. From the comparison with real data from Egypt and Italy, we concluded that the proposed model is described well the real data in Italy compared the real data in Egypt. Specifically, we had successfully applied a kind of Pontryagin’s maximum principle to reduce the threatened population by using three suitable control variables depending on time. NGRK4M and NGRK5M are constructed to solve the optimality system. We can claim from Fig. 6 that the solutions obtained using CVOT2 is the best. Numerical outcomes are introduced to show the validity and applicability of the proposed scheme. In future, the present study can be extended to piecewise differential model.

CRediT authorship contribution statement

Nasser Sweilam: Conceptualization, Data curation, Formal analysis, Funding acquisition, Resources, Methodology, Investigation, Supervision, Writing – review & editing. Seham Al-Mekhlafi: Conceptualization, Data curation, Formal analysis, Funding acquisition, Resources, Software, Validation, Visualization, Writing – original draft, Methodology, Investigation, Supervision. Salma Shatta: Software, Validation, Visualization, Writing – original draft. Dumitru Baleanu: Supervision, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.