key: cord-0663834-jnpysl0s
authors: Sauerteig, Philipp; Esterhuizen, Willem; Wilson, Mitsuru; Ritschel, Tobias K. S.; Worthmann, Karl; Streif, Stefan
title: Model Predictive Control Tailored to Epidemic Models
date: 2021-11-12
journal: nan
DOI: nan
sha: 8cd15101289403cf6152dfc6f6b720e70f7340bc
doc_id: 663834
cord_uid: jnpysl0s

We propose a model predictive control (MPC) approach for minimising the social distancing and quarantine measures during a pandemic while maintaining a hard infection cap. To this end, we study the admissible and the maximal robust positively invariant set (MRPI) of the standard SEIR compartmental model with control inputs. Exploiting the fact that in the MRPI all restrictions can be lifted without violating the infection cap, we choose a suitable subset of the MRPI to define terminal constraints in our MPC routine and show that the number of infected people decays exponentially within this set. Furthermore, under mild assumptions we prove existence of a uniform bound on the time required to reach this terminal region (without violating the infection cap) starting in the admissible set. The findings are substantiated based on a numerical case study.

As the ongoing COVID-19 pandemic demonstrates, it is important to understand how infectious diseases spread and how countermeasures may affect this spread. Evaluating the impact of countermeasures based on mathematical models is an active field of research. In particular, there is a large body of work applying optimal control to epidemiology. Early work on this topic includes [19] , where vaccination policies are derived for an SIR (susceptible-infectiveremoved) model via dynamic programming; and [28] which solves optimal control problems for an SIS (susceptible-infective-susceptible) model. Many papers apply optimal control to SIR models, [16, 12, 20, 1] . See also the survey [30] .

Model predictive control (MPC) has also been applied to epidemiology, such as in the papers [29, 21, 31] which consider stochastic networked epidemic models. Recently, many papers have appeared on the modelling and control of COVID-19 via MPC, see for example [5, 22] . The papers [2, 26, 14, 24] consider MPC of COVID-19 models where the control is limited to non-pharmaceutical interventions, whereas the papers [13, 25] also consider vaccination strategies.

Set-based methods have also been applied to epidemiology with the specific goal of control design that not only eliminates a disease, but also respects hard infection caps. In [23] the authors describe the admissible set (also known as the viability kernel in viability theory [3] ) of a two-dimensional model of a vector-borne disease, such as dengue fever or malaria. The paper [8] extended this study, showing that the malaria model's admissible set as well as the maximal robust positively invariant set (MRPI) may be described via the theory of barriers. In [7, 9] it is shown that parts of the boundaries of these sets are made up of special integral curves of the system that satisfy a minimum-/maximum-like principle. The recent paper [10] uses the theory of barriers to describe both sets for the well-known SIR and SEIR epidemic models with and without plant-model mismatch.

The papers [23, 8, 10] argue that to maintain an infection cap, the control should be chosen based on the location of the state with respect to these sets. If it is located in the MRPI, intervention measures may be relaxed and economic damage may be minimized. However, if the state is located in the admissible set then intervention measures need to be carefully considered to avoid a breach in the infection cap. Some other research concerned with applying set-based methods to epidemiology include work on the set of sustainable thresholds, [4] , and the paper [27] , which describes the viability kernel of an SIR model through the solution of an associated Hamilton-Jacobi equation.

In this paper, we combine the results from the paper [10] with model predictive control (MPC) [15] . We consider the SEIR model and impose a hard cap on the proportion of infective individuals as well as constraints on the contact rate and removal rate. We first show that it is possible to eliminate the disease from any initial state located in the admissible set while always satisfying the infection cap, and that this is always the case if the initial value is contained in the system's MRPI. Based on our findings, we then construct terminal conditions ensuring that the terminal costs used in our MPC implementation are uniformly bounded and prove initial feasibility under mild assumptions.

In Section 2, we introduce the constrained epidemic model and present important facts regarding the system's admissible set and the MRPI. Section 3 presents the MPC implementation that takes advantage of these sets. In this section we also present the paper's main result (Proposition 1) ensuring initial feasibility and eradication of the disease if the prediction horizon is sufficiently long. A numerical case study illustrates our findings in Section 4.

For an interval I ⊆ R and a set Y ⊆ R m we denote the set of all locally integrable functions from I to Y by

In this section, we introduce the SEIR model and the system's admissible set as well as the MRPI.

We consider the SEIR model [18] S(t) = −β(t)S(t)I(t), S(0) = S 0 (1a)

where S(t), E(t), I(t), and R(t) describe the fractions of people who are either susceptible, exposed, infectious, or removed at time t ≥ 0 and (S 0 , E 0 , I 0 , R 0 ) ∈ [0, 1] 4 with S 0 + E 0 + I 0 + R 0 = 1 denotes the initial value. In this context, compartment R accounts for both recovery and death due to a fatality caused by the disease. The parameter η −1 > 0 denotes the average incubation time in days. Furthermore, in the standard formulation of the SEIR model, β ≡ β nom > 0 is the rate of infectious contacts and γ ≡ γ nom > 0 is the removal rate. However, we take countermeasures into account by allowing β and γ to be time-varying control variables. In particular, a value β(t) < β nom reflects a reduction of the rate of infectious contacts, e.g. via contact restrictions or hygiene measures, while γ(t) > γ nom can be interpreted as quarantining, and, thus, removing infectious people. These considerations motivate our control constraints

with β min > 0 and γ max < ∞. Moreover, we aim to maintain a hard infection cap by satisfying the state constraint

for some I max ∈ (0, 1). Note that R does not affect the other compartments. Thus, we may ignore it and analyse the three-dimensional system (1a)-(1c) under input and state constraints (2) and (3). To ease our notation, we use x(t) = (x 1 (t), x 2 (t), x 3 (t)) ⊤ := (S(t), E(t), I(t)) ⊤ ∈ R 3 to denote the state and u(t) := (β(t), γ(t)) ⊤ ∈ R 2 to denote the control input at time t ≥ 0 and writeẋ(t) = f (x(t), u(t)) short for (1a)-(1c). Furthermore, we denote the set of feasible control values by U := [β min , β nom ] × [γ nom , γ max ] ⊂ R 2 and the set of feasible controls as U = L 1 loc ([0, ∞), U ). Thus, for any initial value x 0 and control input u ∈ U there exits a unique solution of the initial value problem (1). We then write x(·; x 0 , u) to denote the trajectory with respect to the three compartments S, E, and I and highlight the dependence on x 0 and u. Based on g : R 3 → R, x → x 3 − I max , the set of feasible states is given by G := x ∈ R 3 g(x) ≤ 0 and its interior by G − . Furthermore, the state trajectory satisfies

In other words, the system is confined to the positively invariant set

Given an initial state, x 0 ∈ G Π := G ∩ Π, we aim to determine a control u = u(x 0 ) that eliminates the disease, i.e., lim t→∞ (E(t) + I(t)) = 0, while maintaining the hard infection cap, i.e., I(t) ≤ I max for all t ∈ [0, ∞). To this end, we make use of the admissible set and the MRPI. Remark 1. For γ > γ nom , i.e., if quarantine measures are active, a certain (additional) proportion of the infectious people are moved to compartment R, see (1c). Consequently, constraint (3) does not reflect a bound on the actual number of infectious people, but rather on those, who are able to infect others. However, it suffices to illustrate our approach. For a more elaborate way to enforce an infection cap, we refer to [14, 13] .

The admissible set contains all states, for which there exists an input such that the state constraints are satisfied for all future time.

Definition 1. The admissible set for the system (1)-(3) is given by

Moreover, we consider the set of states for which any control ensures feasibility.

In [10] , both sets have been characterized using the theory of barriers [7, 9] . In particular, the sets for system (1) with constraints (2)-(3) are compact and never empty with M ⊆ A.

In order to eliminate the disease, we need to drive the state to the set of equilibria given by

Clearly, E ⊂ M. The next lemma summarizes some facts about the SEIR model in terms of A and M that will be useful in the sequel. Lemma 1. The following assertions hold true.

1. For every x 0 ∈ Π and any u ∈ U the epidemic will die out asymptotically, i.e., the compartments satisfy lim t→∞ E(t) = lim t→∞ I(t) = 0 as well as lim t→∞ (S(t)+R(t)) = 1.

For every x 0 ∈ A there exists a u ∈ U such that lim t→∞ E(t) = lim t→∞ I(t) = 0 and I(t) ≤ I max for all t ≥ 0.

Proof. 1) The idea is based on the proof of Theorem 5.1 in [17] , where the assertion has been shown for the SIR model. Let x 0 ∈ Π and u ∈ U be arbitrary. Since S and R are monotonic and bounded, the limits S ∞ := lim t→∞ S(t) and R ∞ := lim t→∞ R(t) exist with lim t→∞Ṡ (t) = lim t→∞Ṙ (t) = 0. Hence,

and, therefore,

Also, noting that A ⊆ Π, the assertion follows immediately from 1).

3) Since x 0 ∈ M, x(t; x 0 ; u) ∈ G Π for all t ≥ 0 for any u ∈ U. Also, noting that M ⊆ Π, the assertion follows immediately from 1).

In Section 3, the set

will be used to define terminal constraints within MPC. We state the following lemma to motivate our choice.

The set X f is contained in M. Moreover, the compartments E and I decay exponentially for all x 0 ∈ X f and all u ∈ U.

Proof. Clearly, any point x 0 with S 0 ∈ [0, 1] and E 0 = I 0 = 0 is in M. Let x 0 ∈ X f with E 0 + I 0 > 0 be given and consider u = u nom . First, note that if I 0 = I max , theṅ

Due to continuity, there exists some (sufficiently small) time δ > 0 such that I(δ) > 0 and x(δ) ∈ X f . Since S decreases monotonically, the compartments E and I are bounded from above by the solution of the initial value probleṁ

In particular,Ė

holds. The solution of (4) with u = u nom is given by

The eigenvalues of A are

Hence, S(δ) < S 0 ≤ γ/β yields exponential stability of A. Consequently, e and i and, thus, E and I decay exponentially.

As in the SIR model [17] , S(t) = γ(t)/β(t) is a threshold in the following sense

i.e., as long as S is sufficiently big, the total number of infected (exposed or infectious) people increases monotonically and decays otherwise. In the following, we use the notation S = γ nom /β nom to denote the threshold for herd immunity, which is also the inverse of the basic reproduction number [18] .

We propose an MPC scheme for determining a feedback control that aims to minimise the required contact restrictions and quarantine measures while maintaining a hard infection cap, i.e., given an initial value x 0 ∈ G Π we are interested in solving the optimal control problem

For an introduction to MPC we refer to [15] .

We present a solution to the stated problem using continuous-time MPC. Thus, we consider the following finite-horizon optimal control problem,

with terminal cost

Here, E nom ( · ;x) and I nom ( · ;x) denote the exposed and infective compartments of the state trajectory obtained with the nominal input u nom ≡ (β nom , γ nom ) ⊤ starting atx. Given an initial state x 0 ∈ G Π and a horizon length T ∈ R >0 ∪ {∞}, the set U X f T (x 0 ) denotes all control functions u ∈ U for which x(t; x 0 , u) ∈ G Π for all t ∈ [0, T ], and x(T ; x 0 , u) ∈ X f . The optimal value function is defined as V T (x 0 ) :

If U X f T (x 0 ) = ∅, i.e., there does not exist a solution to OCP T , we take the convention that V T (x 0 ) = ∞. The continuous-time MPC algorithm is outlined in Algorithm 1.

We assume that if U X f T (x 0 ) = ∅ in step 1), then a minimizer is an element of U X f T (x 0 ), see [15, p. 56 ] for a discussion on this issue. Thus, given an initial state x 0 ∈ A, OCP T is solved with the finite horizon T = N δ, the first portion of the minimizer is applied to the system, the state is measured, and the process is iterated indefinitely. The algorithm implicitly produces the time-varying sampled-data MPC feedback, µ T,δ : [0, δ) × G Π → U , which results in the MPC closed-loop solution, denoted x µ T ,δ (·; x 0 ).

3. Set x 0 ← x(δ; x 0 , u ⋆ ) and go to step 1). J T (x 0 , u), we denotex := x(δ; x 0 , u ⋆ ) and definê u : [0, T ) → U viaû(t) = u ⋆ (δ + t) for t < T − δ and u(t) = u nom otherwise. Then, g(x(t;x,û)) ≤ 0 for all t ≥ 0 and x(T − δ;x,û) ∈ X f , i.e.,û ∈ U X f T (x). Thus, initial feasibility of OCP T yields recursive feasibility.

The following result is essential to prove that the MPC implementation works. Lemma 3. Without countermeasures the infinite horizon cost functional is uniformly bounded on the MRPI, i.e.,

Proof. Let x 0 ∈ M be given. Then, any control u ∈ U yields x(t; x 0 , u) ∈ G Π for all t ≥ 0. Moreover, I ∞ := lim t→∞ I(t) = lim t→∞ E(t) = 0 and S ∞ := lim t→∞ S(t) and R ∞ := lim t→∞ R(t) exist. Hence,

and, therefore,

Analogously, one can show that

In conclusion, we get

Next we show initial (and, thus, recursive) feasibility if the following mild assumption holds. Assumption 1. If there are infected people, i.e., E 0 + I 0 > 0, then we assume the initial values to be bounded from below, i.e., there exists some ε 0 > 0 such that

We further assume there exists some K ∈ R >0 such that

i.e., the fraction of infected people who are infectious is bounded from below. In order to facilitate the proof of Lemma 4, we further assume that

Note that this can always be achieved by allowing stricter social distancing or quarantine measures.

Our numerical simulations in Section 4 show that Assumption (6) is essential to reach the terminal set in finite time, see Figure 3 . However, from a practical point of view, we are only interested in a situation where there already are infected people.

Based on Assumption 1, we are able to show that the terminal set is reached in finite time.

Lemma 4. Let assumption (7) hold. Then, the time required to reach the terminal set is uniformly bounded on the subset A ′ := { x 0 ∈ A | (6) holds }, i.e.

Proof. For x 0 ∈ M, the assertion follows from Lemma 3. Consider an arbitrary x 0 ∈ A ′ \ M.

Then, x 0 / ∈ X f according to Lemma 2. We construct a controlũ ∈ U such that x(t; x 0 ,ũ) ∈ G Π for all t ≥ 0 and x(T ; x 0 ,ũ) ∈ X f for some finite T > 0. To this end, we ensure a decay of S as follows.

Note that since x 0 ∈ A ′ \ M and lim t→∞ I(t) = 0 there exists some control u 1 ∈ U satisfying x(t; x 0 , u 1 ) ∈ G Π for all t ≥ 0 and some time t 1 ≥ 0 such that I(t 1 ; x 0 , u 1 ) = 1/2 · I max . Moreover, based on assumption (8) we may assume that γ nom /β nom < 1 < γ max /β min . since otherwise, the pandemic would simply abate without control, see Remark 2. Therefore, there exists some time t 3 ∈ (t 1 , ∞] and a control u 2 = (β 2 , γ 2 ) ⊤ ∈ U such that

and, hence, the number of infected E+I is constant on the interval [t 1 , t 3 ]. Since lim t→∞ (E(t)+ I(t)) = 0, this implies t 3 < ∞.

Next, we show that there exists some t 2 ∈ [t 1 , t 3 ] and a control u 3 = (β 3 , γ 3 ) ⊤ ∈ U such that E + I is constant on t ∈ [t 2 , t 3 ] and, in addition,

Suppose ηE(t)/I(t) < γ nom for all t > t 1 . Then,Ė(t) > 0 for all t > t 1 , in contradiction to lim t→∞ E(t) = 0. On the other hand, if ηE(t)/I(t) > γ max for all t > t 1 , thenİ(t) > 0 for all t > t 1 , in contradiction to lim t→∞ I(t) = 0. Therefore, t 2 and u 3 as above exist with γ 3 (t) = ηE(t)/I(t) and β 3 (t) = ηE(t 2 )/(S(t)I(t 2 )) for all t ∈ [t 2 , t 3 ]. (The feasibility of β 3 can be argued analogously to the one of γ 3 .) Thus, we have

for all t ∈ [t 2 , t 3 ]. With the maximal t 3 we arrive at S(t 3 ) =S and, thus,

Due to assumptions (6) and (7), the value

for all t ≥ 0 with S(t) > γ nom /β nom .

As argued in the proof of Lemma 4, assumption (7) is not too restrictive. Whenever I(t) is small while E(t) is big at the same time, thenİ(t) > 0, meaning that I is going to increase. Thus, after some time into the pandemic, there will always be a certain fraction of infected people who are infectious.

Our main result states that the MPC feedback generated by Algorithm 1 approximates the solution of (5). (7) hold. There exists a finite prediction horizon T > 0 such that OCP T (5) is initially and, thus, recursively feasible for every x 0 ∈ A ′ with lim t→∞ x µT ,δ (t, x 0 ) ∈ E, and g(x µT ,δ (t, x 0 )) ≤ 0 for all t ≥ 0 under the MPC feedback µ T,δ produced by Algorithm 1.

Proof. From the proof of Lemma 4, the target set X f is reachable from any x 0 ∈ A ′ in finite time. Thus, there exists a finite T > 0 such that for every x 0 ∈ A ′ the problem OCP T is initially feasible. Hence, OCP T is recursively feasible, implying that the infection cap is respected for all t ≥ 0. After the state reaches X f ⊂ M, it approaches E asymptotically according to Lemma 1. 

We now illustrate the paper's theoretical results with simulations in Matlab, using the script provided by [15] to implement the MPC controller. We consider the constrained SEIR model (1) under constraints (2) and (3). We take β nom = 0.44, γ nom = 1/6.5, η = 1/4.6 and I max = 0.05 [13] . Furthermore, we take γ max = 0.5 and β min = 0.22. The point clouds of the boundaries of M and A, indicated in Figure 1 , are obtained as detailed in the paper [10] , that is, via backward numerical integration of the system with a special "extremal input" from particular initial points on the state constraint I = I max . Also shown is a typical run of the MPC closed loop as produced by Algorithm 1 with δ = 1. The prediction horizon is set to T = 25 days and the initial value x 0 = (0.50, 0.18, 0.01) ⊤ is close to the boundary of A. The impact of imperfect predictions and state estimation is left for future research.

In Figure 2 the evolution of I subject to the MPC feedback µ T is depicted. The trajectory follows the construction in the proof of Lemma 4. First, I increases until it reaches some controlled equilibriumİ(t 2 ) = 0 until herd immunity is achieved with S(t 3 ) =S. Then, it decays asymptotically towards zero. Note that the optimal equilibrium in this context is the infection cap, i.e., I ≡ I max since it maximises the descent of S. Furthermore, we weight the cost terms in (5d) equally since tuning weights is out of scope of this paper. As a result, the cost term with respect to γ T is much bigger than the one associated with β T . Keep in mind that in our model, only infectious people are put into quarantine and we do not consider re-infections, i.e., they are completely removed from the system dynamics. Figure 3 motivates assumption (6) . We observe that starting close to S 0 =S and decreasing E 0 and I 0 , the time for a pandemic to break out explodes. Moreover, the overall cost J ∞ (x 0 , u nom ) seems to be independent of the choice of E 0 and I 0 . One reason for that is the fact that E +I increases as long as S(t) >S. Hence, for the pandemic to die out, S 0 −S many people have to get infected first. Consequently, we do not have cost controllability on A. 1 Figure 3 : Impact of initial value for small x 0 = (S + 0.01, ε 0 , ε 0 ) ⊤ without control, i.e., u ≡ u nom and with machine precision ε = 2.2204 · 10 −16 .

In this paper, we considered the SEIR compartmental model with control inputs representing social distancing and quarantine measures. Based on a hard infection cap, we determined a subset in the state space, where the pandemic is contained without enforcing countermeasures. We used this set to define terminal constraints for our MPC algorithm and showed initial and, thus, recursive feasibility under mild assumptions. Our numerical simulations show that the approach is suitable to maintain the infection cap while keeping the total amount of time short where countermeasures have to be enforced. We found that in order for the pandemic to abate, sufficiently many people have to be infected first, i.e., herd immunity needs to be established. Our approach exploits that S is strictly decreasing. Taking re-infectious or births into account requires further investigation.

Optimal control and cost-effective analysis of malaria/visceral leishmaniasis co-infection

A model for COVID-19 with isolation, quarantine and testing as control measures

Viability theory: new directions

Sustainable thresholds for cooperative epidemiological models

Optimal control techniques based on infection age for the study of the COVID-19 epidemic

Model predictive control, cost controllability, and homogeneity

On barriers in state and input constrained nonlinear systems

Maintaining hard infection caps in epidemics via the theory of barriers

On maximal robust positively invariant sets in constrained nonlinear systems

Epidemic management with admissible and robust invariant sets

Recursive feasibility of continuoustime model predictive control without stabilising constraints

A control theory approach to optimal pandemic mitigation

How to Coordinate Vaccination and Social Distancing to Mitigate SARS-CoV-2 Outbreaks

How much testing and social distancing is required to control COVID-19? Some insight based on an age-differentiated compartmental model

Nonlinear model predictive control: Theory and algorithms

Optimal control of epidemics with limited resources

Three basic epidemiological models

The mathematics of infectious diseases

Optimal vaccination schedules in a deterministic epidemic model

Optimal control of the chemotherapy of HIV

Dynamic resource allocation to control epidemic outbreaks a model predictive control approach

Robust and optimal predictive control of the COVID-19 outbreak

Viable control of an epidemiological model

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

A model predictive control approach to optimally devise a two-dose vaccination rollout: A case study on COVID-19 in Italy

Optimal control of the COVID-19 pandemic with nonpharmaceutical interventions

A model for a vector-borne disease with control based on mosquito repellents: A viability analysis

Quantitative guidelines for communicable disease control programs. Biometrics

Dynamic control of modern, networkbased epidemic models

Optimal control in epidemiology

Robust economic model predictive control of continuous-time epidemic processes