key: cord-0471731-o2plh7bl
authors: Armaou, Antonis; Katch, Bryce; Russo, Lucia; Siettos, Constantinos
title: Designing Social Distancing Policies for the COVID-19 Pandemic: A probabilistic model predictive control approach
date: 2022-02-02
journal: nan
DOI: nan
sha: 89397d4f064ef3e579d51a75682d2f659ece236d
doc_id: 471731
cord_uid: o2plh7bl

The effective control of the COVID-19 pandemic is one the most challenging issues of nowadays. The design of optimal control policies is perplexed from a variety of social, political, economical and epidemiological factors. Here, based on epidemiological data reported in recent studies for the Italian region of Lombardy, which experienced one of the largest and most devastating outbreaks in Europe during the first wave of the pandemic, we address a probabilistic model predictive control (PMPC) approach for the modelling and the systematic study of what if scenarios of the social distancing in a retrospective analysis for the first wave of the pandemic in Lombardy. The performance of the proposed PMPC scheme was assessed based on simulations of a compartmental model that was developed to quantify the uncertainty in the level of the asymptomatic cases in the population, and the synergistic effect of social distancing in various activities, and public awareness campaign prompting people to adopt cautious behaviors to reduce the risk of disease transmission. The PMPC scheme takes into account the social mixing effect, i.e. the effect of the various activities in the potential transmission of the disease. The proposed approach demonstrates the utility of a PMPC approach in addressing COVID-19 transmission and implementing public relaxation policies.

The structure of the paper is as follows. In section 2.1, we present the compartmental dynamical model and we justify the selection of the values of the model parameters based on publicly available epidemiological data. Based on the model, in section 2.2, we derive estimations of the basic and effective reproduction numbers. In section 2.3 we perform an analysis of the Google mobility reports for Lombardy, thus providing a model for the effect of the social distancing measures in all activities to the residential activity which is also an important factor in the transmission of infectious diseases. In section 2.4, we present the proposed PMPC scheme constrained to the impact of the social mixing and activity patterns in the potential transmission. In section 3, we report the simulation results, and in section 4 we conclude with the discussion of the findings. Finally, in the Supplementary Information section, we present the sources of the data used in the present study, competing formulations to the designed PMPC one and an analysis of the effect of PMPC parameters on the resulting policies.

As in our probabilistic model predictive control scheme, we consider as constraint a maximum number of hospitalized persons due to COVID19, we extended the model proposed in Russo et al. [17] to include hospitalized cases and confirmed infected cases with mild symptoms that are isolated at home. Thus, it includes two categories of infected cases, namely the confirmed/reported and the unreported (unknown) asymptomatic cases in the total population [17] . Based on observations and studies, our modelling hypothesis is that the confirmed cases of infected are only a (small) subset of the actual number of infected cases in the total population [20, 6] . Regarding the confirmed cases a study conducted by the Chinese CDC which was based on a total of 72,314 cases in China, about 80.9% of the cases were mild and could recover at home, 13.8% severe and 4.7% critical [21] .

On the basis of the above findings, in our modelling approach, the asymptomatic cases recover from the disease relatively soon and without medical care, while the confirmed cases include all the above types, but on average their recovery lasts longer than the asymptomatic cases, they may also be hospitalized or sadly die from the disease [17] .

Considering a well-mixed population of size N , the state of the system at time t is described by the following compartments (see also Figure 1 for a schematic): S(t) representing the number of susceptible persons, E(t) the number of exposed, I A (t) the number of asymptomatic infected persons in the total population that recover relatively soon without any other complications; R(t) is the number of recovered asymptomatic persons. I S (t) denotes the number of infected cases who may develop more severe symptoms and a part of them is hospitalized (in our model denoted by H(t)), a part of them dies (denoted by D(t)) and a part of them recovers (denoted by R H (t)). The other part of confirmed infected cases from the I S (t) compartment (in our model represented by the variable Q(t)) is assumed to experience mild symptoms and is quarantined at home. For this compartment there are no other available data except from its specific number. The confirmed recovered (in our model denoted by R H (t)) are the recovered cases dismissed from hospitals. Furthermore, there are no reported data for the number of deaths from the infected people that are isolated at home. As said this compartment represents the part of cases that experience mild symptoms and thus they aren't expected to die from the disease.

Here, we remark that a wide testing policy may also result in the identification of asymptomatic cases belonging to the compartment I A that would then be assigned to compartment I S . However, as a generally reported rule in Italy, tests were conducted only for those who asked for medical care with symptoms like fever and coughing. Thus, people who did not seek medical attention were tested on rare occasions [22] . Hence, for any practical means the compartment I S reflects the cases with mild to more severe symptoms.

In our model: a) the number of mild to severe symptomatic cases as represented by the compartment I S does not coincide with the reported/confirmed number of infected cases as: (i) not all symptomatic cases are reported as confirmed cases, (ii) the confirmed reported total positive cases is the sum of the hospitalized and guaranteed cases, b) the number of confirmed recovered cases does not coincide with the total number of recovered cases dismissed from the hospital R H ; the reported cases include also a part of the confirmed cases that were isolated at home and after a second test were identified as recovered [23] ; thus we expect that the number of R H to be less that the reported one, c) there are no data for the number of recovered from the compartment Q; since these cases experience mild symptoms, we assume there are no deaths from this compartment. Hence, in our model only the sum of the H and Q compartments reflect the number of actual confirmed cases.

For our analysis, and for such a short period, we assume that the total number of the population remains constant. Based on demographic data, the total population of Lombardy is N = 10m. The rate at which a susceptible (S) becomes exposed (E) to the virus is proportional to the density of asymptomatic infectious persons I A (t) in the total population and the symptomatic ones I S (t) the period preceding their hospitalization or home quarantine. The Figure 1 : A schematic of the proposed compartmental model with hidden compartments. Dashed lines depict the compartments for which there are no available data. The released data report the total number of confirmed cases I S , the total number of hospitalized cases H, the total number of recovered cases dismissed from the hospital R H , the total number of confirmed infected isolated at home Q and the total number of deaths D.

proportionality constant is the "effective" disease transmission rate, say β =cp, wherec is the average number of contacts per day and p is the probability of infection upon a contact between a susceptible and an infected.

The assumption of the model is that only a fraction, ǫ, of the actual number of exposed cases E(t) will experience mild to more severe symptoms and ask for medical treatment; this fraction is represented by the compartment I S (t). Thus, our discrete mean field compartmental SEASQHRD model reads:

The model is defined in discrete time points t = 1, 2, . . ., with the corresponding initial condition at the very start of the pandemic (t 0 =0=DAY-ZERO) being

The values of the epidemiological parameters were fixed in the proposed model at values reported by clinical studies as follows:

• β(d −1 ) is the "effective" transmission rate of the disease; recently published studies have reported that this rate is for any practical means the same for both asymptomatic and symptomatic cases (see e.g. [24] ). The transmission rate β cannot be obtained by clinical studies, but only by mathematical models. Here we have set β = 0.68 based on a previous work for the region of Lombardy [17] .

• σ(d −1 ) is the average per-day "effective" rate at which an exposed person becomes infectious. In many studies, this is set equal to the inverse of the mean incubation period (i.e. the time from exposure to the development of symptoms). A study in China [25] suggests that it may range from 2 to 14 days, with a median of 5.2 days. Another study using data from 1,099 patients with laboratory-confirmed 2019-nCoV ARD from 552 hospitals in 31 provinces/provincial municipalities of China suggested that the median incubation period is 4 days (interquartile range, 2 to 7). However, the incubation period does not generally coincide with the time from exposure to the time that someone starts to be infectious. Regarding COVID-19, it has been suggested that an exposed person can be infectious well before the development of symptoms [26] . Indeed, it has been suggested that for a mean incubation period of 5.2 days infectiousness starts from 2.3 days (95% CI, 0.8-3.0 days) before symptom onset [26] . In our model, as explained above, 1 /σ represents the period from exposure to the onset of the contagious period and according to the above has been set as σ = 1 /3.

• h(d −1 ) is the average per-day rate at which an infected person is hospitalized. The inverse of the rate of hospitalization h, i.e. 1 /h is the mean time from the onset of symptoms to hospitalization for those infected cases that need hospitalization. For Lombardy, the median value of this parameter has been reported to be 7 days [27] . Thus, considering also the fact that the infectiousness periods start 2-3 days before the onset of symptoms during the incubation period as explained above, we have set 1 /h = 4 + 3 and thus h = 1 /7. This coincides with the corresponding number reported by the Istituto Superiore di Sanit, Italy [27] .

is the average per-day "effective" recovery rate within the group of asymptomatic cases in the total population. Here, it was set equal to ≅ 1 /6.6 based on studies reporting that infectiousness declines quickly within ∼7 days [26] ; this coincides with the mean value of the serial interval reported for Lombardy (6.6 days) [28] but also that for China (mean serial interval of 6.3 days (95% CI 5·2-7·6) ) [29] .

• δ H (d −1 ) is the average per-day "effective" recovery rate within the subset of hospitalized cases (H) that finally recover. As reported recently [16] , the overall median length of stay (LoS) in hospitals, which for any practical means it coincides with the recovery period, decreased steadily from 21.4 (20. • δ Q (d −1 ) is the average per-day "effective" recovery rate within the subset of quarantined cases (Q). In a study that is based on 55,924 laboratory-confirmed cases, the WHO-China Joint Mission has reported a median time of 2 weeks from onset to clinical recovery for mild cases [30] . Hence, based on these reports, we have set δ Q = 1 /14 for the quarantined cases.

• γ(d −1 ) is the average per-day "effective" mortality rate within the subset of H(t) hospitalized cases that finally die. The median time from hospitalization to death for Lombardy during the early phase of the first wave of the pandemic (until April) has been reported to be seven days through the first wave, with an interquartile range (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) [31, 16] . We should also note that the time from hospitalization to death subsequently increased as medical personnel gained experience. Thus based on the above the time from hospitalization to death was set as γ = 1 /5 for the period until March 20, γ = 1 /7 for the period March 21-April 10, and δ H = 1 /10 for the period April 11-May 4.

• f 1 is the ratio of all medium to severe infected cases (represented by the I S compartment) that gets hospitalized. For the early phase of the pandemic (until March 8th), a 50% over a total 5,830 confirmed cases reported in Lombardy were hospitalized [28] , this percentage is also confirmed by other studies (see e.g. [32] ). Regarding the next period (March 20 th to May 4), that was characterized by a saturation of the hospitals beds it is reasonable to expect that this ratio would fall significantly in addition due to the fear of a high-risk nosocomial infection. Indeed, for the period March 20 and on this ratio was about 20% as reported also by other studies (see e.g. [32] ). Based on these studies and through the analysis of the provided data, this ratio was set to f 1 = 0.65 for the period until March 8, f 1 = 0.60 for the period March 9 to March 20, f 1 = 0.5 for the period March 21 to April 10, and f 1 = 0.2 for the period April 11 to May 4.

• f 2 is the case fatality ratio of all hospitalized cases. Based on a study considering 38,715 hospitalized COVID-19 patients in Lombardy between 21 February and 21 April 2020, this ratio was relatively high [33, 16] ) mainly due to the saturation of ICU beds. In particular as reported in a recent published study [16] for the region of Lombardy is estimated to have steadily decreased from 34.6% (32.5-36.6%) in February to 29.8% • ǫ(d −1 ) is the fraction of the exposed cases that enter compartment I S , i.e. those that will experience mild to severe symptoms and ask for medical care. For the case of Italy and in particular that of Lombardy during the first wave of the pandemic; the ratio of the confirmed to unreported cases was estimated to be around 0.07 (interquartile range: 0.023 to 0.1), [17] . In our model, we have set ǫ = 0.12 for the period until March 8, ǫ = 0.10 for the period March 9 to March 20, ǫ = 0.09 for the period March 21 to April 10, and ǫ = 0.05 for the period April 11 to May 4.

Regarding DAY-ZERO in Lombardy, i.e. the initial date of the introduction of the disease in the region, that was set to January 15, as estimated based on a methodological framework based on genetic algorithms together with the level of asymptomatic cases in the total population [17] . Another study [34] based on genomic and phylogenetic data analysis, reports the same time period, i.e. that between the second half of January and early February, 2020, as the time when the novel coronavirus SARS-CoV-2 entered northern Italy.

Note, that there are three infected compartments, namely E, I A , I S that determine the fate of the pandemic. Thus, considering the corresponding equations given by Eq.(2.2),(2.3), (2.4) , and that at the very first days of the epidemic S ≈ N , the Jacobian of the system as evaluated at the disease-free state is:

The eigenvalues (that is the roots of the characteristic polynomial of the Jacobian matrix) dictate if the disease-free equilibrium is stable or not, that is if an emerging infectious disease can spread in the population. In particular, the disease-free state is stable, meaning that an infectious disease will not result in an outbreak, if and only if all the norms of the eigenvalues of the Jacobian J of the discrete time system are bounded by one. Jury's stability criterion [35] (the analogue of Routh-Hurwitz criterion for discrete-time systems) can be used to determine the stability of the linearized discrete time system by analysis of the coefficients of its characteristic polynomial. The characteristic polynomial of the Jacobian matrix is:

The necessary conditions for stability read:

Thus, the sufficient conditions for stability are given by the following two inequalities:

It can be shown, that the second necessary condition of (2.15) and the first sufficient condition of (2.16) are always satisfied for the range of values of the epidemiological parameters considered here. The first inequality (2.14) results in the necessary condition:

It can be also shown, that for the range of the parameters considered here, the second sufficient condition of (2.17) is satisfied if the necessary condition of (2.19) is satisfied. Thus, the necessary condition of (2.19) is also a sufficient condition for stability. Hence, the disease-free state is stable, if and only if, condition of (2.19) is satisfied.

Note that in this necessary and sufficient condition of (2.19), the first term reflects the contribution of compartment I A and the second term the contribution of compartment I S in the spread of the disease. Thus, the above expression reflects the basic reproduction number R 0 which is qualitatively defined by R 0 = β /infection time. Hence, our model results in the following expression for the basic reproduction number:

Note that for ǫ = 0, the above expression simplifies to R 0 for the simple SIR model. Based on the values that we have accounted in our model, the above expression results to R 0 = 4.38. This estimation is close to the one reported in Russo et al. [17] (R 0 = 4.56) and greater than the ones reported in other studies. For example, Allieta et al. [36] reported R 0 = 3.88, while in another study regarding Italy, D'Arienzo and Coniglio [37] used a SIR model to fit the reported data in nine Italian cities and found that R 0 ranged from 2.43 to 3.10. Finally, [38] reported an R 0 around 3.4. However, we note that the above studies that provide estimates of R o are based on the reported cases which, may lead to an underestimation of the basic reproduction number (see also the discussion in [39] ), while our model takes into account the compartment of asymptomatic cases which transmit the disease.

Regarding the estimation of the effective reproduction number R e , this corresponds to the average number of secondary infections from a single infectious individual during an epidemic. For its calculation, we derive the next generation matrix G with elements g ij formed by the average number of secondary infections of type i from an infected individual of type j, given by (see also [17] ):

F is the vector containing the transmission rates, and V is the vector containing the transition rates between the infected compartments. The effective reproduction number R e is the spectral radius, i.e. the dominant eigenvalue of G. In our model we have:

.

In the spring of 2020, many nations enacted policies restricting resident movement in hopes of slowing the spread of COVID-19. Although the individual compliance to these policies may vary, there is a strong link between government mobility policies and the subsequent mobility behavior in a population [40, 41] . Furthermore, these mobility patterns are correlated with changes in COVID-19 cases [42] .

There are several dates that demarcate major changes in Lombardy mobility policies. On March 4, 2020 all schools and universities were closed. On March 8, 2020, the first major lockdown was introduced, with the encouragement of working from home and closure of recreational venues. On March 21, 2020, measures were furthered; a stay-at-home order was issued along with the banning of pedestrian activity. These measures continued until May 4, when policies were gradually relaxed; non-essential industrial activities reopened, and other pedestrian activities was allowed. These key dates can be used to divide the initial wave of COVID-19 into four approximately time periods. The first period-February 15 to March 8-starts with the availability of mobility data in Lombardy and ends with the introduction of lockdown. The second period stretches from March 9 to March 20: the period of initial lockdown. The third period represents the first phase of stricter measures, ranging from March 21 to April 19, while the fourth and final period begins April 20 and continues until policy relaxation on May 4.

To describe mobility patterns in Lombardy, we used data provided by the Google Community Mobility Reports [2] . The reports show movement trends within a population for a given region. The mobility data is organized into six categories (i.e. activities): retail & recreation (RR), grocery/pharmacy (G), parks (P), transit (T), workplaces (W), and residential (R). We note that within this period schools & universities (SU) were closed. The mobility values are expressed as percent deviation from the baseline value; the baseline was set as the median value from January 3 to February 6, 2020. All activities except for the residential one determine the percent deviation from the change in total visitors. The residential category is determined by the change in time spent in places of residence.

The relative impact of each mobility activity on disease transmission was analyzed to inform the parameter β in the model, since the amount of personal contacts made in a given location can be used to estimate the likelihood of disease transmission in that location. Based on a survey by Mossong et. al. [18] , home, work, and leisure activities each account for about 20% of all personal contacts reported, transport accounts for the least amount of contacts at about 5% and a 20% is accounted to multiple activities. Despite travel being responsible for the least amount of physical and verbal contacts, diminished travel has still been strongly correlated with a decrease in the reproduction number [43] . Therefore, based on this information, each reported mobility category likely plays a significant role in the transmission of COVID-19. To account for the role of each activity in the transmission of COVID-19, an aggregate value was found by appropriately weighting the six categories for a given day. Thus, the transmission rate is impacted by the change in mobility such that β = (1 − u)(1 − θ C )β 0 , with u a measure of the reduction in mobility, i.e., social distancing. In this formulation, when u = 0, there are no mobility restrictions, and when u = 1, there is assumed to be complete lack of mobility and social interaction in the population beyond the residential. The caution that people exhibit in their daily interactions is captured by the term θ C , with θ C = 0 being no caution. Following the first curtails in mobility on March 8 th , the caution was quantified to the value θ C = 0.2. This is in agreement with what has been reported also in other studies about the effect of the cautiousness in the reduction of the disease transmission rate [44] .

Analyzing the available data, the effect of policy changes for different aspects of life in society can thus be quantified as

Obviously, the mobility activity in the category "Residential" is affected by the other activities, including the transit which includes also traveling to other regions for vacation or work. In order to model and assess this dependence, a linear regression model of the following form was fit:

The parameters of the model were estimated using the available mobility data from Feb 15, to May 4 2020. In Table 2 , we report the values of the parameters and their 95% confidence intervals together with the fit statistics, namely the standard error (SE) of the coefficients, the t-statistic to test the null hypothesis that the corresponding coefficient is zero against the alternative that it is different from zero, given the other predictors in the model, the corresponding p-value for the t-statistic of the hypothesis test that the corresponding coefficient is equal to zero or not. The Root Mean Squared Error was 0.0194, the R-squared 0.979, and the adjusted R-squared was 0.978. Thus, at a statistical significance threshold of α = 0.05, the variables RR,G,P,W and S are significant for describing the residential mobility, while the transit activity is not statistically significant. Note also that by the statistics, the hypothesis that c 0 is different than zero cannot be rejected at the significance level of 0.05, highlighting that the residential activity is purely dependent on the rest of the activities and can be discounted in the following section. Employing the statistical analysis results of Table 1 , the resulting standard deviation for social distancing, u(t), was identified to be σ u = 0.0282 and furthermore it could be considered constant during the period Feb. 24 th to May 4 th . Note that during this period there were specific governmental policies assigning mobility curtails to each activity. From this analysis the adherence by Lombardy to the policies can be estimated, which must be accounted for.

Employing the SEASQHRD model with the identified parameters and considering the variance in social distancing due to adherence (Table 1) , we simulated the progression of the disease under the social distancing policy set by the government for 200 adherence scenarios. Figure 2 presents the simulations and the reported data for Lombardy. As the results show, the predictions of the model fits fair the actual reported data. We note that the lower limit of R e on May 4 is close enough to 1. As discussed above, the difference between our estimates with the reported R e in other studies (see e.g. [38] who reported a value of R e around 0.6 just before the end of the lockdown on May 4) is due to the fact our model takes into account also the asymptotic cases assuming that they also transmit the disease. Table 1 

The developed SEASQHRD model of §2.1 was employed to propose government policy as the pandemic progresses. A systematic mechanism for the proposed policy hinges on the solution of an optimization problem that balances the impact on quality of life due to prolonged social distancing on the one hand and the impact of the pandemic on life expectancy on the other. Considerations include the question of societal adherence to government advisories as well as other unaccounted factors which affect the outcome of predictions. We thus assume that government policies need to be continually evaluated and adapted to reality. A system's approach towards this end involves the design of a probabilistic model predictive control (PMPC) structure. The design formulates the policy question as an optimization problem with a finite number of policy decisions being made and their effect being predicted by the SEASQHRD simulator for a finite time horizon forward, called the prediction horizon, T P , given the current state of the epidemic. Each policy decision is effective for a finite period, called the Action Horizon, T A . The number of decisions, N d , multiplied by T A is called the control horizon, T C N d T A . Once a policy has been identified, the first policy decision is enacted and by the end of period T A the state of the epidemic is measured. A new optimization problem for the period T A to T A + T P is then formulated and solved. This recursion constitutes the MPC, where one of the advantages is that it can adapt to uncertainty in the epidemic's behavior. Note, that the state of the compartments is a prerequisite for the SEASQHRD simulator, however only three compartments are directly measured, leaving seven compartment states unknown. In the previous section the analysis of the activities resulted in a distribution of values which impacted the epidemic's evolution. Recognizing the uncertainty in the announced policy due to adherence and the limitations of the derived model, the underlying optimization problem is defined as a probabilistic one α * (t) = arg min The cost function of (2.25) balances the social cost of the epidemic to the social distancing, via the relative weights ω I and Ω A . Function E[·] denotes the expectation of the distribution for argument ·, while P[·] denotes the probability argument · is true. The cost function is constructed to be convex to facilitate convergence to an optimum. Constraint (2.27) reflects the need for mobility of essential workers and for essential activities to take place in order for society to keep functioning. These upper bounds where computed from the Google mobility report data (presented in Table 1 ), by estimating the maximal expected curtail on mobility during the period Feb-24 to May 4, during which society was in essence completely closed.

The stress on society due to mobility restrictions can be significant. Constraint (2.28) is introduced to prevent a rapid increase in the restriction of each activity. In the original formulation, we impose a constraint of 25% increase in each activity restriction with the exception of schools & universities (SU ) which is left unconstrained (i.e., δα c = [0.25 0.25 0.25 0.25 0.25 1.0]). Note that there is no constraint imposed for lifting the restrictions. The societal effect of the desired policy is also captured by the use of two constraints, (2.31) limiting the overflow of patients to hospitals and (2.32) limiting the effective reproduction number, defined in (2.22) (with the basic reproduction number, R 0 calculated from (2.20) ) so that new infection cases dynamic becomes a contraction. In Figure 2 (f), we observe that during the period February 24-May 4, R e remains above one (even though it is greatly reduced) which hints at the increasing rate of asymptomatic cases that transmit the disease during the same period. The use of probabilities in constraints (2.31) & (2.32) captures the degree of risk aversion in the desired policy, by imposing a limit to the risk of violating them via tunable parameters P c bed ∈ [0, 1] and P c rep ∈ [0, 1]. A value of 0.5 converts this constraint to an expectation constraint. For the specific formulation, as the values of P c bed and P c rep increase the PMPC formulation converges to the robust MPC formulation [45] .

The model predictive control simulation begins on February 24 th , three days after the first case was reported in Lombardy. Within this 70 day period until May 4, when the restrictions were relaxed, the horizon lengths determine the duration and frequency of policy investigation and implementation. The prediction horizon, or how far into the future the simulation investigates at a given time point, was examined at lengths of 1 to 4 weeks. As the prediction horizon increases from 1 to 4 weeks, the action during the 1 st period converges to a set value. A prediction horizon of 3 weeks was chosen to balance this convergence and the uncertainty in model predictions. The control horizon, which determines the number of policies during a given prediction, was chosen as 2 weeks to balance problem complexity and solution convergence. Figure S1 of the Supporting Information showcases the effect of the different horizon choices. The action horizon, that is the actual duration of a given policy, was chosen to be 1 week. This choice of horizon allows for frequent policy changes if necessitated by the model predictions, while not being so frequent as to be unrealistic from a policy implementation standpoint.

The choice of the objective function determines the optimization calculations performed by the controller. A linear objective function of (S.3.19) could be used, but the quadratic objective function of (2.25) resulted in stable control pattern and was used in this work. The probabilistic form of this optimization function was used to account for inherent uncertainty in policy implementation; it is less aggressive than a robust formulation of (S.3.10)-(S.3.18), which assumes the worst-case scenario during optimization. Consequently, the choice of a PMPC formulation allows for less heavy-handed policies while still considering the uncertainty associated with population adherence. Another consideration is the relative cost assigned to mobility restrictions, scaled by the relation between parameters Ω A and ω I & ω R . Focusing first on the relative cost of activity restrictions as reflected by the Google mobility report data during the period Feb-24 to May-4, was used to elucidate the preference given in the restrictions. Specifically, the need for groceries is central for a society where activities such as restaurants are curtailed. Similarly, transit is needed in order for basic social activities to take place. On the other hand, parks are not explicitly needed for society to function during epidemic conditions. Using similar arguments a convenient diagonal matrix Ω A was defined, with the diagonal elements expressing the relative resistance to curtails being ω A = [0.2 1.0 0.5 0.5 1.0 0.5]. For a prediction horizon of 3 weeks, control horizon of 2 weeks, and action horizon of 1 week, a selection of Ω A F → 0 results in government policies of aggressive lockdown. A selection of Ω A F ≈ ω R results in moderate mobility rstrictions. To equally weight the impact on hospitalization and impact on disease progression we employed ω I ≈ 10ω R . Figure  3 shows the simulation results of disease progression for the period Feb-24 to May-4 under the optimal policy as identified by the MPC formulation of (2.25)-(2.33) with ω R = 10, ω I = 10, Ω A as discussed previously, T p = 21, T c = 14, T a = 7 days. The uncertainty in the simulation is captured via presenting 200 population adherence scenaria(i.e., N cases = 200), informed from Table 1 . These parameters ultimately prescribe a government policy that corresponds to approximately a 60% reduction in mobility. The difference in outcome between the implementation of the prescribed restrictions with the historical data from Lombardy for reference by comparing Figure 3 with Figure 2 . The simulation of optimal restrictions leads to less infections, hospitalizations, and deaths than the historical data. A major contributor to the difference in deaths is the date of introduction of major curtails which in our simulation starts on Feb-24 which also has the effect that population caution becomes θ C = 20% on that day. The historical data in contrast show the effect of casual curtails starting on Feb-21 with major (and aggressive) curtails starting on March 8; population caution becomes θ C = 20% on March 8.

One of the more challenging aspects of addressing the COVID-19 crisis is the scale of response, promptness, and fast adaptation to a number of factors ranging from most importantly the strengthening of the public health system both in staff and infrastructures with a particular emphasis on the enhancement of the prevention and primary care facilities network to the design of vaccination policies and from the design of public awareness campaigns to the implementation of appropriate social distancing policies that are needed to combat efficiently this global threat. In an emergent situation in which the inherent uncertainly that characterizes the reported data and generally our knowledge at least at the first stages of an pandemic, mathematical modeling and control theory play an important role towards the simulation and investigation of what if scenarios.

Here, retrospectively, we attempted to assess the impact of social distancing measures on the various activities on the evolution of the pandemic in Lombardy during the first wave of the pandemic (Feb 21-May 4 2020) and its consequence on the fatalities and the pressure on the hospital structures, following a model predictive control approach. Towards this aim, we addressed a compartmental modelling approach that quantifies the uncertainty of the level of asymptomatic cases [46, 17] in the population and also the effect of social contacts in different activities to the spread of infectious diseases, In particular, the proposed model extends the one proposed by Russo et al. [17] for the region of Lombardy, thus considering the hospitalized and guaranteed cases. Furthermore, the effect of the social contacts in various activities is modelled based on the Google mobility reports for the region of Lombardy as well as by metanalyses of the social contacts and their effect in the spread of the infectious diseases as reported for Italy by Mossong et al. [18] . Thus, the proposed PMPC scheme models also the effect of the reported activities, namely, workplaces, grocery and pharmacy, retail and entertainment, parks and transit to the residential change of mobility. This was achieved by fitting a linear regression model based on the corresponding Google mobility reports. We show that the proposed modeling framework approximates quite well the actual evolution of the pandemic in terms of reported infected cases, deaths and hospitalized cases. By activating the PMPC scheme from February 21, 2020, i.e. the day that the first official case was reported in Lombardy, we show that the obtained number of deaths would be significantly less (around 4,000 at May 4 instead of the actual number of 14,000), however the spread of the disease would not have been eradicated by May 4. Indeed our simulations show that due to the presence of a significant portion of asymptomatic cases which transfer the disease but are rarely reported, due to the mild or no symptoms, the effective reproduction number is still above 1 on May 4. However, the pressure in the hospitals would be much less, with a pick at a maximum of 3500 hospitalized cases around the mid of April, instead of the actual 13000 cases in the same period. Of course, as our study is about a retrospective analysis, it should be regarded as indicative of what would happen if we knew a-priori the scale and severity of the threat. Thus, our analysis should not be considered by any means as a "ground truth" of any kind of accurate prediction of what would happen if such measures would have been adopted. However, the results are indicative and highlight the importance of the need of very prompt responses and alertness for analogous emergent and re-emergent infectious diseases that will certainly occur into the future.

On the other hand and importantly, our study pinpoints the urgent need of strengthening the public structures including the epidemiological surveillance mechanisms that will allow faster responses, the enhancement of the (fragmented in most of the countries of the world) primary care facilities networks, the investment in better both qualitatively and quantitatively public hospitals and hiring of medical staff but also in other structures such as schools and mass public transportation.

The proposed framework can be extended at various fronts. One issue at the forefront of current policy making is vaccination, which this model does not address. To accurately forecast policy decisions today, the model would have to be adjusted for vaccine-induced immunity to COVID-19. Other considerations include the probability of reinfection with the emergence of the virus variants, which was a less significant factor in this study due to the time period examined. Furthermore, the proposed approach could be used for the quantification of the uncertainty of the evolving dynamics, taking into account the reported, from clinical studies, distributions of the epidemiological parameters rather than their expected values and a more detailed compartmental model that takes also the age distribution in the population (see e.g. [47] . A critical point that is connected with the above is that with such a small number of infectious cases at the initial stage of the simulations, a stochastic model or a hybrid stochastic model could be more realistic in which uncertainty could be modelled in the form of realistic perturbations (see for example [48] ).

Another challenge posed by the COVID-19 pandemic is the accurate state estimation. Here, in order to set the values of the model parameters, we used the reported values of the various epidemiological parameters. However, to examine longer periods of time, methods like moving horizon estimation and design of stochastic state observers, could be used to more accurately capture and assess the uncertainty in the state of the pandemic. Finally, one could extend the proposed framework by extended as a network of regions and interventions (see [15] . This is an important aspect that we aim to tackle in a future work.

All the relevant data used in this paper are publicly available and accessible at https://lab.gedidigital.it/gedi-visual/2020/coronavirus-i-contagi-in-italia/. The reported cumulative numbers of cases from February 21 to March 19 are listed in Table S1 . The data from February 21 to March 8 have been used for the calibration of the model parameters and the data from March 9 to March 19 have been used for the validation of the model.

All the mobility data used in the present study are publicly available at https://www.google.com/covid19/mobility/. These data were employed to estimate the standard deviation and bias between the government guidelines and population behavior, discussed in §2.3.

Beyond the PMPC formulation of (2.25)-(2.33), deterministic and robust MPC problems, ((S.3.1)-(S.3.9) and (S.3.10)-(S.3.18) respectively, were also formulated and their performance was evaluated. The deterministic problem does not consider population adherence to government guidelines, assuming perfect response. This leads to a simpler formulation at the cost of a potentially optimistic view of the progression of the epidemic that can lead to the violation of constraints (S.3.9) and importantly of (S.3.8).

α * (t) = arg min 

At the other end of the spectrum, the robust formulation considers bounded uncertainties in the state evolution and the enforcement of the control action. It attempts to identify a policy that will not violate the constraints for all investigated scenaria. It is important to note that basic robust MPCs are designed to address problems where the uncertainty is bounded; elaborate RMPC formulations have been designed using Lyapunov arguments [49, 45] . In the presented work, this is addressed by creating a set of "scenaria"; each scenario contains values of the adherence of the population to each activity curtail and the associated expected initial state X 0 that will be used by the SEASQHRD model. A robust optimization problem was also formulated where a non-convex set was defined containing value ranges for each activity and unmeasurable initial state. As the results from these two robust MPC formulations were almost identical, we only present the first RMPC formulation here. 

subject to 

Numerous simulations were performed to investigate the effect MPC parameters have on proposed action. Figure S1 showcases the effect of the MPC horizon parameters on the calculated optimal action. We chose horizons T p = 21, T c = 14, T a = 7 days that balanced problem complexity and accuracy well.

The effect of the choice of weights of Ω A can be visualized for the case when ω I = 0 & ω R = 1, i.e. they are relatively small compared to Ω A F =≈ 1.6. We observe that when Ω A = 0.6I 6×6 stringent curtails are imposed on important actives such as groceries and work (shown in Figure S3 ) compared to when the "variable" weight is employed (shown in Figure S2 ), with the aggregate curtail being u(t) ≈ 0.57. Considering the societal effects of activities the nominal weight for the optimization used in formulation of (2.25)-(2.33) was chosen as 

Johns Hopkins Center for Health Security. Coronavirus COVID-19 Global Cases by Johns Hopkins CSSE

Preliminary estimation of the basic reproduction number of novel coronavirus (2019-nCoV) in china, from 2019 to 2020: A data-driven analysis in the early phase of the outbreak

COVID-19 and italy: what next? The Lancet

Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in wuhan, china: a modelling study. The Lancet

Athanasios Tsakris, and Constantinos Siettos. Data-based analysis, modelling and forecasting of the covid-19 outbreak

Analysis of a model for coronavirus spread

The effect of travel restrictions on the spread of the 2019 novel coronavirus (COVID-19) outbreak. Science, page eaba9757

Covid-19 and Flattening the Curve: A Feedback Control Perspective

Robust and optimal predictive control of the COVID-19 outbreak

Control of COVID-19 system using a novel nonlinear robust control algorithm

Nonlinear model predictive control with logic constraints for COVID-19 management

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

A network model of italy shows that intermittent regional strategies can alleviate the covid-19 epidemic

Decreasing hospital burden of covid-19 during the first wave in regione lombardia: an emergency measures context

Tracing day-zero and forecasting the covid-19 outbreak in lombardy, italy: A compartmental modelling and numerical optimization approach

Social Contacts and Mixing Patterns Relevant to the Spread of Infectious Diseases

Retrospective analysis of the italian exit strategy from covid-19 lockdown

Estimating the scale of covid-19 epidemic in the united states: Simulations based on air traffic directly from wuhan, china. medRxiv

The Novel Coronavirus Pneumonia Emergency Response Epidemiology Team. The epidemiological characteristics of an outbreak of 2019 novel coronavirus diseases (covid-19) in china

Il messaggero. Coronavirus, stop ai test facili: tampone solo a chi ha i sintomi

Perché non stiamo facendo più tamponi?

Comparison of Symptoms and RNA Levels in Children and Adults With SARS-CoV-2 Infection in the Community Setting

Early transmission dynamics in wuhan, china, of novel coronavirus infected pneumonia

Temporal dynamics in viral shedding and transmissibility of covid-19

Characteristics of SARS-CoV-2 patients dying in Italy Report based on available data on

The early phase of the covid-19 epidemic in lombardy

Epidemiology and transmission of COVID-19 in 391 cases and 1286 of their close contacts in shenzhen, china: a retrospective cohort study. The Lancet Infectious Diseases

Clinical progression of patients with COVID-19 in shanghai

Characteristics of COVID-19 patients dying in Italy Report based on available data on

Covid-19 in italy: Considerations on official data

Survival of hospitalized covid-19 patients in northern italy: A population-based cohort study by the ita-covid-19 network

Genomic characterisation and phylogenetic analysis of sars-cov-2 in italy

Inners and stability of dynamic systems

Covid-19 outbreak in italy: estimation of reproduction numbers over 2 months prior to phase 2

Assessment of the sars-cov-2 basic reproduction number, r0, based on the early phase of covid-19 outbreak in italy

Tracking the time course of reproduction number and lockdown's effect on human behaviour during sars-cov-2 epidemic: nonparametric estimation

A new framework and software to estimate time-varying reproduction numbers during epidemics

Do covid-19 policies affect mobility behaviour? evidence from 75 canadian and american cities

Effects of information-induced behavioural changes during the covid-19 lockdowns: the case of italy

Association between mobility patterns and COVID-19 transmission in the USA: a mathematical modelling study. The Lancet Infectious Diseases

Using mobility to estimate the transmission intensity of COVID-19 in Italy: A subnational analysis with future scenarios. medRxiv

Quantifying social distancing arising from pandemic influenza

Robust hybrid predictive control of nonlinear systems

Coronavirus -scientific insights and societal aspects

Effects of latency and age structure on the dynamics and containment of covid-19

Bounded Noises in Physics, Biology, and Engineering

Predictive control of switched nonlinear systems with scheduled mode transitions

The work was supported by the Ministero dell'Università e della Ricerca through the program of Fondo Integrativo Speciale per la Ricerca (FISR).Program: FISR2020IP-02893.

The authors declare there is no conflict of interest.