Optimal control of COVID-19

Article History: Received 21 April 2020 Accepted 10 September 2020 Available 30 January 2021 Coronavirus disease of 2019 or COVID-19 (acronym for coronavirus disease 2019) is an emerging infectious disease caused by a strain of coronavirus called SARS-CoV-22, contagious with human-to-human transmission via respiratory droplets or by touching contaminated surfaces then touching them face. Faced with what the world lives, to define this problem, we have modeled it as an optimal control problem based on the models of William Ogilvy Kermack et Anderson Gray McKendrick, called SEIR model, modified by adding compartments suitable for our study. Our objective in this work is to maximize the number of recovered people while minimizing the number of infected. We solved the problem theoretically using the Pontryagin maximum principle, numerically we used and compared results of two methods namely the indirect method (shooting method) and the Euler discretization method, implemented in MATLAB.


Introduction
COVID-19 disease, which appeared in China in late 2019, is caused by SARS-CoV-2, a virus that belongs to the large family of coronaviruses. Very common, they can cause a simple cold as well as a serious respiratory infection like pneumonia, causing fatal epidemics as it was the case with Sras or Mers and this year, COVID-19. This virus passes into humans via animal secretions, under special conditions, then it is transmitted from human to human. Long-term contact is required for the transmission of this virus. As with many infectious diseases, people with pre-existing chronic diseases such as hypertension, cardiovascular disease, diabetes, liver disease, respiratory disease, etc. Seem more likely to develop severe forms. Note that there are serious forms of the disease, even in young adults. Children under the age of fifteen are unlikely to trigger a severe form of coronavirus [CSSE]. In our work, we used a compartmental model in epidemiology based on the work of Kermack and Mckendrick [1], using ordinary differential equations we modeled this epidemic as an optimal control problem, our objective (and of course the objective of all humanity) is to maximize the number of recovered and minimize the number of infected and then draw the attention of the Algerian people to the seriousness of the situation. Initially the two authors proposed a founding model in which the population was divided between individuals susceptible to contracting the disease (compartment S), and infectious individuals (compartment I), it is noted model SI. Then they added the proportion of healed people, we get the SIR model [2].
In the first part of this work, we consider the SIR model with adding other compartments namely the proportion of individuals most exposed to the disease (chronic diseases and the elderly), people who have recovered as well as those who have died. In this model, compartment S is necessary, since it must initially exist individuals who have not yet been infected, when an individual *Corresponding Author in compartment S is exposed to the disease, he does not necessarily become capable of transmitting it immediately, the virus requires a few days to make the individual sick [3], this is called the latency period, where the existence of compartment E then compartment I. After an individual has been infected, two scenarios can occur. First, the individual can die, in which case he belongs to compartment D, or recovered, in this case he belongs to compartment R. In our work, the command u(.) represents the vaccination rate [4]. Our goal is to minimize the number of infected people and maximize recoveries. The study is made in the case of Algeria. We estimate the states at 70 days and 150 days, from 02 April 2020 when the number of infected reached 986 people. This work is structured as follows: after a brief introduction, in section 2, we defined the model used and explained the importance of each of its compartments. Section 3 is devoted to the theoretical resolution of the modeled problem, using the Pontryagin maximum principle [4][5][6][7][8]. The numerical resolution is evoked in section 4, we solved the considered problem with two numerical methods, namely the direct method (Euler discretization ) and an indirect method (shooting method). We developed an implementations with the MATLAB programming language, then we presented some simulation results. In the numerical application, two scenarios were considered, T = 70 days and T = 150 days, with and without an effective vaccine. In the rest of the work, we considered compartment Q of people quarantined, the obtained problem is solved thoerically then numerically in section 5. In the rest of the work, in the absence of treatment, the only remedy to face this epidemic is quarantine. We discussed and commented the results found with two methods of the previous section, in section 6. Then we finished our document with a conclusion.

Statement of the problem
Starting from the basic model (SIR), we tried to improve the equations so that it better correspond to reality. Indeed, many parameters come into play during the spread of a virus and are not taken into account in the basic model. The letter used to represent a compartment is also used to represent the number of individuals in the compartment. At all epidemic, some susceptible peoples become infected. The study is done over a period of T = 70 days, and T = 150 days. For t ∈ [0, T ], denote α the infection rate, S(.) verify the following relation: As in any epidemic, infected people are vaccinated, note u(t), t ∈ [0, T ] the vaccination rate, 0 ≤ u(.) ≤ u max , in this case we will have: I(.) Is such that: where E(.), t ∈ [0, T ] denotes people witch are exposed to the epidemic, it verify the following differential equation: As in any epidemic, infected people are either recovered (compartment R) after being vaccinated or dead (compartment D), this is governed by the following differential equations: The positive coefficients are such that: α > 0 is infection rate, incubation rate β, is the rate of latent individuals becoming infectious. Given the known average duration of incubation Y , β = 1 Y , γ = 1 K , is determined by the average duration of recovery K, of infection. δ is the mortality rate. The evolution of the considered states is shown in the following figure(see Figure 1):  After this modeling, our goal is to minimize the number of infected people and maximize recoveries . The objective function in our case is: We will therefore solve a following optimal control problem:

Theoretical resolution
Taking the case of Algeria, out of 10.0000 healthy people, 3000 are exposed to the virus (among others doctors, nurses and people at risk). As of 02 April 2020, 986 person are infected in this country, 83 died and 108 were cured(Pasteur Institute results).
Each exposed person (E) is transformed into the infected person (I) with rate β, the incubation period refer to the time of infection to clinical symptoms of the disease. We shall assume that incubation period is a maximum of 7 days then β = 1 7 . An infected person is healed up after 6 days, as well γ = 1 6 [3]. This obtained problem is solved theoretically by using the Pontryagin maximum principle(PMP), then numerically, we implemented two methods on MATLAB namely the indirect method (shooting method) and a direct one.
The Hamiltonian of problem (7) is given for t ∈ [0, T ] as follows: where, x(t) = (S(t), I(t), E(t), R(t), D(t)) and The transversality conditions are such that: That gives: Let us Calculate the optimal controle u(.): The Hamiltonien's maximum is given by: The control which maximizes the Hamiltonian is Then:

Numerical implementation and methods
To compute the solution of the optimality system and maximize the objective functional J(u), two implementations are considered: indirect and direct methods, these implementations were compared in terms of model dynamics for the first 70 days (t ∈ [0, 70]) using both the parameter values α, β, γ, and σ in section 3, and the initial conditions Indirect method to solve optimal control problems are based on the Pontryagin's maximum principle to compute the optimal solution. Therefor, these methods require the state system with initial value [4], and necessary conditions [4,5,7], to find the solution of a given optimal control problem, by convering it into a boundary value problem.
On the other hand, direct methods use discretization with respect to time to solve optimal control problems, in which the cost functionl is directly optimized, treating the optimal control problem as a nonlinear optimization problem (N LP ). Some advantages and drawbacks can be pointed out to both indirect and direct methods. According to Trélat [4], direct methods are more robust and less sensitive to the choice of the initial conditions than indirect methods, being more easier to implement. However, in comparison with indirect methods, reaching to a desirable precision is not so easy when direct methods are employed. In addition, the author states not only the possibility of obtain local minima when the direct discretization of an optimal control problem is employed, but also the necessity of a large amount of memory, which in turn can lead to inefficiencies when, for instance, a large dimension problem is considered. On the other hand, indirect methods provide high levels of numerical accuracy, but their implementation can be quite difficult due to the necessity of computing derivatives and necessary conditions related to the PMP.

Resolution by Shooting method
The shooting method is based on the Pontryagin's maximum principle [9]. It consists in finding a zero of the shooting function associated with the original problem. It is a fast, high-precision method that does not require assumptions about the control structure. The shooting method consists in three main steps [7,8]: • Step 1: Form a boundary value problem using the model equations and the adjoint vectors equations as well as the transversality conditions. • Step 2: Determine the shooting function.
• Step 3: Solve a system of nonlinear equations.
Using shooting method, in the presence of an effective vaccine, after 70 days, the results found are illustrated in the following figure (see Figure2): After 70 days, in the presence of an effective treatment, the number of people recovered is very high, the number of infected decreases considerably.
In the absence of an effective vaccine, we will use the shooting method, after for T = 150 days, the results are illustrate in the following figure(see After 150 days, in the absence of an effective vaccine, the number of susceptible individuals does not almost decrease, the number of deaths increases. The number of susceptible individuals as well as the number of exposed to the virus increases considerably(see Figure 3). The execution time of the shooting method is CP U time equals 3.18s.

Resolution by the Euler discretization method
For a number of subintervals N chosen in advance, we will have the step of discretization h = T N and the following times: The application of the Euler discretization scheme for solving boundary value problems gives us the following nonlinear programming problem: (14) We also solved the nonlinear program (14)   After 70 days, in the presence of an effective vaccine, the number of infected individuals decreases, and the number of deaths stabilizes around 180 dead from 40 days. The number of people recovered is considerable(see Figure 4). It has been noted in the previous sections that in the absence of treatment the number of infected persons increases considerably, the best solution is therefore isolation. The execution time of the shooting method is CP U time equals 26.85s.

Advantage of quarantine
Based on China's experience, the vast majority of people who contracted COVID-19 contracted it from contact with infected relatives, hence the importance of isolation so that the epidemic can be controlled. taking the example of the Ebola epidemic and the case of China with this new virus, peoples had been deprived of their freedom of movement, this has paid off, the number of infected individuals has decreased considerably. Based on this protection strategy, we modeled another problem inspired by the SEIRD model, by adding the compartment of quarantined people (see Figure 5), the following optimal control problem has been obtained: (15) Where: α is the infection rate, β −1 average latente time, κ −1 average isolation time ,γ is cure time and σ is the mortality rate. We noticė S(t) +Ė(t) +İ(t) +Q(t) +Ṙ(t) +Ḋ(t) = 0, therefore S(t)+E(t)+I(t)+Q(t)+R(t)+D(t) = N P , N p ∈ R, and like We can simplify the problem.

Theoretical resolution
The Hamiltonian of problem (16) is given for t ∈ [0, T ] as follows: where, x(t) = (S(t), E(t), I(t), Q(t), R(t)) and is adjoint vector, it's solution of the system: The transversality conditions are such that: That gives: Let us Calculate the optimal controle u(.): The Hamiltonien's maximum is given by where H * = H(x * (t), p * (t), u * (t)). The control which maximizes the Hamiltonian is Then: u * (t) = −u max min(0, sign(p S (t) − 1)).

Numerical implementation
To compute the solution of the optimality system and maximize the objective functional J(u), two implementations are considered: indirect and direct methods, these implementations were compared in terms of model dynamics for the first 70 days (t ∈ [0, 70]) using both the parameter values α, β, γ, and σ in section 3, and the initial conditions:

Numerical comparison with results discussion
The direct method is characterized by the simplicity of its implementation, without a priori knowledge of the solution, this method is not very sensitive to the choice of the initial condition, but it has low numerical precision. The indirect method has very high digital precision, but it is very sensitive to the choice of the initial condition. We tried by this work and with two methods, to evaluate the considered states (in particular the infected and the recovered) with and without effective vaccine against Coronavirus. The results found by those methods are similar, but at least, the results found by shooting method are more accurate and the execution time of the shooting method is much less.

Conclusion
In this work, we have modeled a topical problem (COVID-19 problem) as a nonlinear optimal control problem. We considered the SEIR model by adding a compartment D for dead persons. In the first part of the work, we added compartment D, then we discussed and compared the results, in the second part of the work, we included in the considered model, compartment Q of people under quarantine. First, the Pontryagin's maximum principle has been used which gives a necessary condition of optimality. By implementing two numerical methods (the shooting method and the Euler discretization method) on the MATLAB software, we have traced the evolutions of the states considered in the problem, at a duration of 70 days and 150 days with and without vaccine. The results found show that the shooting method is fast and gives accurate results. We deduce that finding an effective vaccine is more than essential to fight against this virus. if not, the solution is quarantining the infected (see Figure 8).