Optimal control for a SIR epidemic model with limited quarantine.

Social distance, quarantines and total lock-downs are non-pharmaceutical interventions that policymakers have used to mitigate the spread of the COVID-19 virus. However, these measures could be harmful to societies in terms of social and economic costs, and they can be maintained only for a short period of time. Here we investigate the optimal strategies that minimize the impact of an epidemic, by studying the conditions for an optimal control of a Susceptible-Infected-Recovered model with a limitation on the total duration of the quarantine. The control is done by means of the reproduction number [Formula: see text], i.e., the number of secondary infections produced by a primary infection, which can be arbitrarily varied in time over a quarantine period T to account for external interventions. We also assume that the most strict quarantine (lower bound of [Formula: see text]) cannot last for a period longer than a value [Formula: see text]. The aim is to minimize the cumulative number of ever-infected individuals (recovered) and the socioeconomic cost of interventions in the long term, by finding the optimal way to vary [Formula: see text]. We show that the optimal solution is a single bang-bang, i.e., the strict quarantine is turned on only once, and is turned off after the maximum allowed time [Formula: see text]. Besides, we calculate the optimal time to begin and end the strict quarantine, which depends on T, [Formula: see text] and the initial conditions. We provide rigorous proofs of these results and check that are in perfect agreement with numerical computations.

Introduction

The Covid-19 pandemic outbreak raises an unprecedented series of decisions in different countries around the world. Since vaccines and effective pharmaceutical treatments were not initially available, governments had decided to impose non-pharmaceutical interventions like social distance, quarantines and total lock-downs as the most effective tools to mitigate the spread of the disease. Although these kinds of measures are helpful in reducing the virus transmission and giving time to health systems to adapt, they could be extremely stressful in terms of economic and social costs, and, in longer periods, tend to have less compliance with the population.

In this article we consider the classical SIR model introduced by Kermack and McKendrick[1] and widely used in epidemiology[2,3], where the population is divided in compartments of Susceptible, Infected and Recovered (or Removed) individuals. As it is usual in SIR models, we assume that people who have recovered develop immunity and, therefore, would not be able to get infected nor infect others. We consider that infection and recovery rates are allowed to change over time, and that are homogeneous among the population, instead of heterogeneous rates, not depending on age, individual protection measures, or awareness[4,5]. Let us observe that these heterogeneous rates implies the existence of effective rates, obtained as weighted means of the individual rates[6]. However, for age-structured models, better results are obtained by considering integro-differential equations, and different tools from mathematical control theory are needed[7]. We also assume a mean-field hypothesis that implies random interactions between any pair of agents, unlike other works[8-10] where interactions are mediated by an underlying network of contacts. In this case, by weighted means of agents degrees, it is possible to derive a system equivalent to a SIR model[11]. Also, a full proof of a derivation using probabilistic tools was studied recently[12] together with the optimal control problem for vaccination.

Optimal control problems for a system governed by a SIR or a SEIR model (with the addition of the Exposed compartment) with pharmaceutical interventions as vaccination or treatment were widely studied[13-16]. On its part, in the field on optimal control problems, non-pharmaceutical interventions were studied mostly for a control consisting of isolation acting only on the infected[17-20]. Taking into account that there is a window of time when the infected are not detected, in this article we will consider that the quarantine is applied to the whole population. Non-pharmaceutical interventions can range from a mild mitigation policy to a strong suppression policy. As discussed by Ferguson et al.[21], a suppression policy “aims to reverse epidemic growth, reducing case numbers to low levels and maintaining that situation indefinitely”. Suppression can be achieved by restricting travels, closing schools and nonessential businesses, banning social gatherings, and asking citizens to shelter in place. These measures, often referred to as a lockdown, are highly restrictive on social rights and damaging to the economy. In contrast, a mitigation policy “focuses on slowing but not necessarily stopping epidemic spread”. Mitigation measures may involve discouraging air travel while encouraging remote working, requiring companies to provide physical separation between workers, banning large gatherings, isolating the vulnerable, and identifying and quarantining contagious individuals and their recent contacts.

A critical parameter in the SIR model is the basic reproduction number , defined as the expected number of cases directly generated by one case in a population where all individuals are susceptible to infection. At the beginning of the epidemic, when no one in the population is immune, infected individuals will infect other people on average. Let us observe that, for , the number of new cases decline, and when , the number of new cases grows. However, at any time , the effective reproduction number replaces , since the number of contacts between infected and susceptible agents is reduced due to the interactions with recovered individuals that are immune. Hence, the epidemic grows until a sufficient fraction of the population becomes infected, and after reaching a peak starts to gradually decline. Following Ferguson et al.[21], for the Covid-19 the suppression phase can achieve , while the mitigation measures are unlikely to bring below 1. Therefore, the number of new cases are expected to decline during the suppression phase and to start rising again during the mitigation phase, although at a slower rate than in a non-intervention scenario.

In this work we assume that the intervention will occur in a preset period of time T, as proposed by Greenhalgh[22] and recently by Ketcheson[23], since it is unrealistic that interventions can be sustained indefinitely. Also, the lockdown or (strict quarantine) can last at most for a period , the maximum time that the population will adhere. Now, there are several interesting questions related to the implementation of the measures:

When should the suppression policy begin in [0, T]?

Is it convenient to split the maximum time into different intervals?

Is it better to apply a strong lockdown followed by mild mitigation measures or not?

In this article we study the previous questions using optimal control tools and numerical computations[24-27]. The answers clearly depend on the goal, which in our case is to minimize the overall impact of the epidemics in terms of the final number of infected individuals and the social and economic cost of the interventions, which we assume to increase as the quarantine becomes more strict. To account for quarantine measures, we consider a time-dependent reproduction number. Using an optimal control approach we show that the optimal strategy is of a single bang-bang type, that is, the lockdown or strict quarantine is applied in a single interval of time. Moreover, we characterize the time to start and finish the lockdown during the intervention phase. Let us remark that these questions make sense also in SIHR models, which include hospitalized individuals, since it must be necessary to keep the maximum of the hospitalized group below some threshold[28].

Recently, many works have appeared dealing with these and related issues. The optimal time to start the suppression measures that maximizes this type of objective function was studied by Ketcheson[23], where it was proved that a bang-bang control is optimal. However, in that work it is assumed that the lockdown corresponds to a zero reproduction number, something that is impossible to achieve in the real world. Moreover, it is assumed that the strict lockdown can last during the whole intervention, which seems to be impracticable. This problem was also analyzed for a different objective function[29], i.e., minimizing the peak of infected individuals, for which they proved that the optimal policy is not bang-bang. Besides, Kruse and Strack[30] minimize a functional that depends on the number of infectives during the intervention plus a term that measures the social and economic cost of interventions, and prove that the optimal control is bang-bang, but they do not investigate the optimal time to start the suppression policy.

The second question is suggested by the strategy proposed by Ferguson et al.[21]: the lockdown must be turned on and off several times based on the incidence of the virus in the population. A control-theoretic approach was considered in several works[28,31-33], although no time limits for the interventions were imposed. We shall see that the optimal policy is of bang-bang type, which consists on turning the lockdown on only once and turning it off after the maximum allowed time , in agreement with other authors[23,30].

Finally, the third question involves both suppression and mitigation phases, and one of the policies was colorfully characterized as the hammer and the dance in[34]: a strict lockdown, followed by mitigation measures in order to keep under control the propagation of the disease. However, our main results indicate that the best strategy actually depends on the initial condition, determined by the relation between and the initial fraction of susceptible individuals . On the one hand, when is smaller than the optimal strategy consists on applying a strong lockdown right at the beginning of the intervention period [0, T] for the maximum time period , followed by a mild mitigation measure until the end of the intervention (strong-mild strategy). On the other hand, when is larger than the optimal strategy is to apply mild mitigation measures at the beginning of the intervention, followed by a strong lockdown, and then a mild mitigation measure again in some cases (mild-strong-mild strategy). In this case, the optimal time to start the strong lockdown depends non-trivially on the initial condition.

The paper is organized as follows. We start by describing the basic SIR model in “The SIR model” section, and by introducing the SIR model with control in “The SIR control model” section. We then present the main results of the article supported by numerical simulations in “Results” section, followed by a discussion and conclusions including future work in “Discussion and conclusions” section. Finally, in “Methods” section we provide rigorous proofs of the results by applying Pontryagin’s maximum principle to the control problem, and we prove that the optimal control is bang-bang in Lemma 2. The main results that characterize the optimal control are given in Theorem 2, Corollary 1 and Theorem 3, for different case scenarios.

The SIR model

The basic SIR compartmental model of infectious diseases introduced by Kermack and McKendrik[1] considers a population of individuals that is divided in three compartments with homogeneous characteristics: Susceptible (S), Infected (I), and Recovered or Removed (R). The fraction of susceptible, infected and recovered individuals at time t is denoted by x(t), y(t) and , respectively. It is assumed that each infected individual is in contact with an average number of random individuals per unit time, and that infects only those who are susceptible ( transition), generating new infections at an average rate . Besides, each infected individual recovers at a rate ( transition). Births and deaths are neglected, and the recovered population is assumed to no longer infect others and cannot be reinfected. The infection and recovery rates and , respectively, are related to the basic reproduction number by , that is, the mean number of infections produced by a single infected individual in a totally susceptible population () during its mean infectious period . Then, the evolution of the system is governed by the following set of coupled nonlinear ordinary differential equations (see Hethcote[26], section 2.1): with . Here and are short notations for the time derivatives dx/dt and dy/dt, respectively. The region is forward-invariant and there exists a unique solution for all time[26]. Then, since , the proportion of infectious individuals is positive at any time. Even though the temporal dynamics of Eq. (1) depends on both and , the set of system’s trajectories on the space depends only on the basic reproduction number because only affects the overall time scale of the system. In Fig. 1 we depict typical trajectories starting from different initial conditions x(0), y(0), for [panel (a)] and [panel (b)].

Figure 1

Trajectories of the system in the phase space for the SIR model Eq. (1) with basic reproduction number (a) and (b), where x and y are the fractions of susceptible and infected individuals, respectively. Each curve corresponds to a trajectory starting from a given initial condition (x(0), y(0)), as indicated in the legends. Arrows denote the direction of the time evolution of x and y. The vertical dashed line corresponds to the critical value where .

The system of Eq. (1) is at equilibrium if . This equilibrium is stable only if , a condition referred to as herd immunity. If this condition is not satisfied at the initial time (), then y(t) first increases until it reaches its maximum value at a time t for which (dashed vertical lines in Fig. 1), and then decreases and approaches zero asymptotically, i.e., . That is, for , and is y(t) larger than zero for any finite time . The fraction of susceptible individuals x(t) is strictly decreasing, and its value in the long time limit is always positive. Therefore, the state of the system in the long time limit consists only of susceptible and recovered individuals, , where . Also, it is known that (see Theorem 2.1 of the work by Hethcote[26]).

The SIR control model

We now extend the classical SIR model to address the problem of controlling the spread of an epidemics with no access to vaccination, where the only possible control is isolation. We model this non-pharmaceutical intervention via a time dependent reproduction number that can be varied in the interval , where corresponds to a more strict isolation (“strict” quarantine) than (“mild” quarantine), with , and assume that this intervention can only be applied over a finite time interval [0, T]. Here T is the length of the intervention period. After the intervention, the restrictions are removed, thus the disease spreads freely and for all . We think of the control parameter as capturing political measures such as social distancing, and the lockdown of businesses, schools, universities and other institutions. Then, the system evolves according to the following set of coupled nonlinear ordinary differential equations: with , for and for , where . We also assume that during the intervention period [0, T] it is not possible to impose an extremely restrictive isolation for a long time. Thus, we consider that the strict quarantine—corresponding to —can last at most for a fixed time period , with . Once the period of intervention is finished at time T we compute , where (x(t), y(t)) is the solution of the system of Eq. (2) with initial condition (x(T), y(T)) and constant reproduction number for . Note that in this case, from Hethcote[26] we deduce that .

While political measures reduce the spread of the disease, they often come at an important economic and social cost. A long and strict quarantine can be very effective at reducing contagions, but at the expense of having a negative impact on the economy. Our goal is to find the optimal control on the SIR model described above that minimizes the total damage of a pandemic in terms of both, the total number of infections and also the socioeconomic costs. We model this trade-off by considering a global cost capturing the total number of individuals that were infected during the epidemics, i.e., those who are recovered in the long-time limit , and the socioeconomic cost of shutting down society during the intervention on [0, T], which we assume to increase as decreases (more restrictions). In order to find the optimal it proves convenient to work with the fraction of susceptible individuals in the long term instead. Then, given that minimizing is equivalent to maximizing , since , we define the functionalwhere

Our goal is to maximize the functional J, which has the following interpretation. The first term of J is the fraction of individuals that remain susceptible in the long term , and that we want to keep as large as possible subject to the condition of maximizing the second term of J as well, the functional . The functional is taken to be inversely proportional to the socioeconomic cost of the intervention ( increases as the cost decreases), as the function L is assumed to be a monotonously increasing function of . Then, an increase of the socioeconomic cost is achieved by decreasing (more restrictions or stricter quarantine), which leads to decreasing L and consequently . Therefore, we see that there is a non-trivial competition between the two terms of Eq. (3), given that by decreasing the value of increases, while decreases.

In the next section we describe the main results about the optimal control and we test them via numerical simulations.

Results

As mentioned in the last section, the optimal control is given by the shape of that minimizes both, the final number of infected individuals and the socioeconomic costs, which corresponds to maximizing the functional J from Eqs. (3) and (4). From now on we restrict ourselves to the case where the socioeconomic cost of imposing a quarantine is linear in the control . This is a simplified first approach that narrows the analysis of the general problem formulated in ”The SIR control model” section but, as we shall see, has the advantage of providing further insight into the structure of the optimal policy. The assumption that the cost of socioeconomic measures that reduce the transmission rate is linear in can be interpreted in the context of social distancing as given by Kruse and Strack[30]: “Shutting down half of the economy for two days is equally costly as shutting down the whole economy for a single day”. Then, we consider that the function L in Eq. (4) is a linear and increasing function that depends only on the control , that is , which satisfies the condition of being a monotonically increasing function of expressed in the last section. In this case, the parameter could be interpreted as the assessment that a policy maker gives to the socioeconomic impact of the quarantine compared to the final number of infected individuals. In this regard, is a fixed real number that can be chosen small enough by a government that intends to reduce the final number of infected individuals regardless the socioeconomic impact, or can be chosen large enough by a government that can face a large number of final infected and intends to control the socioeconomic impact. Here we mainly focus on the case where is small (see condition Eq. (33) in “Methods” section). Also, when is large enough we prove that the optimal strategy consists in calling off the lockdown and take mild mitigation measures for all the intervention period (see Lemma 6 in “Methods” section).

Under these conditions, we prove in “Methods” section that the optimal control is of the form of a single bang-bang. This consists on turning the strict quarantine on only once and switching it off after the maximum allowed time or, eventually, when the intervention ends at time T, depending on the initial condition and the values of the strict and mild reproduction numbers and , respectively. Then, the problem is reduced to find the optimal time to start the strict quarantine, which we call , and its length called . We also show in “Methods” section that the length of the strict quarantine for the optimal control could be less than the maximum time in some cases, as we shall see below. This means that, surprisingly, sometimes it is more convenient to make a shorter use of the strict quarantine to obtain better results in terms of pandemic costs.

We analyze the optimal control in three different case scenarios: i) , and , ii) , and and iii) and . The optimal times and for each case are given in Corollary 1, Theorems 2 and 3, respectively, of “Methods” section, where the interested reader can find rigorous proofs. The optimal control for the different cases, given by and , is summarized and numerically tested below for specific parameter values. For that, we integrate the system of Eq. (2) using the Adams’ method[35,36], for various time periods of the strict quarantine () in the bang-bang control, starting at time t but ending before T, which is the control period ( for ). That is, the value of adopts the following form, depending on t, and T:where

We start by testing the simplest case , and , and we then test the most general case and . We take the value for the recovery rate, which corresponds to a mean recovery time of 10 days that falls within the range of Covid-19 estimates[23]. This value of sets the time scales of the system. For now on all time scales are given in units of “days”, even though we omit units for the sake of simplicity.

Case , and (Corollary 1)

We first analyze the case with a bang-bang control in the interval [0, T] that consists of a mild quarantine and an extremely strict and unrealistic quarantine () during which there are no infections. The other parameter in the simulations is , together with the initial condition and . We can see from Corollary 1 that the optimum initial time of the strict quarantine is for (), while for () is given bywhere is the unique value, independent from , such that . Likewise, is the unique value, independent from , such that . As , reaches a maximum value when the strict quarantine starts at the optimal time [see Eq. (12)]. In Fig. 2 we plot vs w(0) for , calculated from Eq. (54) (squares) and by estimating the maximum of (circles). We can see that takes values close to zero for . In the rest of this section we consider the case .

Figure 2

Optimal initial time of the strict quarantine vs w(0) for and . The other parameters are and .

The behaviour of from Eq. (7) for () is tested in Fig. 3a, where we compare numerical results (circles) with that obtained from Eq. (7) (squares, Corollary 1). We observe that the agreement between numerical computations and Corollary 1 is very good. Figure 3b is an auxiliary plot that shows how to obtain graphically the optimum times and that define the three different regimes of from Eq. (7). These times are obtained by estimating the values of for which the curve crosses the lines and , which happens at and , respectively.

Figure 3

(a) Optimal initial time vs strict quarantine length for and . (b) Graphical determination of the times and that define the three regions for the different behaviours of . The parameters are and . The initial condition corresponds to (). The optimum times of the first and last regions are and , respectively, determined by and .

Remark 1

The effective reproductive number represents the mean number of individuals that an agent infects during its infectious period, at time t. It is interesting to note that the optimal time from Eq. (7) can be rewritten in terms of as

Here , and , where and are determined from the relations

In Fig. 4 we show the evolution of the system in the phase space for a given and various t (right panels), together with the evolution of (left panels), which describe the three different behaviours of . All curves start at and follow the top curve with mild quarantine () until the strict quarantine starts at t (), vertically falling down up to a lower level curve when the mild quarantine starts, and finally following this curve until the fixed point is asymptotically reached. The vertical trajectory describes the evolution within the strict quarantine where x(t) remains constant, given that in that period. The optimum time that leads to the maximum of corresponds to the time for which y(t) drops to the lowest level curve in the interval (pink curve). For (Fig. 4 top panels) we see that the maximum of is reached starting the strict quarantine at , where the effective reproduction number is , and thus there is no new outbreak when the strict quarantine is released (). In this case the entire quarantine period is used. For (Fig. 4 middle panels) the optimum initial time is , obtained by still using the entire strict quarantine period but starting earlier than . Finally, for (Fig. 4 bottom panels) the optimum is , where it turns more effective to use the strict quarantine for a shorter time . Notice that implementing a shorter but later strict quarantine is more efficient than using a longer and earlier strict quarantine, as we can see by comparing for and in the bottom panels for the case.

Figure 4

Case , and (Corollary 1). System’s trajectory in the phase space (right panels), for , , , and , and the three values of indicated in the legends corresponding to the different regimes of the optimum time (pink lines). Left panels show the time evolution of for three different initial times t of the strict quarantine in each case. The optimum times are for (top panels), for (middle panels) and for (bottom panels).

Case , and (Theorem 2)

We now analyze the case , and , with . This corresponds to a strict quarantine that is softer than in the previous case , and during which there are infections. Initially is and . We can see from Theorem 2 that the optimum initial time of the strict quarantine for iswhere is the unique value, depending on , such that . On the other hand, is the unique value, independent from , such that . The dependence and independence of w(t) on for and , respectively, can be seen from the definition of w(t) in Eq. (42).

In Fig. 5a we compare numerical results (circles) with results from Eq. (10) (squares, Theorem 2), where we see a very good agreement. At , reaches a maximum. Unlike the case, for the optimal time in the interval depends on , that is, , while for is independent of . Figure 5b,c show that the optimal times and are estimated, respectively, as the values of for which the curve crosses the horizontal line 0 and the curve .

Figure 5

(a) Optimal initial time vs strict quarantine length for and . (b) and (c) Graphical determination of the times and , respectively, which define the three regions for the different behaviours of . The parameters are . The initial condition corresponds to . The optimum time for the region is , while has a slight dependence on for .

Figure 6 is analogous to Fig. 4 for the case, and depicts the three different behaviours of . Curves are similar to those of , where the main difference is that for the trajectory of the system within the strict quarantine in the space is described by a diagonal line (see inset of top-right panel), given that is larger than zero and thus x(t) decreases in this period. At the optimum time , y(t) drops to the lowest level curve in the interval (pink curves).

Figure 6

Case , and (Theorem 2). System’s trajectory in the phase space (right panels), for , , , and , and the three values of indicated in the legends corresponding to the different regimes of the optimum time (pink lines). Left panels show the time evolution of for three different initial times t of the strict quarantine in each case. The optimum times are for (top panels), for (middle panels) and for (bottom panels).

General case and (Theorem 3)

In this section we analyze the most general case , with a mild quarantine () together with a strict quarantine () during the control interval , and with for the case of no restrictions after the control period . We take , and the rest of the parameters are the same as those in the previous studied cases. Then, from Theorem 3 the optimum initial time is given bywhere is a unique value that depends on and satisfies , while is a unique value independent of that satisfies . Here is given by Eq. (36), whereas the dependence and independence of w(t) on for and , respectively, is seen in the definition of w(t) in Eq. (35).

Given that we consider here , J reaches a maximum at the optimum time (see Eq. (12)). Figure 7a shows the behaviour of as a function of , where we observe a very good agreement between numerical results (circles) and Theorem 3 (squares). We also see that depends slightly on in the interval, while for . The optimal times and are estimated as the values of for which the curve crosses the horizontal line 0 and the curve , respectively (Fig. 7b,c).

Figure 7

(a) Optimal initial time vs strict quarantine length for and . (b) and (c) Graphical determination of the times and , respectively, which define the three regions for the different behaviours of . The parameters are , , , and . The initial condition corresponds to . The optimum time for the region is , while has a slight dependence on for .

In the right panels of Fig. 8 we show the system’s evolution in the space for three different values of corresponding to the different behaviour of . Unlike the previously studied cases where (Figs. 4 and 6), here we observe that the curves (x(t), y(t)) may exhibit up to three different regimes within the control period T, which is due to the fact that jumps three times in that interval: from to at time t, from to at and from to at T. This can be clearly seen in the curve for (inset of top-right panel of Fig. 8). For in the other two regions ( and 34), the strict quarantine ends at T for , and thus jumps twice and (x(t), y(y)) exhibits two regimes in [0, T] (insets of middle-right and bottom-right panels). As in the previously studied cases, y(t) drops to the lowest level curve in the interval for the optimum time (pink curves).

Figure 8

Case and (Theorem 3). System’s trajectory in the phase space (right panels), for , , , , and , and the three values of indicated in the legends corresponding to the different regimes of the optimum time (pink lines). Left panels show the time evolution of for three different initial times t of the strict quarantine in each case. The optimum times are for (top panels), for (middle panels) and for (bottom panels).

Discussion and conclusions

In this paper, we have studied an optimal control problem on a SIR dynamics, with a control on the reproduction number and a limitation in the duration of the intervention T and strict quarantine. Based on the Pontryagin’s maximum principle, we have given first order necessary conditions with an overall cost of the epidemic that takes into account both the maximization of the susceptible population in the long term (equivalently, a minimization of the ever infected population) and a penalization of the lockdown associated to a social and economic cost of the epidemic. We also point out that we have employed a novel proof to establish our analytical results. Moreover, some numerical examples have been provided to show the validity of our theoretical results.

Given a fixed time of intervention T where control strategies can be applied, and a strict quarantine period that represents the maximum time lapse for the stronger intervention, we proved that the optimal strategy is bang-bang when the term representing the socioeconomic cost of the objective functional is linear with respect to the control . More precisely, the optimal solution consists of switching at most twice between a mild () and a strict () quarantine, where the latter lasts at most a time period .

Although some studies have supported the idea that a too soon or too late intervention may not minimize the total mortality, we found a broader scenario. This is because the optimal solution takes the value corresponding to the lockdown on an interval , with and depending on the initial fractions of susceptible and infected individuals and , respectively, and the parameters and T. In fact, we showed that, in some cases, the optimal strategy consists of taking or (see Theorem 2 items 1 and 3-4 respectively). However, for an initial condition that corresponds to a real-life case scenario in which the percentage of the population that is infected is small when non-pharmaceutical interventions start, we obtained that the optimal strategy consists on delaying the beginning of the lockdown (items 2-4 from Theorem 2). For the case and , this optimum consists in applying a mild mitigation policy () at the beginning of the intervention, followed by a strong suppression policy () and then a mild mitigation again (mild–strict–mild strategy). Here the optimum time to start the strict quarantine corresponds to one that leaves the effective reproduction number at or just below the threshold value 1.0 when the strict quarantine is released, preventing a new outbreak. For the case the optimum corresponds to a mild–strict mitigation strategy, with a strict quarantine that starts late and lasts for a period shorter than . Surprisingly, it turns more effective to implement a short strict quarantine that starts late than a long strict quarantine that starts early.

We remark that these are optimal strategies within the basic SIR model defined in Eq. (2), which describe in an oversimplified manner the spread of an epidemic on an infinitely large population of individuals with homogeneous recovery, contagion and contact rates, where stochastic fluctuations due to finite-size effects are neglected. Then, stochastic fluctuations in the SIR model on finite populations may play a major role at the beginning of the epidemics if the fraction of infected individuals is relatively small, and thus starting with a strict quarantine may prove more effective if we want to drive the epidemic to extinction. However, we expect that the results presented in this article hold in the limit of very large populations.

We have also studied the possibility of implementing intermittent quarantines, and the possibility of applying suppression measures first, followed by mitigation measures. In both cases, if the total duration of measures is limited, we have shown that they are not optimal in order to maximize the fraction of susceptible individuals at the end of the pandemic.

A major concern with respect to the current COVID-19 crisis is the possibility of an overload of available treatment resources. Since the hospitalized individuals are a fraction of the infected population, a natural objective is to keep the number of infected individuals below some threshold for all times. In a future work we intend to extend our analytic results including a running state constraint that takes this restriction into account. It might also be interesting to study the agent-based version of the SIR model, which naturally accounts for finite-size fluctuations, in order to investigate the role played by stochastic fluctuations in the different optimal strategies described above. It would be worthwhile to explore how the results are affected by the heterogeneity in recovery and infection parameters related to age and social stratum. Finally, we also aim to study the role of an underlying network of contacts, and changes in contact rates due to individual measures triggered by fear of contagion.

Methods

Formalization of the optimal control problem

As said in “Results” section, we assume that the function L [integrand of ] depends linearly on the control . This is a simplified first approach which provides further insight into the structure of the optimal policy. In fact, under this linearity assumption we will prove that the optimal control must be bang–bang (Lemma 2), that is, the strict quarantine is turned on and off. Thus, we consider that with . In this case the functional J reads

Moreover, we consider a restriction on the admissible controls that assumes that the control can take the value for at most time, and also takes into account a maximum economic cost that the policy maker can afford. In regard to the latter, we consider a maximum cost for imposing the strict lock down for the entire period . Smaller values of represent stricter measures and thus a larger socioeconomical cost. Therefore, we impose an inferior bound to the average of on [0, T], meaning an upper bound for the socioeconomic cost. Thus, we consider the restriction

Note that any control satisfying Eq. (13) takes the value for a period no longer than . In fact, if is a control that takes the value for a longer period of time than , for instance , then we would have that contradicting the inequality from Eq. (13).

We are now in conditions to formalize the problem of finding the optimal control that maximizes the functional J. Given , fixed satisfying , , we then consider the following optimal control problem with an objective function J :

In what follows we compute the partial derivatives of with respect to x(t) and y(t) in the same way that it is done in[23]. We begin reviewing the solution of the SIR model without control Eq. (1) as done by[23]. It can be shown[1] that x(t) satisfies which combined with the identity implies thatis constant in time for any solution of Eq. (1). The trajectories in Fig. 1 are also contour lines of . Since , we then have that . Then satisfies the equation and therefore where is the principal branch of Lambert’s function[37], and thus for any

From this expression we can compute the partial derivatives of with respect to x(t) and y(t) in the same way that it is done in[23].

In order to solve problem Eq. (14) by means of the Pontryagin’s Maximum Principle, we add a new state variable given by and consider in the class of Lebesgue-measurable functions (so that we have an existence result for optimal solution). Thus, we can study the equivalent optimal control problem:

We will refer to a 4-tuple as an admissible process of the underlying control system if the control is a measurable function and the state (x, y, v) is an absolutely continuous vector function satisfying Eqs. (16b)–(16f). The optimal control problem consists in finding an optimal admissible process that maximizes the cost J. In this case, we refer to the control as optimal control.

Next, we give a result on existence of solution for the optimal control problem Eq. (16).

Proposition 1

The optimal control problem Eq. (16) admits a solution.

Proof

Since is continuous, convex and satisfies that there exists a constant such that for all it holds , the proof follows directly from Theorem 23.11[24]. Using that the control space is closed, solutions of the system of Eqs. (16b), (16c) satisfy that and the application is continuous, it is straightforward to prove conditions (a) to (f) from Theorem 23.11. Moreover, taking , we see that the unique solution of the system of Eqs. (16b)–(16d) together with gives an admissible process for which J is finite completing the hypothesis of Theorem 23.11.

The optimal control is bang-bang

In what follows, we derive the necessary conditions for problem Eq. (16) where J is given by Eq. (12). We consider the Hamiltonian Hwhere , . The necessary conditions for a maximum process on [0, T] are the following[24,25]: There exists a real number , the adjoint variable which is absolutely continuous, and such that for every t and the following conditions hold:

The adjoint variables satisfy a.e. with final time conditions (using the abbreviation for ) and the adjoint variable satisfies obtaining for all .

For a.e. Using that and defining we obtain for a.e. ,

There exists a constant C such that for a.e. thus, for all

We have the following result:

Lemma 1

The optimal control problem is normal (that is, the multiplier ).

Assume . From Eqs. (18a), (18b) with final time conditions , for all . Since the multipliers , then yielding . Therefore, from the optimality condition given in Eq. (23), a.e. contradicting the complementarity condition given in Eq. (21). Thus, we can assume , and the proof is finished.

Lemma 2

Let with and let be an optimal control of problem Eq. (16) . Then is a bang-bang control.

Assume on an interval , then computing its derivative we obtainThus, from Eq. (18a) for all and therefore for all , contradicting the end point conditions. Then, there cannot be singular arcs and the control is given by:

Lemma 3

Let be given and an admissible process. Then for and therefore

See[23].

In the next lemma we will see that the switching function changes sign at most two times, concluding that an optimal control jumps at most twice.

Lemma 4

The switching function given in Eq. (22) changes sign at most twice.

The proof follows by analysing the phase diagram of . We begin by noting that a solution of the system of Eqs. (18a), (18b) cannot cross both semilines and . This is a consequence of condition Eq. (24). In fact, assume there exist such that and . Evaluating Eq. (22) on for we have that and from Eq. (24), if or, using Eq. (25), if , both cases contradicting that and had opposite signs.

Since we have end time conditions on T we go backwards from with and . From Eq. (18a), for , , thus is decreasing and for the semiplane we have that is increasing. Also, from Eq. (18b), for we have that is increasing and for , is decreasing. Finally, from Eqs. (18a) and (18b), for and , both and are increasing.

Thus, since the end time conditions are on the region of the phase diagram with and we have that the solution backwards in time moves to the left where decreases and keeps being negative. At some point in time it could cross the semiline (note that there cannot be touch points). If the solution crosses this line it cannot cross the semiline for a previous time and thus it stays on the region for all previous times.

From the definition of Eq. (22), for , we have . For , and for , could become negative. Since , we see that for such that is on the region the function decreases for such with and increases for . Also, let such that , and , then , reaches the minimum value on at and .

Thus, we conclude that has at most two zeros on [0, T] and the proof is finished.

Theorem 1

Let be an optimal process, then:with .

Since the optimal control must be bang-bang satisfying Eq. (25), from Lemma 4 it has at most two jumps and from Eq. (13), it takes the value at most for time. Thus the proof is completed.

As a consequence of Theorem 1 we have that if is an optimal process, then the optimal control is a piecewise constant function having at most two jumps and therefore its unique associated state is a piecewise continuously differentiable function.

Characterization of the optimal control

In this section we will give the main Theorem of the article, that characterizes the switching times and (where is the beginning of the lockdown and is its duration) for an optimal control.

Let us consider the compact set (see Fig. 9)

Figure 9

Graphic of set R.

Given , for simplicity of notation, we will denote , .

Moreover, given , we will denote the solution of equation Eqs. (16b), (16c) for byand , . Then, if we call , the solution of equation Eqs. (16b), (16c) associated to the control given by equation Eq. (27) with initial data , we have that

From Theorem 1 we need to determine the maximum of the functionon the compact set R.

In order to do that, we need to compute the derivatives of J.

After some computations (see Eqs. (73) and (76) from Appendix) we obtain thatandwhere . Note that for , if and only if

In the results given in the next sections we will assume that for all and therefore in the following remark we analyse the derivatives of J restricted to the superior border of R (see Eq. (37) and (38)) which will be used later.

Remark 2

Assume that for all , then the maximum value of J on R must be attained at the superior border

Thus, in this case, for with we have and for with we have , and the control is as in Fig. 10

Figure 10

Control for In (a) we have a mild-strict-mild quarantine in the intervention interval [0, T]. In (b) we have a mild-strict quarantine in the intervention interval [0, T].

Following, we define a continuous function w for which is the second factor in equation Eq. (31) and thus gives information on the monotonicity of J (see also Eqs. (40) and (41) below).and for , we define

Then we have that for and for with

We consider the continuous function defined as the restriction of to P, that is

From Eq. (37) we have thatand from Eqs. (37) and (38), we getfor .

In the next subsection, we prove the main result of this article (Theorem 2) for the case , and . Then, we derive the result for (Corollary 1) in order to compare our result with the one obtained in[23].

Case and .

For , from Eq. (74) with , we obtain

In addition, from Eqs. (32) and (74) and using we have thatfor all .

In this case, the sign and zeros of w(t) for are the same as those for the functionwhere z(t) can be interpreted as the average on the time window of the difference between the effective reproduction number and the threshold 1.0, weighted by the inverse of . Looking at the phase diagram of Fig. 1 we see that the trajectories travel through the contour line until the strict lockdown is activated, and then descends for the time the lock down lasts (always less than ) to another contour line of with . For , the zero of w(t) () captures the moment when this travel to a lower contour line is faster, in the sense that leaves the trajectory on the lowest possible contour line at the end of the strict lockdown and, therefore, at the maximum of . The function z(t) can also be seen as an external parameter that becomes zero at the optimal time for which J reaches its maximum, i.e., . In the next remark we discuss the sign of z(t) for when .

Remark 3

Given with , assume there exists such that the solution , that is (red line in Fig. 11). Then, for , and therefore for implying . Additionally, assume there exists such that (blue line). In this case, it is clear that and also that for all there exists a unique such that . Moreover we can conclude that for , for all and therefore . If defined before does not exist, that is for all , then we take . Likewise, if does not exist, that is for all , , then we take and conclude that in either case, the sign of z(t) for , is determined in the complement of .

Figure 11

Trajectories for when .

In the next lemma we prove that z is a decreasing function on the interval introduced in Remark 3.

Lemma 5

Let be given by Remark 3. Then the function z defined in Eq. (44) is decreasing on . Moreover, for , for and consequently w changes sign at most once on .

Proof

Since the duration of the strict quarantine is fixed, for simplicity of notation in this proof we neglect the subindex from solutions x and y. Given , we define the auxiliary functions

By computing the derivative for we obtain and we have that

Note that both and have the opposite sign of . First, assume for all , then from Eq. (46a), is a decreasing function. Moreover, since and for , we deduce that . In addition, from Eq. (45c) we obtain is positive.

On the other side, from the fact that is a decreasing function, if we assume that there exists such that we have that for , and for , . Therefore, and attains a global minimum on at and thus and for all . Furhtermore, from Eq. (46b) we have that is also a decreasing function.

Thus, we have proved that for , is decreasing on and for all , yielding from Eq. (47) that for all . Finally, from Remark 3 we deduce that z changes sign at most once on (see Fig. 12) concluding that w changes sign at most once on .

Figure 12

Behaviour of z(t) when .

From Eq. (40), we see that the monotonicity of J on depends on the sign of w. For all the foregoing, if J attains its maximum at , we can interpret w(t) as an external parameter that shifts the extremal point of J from to zero.

In the next theorem we assume that . This condition is always satisfied for .

Theorem 2

Let with , and w be given by Eq. (42). Then the optimal control is unique and is given bywhere

For : and .

For and : and where is the unique value on such that .

For : and .

For : where is the unique value on such that and .

From Eq. (43) and Remark 2 the maximum value of J on R must be attained at the superior border P defined in Eq. (34). Therefore, from Eqs. (40) and (41) we haveand using that for , , thenfor . Note that from Eq. (75) with and , for it holds the identitywhereis a decreasing function in . In fact, using that is a positive function [see Eqs. (68)) and (62a)] it is easy to see thatand therefore, for . From Eqs. (50) and (51) we also have that for

We consider the following cases:

If , then from Lemma 5, for all . Thus, from Eq. (49) we have Moreover, using that , the positivity of and Eq. (51), we obtain that and being h a decreasing function, from Eq. (53) we deduce that Therefore .

Since and , from Lemma 5 there exists an unique such that , for all and for all .

Moreover, since , in the same way as for the previous item, we have and from Eq. (49), we obtain concluding that .

Since from Lemma 5, for all . On the other side, since , from Eq. (52) we have that and using that h is a continuous and decreasing function, we obtain that for all . Thus, Therefore, and .

Since from Lemma 5, we have for all . On the other side, since , then and using that and h is a continuous and decreasing function, there exists a unique such that , for and for . Therefore, Consequently, and .

Case , and .

Let w(t) defined as in Eq. (35). Note that for and , from Eq. (75) with , we obtain

It is easy to observe that in this case the sign of w(t) on is given by . Moreover, w changes sign at most once on , going from positive to negative values.

Corollary 1

Let , and . Then, the optimal control is unique and is given bywhere

For : and .

For and : and , where is the unique value on such that .

For : and .

For : and , where is the unique value on such that .

The proof follows from Theorem 2 using the fact thatand

Note that in this case if we take , the corollary is reduced to only two possible cases: or , obtaining the same result as Ketcheson[23] in Theorem 3.

General case

In this section we study the behaviour of optimal solutions for the general case when and , that is, objective function J includes the term that accounts the running cost of the control and allows us to account for factors like the economic cost of intervention or heightened risks caused by hospital overflow.

Lemma 6

Assume for all , then the optimal control is given by .

From Eq. (32), we have that and therefore the maximum value of J on R is attained at the inferior border of R where and is constant.

In the next theorem we give a general result including both the economic cost of intervention () and a mitigation phase different from the no intervention one, that is . In the next section we give numerical simulations supporting this result. When and we recover Theorem 2.

Theorem 3

Let with , satisfying Eq. (33) for all and let w and defined as in Eqs. (35) and  (36) respectively , then the optimal control is unique and is given bywhere

For : and .

For and : and where is the unique value on such that .

For : and .

For : where is the unique value on such that and .