key: cord-0521480-ptpiu3kl
authors: Morato, Marcelo M.; Pataro, Igor M. L.; Costa, Marcus V. Americano da; Normey-Rico, Julio E.
title: A Parametrized Nonlinear Predictive Control Strategy for Relaxing COVID-19 Social Distancing Measures in Brazil
date: 2020-07-19
journal: nan
DOI: nan
sha: f473e4d76cade514cd640433f5e0637d332db159
doc_id: 521480
cord_uid: ptpiu3kl

In this paper, we formulate a Nonlinear Model Predictive Control (NMPC) to plan appropriate social distancing measures (and relaxations) in order to mitigate the COVID-19 pandemic effects, considering the contagion development in Brazil. The NMPC strategy is designed upon an adapted data-driven Susceptible-Infected-Recovered-Deceased (SIRD) contagion model, which takes into account the effects of social distancing. Furthermore, the adapted SIRD model includes time-varying auto-regressive contagion parameters, which dynamically converge according to the stage of the pandemic. This new model is identified through a three-layered procedures, with analytical regressions, Least-Squares optimization runs and auto-regressive model fits. The data-driven model is validated and shown to adequately describe the contagion curves over large forecast horizons. In this model, control input is defined as finitely parametrized values for social distancing guidelines, which directly affect the transmission and infection rates of the SARS-CoV-2 virus. The NMPC strategy generates piece-wise constant quarantine guidelines which can be relaxed/strengthen as each week passes. The implementation of the method is pursued through a search mechanism, since the control is finitely parametrized and, thus, there exist a finite number of possible control sequences. Simulation essays are shown to illustrate the results obtained with the proposed closed-loop NMPC strategy, which is able to mitigate the number of infections and progressively loosen social distancing measures. With respect to an"open-loop"/no control condition, the number of deaths still could be reduced in up to 30 %. The forecast preview an infection peak to September 2nd, 2020, which could lead to over 1.5 million deaths if no coordinate health policy is enacted. The framework serves as guidelines for possible public health policies in Brazil.

The SARS-CoV-2 virus was first registered in humans in Wuhan, China by December 2019. Since then, it has reached almost all countries around the globe with pandemic proportions and devastating effects. This contagion seems unprecedented and can already be rules as the definite health crisis of the 21 rst century. The SARS-CoV-2 virus causes a severe acute respiratory syndrome, which can become potentially fatal. (SIRD) models Peng et al. (2020) ; Kucharski et al. (2020) . These SIRD models comprise coupled nonlinear differential equations, as originally presented for population dynamics in work of Kermack & McKendrick (1927) . Therefore, the Nonlinear Model Predictive Control (NMPC) Camacho & Bordons (2013) framework shows itself as a rather convenient approach to guide model-based public health policies, given the fact that it can adequately consider the nonlinear dynamic of the virus spread (through these SIRD equations) together with the effect of lockdown/quarantine measures, introduced as constraints of the optimization problem. Recent literature has indeed shown how SIRD-based NMPCs can be applied to plan COVID-19 social distancing measures in three novel works: • Alleman et al. (2020) consider the application of a NMPC algorithm regarding the data from Belgium.

The control input is the actual isolation parameter that is plugged into the SIRD model;

• Morato et al. (2020) consider an optimal On-Off MPC design, which is formulated as mixed-integer problem with dwell-time constraints. Furthermore, the response of the population to social isolation rules with an additional dynamic variable, which is then incorporated to the SIRD dynamics;

• Köhler et al. (2020) applied the technique for the German scenario. The approach considers that the control input is a variable that directly affects the infection and transmission rates (parameters of the SIRD model).

Building upon this previous results, herein we propose an NMPC algorithm formulated on the basis of a SIRD-kind model, which is obtained by the means of an identification algorithm. Following the lines of Morato et al. (2020) , we assumed that the control action is a social isolation guideline, which is enacted and passed on to the population. Then, the population responds in some time to these measures, which can be mathematically described as a dynamic social isolation factor. Furthermore, in accordance with the discussion presented by Köhler et al. (2020) , we also introduce an additional dynamic nonlinear model, which gives the input/output (I/O) relationship between the COVID-19 contagion infection and transmission rates according to the stage of the pandemic, as suggested by . Complementary, we use autoregressive equations for the epidemiological parameters of the disease (transmission factor, infection factor, lethality rate) directly related to social distancing as control input. These time-varying regressions are used to provide better long-term forecasts than the regular SIRD equations. Regarding the NMPC framework, the novelty in this paper resides in the formulation of the control input as a finitely parametrized variable, which makes the NMPC implementation persuadable through search mechanisms, which run much faster than the Nonlinear Programming methods from the prior.

Motivated by the previous discussion, the problem of how optimal predictive control can be used to formulate adequate social distancing policies, regarding Brazil, is investigated in this work. The main contribution of our approach comprises the following ingredients:

• Firstly, an adapted SIRD model is proposed, which incorporates delayed auto-regressive dynamics for the transmission, infection and lethality rates of the SARS-CoV-2 virus, which vary according to the stage of the pandemic. The adapted model also incorporates the dynamic response of the population to a given social distancing guideline, which is the NMPC control input.

• In order to identify the SIRD model parameter auto-regressive dynamics (for the epidemiological parameters), we use a three-layered identification procedure, which concatenates analytical expansions and Least-Squares optimization procedures. The identified models are validated with regard to the available set of data from Brazil, which includes the social distancing factor observed in the country.

• The main innovation of the paper is the proposed Parametrized Nonlinear Predictive Control algorithm, which is based on the identified adapted SIRD model with auto-regressive epidemiological parameters. The control input (social distancing guideline) is finitely parametrized over the NMPC prediction horizon, at each sampling instant. Then, an explicit nonlinear programming solver is simulated for all possible input sequences along time. The NMPC solution, then, is simply found through a search mechanism regarding these simulated sequences, which is rather numerically cost-efficient.

• Finally, we present results considering the application of this NMPC algorithm to the COVID-19 scenario in Brazil. These results are a twofold: a) those that regard the application of the optimal strategy to control the pandemic since its beginning, comparing the simulation results to those seen in practice; and b) applying the control method from 30th of July onward, aiming to mitigate and revert the current health crisis catastrophe.

This paper is structured as follows. Section 2 presents the new SIRD model with auto-regressive epidemiological parameters, the respective identification procedure and validation results. Section 3 discusses the proposed NMPC strategy. Section 4 depicts the obtained control results, regarding the COVID-19 contagion mitigation for Brazil. General conclusions are drawn in Section 5.

In this Section, we describe the used modeling framework for the COVID-19 contagion dynamics in Brazil. We build upon the SIRD models for the SARS-CoV-2 virus from works of Peng et al. (2020) ; Kucharski et al. (2020) ; Ndairou et al. (2020) . The model adaptations are done such that the epidemiological parameters are taken as time-varying piece-wise constant functions of the stage of the pandemic. This adaptation is in accordance with recent immunology results ; Dowd et al. (2020) ; He et al. (2020) , which discuss the transmission and reproduction rate of this virus. It is also taken into account two different sampling periods, to include the issue of the incubation period of the virus in human, reportedly, in average, as 5.1 days Lauer et al. (2020) .

The SIRD model describes the spread of a given disease with respect to a population split into four non-intersecting classes, which stand for:

• The total amount of susceptible individuals, that are prone to contract the disease at a given (discrete) sample of time k, denoted through the dynamic variable S(k);

• The individuals that are currently infected with the disease (active infections at a given sample of time k), denoted through the dynamic variable I(k);

• The total amount of recovered individuals, that have already recovered from the disease, from an initial instant k 0 until the current sample k, denoted through the dynamic variable R(k);

• Finally, the total amount of deceased individuals due to a fatal SARS-CoV-2 infection, from an initial instant k 0 until the current sample k, is denoted D(k).

Due to the evolution of the spread of the disease, the size of each of these classes change over time. Therefore, the total population size N is the sum of the first three classes as follows:

Since the Brazilian government discloses daily samples of total infections and accumulated deaths, we consider that these discrete-time dynamics samples k, given each T 1 = 1 day. Furthermore, to account for the average incubation period of the disease, we consider that the epidemiological parameters vary weekly (each T 2 = 7 days), kept as zero-order-held/piece-wise-constant samples. We will further assess this matter in the sequel.

In the SIRD model, there are three major epidemiological parameters, which express the specific dynamics of the SARS-CoV-2 virus in the population set. The dynamics of these parameters are given in the sparser discrete-time weekly samples, since they change according to the incubation of the virus:

• The transmission rate parameter β stands for the average number of contacts that are sufficient for transmission of the virus from one individual. According to the detailed classes of individuals, then, it follows that T 1 β(k)I(k)/N (k) determines the number of contacts that are sufficient for transmission from infected individuals to one susceptible individual; and (T 1 β(k)I(k)/N (k))S(k) determines the number of new cases (per day) due to the amount of S(k) susceptible individuals (these are the ones "available for infection").

• The infectiousness Poisson parameter γ denotes the inverse of the period of time a given infected individual is indeed infectious. Consequently, γ affects the rate of recovery (or death) of an infected person. This parameter directly quantifies the amount of individuals that "leaves" the infected class, in a given sample.

• The mortality rate parameter ρ stands for the observed mortality rate of the COVID-19 contagion.

We model the amount of deceased individuals due to the SARS-CoV-2 infection following the lines of (Keeling & Rohani, 2011) : the new number of deaths, at each day, can be accounted for through the following expressions:

Considering these explanations, the "SIRD" (Susceptible-Infected-Recovered-Dead) model is expressed through the following nonlinear discrete-time difference equations:

where I S denotes the portion of the infected individuals which in fact display symptoms. This "symptomatic" class has been accounted for in previous papers ; Morato et al. (2020) . We note that only these symptomatic individuals will require possible hospitalization and may die. The remainder (1 − p sym )I are asymptomatic or lightly-symptomatic, which do transmit the virus but do not die or require hospitalization. The class of recovered individuals considers immunized people, encompassing those that displayed symptoms and those that did not. The symptomatic ration parameter p sym is constant and borrowed from previous papers. The cumulative number of cases I c denotes the total number of people that have been infected by the SARS-CoV-2 virus until a given day. In fact, this cumulative variable is equivalent to those that are currently infected I, summed with those that have recovered R and those that have died D.

We must also stress that, in the SIRD model equations given above, parameter ψ represents a transmission rate mitigation factor. This factor accounts for the observed social isolation factor with the population set N . It follows that ψ = 0 denotes the situation where the whole population set has sustained social interactions. As discussed in previous papers (Keeling & Rohani, 2011; , there exists some "natural" ψ = ψ factor, which stands for a normal/nominal conditions, still with "no control" of the viral spread (with no social isolation guidelines, this kind of situation has been observed in Brazil in the first weeks of the contagion, end of February, beginning of March). In contrast, ψ = 1 represents a complete lockdown scenario, where there are no more social interactions (this conditions has not been seen and is, in practice, unattainable). In practice, there is some maximal social isolation factor ψ = ψ that can be put in practice. In this paper, we consider the values for social isolation in Brazil from the work of , which discloses weekly estimates for ψ.

Another essential information in epidemiology theory is the basic reproduction number, usually denoted by R 0 . This factor is able to measure the average potential transmissibility of the disease. In practice, this value is fixed/constant and inherent to a given disease. Through the sequel, a time-varying version of this basic reproduction number is considered, denoted R t , which represents how many COVID-19 cases could be expected to be generated due to a single primary case, in a population for which all individuals are susceptible. From a control viewpoint, R t represents the epidemic spread velocity: if R t > 1, the infection is spreading and the number of infected people increases along time (this typically happens at the beginning of the epidemic), otherwise, if R t < 1, it means that more individuals "leave" from the infected class, either recovering or dying, and, thereby, the epidemic ceases. The reproduction number R 0 is affected by different factors, including immunology of the virus, biological characteristics, but also governments policies to control the number of susceptible population. Since the epidemic parameters change along the time, we process by considering the effective reproduction number R 0 also as a dynamic variable.

The underlying assumption used to calculate R t , is that, at the beginning of the pandemic, S ≈ N . Considering parameters β, γ, ρ and ψ from the SIRD model Eq.

(2), R t can be roughly given by:

Remark 1. Regarding the SIRD model given in Eq.

(2), it follows that N (0) = N 0 is the initial population size (prior to the viral infection). Furthermore, we stress that, in SIRD-kind models, I(k) represents the active infections at a given moment, while D(k) represents the total amount of deaths until this given moment; for this reason, it follows that D(k + 1) − D(k) is proportionally dependent to I(k).

Remark 2. We note that we are not able to use more "complex" descriptions of the COVID-19 contagion in Brazil, such as the "SIDARTHE" model used by Köhler et al. (2020) (that splits the infections into (symptomatic, asymptomatic) detected, undetected, recovered, threatened and extinct) because we have insufficient amount of data. The Ministry of Health only disclosed the total amount of infections (I(k) + D(k) + R(k)) and the total amount of deaths (D(k)) at each day. Due to the absence of mass (sampled) testing, there is no data regarding detected asymptomatic individuals, for instance, as it is available in Germany (where (Köhler et al., 2020) originate). To choose a more complex model may only decrease the truthfulness/validity of the identification results, since the parameters can only represent a singular combination that matches the identification datasets and cannot be used for forecasting/prediction purposes.

To illustrate the dynamics of a contagion like the pandemic outbreak of COVID-19, in Figure 3 we present long-term forecast using a regular SIRD model with constant epidemiological parameters, following the methodology presented by . These predictions were computed on 11/06/2020 and are here shown to demonstrate the qualitative evolution of the I and D curves. The active infections curve I shows a increase-peak-decrease characteristic, while the total number of deaths D shows an asymptotic behavior to some steady-state condition. We note that, as of 30/06, the number of infections had already significantly increased, which means that the most recent forecasts preview an even worse catastrophe.

In this paper, the SIRD dynamics from Eq. (2) are adapted. These adaptations are included for three main reasons:

1. In accordance with the discussions presented by Morato et al. (2020) , we understand that it is imperative to embed to the SIRD model how the population responds to public health policies, that is, to incorporate the dynamics and effects of social distancing. When a government enacts some social distancing policy, the population takes some time to adapt to it and to, in fact, exhibit the expected social distancing factor ψ. This is, to take ψ as a time-varying dynamic map of the enacted social distancing guidelines (which will, later on, be defined by an optimal controller).

Apr 2020 It is reasonable that when more active infections occur, the virus tends to spread more efficiently, for instance. Moreover, the mortality rate should stabilize after the infection curve has decreased. To incorporate this issue in the model, auto-regressive moving-average time-varying dynamics for the viral parameters are proposed, which stabilize according to stage of the pandemic. Furthermore, by doing so, the inherent incubation characteristic of the virus is also taken into account, since the viral transmission, lethality and infection rates should vary with dynamics coherent with the average incubation period. (2020); Bhatia et al. (2020) for the evolution of the virus assuming that the epidemiological parameters remain constant along the prediction horizon. This is not reasonable, due to the fact that such parameters are inherently time-varying, given according to the exhibited pandemic level. Thus, if parameters are held constant for forecasts, these forecasts are only qualitative and allow short-term conclusions. Since dynamic models for the epidemiological parameters are proposed, make more coherent long-term forecasts can be derived. Of course, it must be stressed that these forecasts will still be qualitative, since one cannot account for perfect accurateness regarding the number of infection and deaths due to the absence of mass testing in Brazil. Furthermore, the effect of future unpredicted phenomena cannot be accounted for, such as the early development of an effective vaccine, which would certainly make the infections drop largely.

These model adaptations are discussed individually in the sequel.

In order to design and synthesize effective control strategies for social distancing (public) policies, to be oriented to the population by local governments, the social distancing factor ψ is further exploited.

In what concerns the available data from Brazil (that is used for identification of the model parameters), the social distancing factor ψ is a known variable, given in weekly piece-wise constant samples. Later on, regarding the control procedure (Section 4), this isolation factor will vary according to the enacted social distancing policy u, as defined by a nonlinear optimal predictive control algorithm.

The differential equation that models the response of the susceptible population to quarantine guidelines is taken as suggests Morato et al. (2020) , this is:

where u(k) is the actual control input, the guideline that defines the social isolation factor goal, as regulated by the government (this signal will be later on determined by the proposed optimal controller), ψ = 0.4317 day −1 is a settling-time parameter, which is related to the average time the population takes to respond to the enacted social isolation measures, and K ψ is a time-varying stating gain relationship between ψ and u.

As recommend Morato et al. (2020) , we assume that when more infected cases have been reported, and when the hospital bed occupation surpasses 70%, the population will be more prone to follow the social distancing guidelines, with larger values for the gain relationship K ψ :

being n ICU the total number of ICU beds in the country. We recall that p sym is a parameter which gives the amount of infected individuals that in fact display symptoms and may need to be hospitalized (I s = p sym I).

Remark 3. We note that Brazil has, as of February, 45, 848 ICU beds. This number is estimated to increase in up to 80 % with field hospitals that were built specifically for the COVID-19 contagion (i.e. n ICU = 82, 526). The percentage of symptomatic individuals is taken as 16 %, according to the suggestions of . This value is coherent with the available information regarding this virus, concerning multiple countries Lima (2020) . The percentage of symptomatic considered groups the infections with severe/acute symptoms, which will, indeed, most possibly require treatment.

Remark 4. In practice, Eq. (4) is bounded to the minimal and maximal values for the social distancing factor ψ and ψ, respectively. These values are the same are that used as saturation constraints on u, as discussed in the sequel.

We proceed by considering that u is a finitely parametrized control input. This is: the enacted social distancing guidelines can only be given within a set of pre-defined values. This approach is coherent with possible ways to enforce and put in practice social distancing measures. Pursuing this matter, the actual control signal, that is to be defined by the proposed optimal automatic controller, must abide to a piecewise constant characteristic with a possible increase/decrease of 5% of isolation per week. This percentage is the average increase of social isolation, in the beginning of the pandemic in Brazil, when social distancing measures were strengthened significantly over a small time period. The values for social isolation ψ(k) in Brazil have been estimated by InLoco (2020) and are presented in Figure 4 .

Moreover, u is considered to be defined within the admissibility set [ψ , ψ]. As gives Figure 4 , which shows the observed social isolation factor ψ in Brazil, the "natural" isolation factor is of ψ = 0.3 (which stands for a situation of no-isolation guidelines). Furthermore, ψ is the maximal attainable isolation factor for the country; the maximal observed value in Brazil is of roughly 53 % (see Figure 4) . Anyhow, coordinated "lockdown" measures were not forcefully enacted. For this matter, we consider ψ = 0.7, which would represent that the population can be restricted to, at most, a 70% reduction on the level of social interactions. As reported by , it seems unreasonable to consider larger values for social isolation in the country, due to multiple reasons (hunger, social inequalities, laboring necessities, lack of financial aid from the government for people to stay home, etc.). Remark 5. It must noted that the SIRD model considers some "natural" isolation factor (which is the lower bound on u). As seen in Figure 4 , the social isolation factor of 30 % is the observed lower bound for ψ seen in Brazil.

Mathematically, these constraints are expressed through the following Equation set: 

This constraints given in Eq. (6) are, in fact, very interesting from an implementation viewpoint, as the government could convert the finitely parametrized values for u into actual practicable measures, as illustrates Table 1 . We remark that this Table is only illustrative. Epidemiologists and viral specialists should be the ones to formally discuss the actual implemented measures that ensure the social distancing factor guideline given by u.

The second main modification of the original SIRD model is to consider dynamic models for the epidemiological parameters of the SARS-CoV-2 virus, γ, β and ρ. These models are taken as auto-regressive moving average functions, which converge as the pandemic progresses. Therefore, the following dynamics are considered:

The models given in Eq. (7) are possibly delayed and auto-regressive. Anyhow, despite the parameters β, γ and ρ being time-varying, the model functions f β , f γ and f ρ are constant. The order and number of regressions are found through an optimization procedure, which is further detailed in Section 2.3. 

The complete model used in this work to describe the COVID-19 contagion outbreak in Brazil is illustrated through the block-diagram of Figure 5 . We note that the "COVID-19 Epidemiological Parameter Models", represented through Eq. (7), enable the complete model to offer long-term predictions with more accurateness, as discussed in the beginning of this Section. We denote as the "SIRD+ARIMA" model, henceforth, the cascade of the population response from Eq. (4) to the SIRD model, with auto-regressive time-varying epidemiological parameters as give Eq. (7). (2020). The majority of these papers indicates that SIRD models provide the best model-data fits.

We note, anyhow, that the SIRD model offers three degrees-of-freedom at each instant k (i.e. β(k), γ(k) and ρ(k)), as gives Equation 2. This means that different instantaneous combinations of these parameters can yield the same numerical values for S(k), I(k), R(k) and D(k). Therefore, although mathematical and graphical criteria have been used to validate these dynamic models when compared to real data, the estimated values for these parameters should be coherent with biological characteristics of this viral pandemic. An indicative of badly adjusted SIRD model parameters is the effective reproduction number of the disease R 0 , which should naturally decrease to a threshold smaller than one as the pandemic ceases.

Thus, in order to identify a SIRD model that is indeed able to describe the pandemic contagion accurately (which is especially important for forecast and control purposes), we proposed a two step identification method to determine the time-varying dynamics for β(k), γ(k) and ρ(k), with respect to the dynamic parameter models in Eq. (7).

Regarding datasets, we must comment that the Brazilian Ministry of Health discloses, daily, the values for the cumulative number of confirmed cases I c (k) and the cumulative number of deaths D(k). Furthermore, we recall that, through the sequel, the values of the observed average social isolation in the country, ψ(k), are assumed to be known (as given in Figure 4 ). The symptomatic percentage is assumed as constant along the prediction horizons.

Then, we proceed by depicting the used identification procedure, which is performed in three consecutive steps/layers, as follows:

• The first step resides in analytically solving the SIRD regressions from Equation 2 for a fixed interval of points (days).

• Then, the derived values from the analytical solutions are passed as initial conditions to an optimization layer, which solves a constrained Ordinary Least Square minimization problem, aiming to adjust the parameter values to pre-defined (biologically coherent) sets, so that the identified model yields small error with adequate parameter choices. The output of this second step stands for the time-series vectors regarding the SIRD parameters.

• Finally, these time-series are used to fit auto-regressive models, via Least Squares, as gives Eq. (7).

By following this three-layered procedure, the estimated parameter equations (Eq. (7)) are found in accordance with biological conditions, which are described as the pre-defined optimization sets. We note that we embed to the optimization layer sets that are also in accordance with previous results for SIRD model estimations for Brazil ; Morato et al. (2020) . The complete identification procedure used to estimate the SIRD model dynamics follows the lines illustrated in Figure 6 .

We note that, in order to simplifying the formulation of the optimization algorithm, the difference Equations for I(k + 1) and D(k + 1) are modified in order to decouple parameters related to the number of fatal and related to the number of recovered individuals, with regard to I(k). This is:

where α(k) = ρ(k) γ(k) (1−ρ(k)) . Therefore, instead of identifying β, γ and ρ, we pursue the identification of β, γ and α, and, then, ρ is computed as follows:

The identification procedure starts by collecting the available datasets (regarding the cumulative number of COVID-19 cases I c and deaths D) inside a fixed interval k ∈ [k i , k f ]. Then, the first layer computes "exact", analytical values for the epidemiological parameters, denotedβ,γ andα according to the following discrete analytical expansions:

(13)

It is assumed that the used datasets might be corrupted by series of issues, such as cases that are not reported on given day k and simply accounted for on the following days, or sub-reported cases, as discussed by . Therefore, we adjust these parameters through the optimization layer, in order to improve the reliability of the identification.

Once again, the available data from the same interval discrete [k i , k f ] is used. The optimization procedure minimizes a constrained Ordinary Least Squares problem, whose solution comprises the following vectors: S, I, R and D. These vectors comprise the model-based outputs within the considered discrete interval. The decision variables for the optimization procedure are the degrees-of-freedom of the SIRD model, in this layer denotedβ,γ,α. The Ordinary Least Square criterion ensures that the optimization minimizes the error between the real data and the estimated values from the SIRD model, by choosing coherent values for the decision variables.

The quadratic model-data error variables terms used in the optimization layer are the following:

Er R (i) = (R(i) −R(i,β(1 − ψ),γ,α)) 2 ,

for which the variablesÎ,R,D are estimated according to the SIRD model equations. The complete optimization problem is formulated as follows:

wherein δ and δ are uncertainty interval margins used to define the lower and upper bound of each decision variable on the optimization problem, β 0 , γ 0 and α 0 are the initial conditions of the optimization problem and w 1 , w 2 and w 3 are taken as positive weighting values (tuning parameters), used to normalize the magnitude order of the total cost and adjust the model fit with respect to variables Er I , Er R and Er D . With respect to the discussion regarding the average incubation period of the SARS-CoV-2 virus, the optimization procedure compute piece-wise constant values for β(k), γ(k) and α(k) each T 2 = 7 days. This is also done to improve the model-data fitting results, as suggest Morato et al. (2020) .

Considering the application of this second layer to a complete data-set, with new values for the epidemiological parameter each T 2 days, the output of this layer stands for the epidemiological time-series vectors denoted β opt , γ opt and α opt . These time-series are, then, used to fit the auto-regressive models in Eq. (7), in the third layer. We note that these SIRD parameter dynamics are the ones that can be can be used for forecasting and control purposes (and also to calculate the effective reproduction number R t ).

The third layer of the identification procedure resides in fitting an auto-regressive model to the timeseries derived from the optimization β opt , γ opt and ρ opt (note that ρ opt is given as a function of β opt , γ opt and α opt , as gave Eq. (10)). In this paper, the auto-regressive functions in Eq. (7) are of Auto-Regressive Integrated Moving Average (ARIMA) kind. The ARIMA approach is widely used for prediction of epidemic diseases, as shown by Kirbas et al. (2020) . It follows that, from a time-series viewpoint, the ARIMA model can express the evolving of a given variable (in this case, the epidemiological parameters) based on prior values. Such models, then, are coherent with the prequel discussion regarding the convergence of these parameters to steady-state conditions (Sec. 2.2).

For simplicity, instead of presenting the ARIMA fits for the three epidemiological time-series (β opt , γ opt and ρ opt ), we proceed by focusing on the SARS-CoV-2 transmission factor β. We note that equivalent steps are pursued for the other parameters. The main purpose of this third layer is to model the trends of the SIRD epidemiological parameters (as provided by the two previous layers) and use these trends, in the fashion of Eq. (7), to improve the forecasting of the SIRD+ARIMA model, making it more coherent and consistent for feedback control strategies.

It is worth mentioning that this layer is an innovative and important advantage of the SIRD model identification proposed in this work. As depicted by Fernndez-Villaverde & Jones (2020), using a SIRD model with time-varying epidemiological parameters allows one to provide forecasts that differ at the initial and long-term moments of the COVID-19 contagion.

As exploited in Section 2, the transmission rate parameter β gives an important measure to analyze the pandemic panorama. It has been shown that this parameter varies according to health measures applied to the prone population Fernndez-Villaverde & Jones (2020); Oliveira et al. (2020) ; Caccavo (2020) . The ARIMA expression is given as follows:

where δβ k = [δβ(k) δβ(k − 1) δβ(k − 2) . . . δβ(k − n β )] T is the vector of each incremental difference δβ(k) = β(k) − β(k − 1). Notice that, regarding the total regression order of this model, in accordance with Eq. (7) . Then, the ARIMA fit is performed minimizing an Extended Least-Squares (ELS) Procedure, used to find depolarized estimations for the ARIMA parameter values (a β0 , . . . , b βn β ). This procedure is based on the following steps:

(i) Concatenate the ARIMA the regression term from Eq. (18) as β(k) = ω T (k − 1)Θ + (k), where ω(k − 1) concatenates the input values from the previous layers (in this case, the time-series β opt ) and Θ concatenates the ARIMA coefficients;

(ii) Determine a constrained Least-Squares estimationΘ, w.r.t. Eq. (18);

(iii) Compute the Least-Squares residue = β opt − ωΘ and use it as the initial guess for for the next iteration;

(iv) Iterate until convergence or some (2-norm wise) small residue is achieved.

With respect to the detailed identification procedure, the datasets provided by Brazilian Ministry of Health are used for validation purposes, considering the interval from the first confirmed COVID-19 case in the country, dating February 26th, until the data from June 30th, 2020. The Ministry of Health disclosed daily values for the cumulative number of cases and deaths. The data is sufficient to obtain the behaviors for S, I, R, and D. Note that the total population size of Brazil is used as the initial condition N 0 ≈ 210 million.

Regarding the second layer, the optimization (normalization) weights are presented in Table 2 . Moreover, the uncertainties bounds taken as δ = 0.95 and δ = 1.05. Parameter w 1 w 2 w 3 Value 1 10 2

We proceed the model validation in a twofold: a) first we consider the first 100 data points to tune the SIRD+ARIMA model and demonstrate its validity, globally well representing the following points; and b) we tune the SIRD+ARIMA model by identifying its parameters for the whole dataset; this is the model used for control. We note that the values used for the social isolation factor ψ are those as exhibited in Figure 4 .

Regarding the first validation part, Figure 7 shows the estimated time-series values from the optimization output. Recall that these parameters are piece-wise constant for periods of T 2 = 7 days. As discussed in previous literature, these parameters tend to follow stationary trends as the pandemic progresses. Figure 7 : β, γ, α and ρ parameters identified by the two layer approach -from February 26th, 2020 to June, 04th 2020.

In order to illustrate the ARIMA fits, the identified auto-regressive model for SARS-CoV-2 transmission parameter β is provided in Figure 8 . The ELS procedure yields a regression f β (·) with n β = 21 daily steps (i.e. 3 weeks). Figure 8 demonstrates how the ARIMA is indeed able to describe the time-series β opt , globally catching the time-varying behavior with 99.04 % accuracy. The first 100 data points comprise the period from February 26th to June, 04th 2020. Using the resulting SIRD+ARIMA model identified in this period, we proceed by extending the model forward, making forecasts from June, 05th until June 30th, 2020, in order to demonstrate the validity of the identification. The forecasts are made using Eqs. (2) wherein parameters γ, β and ρ are given by the ARIMA regressions model in Eq. (7). Figures 9, 10 and 11 show the model forecasts compared with real data for active infections, recovered individuals and deaths, respectively, considering a interval of confidence of 95%. Clearly, the global behavior of the COVID-19 contagion is well described by the proposed models. To further illustrate this fact, Figure 12 depicts the model-data error terms (from Eqs. (14)-(16)), given in percentage. The coefficient of determination for these identification results are all above 0.99.

We note that the data-driven model, that is used for the NMPC strategy, is the based on the whole set of data. The simulation of this model is given in Figure 13 .

Based on the prior developments, given the high coefficient of determination and observing Figures 8-12 (in which a 95% confidence interval is considered), it can be noticed that the identified model SIRD+ARIMA describes very well the observed COVID-19 pandemic behavior in Brazil and, in addition, it can forecast the SIRD curve with relatively small error. It is worth mentioning that this model identification approach needs to be updated regularly. From Figure 12 , one notes that 20-days-ahead forecasts can be derived with model-data errors below 10%. This error margins may increase when considering larger forecast horizons. Nevertheless, regarding the control purpose of this work and taking into account prediction horizons of roughly one month, the proposed SIRD+ARIMA is indeed able to be highly representative regarding the COVID-19 spread. For perfectly coherent forecasts, it is recommendable for the identification to be reperformed each week, and the ARIMA fits updated. This recommendation has also been provided by Morato et al. (2020) .

We also remark that the ARIMA fits for the epidemiological parameters are given in weekly samples i.e. k week , while the SIRD variables are given for each day i.e. k day . Regarding this matter, we can represent the weekly-sampled variables, from the viewpoint of the daily-sampled variables, as:

where T 2 = 7 days. Finally, considering this validated SIRD+ARIMA model, the basic reproduction number R t of the SARS-CoV-2 virus in Brazil can be inferred through Eq. (3). In Figure 14 , the evolution of the viral reproduction number R t for each week of the pandemic in show. As it can be seen, this reproduction number presents stronger variations at the beginning stage of, but then tends to converge to a steadier values, which corroborates with the model extensions considering the epidemiological parameters. Moreover, notice that the reproduction number for since June 4th, 2020 is somewhat steady near 1.603, which indicates that the COVID-19 pandemic is still spreading in the country. 

In this Section, we detail the proposed NMPC strategy used to determine public health guidelines (regarding social isolation policies) to mitigate the spread of the COVID-19 contagion in Brazil. This NMPC algorithm is conceived under the feedback structure illustrated in Figure 5 .

We note that the generated control action u(k) represents the input to the population's response to social distancing measures, as gives the differential Eq. (4). For this reason, it follows that the proposed NMPC algorithm operates with a weekly sampling period (T 2 ). Implementing a new social distancing guideline every week is coherent with previous discussions in the literature Morato et al. (2020) . It does not seem reasonable to change the distancing measures every few days.

Before presenting the actual implementation of NMPC tool, we note that its generated control sequence must be given in accordance with the constraints expressed in Eq. (6).

Considering that the NMPC has a horizon of N p steps (given in weekly samples) from the viewpoint of each week k, the generated control signal is 2 :

Since the variations from each quarantine guideline u(k − 1) to the following u(k), denoted δu(k) are equal to ± 0.05 or 0, it follows that all possible control sequences can be described as, from the viewpoint of sample k:

From Eq. (21), we can conclude that, in the considered settling, there are 3 Np possible control sequences.

The main purpose of social isolation is to distribute the demands for hospital bed over time, so that all infected individuals can be treated. It seems reasonable to act though social isolation to minimize the number of active infections (I), while ensuring that this class of individuals remains smaller than the total number of available ICU beds n ICU .

Moreover, it seems reasonable to ensure that social isolation measures are enacted for as little time as possible, to mitigate the inherent economic backlashes.

This trade-off objective (mitigating I against relaxing social isolation) is expressed mathematically as:

where n I is a nominal limit for I (i.e. the initial population size N 0 ) included for a magnitude normalization of the trade-off objective. Note that u is given within [0.3 0.7], so there is no need for normalization.

The proposed NMPC algorithm must act to address the Control Objective given in Eq. (22) while ensuring some inherent process constraints. These are:

1. Ensure that the control signal has the piece-wise constant characteristic, as gives Eq. (6); and 2. Ensuring that the number of infected people with active acute symptoms does not surpass the total number of available ICU hospital beds in Brazil, this is:

Bearing in mind the previous discussions, the NMPC procedure is formalized under the following optimization problem, which is set to mitigate the SARS-CoV-2 viral contagion spread with respect to the trade-off control objective of Eq. (22) and the constraints from Eqs. (6)-(23), for a control horizon of N p (weekly) steps, from the viewpoint of each sampling instant k:

δu(k + i − 1) = −0.05 or 0 or 0.05 ,

The finitely parametrized NMPC methodology has been elaborated by Alamir (2006) . This paradigm has recently been extended to multiple applications Rathai et al. (2018 Rathai et al. ( , 2019 , with successive real-time results.

Given that this control framework offers the tool to formulate (finitely parametrized) social distancing guidelines for the COVID-19 spread in Brazil, with regard to the discussions in the prior, we proceed by detailing how it is implemented at each sampling instant k, ensuring that the Nonlinear Optimization Problem in Eq. (24) is solved.

Basically, the parametrized NMPC algorithm is implemented by simulating the SIRD+ARIMA model with an explicit nonlinear solver, testing it according to all possible control sequences (as gives Eq. (21)). Then, the predicted variables are used to evaluate the cost function J(·). The control sequences that imply in the violation of constraints are neglected. Then, the resulting control sequence is the one that yields the minimal J(·), while abiding to the constraints. Finally, the first control signal is applied and the horizon slides forward. This paradigm is explained in the Algorithm below. Figure 15 illustrates the flow of the implementation of the proposed Social Distancing control methodology for the COVID-19 viral spread in Brazil.

Initialize: N (0), S(0), I(0), R(0) and D(0). Require: Q, n I , n ICU Loop every T 1 days:

• Step (i ): "Measure" the available contagion data (N (k), S(k), I(k), R(k) and D(k));

• Step (ii ): Loop every T 2 days:

-Step (a): For each sequence j ∈ 3 Np : * Step (1): Build the control sequence vector according to Eq. (21); *

Step (2): Explicitly simulate the future sequence of the SIRD variables; *

Step (3): Evaluate if constraints are respected. If they are not, end, else, compute the cost function J(·) value.

-end -Step (b): Choose the control sequence U k that corresponds to the smallest J(·).

-Step (c): Increment k, i.e. k ← k + 1.

• end

• Step (iii ): Apply the local control policy u(k) as gives Eq. (6);

• Step (iv ): Simulate the SIRD+ARIMA model, as give Eqs.

(2)-(4).

•

Step (v ): Increment k, i.e. k ← k + 1.

end Figure 15 : Algorithm Implementation Flowchart.

Bearing in mind the previous discussions, the proposed SIRD+ARIMA model and the finitely parametrized NMPC toolkit, this Section is devoted in presenting the simulation results regarding the COVID-19 contagion spread in Brazil. The following results were obtained using Matlab. The average computation time needed to solve the solution of the proposed NMPC procedure is of 2 ms 3 .

We proceed by depicting two simulation scenarios, which account for the following situations:

(a) What would be the case if the NMPC generated social distancing policy was enacted back in the beginning of the COVID-19 pandemic in Brazil.

(b) How plausible is it to apply the NMPC technique from now on (30/06/2020) and reduce and mitigate the effected of the spread.

In order to provide these scenarios, we consider an "open-loop" comparison condition, which represents the simulation of the SIRD+ARIMA model with the social distancing factor in fact observed in Brazil (see Figure 4 ). The results are based on 126 days of data. From the 127 th sample on, the open-loop simulation takes ψ = 0.3 (no isolation).

The proposed NMPC is designed with a cost function J(·) with tuning weight Q = 0.7 (which means reducing the number of infections is prioritized over relaxing social distancing). Furthermore, the NMPC optimization is set with a prediction horizon of N p = 4 weeks (28 days), which means that the controller makes its decision according to model-based forecasts of the SIRD+ARIMA model for roughly one month ahead of each sample k. For simplicity n I is taken as n ICU psym , since the main goal of the MPC is to ensure p sym I ≤ n ICU .

We recall that the COVID-19 pandemic in Brazil formally "started" 26/02/2020, when the first case was officially registered.

This first simulation scenario considers the application of the NMPC control strategy to guide social distancing from 30 days after the first case (27/03/2020). We choose this date because, in Brazil, the majority of states started social distancing measures (in different levels) around the period of late-March/mid-April. We do not deem it reasonable to consider the application of the NMPC strategy from "day 1", since this was not seen anywhere in Brazil.

Considering this control paradigm, Figure 16 depicts the predicted evolution of severe/acute symptomatic COVID-19 infections over time. These active symptomatic infections stand for those that may need ICU hospitalization. The results shows that the NMPC is able to thorough attenuate the peak of infections, ensuring that it always stays below the ICU threshold. This is a quite significant results, which shows that the proposed feedback framework is able to offer an enhanced paradigm, with time-varying social distancing measures, such that the COVID-19 spread curve is indeed flattened, never posing serious catastrophic difficulties to Brazilian hospitals (and health system overall). The result also indicated that the social distancing should vary with relaxing and strengthening periods over the years of 2020, 2021 and 2022. We note, once again, that health professionals should better qualify the relationship between the social distancing factor and actual economic/social restrictions, as illustrated in Table 1 . The results also indicate that if no coordinated action is deployed by the Brazilian federal government, the number of active symptomatic infections at a given day could be up to 540000 individuals, with this peak forecast to October 16th, 2020. If the NMPC strategy was indeed applied, two peaks would have been seen, with 80 % of total ICU capacity, previewed for September 25th, 2020 and March, 31rst, 2021.

Regarding this scenario, Figure 17 shows the evolution of the total amount of infections and cumulative number of deaths. This Figure also places the real data points against the simulated model. The possibilities to come are catastrophic, as also previewed by Morato et al. (2020) . We note, once again, that the SIRD+ARIMA models offers qualitative results, the magnitude of 1.5 million deaths seems quite alarming, but the model is identified considering a large margin for sub-reports (in number of deaths and confirmed cases). Not all deaths due to COVID-19 are currently being accounted for in Brazil, as discuss The Lancet (2020); Zacchi & Morato (2020) . 

The second simulation scenario that we consider is to control the situation from the current moment (30/06/2020) to avoid a total collapse of the Brazilian health system. We note that as of this date, the country counts over 1.4 million confirmed COVID-19 cases and more than 59500 deaths due to the SARS-CoV-2 virus. Figure 18 shows the main results regarding this second scenario, considering the active symptomatic infections, that may require ICU hospitalization. Even though a partial catastrophe is already under course in Brazil, if social distancing are guided through the proposed NMPC, a total collapse of the public health system can still be avoided. Figure 19 shows the evolution for recovered individuals and cumulative number of cases I c , Figure 20 shows the evolution of all active infections (symptomatic and asymptomatic), while Figure 21 shows the mounting number of deaths. The proposed NMPC, if rapidly put in practice, could still be able to slow the viral spread, saving 25 % of lives w.r.t. the open-loop condition. The peak of infections, if such technique is applied, has its forecast previewed to September 2nd, 2020, being anticipated in 17 days. 

In this paper, an optimal control procedure is proposed for ruling social isolation guidelines in Brazil, in order to mitigate the spread of the SARS-CoV-2 virus. To do so, a new model is proposed, based on extensions of the SIRD equations. The proposed model embeds weekly auto-regressive dynamics for the epidemiological parameters and also takes a dynamic social distancing factor (as seen in previous papers ). The social distancing factor measures and expresses the population's response to quarantine measures, guided by the Nonlinear Model Predictive Control procedure. The NMPC strategy is designed within a finitely parametrized input paradigm, which enables its fast realization.

In this work, some key insights were given regarding the future panoramas for the COVID-19 pandemic in Brazil. Below, some of the main findings of this paper are summarized:

• The presented results corroborate the hypothesis formulated in many of the previous papers regarding the COVID-19 pandemic in Brazil Hellewell et al. (2020) ; The Lancet (2020); Morato et al. (2020) : herd immunity is not a plausible option for the country 4 ; if no coordinates social distancing action is enforced, the ICU threshold will be largely surpassed, which can lead to elevated fatality.

• The prediction of evolution of the viral spread is relatively accurate with the proposed adapted SIRD model for up to 20 days. Larger prediction horizons can be considered, but daily model-updates are recommended.

• The simulation forecasts derived with the NMPC strategy and with an open-loop condition (no social distancing) indicate that social distancing measures should still be maintained for a long time. The strength of these measures will be diluted as time progresses. The forecasts indicate an infection peak of over 600000 symptomatic individuals to late September, 2020, in the current setting. If model-based control is enacted, the peak could be anticipated and the level of infections could be contained below the ICU hospital bed threshold. The NMPC could save over 400000 lives if enacted from now (July, 2020).

• The results also indicate that if such coordinated control strategy was applied since the first month of COVID-19 infections in Brazil, a more relaxed social distancing paradigm would be possible as of late 2020. Since this has not been pursued, the social distancing measures may go up until late 2021 if no vaccine is made available.

These results presented in this paper are qualitative. Brazil has not been testing enough its population (neither via mass testing or sampled testing), which means that the data regarding the number of infections is very inconsistent. As thoroughly details, the uncertainty margin associated to the available data (in terms of case sub-reporting) is very significant. Anyhow, the results presented herein can help guiding long-term regulatory decision policies in Brazil regarding COVID-19. One must note that social distancing measures, in different levels, will be recurrent and ongoing for a long time. Due to this fact, compensatory social aid policies should also be developed in order to reduce the effects of a possibly long-lasting economic turn-down. Recent papers have discussed this matter Ahamed & Gutiérrez-Romero (2020); Zacchi & Morato (2020) .

The authors report no financial disclosure nor any potential conflict of interests.

The simulations driving the worlds response to covid-19. how epidemiologists rushed to model the coronavirus pandemic?

The role of financial inclusion in reducing poverty: A covid-19 policy response

Stabilization of nonlinear systems using receding-horizon control schemes: a parametrized approach for fast systems volume

COVID-19: from model prediction to model predictive control

Modeling and forecasting the early evolution of the COVID-19 pandemic in Brazil

The COVID-19 (SARS-CoV-2) uncertainty tripod in Brazil: Assessments on model-based predictions with large under-reporting

Social distancing and movement constraint as the most likely factors for COVID-19 outbreak control in brazil

A new twenty-first century science for effective epidemic response

Short-term forecasts of COVID-19 deaths in multiple countries

Chinese and italian covid-19 outbreaks can be correctly described by a modified sird model. medRxiv

Model predictive control

Tracking the onset date of the community spread of SARS-CoV-2 in western countries. medRxiv

Demographic science aids in understanding the spread and fatality rates of COVID-19

Estimating and Simulating a SIRD Model of COVID-19 for Many Countries, States, and Cities. NBER Working Papers

Protect indigenous peoples from COVID-19

Temporal dynamics in viral shedding and transmissibility of COVID-19

Feasibility of controlling COVID-19 outbreaks by isolation of cases and contacts. The Lancet Global Health

Social isolation map covid-19

Modeling Infectious Diseases in Humans and Animals

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

A contribution to the mathematical theory of epidemics

Comparative analysis and forecasting of covid-19 cases in various european countries with arima, narnn and lstm approaches

Robust and optimal predictive control of the COVID-19 outbreak

Early dynamics of transmission and control of COVID-19: a mathematical modelling study. The Lancet infectious diseases

The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: estimation and application

Information about the new coronavirus disease (COVID-19)

Coronavirus pandemic (covid-19)

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

Mathematical modeling of COVID-19 transmission dynamics with a case study of wuhan

Evaluating the burden of covid-19 on hospital resources in bahia, brazil: A modelling-based analysis of 14.8 million individuals. medRxiv

Estimation of COVID-19 under-reporting in brazilian states through SARI

Epidemic analysis of COVID-19 in China by dynamical modeling

A parameterized NMPC scheme for embedded control of semi-active suspension system

GPU-based parameterized nmpc scheme for control of half car vehicle with semi-active suspension system

Expected impact of COVID-19 outbreak in a major metropolitan area in Brazil

COVID-19 in latin america: the implications of the first confirmed case in Brazil

Modeling and forecasting the covid-19 pandemic in india

A bayesian analysis of the total number of cases of the COVID-19 when only a few data is available. a case study in the state of Goias

Understanding of COVID-19 based on current evidence

Modeling covid-19 epidemic in heilongjiang province, china

Covid-19 in Brazil: so what

The COVID-19 pandemic in Brazil: chronicle of a health crisis foretold

The COVID-19 pandemic in Brazil: An urge for coordinated public health policies

The Authors acknowledge the financial support of National Council for Scientific and Technological Development