key: cord-0681331-9hetvxmg
authors: Bouhali, Amira; Aribi, Walid Ben; Miled, Slimane Ben; Kebir, Amira
title: Optimal Control applied to SIRD model of COVID 19
date: 2021-09-03
journal: nan
DOI: nan
sha: 3649fc6978b31f21c6bdfecf3be46315db2bed23
doc_id: 681331
cord_uid: 9hetvxmg

In this study, we present an epidemic-controlled SIRD model with two types of control strategies: mask wear and screening. The aim of this study is to minimize the number of deceased keeping a minimal cost of mask advertising and screening. The model is proved to be well-posed and to have an invariant region. Also, a thorough study of the dynamics is effected and the basic reproduction number is used to study the stability of the steady states of the model. As for the optimal control analysis, the existence of an optimal control was checked. Then its characterization was carried out using Pontryagin's minimum principle. Numerical simulations are conducted after that with different values of maximal screening for comparison. The findings of the optimal control analysis and numerical simulations both reveal that the optimal pair of strategies contribute enormously in lowering the number of infected and dead individuals. Although zero infection is not achieved in the population, this study implies that carrying an optimal approach constitutes a major step in controlling the spread of the disease to the barest minimum that can buy time for middle and low-income countries to carry on with their vaccination strategies.

Severe acute respiratory syndrome coronavirus 2 commonly known as SARS-CoV-2 is a novel coronavirus that has caused the global pandemic of COVID-19 first reported in Wuhan China a year ago. Soon after, the virus has spread in a very rapid way to extend to the whole world making the World Health Organization declare a global pandemic on March 11 th , 2020. The virus has proved to be very difficult to contain out of the quarantine measures due to its high contagion that lead to over 180000000 people to be infected worldwide. However its high contagion is not the only problem as it has proven as well to be lethal causing over 3900000 deaths around the world. On the other hand, the economic pressure on the governments has shown how inconvenient the lockdown strategy is on a long term and how much required it is to carry on with a more normal way of life. For that, the problem has been treated not only biologically in the purpose of vaccine and treatment implementation but also mathematically to study its social effects. This is not a first as mathematical modeling has provided a very powerful tool for investigating the dynamics of infectious diseases and controlling them. Previous studies have introduced different models allowing to predict and assess intervention strategies during pandemic spread [1, 2] such as Ebola [3] , Tuberculosis [4] or the current Covid-19 [5, 6] . Some models assumed life-long immunity such as SIR models, others considered the lower limit as a recovered individual is not supposed to be immune for no matter how short a period and they introduced the SIS models. Some also were more realistic and assumed a gained immunity for a period of time revealing the SIRS models. In our case, we consider an SIRD model where we introduce the disease-caused death equation into the model dynamics as our focus, in the second part of this study, is on minimising the number of these deaths. The optimal control efforts serve mainly that end. Thus, the ultimate goal of this study is to minimize the number of deaths with basic strategies only: mask wear and screening at a minimal cost. This presents the possibility of containing the disease without any extreme measures such as lockdown or vaccination which represents a very suitable solution for middle and specially low-income countries. As it allows them to minimize the costs and provides time so that they can carry on with their vaccination strategies. In this work, both mathematical and numerical analysis of a controlled epidemiological model of four sub-populations: susceptible, infectious, recovered and dead are presented. Section 2 is a study of the dynamics of the SIR model, its equilibria and their stability. Section 3 focuses on the optimal control problem that aims to reduce the number of the deceased keeping a minimal screening cost. Section 4 is dedicated to the numerical simulations and the discussion. Then, a conclusion was drawn in the last section.

This section outlines the formulation of a deterministic SIRD model for COVID-19. The total population at time t is divided into four subpopulations: Susceptible, S(t) ; Infectious, I(t) ; Recovered, R(t) and Dead, D(t). Two types of control u 1 (t) and u 2 (t) are used where 1 − u 1 (t) is the probability of mask wear and u 2 (t) is the screening rate. In the Susceptible compartment, S(t), people are recruited into the population at a constant rate, Λ, through migration/birth. They exit this compartment either through infection induced by the disease with the force of infection, u 1 β I(t) or natural mortality. The infectious compartment, I(t), gains population through infection induced by the disease at the rate of u 1 (t) β S(t). A proportion, α, exits this compartment through recovery at a rate u 2 (t) + δ after screening or end of incubation period, the remaining proportion, 1 − α, of the infectious individuals leaves this compartment at a rate u 2 (t) + δ towards the dead compartment through disease induced death, D(t). Recovered individuals are assumed to develop permanent immunity to COVID-19, and compartments, S , I and R are assumed to have a natural mortality rate, µ. Therefore, the epidemic model is given by the following system:

subject to the following initial conditions

All parameters of this model are considered positive .

In what follows, we will study the dynamic of the sub-model, susceptible, infected and recovered (SIR) model, in the case where controls are constants.

The SIR model corresponds to the first three equations of the system (1):

We aim here to understand the impact of time independent control parameters, i.e., u 1 (t) = u 1 and u 2 (t) = u 2 , on the transmission dynamics of the COVID-19.

By the following, we prove that the solutions are uniformly bounded in a positive invariant region,

Theorem 2.1. For any non-negative initial condition, the solution of system (2) remains non-negative and positively bounded. In addition, the set Ω is positively invariant for the epidemic model (2).

Proof 1. The positivity of the solutions of the system (2) can be verified by examining the direction of the vector field dS(t) dt , dI(t) dt , dR(t) dt T of (2) on each coordinate plane. In the {IR} hyper-plane, for S = 0, one has,

This shows that the vector field points to the interior of R 3 + . Therefore, no trajectory can leave the positive octant by crossing the boundary face S = 0. Solutions starting from the {IR} hyper plane remain in the same hyper plane. Similarly, if R = 0 and I > 0, one has dR(t) dt | R=0 = α (u 2 + δ) I(t) > 0. Thus, no trajectory can leave through the boundary face R = 0. Likewise, if I = 0 and R > 0 one has dI(t) dt | I=0 = 0. This indicates that once a trajectory enters this boundary face, it will remain there. In addition, we have:

In addition, let N (t) = S(t) + I(t) + R(t) be the total population number at time t, then ∀t ∈ R + :

Consequently, according to Gronwall's lemma, one has

and then,

Existence and global stability of equilibrium points In this section the existence and the stability of disease-free equilibrium and the endemic equilibrium states of model (2) are examined.

First, we need to define the basic reproduction number, R 0 . This quantity predicts the spread of a disease in the population. It is defined as the average number of secondary infections generated when an infected person is introduced into a host population where everyone is susceptible and it is given by :

where F (S, I, R) = u 1 (t)βS(t)I(t) and V (S, I, R) = (u 2 (t) + µ+ δ)I(t) denote respectively the rates of the transfer in and out of the infected compartment. Then, It is easy to show that the system (2) has two steady states: a disease-free equilibrium (DFE) given by E * 0 = ( Λ µ , 0, 0) that exists for any value of the parameters and an endemic equilibrium E * 1 = (S * , I * , R * ) in the interior of Ω that exists if and only if R 0 > 1 and where,

For the global stability of equilibrium we use popular types of Lyapunov functions i.e, the common quadratic and Volterra-type functions.

Theorem 2.2. If R 0 ≤ 1, then the DFE, E * 0 , is globally asymptotically stable on Ω. If R 0 > 1, then the endemic equilibrium, E * 1 , is globally asymptotically stable.

Proof 2. Since R is not present in the first two equations of the system (2), then by theorem 3.1 of [7] , to study the stability, it is sufficient to analyze, the following isolated subsystems:

It is obvious that equation (6) is globally exponentially stable for I = 0 or I = I * .

As for the global asymptotic stability of (5), one can use the Lyapunov function and LaSalle's theorem [8, 9] .

For the DFE, we define by

It follows from the LaSalle's invariance principle [9] that ( Λ µ , 0) is globally asymptotical stable for (5). Therefore, using Theorem3.1 of [7] , we conclude that E * 0 is an asymptotically stable equilibrium point of (2).

If now R 0 > 1, the endemic equilibrium exists, and then we define the Lyapunov function L : {(S, I) ∈ Ω 1 : S > 0, I > 0} → R by:

L is C 1 on the interior of Ω 1 , E * is the global minimum of L on Ω 1 , and L(S * , I * ) = 0.

The time derivative of L computed along solutions of (2) iṡ

Knowing that, Λ = u 1 βS * I * + µS * , then we havė

ThenL is negative definite andL = 0 if and only if S = S * . Therefore the largest compact invariant set in (S, I) ∈ Ω 1 :L = 0 is the singleton {(S * , I * )}. Indeed if S = S * we getṠ = 0 and then I = I * . By LaSalle's invariance principle, we conclude that (S * , I * ) is globally asymptotically stable for the system (5). Therefore, using Theorem3.1 of [7] , we conclude that E * 1 is an asymptotically stable equilibrium point of (2).

In this section, we investigate the effect of mask wear and screening on the spread of the disease. The two strategies are still assumed constant.

First it is clear that for a 100% of mask wear i.e. u 1 = 0, one has R 0 = 0 and thus there is no transmission of the disease. Screening at this point is not needed . On the other hand, for 0% of mask wear i.e. u 1 = 1 one has R 0 = βΛ µ(u 2 +µ+ δ) . And since stopping the spread of the disease requires that R 0 < 1, it becomes necessary for the screening value to exceed Λ β µ − (δ + µ) in order to keep the spread of the disease under control. To investigate the effect of a combination of mask wear and screening, R 0 is plotted for different values of mask wear 1 − u 1 ∈ {0; 0.2; 0.4; 0.6; 0.9}. The curve of R 0 is decreasing and it crosses the value of one at some point. We noticed that the more people wear masks, the less screening is needed and the faster R 0 reaches the value 1, see figure (1). Hence, it is obvious that wearing masks, even without screening, helps in flattening the infection curve and increasing the number of susceptible, see figure (2) . 

As optimal control techniques are of great use in developing optimal strategies to control various kinds of diseases we use them in this section with the aim of reducing the number of deceased individuals at a finite time, D(t f ), keeping a minimal cost of screening t f 0 u 2 2 (t)dt. Therefore, the objective function that we seek to minimize over a finite time horizon [0, t f ] is given by:

Where the set of admissible controls U is given by

Note that the controls are no longer considered constant.

Theorem 3.1. There exists an optimal control u * and a corresponding state variables vector (S 0 , I 0 , R 0 , D 0 ) that minimizes the objective function.

Proof 3. The existence of the optimal control pair can be obtained using a result by Fleming and Rishel (1975) and by Lukes (1982) [1, 10] .

In fact, one can easily verify that:

1. The set of controls and corresponding state variables is nonempty. 2. The admissible set U is convex and closed.

3. The right hand side of the state system 1 is bounded by a linear function in the state and control variables. 4. The integrand of the objective functional L is convex on U and there exists constants ω 1 > 0, ω 2 > 0 and ρ > 1 such that

In order to determine the optimal control, Pontryagin's Minimum Principle was used [10] . The latter changes the optimality system into a study of the Hamiltonian variations through the use of adjoint functions . The Hamiltonian is given by

is the vector of state variables and λ = (λ 1 (t), λ 2 (t), λ 3 (t), λ 4 (t)) is the vector of adjoint variables and < ., . > is the scalar product.

Theorem 3.2. Given optimal controls u * 1 (t) , u * 2 (t) and the corresponding solution S 0 (t) , I 0 (t) , R 0 (t) and D 0 (t) of the corresponding state system (1) -(8), there exists adjoint variables λ 1 , λ 2 , λ 3 and λ 4 that satisfẏ

with transversality conditions:

Furthermore, the optimal control is given by

Proof 4. According to Pontryagin's minimum principle,

Thus, the adjoint functions (λ 1 (t), λ 2 (t), λ 3 (t), λ 4 (t)) have the following dynamicṡ

with the final conditions λ(t f ) = (0, 0, 0, 0).

From the third and fourth equations we can deduce that λ 3 ≡ 0 and λ 4 ≡ 0.

(11) and the adjoint variables dynamics is reduced tȯ

Also, the Pontryagin's Minimum Principle states that the optimal control u * minimizes the Hamiltonian, hence we should seek the minimum of H . So we need to study the critical points of the Hamiltonian. A critical point of

The equation ∂H ∂u 2 = 0 implies that

The first equation however, shows that the minimum is either reached at

according to the sign of λ 2 − λ 1 . In fact when u 1 is supposed constant; H would depend on u 2 only and therefore u * 2 is a minimum to H since A 2 > 0. In that case, one has

Since u min

then two scenarios are possible

and consequently,

and consequently,

Note that u * 2 must satisfy u min 2 < u * 2 < u max 2 to be taken into consideration. Otherwise,

• min

Since A 2 > 0, we deduce that it is not possible to have ∂H ∂u 2 = 0 and therefore we cannot discuss the case of singular control in the usual terms. However, it is possible to have ∂H ∂u 1 = 0 which implies that βSI(λ 2 −λ 1 ) = 0. Consequently, either S.I = 0 or λ 2 − λ 1 = 0. As the first case does not present quite an interesting case of study, we move to the latter that yields

However, according to the co-state variables dynamics, one hasλ 1 = µλ 1 which implies that λ 1 (t) = λ 1 (t 0 )e µ(t−t 0 ) .Consequently,

. This equality is absurd except for one particular case α = 1 and λ 1 (t 0 ) = 0. Therefore, the existence of an interval [t 0 ; t 1 ] such that ∂H ∂u = 0 ∀t ∈ [t 0 ; t 1 ] is not possible.

In this section, the system (1) is solved numerically, and the results obtained are presented below. The numerical simulations were carried out by implementing a 4th order Runge-Kutta Method (see, for example [10] ). This iterative method consists in solving the system of equation (1) . Details of the application of this method are developed in [11] . Then, the adjoint variables equations are solved by a reverse fourth order Runge-Kutta scheme using the current iteration solution of equation (4) . The iteration stops if the values of the unknowns at the previous iteration are very close to those of the current iteration. The parameters used are presented in the table 1.

References α

The rate at which infected individuals become cured ≈ 0.99

The total size of the population 11172177

The disease transmission coefficient 0.15473652/N (0) Fitted 1/δ

The mean duration of infection 10.14 days Fitted µ

The death rate 0.000017534

The birth rate 510.5937

The balancing factor associated to the cost component 30

The balancing factor associated to the cost component 10

Mask wear rate per unit of time 0.4 < u 1 < 1 Assumed u 2

Screening rate per unit of time 0 < u 2 < 0.5 Assumed To start, the system is solved using the set of parameters listed above and the following initial conditions We introduced the control and solved the optimality system. With the use of these parameters, and the adjoint variables dynamics, the following solutions for λ 1 and λ 2 were obtained. For this set of parameters, λ 2 − λ 1 is always positive (see figure 3 ). According to the optimal control study conducted above, this results in

For that value of u 1 , one has maximal constant mask wear while the screening rate starting at a value near 0.5 remains constant during the first 50 days then starts decreasing until it reaches 0 (see figure 4 ). Then, the state variables were plotted in two cases: controlled and uncontrolled. In the absence of any form of control, the susceptible curve starts decreasing at the 50 th day until it reaches a value near zero. On the other hand, the curve of the infected reaches a peak that exceeds 2. 10 6 and the number of deaths reaches 3. 10 5 . However, once the system is controlled, a huge difference in the dynamics is observed. The susceptible number is increasing as opposed to the the infections that start at a maximal value of 10 3 then decrease to a value near 100. The dead curve is still increasing; however, to a maximal value less than 50 ( see figure 5 ).

The coefficients, A 1 and A 2 , are balancing cost factors due to the size and importance of the two parts of the objective function. We assume that the coefficient A 1 associated with D(t f ) is greater than or equal to A 2 , which is associated with the control u 2 . This assumption is based on the following facts: The cost associated with D(t f ) will include the cost of the dead person, and the cost associated with u 2 will include the cost of screening.The fractions of the weighting factors, A 1 /A 2 = 1, 3, 10 and 100, are presented in Figure  (6 ). And to illustrate the optimal strategy we have chosen the weighing factor, A1/A2 = 3 since the only change observed was in the values of the controls rather than the dynamics. In order to present the importance of maximal mask wear and screening values, the state variables were represented for two values of maximal screening u max 2 ∈ {0.2; 0.5} and four different values of mask wear 0.2, 0.4, 0.6, 0.9 (See figure 8) . In the first case, with maximal screening at 0.2, and for all values of mask wear, the susceptible and recovered were increasing and the dead didn't exceed 60. However, the infected curve soon showed an increase after having decreased for low mask wear 1 − u 1 ∈ {0.2, 0.4}. For 40% of mask wear, the number of infections decreases to 400 and remains constant after day 50. For 90% of mask wear, the infected keeps decreasing. On the other hand, with u max 2 = 0.5, the susceptible and recovered curves reach their maximal values much faster. The dead curve has a maximal value lower than the first case. And interestingly, the curve of the infected does not show the increase registered in the first case. It also stabilizes at a very low value that is almost equal to zero.

It is important to mention that the values used for the transmission rate β and the mean duration of infection 1 δ are not assumed but are rather estimated based on the real values registered by the Tunisian authorities as presented bellow.

The root mean square error (RMSE) [14] is a frequently used method to measure the difference between the values predicted by a model and the values observed in reality. Let X obs be the vector of the observed values and X model the vector of modeled ones. The RMSE of a prediction model with respect to the estimated variable X model is defined as follows

Hence, to obtain optimal parameters {β, δ} for our model, one should solve the following problem : min RM SE

Here, the fit is measured by computing the value of the RMSE function using data of deaths for the beginning of the second wave in Tunisia which is calibrated from September 2021, provided by 1 as X obs data. X model is the death data obtained by the SIRD model (1) subject to the following initial condition In addition, to minimize the RMSE function, we used the genetic method 2 to update the parameters β and δ. Figure 9 shows the result of the fitted values using the optimal parameters β and δ.

In this work, a study of COVID-19 transmission for the case of Tunisia was carried out. A four compartmental mathematical model with mask wear and screening as time-dependent control measures is developed. The model is proven to have an invariant region where it is well-posed and makes biological sense to be studied for human population. Different properties of (a) Infected (b) Death Fig. 9 : The fitted value using the optimal parameters β = 1.3201448967113115e − 08 and δ = 0.098575 the model including global stability analysis of the disease-free and endemic equilibrium points have been studied. Some of the parameter estimates were taken from literature and the remaining parameters were computed based on real daily data of COVID-19 confirmed cases of Tunisia. The basic reproduction number was also checked. An optimal analysis of the model for the purpose of assessing the effect of screening companions was conducted. The result showed that the optimal practice of combination of these two strategies significantly reduces the number of infections and deaths(see figure 8 ). In fact, the usage of any form of prevention strategy such as personal mask wear alone led to a decrease in the number of cases in the infected and dead compartments. The screening applied alone also contributed in lowering the number of infections by lowering the basic reproductive number. It is also found that combining control strategies mask wear and screening is even better at combating the deadly COVID-19 pandemic in Tunisia . The optimal application of the control measures (mask wear and screening) though led to much better and faster results. However, for quicker results, governments are required to set higher maximal values of screening and mask wear (see figure 8 ).

Optimal control of a delayed sirs epidemic model with vaccination and treatment

Optimal control of an epidemic model with a saturated incidence rate

Modeling the spread of ebola with seir and optimal control

Optimal control of treatment in a two-strain tuberculosis model, Discrete and continuous dynamical systems-seriesB

Optimal control applied to a seir model of 2019-ncov with social distancing

Optimal quarantine strategies for covid-19 control models

Decomposition techniques for large-scale systems with nonadditive interactions : Stability and stabilizability

The stability of Dynamical Systems

Some extensions of lyapunov's second method

Numerical optimal control of hiv transmission in octave/matlab

WHO), Novel coronavirus (2019-ncov) situation reports

Indicateurs du world-factbook [archive] publié par la cia

Bayesian estimation and tracking: a practical guide

This work was supported in part by the French Ministry for Europe and Foreign Affairs via the project "REPAIR Covid-19-Africa" coordinated by the Pasteur International Network association and by European Union's Horizon 2020 research and innovation program under grant agreement No. 883441 (STAMINA).