Amira Bouhali
email: amira.bouhali.2893@gmail.com

Slimane Walid Ben Aribi

Ben Miled

Amira Kebir

Amira Kebir Optimal

Slimane Ben Miled

Optimal control for a COVID-19 problem through a SIRD model

Keywords: Optimal control, Structured models, COVID 19, Basic reproduction number

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I INTRODUCTION

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. The virus has proved to be very difficult to contain out of the quarantine measures due to its high contagion and lethalness. On the other hand, the economic pressure on the governments has shown how inconvenient the lockdown strategy is and how much required it is to carry on with a normal way of life. The problem has been treated biologically in the first place by trying to develop vaccine and treatment. However, despite the implementation of several vaccines and their use, many countries kept registering high numbers of deaths and infections. For that, the problem was also considered from a mathematical point of view. 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 [START_REF] Laarabi | Optimal control of an epidemic model with a saturated incidence rate[END_REF][START_REF] Laarabi | Optimal Control of a Delayed SIRS Epidemic Model with Vaccination and Treatment[END_REF] such as Ebola [START_REF] Boujakjian | Modeling the spread of Ebola with SEIR and optimal control[END_REF], Tuberculosis [START_REF] Jung | Optimal control of treatment in a two-strain Tuberculosis model[END_REF] or the current Covid-19 [START_REF] Mallela | Optimal Control applied to a SEIR model of 2019-nCoV with social distancing[END_REF][START_REF] Grigorieva | Optimal quarantine strategies for covid-19 control models[END_REF]. In our case, we consider a 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 study was conducted over an eight-month period as it is the vaccine-induced immunity time interval. During that time, the population is supposed to be completely vaccinated. Thus, the ultimate goal of this study is to minimise the number of deaths among a vaccinated population with two basic strategies: mask wear and screening at a minimal cost. This presents the possibility of containing the disease without any extreme measures such as lockdown. Modeling such a situation represents a very good opportunity as it gives a larger view of the situation and offers the chance of setting a good vaccination schedule that can lead to a total containment of the disease. In this work, both the mathematical and numerical analysis of a controlled epidemiological model of four subpopulations: 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 while 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.

II MODEL DESCRIPTION AND ANALYSIS

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 birth. They exit this compartment either through natural mortality, µ, or infection induced by the disease with the force of infection, u 1 (t) β I(t) where β denotes the disease transmission coefficient. 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 of u 2 (t) + δ after screening or end of the incubation period, the remaining proportion, 1 -α, of the infectious individuals leaves this compartment at a rate of u 2 (t) + δ towards the dead compartment through disease induced death, D(t). Recovered individuals are assumed to develop 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:

         dS(t) dt = Λ -u 1 (t)βS(t)I(t) -µS(t) dI(t) dt = u 1 (t)βS(t)I(t) -(u 2 (t) + µ+ δ)I(t) dR(t) dt = α (u 2 (t) + δ) I(t) -µR(t) dD(t) dt = (1-α) (u 2 (t) + δ) I(t) (1)
subject to the following initial conditions

S(0) ≥ 0, I(0) ≥ 0, R(0) ≥ 0, D(0) ≥ 0
The whole population is assumed to be vaccinated, and all parameters of the model are positive. Moreover, parameter α is less than 1.

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

Analysis of the SIR model with constant controls

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

     dS(t) dt = Λ -u 1 βS(t)I(t) -µS(t) dI(t) dt = u 1 βS(t)I(t) -(u 2 + µ+ δ)I(t) dR(t) dt = α (u 2 + δ) I(t) -µR(t) (2)
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

Ω = (S, I, R) ∈ R 3 + : S + I + R ≤ Λ µ (3)
Theorem 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).

Existence and global stability of equilibrium points

In this section, the existence and the stability of both the disease-free and the endemic equilibria 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:

R 0 = ∂ I F (S, I, R) ∂ I V (S, I, R) | ( Λ µ ,0,0) = u 1 βΛ µ (u 2 + µ+ δ) (4)
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,

S * = Λ µR 0 , I * = Λ u 2 + µ + δ 1 - 1 R 0 , R * = α(u 2 + δ)(R 0 -1) u 1 β .
For the global stability of the endemic equilibrium we use popular types of Lyapunov functions i.e, the common quadratic and Volterra-type functions.

Theorem 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.

III OPTIMAL CONTROL

In this section, we aim to reduce the number of deceased individuals while keeping a minimal cost of screening. The cost of mask wear advertising campaign is treated in two cases: constant in the first part and variable in the second. Note that the controls in this section are no longer considered constant.

Constant mask wear cost

As aforementioned, the objective is to reduce the number of deceased individuals at a finite time, D(t f ), with a minimal cost of screening

t f 0 u 2 2 (t)dt.
The constant cost of mask wear advertising campaign has no effect on the objective function and is assumed equal to zero. Therefore, the objective function that we seek to minimise over a finite time horizon [0, t f ] is given by:

J(u 1 , u 2 ) = A 1 D(t f ) + A 2 t f 0 u 2 2 (t)dt = t f 0 (A 1 (1-α) (u 2 (t)+ δ)I(t) + A 2 u 2 2 (t)) dt (5)
Where the set of admissible controls U is given by:

U = {u = (u 1 , u 2 ) ∈ (L ∞ (0; t f )) 2 | 0 ≤ u min i ≤ u i (t) ≤ u max i ≤ 1 , for i = 1, 2} Theorem 3:
There exist an optimal control u * and a corresponding state variables vector (S 0 , I 0 , R 0 , D 0 ) that minimize the objective function.

Theorem 4:

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) -( 5), there exist adjoint variables λ 1 , λ 2 , λ 3 and λ 4 that satisfy

λ1 = βI(λ 1 -λ 2 )u 1 + λ 1 µ λ2 = βS(λ 1 -λ 2 )u 1 + λ 2 (u 2 + µ + δ) -A 1 (1 -α)(u 2 + δ) λ3 = µλ 3 λ4 = 0 (6)
with transversality conditions:

λ i (t f ) = 0 , i = 1, 2, 3, 4. (7)
Furthermore, the optimal control is given by u

* = (u * 1 , u * 2 )
where

u * 1 = u min 1 , if λ 2 -λ 1 > 0 u max 1 , if λ 2 -λ 1 < 0 u * 2 =      (λ 2 -(1-α)A 1 )I 2A 2 , if u min 2 < (λ 2 -(1-α)A 1 )I 2A 2 < u max 2 u min 2 , if (λ 2 -(1-α)A 1 )I 2A 2 < u min 2 u max 2 , if (λ 2 -(1-α)A 1 )I 2A 2 > u max 2

Variable mask cost

In this section we add the mask advertising campaign cost as a quadratic term -t f 0 A 3 u 1 (t) to the previous objective function. The latter becomes:

J(u 1 , u 2 ) = t f 0 A 1 (1-α) (u 2 (t)+ δ)I(t) + A 2 u 2 2 (t) -A 3 u 2 1 (t) dt (8)
For the same set of admissible controls aforementioned, one has the following results:

Theorem 5:

There exist an optimal control u * and a corresponding state variables vector (S 0 , I 0 , R 0 , D 0 ) that minimize the objective function [START_REF] Laarabi | Optimal Control of a Delayed SIRS Epidemic Model with Vaccination and Treatment[END_REF] if and only if

A 3 A 2 ≤ ( u 2 -v 2 u 1 -v 1
), for all u, v ∈ U .

Theorem 6:

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 exist adjoint variables λ 1 , λ 2 , λ 3 and λ 4 that satisfy

λ1 = βI(λ 1 -λ 2 )u 1 + λ 1 µ λ2 = βS(λ 1 -λ 2 )u 1 + λ 2 (u 2 + µ + δ) -A 1 (1 -α)(u 2 + δ) λ3 = µλ 3 λ4 = 0 (9)
with transversality conditions:

λ i (t f ) = 0 , i = 1, 2, 3, 4. (10)
Furthermore, the optimal control is given by u * = (u * 1 , u * 2 ) where

u * 1 = u min 1 , if u min 1 < βSI(λ 2 -λ 1 ) 2A 3 < 1 and H(u min 1 ) < H(u max 1 ) u max 1 , if 0 < βSI(λ 2 -λ 1 ) 2A 3 < u max 1 and H(u min 1 ) > H(u max 1 ) u * 2 =      (λ 2 -(1-α)A 1 )I 2A 2 , if u min 2 < (λ 2 -(1-α)A 1 )I 2A 2 < u max 2 u min 2 , if (λ 2 -(1-α)A 1 )I 2A 2 < u min 2 u max 2 , if (λ 2 -(1-α)A 1 )I 2A 2 > u max 2
where H is the Hamiltonian of the model.

IV NUMERICAL SIMULATIONS AND DISCUSSION

In this section, system (1) is solved numerically for the constant mask wear case, and the results obtained are presented below. The numerical simulations were carried out by implementing a 4 th order Runge-Kutta Method (see, for example [START_REF] Lenhart | Optimal Control Applied to Biological Models[END_REF]). This iterative method consists of solving the system of equations [START_REF] Lasalle | Some extensions of Lyapunov's second method[END_REF]. Details of the application of this method are developed in [START_REF] Campos | Numerical Optimal Control of HIV Transmission in Octave/MATLAB[END_REF].

The parameters used are presented in table 1. To start, the system is solved using the set of parameters listed above and the following initial conditions S(0) = 11718548; I(0) = 2629; R(0) = 0; D(0) = 0. We introduced 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 1). According to the optimal control study conducted above, this results in u * 1 = u min 1 . For that value of u 1 , one has maximal constant mask wear while the screening rate starting at 0.2 remains almost constant during the first 150 days then decreases to 0 (see figure 2). 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 almost from the start until it reaches a value near zero. On the other hand, the curve of the infected reaches a peak that exceeds 2. 10 5 and the number of deaths reaches 4. 10 4 . However, once the system is controlled, a huge difference in the dynamics is observed. The susceptible number is increasing starting from day 100 as opposed to the infections that start at a maximal value of 3. 10 4 then decrease to 0. The dead curve is still increasing; however, to a maximal value of less than 450 ( see figure 3).

The coefficients, A 1 and A 2 , are balancing cost factors. We assume that A 1 associated with the number of deaths D(t f ) is greater than or equal to A 2 , associated with the screening u 2 . The fractions of the weighing factors, A 1 A 2 = 1, 3, 10 and 100, are presented in Figure [START_REF] Jung | Optimal control of treatment in a two-strain Tuberculosis model[END_REF]. And to illustrate the optimal strategy we have chosen the weighing factor, A 1

A 2 = 3 since the only change observed was in the values of the controls rather than the behaviour. 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 5). In both cases, the same behaviour is observed: infections decrease, both the susceptible and recovered increase and deaths have a maximal threshold. However, a difference in the pace is noticed as infections decrease much faster for high screening values and the maximal value of deaths is lower: less than 450 with high screening as opposed to over 600 for low screening.

V CONCLUSION

In this work, a study of COVID-19 transmission for the case of Tunisia was carried out. A fourcompartmental mathematical model of a vaccinated population with mask wear and screening as time-dependent control measures has been developed. The model is proven to have an invariant region where it is well-posed and makes biological sense to be studied for the human population. Different properties of the model, including global stability analysis of the equilibria, have been studied. Some of the parameter estimates were taken from literature, and the remaining parameters were estimated based on real daily data of COVID-19 confirmed cases in Tunisia. An optimal analysis of the model for the purpose of assessing the effect of mask wear by the individuals and screening campaigns was conducted. The results showed that the optimal practice of a combination of these two strategies in a vaccinated population significantly reduces the number of infections and deaths. And for quicker results, it is required to set higher maximal values of screening and mask wear (see figure 5 ).

A ANNEX 1

Proof of Theorem 3. The existence of the optimal control pair can be obtained using a result from [START_REF] Lenhart | Optimal Control Applied to Biological Models[END_REF][START_REF] Laarabi | Optimal Control of a Delayed SIRS Epidemic Model with Vaccination and Treatment[END_REF]. 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 L(u) ≥ ω 2 + ω 1 (|u 1 | 2 + |u 2 | 2 ) ρ 2 .
Proof of Theorem 4. In order to determine the optimal control, Pontryagin's Minimum Principle was used [START_REF] Lenhart | Optimal Control Applied to Biological Models[END_REF]. The latter changes the optimality system into a study of the Hamiltonian variations through the use of adjoint functions. The Hamiltonian is given by

H(t, u, X, λ) =< λ(t), Ẋ(t) > +A 1 (1 -α) (u 2 (t) + δ) I(t) + A 2 u 2 2 (t)
where X = (S, I, R, D) 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. According to Pontryagin's minimum principle, the adjoint functions (λ 1 (t), λ 2 (t), λ 3 (t), λ 4 (t)) have the following dynamics

λ1 = βI(λ 1 -λ 2 )u 1 + λ 1 µ λ2 = βS(λ 1 -λ 2 )u 1 + λ 2 (u 2 + µ + δ) -λ 3 α(u 2 + δ) -(1 -α)(u 2 + δ)λ 4 -A 1 (1 -α)(u 2 + δ) λ3 = µλ 3 λ4 = 0
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.

Consequently, the Hamiltonian becomes

H = (Λ -u 1 (t)βS(t)I(t) -µS(t)) λ 1 + (u 1 (t)βS(t)I(t) -(u 2 (t) + µ+ δ)I(t)) λ 2 +A 1 (1 α) (u 2 (t)+ δ)I(t) + A 2 u 2 2 (t). (11)
and the adjoint variables dynamics is reduced to

λ1 = βI(λ 1 -λ 2 )u 1 + λ 1 µ λ2 = βS(λ 1 -λ 2 )u 1 + λ 2 (u 2 + µ + δ) -A 1 (1 -α)(u 2 + δ)
B ANNEX 2

Parameter estimation

The root mean square error (RMSE) [START_REF] Haug | Bayesian estimation and tracking: a practical guide[END_REF] 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:

RM SE = 1 n n j=1
(X model,j -X obs,j )

2

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 June 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 S(0) = 1686692, I(0) = 34485, R(0) = 0, D(0) = 0

In addition, to minimize the RMSE function, we used the genetic method 2 to update the parameters β and δ.

Figure 6 shows the result of the fitted values using the optimal parameters β and δ.

Figure 1 :

 1 Figure 1: Adjoint variables λ 1 and λ 2

Figure 2 :Figure 3 :

 23 Figure 2: Optimal mask wear rate 1 -u * 1 (left) and screening rate u * 2 (right) per unit of time

Figure 4 :Figure 5 : 2 = 0 . 2 ( 1 ∈

 452021 Figure 4: The screening rate u * 2 per unit of time

Figure 6 :

 6 Figure 6: The fitted value using the optimal parameters β = 1.3201448967113115e -08 and δ = 0.098575

Table 1 :

 1 Description and values of the parameters.

	Parameters	Description	Values	References
	α	The rate at which infected individuals become cured	≈ 0.99	[12]
	N (0)	The total size of the population	11172177	[14]
	β	The disease transmission coefficient	0.24032955/N (0)	Fitted
	1/δ	The mean duration of infection	5.073 days	Fitted
	µ	The death rate	0.000017534	[14]
	Λ	The birth rate	510.5937	[14]
	A 1	The balancing factor associated with cost component	30	Assumed
	A 2	The balancing factor associated with cost component	10	Assumed
	1 -u 1	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.2	Assumed

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. In fact, when u 1 is considered 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

Assume now that there exists a subset

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

However, according to the co-state variables dynamics, one has λ1 = µλ 1 which implies that λ 1 (t) = λ 1 (t 0 )e µ(t-t0) . Consequently, λ 1 (t 0 )e µ(t-t0) = A 1 (1 -α) , ∀t ∈ [t 0 , t 1 ]. 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.

Proof of Theorem 5. The two first conditions for the existence of optimal control are checked in the first case and remain unchanged.

Moreover, one has

where

Proof of Theorem 6. Using Pontryagin's minimum principle, one has the following dynamics for the adjoint functions (λ

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

Similarly to the first case, one can deduce from the third and fourth equations that λ 3 (t) = 0 and λ 4 (t) = 0. Consequently, the Hamiltonian becomes

and the adjoint variables dynamics are reduced to

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 H,

The first equation, shows that the critical value is u crit 1 = βSI(λ2-λ1)

2A3

. However, it is easy to verify that this value is a maximum to the parabola H(u 2 ). In fact, one has three possible cases • If u crit