A Parametrized Nonlinear Predictive Control Strategy for Relaxing COVID-19 Social Distancing Measures in Brazil
Marcelo M. Morato, Igor M. L. Pataro, Marcus V. Americano da Costa, Julio E. Normey-Rico
AA Parametrized Nonlinear Predictive Control Strategy forRelaxing COVID-19 Social Distancing Measures in Brazil
Marcelo M. Morato a , Igor M. L. Pataro c , Marcus V. Americano da Costa b , Julio E. Normey-Rico a a Renewable Energy Research Group (
GPER ), Department of Automation and Systems (
DAS ),Federal University of Santa Catarina (UFSC), Florian´opolis, Brazil. b Department of Chemical Engineering (
DEQ ), Federal University of Bahia (UFBA),02 Professor Aristides Novis St., Salvador, BA-40210910, Brazil. c CIESOL, Department of Informatics, University of Almer´ıa, Ctra. Sacramento s/n 04120, Almer´ıa, Spain
Abstract
The SARS-CoV-2 virus was first registered in Brazil by the end of February, 2020. Since then, the countrycounts over 70000 deaths due to COVID-19, and profound social and economical backlashes are felt allover the country. The current situation is an ongoing health catastrophe, with the majority of hospitalbeds occupied with COVID-19 patients. In this paper, we formulate a Nonlinear Model Predictive Control(NMPC) to plan appropriate social distancing measures (and relaxations) in order to mitigate the pandemiceffects, considering the contagion development in Brazil. The NMPC strategy is designed upon an adapteddata-driven Susceptible-Infected-Recovered-Deceased (SIRD) contagion model, which takes into accountthe effects of social distancing. Furthermore, the adapted SIRD model includes time-varying auto-regressivecontagion parameters, which dynamically converge according to the stage of the pandemic. This new modelis identified through a three-layered procedures, with analytical regressions, Least-Squares optimization runsand auto-regressive model fits. The data-driven model is validated and shown to adequately describe thecontagion curves over large forecast horizons. In this model, control input is defined as finitely parametrizedvalues for social distancing guidelines, which directly affect the transmission and infection rates of theSARS-CoV-2 virus. The NMPC strategy generates piece-wise constant quarantine guidelines which can berelaxed/strengthen as each week passes. The implementation of the method is pursued through a searchmechanism, since the control is finitely parametrized and, thus, there exist a finite number of possiblecontrol sequences. Simulation essays are shown to illustrate the results obtained with the proposed closed-loop NMPC strategy, which is able to mitigate the number of infections and progressively loosen socialdistancing measures. With respect to an “open-loop”/no control condition, the number of deaths still couldbe reduced in up to 30 %. The forecast preview an infection peak to September 2nd, 2020, which could leadto over 1 . Keywords:
Nonlinear Model Predictive Control, COVID-19, Social isolation, SIRD Model, SystemIdentification.
1. Introduction
The SARS-CoV-2 virus was first registered in humans in Wuhan, China by December 2019. Since then,it has reached almost all countries around the globe with pandemic proportions and devastating effects. Thiscontagion seems unprecedented and can already be rules as the definite health crisis of the 21 rst century.The SARS-CoV-2 virus causes a severe acute respiratory syndrome, which can become potentially fatal.
Email address: [email protected] (Marcelo M. Morato)
Preprint submitted to ISA Transactions July 21, 2020 a r X i v : . [ q - b i o . P E ] J u l his virus spreads rapidly and efficiently: by mid-June the virus had already infected over 6 million people.As of mid-July, over 13 million COVID-19 cases have been confirmed.To tackle and extenuate the effects caused by this virus, global scientific efforts have been called out, beingurgently necessary Bedford et al. (2019). Unfortunately, vaccines take quite some time to be procedure andare, a priori , previewed to be ready only by mid-2021. To retain the diffuse of this disease, most countrieshave chosen to adopt social distancing measures (in different levels and with diverse strategies) since March,2020 Adam (2020). This has been explicitly pointed out as the most pertinent control option for the COVID-19 outbreak, including the cases of countries with large social inequalities, as Brazil Baumgartner et al.(2020). It should be mentioned that the concept underneath social distancing is to impede the saturationof health systems due to large amounts of active COVID-19 infections, which would require treatment atthe same time. In this way, when social distancing policies are enacted, the demands for treatment becomediluted over time, and the health systems do not have to deal with hospital bed shortages associated witha large peak of active infections.Brazil has shown itself as quite a particular case regarding COVID-19 Werneck & Carvalho (2020): thecountry is very large, with 26 federated states, and each federated state has had autonomy to choose its ownhealth policies, which has lead to different levels of social distancing measures for each state. Furthermore,there has been no coordinated nation-wide public health policy to address the viral spread by the federalgovernment, which is very reluctant to do so, claiming that the negative economic effects are too steep andthat social distancing is an erroneous choice The Lancet (2020); Zacchi & Morato (2020). The countrycurrently ranks as second with respect to numbers of cases and deaths. The expectations disclosed on therecent literature suggest catastrophic scenarios for the next few months Rocha Filho et al. (2020); Moratoet al. (2020), which might pursue until mid-2021.Regarding this discussion, we consider the COVID-19 contagion data from Brazil in this paper. In orderto better illustrate the Brazilian scenario with respect to other countries, Figure 1 depicts the evolutionof the SARS-CoV-2 contagion curves in different countries of the world, considering the cumulative casesconfirmed according to the pandemic period of each country (1a). We note that, as of July 7th, 2020, Brazilcounts over 1 . (a) The COVID-19 contagion curve worldwide (in days,since the 100th confirmed case). (b) World map COVID-19 case concentration level. Figure 1: COVID-19 in the world. Data published by the European Center for Disease Prevention and Control (ECDC),Max Roser & Hasell (2020).
Brazil is currently facing many issues due to the SARS-CoV-2 contagion, despite having a strong univer-sal health system. The current situation is near-collapsing, since the majority of Intense Care Unit (ICU)hospital beds are occupied with COVID-19 patients, all around the country. In addition, the virus is pro-2 igure 2: Brazilian map w.r.t. to “hot-spots” of COVID-19 cases, from Google, dating 30 / / gressing to farthest western cities of the country, away from urban areas, where medical care is somehow lesspresent. It can also be noted that SARS-CoV-2 is currently posing great threat to indigenous communities,such as the Yanomami and Ye’kwana ethnicities . Overall, the current situation in the country is verycritical.The first official death due to the SARS-CoV-2 virus in Brazil was registered in March, 17th 2020.Anyhow, through inferential statistics, Delatorre et al. (2020) acknowledge the fact that community trans-mission has been ongoing in the state of S˜ao Paulo since the beginning of February (over one month beforethe first official reports). This points to empirical evidence that the true amount of infected individuals, andpossibly registered deaths, are actually very under-reported. Moreover, due to the absence of mass testing,the country is only accounting for COVID-19 patients with moderate to severe symptoms. People withmild or no symptoms are not being accounted for, following guidelines from the federal Ministry of Health.Additionally, we know that the SARS-CoV-2 virus is, for a large number of individuals, an asymptomaticdisease (yet transmissible). For these reasons, the scientific community has been warning for a possibly hugemargins of underestimated cases in Brazil Silva et al. (2020); Rocha Filho et al. (2020); Delatorre et al.(2020); Paixo et al. (2020); Bastos et al. (2020).Regarding this issue, it becomes fundamental to perform coordinated social distancing interventions, assoon as possible. Furthermore, these public health interventions should be put in practice at the correct timeand for the correct duration. Well-designed social distance guidelines should help mitigating the contagionand thus avoiding the saturation of the heath systems, while, at the same time, trying to balance socialand economic side-effects by releasing/relaxing the quarantine measures as soon as possible. Modeling andidentifying the COVID-19 epidemiological process, and designing an optimal controller that addresses thisissue are the main focuses of this paper.Regarding model frameworks for this contagion, recent literature has demonstrated how the SARS-CoV-2 viral contagion dynamics can be very appropriately described by Susceptible-Infected-Recovered-Deceased The Brazilian Socioenvironmental Institute (ISA,
Instituto Socioambiental , see ) hasreleased a technical note Ferrante & Fearnside (2020) which warns for the contagion of COVID-19 of up to 40 % of YanomamiIndigenous Lands, amid the states of Amazonas and Roraima and a long the border between Brazil and Venezuela, due tothe presence of approximately 20000 illegal mining prospectors. Datasets regarding the COVID-19 spread amid indigenouscommunities are available in https://covid19.socioambiental.org . • Alleman et al. (2020) consider the application of a NMPC algorithm regarding the data from Belgium.The control input is the actual isolation parameter that is plugged into the SIRD model; • Morato et al. (2020) consider an optimal On-Off MPC design, which is formulated as mixed-integerproblem with dwell-time constraints. Furthermore, the response of the population to social isolationrules with an additional dynamic variable, which is then incorporated to the SIRD dynamics; • K¨ohler et al. (2020) applied the technique for the German scenario. The approach considers that thecontrol input is a variable that directly affects the infection and transmission rates (parameters of theSIRD model).Building upon this previous results, herein we propose an NMPC algorithm formulated on the basis ofa SIRD-kind model, which is obtained by the means of an identification algorithm. Following the lines ofMorato et al. (2020), we assumed that the control action is a social isolation guideline, which is enactedand passed on to the population. Then, the population responds in some time to these measures, whichcan be mathematically described as a dynamic social isolation factor. Furthermore, in accordance with thediscussion presented by K¨ohler et al. (2020), we also introduce an additional dynamic nonlinear model, whichgives the input/output (I/O) relationship between the COVID-19 contagion infection and transmission ratesaccording to the stage of the pandemic, as suggested by Bastos et al. (2020). Complementary, we use auto-regressive equations for the epidemiological parameters of the disease (transmission factor, infection factor,lethality rate) directly related to social distancing as control input. These time-varying regressions are usedto provide better long-term forecasts than the regular SIRD equations. Regarding the NMPC framework,the novelty in this paper resides in the formulation of the control input as a finitely parametrized variable,which makes the NMPC implementation persuadable through search mechanisms, which run much fasterthan the Nonlinear Programming methods from the prior.Motivated by the previous discussion, the problem of how optimal predictive control can be used toformulate adequate social distancing policies, regarding Brazil, is investigated in this work. The maincontribution of our approach comprises the following ingredients: • Firstly, an adapted SIRD model is proposed, which incorporates delayed auto-regressive dynamics forthe transmission, infection and lethality rates of the SARS-CoV-2 virus, which vary according to thestage of the pandemic. The adapted model also incorporates the dynamic response of the populationto a given social distancing guideline, which is the NMPC control input. • In order to identify the SIRD model parameter auto-regressive dynamics (for the epidemiologicalparameters), we use a three-layered identification procedure, which concatenates analytical expansionsand Least-Squares optimization procedures. The identified models are validated with regard to theavailable set of data from Brazil, which includes the social distancing factor observed in the country. • The main innovation of the paper is the proposed Parametrized Nonlinear Predictive Control algorithm,which is based on the identified adapted SIRD model with auto-regressive epidemiological parameters.The control input (social distancing guideline) is finitely parametrized over the NMPC predictionhorizon, at each sampling instant. Then, an explicit nonlinear programming solver is simulated forall possible input sequences along time. The NMPC solution, then, is simply found through a searchmechanism regarding these simulated sequences, which is rather numerically cost-efficient.4
Finally, we present results considering the application of this NMPC algorithm to the COVID-19scenario in Brazil. These results are a twofold: a) those that regard the application of the optimalstrategy to control the pandemic since its beginning , comparing the simulation results to those seen inpractice; and b) applying the control method from 30th of July onward, aiming to mitigate and revertthe current health crisis catastrophe.This paper is structured as follows. Section 2 presents the new SIRD model with auto-regressive epidemi-ological parameters, the respective identification procedure and validation results. Section 3 discusses theproposed NMPC strategy. Section 4 depicts the obtained control results, regarding the COVID-19 contagionmitigation for Brazil. General conclusions are drawn in Section 5.
2. Model, Identification and Validation
In this Section, we describe the used modeling framework for the COVID-19 contagion dynamics in Brazil.We build upon the SIRD models for the SARS-CoV-2 virus from works of Peng et al. (2020); Kucharski et al.(2020); Ndairou et al. (2020). The model adaptations are done such that the epidemiological parametersare taken as time-varying piece-wise constant functions of the stage of the pandemic. This adaptation is inaccordance with recent immunology results Sun et al. (2020); Dowd et al. (2020); He et al. (2020), whichdiscuss the transmission and reproduction rate of this virus. It is also taken into account two differentsampling periods, to include the issue of the incubation period of the virus in human, reportedly, in average,as 5 . The SIRD model describes the spread of a given disease with respect to a population split into fournon-intersecting classes, which stand for: • The total amount of susceptible individuals, that are prone to contract the disease at a given (discrete)sample of time k , denoted through the dynamic variable S ( k ); • The individuals that are currently infected with the disease (active infections at a given sample of time k ), denoted through the dynamic variable I ( k ); • The total amount of recovered individuals, that have already recovered from the disease, from an initialinstant k until the current sample k , denoted through the dynamic variable R ( k ); • Finally, the total amount of deceased individuals due to a fatal SARS-CoV-2 infection, from an initialinstant k until the current sample k , is denoted D ( k ).Due to the evolution of the spread of the disease, the size of each of these classes change over time.Therefore, the total population size N is the sum of the first three classes as follows: N ( k ) = S ( k ) + I ( k ) + R ( k ) . (1)Since the Brazilian government discloses daily samples of total infections and accumulated deaths, weconsider that these discrete-time dynamics samples k , given each T = 1 day. Furthermore, to account forthe average incubation period of the disease, we consider that the epidemiological parameters vary weekly(each T = 7 days), kept as zero-order-held/piece-wise-constant samples. We will further assess this matterin the sequel.In the SIRD model, there are three major epidemiological parameters, which express the specific dynamicsof the SARS-CoV-2 virus in the population set. The dynamics of these parameters are given in the sparserdiscrete-time weekly samples , since they change according to the incubation of the virus:5 The transmission rate parameter β stands for the average number of contacts that are sufficient fortransmission of the virus from one individual. According to the detailed classes of individuals, then, itfollows that T β ( k ) I ( k ) /N ( k ) determines the number of contacts that are sufficient for transmissionfrom infected individuals to one susceptible individual; and ( T β ( k ) I ( k ) /N ( k )) S ( k ) determines thenumber of new cases (per day) due to the amount of S ( k ) susceptible individuals (these are the ones“available for infection”). • The infectiousness Poisson parameter γ denotes the inverse of the period of time a given infectedindividual is indeed infectious. Consequently, γ affects the rate of recovery (or death) of an infectedperson. This parameter directly quantifies the amount of individuals that “leaves” the infected class,in a given sample. • The mortality rate parameter ρ stands for the observed mortality rate of the COVID-19 contagion.We model the amount of deceased individuals due to the SARS-CoV-2 infection following the lines of(Keeling & Rohani, 2011): the new number of deaths, at each day, can be accounted for through thefollowing expressions: T ρ ( k )1 − ρ ( k ) γ ( k ) I ( k )Considering these explanations, the “ SIRD ” (Susceptible-Infected-Recovered-Dead) model is expressedthrough the following nonlinear discrete-time difference equations: S ( k + 1) = S ( k ) − T (1 − ψ ( k )) β ( k ) I ( k ) S ( k ) N ( k ) I ( k + 1) = I ( k ) + T (1 − ψ ( k )) β ( k ) I ( k ) S ( k ) N ( k ) − T γ ( k ) I ( k )1 − ρ ( k ) I S ( k + 1) = p sym I ( k + 1) I c ( k + 1) = ( I ( k + 1) + R ( k + 1) + D ( k + 1)) R ( k + 1) = R ( k ) + T γ ( k ) I ( k ) D ( k + 1) = D ( k ) + T ρ ( k )1 − ρ ( k ) γ ( k ) I S ( k ) [SIRD] , (2)where I S denotes the portion of the infected individuals which in fact display symptoms. This “symptomatic”class has been accounted for in previous papers Bastos & Cajueiro (2020); Morato et al. (2020). We notethat only these symptomatic individuals will require possible hospitalization and may die. The remainder(1 − p sym ) I are asymptomatic or lightly-symptomatic, which do transmit the virus but do not die or requirehospitalization. The class of recovered individuals considers immunized people, encompassing those thatdisplayed symptoms and those that did not. The symptomatic ration parameter p sym is constant andborrowed from previous papers.The cumulative number of cases I c denotes the total number of people that have been infected by theSARS-CoV-2 virus until a given day. In fact, this cumulative variable is equivalent to those that are currentlyinfected I , summed with those that have recovered R and those that have died D .We must also stress that, in the SIRD model equations given above, parameter ψ represents a transmissionrate mitigation factor. This factor accounts for the observed social isolation factor with the population set N .It follows that ψ = 0 denotes the situation where the whole population set has sustained social interactions.As discussed in previous papers (Keeling & Rohani, 2011; Bastos et al., 2020), there exists some ”natural” ψ = ψ factor, which stands for a normal/nominal conditions, still with “no control” of the viral spread(with no social isolation guidelines, this kind of situation has been observed in Brazil in the first weeks ofthe contagion, end of February, beginning of March). In contrast, ψ = 1 represents a complete lockdownscenario, where there are no more social interactions (this conditions has not been seen and is, in practice,unattainable). In practice, there is some maximal social isolation factor ψ = ψ that can be put in practice.In this paper, we consider the values for social isolation in Brazil from the work of Bastos et al. (2020),which discloses weekly estimates for ψ . 6nother essential information in epidemiology theory is the basic reproduction number, usually denotedby R . This factor is able to measure the average potential transmissibility of the disease. In practice, thisvalue is fixed/constant and inherent to a given disease. Through the sequel, a time-varying version of thisbasic reproduction number is considered, denoted R t , which represents how many COVID-19 cases couldbe expected to be generated due to a single primary case, in a population for which all individuals aresusceptible. From a control viewpoint, R t represents the epidemic spread velocity: if R t >
1, the infection isspreading and the number of infected people increases along time (this typically happens at the beginningof the epidemic), otherwise, if R t <
1, it means that more individuals “leave” from the infected class, eitherrecovering or dying, and, thereby, the epidemic ceases. The reproduction number R is affected by differentfactors, including immunology of the virus, biological characteristics, but also governments policies to controlthe number of susceptible population. Since the epidemic parameters change along the time, we process byconsidering the effective reproduction number R also as a dynamic variable.The underlying assumption used to calculate R t , is that, at the beginning of the pandemic, S ≈ N .Considering parameters β , γ , ρ and ψ from the SIRD model Eq. (2), R t can be roughly given by: R t ( k ) ≈ (1 − ψ ( k )) β ( k )(1 − ρ ( k )) γ ( k ) . (3) Remark 1.
Regarding the SIRD model given in Eq. (2), it follows that N (0) = N is the initial populationsize (prior to the viral infection). Furthermore, we stress that, in SIRD-kind models, I ( k ) represents the active infections at a given moment, while D ( k ) represents the total amount of deaths until this givenmoment; for this reason, it follows that D ( k + 1) − D ( k ) is proportionally dependent to I ( k ). Remark 2.
We note that we are not able to use more “complex” descriptions of the COVID-19 contagionin Brazil, such as the “SIDARTHE” model used by K¨ohler et al. (2020) (that splits the infections into(symptomatic, asymptomatic) detected, undetected, recovered, threatened and extinct) because we haveinsufficient amount of data. The Ministry of Health only disclosed the total amount of infections ( I ( k ) + D ( k ) + R ( k )) and the total amount of deaths ( D ( k )) at each day. Due to the absence of mass (sampled)testing, there is no data regarding detected asymptomatic individuals, for instance, as it is available inGermany (where (K¨ohler et al., 2020) originate). To choose a more complex model may only decreasethe truthfulness/validity of the identification results, since the parameters can only represent a singularcombination that matches the identification datasets and cannot be used for forecasting/prediction purposes.To illustrate the dynamics of a contagion like the pandemic outbreak of COVID-19, in Figure 3 wepresent long-term forecast using a regular SIRD model with constant epidemiological parameters, followingthe methodology presented by Bastos et al. (2020). These predictions were computed on 11 / / I and D curves. The active infections curve I shows a increase-peak-decrease characteristic, while the total number of deaths D shows an asymptoticbehavior to some steady-state condition. We note that, as of 30 /
06, the number of infections had alreadysignificantly increased, which means that the most recent forecasts preview an even worse catastrophe.
In this paper, the SIRD dynamics from Eq. (2) are adapted. These adaptations are included for threemain reasons:1. In accordance with the discussions presented by Morato et al. (2020), we understand that it is im-perative to embed to the SIRD model how the population responds to public health policies, that is,to incorporate the dynamics and effects of social distancing. When a government enacts some socialdistancing policy, the population takes some time to adapt to it and to, in fact, exhibit the expectedsocial distancing factor ψ . This is, to take ψ as a time-varying dynamic map of the enacted socialdistancing guidelines (which will, later on, be defined by an optimal controller).7 ctive Infections Apr 2020 Jul 2020 Oct 2020 Jan 2021 Apr 2021 Jul 2021 Oct 2021 Jan 20220246810121416 I ( t ) Total Deaths
Jul 2020 Jan 2021 Jul 2021 Jan 202200.20.40.60.811.21.41.61.82 D ( t ) Mar 25 Apr 08 Apr 22 May 06 May 20 Jun 032020 00.511.522.53 10 Apr 28 May 12 May 26 Jun 092020 05 10 Figure 3: Long-term SIRD Forecasts for Brazil (from 11 / / I ( t ) and D ( t )): a) No-ControlSituation ( − ), b) Hard social distancing ( ψ = 0 .
6) from 11 / / − ), c) ψ = 0 . − ) and d) Real data ( × ). Shades represent total variation over a 95% confidence interval, solid lines representmean values.
2. Many results seen in the literature Bastos & Cajueiro (2020); Bhatia et al. (2020); Ndairou et al.(2020) consider constant parameter values for γ , β and ρ , understanding that these factor are inherentto the disease. Anyhow, recent results Bastos et al. (2020); Kucharski et al. (2020) indicate that theseparameters are indeed time-varying, achieve steady-state values at ending stages of the pandemic.It is reasonable that when more active infections occur, the virus tends to spread more efficiently,for instance. Moreover, the mortality rate should stabilize after the infection curve has decreased.To incorporate this issue in the model, auto-regressive moving-average time-varying dynamics for theviral parameters are proposed, which stabilize according to stage of the pandemic. Furthermore,by doing so, the inherent incubation characteristic of the virus is also taken into account, since theviral transmission, lethality and infection rates should vary with dynamics coherent with the averageincubation period.3. Recent papers provide forecasts Bastos & Cajueiro (2020); Bhatia et al. (2020) for the evolution of thevirus assuming that the epidemiological parameters remain constant along the prediction horizon. Thisis not reasonable, due to the fact that such parameters are inherently time-varying, given accordingto the exhibited pandemic level. Thus, if parameters are held constant for forecasts, these forecastsare only qualitative and allow short-term conclusions. Since dynamic models for the epidemiologicalparameters are proposed, make more coherent long-term forecasts can be derived. Of course, it must bestressed that these forecasts will still be qualitative, since one cannot account for perfect accuratenessregarding 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 ofan effective vaccine, which would certainly make the infections drop largely.These model adaptations are discussed individually in the sequel. In order to design and synthesize effective control strategies for social distancing (public) policies, to beoriented to the population by local governments, the social distancing factor ψ is further exploited.8n what concerns the available data from Brazil (that is used for identification of the model parameters),the social distancing factor ψ is a known variable, given in weekly piece-wise constant samples. Later on,regarding the control procedure (Section 4), this isolation factor will vary according to the enacted socialdistancing policy u , as defined by a nonlinear optimal predictive control algorithm.The differential equation that models the response of the susceptible population to quarantine guidelinesis taken as suggests Morato et al. (2020), this is: (cid:8) ψ ( k + 1) = ψ ( k ) + T (cid:37) ψ ( K ψ ( k ) u ( k ) − ψ ( k )) , (4)where u ( k ) is the actual control input, the guideline that defines the social isolation factor goal, as regulatedby the government (this signal will be later on determined by the proposed optimal controller), (cid:37) ψ =0 . − is a settling-time parameter, which is related to the average time the population takes torespond to the enacted social isolation measures, and K ψ is a time-varying stating gain relationship between ψ and u .As recommend Morato et al. (2020), we assume that when more infected cases have been reported, andwhen the hospital bed occupation surpasses 70%, the population will be more prone to follow the socialdistancing guidelines, with larger values for the gain relationship K ψ : K ψ ( k ) = max (cid:26) p sym I ( k )0 . n ICU (cid:27) , (5)being n ICU the total number of ICU beds in the country. We recall that p sym is a parameter which gives theamount of infected individuals that in fact display symptoms and may need to be hospitalized ( I s = p sym I ). Remark 3.
We note that Brazil has, as of February, 45 ,
848 ICU beds. This number is estimated toincrease in up to 80 % with field hospitals that were built specifically for the COVID-19 contagion (i.e. n ICU = 82 , Remark 4.
In practice, Eq. (4) is bounded to the minimal and maximal values for the social distancingfactor ψ and ψ , respectively. These values are the same are that used as saturation constraints on u , asdiscussed in the sequel.We proceed by considering that u is a finitely parametrized control input. This is: the enacted socialdistancing guidelines can only be given within a set of pre-defined values. This approach is coherent withpossible ways to enforce and put in practice social distancing measures. Pursuing this matter, the actualcontrol signal, that is to be defined by the proposed optimal automatic controller, must abide to a piece-wise constant characteristic with a possible increase/decrease of 5% of isolation per week. This percentageis the average increase of social isolation, in the beginning of the pandemic in Brazil, when social distancingmeasures were strengthened significantly over a small time period. The values for social isolation ψ ( k ) inBrazil have been estimated by InLoco (2020) and are presented in Figure 4.Moreover, u is considered to be defined within the admissibility set [ ψ , ψ ]. As gives Figure 4, which showsthe observed social isolation factor ψ in Brazil, the “natural” isolation factor is of ψ = 0 . ψ is the maximal attainable isolation factor forthe country; the maximal observed value in Brazil is of roughly 53 % (see Figure 4). Anyhow, coordinated“lockdown” measures were not forcefully enacted. For this matter, we consider ψ = 0 .
7, which wouldrepresent that the population can be restricted to, at most, a 70% reduction on the level of social interactions.As reported by Bastos et al. (2020), it seems unreasonable to consider larger values for social isolation inthe country, due to multiple reasons (hunger, social inequalities, laboring necessities, lack of financial aidfrom the government for people to stay home, etc. ). 9 ar Apr May Jun
Time [days] S o c i a l D i s t an c i ng [ - ] Figure 4: Social distancing factor in Brazil from February 26th, 2020 to June, 30th 2020 InLoco (2020).
Remark 5.
It must noted that the SIRD model considers some “natural” isolation factor (which is thelower bound on u ). As seen in Figure 4, the social isolation factor of 30 % is the observed lower bound for ψ seen in Brazil.Mathematically, these constraints are expressed through the following Equation set: u ( k ) = u ( k −
1) + δu ( k ) δu ( k ) = {− .
05 or 0 or 0 . } u ( k ) ∈ (cid:2) .
30 0 .
35 0 . . . . .
60 0 .
65 0 . (cid:3) . (6)This constraints given in Eq. (6) are, in fact, very interesting from an implementation viewpoint, asthe government could convert the finitely parametrized values for u into actual practicable measures, asillustrates Table 1. We remark that this Table is only illustrative. Epidemiologists and viral specialistsshould be the ones to formally discuss the actual implemented measures that ensure the social distancingfactor guideline given by u . The second main modification of the original SIRD model is to consider dynamic models for the epi-demiological parameters of the SARS-CoV-2 virus, γ , β and ρ . These models are taken as auto-regressivemoving average functions, which converge as the pandemic progresses. Therefore, the following dynamicsare considered: β ( k ) = f β ( β ( k − . . . , β ( k − n β )) γ ( k ) = f γ ( γ ( k − . . . , γ ( k − n γ )) ρ ( k ) = f ρ ( ρ ( k − . . . , ρ ( k − n ρ )) . (7)The models given in Eq. (7) are possibly delayed and auto-regressive. Anyhow, despite the parameters β , γ and ρ being time-varying, the model functions f β , f γ and f ρ are constant. The order and number ofregressions are found through an optimization procedure, which is further detailed in Section 2.3.10 able 1: Illustrative Example of Finitely Parametrized Social Distancing Measures. Control Signal / Social Distancing Guideline ( u ) Implemented measures Infection Risk u = 0 . u = 0 .
35 Low restriction levels. Use ofmasks to go outside. Publictransport functioning. Lim-ited opening of shops andsmall public spaces. Low Risk... ... u = 0 . u = 0 .
45 Very restrictive policies.Reduced public transport.Urge for people to stay homeat all times. Very restrictiveopenings only. High Risk... ... u = 0 . .2.3. The Complete COVID-19 Model The complete model used in this work to describe the COVID-19 contagion outbreak in Brazil is illus-trated through the block-diagram of Figure 5. We note that the “COVID-19 Epidemiological ParameterModels”, represented through Eq. (7), enable the complete model to offer long-term predictions with moreaccurateness, as discussed in the beginning of this Section. We denote as the “SIRD+ARIMA” model,henceforth, the cascade of the population response from Eq. (4) to the SIRD model, with auto-regressivetime-varying epidemiological parameters as give Eq. (7).
Figure 5: The “SIRD+ARIMA” model for COVID-19 in Brazil.
Recent numerical algorithms have been applied to estimate the model parameters of the COVID-19pandemic Bastos et al. (2020); Sarkar et al. (2020); Sun & Wang (2020); Kennedy et al. (2020); Caccavo(2020); Oliveira et al. (2020); Fernndez-Villaverde & Jones (2020). The majority of these papers indicatesthat SIRD models provide the best model-data fits.We note, anyhow, that the SIRD model offers three degrees-of-freedom at each instant k (i.e. β ( k ), γ ( k )and ρ ( k )), as gives Equation 2. This means that different instantaneous combinations of these parameterscan yield the same numerical values for S ( k ), I ( k ), R ( k ) and D ( k ). Therefore, although mathematical andgraphical criteria have been used to validate these dynamic models when compared to real data, the estimatedvalues for these parameters should be coherent with biological characteristics of this viral pandemic. Anindicative of badly adjusted SIRD model parameters is the effective reproduction number of the disease R ,which should naturally decrease to a threshold smaller than one as the pandemic ceases.Thus, in order to identify a SIRD model that is indeed able to describe the pandemic contagion accurately(which is especially important for forecast and control purposes), we proposed a two step identificationmethod to determine the time-varying dynamics for β ( k ), γ ( k ) and ρ ( k ), with respect to the dynamicparameter models in Eq. (7).Regarding datasets, we must comment that the Brazilian Ministry of Health discloses, daily, the valuesfor the cumulative number of confirmed cases I c ( k ) and the cumulative number of deaths D ( k ). Furthermore,we recall that, through the sequel, the values of the observed average social isolation in the country, ψ ( k ),are assumed to be known (as given in Figure 4). The symptomatic percentage is assumed as constant alongthe prediction horizons.Then, we proceed by depicting the used identification procedure, which is performed in three consecutivesteps/layers, as follows: • The first step resides in analytically solving the SIRD regressions from Equation 2 for a fixed intervalof points (days). • Then, the derived values from the analytical solutions are passed as initial conditions to an optimizationlayer, which solves a constrained Ordinary Least Square minimization problem, aiming to adjust the12arameter values to pre-defined (biologically coherent) sets, so that the identified model yields smallerror with adequate parameter choices. The output of this second step stands for the time-seriesvectors regarding the SIRD parameters. • Finally, these time-series are used to fit auto-regressive models, via Least Squares, as gives Eq. (7).By following this three-layered procedure, the estimated parameter equations (Eq. (7)) are found inaccordance with biological conditions, which are described as the pre-defined optimization sets. We notethat we embed to the optimization layer sets that are also in accordance with previous results for SIRDmodel estimations for Brazil Bastos & Cajueiro (2020); Morato et al. (2020).
Figure 6: Parameter Identification Procedure.
The complete identification procedure used to estimate the SIRD model dynamics follows the linesillustrated in Figure 6.
We note that, in order to simplifying the formulation of the optimization algorithm, the differenceEquations for I ( k + 1) and D ( k + 1) are modified in order to decouple parameters related to the number offatal and related to the number of recovered individuals, with regard to I ( k ). This is: I ( k + 1) = I ( k ) + T (1 − ψ ( k )) β ( k ) I ( k ) S ( k ) N ( k ) − T γ ( k ) I ( k ) − T α ( k ) I ( k ) (8) D ( k + 1) = D ( k ) + T α ( k ) I ( k ) (9)where α ( k ) = ρ ( k ) γ ( k )(1 − ρ ( k )) . Therefore, instead of identifying β , γ and ρ , we pursue the identification of β , γ and α , and, then, ρ is computed as follows: ρ ( k ) = α ( k ) γ ( k ) α ( k ) γ ( k ) . (10)The identification procedure starts by collecting the available datasets (regarding the cumulative numberof COVID-19 cases I c and deaths D ) inside a fixed interval k ∈ [ k i , k f ]. Then, the first layer computes13exact”, analytical values for the epidemiological parameters, denoted ˜ β , ˜ γ and ˜ α according to the followingdiscrete analytical expansions:˜ γ = R ( k f ) − R ( k i ) (cid:80) i = k f i = k i I ( i ) , (11)˜ α = D ( k f ) − D ( k i ) (cid:80) i = k f i = k i I ( i ) , (12)˜ β = 1(1 − ψ ( k f )) I ( k f ) − I ( k i ) + (˜ α + ˜ γ ) (cid:80) i = k f i = k i I ( i ) (cid:80) i = k f i = k i S ( i ) I ( i ) /N ( i ) . (13) It is assumed that the used datasets might be corrupted by series of issues, such as cases that are notreported on given day k and simply accounted for on the following days, or sub-reported cases, as discussedby Bastos et al. (2020). Therefore, we adjust these parameters through the optimization layer, in order toimprove the reliability of the identification.Once again, the available data from the same interval discrete [ k i , k f ] is used. The optimization procedureminimizes a constrained Ordinary Least Squares problem, whose solution comprises the following vectors: S , I , R and D . These vectors comprise the model-based outputs within the considered discrete interval.The decision variables for the optimization procedure are the degrees-of-freedom of the SIRD model, in thislayer denoted ˆ β , ˆ γ , ˆ α . The Ordinary Least Square criterion ensures that the optimization minimizes theerror between the real data and the estimated values from the SIRD model, by choosing coherent values forthe decision variables.The quadratic model-data error variables terms used in the optimization layer are the following: Er I ( i ) = ( I ( i ) − ˆ I ( i, ˆ β (1 − ψ ) , ˆ γ, ˆ α )) , (14) Er R ( i ) = ( R ( i ) − ˆ R ( i, ˆ β (1 − ψ ) , ˆ γ, ˆ α )) , (15) Er D ( i ) = ( D ( i ) − ˆ D ( i, ˆ β (1 − ψ ) , ˆ γ, ˆ α )) , (16)for which the variables ˆ I, ˆ R, ˆ D are estimated according to the SIRD model equations. The complete opti-mization problem is formulated as follows:min β,γ,α J = i = k f (cid:88) i = k i ( w Er I ( i ) + w Er R ( i ) + w Er D ( i )) ,s.t.: δ ˜ β (1 − ψ ( i )) ≤ β (1 − ψ ( i )) ≤ δ ˜ β (1 − ψ ( i )), (17) δ ˜ γ ≤ γ ≤ δ ˜ γ , δ ˜ α ≤ α ≤ δ ˜ α . β γ α = ˜ β ˜ γ ˜ α .wherein δ and δ are uncertainty interval margins used to define the lower and upper bound of each decisionvariable on the optimization problem, β , γ and α are the initial conditions of the optimization problem and w , w and w are taken as positive weighting values (tuning parameters), used to normalize the magnitudeorder of the total cost and adjust the model fit with respect to variables Er I , Er R and Er D .With respect to the discussion regarding the average incubation period of the SARS-CoV-2 virus, theoptimization procedure compute piece-wise constant values for β ( k ), γ ( k ) and α ( k ) each T = 7 days. Thisis also done to improve the model-data fitting results, as suggest Morato et al. (2020).14onsidering the application of this second layer to a complete data-set, with new values for the epidemi-ological parameter each T days, the output of this layer stands for the epidemiological time-series vectorsdenoted β opt , γ opt and α opt . These time-series are, then, used to fit the auto-regressive models in Eq. (7),in the third layer. We note that these SIRD parameter dynamics are the ones that can be can be used forforecasting and control purposes (and also to calculate the effective reproduction number R t ). The third layer of the identification procedure resides in fitting an auto-regressive model to the time-series derived from the optimization β opt , γ opt and ρ opt (note that ρ opt is given as a function of β opt , γ opt and α opt , as gave Eq. (10)). In this paper, the auto-regressive functions in Eq. (7) are of Auto-RegressiveIntegrated Moving Average (ARIMA) kind. The ARIMA approach is widely used for prediction of epidemicdiseases, as shown by Kirbas et al. (2020). It follows that, from a time-series viewpoint, the ARIMA modelcan express the evolving of a given variable (in this case, the epidemiological parameters) based on priorvalues. Such models, then, are coherent with the prequel discussion regarding the convergence of theseparameters to steady-state conditions (Sec. 2.2).For simplicity, instead of presenting the ARIMA fits for the three epidemiological time-series ( β opt , γ opt and ρ opt ), we proceed by focusing on the SARS-CoV-2 transmission factor β . We note that equivalentsteps are pursued for the other parameters. The main purpose of this third layer is to model the trends ofthe SIRD epidemiological parameters (as provided by the two previous layers) and use these trends, in thefashion of Eq. (7), to improve the forecasting of the SIRD+ARIMA model, making it more coherent andconsistent for feedback control strategies.It is worth mentioning that this layer is an innovative and important advantage of the SIRD modelidentification proposed in this work. As depicted by Fernndez-Villaverde & Jones (2020), using a SIRDmodel with time-varying epidemiological parameters allows one to provide forecasts that differ at the initialand long-term moments of the COVID-19 contagion. β As exploited in Section 2, the transmission rate parameter β gives an important measure to analyze thepandemic panorama. It has been shown that this parameter varies according to health measures appliedto the prone population Fernndez-Villaverde & Jones (2020); Oliveira et al. (2020); Caccavo (2020). TheARIMA expression is given as follows: (cid:16) a β + a β β ( k −
1) + a β β ( k −
2) + · · · + a β nβ β ( k − n β + (cid:15) ( k ) (cid:17) δβ k = (18) (cid:16) b β (cid:15) ( k ) + b β (cid:15) ( k −
1) + b β (cid:15) ( k −
2) + · · · + b β nβ (cid:15) ( k − n β ) (cid:17) ,where δβ k = [ δβ ( k ) δβ ( k − δβ ( k − . . . δβ ( k − n β )] T is the vector of each incremental difference δβ ( k ) = β ( k ) − β ( k − a β , . . . , b β nβ ). This procedure is based onthe following steps:(i) Concatenate the ARIMA the regression term from Eq. (18) as β ( k ) = ω T ( k − (cid:15) ( k ), where ω ( k −
1) concatenates the input values from the previous layers (in this case, the time-series β opt ) andΘ concatenates the ARIMA coefficients;(ii) Determine a constrained Least-Squares estimation ˆΘ, w.r.t. Eq. (18);(iii) Compute the Least-Squares residue (cid:15) = β opt − ω ˆΘ and use it as the initial guess for (cid:15) for the nextiteration;(iv) Iterate until convergence or some (2-norm wise) small residue is achieved.15 .7. Model Validation Results With respect to the detailed identification procedure, the datasets provided by Brazilian Ministry ofHealth are used for validation purposes, considering the interval from the first confirmed COVID-19 case inthe country, dating February 26th, until the data from June 30th, 2020. The Ministry of Health discloseddaily values for the cumulative number of cases and deaths. The data is sufficient to obtain the behaviorsfor S , I , R , and D . Note that the total population size of Brazil is used as the initial condition N ≈ δ = 0 .
95 and δ = 1 . Table 2: Optimization Weights.
Parameter w w w Value 1 10 2We proceed the model validation in a twofold: a) first we consider the first 100 data points to tune theSIRD+ARIMA model and demonstrate its validity, globally well representing the following points; and b)we tune the SIRD+ARIMA model by identifying its parameters for the whole dataset; this is the modelused for control. We note that the values used for the social isolation factor ψ are those as exhibited inFigure 4. Regarding the first validation part, Figure 7 shows the estimated time-series values from the optimizationoutput. Recall that these parameters are piece-wise constant for periods of T = 7 days. As discussed inprevious literature, these parameters tend to follow stationary trends as the pandemic progresses. Mar Apr May Jun
Time [days] [ / d ay ] , , [ / d ay ] (x10) (x10) Figure 7: β , γ , α and ρ parameters identified by the two layer approach - from February 26th, 2020 to June, 04th 2020. In order to illustrate the ARIMA fits, the identified auto-regressive model for SARS-CoV-2 transmissionparameter β is provided in Figure 8. The ELS procedure yields a regression f β ( · ) with n β = 21 dailysteps (i.e. 3 weeks). Figure 8 demonstrates how the ARIMA is indeed able to describe the time-series β opt ,globally catching the time-varying behavior with 99 .
04 % accuracy.16
Time [days] opt
ARIMA model
Figure 8: Comparison between ARIMA model and previous identified β . FIT - 99.04% The first 100 data points comprise the period from February 26th to June, 04th 2020. Using the resultingSIRD+ARIMA model identified in this period, we proceed by extending the model forward, making forecastsfrom June, 05th until June 30th, 2020, in order to demonstrate the validity of the identification. Theforecasts are made using Eqs. (2) wherein parameters γ , β and ρ are given by the ARIMA regressions modelin Eq. (7). Figures 9, 10 and 11 show the model forecasts compared with real data for active infections,recovered individuals and deaths, respectively, considering a interval of confidence of 95%. Clearly, theglobal behavior of the COVID-19 contagion is well described by the proposed models. To further illustratethis fact, Figure 12 depicts the model-data error terms (from Eqs. (14)-(16)), given in percentage. Thecoefficient of determination for these identification results are all above 0 . + ARIMA Model
We note that the data-driven model, that is used for the NMPC strategy, is the based on the whole setof data. The simulation of this model is given in Figure 13.Based on the prior developments, given the high coefficient of determination and observing Figures 8-12(in which a 95% confidence interval is considered), it can be noticed that the identified model SIRD+ARIMAdescribes very well the observed COVID-19 pandemic behavior in Brazil and, in addition, it can forecastthe SIRD curve with relatively small error. It is worth mentioning that this model identification approachneeds to be updated regularly. From Figure 12, one notes that 20-days-ahead forecasts can be derived withmodel-data errors below 10%. This error margins may increase when considering larger forecast horizons.Nevertheless, regarding the control purpose of this work and taking into account prediction horizons ofroughly one month, the proposed SIRD+ARIMA is indeed able to be highly representative regarding theCOVID-19 spread. For perfectly coherent forecasts, it is recommendable for the identification to be re-performed each week, and the ARIMA fits updated. This recommendation has also been provided byMorato et al. (2020).We also remark that the ARIMA fits for the epidemiological parameters are given in weekly samples i.e. k week , while the SIRD variables are given for each day i.e. k day . Regarding this matter, we can represent17 ar Apr May Jun Time [days] A c t i ve I n f ec t e d C ases SIRD+ ModelOfficial data - MS Brazil
June 4thStart Validation
Figure 9: Validation of the SIRD+ARIMA model using estimated parameters with official data - Active Infected Curve. the weekly-sampled variables, from the viewpoint of the daily-sampled variables, as: β ( k day ) = β ( k week T ) ∀ k day ∈ [ k week T , ( k week + 1) T ) , γ ( k day ) = γ ( k week T ) ∀ k day ∈ [ k week T , ( k week + 1) T ) , ρ ( k day ) = ρ ( k week T ) ∀ k day ∈ [ k week T , ( k week + 1) T ) , (19) ψ ( k day ) = ψ ( k week T ) ∀ k day ∈ [ k week T , ( k week + 1) T ) , u ( k day ) = u ( k week T ) ∀ k day ∈ [ k week T , ( k week + 1) T ) ,where T = 7 days.Finally, considering this validated SIRD+ARIMA model, the basic reproduction number R t of the SARS-CoV-2 virus in Brazil can be inferred through Eq. (3). In Figure 14, the evolution of the viral reproductionnumber R t for each week of the pandemic in show. As it can be seen, this reproduction number presentsstronger variations at the beginning stage of, but then tends to converge to a steadier values, which cor-roborates with the model extensions considering the epidemiological parameters. Moreover, notice thatthe reproduction number for since June 4th, 2020 is somewhat steady near 1 . ar Apr May Jun Time [days] R ec o ve r e d C ases SIRD+ ModelOfficial Data - MS Brazil
June 4thStart Validation
Figure 10: Validation of the SIRD+ARIMA model using estimated parameters with official data - Recovered Curve.
Mar Apr May Jun
Time [days] F a t a l C ases SIRD+ ModelOficial Data - MS Brazil
June 4thStart Validation
Figure 11: Validation of the SIRD+ARIMA model using estimated parameters with official data - Fatal Cases Curve. ar Apr May Jun Time [days] E rr o r D a t a - M od e l [ % ] Infected cases: R = 0.9904Recovered cases: R = 0.9976Fatal cases: R = 0.9954 June 4thStart Validation
Figure 12: Error between the simulated model and official data (in percentage).
Mar Apr May Jun Jul
Time [days] I nd i v i du a l [-] I nd i v i du a l [-] Active Infections Recovered cases Fatal cases (x10) Susceptibles
Figure 13: SIRD model curves used for NMPC framework. Right axis for Susceptible individuals (S). Left axis for ActiveInfections (I), Recovered (R) and Fatal (D) Cases. ar Apr May Jun Time [days] B as i c R e p r odu c t i on N u m b e r ( R t ) June 04thR t = 1.603 Figure 14: Basic reproduction number R t calculated by the multi-layered identification framework. . The NMPC strategy In this Section, we detail the proposed NMPC strategy used to determine public health guidelines(regarding social isolation policies) to mitigate the spread of the COVID-19 contagion in Brazil. ThisNMPC algorithm is conceived under the feedback structure illustrated in Figure 5.We note that the generated control action u ( k ) represents the input to the population’s response to socialdistancing measures, as gives the differential Eq. (4). For this reason, it follows that the proposed NMPCalgorithm operates with a weekly sampling period ( T ). Implementing a new social distancing guidelineevery week is coherent with previous discussions in the literature Morato et al. (2020). It does not seemreasonable to change the distancing measures every few days. Before presenting the actual implementation of NMPC tool, we note that its generated control sequencemust be given in accordance with the constraints expressed in Eq. (6).Considering that the NMPC has a horizon of N p steps (given in weekly samples) from the viewpoint ofeach week k , the generated control signal is : U k = (cid:2) u ( k | k ) u ( k + 1 | k ) · · · u ( k + N p − | k ) (cid:3) . (20)Since the variations from each quarantine guideline u ( k −
1) to the following u ( k ), denoted δu ( k ) areequal to ± .
05 or 0, it follows that all possible control sequences can be described as, from the viewpointof sample k : U k = ( u ( k −
1) + δu ( k )) · · · u ( k − N p times (cid:122) (cid:125)(cid:124) (cid:123) + δu ( k ) · · · + δu ( k + N p − . (21)From Eq. (21), we can conclude that, in the considered settling, there are 3 N p possible control sequences. The main purpose of social isolation is to distribute the demands for hospital bed over time, so thatall infected individuals can be treated. It seems reasonable to act though social isolation to minimize thenumber of active infections ( I ), while ensuring that this class of individuals remains smaller than the totalnumber of available ICU beds n ICU .Moreover, it seems reasonable to ensure that social isolation measures are enacted for as little time aspossible, to mitigate the inherent economic backlashes.This trade-off objective (mitigating I against relaxing social isolation) is expressed mathematically as: || I ( k ) || Qn I + || u ( k ) || (1 − Q ) = I ( k ) T Qn I I ( k ) + u ( k ) T (1 − Q ) u ( k ) , (22)where n I is a nominal limit for I (i.e. the initial population size N ) included for a magnitude normalizationof the trade-off objective. Note that u is given within [0 . . The proposed NMPC algorithm must act to address the Control Objective given in Eq. (22) whileensuring some inherent process constraints. These are:1. Ensure that the control signal has the piece-wise constant characteristic, as gives Eq. (6); and2. Ensuring that the number of infected people with active acute symptoms does not surpass the totalnumber of available ICU hospital beds in Brazil, this is: p sym I ( k ) ≤ n ICU . (23) Notation χ ( k + j | k ) denotes the predicted values for variable χ ( k + j ), computed at the discrete instant k . .4. The Complete NMPC Optimization Bearing in mind the previous discussions, the NMPC procedure is formalized under the following op-timization problem, which is set to mitigate the SARS-CoV-2 viral contagion spread with respect to thetrade-off control objective of Eq. (22) and the constraints from Eqs. (6)-(23), for a control horizon of N p (weekly) steps, from the viewpoint of each sampling instant k :min U k J ( · ) = min U k N p (cid:88) i =1 (cid:18)(cid:18) I ( k + i ) T Qn I I j ( k + i ) (cid:19) + (cid:0) u ( k + i − T (1 − Q ) u ( k + i − (cid:1)(cid:19) , (24)s.t.: SIRD+ARIMA Model ∀ i ∈ N [1 , N p ] , ψ ≤ u ( k + i − ≤ ψ , u ( k + i −
1) = u ( k + i −
2) + δu ( k + i −
1) , δu ( k + i −
1) = − .
05 or 0 or 0 .
05 , p sym I ( k + i ) ≤ n ICU . The finitely parametrized NMPC methodology has been elaborated by Alamir (2006). This paradigmhas recently been extended to multiple applications Rathai et al. (2018, 2019), with successive real-timeresults.Given that this control framework offers the tool to formulate (finitely parametrized) social distancingguidelines for the COVID-19 spread in Brazil, with regard to the discussions in the prior, we proceedby detailing how it is implemented at each sampling instant k , ensuring that the Nonlinear OptimizationProblem in Eq. (24) is solved.Basically, the parametrized NMPC algorithm is implemented by simulating the SIRD+ARIMA modelwith an explicit nonlinear solver, testing it according to all possible control sequences (as gives Eq. (21)).Then, the predicted variables are used to evaluate the cost function J ( · ). The control sequences that implyin the violation of constraints are neglected. Then, the resulting control sequence is the one that yields theminimal J ( · ), while abiding to the constraints. Finally, the first control signal is applied and the horizonslides forward. This paradigm is explained in the Algorithm below. Figure 15 illustrates the flow of theimplementation of the proposed Social Distancing control methodology for the COVID-19 viral spread inBrazil. 23 initely Parametrized NMPC for Social Distancing GuidelinesInitialize : N (0), S (0), I (0), R (0) and D (0). Require : Q , n I , n ICU
Loop every T days : • Step ( i ): ”Measure” the available contagion data ( N ( k ), S ( k ), I ( k ), R ( k ) and D ( k )); • Step ( ii ): Loop every T days : – Step (a):
For each sequence j ∈ N p : ∗ Step (1): Build the control sequence vector according to Eq. (21); ∗ Step (2): Explicitly simulate the future sequence of the SIRD variables; ∗ Step (3): Evaluate if constraints are respected. If they are not, end, else, compute the costfunction J ( · ) value. – end– Step (b): Choose the control sequence U k that corresponds to the smallest J ( · ). – Step (c): Increment k , i.e. k ← k + 1. • end • Step ( iii ): Apply the local control policy u ( k ) as gives Eq. (6); • Step ( iv ): Simulate the SIRD+ARIMA model, as give Eqs. (2)-(4). • Step ( v ): Increment k , i.e. k ← k + 1. end igure 15: Algorithm Implementation Flowchart. . Simulation Results Bearing in mind the previous discussions, the proposed SIRD+ARIMA model and the finitely parametrizedNMPC toolkit, this Section is devoted in presenting the simulation results regarding the COVID-19 conta-gion spread in Brazil. The following results were obtained using Matlab. The average computation timeneeded to solve the solution of the proposed NMPC procedure is of 2 ms . We proceed by depicting two simulation scenarios, which account for the following situations:(a) What would be the case if the NMPC generated social distancing policy was enacted back in thebeginning of the COVID-19 pandemic in Brazil.(b) How plausible is it to apply the NMPC technique from now on (30 / / th sample on, the open-loop simulationtakes ψ = 0 . J ( · ) with tuning weight Q = 0 . N p = 4 weeks (28 days), which means that the controllermakes its decision according to model-based forecasts of the SIRD+ARIMA model for roughly one monthahead of each sample k . For simplicity n I is taken as n ICU p sym , since the main goal of the MPC is to ensure p sym I ≤ n ICU . We recall that the COVID-19 pandemic in Brazil formally “started” 26 / / / / The implementation was performed on a i5 2 . . Figure 16: Scenario (a): Symptomatic Active Infections p sym I and Control Input u (Social Distancing, ψ ). igure 17: Scenario (a): Total Active Infections I and Total Deaths D . The second simulation scenario that we consider is to control the situation from the current moment(30 / / . I c , Figure 20 shows the evolution of all active infections (symptomatic and asymptomatic), whileFigure 21 shows the mounting number of deaths. The proposed NMPC, if rapidly put in practice, could stillbe able to slow the viral spread, saving 25 % of lives w.r.t. the open-loop condition. The peak of infections,if such technique is applied, has its forecast previewed to September 2nd, 2020, being anticipated in 17 days.28 igure 18: Scenario (b): Symptomatic Active Infections p sym I and Control Input u (Social Distancing, ψ ).Figure 19: Scenario (b): Cumulative Cases I c and Recovered Individuals R . igure 20: Scenario (b): Total Active Cases I .Figure 21: Scenario (b): Total Deaths D . . Conclusions In this paper, an optimal control procedure is proposed for ruling social isolation guidelines in Brazil,in order to mitigate the spread of the SARS-CoV-2 virus. To do so, a new model is proposed, based onextensions of the SIRD equations. The proposed model embeds weekly auto-regressive dynamics for theepidemiological parameters and also takes a dynamic social distancing factor (as seen in previous papers(Morato et al., 2020; Bastos et al., 2020)). The social distancing factor measures and expresses the popula-tion’s response to quarantine measures, guided by the Nonlinear Model Predictive Control procedure. TheNMPC strategy is designed within a finitely parametrized input paradigm, which enables its fast realization.In this work, some key insights were given regarding the future panoramas for the COVID-19 pandemicin Brazil. Below, some of the main findings of this paper are summarized: • The presented results corroborate the hypothesis formulated in many of the previous papers regardingthe COVID-19 pandemic in Brazil Hellewell et al. (2020); The Lancet (2020); Morato et al. (2020):herd immunity is not a plausible option for the country ; if no coordinates social distancing action isenforced, the ICU threshold will be largely surpassed, which can lead to elevated fatality. • The prediction of evolution of the viral spread is relatively accurate with the proposed adapted SIRDmodel for up to 20 days. Larger prediction horizons can be considered, but daily model-updates arerecommended. • The simulation forecasts derived with the NMPC strategy and with an open-loop condition (no socialdistancing) indicate that social distancing measures should still be maintained for a long time. Thestrength of these measures will be diluted as time progresses. The forecasts indicate an infection peakof over 600000 symptomatic individuals to late September, 2020, in the current setting. If model-basedcontrol is enacted, the peak could be anticipated and the level of infections could be contained belowthe ICU hospital bed threshold. The NMPC could save over 400000 lives if enacted from now (July,2020). • The results also indicate that if such coordinated control strategy was applied since the first month ofCOVID-19 infections in Brazil, a more relaxed social distancing paradigm would be possible as of late2020. Since this has not been pursued, the social distancing measures may go up until late 2021 if novaccine 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 infectionsis very inconsistent. As Bastos et al. (2020) thoroughly details, the uncertainty margin associated to theavailable data (in terms of case sub-reporting) is very significant. Anyhow, the results presented hereincan help guiding long-term regulatory decision policies in Brazil regarding COVID-19. One must note thatsocial distancing measures, in different levels, will be recurrent and ongoing for a long time. Due to thisfact, compensatory social aid policies should also be developed in order to reduce the effects of a possiblylong-lasting economic turn-down. Recent papers have discussed this matter Ahamed & Guti´errez-Romero(2020); Zacchi & Morato (2020).
Acknowledgments
The Authors acknowledge the financial support of National Council for Scientific and TechnologicalDevelopment (CNPq, Brazil) under grants 304032 / − / − Previous paper have also elaborated on the fact that vertical isolation is also not an option for the time being, since Brazildoes not have the means to pull of efficient public policies to separate the population at risk from those with reduced risk,due to multiple social-economical issues of the country Silva et al. (2020); Rocha Filho et al. (2020); Rodriguez-Morales et al.(2020). otes The authors report no financial disclosure nor any potential conflict of interests.
References
Adam, D. (2020). The simulations driving the worlds response to covid-19. how epidemiologists rushed to model the coronaviruspandemic?
Nature , April .Ahamed, M., & Guti´errez-Romero, R. (2020). The role of financial inclusion in reducing poverty: A covid-19 policy response.Alamir, M. (2006).
Stabilization of nonlinear systems using receding-horizon control schemes: a parametrized approach forfast systems volume 339. Springer.Alleman, T., Torfs, E., & Nopens, I. (2020). COVID-19: from model prediction to model predictive control. https://biomath.ugent. e/sites/default/files/2020-04/Alleman etal v2.pdf , , 2020.Bastos, S. B., & Cajueiro, D. O. (2020). Modeling and forecasting the early evolution of the COVID-19 pandemic in Brazil. arXiv:2003.14288 .Bastos, S. B., Morato, M. M., & anda Julio E Normey-Rico, D. O. C. (2020). The COVID-19 (SARS-CoV-2) uncertaintytripod in Brazil: Assessments on model-based predictions with large under-reporting. arXiv:2006.15268 .Baumgartner, M. T., Lansac-Toha, F. M., Coelho, M. T. P., Dobrovolski, R., & Diniz-Filho, J. A. F. (2020). Social distancingand movement constraint as the most likely factors for COVID-19 outbreak control in brazil. medRxiv , .Bedford, J., Farrar, J., Ihekweazu, C., Kang, G., Koopmans, M., & Nkengasong, J. (2019). A new twenty-first century sciencefor effective epidemic response. Nature , , 130–136.Bhatia, S., Cori, A., Parag, K. V., Mishra, S., Cooper, L. V., Ainslie, K. E. C., Baguelin, M., Bhatt, S., Boonyasiri, A.,Boyd, O., Cattarino, L., Cucunub, Z., Cuomo-Dannenburg, G., Dighe, A., Dorigatti, S., Ilaria van-Elsland, FitzJohn,R., Fu, H., Gaythorpe, K., Green, W., Hamlet, A., Haw, D., Hayes, S., Hinsley, W., Imai, N., Jorgensen, D., Knock,E., Laydon, D., Nedjati-Gilani, G., Okell, L. C., Riley, S., Thompson, H., Unwin, J., Verity, R., Vollmer, M., Walters,C., Wang, H. W., Walker, P. G., Watson, O., Whittaker, C., Wang, Y., Winskill, P., Xi, X., Ghani, A. C., Donnelly,C. A., Ferguson, N. M., & Nouvellet, P. (2020). Short-term forecasts of COVID-19 deaths in multiple countries. URL: https://mrc-ide.github.io/covid19-short-term-forecasts/index.html [Online; Accessed 29-April-2012].Caccavo, D. (2020). Chinese and italian covid-19 outbreaks can be correctly described by a modified sird model. medRxiv , .doi: .Camacho, E. F., & Bordons, C. (2013). Model predictive control . Springer Science & Business Media.Delatorre, E., Mir, D., Graf, T., & Bello, G. (2020). Tracking the onset date of the community spread of SARS-CoV-2 inwestern countries. medRxiv , .Dowd, J. B., Andriano, L., Brazel, D. M., Rotondi, V., Block, P., Ding, X., Liu, Y., & Mills, M. C. (2020). Demographicscience aids in understanding the spread and fatality rates of COVID-19.
Proceedings of the National Academy of Sciences , , 9696–9698.Fernndez-Villaverde, J., & Jones, C. I. (2020). Estimating and Simulating a SIRD Model of COVID-19 for Many Countries,States, and Cities . NBER Working Papers 27128 National Bureau of Economic Research, Inc. URL: https://ideas.repec.org/p/nbr/nberwo/27128.html . doi: .Ferrante, L., & Fearnside, P. M. (2020). Protect indigenous peoples from COVID-19.
Science , , 251–251.He, X., Lau, E. H., Wu, P., Deng, X., Wang, J., Hao, X., Lau, Y. C., Wong, J. Y., Guan, Y., Tan, X. et al. (2020). Temporaldynamics in viral shedding and transmissibility of COVID-19. Nature medicine , , 672–675.Hellewell, J., Abbott, S., Gimma, A., Bosse, N. I., Jarvis, C. I., Russell, T. W., Munday, J. D., Kucharski, A. J., Edmunds,W. J., Sun, F. et al. (2020). Feasibility of controlling COVID-19 outbreaks by isolation of cases and contacts. The LancetGlobal Health , .InLoco (2020). Social isolation map covid-19 (in portuguese). https://mapabrasileirodacovid.inloco.com.br/pt/ .Keeling, M. J., & Rohani, P. (2011).
Modeling Infectious Diseases in Humans and Animals . Princeton University Press.Kennedy, D. M., Zambrano, G. J., Wang, Y., & Neto, O. P. (2020). Modeling the effects of intervention strategies on covid-19transmission dynamics.
Journal of Clinical Virology , , 1–7. doi: https://doi.org/10.1016/j.jcv.2020.104440 .Kermack, W. O., & McKendrick, A. G. (1927). A contribution to the mathematical theory of epidemics. Proceedings of theRoyal Society A , , 700–721.Kirbas, I., Szen, A., Tuncer, A. D., & Kazancolu, F. (2020). Comparative analysis and forecasting of covid-19 cases in variouseuropean countries with arima, narnn and lstm approaches. Chaos Solitons & Fractals , (p. 110015). doi: .K¨ohler, J., Schwenkel, L., Koch, A., Berberich, J., Pauli, P., & Allg¨ower, F. (2020). Robust and optimal predictive control ofthe COVID-19 outbreak. arXiv preprint arXiv:2005.03580 , .Kucharski, A. J., Russell, T. W., Diamond, C., Liu, Y., Edmunds, J., Funk, S., Eggo, R. M., Sun, F., Jit, M., Munday, J. D.et al. (2020). Early dynamics of transmission and control of COVID-19: a mathematical modelling study.
The Lancetinfectious diseases , .Lauer, S. A., Grantz, K. H., Bi, Q., Jones, F. K., Zheng, Q., Meredith, H. R., Azman, A. S., Reich, N. G., & Lessler, J.(2020). The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: estimationand application.
Annals of internal medicine , , 577–582.Lima, C. M. A. d. O. (2020). Information about the new coronavirus disease (COVID-19). Radiologia Brasileira , , V–VI.Max Roser, E. O.-O., Hannah Ritchie, & Hasell, J. (2020). Coronavirus pandemic (covid-19). Our World in Data , .Https://ourworldindata.org/coronavirus. orato, M. M., Bastos, S. B., Cajueiro, D. O., & Normey-Rico, J. E. (2020). An optimal predictive control strategy forCOVID-19 (SARS-CoV-2) social distancing policies in Brazil. arXiv:2005.10797 .Ndairou, F., Area, I., Nieto, J. J., & Torres, D. F. (2020). Mathematical modeling of COVID-19 transmission dynamics witha case study of wuhan. Chaos, Solitons & Fractals , (p. 109846).Oliveira, J. F., Jorge, D. C. P., Veiga, R. V., Rodrigues, M. S., Torquato, M. F., da Silva, N. B., Fiaconne, R. L., Castro, C. P.,Paiva, A. S. S., Cardim, L. L., Amad, A. A. S., Lima, E. A. B. F., Souza, D. S., Pinho, S. T. R., Ramos, P. I. P., Andrade,R. F. S., & (2020). Evaluating the burden of covid-19 on hospital resources in bahia, brazil: A modelling-based analysis of14.8 million individuals. medRxiv , . doi: .Paixo, B., Baroni, L., Salles, R., Escobar, L., de Sousa, C., Pedroso, M., Saldanha, R., Coutinho, R., Porto, F., & Ogasawara,E. (2020). Estimation of COVID-19 under-reporting in brazilian states through SARI. arXiv:2006.12759 .Peng, L., Yang, W., Zhang, D., Zhuge, C., & Hong, L. (2020). Epidemic analysis of COVID-19 in China by dynamical modeling. arXiv preprint arXiv:2002.06563 , .Rathai, K. M. M., Alamir, M., Sename, O., & Tang, R. (2018). A parameterized NMPC scheme for embedded control ofsemi-active suspension system.
IFAC-PapersOnLine , , 301–306.Rathai, K. M. M., Sename, O., & Alamir, M. (2019). GPU-based parameterized nmpc scheme for control of half car vehiclewith semi-active suspension system. IEEE Control Systems Letters , , 631–636.Rocha Filho, T. M., dos Santos, F. S. G., Gomes, V. B., Rocha, T. A., Croda, J. H., Ramalho, W. M., & Araujo, W. N. (2020).Expected impact of COVID-19 outbreak in a major metropolitan area in Brazil. medRxiv , .Rodriguez-Morales, A. J., Gallego, V., Escalera-Antezana, J. P., Mendez, C. A., Zambrano, L. I., Franco-Paredes, C., Su´arez,J. A., Rodriguez-Enciso, H. D., Balbin-Ramon, G. J., Savio-Larriera, E. et al. (2020). COVID-19 in latin america: theimplications of the first confirmed case in Brazil. Travel medicine and infectious disease , .Sarkar, K., Khajanchi, S., & Nieto, J. J. (2020). Modeling and forecasting the covid-19 pandemic in india.
Chaos, Solitons &Fractals , , 1–16. doi: https://doi.org/10.1016/j.chaos.2020.110049 .Silva, R. R., Velasco, W. D., da Silva Marques, W., & Tibirica, C. A. G. (2020). A bayesian analysis of the total number ofcases of the COVID-19 when only a few data is available. a case study in the state of Goias, Brazil. medRxiv , .Sun, P., Lu, X., Xu, C., Sun, W., & Pan, B. (2020). Understanding of COVID-19 based on current evidence. Journal of medicalvirology , , 548–551.Sun, T., & Wang, Y. (2020). Modeling covid-19 epidemic in heilongjiang province, china. Chaos, Solitons & Fractals , ,1–5. doi: https://doi.org/10.1016/j.chaos.2020.109949 .The Lancet (2020). Covid-19 in Brazil: so what?. The Lancet , , 1461. doi: https://doi.org/10.1016/S0140-6736(20)31095-3 .Werneck, G. L., & Carvalho, M. S. (2020). The COVID-19 pandemic in Brazil: chronicle of a health crisis foretold.Zacchi, L. L., & Morato, M. M. (2020). The COVID-19 pandemic in Brazil: An urge for coordinated public health policies.URL: https://hal.archives-ouvertes.fr/hal-02881690 preprint.preprint.