key: cord-0803420-yldlknad
authors: Borri, Alessandro; Palumbo, Pasquale; Papa, Federico; Possieri, Corrado
title: Optimal design of lock-down and reopening policies for early-stage epidemics through SIR-D models()
date: 2020-12-23
journal: Annu Rev Control
DOI: 10.1016/j.arcontrol.2020.12.002
sha: 9dc336e0177199b567e6d622169c3fd4b9ce3946
doc_id: 803420
cord_uid: yldlknad

The diffusion of COVID-19 represents a real threat for the health and economic system of a country. Therefore the governments have to adopt fast containment measures in order to stop its spread and to prevent the related devastating consequences. In this paper, a technique is proposed to optimally design the lock-down and reopening policies so as to minimize an aggregate cost function accounting for the number of individuals that decease due to the spread of COVID-19. A constraint on the maximal number of concomitant infected patients is also taken into account in order to prevent the collapse of the health system. The optimal procedure is built on the basis of a simple SIR model that describes the outbreak of a generic disease, without attempting to accurately reproduce all the COVID-19 epidemic features. This modeling choice is motivated by the fact that the containing measurements are actuated during the very first period of the outbreak, when the characteristics of the new emergent disease are not known but timely containment actions are required. In fact, as a consequence of dealing with poor preliminary data, the simplest modeling choice is able to reduce unidentifiability problems. Further, the relative simplicity of this model allows to compute explicitly its solutions and to derive closed-form expressions for the maximum number of infected and for the steady-state value of deceased individuals. These expressions can be then used to design static optimization problems so to determine the (open-loop) optimal lock-down and reopening policies for early-stage epidemics accounting for both the health and economic costs.

The fast diffusion of the COVID-19 outbreak, as well as the severe consequences that such a pandemic disease produces on the national health system, are common problems that the whole world is currently coping with. In the last two decades, the frequent emergence of new diseases, like the respiratory syndromes SARS (2002) and MERS (2012) , and the re-emergence of new foci of older diseases, like ebola (2007, 2008, 2012) and zika (2010) viruses, or cholera (2010), have been real threats for the entire world because of the high risk of a fast disease spread due to the globalization and the highly interconnected network among people. However, those epidemics have not reached such an extensive spread like COVID-19, that can be considered as the next world wide pandemic after the Spanish flu of 1918.

Nowadays, after almost six months since the first notified cases coming from the Wuhan region in China (Surveillances, 2020) , the governments of the entire world are still adopting containment measures to stop the pandemic spread and facing its devastating consequences on the health and economic systems of their countries. To this end, these efforts have aroused the attention of the scientific community on the study of mathematical models of pandemic evolution and their ✩ The material in this paper has not been presented at any conference. * Corresponding author.

to properly identify the model parameters in a model-based design of efficient control strategies. In summary, dealing with an early-stage epidemics and the necessity of readily identifying the main features of the epidemic model, the simplest choice provides smaller confidence intervals for the model parameters, especially in the present COVID-19 outbreak (Bertozzi, Franco, Mohler, Short, & Sledge, 2020; Roda, Varughese, Han, & Li, 2020) . Indeed, an important aspect in the mathematical modeling of COVID-19 is also the inclusion in the model structure of the interventions adopted for preventing and containing the outbreak, because actually some medical, social and political actions were progressively actuated since the first notifications (Gatto et al., 2020; Giordano et al., 2020) .

Many other papers, especially published at the beginning of the outbreak, are focused on the evaluation of the main features of the disease evolution and on the estimation of the transmission risk and of the basic reproduction number of the disease, on the basis of the initial available epidemiological data (Li, Guan, Wu, Wang, Zhou, et al., 2020; Tang, Bragazzi, Li, Tang, Xiao, & Wu, 2020; Tang, Wang et al., 2020a; Wu, Leung, & Leung, 2020; Zhang, Diao, Yu, Pei, Lin, & Chen, 2020; Zhao et al., 2020) .

As previously stated, on spite of the high degree of accuracy reached by complex modeling frameworks in reproducing the pandemic data, things look different from the perspective of model-based feedback design, in which most of these models result to be too complex in view of the construction of control/optimization policies, as well as for the study of the qualitative properties of the resulting system. This imposes a trade-off between simplicity and accuracy.

It is worth pointing out that, in the literature, several techniques have been proposed to optimally design open-loop and closed-loop lock-down policies. Some examples in the rapidly expanding literature include (Köhler et al., 2020) , where the SIDARTHE model has been used to design a feedback strategy, based on the model predictive control (MPC) method, to minimize the fatalities due to COVID-19, (Carli, Cavone, Epicoco, Scarabaggio, & Dotoli, 2020) where an optimal control approach based on the SIRQTHE model has been proposed to determine the most effective strategies to be adopted during postlockdown mitigation phases in a multi-region scenario, (Morato, Bastos, Cajueiro, & Normey-Rico, 2020) where an MPC policy is formulated to mitigate the COVID-19 contagion in Brazil, designed as optimal onoff social isolation strategy, and (Gonzàlez, Torfs, Nopens, & Alleman) where the SEIR model has been used to predict and control the spread of COVID-19. Differently from these works wherein complex models have been used to predict the evolution of the COVID-19 pandemic, in this paper the transmission of COVID-19 is modeled by means of a simple SIRD model. The rationale behind this choice is that such a model has few parameters that can be identified by using a short interval of observation, so to immediately take the necessary actions to limit the spread of COVID-19. Further, thanks to its simple structure, the solution to the SIRD model can be computed exactly, thus allowing us to exploit explicit expressions to carry out the optimization task.

In the present paper an optimization of the lock-down strategy, the common measure adopted by the majority of the countries, is proposed. The optimization approach aims at containing the outbreak during the very first spreading period, when the characteristics of a new emergent disease are not known but timely containment actions are required. For this reason, we formulated the problem on the basis of a simple dynamical ODE model describing the single outbreak of a general disease, without attempting to accurately reproduce all the features of the new COVID-19. However, the model parameters are assessed on the basis of official epidemiological data related to the first 14 days of the Italian COVID-19 pandemic case (Rosini, 2020) . The chosen model is a modified version of the basic SIR model that explicitly distinguishes between recovered and death (it will be referred to as SIRD). The modeling choice allows to get a good fitting of the initial outbreak, according to which the model-based strategy will be applied. Lock-down and release events act in terms of reset of the susceptible population, so that the resulting system is hybrid (Goebel, Sanfelice, & Teel, 2009) .

The paper is organized as follows. Section 2 introduces the modeling framework and studies some fundamental properties of the SIRD model, which are useful in the sequel. Section 3 addresses the optimal design of lock-down, while Section 4 addresses the optimal design of joint lock-down and reopening. Section 5 illustrates some simulations of the proposed methods based on real data. Section 6 offers concluding remarks.

Consider the following modified version of the McKendrick and Kermack (compartmental) model (Bailey, 1975; Kermack & McKendrick, 1927) accounting for the number of deceased individuals due to the disease:

with initial conditions (0) = 1 > 0, (0) = 2 > 0, (0) = 3 ≥ 0, (0) = 4 ≥ 0, where > 0 is the infection rate, > 0 is the recovery rate, and > 0 is the mortality rate. The compartments used for this model consist of four classes:

• individuals ( ) not yet infected, at time , and susceptible to be infected; • individuals ( ) that have been infected at time and are capable of spreading the disease to those in the susceptible category; • individuals ( ) not susceptible to the disease, e.g. because already infected and recovered; • the individuals ( ) that deceased due to the disease.

In this Section, we compute three first integrals for the SIRD model given in Eq. (1). These functions are used in the subsequent Section 2.2 to compute its solution and in Section 2.3 to determine the limit of such a solution for → +∞.

The following proposition introduces the first integrals of the SIRD model.

be the vector field governing its dynamics. Then, functions 1 , 2 , and 3 defined as

are first integrals of the SIRD model (1), since their gradient

is a vector orthogonal to (i.e. ∇ ⋅ = 0, = 1, 2, 3). Hence, they assume a constant value along the solutions of system (1), that is Note that has the clear meaning of the number of individuals of the overall population.

According to Proposition 1, one has that susceptibles and infected are constrained to evolve onto the following trajectory of the ( , ) phase plane:

( ) + ( ) − + log( ( )) = .

(4)

The first integrals given in Proposition 1 can be used to compute the solution to the SIRD model (1) as we show in the following theorem.

Then, the solution to the SIRD model (1) is given by

Proof. By considering Eq. (1a), we have that

By substituting Eq. (7) in Eq. (4), we obtain the following first order ordinary differential equation in :

Letting ( ) = log( ( )) and substituting it into Eq. (8), we obtain the following first order ordinary differential equation in :

According to Eq. (5), it comes that = −1 ( ) is the solution to Eq. (9), sincė

Hence, by considering that ( ) = exp( ( )), and using the first integrals given in Proposition 1, we obtain the expressions given in Eq. (6). □

It is worth pointing out that an alternative approach to compute the solution to system (1) has been given in Harko, Lobo, and Mak (2014) by using different arguments involving the solution to a second order differential equation and a slightly different parameterization for ( ).

( ) is the following:

• if 1 > ( + )∕ , then ( ) will initially increase at = 0, then come to a maximum for

and eventually converge to zero for ↦ +∞. • if 1 ≤ ( + )∕ , then ( ) will be monotonically decreasing, eventually vanishing for ↦ +∞.

Proof. As a preliminary remark, it is apparent that the state evolution of the 4 state variables is kept non-negative for any non-negative initial conditions. According to the constraint provided by the first integral in Eq. (4), we can write as the following function of

and a gradient ∕ providing a monotonically increasing phase for ∈ (0, ( + )∕ ) and monotonically decreasing phase for ∈ (( + )∕ , +∞). Because of the non-negative initial conditions of both susceptible and infected individuals, the unique maximum for = ( ), achieved for = ( + )∕ , is ensured to be positive. Besides, from Eq. (1a) it comes that ( ) has a monotonically decreasing fashion, so that the relative initial position of (0) = 1 with respect to the argument of the maximum of ( ) decides whether evolution has an initial increasing period or not, according to the Proposition theses. By inverting Eq. (6a), with = ( + )∕ , the formula readily comes. □ In summary, from Proposition 2, the maximum value attained by the state response in the variable of the SIRD model (1) is

The first integrals given in Proposition 1 can be used also to characterize the limit for → +∞ of the solution to the SIRD model (1), as shown in the following proposition.

with W 0 (⋅) being the principal branch of the Lambert omega function, defined as the inverse of ( ) = (Rudin, 1964) .

Proof. By setting equal to zero the dynamics of the SIRD model (1), it readily comes that stationary points are provided by a manifold, sharing the common feature that ∞ = 0. Therefore, the limit for ↦ +∞ of the other 3 state variables depends on the initial conditions. Substituting ∞ = 0 in Eq. (4) and taking into account the monotonically decreasing behavior of , it comes out that the remaining number of susceptible individuals at the end of the epidemics ∞ is given by the smallest solution to the transcendental equation

By further manipulating Eq. (14) it comes that

so that, by exploiting the Lambert omega function notation, the remaining number of susceptible individuals at the end of the epidemics is the one given in Eq. (13a). Hence, using the remaining first integrals given in Proposition 1, we obtain the expressions given in Eq. (13). □

In this Section, we use the results gathered in Section 2 to optimally design a lock-down. A lock-down may be modeled by removing a number of individuals from the susceptible compartment, with the goal of reducing the number of people that get infected by the disease: supposed to be confined at home, they become no longer susceptible to be infected. Our objective is to design the number together with 

Behavior of ( ) if 0 < ≤ . Descending arrows represent a decaying behavior for ( ), whereas ascending arrows represent an increasing behavior for ( ). The same notation is used in all the figures reported hereafter. the time ≥ 0 at which the lock-down takes place so to minimize the number of deaths due to the disease and, at the same time, to keep the number of infected individuals below a given threshold value. By Theorem 1, the effect of removing individuals from the susceptible compartment at time ≥ 0 results in the new (post-removal) initial conditions ( 1 ),

From these new initial conditions, the system evolves in free response with the updated first integrals ′ ( , ) = 1 ( , ) + 2 ( ) − + log( 1 ( , )),

Hence, by using Proposition 3 with the updated constant, the number of individuals that decease due to the epidemics is

We now characterize the maximum number of infected individuals considering the effects of the lock-down. Since the evolution of susceptible individuals is monotonically decreasing, we have 1 ( , ) ≤ exp( −1 ( )) ≤ 1 for all ≥ 0 and ≥ 0. Hence, if the infected evolution is decreasing at the lock-down time , then it will keep decreasing also after the lock-down. This straightforwardly comes from the computation of the time derivative of written according to Eq. (11):

and from the consideration thaṫ( ) ≤ 0 for any ≥ 0 and that ( + ) < ( − ). In order to streamline the exposition, denote

( 1 , 2 ) ∶= 1 + 2 − + + + log (

where, by the construction given in Section 2, Eqs. (17) and (18) correspond to the maximum value attained by the state response of system (1) in the variable with initial condition (0) = 1 , (0) = 2 , and to the time at which such a maximum is attained, respectively, if the lock-down does not take place; Eq. (19) corresponds to the updated time at which the maximum is attained with the post-removal initial conditions. Hence, consider the following cases.

• If ≤ 0 (that is 1 ≤ + ) then, by the construction given at the end of Section 2.2, the time behavior of the state response in the variable ( ), before and after the lock-down, is the one given in Fig. 1 .

Thus, in such a case, the maximum number of infected individuals is given by

• If 0 < ≤ , then the time behavior of ( ) is the one given in Fig. 2 . Thus, in such a case, the maximum number of infected individuals is given by

• If > and ′ ( , ) < 0, then the time behavior of ( ) is the one given in Fig. 3 . Thus, in such a case, the maximum number of infected individuals is given by

• If > and ′ ( , ) ≥ 0 then, the time behavior of ( ) is the one given in Fig. 4 . Thus, in such a case, the maximum number of infected individuals is given by

Hence, by the analysis carried out above, the maximum number of infected individuals is given by

Therefore, by considering that, in view of Theorem 1, the maximum number of individuals that can be removed at time ≥ 0 from the susceptible compartment is

letting ( , ) be the economic cost of removing individuals at time (an example for such a cost is given in the subsequent Section 5) and letting be the maximum admissible number of individuals in the infected compartment, the problem of optimally designing the lockdown is equivalent to determine a solution to the following constrained optimization problem

The constraint on the maximum number of infected is motivated by the economic burden on public (and private) health services and hospitals, in order to prevent the collapse of the health system.

It is of interest to determine the solution to the optimization problem (20) in the case that ( , ) = 0 and = +∞, i.e. without accounting for the economic cost of the lock-down and for the burden on the health care. By Eq. (1a) and Proposition 3, the function ↦ exp( −1 ( )) is monotonically decreasing with

Hence, letting ∞ = − + 0 (− 1 + exp(− ( 1 + 2 ) + )) and defining the change of variables

the problem of minimizing ( , ) subject to ≤ ( ) is equivalent to solve the following problem

wherē( , ) is obtained from ( , ) using the substitutions in Eq. (21). By computinḡ( , ) , one obtains that such a derivative is non-positive within the feasible domain. Using a wholly similar reasoning with respect to the variable , it can be easily derived that the solution to the problem given in Eq. (22) is = 1 and = 1. Not surprisingly, it corresponds to = 0 and = 1 , i.e., the optimal lock-down policy that minimizes the number of deaths without accounting for the welfare cost is to lock-down all the susceptible individuals at the beginning of the epidemic. Note that, with such a choice, the maximum number of infected individuals is 2 . Clearly, this analysis is not realistic since one has to consider also the economic cost and the number of beds available in intensive care units when planning a lock-down. These real cases will be discussed in details in Section 5.

One of the key assumption in Section 3 for the optimal design of the lock-down is that individuals are removed from the susceptible compartment at time ≥ 0 up to the end of the epidemics. The main goal of this Section is to also design the number of individuals 2 and the time 2 at which such individuals return in the susceptible compartment due to the reopening.

The construction given in Section 3 to determine the post lock-down initial condition still holds in the setting of this Section. Hence, letting By Theorem 1, the response of system (1) at ∈ ( 1 , 2 ) is

and the initial (post-reopening) conditions are given by

( 2 − 1 ) + log( 1 ( 1 , 1 )), ( + 2 ) = ℎ 4 ( 1 , 1 , 2 ) ∶= 4 ( 1 ) − ′−1 1 , 1

( 2 − 1 ) + log( 1 ( 1 , 1 )).

From these initial conditions, the system evolves freely up to the end of the epidemics, with the updated constants ′′ = ℎ 1 ( 1 , 2 , 1 , 2 ) + ℎ 2 ( 1 , 1 , 2 ) − + log ( ℎ 1 ( 1 , 2 , 1 , 2 ) ) , ′′ = ℎ 3 ( 1 , 1 , 2 ) − ℎ 4 ( 1 , 1 , 2 ), ′′ = ℎ 1 ( 1 , 1 , 2 ) + ℎ 2 ( 1 , 1 , 2 ) + ℎ 3 ( 1 , 1 , 2 ) + ℎ 4 ( 1 , 1 , 2 ).

Hence, by Proposition 3, the total number of people that decease due to the disease is

( 1 , 2 , 1 , 2 ) ∶= ℎ 4 ( 1 , 1 , 2 ) + + (ℎ 1 ( 1 , 2 , 1 , 2 ) + ℎ 2 ( 1 , 1 , 2 ))

) ) .

On the other hand, following the same construction used in Section 3, we have thaṫ

In order to streamline the exposition, define ( 1 , 2 ) and as in Section 3, and

where, by the construction given in Section 2, Eqs. (24) and (25) correspond to the time at which the state response of system (1) in the variable attains its maximum starting from the post lock-down and reopening initial conditions, respectively. Hence, consider the following cases.

• If ≤ 0 and ′′ ( 1 , 2 , 1 , 2 ) ≤ 0, then, by the construction given in Section 3, the time behavior of the state response in the variable ( ) is the one given in Fig. 5 . Thus, in such a case, the maximum number of infected individuals is given by

• If ≤ 0 and ′′ ( 1 , 2 , 1 , 2 ) > 0, then the time behavior of ( ) is the one given in Fig. 6 . Thus, in such a case, the maximum number of infected individuals is given by

( 1 , 2 , 1 , 2 ) = max{ 2 , (ℎ 1 ( 1 , 2 , 1 , 2 ), ℎ 2 ( 1 , 1 , 2 ))}.

• If 0 < ≤ 1 and ′′ ( 1 , 2 , 1 , 2 ) ≤ 0, then the time behavior of ( ) is the one given in Fig. 7 . Thus, in such a case, the maximum number of infected individuals is given by

• If 0 < ≤ 1 and ′′ ( 1 , 2 , 1 , 2 ) > 0, then the time behavior of ( ) is the one given in Fig. 8 . Thus, in such a case, the maximum number of infected individuals is given by ( 1 , 2 , 1 , 2 ) = max{ ( 1 , 2 ), (ℎ 1 ( 1 , 2 , 1 , 2 ), ℎ 2 ( 1 , 1 , 2 ))}.

• If > 1 , ′ ( 1 , 1 ) ≤ 0 and ′′ ( 1 , 2 , 1 , 2 ) ≤ 0, then the time behavior of ( ) is the one given in Fig. 9 . Thus, in such a case, the maximum number of infected individuals is given by ( 1 , 2 , 1 , 2 ) = 2 ( 1 ).

• If > 1 , ′ ( 1 , 1 ) ≤ 0 and ′′ ( 1 , 2 , 1 , 2 ) > 0, then the time behavior of ( ) is the one given in Fig. 10 . Thus, in such a case, the maximum number of infected individuals is given by ( 1 , 2 , 1 , 2 ) = max{ 2 ( 1 ), (ℎ 1 ( 1 , 2 , 1 , 2 ), ℎ 2 ( 1 , 1 , 2 ))}.

• If > 1 , 0 < ′ ( 1 , 1 ) ≤ 2 − 1 and ′′ ( 1 , 2 , 1 , 2 ) ≤ 0, then the time behavior of ( ) is the one given in Fig. 11 . Thus, in such a case, the maximum number of infected individuals is given by

( 1 , 2 , 1 , 2 ) = ( 1 ( 1 , 1 ), 2 ( 1 )).

• If > 1 , 0 < ′ ( 1 , 1 ) ≤ 2 − 1 and ′′ ( 1 , 2 , 1 , 2 ) > 0, then the time behavior of ( ) is the one given in Fig. 12 .

Thus, in such a case, the maximum number of infected individuals is given by

( 1 , 2 , 1 , 2 ) = max{ ( 1 ( 1 , 1 ), 2 ( 1 )), (ℎ 1 ( 1 , 2 , 1 , 2 ), ℎ 2 ( 1 , 1 , 2 ))}.

• If > 1 , ′ ( 1 , 1 ) > 2 − 1 , then the time behavior of ( ) is the one given in Fig. 13 . Thus, in such a case, the maximum number of infected individuals is given by

( 1 , 2 , 1 , 2 ) = (ℎ 1 ( 1 , 2 , 1 , 2 ), ℎ 2 ( 1 , 1 , 2 )).

Hence, using the function defined above, considering that the maximum number of individuals that can be removed at time 1 ≥ 0 from the susceptible compartment is ( 1 ), letting ( 1 , 2 , 1 , 2 ) be the economic cost of removing 1 individuals at time 1 and releasing 2 ≤ 1 individuals at time 2 ≥ 1 , and letting be the maximum admissible Fig. 7 . Behavior of ( ) if 0 < ≤ 1 and ′′ ( 1 , 2 , 1 , 2 ) ≤ 0. Fig. 8 . Behavior of ( ) if 0 < ≤ 1 and ′′ ( 1 , 2 , 1 , 2 ) > 0.

Behavior of ( ) if > 1 , ′ ( 1 , 1 ) ≤ 0 and ′′ ( 1 , 2 , 1 , 2 ) ≤ 0. Fig. 10 . Behavior of ( ) if > 1 , ′ ( 1 , 1 ) ≤ 0 and ′′ ( 1 , 2 , 1 , 2 ) > 0. Fig. 11 . Behavior of ( ) if > 1 , 0 < ′ ( 1 , 1 ) ≤ 2 − 1 and ′′ ( 1 , 2 , 1 , 2 ) ≤ 0.

Behavior of ( ) if > 1 , 0 < ′ ( 1 , 1 ) ≤ 2 − 1 and ′′ ( 1 , 2 , 1 , 2 ) > 0. number of individuals in the infected compartment, the problem of optimally designing the lock-down and the reopening is equivalent to determine a solution to the following optimization problem

Remark 2. Following the same construction used in Remark 1, it can be easily derived that the optimal lock-down and reopening policy that minimizes the number of deaths without accounting for the welfare cost is to lock-down all the susceptible individuals at the beginning of the epidemic and release them when it ends.

In this Section, we show some simulations of the proposed containment strategy and solve the optimization problems described in the previous Sections based on real data.

In order to obtain realistic parameters for the optimization problems, the SIRD model (1) The identified parameters have been used to determine the best and worst case epidemic scenarios in terms of number of deceased individuals. Namely, the formula given in Eq. (13d) has been used to find the best-case and worst-case set of parameters , , , leading to the lowest and highest number of deaths, respectively, if lock-down measurements were not taken. Fig. 15 depicts the number of deceased individuals due to the epidemics letting the parameters , , and vary within the 99% confidence intervals identified above.

As shown by such a figure, thanks to the reduced confidence interval of parameter estimates, the uncertainty related to the number of asymptotic deceased individuals is mainly due to the and uncertainties. In more detail, without taking any control action, the parameters that led to the lowest number of deceased individuals are 

These parameters and the fitted initial conditions have been used to solve the optimization problem (20) in the best and worst case scenarios to optimally design the lock-down letting = 7 ⋅ 10 4 , that is the maximum admissible number of infected individuals, and letting 

where ∈ (0, 1] is a positive dimensionless scalar used to account for the percentage of individuals that are negatively affected by the lock-down from an economic viewpoint (in the following, we have assumed that = 0.75). The Hill function modulates such economic suffering with respect to the delay , making it lighter and lighter by increasing the lock-down delay , with the halving time set at = 1 = 4 days: if the lock-down is delayed, then there is a lower economic and psychological cost thanks to the preparedness of the population. 1 regulates the steepness of the curve (in the following, 1 = 6). Note that both the health and economic costs building the optimization function are expressed in number of individuals.

Parameter has been designed considering the number of intensive care units in Italy (over 7000 units, (Ministero della Salute)), and the fact that after 2 weeks of the COVID-19 pandemics in Italy the number of infected people treated in intensive care was around the 10% of the total number of infected people.

The solution to the optimization problem (20) in the worst case scenario. As a matter of fact, the two optimal solutions, though achieved according to opposite frameworks, are very similar, especially with respect to the number of quarantined people . Fig. 16 depicts the results of two numerical simulations showing the time behavior of the number of susceptible, infected, recovered, and dead individuals according to the proposed optimal solution designed for both the best and the worst scenario.

As shown by Fig. 16 , the determined lock-down policies guarantee that the number of infected individuals is far below the threshold , i.e. the constraint on the maximum number of infected is not A. Borri et al. active for the optimal policy. Furthermore, note that both policies are rather similar to the one that has been implemented by the Italian government, where substantially all individuals have been forced to lock-down. Also the beginning of the lock-down is actually close to the real intervention time in Italy.

In order to evaluate the balance between the health and the economic cost in the considered problem, we also investigated what happens by varying the weight of the aforementioned terms of the cost function, resulting in a modified cost function (1 − ) ( , ) + ( , ), ∈ (0, 1), subject to the same constraints considered in problem (20). Fig. 17 depicts the optimal values of and for different values of ∈ (0, 1).

On the other hand, Fig. 18 depicts the health and economic component of the optimal cost corresponding to different values of .

As shown by Fig. 17 , the optimal value ⋆ is anchored to a very high value that substantially resembles a total lock-down; on the other hand, the optimization algorithm modulates more finely the timing ⋆ of the lock-down, which reasonably increases with . However, the difference between the two values of ⋆ , related to the best and worst scenarios, is almost constant with showing a delay of almost 5 days between the two scenarios. This result confirms the robustness of the proposed procedure that indicates a clear time of intervention. Furthermore, in both the best and worst case scenarios the health component of the overall cost increases with . Clearly, the health and economic costs are larger in the worst case scenario.

A. Borri et al. It is worth noticing that the analysis carried out in Section 3 can be used also to evaluate the effects of a soft lock-down. In fact, by inspecting Figs. 3 and 4, one has that, according to the non-trivial case of > , the spread of the epidemics is suddenly stopped (with the infected starting to definitively decrease at the very beginning of the lockdown) only if ′ ( , ) ≤ 0. For a given parameter setting and lock-down time , such a condition is satisfied if and only if ≥̃, wherẽ = ( ) − ( + )∕ .

The above relations can be easily inferred from the integration limits of Eq. (19), recalling that its integrating function is non-positive (aṡ =̇∕ ≤ 0, and susceptibles can only decrease).

In particular, by considering the parameter setting of the best case scenario and the related optimal time ( = ⋆ ), we obtaiñ= 50.936⋅10 6 individuals (corresponding to 84.39% of the total population). On the other hand, considering the worst case scenario, we havẽ= 51.960⋅10 6 individuals (corresponding to 86.08% of the total population). In other words, even according to the best scenario framework, any soft lockdown consisting of less than 51 million people in quarantine would not prevent the infected population to keep on increasing soon after the lock-down policy.

We also tested the optimal lock-down and reopening policy designed according to the analysis carried out in Section 4.

A. Borri et al. Fig. 19 . Results of two numerical simulations in which the optimal lock-down and reopening policies corresponding to the worst and best case scenarios have been implemented. Fig. 19 depicts the results of two numerical simulations showing the time behavior of the number of susceptible, infected, recovered, and dead individuals according to the proposed optimal solution designed for both the best and the worst scenarios.

As shown by Fig. 19 , the determined lock-down and reopening policies guarantee that the number of infected individuals is far below the threshold , that means, even in this case, the constraint on the maximum number of infected individuals is not active. Furthermore, note that the optimal number of individuals removed from the susceptible compartment and the timing of the lock-down is rather similar to that obtained solving problem (20).

Similarly to the lock-down policy, what comes out by the inspection of the optimal solutions provided in the best and worst scenarios of the lockdown-and-reopening policy is a clear suggestion of an extraordinarily strict lock-down (quite the whole population for both cases) with a slightly different delay in the onset of the lock-down policy (again, about 4 days); similarly, the reopening is suggested to be very mild for both scenarios, ranging from 2.5 to 3 million individuals released after about 3.5 months (with a displacement of about 5 days).

It could be interesting to check out what happens when applying a chosen policy to the wrong scenario. For instance, what happens if we apply the optimal policy computed for the best case to the worst scenario? Do few days of delay really count providing a substantially different evolution of the disease? Results are reported in Figs. 20-21. The histograms show the cost function computed for the worst case according to the optimal policy and to the wrong policy (the one achieved assuming the best epidemic scenario). Results tell us that even a few days of delay in the application of the lock-down policy do count in increasing the overall number of deceased individuals, since such a number is more than twice as much as in the optimal case. Dealing with the SIR mathematical model, this is not surprising: the mortality rate showed to have the largest relative uncertainty within its confidence interval; besides, according to the approximately exponential increase of the early-stage and to the worst-case scenario parameters, 4 days may well more than double the number of infected population. On the contrary, in the best scenario, using the wrong policy (i.e., the one determined considering the worst scenario) leads to an economic cost that is more than three times the optimal one. Also in the lockdown-and-reopening control problem, we investigated what happens by varying the weight of the healthy and economic terms of the cost function, by properly varying parameter ∈ (0, 1) in the modified cost (1 − ) ( 1 , 2 , 1 , 2 ) + ( 1 , 2 , 1 , 2 ), subject to A. Borri et al. the same constraints considered in problem (26). Fig. 22 depicts the optimal values of 1 , 2 , 1 , and 2 for different values of ∈ (0, 1). On the other hand, Fig. 23 depicts the health and economic component of the optimal cost corresponding to different values of .

As shown by Fig. 22 , the optimal value ⋆ 1 is increasing with (while the difference between the best and worst case is almost constant), the optimal values ⋆ 1 and ⋆ 2 are decreasing with , whereas ⋆ 2 is essentially constant with . Note that, according to the optimal procedure, the reopening is allowed after about four months for the majority of the values of (see the bottom panel of Fig. 22 , where almost all the dots are under about 120 days, except for low alpha values of the worst case). Further, as in Fig. 18 , in both the best and worst case scenarios the health component of the overall cost increases with . Clearly, the health and economic costs are larger in the worst case scenario.

The optimal solution obtained solving the optimization problem (26) shows a critical point of the problem formulation: substantially all individuals are forced to lock-down and the reopening allows only the reintroduction of a minimal part (about 5%) of the whole population. This ''strong'' solution is mainly due to two main factors: the modeling assumption of a constant relative infection rate and the strong dependence of the infected population on such a crucial parameter, which straightforwardly implies a strong impact on the maximal number of concomitant infected patients and on the final number of deaths. So the reintroduction of a substantial number of susceptibles in would produce a second wave that could violate the constraint on the maximal number of concomitant infectives and would produce a very high number of deaths at the end of the pandemic. This behavior is caused by the absence of a time-varying dependence of in the model formulation. In other words, we do not take into account for the ''realistic'' possibility of the susceptible population to learn ''something'' on the transmission mechanisms of the new infection during the lock-down time interval. Conversely, a realistic scenario would reasonably imply a change of during the lock-down, due to the increasing awareness of the population on the transmission risk and on the safety behaviors. Such a critical aspect of the problem is highlighted in Fig. 24 , where the strong dependence of the number of deaths on both the parameter and the number of released susceptibles 2 has been shown for both the best (left panel) and worst case (right panel).

Based on such observations, it is now interesting to evaluate how many people would be allowed to be released assuming a different awareness among the susceptible population, that is a different . Fig. 25 shows this interesting aspect by reporting the behavior of the function 2 = ( ) obtained solving Eq. ( * 1 , 2 , * 1 , * 2 ) = * w.r.t. , 2 , where * 1 , * 1 , * 2 are the optimal quantities obtained solving problem (26) (for both the best and worst scenarios) and * is the number of deaths as resulting from the optimization procedure. In other words, Fig. 25 tells us about how much can we improve the strength of the release if a reduction of the infection rate is allowed. It is important to notice that a change in the population behavior would produce a strong impact on the policy planning, which is a huge impact on the number of released people 2 . Indeed, a reduction of of almost 10 A. Borri et al. , 2 , * 1 , * 2 ) = * ; * 1 , * 1 , * 2 , * : optimal values obtained solving problem (26) for the best scenario ( * 1 , * 1 , * 2 , * ) or the worst scenario ( * 1 , * 1 , * 2 , * ).

times would produce an increase of 2 up to about 50-60 millions of released people. Fig. 26 shows the result of the optimal problem (26), for both best and worst scenarios, with respect to the decision variables 2 , 2 when varies over time as: ( ) = , , for < 2 , and ( ) = ∕10, ∕10, for ≥ 2 . The simulations are obtained by fixing also the remaining decision variables to the optimal values of Fig. 19 ( 1 = * 1 , * 1 and 1 = * 1 , * 1 ). The figure shows that, assuming a change in the population behavior that produces a reduction of the initial by a factor of 10, the reintroduction of almost all the population in is allowed.

It is worth noticing that the analysis carried out in Section 4 can be used also to evaluate the admissible values of the infection rate that avoid a second epidemic wave due to the reopening. In fact, by inspecting Figs. 6, 8, 10, 12, and 13, one has that there is a second wave due to the reopening only if ′′ ( 1 , 2 , 1 , 2 ) ≥ 0. More in details, denoting by 2 the relative infectivity rate reflecting the population behavior after the reopening, i.e. for ≥ 2 , we want to find 2 such that ′′ ( 1 , 2 , 1 , 2 ; 2 ) ≤ 0. So, with a slight abuse of notation, rewriting Eq. (25) as ′′ ( 1 , 2 , 1 , 2 ; 2 ) = ∫ log( + 2 ) log(ℎ 1 ( 1 , 2 , 1 , 2 ; )) 1 2 − 2 ′′ ( 1 , 2 , 1 , 2 ; 2 )− − d ,

it is easy to conclude that ′′ ( 1 , 2 , 1 , 2 ; 2 ) ≤ 0 if and only if 2 ≤̃2, wherẽ

In this work, we proposed a method for the joint optimal lock-down and release design in a pandemic. The approach is based on a SIRD formulation, which provides a satisfactory compromise between accuracy and simplicity and is able to account for the basic qualitative features of a pandemic, in particular in its early stages. Optimal decisions are taken in terms of optimal timing of lock-down/release, and optimal percentage of quarantined/released susceptibles. The technique is then applied in a realistic simulation scenario based on the data of COVID-19 evolution in Italy.

The analysis carried out in this paper suggests that when dealing with a pandemic such as the one caused by the novel coronavirus it may be necessary to take immediate measurements by locking-down a vast majority of the population at the very beginning of the epidemic spread. In fact, the analysis carried out in Section 5 suggests to take measurements in the first 20 days, which may further reduce in the case that there is a delay between measurements and their effects, as it has been evinced in the COVID-19 spread in Italy. Moreover, the results suggest that a strong lock-down, removing about 50-60 millions of individuals from the susceptible population, is actually required to avoid a large number of deaths.

The results also show that the number of released individuals can deeply change on the basis of the relative infectivity at the realizing time: such a value can be different from the initial one because of the increasing awareness (during the lock-down) on the transmission risk and the actualization of some precautionary measures (mask wearing, social distancing, hygienic measures, etc.). In future research efforts we aim at extending the framework by considering the effect of stochastic fluctuations in the pandemic evolution, which will be modeled by means of properly defined continuoustime Markov chains (CTMCs) (Van Kampen, 1992) or stochastic differential equations (SDEs) (Langevin, 1908) . Further, note that the lock-down and reopening policies proposed in this work are open loop. Another possible extension of the results given in this paper is to envision an MPC-like strategy to transform these policies into feedback control laws.

In fact, if the precautionary measures used to reduce the infection rate are not sufficient to limit the post lock-down spread of the epidemics, one may envision an iterative procedure to design periodic lock-downs and reopenings based on the machinery proposed in this work. Namely, the OLS approach outlined at the beginning of Section 5 can be used to estimate the current model parameters using the new available data over a moving time window. This updated model can be then employed as outlined in Section 5 to periodically design lockdown and reopening strategies with the aim of balancing health and economic costs.

The mathematical theory of infectious diseases and its applications

The challenges of modeling and forecasting the spread of COVID-19

Model predictive control to mitigate the COVID-19 outbreak in a multi-region scenario

Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

Hybrid dynamical systems

COVID-19: from model prediction to model predictive control

Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic model and of the SIR model with equal death and birth rates

Demographic balance for the year 2018 and resident population in Italy

A contribution to the mathematical theory of epidemics

Sur la théorie du mouvement brownien

Early transmission dynamics in Wuhan, China, of novel coronavirus-infected pneumonia

An optimal predictive control strategy for COVID-19 (SARS-CoV-2) social distancing policies in Brazil

The effect of control strategies to reduce social mixing on outcomes of the COVID-19 epidemic in wuhan, China: a modelling study

Why is it difficult to accurately predict the COVID-19 epidemic?

COVID-19 data in Italy

Principles of mathematical analysis

The epidemiological characteristics of an outbreak of 2019 novel coronavirus diseases (COVID-19)-China

An updated estimation of the risk of transmission of the novel coronavirus (2019-nCov)

Estimation of the transmission risk of the 2019-nCoV and its implication for public health interventions

Stochastic processes in physics and chemistry

Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study

Estimation of the reproductive number of Novel Coronavirus (COVID-19) and the probable outbreak size on the Diamond Princess cruise ship: A data-driven analysis

Preliminary estimation of the basic reproduction number of novel coronavirus

A data-driven analysis in the early phase of the outbreak

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.