OPTIMAL CONTROL AND FREE OPTIMAL TIME PROBLEM FOR A COVID-19 MODEL WITH SATURATED VACCINATION FUNCTION

In this paper, we investigate a free terminal time optimal control applied to 6 ordinary differential equations which describe the spread of COVID-19 infection We propose an extension of the classical Susceptible-Exposed-Infectious-Recovered (SEIR) model, where the infectious patients are divided into unreported (U) and reported cases (I) To have a more realistic model, we estimate the parameters of our model using real Moroccan data We use Bootstrap as a statistical method to improve the reliability of the parameters estimates The main goal of this work is to find the optimal control strategy and to determine the optimal duration of a vaccination campaign adequate to eradicate the infection in Morocco For this, we introduce into the model a saturated vaccination function, which takes into account the limited resources on the COVID-19 vaccine, and we formulate a minimization problem where the final time is considered to be free The existence of optimal control is investigated The characterization of the sought optimal control and optimal final time is derived based on Pontryagin's maximum principle Using Matlab, we solve the optimality system with an iterative method based on the iterative Forward-Backward Sweep Method (FBSM) The numerical simulation results show the efficiency of a vaccination strategy on reducing the number of infectious individuals within an optimal period time which is approximately equal to 44 days


INTRODUCTION
The infectious COVID-19 spreads via direct contact or droplets and results for the majority of cases in initial symptoms including fever, cough, dyspnea, myalgia, or fatigue. Emergency symptoms can occur, they include intestinal symptoms like diarrhea, cardiac injury, chest pain or pressure, confusion, difficulty walking, and kidney injury. The complication of the infection would lead to death (see Chan et al. [1] and Sohrabi et al. [2]). Epidemic models are of very relevant interest in the background literature related to understanding the propagation, extinction, oscillatory behaviours and convergence properties of the state-trajectory solutions to their equilibrium points see Wang et al. [3], Soares and Bassanezi [4], and Köhler-Rieper et al. [5]. Many mathematical models have been developed in recent decades to gain insight into the transmissible nature of infectious diseases, such as SIR models, SIS models and SEIR models with or without time delays, see for example Beretta [6], Li et al. [7], Yew Ng and Mei Gui [8] and De Luca and Romeo [9]. When the pandemic start, we are interested to understand where and how did the pandemic start, what is the risk of its spread in the region and importation in other regions of the world, every one try to understand the pathogen and its epidemiological characteristics. At the stage where the pandemic takes hold, researchers investigating various intervention and control strategies; usually pharmacological interventions do not work in the event of a pandemic and thus non-pharmacological interventions are most appropriate. Newbold et al. [10] develop an epidemiological-economic model to examine the optimal duration and intensity of physical distancing measures aimed to control the spread of COVID-19. Their model was applied to the United States where they consider the trade-off between the lives saved by physical distancing both directly from stemming the spread of the virus and indirectly from reductions in air pollution during the period of physical distancing and the short and long run economic costs that ensue from such measures.
Optimal control of epidemics with purely non-pharmaceutical procedures was treated in Kantner and Koprucki [11]. They consider an extension of SEIR model and continuous-time optimal control theory to compute the optimal non-pharmaceutical intervention strategy, to minimize disease-related deaths and to establish a sufficient degree of natural immunity in order to exclude a second wave, when the eradication of the epidemic isn't possible. Their obtained solution was close to the stability boundary of the system, and they calibrate their model to reproduce the initial exponential growth phase of the COVID-19 pandemic in Germany. Chen et al. have developed a numerical solution for uncertain SIR model with application to COVID-19 (see Chen et al. [12]). They presented an α-path-based approach, that can handle the highdimensional SIR model, to calculate the uncertainty distributions and related expected values of solutions. Furthermore, they estimated the parameters by using the method of moments and designed a numerical algorithm to solve them. Their model was applied to the infected and recovered data of the province of Hubei to describe the development trend of COVID-19.
They affirmed that the lockdown policy achieves almost 100% efficiency after February 13, 2020. Chowdhury et al. [13] employed a multivariate prediction model, based on up-to-date transmission and clinical parameters, to simulate outbreak trajectories in 16 countries, from diverse regions and economic categories. In each country, they modelled the impacts on intensive care unit (ICU) admissions and deaths over an 18-month period. This multi-country analysis demonstrates that intermittent reductions of R 0 below 1 through a potential combination of suppression interventions and relaxation can be an effective strategy for COVID-19 pandemic control. Global stability for delay SIR and SEIR epidemic models with non-linear incidence rate was developed by Huanga et al. to show that the global properties of equilibria only depend on the basic reproductive number and that the latent period in a vector does not affect the stability, but the latent period in an infected host plays a positive role to control disease development (see Huanga et al. [14]). They incorporate time delays into the ordinary differential equation models and consider two delay differential equation models in which delays are caused by the latency of the infection in a vector, and by the latent period in an infected host. SEIR epidemic model for the spread of COVID-19 using the Caputo fractional derivative was investigated by Shahram [15]. By using the fixed point theory, the existence of a unique solution for the model was treated. Also, by using the fractional Euler method, they get an approximate solution to the model. The generalized fractional order SEIR model was employed by Xu et al. [16] to forecast analysis of the epidemics trend of COVID-19 in the USA, two modified models SEIQRP and SEIQRPD successfully capture the development process of COVID-19, which provides an important reference for understanding the trend of the outbreak. To control the diffusion of infectious disease such as COVID-19, Park and Kim analysed the basic reproduction number in South Korea that allow identification of a necessary level of vaccine stockpile in order to achieve herd immunity(see Park and Kim [17]. They adopted an susceptible-infected-susceptible model that allows a stochastic diffusion. Their results shows that the basic reproduction number of South Korea is substantially lower than those of the other regions and is approximately equal to 2. The herd immunity suggests that at least 62% of the susceptible population be vaccinated when the vaccine becomes available. Other works on the optimal control problem of COVID-19 can be found in [18,19,20,21] and the references therein.
With the lack of effective treatment for COVID-19, vaccination remains one of the possible solutions to control this infectious disease. In order to evaluate the impact of a vaccination campaign on the spread of COVID-19, we propose in this work an optimal control problem based on a SEUIRD model adapted to this epidemic and which takes into consideration real data from the epidemiological situation in Morocco. Our contribution essentially consists of estimating the parameters of our model from the number of cases communicated by the Ministry of Health in Morocco; the aim being to build a model that reflects, as possible as, the reality. Solving an optimal control problem where we considered a saturated vaccination function uS 1 + bS where u represents the percentage of susceptible people vaccinated. This function is the most appropriate for the current context, which is characterised by high international demand on available vaccines and limited manufacturing capacities. Also, through this work we propose an answer to a fundamental question concerning the optimal duration of the vaccination campaign in order to reduce the number of the infection cases.
The paper is organised as follows: in Section 2, we present our mathematical model and the parameters estimates. In Section 3, we formulate the optimal control problem, we prove the existence of a solution and put forward the control expression. In Section 4, we propose numerical simulations to show the efficiency of our control strategy. Finally, some conclusions are given in Section 5.

COVID-19 MATHEMATICAL MODEL
2.1. Model construction. In this paper, we consider a model of six non-linear differential equations which describe the dynamic of the COVID-19 disease. The total population, denoted N, is divided into six different compartments: • Susceptible individuals (S), people who may be infected by the virus; • Exposed individuals and asymptomatic individuals (E), infected with the virus but without typical symptoms of infection; • Unreported infection cases (U), who are infectious but not yet confirmed by a test; • Reported infection cases (I), people who are diagnosed as COVID-19 positive patients and are hospitalized; • Recovered individuals (R); • Individuals that deceased due to the disease (D).
We assume that the transmission of COVID-19 occurs following adequate contact between susceptible and infected people in E and U compartments. As the positive diagnosed people in I are isolated, we assume that they have no contact with susceptible and do not contribute to the disease spread. Due to the non-linear contact dynamics in the population, we use the incidence function β E +U N S to indicate successful transmission of COVID-19, where β denotes the infection contact rate. We assume that all newly infected individuals enter to the exposed compartment E for k −1 days (k is the rate at which individuals leaves the latent class by becoming infectious). A proportion p of the infectious individuals are diagnosed and enter the compartment I. While the remaining infectious patients are considered as free infectious people and they are regrouped in the unreported compartment U. Among infectious patients who are not yet detected and isolated, some of them are diagnosed at a rate γ. In our model, we use a time-dependent decreasing mortality rate m(t) given by the following formula From the clinical data, we have observed that there was an increasing change within time for the number of individuals who progress from I to the class R; to reflect this, we assume that people in I enter into the removed class with a time-dependent variable r(t), which is derived from the following logistic function: where t 1/2 represents the time at which r(t) reaches its half value, ∆ determines the width of r(t). The parameters r 0 and r f model asymptotic values.
The transfer diagram of our model is given in Figure 1.
The dynamics of the model is governed by the following differential equations: with initial data,

Positivity and boundedness of solutions. As the model (3) describes the temporal evo-
lution of a human populations, we shall show that the solutions remain non-negative and bounded.
From the system (3), we have We note that all above rates are non-negative on bounding planes of the non-negative cone of

Parameters estimation.
In this section, we estimate parameters of the model (3), as well as the unknown initial conditions E 0 and U 0 . For this, we will use the least squares data fitting method and the cumulative daily data of confirmed, removed and death cases form June 10,2020 to November 14,2020, reported by the Ministry of Public Health in Morocco and available in [22].
Firstly, we start by computing parameters m 0 and m 1 in (1). We assume that the cumulative numbers of reported cases and deceased people are respectively given by Combining equations (1), (6) andḊ(t) = m(t)I(t), we set Using real data of reported cases I(t) and deaths D(t) from [22] and the Matlab curve toolbox, we obtain the estimates of a 0 , a 1 , b 0 and b 1 (see Table 1). Thus, we get m 0 = 3.7429 × 10 −04 and m 1 = −6.9 × 10 −04 . Then, as the progression rate, k, form E to the infectious stage is equal to the inverse of the mean incubation period [23], and knowing that several studies show that the mean incubation period of COVID-19 is around 5 days (see for example [24,25,26]); we fix k = 1/5 per days.
Finally, the remaining parameters and the unknown initial conditions E 0 and U 0 are estimated by fitting in the model (3) to the observed data using the MATLAB routine lsqcurvefit. The resulting estimations are listed in Table 2 and the best model fit to the real data is shown in  Bootstrap sampling has been used to improve the reliability of the estimates of the parameters β , p, γ, r f , r 0 , t 1/2 , ∆ as well as the initial condition E 0 and U 0 . The distributions of the estimated parameters obtained from 1000 generated data sets are given in Figure 5, and the mean, the standard deviation, and 95% confidence interval are listed in Table 2. The results provide a good confidence level of the set of our estimated parameters.   we can limit our study to the following system with initial data, The control strategy we adopt consist of a vaccination campaign that aims to limit the spread of the COVID-19 disease among Moroccan population. We consider a control variable u(t) which represents the percentage of individuals being vaccinated per unit of time. As there is a limited resources and the whole population can not be vaccinated, we chose a saturated vaccination function given by u 1 + bS where b is the saturation factor that measures the effect of the susceptible being delayed for vaccination. A positive number, denoted ω, is used to model the efficiency of the vaccine. We assume that the vaccinated people will join the recovered individuals compartment. Our model with control is given as follow with initial data, Mathematically, we formulate an optimal control problem with free terminal time based on the objective functional where A and B are the weight constants of the control and time respectively.
In other words, we seek optimal control u * and optimal terminal time t * f such that (11) J(u * ,t * f ) = min J(u,t f ) | u ∈ U and t f ∈ R + , subject to system (9), with U ad is the set of the admissible controls, defined by

For the existence of a solution of the system (9), we note X(t) = [S(t), E(t), I(t),U(t), R(t)] T
and we write this system as follows The function F satisfies Furthermore, one has Thus, it follows that the function G is uniformly Lipschitz continuous. From the boundedness of the control u and the restriction on the state variables, we conclude that there exist a solution of the system (9) (see [29]). In the rest of this section, we will be interested in proving the existence of the optimal control and deriving its characterization. Theorem 1. Consider the control problem with system (9). There exist an optimal control u * such that (1) The set of controls and corresponding solutions to the system (9) is non-empty.

Existence
(2) The control set U ad is convex and closed. ( 3) The function f is bounded by a linear function in the state and control variables.
(5) The pay-off term at the terminal time in the objective functional φ X(t f ),t f = Bt 2 f is continuous. Proof.
• Condition 1 To prove that the set of controls and corresponding state variables is nonempty, we use a simplified version of an existence result [31] (see Theorem 7.1.1).
• Condition 2 The set U ad is closed and convex by definition.
• Condition 3 We have Using the boundness of the state variables, system (20) can be rewrite in matrix form: which gives a linear function of controls vector and state variables vector. Thus, we can Hence, we see that the right hand side of the state system is bounded by a sum of state and control vector. Therefore, condition 3 is satisfied.
• Condition 4 To verify condition 4, we note that the integrand in the objective functional (10) is quadratic in the control, it is convex on U ad . As the state variables are bounded by N, and by considering c 2 = N, c 1 = A 2 and ρ = 2 it is easy to check that the constants c 1 , c 2 > 0, and ρ > 1 satisfies L(X, u) ≥ c 1 |u| 2 β 2 − c 2 .
• Condition 5 It is obvious that φ X(t f ),t f is continuous.

Characterization of the optimal control.
In order to derive the necessary conditions for our optimal control problem, we use the Pontryagin's Maximum Principale given in [32].
This principal converts the problem (9)-(11) into a problem of minimizing a Hamiltonian, H, defined by where f i for i = 1, . . . , 5 is the right side of the differential equation of the i th state variable.
By applying the Pontryagin's maximum principle [32] and the existence result for optimal control and corresponding optimal states from [30], we obtain the following theorem: Theorem 2. Given an optimal control u * , an optimal terminal time t * f and solutions S * , E * , I * and R * of the corresponding state system (9), there exists an adjoint vector λ = [λ 1 , λ 2 , λ 3 , λ 4 , λ 5 ] that satisfies with the transversality condition λ i (t f ) = 0 for i = 1, . . . , 5.
Furthermore, the optimal control u * is given by and the optimal final time is given by .
Proof. The adjoint equations and transversality conditions can be obtained by using Pontryagin's Maximum Principle such that The optimal control u * can be solved from the optimality condition, with the control bounds in U ad , it is easy to obtain u * in the form (26).
The transversality condition for t f to be the optimal terminal time can be stated as where φ (t, X(t)) = Bt 2 that is

NUMERICAL SIMULATION
In this section, we illustrate the numerical simulations associated to the resolution of the optimality system. This system consists of the state system, adjoint system, the control characterization and separated boundary conditions at times t 0 = 0 and t f . To solve this two-point boundary value system, we use an iterative method based on a combination of forward and backward difference approximation, which converges when a tolerance criterion is reached.
The numerical scheme we propose to compute the solution of our optimality system and the optimal final time is summarized in the algorithm given below. Let X the state variable and λ the adjoint variable.

Algorithm
Step 0 • Make initial guess for the terminal time t f and the control u 0 .
Step 1 • Solve the state system (9) with initial condition X 0 forward in time using the stored values for X and u.
• Solve the adjoint system (25) with transversality condition λ (t f ) = 0 backward in time using the stored value for the controls and the state variable.
Step 2 • Update the control using new X, λ and the formula (26).
• Update the free terminal time using the formula (27).
Step 3 • Testing the convergence : if the difference of values of these variables in this iteration and the last iteration is sufficiently small, output the obtained current values as solutions.
If the difference is not considerably small, return to Step 1. In figure 6, we remark that the number of susceptible individuals decreases more rapidly during the vaccination campaign, it reach 160 at the day 157.
The vaccination contains the infection to the lower values on exposed individuals represented in the figure 7, it goes from 2928 at the day 14 to 0 at the day 157, while it grows to 5, 7 × 10 4 without control at the day 157.
In figure 8, the evolution of the infected individuals takes more time to decrease than the exposed one, it goes from 1, 975 × 10 4 at the day 48 to 4345 at the day 157, while without vaccination, the infected compartment of the epidemic is uncontrolled and goes to high levels represented on the value 2, 84 × 10 5 at the day 157. Figure 9 shows that the vaccination strategy gets the number of U to 0 at the day 157, while, without control, it still reaches high levels as the previous states.
The number of death is increasing in both the parts with and without control, but the figure 10 shows that we can limit it to 926 at the day 157, while it's 4864 at the same day without control.
The figure 11 represents the optimal control where we see that the full effort on u is applied during almost the entire control period.
Finally, we can conclude that the vaccination control strategy has an efficient impact on reducing the number of susceptible, exposed and infectious people in an optimal period time, which approximates 44 days as shown in figure 12.

CONCLUSION
Recent studies reported the severe consequences of COVID-19 on public health and the economy. The number of infection and deaths around the world motivates the researchers to work for a vaccine that can help to eradicate this pandemic. The concern of this paper is to contribute to the world efforts that aim at understanding the dynamic of this infectious disease and limiting its spread. By setting an optimal control problem based on a model that reflects the epidemiological situation of Morocco, and includes a saturated vaccination function; we investigate the impact of a vaccination campaign on reducing the number of infectious individuals as well as the optimal duration of this campaign. Our results show that vaccination is an effective strategy that can help to control the spread of COVID-19 in Morocco within an optimal time that approximates 44 days. Note that the model and the optimal approach used in this study can be applied to other countries. Also, taking into account the main role of the spatial diffusion of the COVID-19 disease, one possible extension of this work is to consider a regional control approach based on a reaction-diffusion model.

CONFLICT OF INTERESTS
The author(s) declare that there is no conflict of interests.