key: cord-0978160-lvp1jcil
authors: Coraggio, Marco; Xie, Shihao; Lellis, Francesco De; Russo, Giovanni; Bernardo, Mario di
title: Intermittent non-pharmaceutical strategies to mitigate the COVID-19 epidemic in a network model of Italy via constrained optimization
date: 2021-03-25
journal: 60th IEEE Conference on Decision and Control, CDC 2021
DOI: 10.1109/cdc45484.2021.9683420
sha: de6d18077bd47ea2738ee1b3e172bc87988d1350
doc_id: 978160
cord_uid: lvp1jcil

This paper is concerned with the design of intermittent non-pharmaceutical strategies to mitigate the spread of the COVID-19 epidemic exploiting network epidemiological models. Specifically, by studying a variational equation for the dynamics of the infected in a network model of the epidemic spread, we derive, using contractivity arguments, a condition that can be used to guarantee that, in epidemiological terms, the effective reproduction number is less than unity. This condition has three advantages: (i) it is easily computable; (ii) it is directly related to the model parameters; (iii) it can be used to enforce a scalability condition that prohibits the amplification of disturbances within the network system. We then include satisfaction of such a condition as a constraint in a Model Predictive Control problem so as to mitigate (or suppress) the spread of the epidemic while minimizing the economic impact of the interventions. A data-driven model of Italy as a network of three macro-regions (North, Center, and South), whose parameters are identified from real data, is used to illustrate and evaluate the effectiveness of the proposed control strategy.

Because of recent events, the design of control strategies aimed at mitigating the spread of epidemic processes has suddenly become a major pressing research problem with clear societal and economic impacts [1] , [2] . Among these strategies, nonpharmaceutical interventions (NPIs) are essential to mitigate or suppress the epidemic (see [3] , [4] for further details on the definition and difference between mitigation and suppression). In this context, two key challenges are: (i) devising proxies to detect the onset of the epidemic spreading; (ii) synthesizing control strategies that balance the need of containing the epidemic with the need of reducing its economic impact.

Related work: A fast intermittent lockdown strategy was proposed in [1] with the goal of alleviating the effects of measurement uncertainty. Following a similar argument, a cyclic periodic lockdown strategy was studied in [5] , whereas a stepping down policy was presented in [6] . The use of feedback to dynamically adapt the intervention strategies was recently suggested in several papers. For example, in [7] , several feedback intervention strategies are presented, whereas in [8] the benefits were shown of adopting feedback intermittent regional control strategies triggered by occupancy levels of the intensive care units (ICUs) in each region. Furthermore, several epidemic control strategies were recently proposed, based on model predictive control (MPC). These include [9] , where MPC is used to decide whether to enforce social distancing in a region/country of interest, using a SIRASDC model ( : susceptible, : infected, : recovered, : asymptomatic, : symptomatic, : dead, : with control) and including hard constraints on ICUs occupancy and minimal dwell time of the control. MPC was also implemented on a number of different models, e.g., a SIDARTHE model ( : susceptible, : infected, : diagnosed, : ailing, : recognized, : threatened, : healed, : extinct) in [10] , a SEAIR model ( : susceptible, : exposed, : asymptomatic infected, : confirmed infections, : recovered) in [11] , and a multi-region SIRQTHE network model ( : susceptible, : infected, : recovered, : quarantined, : threatened, : healed, : extinct) in [12] .

Contribution: We consider network epidemic models, and give a sufficient condition, directly dependent on the model parameters, ensuring that the epidemiological effective reproduction number R is less than one and that the infected population decreases to zero. Then, we consider a network model of Italy parametrized on real data, and obtain a set of intermittent non-pharmaceutical intervention strategies at the regional level by formulating a constrained optimization problem that is effectively solved via MPC. The formulation (i) embeds our novel stability condition as a constraint, (ii) is based on the minimization of an economic cost induced by both social distancing and testing and (iii) considers a lower bound on the control dwell time. To the best of our knowledge, we are the first to implement the MPC on a regional model including all these features.

We denote by I the identity matrix of appropriate dimensions, by 0 and 1 column vectors of appropriate dimensions consisting of all zeros and ones, respectively; · denotes the vector ℓ -norm and its induced matrix norm. R 0 is the positive orthant in R (including the axes); N ≥0 and N >0 are the sets of semipositive and positive natural numbers, respectively. Letting A be a matrix, A 0 indicates that all its elements ≥ 0 (analogously for , , ≺).

Lemma 1 ( [13] , p. 368). Let B 1 , B 2 ∈ R × , with ∈ N >0 . If B 1 0, B 2 0, and B 1 B 2 , then B 1 ≤ B 2 .

To model the dynamics of the epidemic, we adapted and updated a network model that was originally presented and validated in [8] . We re-parameterized and updated the model with real data related to the spread of the epidemic in Italy collected by the National Civil Protection Authority [14] , between February 24, 2020 and February 25, 2021. The model comprises = 3 macro-regions corresponding to the North, the Center, and the South of Italy. 1 The dynamics of the disease in each region ∈ {1, . . . , } is captured by a SIQHDR model, describing the discrete-time evolution of the number of susceptible ( ), (undetected) infected ( ), quarantined ( ), hospitalized ( ), deceased ( ), and recovered ( ) as a function of time ∈ N ≥0 (in days). We also define S

Moreover, for each pair ( , ) of macro-regions, ( ) ∈ [0, 1] describes the fraction of people in region that commute and interact with people in region , returning to their region of origin at the end of each day (see also [15] ). The resulting network model is

and the initial conditions are given at time 0 ∈ N ≥0 as x( 0 ) = x 0 ∈ R 6 0 . Following [8] , in (1) the infection rate and the recovery rate are taken to be equal for all the macro-regions, as they have not been shown to strongly depend on regional factors. The other parameters vary among different macro-regions to capture the heterogeneity of containment strategies, different regional health care systems, etc. In particular, note that

as in (1b) the number of people leaving compartment can never be greater than . Moreover, as done in [8] , we also estimate the number of patients requiring treatment in intensive care as 10% of those who are hospitalized. Thus, letting be the maximum number of beds in the ICUs of region [16] , [17] , following [8] we express the mortality rate as ( ( )) = 0 + b min 0.1 ( ) , 1 . Finally, p ( ) is the population in region that is free-to-move at time , whereas is the total population in region . All parameters' values are identified from data reported in [14] and are contained in Table  I .

Social distancing: The parameter ( ) ∈ [0, 1] models the effect of social distancing in region at time . Namely, ( ) = 1 corresponds to the absence of social distancing measures, whereas ( ) = 0 means that contact between people is forbidden. Intermediate values indicate different levels of social distancing.

Travel restrictions: In what follows we assume [18] where ( ) ∈ [0, 1] is a parameter modelling the strength of travel restrictions to/from region on day , and 0 is the influx from region to region in the absence of travel restrictions. In particular, ( ) = 1 means that travel to/from region is possible, whereas ( ) = 0 implies it is forbidden.

Testing: The effect of testing-and-tracing is captured by the parameter ( ) in (1), which models the rate at which undetected infected are moved to the quarantined compartment at time . To describe the ability of a region to modulate the amount of testing carried out within its territory, we set ( ) = 0 + ( )˜, where ( ) ∈ [0, 1] is a parameter denoting the testing capacity of region on day . In particular, testing is set to the minimal level, 0 ∈ R >0 , when ( ) = 0, or to its maximal value, 0 +˜, when ( ) = 1.

The (1)) over the fraction of population tested daily between April 22, 2020 and February 25, 2021, taken from [14] .

A crucial quantity to assess the rate of propagation of an epidemic disease is the reproduction number, defined in epidemiology as "the number of secondary infections generated by one infected person" [19, Ch. 27] . Typically, the basic reproduction number R 0 is computed when the entire population is assumed to be susceptible [20] , while the effective reproduction number R is estimated from data [19] , [21] . Specifically, the effective reproduction number of region , denoted by R , , can be estimated via the method described in [22, §8] , that is, by computing

Note that while an analytical expression of R 0 for network epidemic models can be found in the existing literature, e.g., [15] , [23] , at the moment an algebraic expression for R does not exist in a closed form for network models. Later in Section III, we derive a sufficient condition to ensure that R (for the network) is below the critical unitary threshold.

We let ( )

, which will be used later in Sections IV and V as the vector of control inputs to mitigate the epidemic.

We focus on the dynamics of the infected recasting (1b) in a more compact form as

where we stressed the dependency of on the state vector x and the control input u, and where

In what follows, sayx( )

T the state evolution of (1) starting from some initial conditionx( 0 ) =x 0 ∈ R 6 0 , 0 ∈ N ≥0 , when u( ) is equal to someū( ). Additionally, following the epidemiological definition of R we obtain an expression that yields the variations of I( ), say I( ), aroundĪ( ), due to small perturbations in the infected only. Formally, this can be done by considering in (5) all the states (x( )) except I( ) as inputs to the dynamics of the infected compartment. After some algebraic manipulations, this dynamics can be given as 2 I( + 1) = (I + (x( ),ū( )) + (x( ),ū( ))) I( ),

where

Next, we use (5) and (7) to give a simple yet effective condition to (i) enforce suppression of the epidemic and (ii) to ensure that the number of secondary infections generated by perturbing the number of infected along a trajectory decays exponentially to zero. Proposition 1. Assume that there exists some matrix norm, say · , and some constant < 1 such that

Then, (i) Ī ( ) ≤ − 0 Ī ( 0 ) for all ≥ 0 , and (ii) I( ) ≤ − 0 I( 0 ) for all ≥ 0 .

Proof. For the sake of brevity, we denote (x( ),ū( )) by ( ) and (x( ),ū( )) by ( ). Part (i) follows directly from the fact that, from (5) 

. Analogously, to prove part (ii), it suffices to exploit (7) and show that (9) implies

To show this, we let B 1 ( ) I + ( ) + ( ) and B 2 ( ) I + ( ), and note the following facts: (a) recalling (3) and (6), we have B 2 ( ) 0, ∀ ∈ N ≥0 ; (b) as (x, u) 0, ∀x ∈ R 6 0 , ∀u ∈ [0, 1] 3 , it holds that B 1 ( ) B 2 ( ), ∀ ∈ N ≥0 ; (c) we can rewrite the term in (8), using (2), as

and, recalling (3), (6), (8), we get that B 1 ( ) 0, ∀ ∈ N ≥0 . Finally, exploiting the facts (a), (b), (c) we can apply Lemma 1 and, considering (9), we have B 1 ( ) ≤ B 2 ( ) ≤ , ∀ ≥ 0 , that is (10), which proves (ii).

Next, we give an interpretation of Proposition 1 in terms of R . Recall that R is defined as the average number ≥ 0 of secondary infections generated by = 1 infected people. The well-known condition R < 1 means that one infected person will infect on average less than one person, and the epidemic will disappear, with the total number of infected (that is I( ) 1 = =1 | ( )|) decreasing at all time. Indeed, this corresponds to part (i) of the thesis of Proposition 1, when considered with = 1. Moreover, noting that in our framework corresponds to I( ) 1 , it is obvious that each person infects on average less than one person if and only if I( ) 1 decreases monotonically over time, which corresponds to part (ii) of Proposition 1.

Additionally, as recently noted in [24, Remark 6] , the use of the ∞-norm prohibits the amplification of bounded deterministic disturbance on the infected dynamics within the network system. Its application in condition (9) yields max ∈ {1,..., } | ( + 1)| < max ∈ {1,..., } | ( )|, which means that the maximum variation of infected across all regions decreases over time. We argue that this condition might be effectively used to assess quantitatively the trend of an epidemic spread.

As proven in the Appendix, Proposition 1 can also be related to a condition independently obtained in [25] when specialized to a SIS network model (see the Appendix for the formal statement and further details).

We now formulate our control problems, aimed at either mitigating or suppressing the epidemic, using different combinations of control actions, while minimizing the economic impact of the measures.

Suppression of the epidemic: First, we consider the problem of finding the cheapest strategies in terms of social distancing ( ) and travel restrictions ( ) that enforce suppression of the epidemic, assuming testing ( ) is not increased beyond its nominal value. To this aim, we define-omitting time dependence for the sake of brevity,

Clearly, if ≤ < 1, ∀ , ∀ ≥ 0 , then (9) in Proposition 1 holds with the infinity norm, and the infected population decreases to zero. Thus, the control goal can be formulated as the constrained optimization problem

where [1, ] is the economic cost function over the time window [1, ] , and will be presented in Section IV-A, U [1, −1] (u(1), u(2), . . . , u( − 1)) is the set of all decision variables, S , S are sets of possible values of and , which will be selected in Section V, and u is the minimum dwell time on u. In loose terms, (12c) means that the value of u cannot be changed before at least * ∈ N >0 days have passed since the last change (similarly to what done in [9] , [12] ); a formal implementation of (12c) is detailed in Section V.

Mitigation of the epidemic: To better preserve the economy, and as currently done in many countries affected by the epidemic, we also study a different approach, where some of the restrictions can be lifted when certain conditions are met.

To this aim, we obtain a relaxation of condition (9) by considering the condition

where ( ) is a Boolean condition associated to the "criticality" of the epidemic in region , that we choose as

linking it to the level of saturation of the ICUs in each region and the value of R , in each region, estimated using (4), with ∈ (0, 1], R ∈ R >0 . The resulting control problem can be then formulated as (12) with (12b) substituted by (13) . The effect of testing: Finally, we explore if and how increasing the amount of testing beyond the nominal value in each region can further reduce costs by formulating the problem as (12) where the constraint on ( ) is chosen as ( ) ∈ S , ∀ , ∀ , with S being a set of possible values, selected later in Section V.

As a cost function for the optimization problems above, we consider the total economic cost of the measures undertaken by each region, computed accounting for the cost of people not being able to work and the cost of testing.

Let active, ( ) ( ) + ( ) + ( ); m = e83.992 be the mean daily GDP pro-capita in Italy [26] , w = 0.617 be the fraction of people that are not able to work when social distancing is in place [27] , and be the daily cost of carrying out maximal testing, normalized over the population. Then, the economic cost at time in region is given by ( ) = ,1 ( ) + ,2 ( ) + ,3 ( ) + ,4 ( ) + ,5 ( ), where the addends are defined as follows: 1. ,1 ( ) is the loss due to people that live and work in region , but cannot work because of social distancing: 5. ,5 ( ) is related to the cost of carrying out increased testing in region and is computed as ,5 ( ) ( ). The expressions in items 1. and 2. were obtained by noticing that in (1)-omitting subscripts for notational easiness-the amount of new infected at each time step is proportional to = ( √ ) ( √ ). Therefore, the portion of people that is not restrained is proportional to √ , and consequently the fraction of population that is restrained is proportional to (1 − √ ).

In item 5., =˜t est = 0.3583, where test = e59 is the economic cost of a single test [28] , = 1.74 is the average number of tests carried out per tested person [14] , and˜= 3.49 × 10 −3 is the fraction of people tested daily when testing is maximal. Letting the regional terminal cost be¯( ) = 0, and assuming the presence of a discount factor 0 < ≤ 1, the total cost over the time window [ 0 , ] can then be given as

).

We solved the control problems described in the previous section by using a Model Predictive Control approach with prediction horizon p [29] .

We chose the sets S , S , S of possible values of the control inputs as follows. For ( ), we approximated the values of this parameter estimated in [8] for each of the Italian regions when different levels of social distancing were enforced, obtaining S = {0.3, 0.4, 0.5, 0.6, 0.7}. For ( ), we replicated the decrease of 70% in mobility fluxes among regions as registered in [30] during the first months of the Italian national lockdown (February-May 2020), obtaining S = { √ 0.3, 1}. As for ( ), we chose three values to represent low, medium, and intense testing levels, i.e., S = {0, 0.5, 1}.

To implement the dwell time constraint (12c), we prescribe that the control input u on day can be changed only if ( ) = * , where ( ) is the number of days (up to * ), counting back from (excluding ), such that the control variables stayed equal to the value they had in − 1. To compute ( ), we initialize (1) = 0, and update according to the law

Moreover, to simplify computations, when ( ) = * , if the MPC finds that u( ) can be applied continuously for days, and that u( ) = u( − 1), we apply u( ) for min{ , MPC } days (rather than just one day), for some MPC ∈ N >0 .

The optimization problems are solved using a genetic algorithm implemented in MATLAB 2020 [31] , . Plots are shaded in red when condition (9) is applied. In the middle panels, the control variables (blue), (orange), and (yellow) are depicted, with the plots being shaded in grey when ( ) < * (and the control values cannot be changed). In the bottom panels, the current estimate of R , estimated using (4) is plotted.

A. Numerical results 1) Epidemic suppression: We can observe the solution to (12) (with constraint (12b)) in Figure 2 . As expected, the epidemic is effectively eradicated by the control strategy, as the number of infected people asymptotically goes to zero. However, this strategy is particularly expensive with a total cost of e421.989 · 10 9 (computed with = 1).

2) Mitigation: The results of the control strategy when applied to solve (12) with constraint (13) is reported in Figure 3 . In this case, the total cost is e337.172 · 10 9 , significantly lower than the previous scenario. This is because in each region NPIs interventions are taken only when some safety thresholds are met and only in that specific region. At the same time, the number of people in ICUs remains within capacity, and all regional effective reproduction numbers R , remain within acceptable values.

3) Effect of additional testing: Finally, we report the solution to the problem when ( ) ∈ S in Figure 4 . The total cost is e262.379 · 10 9 , even lower than the previous case. Clearly, the use of additional testing allows to further reduce costs by relying less on social distancing to mitigate the epidemic. Indeed, in Figure 4 , it is possible to note that travel restrictions are never used in this case, and social distancing is significantly less severe when compared to that in Figures 2 and 3 .

We considered the problem of designing localized intermittent non-pharmaceutical intervention strategies to mitigate the COVID-19 epidemic in a regional model of the Italian epidemic. We first gave a sufficient condition to determine the stability of the infected population, in a network epidemic model. We then exploited this condition, or a relaxed version of it, as (13) . Top panels are shaded in red when (x( ), u( )) ≤ (see (11) ), while green is used when (x( ), u( )) > . In the bottom panels, the threshold value R for R , is plotted in a dash-dotted line. a constraint to formulate the control problem as an optimization problem in which the levels of social distancing, travel restrictions and additional testing must be selected to avoid saturation of the regional health systems and containment (or suppression) of the epidemic spread, while minimizing the overall cost to the economy. Our results show the effectiveness of the proposed approach. In the future, we will study robustness of our control strategy to model uncertainties, and measurement error.

The code and data to run the model identification ( § II), the linear regression ( § II), and the simulations ( § V-A) are available at https://github.com/diBernardoGroup/MPC-to-mitigate-COVID-19. 

Post-lockdown abatement of COVID-19 by fast periodic switching

Projecting the transmission dynamics of SARS-CoV-2 through the postpandemic period

Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand

Inferring the effectiveness of government interventions against COVID-19

Cyclic exit strategies to suppress COVID-19 and allow economic activity

Modeling the effects of intervention strategies on COVID-19 transmission dynamics

How control theory can help us control COVID-19

A network model of Italy shows that intermittent regional strategies can alleviate the COVID-19 epidemic

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

Robust and optimal predictive control of the COVID-19 outbreak

Modeling, state estimation, and optimal control for the US COVID-19 outbreak

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

Matrix Analysis

Dipartimento della protezione civile-presidenza del consiglio dei ministri official COVID-19 data repository

SIR-network model and its application to dengue fever

Ministero della Salute. (2020) Posti letto per regione e disciplina

Linee di indirizzo organizzative per il potenziamento della rete ospedaliera per emergenza COVID-19

Ricoverati e posti letto in area non critica e terapia intensiva

Modern epidemiology

Perspectives on the basic reproductive ratio

The Effective Reproduction Number as a Prelude to Statistical Estimation of Time-Dependent Epidemic Trends

Demystifying the math of the coronavirus

On the definition and the computation of the basic reproduction ratio 0 in models for infectious-diseases in heterogeneous populations

Scalability in nonlinear network systems affected by delays and disturbances

Optimal lockdown for pandemic control

World Economic Outlook Database

Who and how many can work from home? evidence from task descriptions

Analisi dei modelli organizzativi di risposta al Covid-19

Model predictive control: past, present and future

Google mobility data

Global Optimization Toolbox

For the sake of completeness, we explore here the relationship between the stability condition in Proposition 1 and a recent condition derived independently in [25] , when both are applied to a discrete-time SIS model (adapted from the continuous-time model adopted in [25] ).Namely, the dynamics of the infected we consider is given bywhere ∈ [0, 1] is the portion of infected in region , with ∈ {1, . . . , }, I [ 1 · · · ] T , 0 ≤ ≤ 1 is the infection rate, 0 ≤ ≤ 1 is the recovery rate, A ∈ R × captures the rate at which infection flows between regions (also accounting for social distancing), with A 0, and diag(v) is a diagonal matrix having on its main diagonal the vector v.Next, we derive and compare two different upper bounds on the maximum eigenvalue of A, that guarantee the stability of the infected dynamics in the origin. The first bound, say (B1), is derived following the arguments in Proposition 1, whereas the second, say (B2), is computed following the steps in [25] .Bound derived from Proposition 1: By using similar arguments as those used in Proposition 1, it can be shown that, if there exists some matrix norm and some < 1 such thatthen I( ) ≤ − 0 I( 0 ) for all ≥ 0 . Letting m ( ) min ( ), and applying the triangle inequality and Lemma 1, we find that (17) is satisfied ifwhere we used the fact that the variables ( ) are bounded between 0 and 1. As A ≥ max (A), then there exists ∈ (0, 1] such that A = max (A). Thus, we rewrite (18) and obtain the first bound asBound derived from [25] : Following the arguments in [25] , we let P (1 − )I + A, noting that P 0. Applying the Perron-Frobenius theorem [13, Theorem 8.3 .1] to P T yields the existence of an eigenvalue max (P) ≥ 0, with associated left eigenvector v max 0. Moreover, max (P) is such that max (P) ≥ | (P)|, for all other eigenvalues of P. Now, if(which constitutes the second bound), we consequently get max (P) ≤ < 1, and it is then immediate to show that this implies that 0 is a stable equilibrium for the dynamics I( ) given in (16) . Indeed, we have v T max I(which means stability of I( ) in 0 recalling that I( ) 0 and v max 0. Note that (B2) is equivalent to the stability condition in Proposition 3.1 in [25] when applied to a discrete-time SIS network model.Comparison: Finally, we relate the two stability conditions (B1) and (B2). Clearly, when m ( ) ∈ (1 − , 1), 1 > 2 and our condition (B1) is less conservative than (B2), while it is more conservative ( 1 < 2 ) when m ( ) ∈ [0, 1 − ).