The COVID-19 pandemic began with a discovery of first cases of pneumonia of unknown origin in Wuhan, China, in late December of 2019. The city of Wuhan has been closed for quarantine since the 22nd of January of 2020. On the 30th of January, World Health Organization (WHO) recognized the outbreak of a new coronavirus (SARS-CoV-2) as a public health emergency of international concern; on the 11th of March it announced that the outbreak had become a pandemic and on the 13th of March that Europe had become its center. A hallmark of the current pandemic, which is making it difficult to deal with, is a varying in a wide range incubation period (in some cases over 14 days) and a very large number of asymptomatic patients who are contagious but demonstrate no clinical manifestations or very mild manifestations.
As no vaccine against COVID-19 is currently (the spring of 2020) available, the quarantine (or a regime of massive self-isolation) remains to be the only accessible control policy. The Chinese experience of the February-March of 2020 shows that the quarantine can effectively stop the spread of the infection and annihilate the virus. However, at the same time, the large-scale quarantine is also extremely expensive policy inflicting huge economic losses. Moreover, while some groups of a population, such as children and the retired, can be quarantined at a comparatively low cost, the cost quickly grows as more and more of the economically active people have to be isolated.
In order to investigate the feasibility of reducing the cost of the policy aimed at the control of the epidemics and annihilating the virus, in this paper we formulate and analytically and numerically consider two problems of optimal control of COVID-19 spread. Our results indicate that the epidemic can be stopped and suggest the cost-effective scenarios.
Let us consider the spread of COVID-19 in a human population of size N(t). We postulate that the population is divided into the following five classes:
• S(t) -the number of susceptible individuals;
• E(t) -the number of exposed individuals (in the latent state) who exhibit no symptoms and are no contagious; • I(t) -the number of infected individuals who have very mild symptoms or no symptoms at all (it is believed that for COVID-19 those are about 80% of the infectious individuals); • J(t) -the number of infected individuals who are seriously ill.
• R(t) -the number of recovered individuals. Hence, the natural equality
holds. Please note that the proposed population structure postulate that there are two parallel pathways of the disease progression, namely, symptomatic, S → E → J → R, and asymptomatic, S → E → I → R.
In this paper, we are only interested in a single epidemic which we assume to be reasonably short. Therefore, we ignore the demographic processes (that is, the new births and natural deaths) assuming that they occur at a considerably slower time scale. The the spread of COVID-19 can be described by the following system of differential equations: System (2.2) postulate that the susceptible individuals S(t) are infected through contact with asymptomatically and symptomatically infected people I(t) and J(t) at incidence rates f 1 (S(t), I(t), N(t)) and f 2 (S(t)J(t), N(t)), respectively. After an instance of infection the infected individual moves into the exposed compartment E(t) where remain for 1/γ days. (Hence, γ is the rate with which the exposed individuals move into the infectious groups.) Moreover, σ 1 and σ 2 are the fractions of the exposed individuals that moves into the classes of the asymptomatically infected individuals I(t) and the symptomatically infected individuals J(t), respectively. In this paper, we assume that σ 1 + σ 2 = 1. It is currently assumed that for COVID-19 σ 1 = 0.8 and σ 2 = 0.2; we use these figures in our computations. Parameters ρ 1 and ρ 2 are the removal (recovery + death) rates for the I(t) and J(t) groups, respectively. We assume that ρ 1 I(t) and (1 − q)ρ 2 J(t) are the recovery rates of individuals in the I(t) and J(t) groups, respectively, and that qρ 2 J(t) is the disease-induced death rate. (Hence, value q ∈ [0, 1] is the death probability.) Precise forms of the incidence rates f 1 (S(t), I(t), N(t)) and f 2 (S(t)J(t), N(t)), and, in particular, their dependence on the population size N(t) are not known. However, as we show below, the dependence of incidence rate on population size is very import for quarantine. In this paper we consider two rates, namely, bilinear rates
and rates
When the population size is constant, or its variations are reasonably small, these two incidence rate types are identical or lead to a qualitatively similar outcome. However, as we show below, when the quarantine is imposed, the corresponding models are principally different.
In this paper, we assume that the individuals in J(t) class exhibit the characteristic symptoms of the disease and, hence, are isolated, either at homes or in a hospital. We assume that this implies that their ability to spread the virus is significantly limited and postulate that β 1 > β 2 and β 1 > β 2 (and even β 1 ≫ β 2 , β 1 ≫ β 2 ) hold.
System (2.2) should be complemented by the corresponding initial conditions
We assume that at t = 0 the values S 0 , E 0 , I 0 , J 0 , N 0 are positive and R 0 ≤ 0. Moreover, the equality:
where N 0 is the initial population size, holds. We assume that model (2.2) is to be a object of a control, and, hence, we consider this model at a given time interval [0, T ].
2.1. Model 1. Model (2.2) combined with incidence rate (2.3) and defined on time interval [0, T ] yields the following system of differential equations:
with initial conditions (2.5) It is easy to see that the equations of system (2.7) together with initial conditions (2.5) and equality (2.6) imply relationship (2.1) and that the value of N(t) naturally varies (decreases due to disease-induced mortality).
For the sack of simplicity, let us perform the normalization of system (2.7) and initial conditions (2.5), using the following formulas:
The substitution yields system of differential equations
with initial conditions
Here s 0 , e 0 , i 0 , j 0 are positive, r 0 ≥ 0 (further in this paper we will use r 0 = 0 assuming that we study an initial stage of the epidemic), and equality
following from (2.6), holds. Please note that normalization (2.8) converts equality (2.1) into equality
The important properties of solutions for system (2.9) are established by the following lemma. The proof of Lemma 1 is deferred to Appendix A. Lemma 1 implies that all solutions of system (2.9) with the initial conditions (2.10) retain their biological meanings for all t ∈ [0, T ].
Let us introduce a control function u(t) into system (2.9). We assume that this control reflects the intensity of the quarantine. We assume that the quarantine implies "effective isolation" of equal fraction u(t) in groups S(t), I(t) and R(t) (or reducing the number of contacts of the individuals in these groups by factor (1 − u(t))), whereas the group J(t) (the symptomatically ill individuals) is assumed to be already isolated and, therefore, its status is not affected by the quarantine. The control satisfies the restrictions:
These considerations lead to the following control system:
with the corresponding initial conditions (2.10).
2.2. Model 2. For the incidence rates (2.4), the change in the size of the compartments is described by the following system of differential equations:
The corresponding initial conditions are defined by (2.5) and (2.6). Here β 1 and β 2 are the per capita rates of virus transmission from asymptomatic infected individuals I(t) and symptomatic infected individuals J(t), respectively. Substituting (2.8) into system (2.15) with initial conditions (2.5), we obtain the following normalized system of differential equations:
with the initial conditions (2.10).
The positiveness, boundedness, and continuation of the solutions s(t), e(t), i(t), j(t), r(t), n(t) on the entire interval [0, T ] for system (2.16) are established by arguments similar to those presented in Lemma 1. Now, we introduce into system (2.16) the control function (the "quarantine intensity") u(t). We have to stress that, due to the difference of incidence rates, the effect of quarantine on model (2.16) is different (smaller) compared to model (2.9).
As above, we assume that the quarantine implies isolation of fraction u(t) in groups S(t), I(t) and R(t), whereas the group J(t) is assumed to be already isolated. Furthermore, group J(t) is comparatively small, and, hence, we can assume that
and that the fraction u(t) of the population N(t) is isolated. These considerations leads to the following control system:
with the corresponding initial conditions (2.10). As above, we assume that control u(t) satisfies restrictions (2.13). Please note that for u(t) = 0 (in the absence of quarantine) systems (2.14) and (2.18) become systems (2.9) and (2.16), respectively. The different manifestation of control u(t) in models (2.14) and (2.18) shows why the dependence of an incidence rate on population size N(t) is so important for studying impacts of quarantine.
To formulate the optimal control problems, let us introduce set Ω(T ) of all admissible controls. The set is formed by all possible Lebesgue measurable functions u(t) that for almost all t ∈ [0, T ] satisfy restrictions (2.13). For control systems (2.14) and (2.18), on the set Ω(T ) of all admissible controls, we consider objective function
Here, α 1 , α 2 are non-negative and α 3 is positive weighting coefficients. For convenience of notation, in function (2.19) we below denote the terminal part by P .
In function (2.19) , the first two terms reflect the COVID-19 level at the end of quarantine period [0, T ] and the cumulative level over the entire quarantine period. The last term determines the total cost of the quarantine. We have to note that the actual cost function is unknown. However, it is obvious that quarantining some groups, such as children or the retired, is comparatively cheap, and that the cost rapidly grows as wider groups of economically active individuals have to be isolated. The quadratic cost function that we use in this paper qualitatively reflects this situation and, at the same time, is reasonably simple.
Please note that phase variables e(t), i(t) and j(t) are only present in the objective function (2.19) . Moreover, the first four equations of system (2.14 are independent of variables r(t) and n(t). Therefore, the last two differential equations can be omitted from system (2.14. Likewise, the fifth differential equation can be excluded from system (2.18).
As a result, we state the first optimal control problem (OCP-1) as a problem of minimizing the objective function (2.19) on the set Ω(T ) of all admissible controls for system
with initial conditions
where s 0 , e 0 , i 0 , j 0 are positive and satisfy equality (2.11). Likewise, we state the second optimal control problem (OCP-2) as a problem of minimizing the objective function (2.19) on the set Ω(T ) of all admissible controls for system
where s 0 , e 0 , i 0 , j 0 are positive and satisfy equality (2.11). It is easy to see that for OCP-1 and OCP-2 the hypotheses of of Theorem 4 (chapter 4, [8] ) are correct. This Theorem and Lemma 1 guarantee the existence of an appropriate optimal solution, which consists of • the optimal control u 1 * (t) and the corresponding optimal solutions s 1 * (t), e 1 * (t), i 1 * (t), j 1 * (t) to system (2.20) for OCP-1; • the optimal control u 2 * (t) and the corresponding optimal solutions s 2
3. Analysis of OCP-1
To study OCP-1 analytically, we apply the Pontryagin maximum principle ( [13] ). Firstly, we write down the Hamiltonian of this problem:
where ψ 1 , ψ 2 , ψ 3 , ψ 4 are the adjoint variables. The Hamiltonian satisfies:
By the Pontryagin maximum principle, for the optimal control u 1 * (t) and the corresponding optimal solutions s 1
With respect to u, this Hamiltonian is a quadratic function of the form
. Therefore, the following relationship holds:
Here function λ 1 * (t) is the so-called indicator function [14] , which for A * (t) = 0 is defined as
It determines the behavior of the optimal control u 1 * (t) according to (3.5). By (3.2), (3.4) and (3.6) we can see that A * (T ) > 0 and B * (T ) > 0, and, therefore, inequality λ 1 * (T ) > 0 holds. According to (3.5) , this means that the following lemma is true.
Lemma 2. At t = T , the optimal control u 1 * (t) is positive and takes either value λ 1 * (T ) or value u max . Now we can see that the following lemma is valid. 
Equality (3.7) shows that, for all values of t ∈ [0, T ], the maximum of Hamiltonian (3.3) is reached with a unique value u = u 1 * (t). Therefore, the following lemma immediately follows from Theorem 6.1 in [4] . Moreover, arguing as in [7, 10, 15] , due to the boundedness of the state and adjoint variables and the Lipschitz properties of systems (2.20) and (3.1) and relationship (3.7), one can establish the uniqueness of this control.
To analytically study OCP-2, we also use the Pontryagin maximum principle. The Hamiltonian of OCP-2 is
where φ 1 , φ 2 , φ 3 , φ 4 , φ 5 are the adjoint variables. The Hamiltonian satisfies
By the Pontryagin maximum principle, for the optimal control u 2 * (t) and the corresponding optimal solutions s 2
is the nontrivial solution of the adjoint system:
• the control u 2 * (t) maximizes the Hamiltonian (4.3)
H(s 2 * (t), e 2 * (t), i 2 * (t), j 2 * (t), n 2
with respect to u ∈ [0, u max ] for almost all t ∈ [0, T ]. Therefore, the following equalities hold:
. is the indicator function that determines the behavior of the optimal control u 2 * (t) according to (4.4) . Now, (4.2) and (4.5) yield
which, by Lemma 1, implies inequality λ 2 * (T ) > 0. According to (4.4) , this implies that the following lemma that is similar to Lemma 2 is true.
Lemma 7. At t = T , the optimal control u 2 * (t) is positive and takes either value λ 2 * (T ), or value u max .
with a unique value u = u 2 * (t). Therefore, the following lemma, which is similar to Lemma 5, immediately follows from Theorem 6.1 in [4] .
Lemma 8. The optimal control u 2 * (t) is a continuous function on the interval [0, T ]. Finally, the arguments in Remark 6 are also applicable here. Specifically, systems (2.22) and (4.1) with initial conditions (2.23), (4.2), relationship (4.4) and equality (4.5) form a two-point boundary value problem for the maximum principle. The optimal control u 2 * (t) satisfies this boundary value problem together with the corresponding optimal solutions s 2 * (t), e 2 * (t), i 2 * (t), j 2 * (t), n 2 * (t) for system (2.22). Moreover, due to the boundedness of the solutions for the state and adjoint variables and the Lipschitz properties of systems (2.22) and (4.1), and equation (4.4) that defines the control, it is possible to establish the uniqueness of this control.
The basic reproduction number ℜ 0 typically characterizes the ability of an infection to spread: it is usually assumed ( [16] ) that an epidemic occurs if ℜ 0 > 1. If ℜ 0 < 1, then the epidemic gradually fades ( [11] ). To find ℜ 0 for systems (2.9) and (2.16), we apply the Next-Generation Matrix Approach [17] . For both these systems, the basic reproduction number is the same:
Further in this paper we use the following values of parameters: We also assume that
For these values of parameters, (5.4) ℜ 0 = 11.62 · β 1 . Table 1 shows the relationships between ℜ 0 and coefficients β 1 and β 2 that are assumed to be related via (5.3) and (5.4) . Further in this paper we use a value of ℜ 0 from {2.5; 3.0; 4.0; 6.0} (see [2, 3, 6, 9, 12, 18, 19] ). Table 1 . Values of parameters β 1 and β 2 depending on ℜ 0 .
For control systems (2.14) and (2.18) , the corresponding reproduction numbers ℜ 1 (u) and ℜ 2 (u) calculated under assumption of the constant control are
Parameters from (5.2) yield ℜ 1 0 (u max ) = 0.01β 1 (0.8 · 14 + 0.2 · 21) = 0.154β 1 , ℜ 2 0 (u max ) = 0.1β 1 (0.8 · 14 + 0.2 · 21) = 1.54β 1 . It is easy to see that for all values of β 1 from Table 1 inequalities
This means that quarantine with the maximum intensity u max = 0.9 should eventually stop the epidemic. It is also clear that for both these systems u max can be reduced.
We conduct numerical study of OCP-1 and OCP-2 using BOCOP 2.0.5 package [1] . BOCOP 2.0.5 is an optimal control interface implemented in MATLAB and used for solving optimal control problems with general path and boundary constraints and with free or fixed final time. By a time discretization, such optimal control problems are approximated by finite-dimensional optimization problems, which are then solved by the well-known software package IPOPT using sparse exact derivatives that are computed by ADOL-C. IPOPT is an open-source software package designed for large-scale nonlinear optimization. In the computations, we set the number of time steps to 5000 and the tolerance to 10 −14 and use the sixth-order Lobatto III C discretization rule [1] .
In the computations we use parameters from (5.2) and Table 1 The following figures present some results of computations for OCP-1 and OCP-2. Though we investigated these problems for ℜ 0 ∈ {2.5; 3.0; 4.0; 6.0} and for time intervals These figures demonstrate that even for this comparatively low basic reproduction number the 15-days quarantine is insufficient to eliminate the epidemic, and in particular for OCP-2. One can see that the graphs of i 1 * (t) and i 2 * (t) do not even start to decrease. (Please, note the sums of the appropriate graphs i 1 * (t) + j 1 * (t) and i 2 * (t) + j 2 * (t) represent the dynamics of the infectious at moment t.) We have to noted that the results for ℜ 0 = 4.0 are very similar to those for ℜ 0 = 3.0. If the virus is more contagious, that is for ℜ 0 = 6.0, then the situation is worse: as one can The results for for the 30-day policy for ℜ 0 = 3.0 are presented in Figures 6.9 and 6.10. One can see in these figures, that for OCP-1 the total infectious level (the graph of i 1 * (t)) passes its maximum and starts to decrease. However, it is noteworthy that for OCP-2 the level of infected i 2 * (t) in increasing even under the strongest quarantine measures. One conclusion that has to be withdrawn from these results is that the importance of actual dependency of the incidence rate on the population size N(t). 
Figures 6.11 and 6.12 that are given for the high infectivity level ℜ 0 = 6.9 show that the longer 60-days quarantine give considerably better results for both OCP-1 and OCP-2. For both these problems, even in the case of the high communicability of the infection, corresponding graphs for the infection level i 1 * (t) and i 2 * (t) pass their maximums and start decreasing. One can expect that for a lower ℜ 0 the outcomes would be considerably better, and our computation confirm this assumption. 
For the both problems, the optimal control u * should be kept as high as possible from begin of the policy and till a considerable decrease of the infection level is reached. Then it can be slowly relived.
As the impact of quarantine in OPC-2 is smaller, the corresponding optimal control u 2 * (t) is to keep the value of u max during the entire isolation period. The optimal control u 1 * (t) should be maximal during the first month and then it is slowly decreasing. u (quarantine intensity) Figure 6 .12. OCP-2: optimal solutions and optimal control for ℜ 0 = 6.0 and T = 60 days: upper row: i 2 * (t), j 2 * (t), i 2 * (t) + j 2 * (t); lower row: e 2 * (t) + i 2 * (t) + j 2 * (t), r 2 * (t), u 2 * (t); s 2 * (T ) = 0.999550, n 2 * (T ) = 0.999982.
In this paper, two SEIR type models that describe the spread of COVID-19 virus in a human population of variable size are considered. The models differ by the incidence rates that describe the virus transmission. A bounded control function that reflects the intensity of quarantine measures in the population was introduced into each of these models. This control reflects all sorts of the direct and indirect measures (quarantine, mask-wearing, various educational and information campaigns) aimed at reducing the possibility of transmission of the virus from infected to healthy individuals. The very first observation that can be done is that for these models the impacts of the quarantine is very different. This observation dignifies the importance of the actual form of the dependency of the incidence rate on the population size.
For each of the control models, the optimal control problem, consisting in minimizing the Bolza type objective function, was stated. Its terminal part determined the level of disease in the population caused by COVID-19 at the end of the quarantine period, and its integral part was a weighted sum of the cumulative level of disease over the entire quarantine period with the total cost of this quarantine. A detail analysis of the optimal solutions for these optimal control problems (OCP-1 and OCP-2) were made using the Pontryagin maximum principle. The properties of the corresponding optimal controls were established. Then, the values for the control models parameters, based on the knowledge of their basic reproductive ratios, were estimated. For these parameters, the results of computations performed using BOCOP 2.0.5 software were presented and discussed.
In our study we assumed that there is neither vaccine, nor drug available for the disease treatment. By the term "Quarantin" we mean all direct (isolation) and indirect protective measures during a specified period. Based on our analysis and computations, we can make the following conclusions:
It looks like for both, OCP-1 and OCP-2, the quarantine time T = 60 works better. However, the type of the optimal control seems to be more realistic for OCP-2 (especially when the quarantine is short), because it is not clear how the policy presented by Figures 6.1, 6 .3 or 6.5 would be implemented. The optimal control, on the other hand, presented by the graphs on Figures 6.2, 6.4 or 6.6 is understood as follows:
• The COVID-19 epidemic can be stopped by quarantine measures, and the virus can be eliminated by decreasing the level of infection below its survival level (e.g., less than 1 infected individual in each location). This outcome can be reached even without vaccination, by quarantine measures only. • However, one cannot expect that a short-term quarantine (e.g. 2-week, or even 1-month quarantine) would bring a decisive outcome. Even a 1-month quarantine can be insufficient to decrease the infection below its survival level. • It is necessary to keep the strongest quarantine for the most time of the planned period. The quarantine intensity can be gradually made easier only when the infection level reaches a certain low level.