A parametrized nonlinear predictive control strategy for relaxing COVID-19 social distancing measures in Brazil.

The SARS-CoV-2 virus was first registered in Brazil by the end of February 2020. Since then, the country counts over 150000 deaths due to COVID-19 and faces a profound social and economic crisis; there is also an ongoing health catastrophe, with the majority of hospital beds in many Brazilian cities currently occupied with COVID-19 patients. Thus, a Nonlinear Model Predictive Control (NMPC) scheme used to plan appropriate social distancing measures (and relaxations) in order to mitigate the effects of this pandemic is formulated in this paper. The strategy is designed upon an adapted data-driven Susceptible-Infected-Recovered-Deceased (SIRD) model, which includes time-varying auto-regressive immunological parameters. A novel identification procedure is proposed, composed of analytical regressions, Least-Squares optimization and auto-regressive model fits. The adapted SIRD model is validated with real data and able to adequately represent the contagion curves over large forecast horizons. The NMPC strategy is designed to generate piecewise constant quarantine guidelines, which can be reassessed (relaxed/strengthened) each week. Simulation results show that the proposed NMPC technique is able to mitigate the number of infections and progressively loosen social distancing measures. With respect to a "no-control" condition, the number of deaths could be reduced in up to 30% if the proposed NMPC coordinated health policy measures are enacted.

Introduction

The SARS-CoV-2 virus causes a severe acute respiratory syndrome, which can become potentially fatal. This virus has spread rapidly and efficiently, becoming a worldwide pandemic. By mid-June , SARS-CoV-2 had already infected over million people. As of mid-August, over 21.5 million COVID-19 cases had been confirmed. In order to illustrate the severity of the COVID-19 contagion, Fig. 1 displays its evolution in different countries and a world map of concentration levels.

Fig. 1

COVID-19 in the world. Data published by the European Center for Disease Prevention and Control (ECDC) [1].

Unfortunately, vaccines for this virus are not expected to be ready until mid-. Therefore, in order to control this disease, most countries have adopted social distancing measures since March (in different levels and with diverse strategies) [2]. This tactic has been explicitly pointed out as the most pertinent option regarding the COVID-19 outbreak, even in the case of countries with great social inequalities, such as Brazil [3]. The concept underneath social distancing is to prevent the saturation of health systems due to large amounts of COVID-19 infections that would require treatment at the same time. Thus, when social distancing policies are enacted, the demands for treatment become diluted over time, and health systems do not have to deal with hospital bed shortages associated with a concentrated peak of infections.

Brazil has shown itself to be a very particular case regarding COVID-19, as discuss Croda et al. [4]: the country is very large, with federated states that have chosen independent tactics to address the contagion. Furthermore, there has been no coordinated national public health policy; the federal government has openly expressed its unwillingness to determine any strict quarantine measure, claiming that the negative economic effects would be too steep [5], [6], [7]. The country currently ranks as second w.r.t. COVID-19 cases and deaths. The expectations disclosed on the recent literature suggest catastrophic scenarios for the coming months, which may last until mid-, or even later [8], [9], [10].

Brazil has displayed near-collapse situations in many hospitals, since the majority of Intense Care Unit (ICU) beds are currently occupied with COVID-19 patients, all around the country.1 Furthermore, the scientific community has been warning that the true amount of infected individuals in the country is very under-reported w.r.t. to the data disclosed by the Ministry of Health [12], [13]; due to the absence of mass testing, only patients with moderate to severe symptoms are being accounted for. As of October , , the country counts over 5.1 million confirmed COVID-19 cases and deaths. The overall situation is very critical.

Recent literature has demonstrated how the SARS-CoV-2 viral contagion dynamics can be appropriately described by Susceptible–Infected–Recovered–Deceased (SIRD) models, as discussed by Wang et al. [14] and Kucharski et al. [15]. Therefore, regarding the previous discussion, we consider the Brazilian COVID-19 contagion, modeled through SIRD equations which incorporate a social distancing factor as an additional input variable. Regarding this context, we investigate how to perform optimal coordinated social distancing interventions to address the ongoing pandemic scenario in Brazil. Well-designed social distancing guidelines, with consistent duration and implementation, help mitigate the contagion spread which, thus, avoids the saturation of heath systems, while, at the same time, balances social and economic side effects, by releasing/relaxing the quarantine measures as soon as possible.

Since SIRD models comprise coupled nonlinear differential equations, the Nonlinear Model Predictive Control (NMPC) [16] framework is a rather convenient approach to guide optimal model-based public health policies for pandemics, given that it can adequately consider the nonlinear dynamic of the virus spread together with the effect of lockdown/quarantine measures, which are introduced as constraints of the optimization problem. Indeed, the optimal control of the COVID-19 outbreak has become a very active research area, since global scientific efforts have been addressing this topic [2]. This paper lies within the context of this evolving field. Below, we highlight some of correlated studies:

Djidjou-Demasse et al. [17] consider moving-window optimal control for COVID-19, taking into account piecewise constant parameter estimation. The parameter variation framework provides interesting results in terms of prediction fidelity;

Kantner [18] also formulates an optimal moving-window, but integrates some notions of uncertainty regarding the model, which may be corrupted due to under-reported data;

Köhler et al. [19] apply the MPC technique for the German COVID-19 context. The approach considers that the control input directly affects the infection and transmission rates (parameters of the SIRD model);

Alleman et al. [20] consider the application of a NMPC algorithm regarding the data from Belgium. The control input is the actual isolation parameter which attenuates the transmission factor in the SIRD model;

Morato et al. [10] 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 is modeled with an additional dynamic variable, which is incorporated to the SIRD dynamics.

Building upon these previous results, our main contribution is a novel NMPC algorithm synthesized upon an adapted SIRD model, which is identified and validated for the Brazilian context. Following the lines of Morato et al. [10], we assume that the NMPC action is a social isolation guideline, which is enacted and passed on to the population. Consequently, 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. [19] and Bastos et al. [13], we upgrade the SIRD model with an additional dynamic nonlinear variable, which gives the input/output (I/O) relationship between the COVID-19 contagion infection and transmission rates according to the stage of the pandemic. The SIRD model parameters (transmission factor, infection factor, lethality rate) are weekly-varying, given through auto-regressive equations. These regressions are used to provide better long-term COVID-19 forecasts than the regular SIRD equations.

The main novelty of this paper is a twofold: (a) the development of an adapted SIRD model, which offers consistent long-term forecasts, and (b) the NMPC algorithm formulation through finitely parametrized control inputs, which makes its implementation persuadable through search mechanisms, which run much faster than the Nonlinear Programming methods from the previously referenced works. Furthermore, the finitely parametrized control input values can be “translated” into practical social distancing guidelines, such as closing schools, obliging the use of masks in public areas, closing public transport and other measures. We highlight the main ingredients of our contribution:

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. The adapted model also incorporates the dynamic response of the population to a given social distancing guidelines. This model is identified through a novel identification procedure, based on analytical expansions, Least-Squares optimization and auto-regressive fits. The method is able to offer long-term forecasts with relatively small error and parameters that are biologically coherent.

Then, a Parametrized Nonlinear Predictive Control algorithm is proposed. 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 the prediction horizon. The NMPC solution is found through a simple search mechanism regarding these simulated sequences, which is rather numerically cost-efficient.

Finally, we present consistent simulation results considering the application of this NMPC algorithm to the Brazilian COVID-19 scenario. These results are a twofold: (a) those that regard the application of the optimal strategy to control the pandemic since its beginning, with comparisons to what, in fact, happened in Brazil; and (b) the application of the method from July 30th onwards, aiming to mitigate and revert the current health crisis.

This paper is structured as follows. Section 2 presents the new SIRD model with auto-regressive epidemiological parameters, its 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.

Model, identification and validation

The COVID-19 dynamic model is detailed in this Section. In this paper, we use an adapted SIRD model, building upon from works of Wang et al. [14], Kucharski et al. [15] and Ndairou et al. [21]. The “SIRD” (Susceptible–Infected–Recovered–Dead) model is adapted in order to incorporate the social distancing factor and piecewise time-varying epidemiological parameters, which vary according to the stage of the pandemic. We note that the use of time-varying epidemiological parameters is in accordance with recent immunology discussion and theoretical results, e.g. [22], [23], [24]; these papers have highlighted how the transmission and reproduction rates of the SARS-CoV-2 virus indeed vary over time.

SIRD Epidemiological model

The SIRD model describes the spread of a given disease w.r.t. a population that is split into four non-intersecting classes:

The total amount of susceptible individuals, which are prone to contract the disease at a given sample of time , denoted ;

the individuals that are currently infected with the disease (active infections at a given sample of time ), denoted ;

the total amount of recovered individuals, that have already recovered from the disease, from an initial instant until the current sample , denoted ;

and, finally, the total amount of deceased individuals, from an initial instant until the current sample , denoted .

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

Since the Brazilian government discloses daily samples of total confirmed infections and deaths due to COVID-19, we consider that the SIRD discrete-time dynamics samples are given each day. Furthermore, in order to account for the average incubation period of the virus (of, in average, 5.1 days, see e.g. [25]), we consider that the epidemiological parameters vary weekly (each days).

There are three major epidemiological parameters which are accounted for in the SIRD model. These parameters express the specific characteristics of the modeled contagion within a given population set:

The transmission rate parameter , which represents the average number of contacts that are sufficient for transmission of the virus from one individual to another. It follows that gives the new contagion infections at a given day .

The infectiousness Poisson parameter , which stands for the inverse of the period of time for which a given infected individual is indeed infectious. This parameter directly quantifies the amount of individuals that “leaves” the infected class.

The mortality rate parameter , which gives the ratio of infected individuals that die. The amount of deceased individuals due to the contagion, at a given moment , is given by .

The SIRD model is expressed through the following nonlinear discrete-time difference equations: where denotes the portion of the infected individuals which in fact display symptoms. The cumulative number of cases stands for the total number of people that have been infected by the SARS-CoV-2 virus since the beginning of the pandemic. This variable is analytically equivalent to the sum of the individuals that are currently infected with those that have already recovered and those that have deceased .

The class of active infections is split into two compartments: active symptomatic infections and active asymptomatic infections , as accounted for in [10], [26]. The “asymptomatic class” encompasses those individuals that are asymptomatic or lightly-symptomatic, that may transmit the virus but will not die or require hospitalization. Accordingly, the “symptomatic class” denotes the infections that will, in fact, require hospital treatment. We consider, for simplicity, that the symptomatic compartment is given by , whereas the asymptomatic remainder is , being the ratio between those that develop acute symptoms over the whole class of infections; this symptomatic ratio parameter is constant and borrowed from [10], [26]. We note that the class of recovered individuals stands for all individuals who became immunized against SARS-CoV-2, encompassing those that displayed symptoms and those that did not.

We emphasize that represents a transmission rate mitigation factor. This factor expresses the observed social isolation ratio within the susceptible population set . It follows that stands for a situation in which the whole population set has sustained social interactions. As discussed in previous papers [13], [27], there exists a “natural” factor, which stands for normal conditions, still with “no control” of the viral spread (with no social isolation guidelines or measures2 ). In contrast, represents a complete lockdown conditions, for which there are no social interactions (this is, in practice, unattainable). In this paper, we consider the bounds for the social isolation variable from the work of Bastos et al. [13].

Another essential epidemiological property of a viral spread is its basic reproduction number, usually denoted . This parameter quantifies the average potential of viral transmission and represents how many COVID-19 cases are expected to be generated due to a single primary case, within a population for which all individuals are susceptible. Theoretically, for stabilized viral contagions (such as common flu), this parameter is constant and inherent to the disease. However, in the case of spreading contagions, the number of susceptible individuals changes over time. Thus, one is only able to quantify the effective reproduction number of the contagion, denoted . From a control viewpoint, represents the epidemic spread velocity: if , it follows that the infection is spreading and the number of infected people increases along time (this typically happens at the beginning of the epidemic), otherwise, if , it means that there are more individuals “leaving” the infected class, either recovering or dying, than “entering” (new infections), which means that the epidemic is ceasing. We note that the viral effective reproduction number is affected by a series factors, including the immunology of the virus, biological characteristics and health policies. Since the epidemic parameters change over time, we take the effective reproduction number as a dynamic variable, i.e. .

The underlying assumption used to calculate the viral effective reproduction parameter, is that, at the beginning of the pandemic, . Considering parameters , , and from the SIRD model Eq. (2), is given by:

The SIRD model given in Eq. (2) accounts for as the initial population size (prior to the viral infection). We stress that represents the active infections at a given moment, while and represent the total amount of deaths and recovered individuals until this given moment, respectively. For this reason, it follows that and are proportionally dependent to .

We note that we are not able to use more “complex” descriptions of the COVID-19 contagion for the Brazilian context, such as the “SIDARTHE” model3 used by Köhler et al. [19] because we have insufficient amount of data. The Ministry of Health only discloses the total amount of infections () and the total amount of deaths (), per 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. [19] originate. If we considered more complex models, the truthfulness/validity of the identification and simulation results may be corrupted, since the identified parameters may represent a singular combination that matches the identification datasets, but that cannot be used for forecasting/prediction purposes.

SIRD model solution

Before further development, we must note that the SIRD model, as given by Eq. (2), has an inherent positivity property. All variables and parameters involved in these equations are positive defined. Furthermore, as discuss and thoroughly investigated by Harko et al. [28] and Bohner et al. [29], the exact analytical solutions to SIRD equations (with constant epidemiological parameters) are time-decaying exponentials, integrative exponentials and logarithmic curves. In a “no-control” condition, with constant parameters and constant social distancing measures, it follows that:

The susceptible curve has an exponentially convergent decay pattern, moving from an initial condition in the direction of an asymptotic equilibrium ;

the active infections curve has an exponential derivative pattern, with an increase, a peak and a decrease (as the contagion ceases); nominally, it follows that ;

the recovered and deceased curves and show exponential convergent increases, departing from the origin to the direction of asymptotically stable equilibria .

The SIRD model in Eq. (2) displays the positiveness property, which means that all variables are non-null and greater than , .

Assume an initial condition for which all individual are susceptible and in constant contact, i.e.  and . Assume, also, that there is a non-null initial number of active infections, i.e.  and that initially there are none deceased or recovered, i.e. .

Note that all virus-related parameters are positive-defined, due to immunological implications, i.e. , and , . Furthermore, we note that is the sampling period (given in number of days). The transmission parameter is necessarily contained within . The Poisson parameter gives the inverse of the average infection period (in days) of the virus, which implies that (at least two days of infection). The lethality parameter is always positive and non-null; it can also be considered as upper-bounded by 0.1, regarding the SARS-CoV-2 case. Lastly, we stress that and, thus, .

We note that the dynamics of , , and are necessarily positive if displays the positivity property, since they are proportional or proportional–integrative to . Thus, it remains to demonstrate the positiveness of and .

Regarding the susceptible curve, we note that, departing from an initial condition , it follows that for positive values of . We can re-arrange the susceptible dynamics as follows: Since , we have that and, thus, .

It only remains to demonstrate that is positive. We can re-arrange the active infection dynamics as follows: Now, we can evaluate the bounds on the time-varying parameter . Based on the immunology properties of the parameters, we have that . Furthermore, since and , we can conclude that . In this case, for any positive initial condition , it follows that . This concludes the proof.  □

In order to illustrate the dynamics of the COVID-19 pandemic outbreak through a SIRD model, Fig. 2 presents forecasts for the COVID-19 contagion spread in Brazil, considering constant epidemiological parameters and following the methodology presented by Bastos et al. [13]. These predictions, computed on June , , are here shown to qualitatively demonstrate the behaviors of and . The active infections curve shows an increase-peak-decrease characteristic, while the total number of deaths shows an asymptotic behavior to a steady-state condition. We note that, as of June , the number of infections had already significantly increased, which means that the most recent forecasts preview even worse scenarios. For the full exact analytic expressions for SIRD model, we invite the Reader to refer to [29].

Fig. 2

Long-term SIRD Forecasts for Brazil (from June , ). Consolidated datasets ( and ): (a) No-Control Situation (), (b) Hard social distancing () from June onwards (), (c) Hard social distancing applied from March () and (d) Real data (). Shades represent total variation over a 95% confidence interval, solid lines represent mean values.

Model extensions

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

In accordance with the discussions presented by Morato et al. [10], we understand that it is essential to incorporate the population’s response to public health policies into the SIRD model. This is, to take into account the dynamic effects of social distancing. When a government enacts a given social distancing policy, the population takes some time to adapt to it and to, in fact, exhibit the expected social distancing factor . Therefore, we take as a time-varying dynamic map of the enacted social distancing guidelines (which will be, later on, defined by an optimal controller).

Many works seen in the recent literature [3], [21], [26] simply opt to consider constant values for , and , acknowledging that these factors are inherent to the disease. Anyhow, recent results [13], [15] indicate that these parameters can be considered as time-varying, achieving steady-states at ending stages of the pandemic. In order to illustrate these time-varying characteristics of the transmission rate parameter , for instance, it is reasonable to consider that, when more active infections occur, the virus tends to spread more efficiently. Likewise, the mortality rate should exhibit an initial peak, when deaths occur for a smaller number of confirmed cases, stabilizing after the infection curve has decreased. With the aim to incorporate this issue to the SIRD model, auto-regressive moving-average time-varying dynamics for the viral parameters are proposed, which stabilize according to evolution 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 this incubation period.

If parameters are held constant for forecasting purposes, only qualitative predictions are offered, allowing short-term conclusions. Since dynamic models for the epidemiological parameters are proposed, more coherent long-term forecasts can be derived. We must stress 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.

Social distancing guidelines/the control input

In order to design and synthesize effective control strategies for social distancing policies, that are to be passed as guidelines to the population by the local government, 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 piecewise constant samples. Regarding the application of the control procedure (Section 4), this isolation factor will vary according to the enacted social distancing policy , as defined by a nonlinear optimal predictive control algorithm.

The differential equation that models the response of the susceptible population w.r.t. quarantine guidelines is taken as suggests [10], this is: being the actual control input: a 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. Furthermore, is a settling time parameter, which is related to the average time the population takes to respond to the enacted social isolation measures, and is a time-varying stating gain relationship between and .

As recommend Morato et al. [10], we assume that when more 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 . This is mathematically expressed as: being the total number of ICU beds in the country. We recall that is a parameter which gives the amount of infected individuals that in fact display symptoms and may need to be hospitalized (). We stress that is typically greater than in situations when the intense care units are saturating, with over 70% occupancy rate, which means that the populations responds to the social distancing guidelines more intensely.

Brazil has 45,848 ICU beds available. This number is estimated to increase in up to with field hospitals that were built specifically for the COVID-19 contagion (i.e. ). The percentage of symptomatic individuals is taken as of the total amount of infections, according to the suggestions of Bastos and Cajueiro [26]. This value is coherent with the available information regarding this virus, concerning multiple countries [30]. The percentage of symptomatic individuals corresponds to those with severe/acute symptoms, which will most possibly require treatment.

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

We proceed by considering that is a finitely parametrized control input. This is: the enacted social distancing guidelines can only be given within a set of predefined values. This approach is coherent with possible ways to enforce and put in practice social distancing measures: the set of predefined values for can be translated into a set of enforcement acts, such as closing public transports, closing shops, enforcing the use of masks, etc.

Pursuing this matter, the control input is defined as a “smooth” and piecewise constant signal. The maximal possible increase/decrease on is defined w.r.t. the historical data for the observed social isolation factor in Brazil, estimated in previous works [31], [32] and presented in Fig. 3. Accordingly, neglecting the strict increase on seen by mid-March (which stands for a moment of many strict isolation measures enforced at the same time), we assume that may vary by 5% of isolation per week. This percentage is the historical average variation of social isolation seen in Brazil (from April to July ), when social distancing measures were progressively strengthened over a time.

Fig. 3

Social distancing factor in Brazil from February , to June ,  [31].

Moreover, is defined within the admissibility set . As gives Fig. 3, which shows the observed social isolation factor in Brazil, the “natural” isolation factor is of , 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 Fig. 3). Anyhow, coordinated “lockdown” measures were not forcefully enacted. For this matter, we consider , which would represent that the population can be restricted to, at most, a 70% reduction in the level of social interactions. As reported by Bastos et al. [13], it seems unreasonable to consider larger values for social isolation in the country, due to multiple reasons (hunger, social inequalities, labor needs, lack of financial aid from the government for people to stay home, etc.).

Mathematically, these constraints are expressed as:

The “natural” isolation factor, which is the lower bound on , can be seen in Fig. 3, which shows the observed social isolation factor in Brazil and its lower bound of .

The constraints given in Eq. (6) are, in fact, very interesting from an implementation viewpoint, since governments could “translate” the finitely parametrized values for into actual practicable measures, as those depicted in Table 1. We remark that this Table is presented for illustrative purposes only; epidemiologists and viral specialists are the ones who could formally discuss which measures should be put in practice to ensure the social distancing factor guideline given by . We stress that changes in these guidelines do not affect the conceptual essence of this work, since the proposed methodology is general and can be applied w.r.t. the epidemic reality of any location.

Table 1

Illustrative example of finitely parametrized social distancing measures.

Control signal/Social distancing guideline (u)	Implemented measures	Infection risk
u=0.3	No public health emergency. All economy sectors can return to their normal activities.	Controlled contagion
⋮	⋮
u=0.35	Low restriction levels. Use of masks to go outside. Public transport functioning. Limited opening of shops and small public spaces.	Low risk
⋮	⋮
u=0.4	Moderate restrictions. Use of masks to go outside. Closed public spaces. Restricted openings only.	Moderate risk
⋮	⋮
u=0.45	Very restrictive policies. Reduced public transport. Urge for people to stay home at all times. Very restrictive openings only.	High risk
⋮	⋮
u=0.7	Severe restrictive policies. No public transport. Urge for people to stay home at all times. Only basic services may open, with reduced capacities.	Very high risk

Dynamic epidemiological parameter models

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

The models given in Eqs. (7)–(9) are possibly delayed and auto-regressive. Anyhow, despite the parameters , and being time-varying, the model functions , and are constant. The order and number of regressions are found by means of an optimization procedure, which is further detailed in Section 2.4.

The complete COVID-19 model

The complete model used in this work to describe the COVID-19 contagion outbreak in Brazil is illustrated by the block diagram of Fig. 4. Note that COVID-19 epidemiological parameter models given through Eqs. (7)–(9) offer accurate long-term predictions, as discussed in the beginning of this Section. Thus, we denote henceforth the cascade of the population response dynamics (Eq. (4)) upon the SIRD equations with auto-regressive time-varying epidemiological (Eqs. (7)–(9)) as the “SIRD＋ ARIMA” model.

Fig. 4

The “SIRD＋ARIMA” model for COVID-19 in Brazil.

Identification procedure

Recent numerical algorithms have been applied to estimate the model parameters of the COVID-19 pandemic, e.g. [13], [32], [33]. The SIRD model offers three degrees-of-freedom at each sampling instant (i.e. , and ), which means that different instantaneous combinations of these parameters can yield the same values for , , and . 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 immunology characteristics of the SARS-CoV-2 virus. An indicative of badly adjusted SIRD model parameters is the effective reproduction number of the disease , which should naturally decrease as the pandemic ceases.

We proceed by depicting the proposed identification procedure, which is performed in three consecutive steps/layers. This procedure, illustrated in Fig. 5, comprises the following steps:

Fig. 5

Parameter identification procedure.

The first step resides in analytically solving the SIRD regressions from Eq. (2) for a fixed interval of data samples, given in number of 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 w.r.t. the SIRD model structure, aiming fine adjustments of the parameter values according to the predefined (biologically coherent) sets, so that the identified model provides adequate parameter estimates. The output of this second step stands for time series vectors regarding the SIRD parameters.

Finally, these time series are used to fit auto-regressive models in the form of Eq. (7), via Extended Least-Squares minimization.

By pursuing this three-layered procedure, the estimated parameter curves are given in accordance with biological conditions. The optimization layer embeds feasibility sets for the epidemiological parameters that are in accordance with previous results for SIRD model estimations for Brazil [10], [26].

Through the sequel, the nonlinear difference equations for and are modified in order to decouple the parameters related to the number of deaths and the total number of recovered individuals, regarding . This is: where . Therefore, instead of identifying , and , we pursue the identification of , and , and, then, is computed as follows:

First layer: Analytical solutions

The identification procedure starts by collecting the available datasets regarding the cumulative number of COVID-19 cases and deaths , inside a fixed interval of points .

The first layer computes “exact” analytical values for the epidemiological parameters, denoted , and , which are found according to the following discrete analytical expansions:

Second layer: Ordinary least squares optimization

Assuming that the used datasets might be corrupted by a number of issues, such as cases that are not reported on given day and simply accounted for on the following days, or sub-reported cases, as discussed by Bastos et al. [13], we adjust the parameter estimated from the first layer through an optimization layer, which is used in order to improve the reliability of the identification.

The optimization procedure minimizes a constrained Ordinary Least Squares problem, whose solution comprises the following vectors: , , and . These vectors concatenate the model-based outputs for the considered discrete interval. The available data from the same discrete interval is used once again. The decision variables for the optimization procedure are the degrees-of-freedom of the SIRD model, denoted as , , . 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 biologically coherent values for these decision variables. The quadratic model-data error terms used in the optimization layer are: being the model-based estimations. The complete optimization problem is formulated as follows: where and are uncertainty interval margins used to define the lower and upper bound of each decision variable of the optimization problem. , and are the initial conditions of the optimization problem and , and are positive weighting values (tuning parameters), used to normalize the magnitude order of the total optimization cost .

Considering the application of this second layer to the complete data-set, we note that different values for the epidemiological parameter are given each days. The output of this second layer stands for the epidemiological time series vectors denoted , and . 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 used for forecasting and control purposes, and also to calculate the effective reproduction number .

Third layer: ARIMA fits

The third layer of the identification procedure resides in fitting an auto-regressive model to the time series derived from the optimization , and . Note that is given as a function of , and , as in Eq. (12).

Then, the auto-regressive functions in Eq. (7) are of “Auto-Regressive Integrated Moving Average” (ARIMA) kind. The ARIMA framework is widely used for prediction of epidemic diseases, as shown by Kırbaş et al. [34]. 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 (Section 2.3).

For presentation simplicity, instead of presenting the ARIMA fits for the three epidemiological time series (, and ), 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), in order 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 Lalwani et al. [35], the use of SIRD model with time-varying epidemiological parameters allows one to provide forecasts coherent with the evolution of the COVID-19 contagion.

ARIMA fit for the viral transmission rate

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. The used ARIMA expression is given as follows:

The ARIMA fit is performed minimizing an Extended Least-Squares (ELS) Procedure, used to find unbiased estimations for the ARIMA parameter values (, , ). This procedure is based on the following steps:

Consider as a variable which comprises the model residue, initialized as null. Concatenate the ARIMA the regression term from Eq. (20) as , where concatenates the input values from the previous layers (in this case, the time series ) and concatenates the ARIMA coefficients;

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

Compute the Least-Squares residue and use it as the initial guess for for the next iteration;

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

Validation results

Regarding 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 , until the data from June , . Note that the total population size of Brazil is used as the initial condition million. The following results were obtained using Matlab software with Yalmip toolbox and fmincon solver. The solver operation is unconstrained. The ARIMA model fitting was derived with Matlab System Identification toolbox. All results were performed in a Intel i Macintosh.

Regarding the second layer, the optimization (normalization) weights are presented in Table 2. Moreover, the uncertainties bounds taken as and .

Table 2

Optimization weights.

Weight	w_1	w_2	w_3
Value	1	10	2

We proceed the model validation in a twofold: (a) first we consider the first data points to tune the SIRD＋ ARIMA model and demonstrate its validity; and (b) we tune the SIRD＋ ARIMA model by identifying the time-varying parameters considering the complete available datasets, which gives the model that is used for control purposes. We note that the values used for the social isolation factor are those as exhibited in Fig. 3.

Regarding the first validation part, Fig. 6 shows the estimated time series from the optimization output. Recall that these parameters are piecewise constant for periods of days. In accordance with the discussion in the literature, we can see that these parameters tend to follow stationary trends as the pandemic progresses.

Fig. 6

Parameters , , and identified by the two layer approach — from February , to June , .

The identified auto-regressive model for SARS-CoV-2 transmission parameter is provided in Fig. 7. The ELS procedure yields a regression with daily steps (i.e.  weeks). Fig. 7 demonstrates how the ARIMA model is indeed able to describe the time series , globally catching the time-varying behavior for this epidemiological parameter with 99.04% accuracy.

Fig. 7

Comparison between ARIMA model and previous identified . FIT — 99.04%.

The first data points comprise the period from February to June , . Using the resulting SIRD＋ ARIMA model identified from this period, we proceed by showing the model-based forecasts for the period from June until June , . The forecasts are made using Eqs. (2) with parameters , and given by the ARIMA model in Eq. (7). Fig. 8, Fig. 9, Fig. 10 show the model forecasts compared with real data for active infections, recovered individuals and deaths, respectively, considering a interval of confidence of %. It is possible to notice that the global behavior of the COVID-19 contagion is well described by the proposed models. Complementary, Fig. 11 depicts the model-data error terms (from Eqs. (16)–(18)), given in percentage. The coefficient of determination for these identification results are all above 0.99, which is a clear indicative of good validation.

Fig. 8

Validation of the SIRD＋ ARIMA model using estimated parameters with official data — Active Infected Curve.

Fig. 9

Validation of the SIRD＋ ARIMA model using estimated parameters with official data — Recovered Curve.

Fig. 10

Validation of the SIRD＋ ARIMA model using estimated parameters with official data — Fatal Cases Curve.

Fig. 11

Error between the simulated model and official data (in percentage).

As previously stated, we note that the data-driven model used for the NMPC strategy is based on the identification results for the whole set of available data. Given the high coefficient of determination of the identified SIRD＋ ARIMA model, we can conclude that it describes the COVID-19 pandemic dynamics very well, being able to provide forecasts with relatively small prediction error.

The main innovation regarding the proposed SIRD＋ ARIMA model is that it provides better forecasts than constant parameter SIRD models, which means that the derived control results present more fidelity. Fig. 11 shows that -days-ahead forecasts can be derived with the proposed SIRD＋ ARIMA, whose modeling error is below %. This is a significant prediction horizon, since SIRD models with constant parameters only yield good forecast for horizons of roughly one week [3], [36].

Regarding the control purpose of this work, thus, the proposed SIRD＋ ARIMA enables highly representative prediction estimates for the COVID-19 contagion spread in Brazil. We also note that the proposed NMPC law is designed to generate new social distancing guidelines each week, which makes the application of the SIRD＋ ARIMA model for predictions even more consistent. This topic is further discussed in Section 3.

For perfectly coherent forecasts, it is recommendable for the identification to be re-performed regularly, and the ARIMA fits updated. This advice is given because the COVID-19 pandemic has inherent time-varying properties that may also be influenced by external phenomena, such as the distribution of efficient anti-viral drugs or vaccines. We note that this recommendation has also been provided by Morato et al. [10].

We also remark that the ARIMA fits for the epidemiological parameters are given in weekly samples i.e. , while the SIRD variables are given for each day i.e. . Note that we can represent the weekly-sampled variables, from the viewpoint of the daily-sampled variables, as follows, with days:

Finally, considering this validated SIRD＋ ARIMA model, the effective reproduction number of the SARS-CoV-2 virus in Brazil can be inferred through Eq. (3). In Fig. 12, the evolution of the viral reproduction number for each week of the pandemic in shown. As it can be seen, this reproduction number presents stronger variations at the beginning stage, but, then, tends to converge to steadier values; this corroborates with epidemiological properties. Moreover, we note that the reproduction number for since June , is somewhat steady near 1.603, which indicates that the COVID-19 pandemic is still spreading in the country.

Fig. 12

Effective reproduction number according to the identified SIRDARIMA model.

The NMPC strategy

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

The generated control action represents the input to the population’s response to social distancing measures, as gives Eq. (4). For this reason, it follows that the proposed NMPC algorithm operates with a weekly sampling period (). Implementing a new social distancing guideline every week is coherent with previous discussions in the literature [10]. We stress that it is not reasonable to change the social distancing rules every few days.

Possible control sequences

Before the actual implementation of NMPC tool is presented, one must note that the derived control sequence must be given in accordance with the constraints expressed in Eq. (6). Considering that the NMPC has a horizon of steps (given in weekly samples), from the viewpoint of each week , the control sequence is4 :

Since the variations from each quarantine guideline to the following , denoted are equal to or , it follows that all possible control sequences can be described, from the viewpoint of sample , as:

Therefore, based on Eq. (23), we can conclude that, in the considered settling, there are possible control sequences at each sampling instant.

Control objectives

The main purpose of social isolation is to dilute the demands for hospital bed over time so that all infected individuals can be treated. It seems reasonable to act by the means of social isolation in order to minimize the number of active infections (), while ensuring that the symptomatic compartment remains smaller than the total number of available ICU beds .

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

This trade-off objective (mitigating against relaxing social isolation) is expressed mathematically as: being a nominal limit for (i.e. the initial population size ) included for a magnitude normalization of the trade-off objective. Note that is given within and, thus, there is no need for normalization.

Process constraints

The proposed NMPC algorithm must act in order to address the control objective given by Eq. (24) while ensuring some inherent process constraints. These are:

That the control signal must exhibit a piecewise constant characteristic, as gives Eq. (6); and

that the number of infected people with active acute symptoms does not surpass the total number of available ICU hospital beds in Brazil:

The complete NMPC optimization

Considering previous discussions, the NMPC procedure is formalized through the following optimization setup. This program focuses on the mitigation of the SARS-CoV-2 viral contagion spread in Brazil given the trade-off control objective and the constraints in Eqs. (6)–(25) and the control horizon of (weekly) steps. From the viewpoint of each sampling instant , this optimization is as follows:

Finitely parametrized NMPC algorithm

The finitely parametrized NMPC methodology has been elaborated by Alamir [37]. This paradigm has recently been extended to multiple applications [38], [39], 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, we proceed by detailing how it is implemented at each sampling instant , ensuring that the Nonlinear Optimization Problem in Eq. (26) 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. (23)). Thus, the predicted variables are used to evaluate the cost function . The control sequences that imply in the violation of constraints are neglected. Then, the resulting control sequence is the one that yields the minimal , 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. Fig. 13 illustrates the flow of the implementation of the proposed Social Distancing control methodology for the COVID-19 viral spread in Brazil. We note that this methodology ensures the optimality of the solution regarding the control objective , as states Rathai et al. [38].

Fig. 13

Algorithm Implementation Flowchart.

Simulation results

Considering 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 control results were obtained using Matlab, implemented in a 2.4 GHz Intel i Macintosh computer. The implementation of the parametrized NMPC algorithm follows the lines (and details) of the work by Rathai et al. [38]. The average computation time needed to solve the solution of the proposed NMPC procedure was of .

Simulation scenarios

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

What would be the case if the social distancing policy generated by the NMPC method was enacted since the beginning of the COVID-19 pandemic in Brazil.

What is the possibility of applying the NMPC technique from now on (June , ) and reduce and mitigate the effect of the spread.

In order to provide these scenarios, we consider a “no-control” comparison condition, which represents the simulation of the SIRD＋ ARIMA model with the social distancing factor in fact observed in Brazil taken as input (see Fig. 3). The results are based on days of data. From the 127th sample onward, the “no-control” simulation takes a constant input of (no isolation). The curves corresponding to this condition are denoted as “open-loop”, because no feedback action is involved in this case.

The proposed NMPC is designed with a cost function with tuning weight , which means reducing the number of infections is prioritized over relaxing social distancing. Furthermore, the NMPC optimization is set with a prediction horizon of weeks ( 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 . For simplicity is taken as , since the main goal of the MPC is to ensure .

Scenario (a): NMPC from the beginning

We recall that the COVID-19 pandemic in Brazil formally “started” February , , when the first case was officially registered.

This first simulation scenario considers the application of the NMPC control strategy to guide social distancing starting days after the first case (March , ). We choose this date because, in Brazil, the majority of states started formal 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 ”, since this was not seen anywhere in Brazil. Furthermore, we note that little was known regarding the characteristic of the COVID-19 contagion before April . As detailed in the previous Sections, with daily samples, the proposed SIRD＋ ARIMA model (and identification methodology) is already able to gather consistent parameter estimates, which is essential for coherent implementation of the proposed control strategy.

Considering this control paradigm, Fig. 14 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 show that the NMPC is able to thoroughly attenuate the peak of infections, ensuring that it always stays below the ICU threshold. This is a quite significant result, 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 difficulties to Brazilian hospitals (and health system overall). The results also indicate that the social distancing should vary with relaxing and strengthening periods over the years of , and . 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. We can conclude 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 individuals, with this peak forecast to October , . If the NMPC strategy was indeed applied, two peaks would have been seen, with of total ICU capacity, previewed for September , and March , .

Fig. 14

Scenario (a) Symptomatic Active Infections and Control Input (Social Distancing, ).

We also stress that, in Fig. 14, the observed population response for social distancing is slightly more intense than the actual guideline , for most of the simulation. This is due to the fact that since the ICU hospital beds are saturated, i.e. over 70% of ICU occupancy rate.

Regarding this scenario, Fig. 15 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. [10]. Once again, one should note that the SIRD＋ ARIMA model offers qualitative results. Moreover, while the magnitude of 1.5 million deaths may seem quite alarming, it must be recalled that 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 discussed by THE [5] and Zacchi and Morato [6].

Fig. 15

Scenario (a) Total Active Infections and Total Deaths .

Scenario (b): NMPC from the now on

The second simulation scenario considered refers to the control of the COVID-19 contagion starting at the current moment (), in order to avoid a total collapse of the Brazilian health system. As of this date, the country counts over 1.4 million confirmed COVID-19 cases and more than deaths due to the SARS-CoV-2 virus.

Fig. 16 presents the main results of this second simulation scenario, regarding the active symptomatic infections, which may require ICU hospitalization. Even though a partial catastrophe is already under course in Brazil, if social distancing were guided using the proposed NMPC, a total collapse of the public health system could still be avoided. Fig. 17 shows the evolution for recovered individuals and cumulative number of cases , Fig. 18 gives the evolution of all active infections (symptomatic and asymptomatic), while Fig. 19 presents the mounting number of deaths. The proposed NMPC, if rapidly put in practice, could still be able to slow the viral spread, saving of lives w.r.t. the open-loop/no-control condition. The peak of infections, if such technique is applied, has its forecast previewed to September , , being anticipated in days.

Fig. 16

Scenario (b) Symptomatic Active Infections and Control Input (Social Distancing, ).

Fig. 17

Scenario (b) Cumulative Cases and Recovered Individuals .

Fig. 18

Scenario (b) Total Active Cases .

Fig. 19

Scenario (b) Total Deaths .

Conclusions

In this paper, an optimal control procedure was proposed for the generation of social isolation guidelines used to mitigate the spread of the SARS-CoV-2 virus in Brazil. In order to develop this procedure, a new contagion model was proposed, based on extensions of the SIRD equations. This model embeds weekly auto-regressive dynamics for the epidemiological parameters and also takes a dynamic social distancing factor. 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 was designed within a finitely parametrized input paradigm, which enables its fast implementation.

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

The presented results corroborate the hypothesis formulated in many of the previous papers regarding the COVID-19 pandemic in Brazil [5], [10], [40]: herd immunity is not a plausible option for the country5 ; if no coordinates social distancing action is enforced, the ICU threshold will be largely surpassed, which can lead to elevated fatality.

The prediction of the viral spread evolution is relatively accurate with the proposed adapted SIRD model for horizons of roughly days, with prediction errors below 10%. Larger prediction horizons can be considered, but daily model-updates are recommended.

The simulation forecasts derived with the NMPC strategy and with a no-control condition (open-loop, with no feedback social distancing guideline) 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 symptomatic individuals to late September, , 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 lives if enacted from now (July, ).

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 . Since this has not been pursued, the social distancing measures may go up until late 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 Bastos et al. [13] 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 [6].

Declaration of Competing Interest

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