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

This paper formulates a Model Predictive Control (MPC) policy to mitigate the COVID-19 contagion in Brazil, designed as optimal On-Off social isolation strategy. The proposed optimization algorithm is able to determine the time and duration of social distancing policies in the country. The achieved results are based on data from the period between March and May of 2020, regarding the cumulative number of infections and deaths due to the SARS-CoV-2 virus. This dataset is assumably largely sub-notified due to the absence of mass testing in Brazil. Thus, the MPC is based on a SIR model which is identified using an uncertainty-weighted Least-Squares criterion. Furthermore, this model includes an additional dynamic variable that mimics the response of the population to the social distancing policies determined by the government, which affect the COVID-19 transmission rate. The proposed control method is set within a mixed-logical formalism, since the decision variable is forcefully binary (existence or the absence of social distance policy). A dwell-time constraint is included to avoid too frequent shifts between these two inputs. The achieved simulation results illustrate how such optimal control method would operate in practice, pointing out that no social distancing should be relaxed before mid August 2020. If relaxations are necessary, they should not be performed before this date and should be in small periods, no longer than 25 days. This paradigm would proceed roughly until January/2021. The results also indicate a possible second peak of infections, which has a forecast to the beginning of October. This peak can be reduced if the periods of days with relaxed social isolation measures are shortened.

Introduction

The COVID-19 pandemic seems to be the global health crisis of our time. Scientist first identified this virus (SARS-CoV-2) in humans in Wuhan, in the province of Hubei, China by December 2019. It causes severe acute respiratory syndrome which can become potentially fatal. The WHO estimated by the end of March that the number confirmed cases was reaching the order 70,000, with more than 33,000 confirmed deaths. Now, by the end of April, it has already spread to almost every country of the world, infecting 3,019,246 and killing more than 208,112 people (World Heath Organization, 2020). Its spread is rapid and efficient and it seems that to tackle this pandemic, global scientific efforts are necessary (Bedford et al., 2019). Since vaccines have not yet been developed, most countries have adopted measures to ensure social distancing, aiming to avoid the spread (Adam, 2020). It seems that the COVID-19 has posed an unique question regarding what are the viable public policies necessary to handle its spread. We note that countries have chosen different strategies to determine social distancing measures, as reviewed by Chudik, Pesaran, and Rebucci (2020) and Cohen Kupferschmidt (2020). Up to the Authors’ best knowledge no governmental policy regarding COVID-19 has been based or driven by advanced feedback control strategies, that account for optimality concerns, input saturation and constraints. There are only theoretical results regarding this topic.

The idea behind social distance is to prevent health systems from becoming saturated due to large amounts of COVID-19 patients being treated at the same time. Therefore, with social distancing policies, the health systems do not have to deal with hospital bed shortages associated with a large peak of infections, since the demands for treatment become distributed over time. Fig. 1 illustrates the evolution of symptomatic individuals due to the SARS-CoV-2 virus with respect to no isolation and hard social isolation policies. The threshold represents an estimate for the number of available Intense Care Unit (ICU) hospital beds.

Fig. 1

Necessity of Social Isolation.

However, social distancing measures exhibit at least three ambiguous side-effects, which policy makers should take into account:

If the governments interrupt the social distancing policy before the correct time, the measure is only able to shift the pattern of contamination to the future, which does not help diminishing the problem of saturating the health system due to excessive demand for ICU beds. This topic has been throughly discussed by Hellewell et al. (2020);

Countries that implement rigid social distancing measures have seen devastating economic effects. The recent papers by Eichenbaum, Rebelo, and Trabandt (2020) and Gormsen and Koijen (2020) elaborate on this issue;

A large part of the population may not be immunized and might suffer from future waves of COVID-19 infections, after the social distancing measures take place. Furthermore, recurrent wintertime outbreaks of SARS-CoV-2 virus will probably occur after this initial pandemic dissemination (Kissler, Tedijanto, Goldstein, Grad, & Lipsitch, 2020).

Therefore, it becomes of fundamental importance to predict the correct time and the duration of social distancing interventions. Well-designed social distance policies may help to control the evolution of the disease, to avoid the saturation of the heath systems and to minimize the economic side effects caused by them.

In this paper, the Brazilian context is taken into account (Werneck & Carvalho, 2020). Brazil is a continent-sized tropical1 country and it has already been facing many issues due to the COVID-19 pandemic. The country has 26 federated states, which have been choosing different social distancing measures since mid-March2 The federal government is reluctant to implement nation-wide policies, claiming that the negative economic effects are too steep and that social distancing is an erroneous choice (The Lancet, 2020); the government suggests that the economy cannot stop and that herd immunity could be a solution to this pandemic. The expected impacts of the disease in Brazil are catastrophic (Ismael, Silva, da Silva, de Melo, Neto, de Melo, Macario, de Souza, Moura, & Domenge, Rocha Filho, dos Santos, Gomes, Rocha, Croda, Ramalho, Araujo, 2020). Moreover, due to lack of testing, Brazil is only accounting for patients with severe symptoms or those how have died; therefore, a huge percentage of sub-notification (over 90%) has been reported (Delatorre, Mir, Graf, Bello, 2020, Rocha Filho, dos Santos, Gomes, Rocha, Croda, Ramalho, Araujo, 2020, Silva, Velasco, da Silva Marques, Tibirica, 2020). The daily reports (”measurements”) delivered by the Ministry of Health, collecting number of infected and deceased patients, only gives an impression of the virus spread of past moment, since, on average, a person will exhibit acute symptoms only 20 days of the infection.

The first death due to the SARS-CoV-2 virus in Brazil was registered in March 17, while the first case was officially notified in February 26. Nonetheless, recent papers (Delatorre, Mir, Graf, Bello, 2020, Rodriguez-Morales, Gallego, Escalera-Antezana, Mendez, Zambrano, Franco-Paredes, Suárez, Rodriguez-Enciso, Balbin-Ramon, Savio-Larriera, et al., 2020) point out that the virus was already present in Brazil since the end of January, before the Carnival festivals. Through inferential statistics, Delatorre et al. (2020) acknowledge that community transmission has been ongoing in the state of São Paulo since the beginning of February (over one month before the official reports).

Even though a strong public health system is available in Brazil, as of April 30, many states were already exhibiting a near-collapse situation, with over 95% of ICU beds occupied with COVID-19 patients. This is illustrated in Fig. 2 , which show the ICU occupancy rate and the number of available in many states (Brazilian Federal Medicine Council, 2020). Clearly, the situation is already border-lining.

Fig. 2

ICU Beds and occupancy rate in Brazil, per state (20/04/30).

Motivated by the previous discussion, in this paper we investigate the problem of controlling the evolution of the COVID-19 pandemic in Brazil using optimal social distancing policies, which are designed through a Model Predictive Control (MPC) framework.

The MPC framework is a widespread optimal control method for the control of processes subject to constraints (Camacho & Bordons, 2013). MPC allows to explicitly consider the effect of input, output and state constraints in the control design procedure, which is rather convenient. As any standard discrete control method (Skogestad & Postlethwaite, 2007), MPC-formulated laws stand for piece-wise-constant signals for sampled-data systems, which is clearly the scope of the COVID-19 dissemination process, since it is measured daily through the number of infected and deceased individuals.

The frameworks to predict and control the complex dynamics of the SARS-CoV-2 virus spread are definitely mixed discrete-continuous problems under multiple objectives due to the nature of the problem (daily measurements and piece-wise constant control). MPC certainly fits this context.

We must remark that many research papers have demonstrated the application and validation of this control tool for health-related, biological and ecological regulation purposes (de Ávila-Simas, Morato, Reynalte-Tataje, Silveira, Zaniboni-Filho, Normey-Rico, 2019, Ionescu, De Keyser, Torrico, De Smet, Struys, Normey-Rico, 2008, Moscoso-Vásquez, Colmegna, Sánchez-Peña, 2016, Zurakowski, Messina, Tuna, Teel, 2004). Thus, MPC fits naturally to COVID-19 social distancing control problem. A possible control approach regarding the COVID-19 contagion is to take into account social distancing constraints and, at the same time, minimize the peak of active infections, ensuring this variable stays below than the available amount of ICU beds. Such strategy is indeed viable through the MPC paradigm. The main contributions of this paper are the following:

We present two modified versions of the Susceptible-Infected-Recovered (SIR) model (Kermack & McKendrick, 1927), embedding the effects of social distancing measures in the evolution of the disease;

We propose an additional dynamic state variable, which models the response of the population to social distancing measures enforced by the government. We also use this state variable to forecast the reduction of the speed of transmission of the virus, with respect to enacted distancing policy;

Due to the fact that large error margins have been reported regarding the available COVID-19 data and statistics in Brazil (see (Bhatia et al., 2020)), we perform an uncertainty-weighted Least-Squares criterion to estimate the parameters of the virus infection/spread model, considering both nominal and inconsistent-data conditions.

Based on these uncertain models, we propose an MPC-based control framework to determine in real time whether to apply or not the social distancing policy. This control strategy resides in the solution of a Mixed-Integer Dynamic Programming Problem, at each sampling instant (day), according to new available datasets (number of infected and deaths). The constraints of the MPC procedure are given with respect to the number of available ICU hospital beds in the country. The MPC also accounts for a minimal dwell-time on each control action (no isolation, complete lock-down), so that frequent social distancing policy shifting does not happen.

Note that, for the sake of practical purposes, we estimate model parameters and distancing policies using real data from the Brazilian Ministry of Health, in a fashion similar to the estimation scheme presented by Bastos and Cajueiro (2020).

Our paper relates to some other recent papers that investigate the COVID-19 pandemic from a control viewpoint. There are some recent works that have also inserted a control variable in the available epidemiological model in order to emulate the control of infections using strategies such as vaccination, isolation, culling and self-isolation (Bolzoni, Bonacini, Soresina, Groppi, 2017, Piguillem, Shi, et al., 2020, Piunovskiy, Plakhov, Tumanov, 2020). In particular, we state that the recent paper (Piguillem et al., 2020) has also addressed the issue of optimal social distancing COVID-19 policies. A case study regarding Belgium has been presented in (Alleman, Torfs, & Nopens, 2020), considering a continuous MPC policy; a robust MPC method, regarding the COVID-19 spread in Germany, has been investigated by Köhler et al. (2020).

We must stress that this paper differs from Piguillem et al. (2020), Alleman et al. (2020) and Köhler et al. (2020) in three main points:

We consider the possibility of large amounts of uncertainty regarding the infected/deaths measurements (which stands true for the Brazilian case, which is ours). The previous papers considered relatively small parametric uncertainties added to the model (identified with real data). In our case, we consider large amounts of uncertainty and use uncertainty-embedded models derived through a series of identification runs.

As in (Alleman et al., 2020) and Köhler et al. (2020), we design and synthetize an optimal control strategy with a Model Predictive Control formalism, which has has wide industrial practice and easy implementation. The MPC design in this paper follows a mixed-integer approach and embeds a dwell-time constraint, which had not been tested in previous papers. Furthermore, the chosen control input differs: in Alleman et al. (2020), the control input is the actual isolation parameter, while in Köhler et al. (2020) it directly affects the infection and transmission rates. In this paper, the control input stands for a new model variable which indicates the government enacted social distancing policy, which affects the average isolation observed in the population, that then influences with the transmission and infection dynamics.

We model the response of the population to these social isolation rules with an additional dynamic variable (and fit our variations of the SIR models accordingly).

The considered dwell-time constraint ensures that the social distancing policy remains constant for at least N days, avoiding frequent shifting between isolation and non-isolation states (which obviously would cause ambiguity and confusion to the population).

This work is organized as follows: In Section 2, we introduce the SIR models used in this work and the modifications we make in order to model people’s response to the social distancing policies. In Section 3, we present the approach used to estimate the parameters of the models; therein, we also present the model-data fitting results. Section 4 presents the optimal control scheme used to mitigate the evolution and the spread of the disease along time. Finally, Section 5 presents the main conclusions of the work.

Epidemiological models

Recent literature (Kucharski, Russell, Diamond, Liu, Edmunds, Funk, Eggo, Sun, Jit, Munday, et al., 2020, Peng, Yang, Zhang, Zhuge, & Hong) shows that the infection rate and evolution dynamics of the SARS-CoV-2 virus can be adequately described by Susceptible-Infected-Recovered (SIR) models, which were originally presented by Kermack and McKendrick (1927).

In this Section, we present two modified versions of the SIR model that take into account the effects of social distancing measures and embedded them to the evolution dynamic of the disease. The new variables account for the dynamics of the population response to such social distancing measures (enacted by local governments).

For reading purposes, we present a thorough Nomenclature and Symbology paper in Appendix A.

Original SIR model and modifications regarding COVID-19

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

Susceptible individuals (S), who are prone to contract the disease;

Infected individuals (I), which currently have the disease;

Recovered individuals (R), who have already recovered from the disease.

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

In the SIR model, the parameter β stands for the average number of contacts that are sufficient for transmission of the virus from one individual, per unit of time t. Therefore, βI(t)/N(t) determines the average number of contacts that are sufficient for transmission from infected individuals, per unit of time, to one susceptible individual; and (βI(t)/N(t))S(t) determines the number of new cases per unit of time due to the amount of S(t) susceptible individuals (they are “available for infection”).

Furthermore, the parameter γ stands for the recovery rate, which is the rate that each infected individual recovers (or dies). This parameter characterizes the amount of individuals that “leaves” the infected class, considering a constant probability quota per unit of time.

Based on these definitions, the SIR dynamics are:

Since the SIR model is used herein to describe a short-term pandemic outbreak, we do not consider the effects of demographic variations. Despite recent discussion regarding the possibilities of reinfection(Del Rio & Malani, 2020), we assume that the recovered individuals will not be reinfected (at least for simplicity purposes), i.e. an individual does not contract the disease twice. We will not implement this first SIR model, it is only included for the sake of referencing.

Indeed, the model should also include the dynamic relationships that appear due to the fraction of people that unfortunately die from the disease. Thus, we include a parameter ρ, which stands for the probability of an individual form the infected class I(t) dying from infection before recovering, as suggested in Keeling and Rohani (2011). In this case, the following set of Equations arise:where stands for the number of people from the population that die due to the disease, per unity of time; and D(t) is the number of people that die due to the disease. Note that, in this case, the number of individuals in the population reduces due to the infection according toFor the ease of reference, this adaptation of the SIR model is named hereafter as the “SIRD” (Susceptible-Infected-Recovered-Dead) model.

Since, in the case of the SARS-CoV-2 virus, there is a relevant percentage of the infected individuals that are asymptomatic, we split the class of infected individuals into the classes of symptomatic and asymptomatic individuals, as suggested in Robinson and Stilianakis (2013), Arino, Brauer, van-den Driessche, Watmough, and Wu (2008) and Longini-Jr., Halloran, Nizam, and Yang (2004):where I (R) is the number of asymptomatic infected (recovered) individuals, I (R) is the number of symptomatic infected (recovered) individuals and p is the proportion of individuals who develop symptoms. For ease of reference, this latter model is named hereafter as the “SIRASD” (Susceptible-Infected-Recovered-Asymptomatic-Symptomatic-Dead) model. The SIRASD model has been widely used in the recent literature to describe the COVID-19 pandemics (Piguillem et al., 2020). Just like with SIRD model, the original condition that N(t) is constant over time form the SIR model no longer holds. Therefore, to evaluate the variation of the population size N(t) over time, one needs to integrate

We note that literature presents more “complex” descriptions of the pandemic diseases, such as the “SIDARTHE” model used in Köhler et al. (2020) (that splits the infections into (symptomatic, asymptomatic) detected, undetected, recovered, threatened and extinct). Anyhow, we cannot account for these models regarding the COVID-19 contagion in Brazil, because we have an insufficient amount of data to identify such models. The Ministry of Health only discloses the total amount of infections and deaths due to the SARS-CoV-2 virus at each day. Due to the absence of mass (sampled) testing, no data is available regarding detected asymptomatic individuals or exact ICU bed occupancy rates, as it is available in Germany, for instance, where the paper by Köhler et al. (2020) originates. To choose more complex models than the ones considered in this paper may only decrease the truthfulness/validity of the identification results, since the estimated parameter values may simply represent a singular solution that matches the identification datasets, but cannot be used truthfully for forecasting/prediction purposes.

The new control-embedded models

In order to design and synthesize effective control strategies for social distancing (public) policies, to be oriented to the population by local governments, these previous SIRD and SIRASD models must be adapted to include the dynamics and effects of social distancing. Therefore, a new differential equation is proposed to model the response of population to these social distancing guidelines. This new dynamic equation is appended to those of S(t), I(t) and D(t), which also become affected by the amount of social distancing at a present time.

We will denote ψ(t) as this varying parameter which account for the population response: stands for the case of a complete lock-down (in a hypothetical and unachievable sense that there would not be any more contact between individuals in the population), while stands for a no-isolation state, where the population is behaving “as usual”, with regular activities as before the pandemic. The differential equation which models the evolution of ψ(t) is the following:where u(t) is a binary variable which determines the social isolation policy regulated by the government (this signal will be later on determined by the proposed optimal controller): for the cases when social isolation is determined, an “On” state is set, with ; when the government does not enact an isolation measure, an “Off” state is set, with . Note that α On and α Off are settling-time parameters which relate to the average time the population takes to respond to the enacted social isolation measures.

The dynamic equation above [C] accounts for these two “On”/“Off” possibilities:

“Off”) if the government determines no isolation is necessary, it follows that and meaning that roughly after days, ψ converges3 to 1. In this case, there is no mitigation of the COVID-19 spread;

“On”) if the government determines that a hard isolation is needed, it follows that and . For the case of a constant static gain it follows that roughly after days, the social isolation variable ψ converges4 to where ψ inf stands for the ”hardest” isolation observed in practice. The time-varying gain K gives the relationship between this “hardest” isolation factor and the social distancing guideline (control input u), as gives the following law:with some positive γ. This relationship implies that K(t)ψ inf decreases (yielding stronger social isolation) as more active infections occur.

Since [C] represents a first-order differential system, the settling-time constants α On and α Off determine the convergence speed of ψ(t). Notice that, for an arbitrary system x(t) reaches 0.99x in 5/α units of time, since the solution for this differential equation is with α > 0, for which .

Since ψ(t) is the average people’s response to public policies to reduce the spread of the virus (such as isolation measures or incentive to wear masks), it affects the transmission factors and, thus, we replace β in the SIRD model by ψ(t)β and β and β in the SIRASD model, respectively, by ψ(t)β and ψ(t)β. It is worth mentioning that the population response to the On-Off isolation policy control may depend on factors such as the incremental number of deaths and the amount of information they known about the disease. Anyhow, since this factor depends on people’s choices, and thus becomes rather difficult to quantify, we opt for the simplicity of the previous differential equation [C], since we may partially estimate its parameters with the available data.

Finally, we include these new dynamics into the SIRD and SIRASD models as discussed. In this case, we get the following models that we denote SIRDC (Susceptible-Infected-Dead-with Control) model:and SIRASDC (Susceptible-Infected-Asymptomatic-Symptomatic-Dead-with Control):which will be used for identification and control purposes.

We must stress that all parameters and variables of the SIRDC and SIRASDC models, used for control synthesis and simulation, are positive.

Tuning of the SIRDC / SIRASDC models

In this Section, we present the results concerning the estimation of the epidemiological parameters of Eq. (7).

We must, at first, re-affirm that recent literature regarding the Brazilian COVID-19 context has raised attention to the large amount of sub-notification in the country (Bhatia, Cori, Parag, Mishra, Cooper, Ainslie, Baguelin, Bhatt, Boonyasiri, Boyd, Cattarino, Cucunub, Cuomo-Dannenburg, Dighe, Dorigatti, Ilaria van-Elsland, FitzJohn, Fu, Gaythorpe, Green, Hamlet, Haw, Hayes, Hinsley, Imai, Jorgensen, Knock, Laydon, Nedjati-Gilani, Okell, Riley, Thompson, Unwin, Verity, Vollmer, Walters, Wang, Walker, Watson, Whittaker, Wang, Winskill, Xi, Ghani, Donnelly, Ferguson, & Nouvellet, Rocha Filho, dos Santos, Gomes, Rocha, Croda, Ramalho, Araujo, 2020, Rodriguez-Morales, Gallego, Escalera-Antezana, Mendez, Zambrano, Franco-Paredes, Suárez, Rodriguez-Enciso, Balbin-Ramon, Savio-Larriera, et al., 2020, Silva, Velasco, da Silva Marques, Tibirica, 2020, The Lancet, 2020). Therefore, we progress by embedding the uncertainty regarding the available datasets to the identification procedure.

The following identification is based the real data provided by the Ministry of Health of Brazil from February 25, 2020 to May 8, 2020. Specifically, we consider the cumulative number of infected individuals () and the number of deaths D(t).

To incorporate the issue of sub-notification, we assume that the data provided by the Ministry of Health is corrupted. Instead of using the original data, we assume: where the instances of Z nominal(t) and D nominal(t) represent the data provided by Ministry of Health. The parameters q ∈ [0, 1] and q ∈ [0, 1] are uncertainty measures that provide a relationship between the nominal (observed) data and the real (latent) variable. We estimate all parameters of our models by minimizing the square-error between the integrated variables and their real values, according to regular identification methodologies (Bard, 1974, Brauer, Castillo-Chavez, Feng, 2019). We proceed by following an hierarchical procedure as done in Bastos and Cajueiro (2020).

Our identification procedure also includes a limit5 to the parameter values, as follows: β, β ∈ [1/20, 2], γ, γ ∈ [1/14, 1/2], ρ, ρ ∈ [0.001, 0.2], α ∈ [0, 1] and ψ inf ∈ [0.3, 0.7]. These limits are in accordance with those presented by Bastos and Cajueiro (2020); Werneck and Carvalho (2020).

Least square procedures

Firstly, we consider the estimation of the infection, transmission and death probability, β, γ and ρ, respectively, of the SIRD model in Eq. (3). This procedure is followed by taking into account that no social distancing policy was enacted in Brazil from February 25, 2020 to March 22, 2020, i.e. for this period; Bastos and Cajueiro (2020) present a through discussion on this matter. The transmission rate parameter γ is then fixed for all the subsequent identification procedures, since this parameter is characteristic of the disease. We refer to the β value of this step as b 1.

Secondly, we assume that during this first period (no social distancing) and estimate the parameters of the SIRDC model from Eq. (6) by minimizing the following square-error:where Z(t) and D(t) represent the data provided by the Ministry of Health of Brazil embedded with uncertainty, as gave Eqs. (8)-(9), and are estimated parameter values using the SIRD model. Note that we use the nonlinear function to correct the exponential characteristic of the series so that the errors of the last values of the series do not dominate the minimization. We assume that there are no recovered individuals at the beginning of the series, and also that (no-isolation state). We refer to the parameters values of this step as and . We limit β in with . For the models with uncertainty, α and ψ inf are also limited in and with and that is, these values under uncertainty are limited by a range defined by the simulation without uncertainty.

We stress that the nonlinear map f(z) is increasing in its domain, and so is f(z)2. This function further weights the last values of the series w.r.t. to the first data steps. Moreover, the formulation in Eq. (10) is coherent with regular Least-Square identification procedures, which minimize the squared difference between model and data.

Thirdly, we estimate the complete SIRASD model. For such, we assume that asymptomatic infected comprise all individuals without symptoms and also those with mild symptoms (that do not need ICU beds), and that symptomatic infected are individuals with moderate to severe symptoms. This line of thought follows the orientation given by the Brazilian Ministry of Health, that incentives people to only seek medical attention if symptoms are moderate or severe, and to stay home otherwise. Therefore, we suppose that and which corresponds to the simulation with uncertainty only upon the number of deaths (i.e. q ≠ 1 and . Since the uncertainty for the number of infected is also related to the asymptomatic individuals, we use . We use this value for the infection probability as constant, so that the initial condition for the asymptomatic infected is . In this context, we estimate the parameters β and ψ inf in order to minimize the following square-error:where and represent the values obtained with the SIRASD model.

Obtained models

The previous identification procedure was realized for a large number of possibilities of uncertainty (q and q). We account for three models:

The ”Nominal” SIRDC model, which is tuned for ;

The ”Uncertain 1” SIRASD model, which is derived from the identification procedure considering 50% more sub-notified deaths and 30 times more infected individuals than reported ( and ;

And the ”Uncertain 2” SIRASD model, which is likewise found through the identification procedure for 50% more sub-notified deaths and 15 times more infected individuals than reported ( and .

The respective model parameters, found through the identification procedure detailed in Section 3, are presented in Table 1 . Regarding the population response to social distancing guidelines, the settling-time parameters and maximal isolation fator values are given in Table 2 .

Table 1

Obtained [SIRD]/[SIRASD] model parameters.

Model	βA	γA	βS	γS	p	ρ
Nominal	–	–	0.4230	0.1395	–	0.09917
Uncertain 1	0.3689	0.0952	0.4307	0.1395	0.0625	0.1462
Uncertain 2	0.4307	0.1395	0.4307	0.1395	0.0322	0.1461

Table 2

Obtained [C] model parameters.

Model	αOn	αOff	ψinf
Nominal	0.1483	0.2966	0.4914
Uncertain 1	0.1462	0.2924	0.5721
Uncertain 2	0.1400	0.2801	0.5664

We must note that the minimization procedures detailed above could be likewise applied to being z the data for a given variable and its model equivalency. In such case, one should consider to exclude possible parameter BIAS. The parameters estimated with both approaches are roughly the same.

Some results

To conclude this Section, we show some short-term simulation results in Fig. 3 . This Figure represents the SIRASDC model running with respect to known data, from 20/03/27 until 20/05/08. The differences between the speed of the COVID-19 spread considering nominal and uncertain conditions are considerable. In Fig. 4 , the data is shown in a logarithmic scale.

Fig. 3

Short-term simulation for the SIRASD model: Nominal vs. Uncertain.

Fig. 4

Short-term simulation for the SIRASD model: Nominal vs. Uncertain (Log-scale).

These Figures depict the differences between the models and also marks how the nominal model perfectly fits the dataset from the Ministry of Healthy. Clearly, the amount of sub-notification plays a significant role in the analysis of this pandemic. If one does not take the uncertainty into account, wrong and rash decisions may be performed. We must stress that an increase on the number of infected individuals (larger sub-notification) means that the mortality rate of the disease decreases, while it increases for a larger sub-notification margin with respect to the number of deaths.

The predictive control strategy

In this Section, we develop the second main contribution of this work, i.e. the optimal control strategy aiming to mitigate the impacts of the COVID-19 pandemic. In practice, the resulting control law which should be implemented by the means of social distancing policies, conducted through orientation by the local government.

For this goal, we consider the control-appended models (SIRDC and SIRASDC), which are now regulated under a closed-loop scheme. In fact, the proposed control strategy is formulated with respect to SIRASDC model, as explained in the sequel.

We must consider that the predictive control strategy is to be synthesized based on a model that a priori encompasses the dataset mismatches (uncertainty in terms of sub-notification). Taking into account a model that considered that all measures (deaths, infections) are not exact can further improve the outcome of the control policy, since it will be more conservative with respect to social isolation measures. In some sense, this kinds of strategy is a robust MPC feedback procedure, because it is based on a worst-case pandemic level. The application of the MPC feedback without the uncertain realization of the model can lead to performances which may borderline the use of the available ICU beds in the country, i.e. the nominal prediction model may give a result in terms of infections which is less than what is observed in practice. This situation would lead to necessity of possibly longer periods of social isolation. Then, such robust MPC procedure is able to avoid this kind of situation and it can be able to significantly reduce the number of fatalities.

Control objectives

As previously discussed, social isolation measures are necessary since they are able to ”flatten” the COVID-19 spread curve. Fig. 1 illustrates this issue, which shows a simulation for the evolution of symptomatic individuals I(t) over time with no social isolation and hard social isolation ( which represents a complete lock-down situation). Clearly, when social isolation policies are implied (mathematically, for ψ < 1), the infection peak is postponed and reduced.

Therefore, the control objectives regarding COVID-19 isolation policies are the following:

To reduce the peak of symptomatic individuals, as much as possible;

To ensure that the peak is smaller than the number of available ICU beds (a tolerance factor can be included);

To determine isolation policies for as little time as possible (in order to mitigate the effects of isolation on economy);

And to avoid shifting between states (isolation or not), maintaining a minimal period of N days in each condition.

These objectives can be mathematically expressed, respectively, as follows:

By minimizing the amount of symptomatic individuals I;

By ensuring that for all t, where n represents the amount of ICU beds and ξ ∈ (0, 1) a tolerance factor (which should also be minimized);

By minimizing the isolation policy u(t), which should be 0 (no isolation) for as long as possible;

To ensure that u(t) is piece-wise constant and maintained at each state (0 or 1) for N samples, where stands for the sampling unit and T the sampling period. We denote this as a dwell-time constraint on u.

Control law

As remarked by Silveira and Pagano (2005), the implementation of any feedback control policy derived from u(t) to real biological system may light upon two complicating issues:

Feedback control in the continuous time-domain t would require the measurement of the involved variables at every instant of time, which is obviously not possible. Anyhow, since the COVID-19 spread has ”slow” dynamics, discretized versions of the SIRDC / SIRASDC continuous models can be used. Given that new measurements are available each day, the discrete sampling is of ;

The control signal should, essentially, model the human action on the studied ecosystem. As previously argued, the human correspondence to the isolation policies have already been included to the SIRDC and SIRASDC models, through the dynamics of ψ(t).

The discretization method used in this paper is very usual: all derivative functions are approximated by the deviance along on sampling period, this is: . This is usually referred to as Euler/forward discretization (Skogestad & Postlethwaite, 2007).

Henceforth, this paper considers a piece-wise constant control signal u(t), generated from periodic measurements, every day. Considering control increments denoted Δu(k), it is implied that:

The actual control law that the controller applies to the system is for the whole time during each sampling period interval, i.e. . This is clearly a piece-wise constant signal due to Eq. (12).

Apart from being piece-wise constant, the control law u must obey another main restriction, in order to be implementable in practice: it must depart from which is the last observed social distancing policy (the country is still in an isolation condition). It must be noted that stands for the instant corresponding to the last sampled field data (infected, deaths dating 29 April).

The optimal On-Off control framework

To design an optimal controller which determines the On-Off control policy u(t) according to the previous discussion, we will follow the Model Predictive Control formalism.

Sometimes named moving/sliding horizon control, the MPC concept is quite straightforward for the optimal control of constrained systems. The basic MPC formulation resides in the solution of an optimization problem with respect to a sequence of control actions U, at each discrete instant. This optimization is written in terms of a process prediction model, performance goals and constraints, which are handled explicitly. The essence of the MPC framework resides in solving this optimization procedure at each discrete time step k and applying the first entry of the control sequence solution U; then, the horizon “rolls forward”.

Therefore, we consider an well-posed quadratic function J which is minimized seeking to accomplish objectives (a) and (c) (reducing infections and reducing the social isolation periods) presented in Section 4.1. This cost function is analytically expressed as:where N is a given prediction horizon, and Q ∈ (0, 1] is a weighting scalar. This weighting is introduced to induce a trade-off between strengthening social isolation measures (Q closer to 1) to mitigate the number of active symptomatic infections and relaxing it, to reduce economic backlashes (Q closer to 0). The notation stands for a model-based prediction for instant made at instant k. The constant stands for the maximal possible value of I, with respect to open-loop simulations (as in Fig. 1); we include this constant to ensure that the magnitude of the first and second term of J are the same (i.e. normalization).

The vector of control efforts inside the prediction horizon U (to be optimized) is also presented:From this control sequence, at each sampling instant k, one takes the first entry u(k) and applies the control signal according to Eq. (12) to the controlled COVID-19 spread process (SIRDC / SIRASDC models).

Notice that if one simply minimized the previous cost function J, at each sampling instant k, with respect to a control sequence U, the results would be a control sequence U which provides a trade-off (according to weighting scalar Q) between the minimization of infected individuals and control effort (u gets closer to 0 and, thus, less isolation is implied).

Anyhow, for an appropriate application of this paradigm, the constraints of each as given in Section 4.2, should be taken into account by the optimization procedure. It follows that:

We stress that:

Eq. (15) ensures that the peak of infections is reduced and does not surpass ;

Eq. (16) ensures the control signal is bounded within the social isolation limits, Eq. (17) ensures that this law is piece-wise constant and Eq. (18) implied that the variation is binary (so that the control is, in fact, ”On”/”Off”);

Eq. (19) implies that a minimal dwell-time of N samples must be accounted for, i.e. the control determines that u stays at a given state ”On” or ”Off” for a minimal period of N days. It is implied that N > N.

Notice that the slack/tolerance variable ξ defined in Section 4.1 should also be minimized so that it does not become a freely-moving degree-of-freedom of the optimization procedure, but a constrained variable. Therefore, a ξ penalty is included to the cost function J, as follows:

We note that to use a quadratic formulation of the cost function J, as given in Eq. (20), is seen in quite a few MPC applications (Camacho & Bordons, 2013). For convenience, we proceed with a quadratic weight on the minimized variables. Anyhow, we must note that, since the considered epidemic control problem has only positive variables (I and u are inherently positive due to the nature of the mathematical problem that describes the COVID-19 contagion), the cost function J could be alternatively formulated as a linearly dependent cost on I and u, which would lead to similar results. Using linear cost functions in MPC applied for epidemic control has been previously seen in Nowzari, Preciado, and Pappas (2016), Köhler, Enyioha, and Allgöwer (2018) and Watkins, Nowzari, and Pappas (2019).

Therefore, bearing in mind this previous discussion, the MPC approach to mitigate the effect of COVID-19 spread consists in minimizing the cost function J at every discrete-time step k, with respect to the previously discussed constraints and taking into account a discretized version of the SIRASDC model, with day. One can mathematically express this problem as follows:

We note that the “Discrete SIRASDC model” used as a prediction model in the MPC optimization is simply found by applying an Euler discretization method (as detailed in Remark 6) to Eq. (7). This discretized model is then extended as a prediction model, describing future instants with respect to the information at (k).

Simulation results, forecasts and discussion

In this Section, we present simulation forecasts using the SIRDC/SIRASDC with parameters identified in Section 3. The following results were obtained with the aid of Matlab software, Yalmip toolbox and Branch-and-Bound (BNB) solver.

In the sequel, the baseline threshold in the I curves represent the number of available ICU beds in the country (see Fig. 2). The maximal threshold stands for an incremented number of ICUs (twice the baseline value), accounting for field hospitals and emergency ICUs that have been made specifically for the COVID-19 pandemic.

The following control results were obtained considering the SIRASD models (Nominal, Uncertain 1 and Uncertain 2), as presented in Section 3.

Before presenting the actual results, we must affirm that the forecasts and arguments that we present in the sequel should not understood by the reader as incontrovertible truths. These forecast are model-based simulations which depend on a number of factors and initial conditions. Furthermore, we must stress that we have aggregated the whole set of Brazilian data in order to provide a general view of the country. However, if anyone intends to use the proposed method to help the formulation of public health policies, we suggest its application to datasets of smaller regions, that share the same hospital chain. We note, as illustrated in Fig. 2, that different regions of the country are facing different levels of the pandemic.

As evidenced in Section 3, the used models grasp the behaviors the SARS-CoV-2 virus dynamics quite accurately, but this does not means that the future predictions are unmistakable. For a fact, we cannot ensure that the social isolation measures will be strictly followed by the population, as we cannot ensure that other factors may come to help ease the spread of the disease (such as vaccines). What we mean by this is that the goal of this work is guide public policies regarding social isolation specially by taking into account the role of uncertainty and sub-notification.

Due to the fact that the proposed model cannot exactly predict the pandemic dynamics, to apply some control polcy conceived based on a nominal model may lack conservatism. This could lead to catastrophic results, risking high levels of mortality. Any possible control policy that the government implements through social distancing measures must be based on recurrent (worst-case) model parameter estimations and recalculations of the optimization problem. One cannot use the models derived with the parameters presented in Table 1 as if they would not change along time. The correct measure is to take into account uncertainty-embedded models, performing the identification procedure detailed in Section 3 every day (when new datasets are available). Such adaptive control procedure (with model tuning and model-based control optimization) would be much more prudent, requiring constant measuring, monitoring, parameter estimation and control computations. As discussed by Köhler et al. (2020), feedback is utterly necessary to ensure a reliable handling of the SARS-CoV-2 outbreak. This is especially critical in Brazil, due to the high level of uncertainty on the datasets.

Through the sequel, the dashed lines represent the results with uncertainty (solid dash, Uncertain 1 model; dot dash, Uncertain 2 model), while the solid lines account for nominal conditions.

Considering these models, Fig. 5 shows the simulation for roughly one year dating from the last data sample (20/03/17), considering a total lock-down condition () and a no-isolation () case. The first 52 samples represent the known dataset, whereas the following samples stand for predicted data. Clearly, even if a hard isolation is enacted, the Brazilian health system will still face issues with large amounts of COVID-19 patients, with a nominal peak forecast to 26th of May. The nominal collapse of the healthy system (threshold) dates very soon, May 23. The amount of deaths expected with the uncertain model is unprecedented. Of course, each life matters and 2 million deceased individuals is a lot to bare. Psychological and social traumas will mark the country. A hard isolation could be able to save more than one million lives, taking into account the results achieved with the worst-case uncertainty scenario.

Fig. 5

Necessity of Social Isolation 2: Model-based Forecasts.

With respect to the forecasts presented in Fig. 5, we must also stress that the possibility of herd immunity must be discarded. These results corroborate the conclusions presented by Köhler et al. (2020), which indicate that neither a complete eradication of the virus nor herd immunity are possible options to attenuate the COVID-19 pandemics without the availability of a vaccine. These results also go along the lines of (Hellewell et al., 2020).

Considering control results, the MPC optimization procedure from Eq. (21) is solved for different cases of N (minimal amount of days in each state: isolation, no isolation). For such, the weighting scalar Q is chosen aiming to imply an adequate trade-off between peak reduction and social isolation. Since the occupancy rate of ICU beds in the country dating 20/05/08 is considerably high, we chose which means that the MPC makes ”more effort” to reduce the amount of infected individuals then to restore a no-isolation policy, which is reasonable considering the observed situation. The control horizon is fixed as days (the MPC makes predictions for two-months ahead of each sampled k, day).

Firstly, we show the results for 5 and 7 days of the minimal days in each state condition. The decision by the MPC optimization and the resulting enacted social isolation measure (ψ) are shown in Fig. 6 . The resulting effects on the amount of symptomatic individuals is shown in Fig. 7 . We must state that, for N ≤ 7 days, the amount of shifting in the observed social isolation variable is quite intense. Furthermore, the obtained results with these values for N were not enough to reduce the infection peak, as observed. Note that this kind of policy would hardly ever be implementable, since a confusing message would be passed to the society. The frequent changing between isolation or no-isolation would not be strictly followed, which is certainly unwanted. Therefore, these results are not considered as practical or viable.

Fig. 6

Control Policy, 5 and 7 Days, Excessive Shifting.

Fig. 7

5 and 7 days, Resulting regulation.

Then, Fig. 8 presents the obtained control results for and 60 days, regarding the I curve. This Figure contains a lot of information, which we try to explain by parts:

Fig. 8

and 60 Days: Infected with Symptoms.

If social isolation is not maintained until roughly August 14, even the shortest “openings” (days in reduced isolation policies) could be catastrophic. Any possible reduction of the hard social isolation measures, before this date, would result in an infection peak which would surpass the amount of available ICU beds in the country in over seven times (considering the worst-case uncertainty). This is very significant and thus, any possible social isolation reduction should not proceed before the initial infection curve starts to decay (in both nominal and uncertain conditions);

Therefore, smallest peaks of infections occurs if the isolation measure is kept at least until August 14. After this date, when relaxations in these measures are enacted, a second infection peak will certainly appear. This second peak is due to the fact that social isolation is reduced by the MPC law after the decay of I running from first peak (August 14). This second peak dates roughly October 3.

The second peak of infection is reduced with a smaller number of days in a no-isolation mode. This means that, after August 14, the MPC control action which results in the smallest values for I are those with, at most, periods of 25 days in the no-isolation mode. The smallest amount of ”open” periods, better the results, as expected.

Note that as N increases, this second peak of infection also increases because the amount of minimal days determined for a no-isolation (or reduced-isolation) policy forces a peak increase, which is later treated by a total isolation after its decay. The amount of deaths are given in Fig. 10. Depending on the amount of days in a no-isolation condition, the amount of deaths may range from 0.13 to 1.88 million individuals, according to the uncertain models.

Fig. 10

and 60 Days: Deaths.

Regarding these obtained results, it seems reasonable to us to ponder the following issues:

The uncertain model forecasts quite harsh infection scenarios. Even though the considered uncertainty is quite high (15, 30 times more cases), it offers us a worst-case forecast to determine public policies. With such uncertain model in mind, it seems evident why social distancing measures are so important right now and why they should not be dropped, despite their possible economical side-effects. It seems extremely necessary for public policies to offer alternative solutions to those without jobs or economically suffering due to the social isolation.

No social distancing measures should be relaxed before mid-August (20/08/14). This would definitely help in avoiding the collapse of the Brazilian health system.

If social distancing is to be relaxed, this should not be done before the first infection peak starts to decay (beginning of August) and the no-isolation periods should be the minimal amount of days possible. To ensure heath safety, a conservative measure indicates that such paradigm of recurrent short periods of reduced isolation, followed by hard isolation periods would proceed until roughly 2021/01/21. This paradigm would be helpful to ensure that the SARS-CoV-2 virus does not cause further infection peaks and to mitigate the amount of deaths.

The resulting control policies from the MPC procedure, for the different values for N, are shown in Fig. 9 . These curves indicate, roughly, when to determine social isolation measures and when to set them off. In fact, the actual implemented policy would depend on a daily update of the MPC results with measured datasets. Anyhow, these results indicate a forecast of roughly when to determine or call off these measures. The best result, in terms of infections, would be to follow the Social Isolation state until mid-August, an then relax this measure with small periods (that should definitely not surpass 25 days).

Fig. 9

Control Policy: and 60 Days of Social Isolation.

We note, that results may also differ according to the chosen trade-off weight Q in the MPC cost function J. These results were elaborated with a tuning weight that “forces” stronger isolation measures than “relaxations”. Anyhow, this choice seems coherent with the alarming scenario in Brazil, with rapidly increasing contagion curves.

It seems, mathematically speaking, that even if the SIRASD model has a new degree-of-freedom (which is the decision variable u, to determine when to determine social isolation), the resulting optimization points out that the best option is to maintain isolation for as long as possible (considering ). Even if allowing social contact for a while, the optimization finds minima solutions of J for the smallest number of days with contact and, then, once again determines isolation. Even though different curves can be retrieved with alternative values for Q, it is very important to bare in mind that situations with I > n may cause an unprecedented health crisis in the country.

The COVID-19 is quite worrisome and presents devastating social and economic effects. Biology literature points out that social isolation is necessary. Using mathematical models and optimization, the answer is the same.

Conclusions

In this paper, we investigate an optimization-based solution for social isolation measures of the COVID-19 spread for the Brazilian context. Since recent works have warned against the large order of sub-notification in Brazil, we take uncertainty into account to determine nominal and uncertainty dynamic models of the COVID-19 pandemics. Such uncertainty-embedded models are SIR-kind equations which also consider a new variable, which accounts for the average response of the population to social distancing measures (as determined by the government). A robust Model Predictive Control framework is designed for the regulation of the COVID-19 through the means of such social isolation policies. The MPC is derived as an optimal On-Off social distancing planner.

In this paper, we have tried to expose some essential insights regarding sub-notification and how possible relaxations of social distancing can be performed in the future. Below, we summarize the main findings of this paper, enlightening the key points:

The presented results corroborate the hypothesis formulated in Hellewell et al. (2020) and also discussed in The Lancet (2020), with respect to the Brazilian scenario: herd immunity cannot be considered a plausible solution, offering great risk and leading to elevated fatality. Furthermore, as illustrate (Rocha Filho, dos Santos, Gomes, Rocha, Croda, Ramalho, Araujo, 2020, Rodriguez-Morales, Gallego, Escalera-Antezana, Mendez, Zambrano, Franco-Paredes, Suárez, Rodriguez-Enciso, Balbin-Ramon, Savio-Larriera, et al., 2020, Silva, Velasco, da Silva Marques, Tibirica, 2020), vertical isolation is also not an option for the time being, since we do not have the means to formulate an efficient public policy to separate the population at risk from those with reduced risk, due to multiple social-economical issues of the country.

Since the spread of the SARS-CoV-2 virus is inherently complex and varies according to multiple factors (some which are possibly unmodelled and external), exact prediction of the pandemic dynamics is not possible. Therefore, the correct control procedure should be based on a recurrent (daily) model tuning and re-calculation of the control law, always taking into account the uncertainty margins.

The simulation forecasts found through the MPC optimization procedure, which accounts for the uncertainty in the spread of the disease, indicate that, at least for now, only one answer is available: maintain social isolation for as long as possible, without relaxing it before mid August 2020. This is a rather strict suggestion, but seems to be the sole possible way to attenuate the (already high) levels of the virus in Brazil. The forecasts also indicate a prediction for the infection peak in the country dating very soon, May 26, with a second (and larger) peak possibly arising in October. The control policy, in terms of social isolation, shows that relaxations (loosening the isolation measures) should be performed in, at most, periods of 25 days of reduced-isolation, after the first infection peaks has passed, until roughly January 2021.

The experience and results provided in this paper could also serve for future pandemic outbreaks, for which social distancing policies can be formulated at the beginning of the spread following the mathematical feedback control methods herein proposed.

Synthetically, we must stress that this paper presents only qualitative results of how an optimization-based On-Off strategy can be formulated regarding the COVID-19 spread, regarding the Brazilian context. Since the country as been experiencing an unwillingness to formally start harder social isolation measures (The Lancet, 2020), the social and economic costs of the pandemic might be brutal.

The qualitative results presented in this paper can definitely serve for long-term regulatory decision policies in Brazil. These results indicate that social distancing measures will be recurrent and ongoing until 2021 if no vaccines are developed by then. Therefore, compensatory social aid policies should be deployed to mitigate a long-lasting economic downfall. We note that some recent findings also discuss this matter, suggesting the need for long lasting income aids during the periods of social isolation (Ahamed, & Gutierrez-Romero, Zacchi, & Morato). The Authors truly hope that the proposition herein formalised can serve to help determining adequate public (health and social) policies from now on.

Finally, we note that Brazil is a very large country, many different demographical and geographical differences. To further investigate the issue of uncertainty and sub-reports regarding the COVID-19 contagion in the country, local data and conditions should be analysed (as per state). The Authors plan on investigating this topic as future works. In these studies, the proposed MPC strategy could be designed with regard to local (uncertainty-embedded) models for each location, determining local social distancing guidelines. The approach in this paper gives an average guideline for the whole country.

Declaration of Competing Interest

There is no conflict of interest derived from this manuscript. Authors state that this work has not been published previously, that it is not under consideration for publication elsewhere, that its publication is approved by all authors and tacitly or explicitly by the responsible authorities where the work was carried out. We also accept that if this submission is accepted, it will not be published elsewhere in the same form, in English or in any other language, without the written consent of the Publisher.

Variable	Meaning	Unit
t	Continuous time variable	days
T_s	Sampling period	day
k	Discrete time variable	day
S	Susceptible individuals	number of people
I	Active infected individuals	number of people
I_S	Active infected individuals displaying symptoms	number of people
I_A	Active infected individuals that do not display symptoms	number of people
Z	Cumulative number of infected individuals in the country	number of people
R	Recovered individuals	number of people
R_S	Recovered individuals that had symptom	number of people
R_I	Recovered individuals that did not display symptoms	number of people
D	Deceased individuals	number of people
N	Total population size (Brazil)	number of people
β	Average number of contacts sufficient for viral transmission	1day
β_A	Transmission rate for the symptomatic class	1day
β_SA	Transmission rate for the asymptomatic class	1day
p	Symptomatic/Asymptomatic individual rate	–
γ	Recovery rate	–
γ_S	Recovery rate for the symptomatic class	–
γ_A	Recovery rate for the asymptomatic class	–
ρ	SARS-CoV-2 lethality rate	–
ψ	Response of the Population to Social Isolation Guidelines	–
ψ_inf	“Hardest” Social Isolation factor	–
u	Social Isolation Guidelines / Control Input	–
Δu	Control increment	–
α_Off	Settling time parameter for “No Isolation” mode	1day
α_On	Settling time parameter for “Social Isolation” mode	1day
K_ψ	Time-varying social isolation gain w.r.t. active infections	–
γ_K	Parameter this time-varying relationship	–
q_D	Sub-report uncertainty w.r.t. to deaths	–
q_I	Sub-report uncertainty w.r.t. to cumulative infections	–
f(·)	Nonlinear identification function	–
n_ICU	Total number of available ICU beds	–
J	MPC Cost function	–
N_p	MPC Prediction Horizon	days
Q_I	MPC trade-off weight	–
ξ	Slack variable	–