Epidemic control via stochastic optimal control

Epidemic control via stochastic optimal control

Introduction

The classic SIR model, originally proposed in KM27, is the standard tool of mathematical epidemiology BCF19 for quantitative analysis of the spread of an epidemic. It describes the state of the affected population in terms of three state variables, traditionally denoted by $S,I$, and $R$:

The model dynamics is given as the following dynamical system:

where the constant parameters $\beta>0$ and $\gamma>0$ are called the infection rate and recovery rate, respectively. The evolution above is subject to the initial conditions:

Notice that this dynamics obeys the conservation law

consistent with the assumption that the variables $S,I$, and $R$ represent population fractions. This means that the variable $R$ is, in a way, redundant, as its current value does not affect the dynamics of $S$ and $I$, and it can be computed in a straightforward manner from (ref). It is thus natural to consider the two dimensional system defined by $S,I$ only.

Eventually every epidemic comes to a natural halt\footnote{In the SIR model (ref) this is reflected by the fact that the solution $(S(t),I(t))=(0,0)$ is globally asymptotically stable.}, but its impact on the population may be very serious. The overall objective of epidemic control is to slow down the spread of infection in a way that it does not overwhelm the healthcare system and it allows the economy to function. All of this should be done within the limits of available resources.

In this note we study the problem of optimal control of an epidemic modeled by means of a stochastic extension of the SIR model (see Section (ref) for definition). We assume that the controlling agent (“government”) has the ability to impact the spread of the disease through one of the policies (or the combination of both):

The controlled dynamics of the SIR model reads AM79, HD11:

For efficiency, we will be using the notation $X_1=S$ and $X_2=I$ throughout the remainder of the paper.

Mathematical models are only as good as (i) their analytic specifications, and (ii) the data that fuel them. During initial phases of an epidemic the available data tend to be of limited usefulness: because of the lack of reliable large scale testing, it is not really known what fractions of the population fall into the different compartments $S$, $I$, and $R$. This may lead to a panic reaction of the population and a chaotic and economically devastating public health response to the epidemic. In the absence of an effective vaccine (which would allow to immunize a portion of the population) the optimal policy is to isolate at least a significant fraction of the infected individuals so that the basic reproduction ratio $R_0$ can be brought significantly below one. Since we lack the knowledge who is infected and who is not, the public health response is to try to isolate everyone, whether susceptible, infected or immune.

These circumstances impose a serious limitation on practical applicability of the approach to optimal epidemic control discussed in this paper, as well as other quantitative approaches. Unless the inputs to the model ($\beta$, $\gamma$, and the current snapshots of $S$ and $I$) can reliably be estimated from the available data, the model's output is unreliable\footnote{“This is indeed a mystery,” Watson remarked. “What do you imagine that it means?”\\ “I have no data yet,” Holmes replied. “It is a capital mistake to theorise before one has data. Insensibly one begins to twist facts to suit theories, instead of theories to suit facts.” ACD}.

The paper is organized as follows. In Section (ref) we review the formulation of the continuous time stochastic SIR model. The optimal control problem for the stochastic SIR model is formulated in Section (ref). The optimal control problem is recast as the stochastic minimum principle problem and formulated in terms of a system of forward backward stochastic differential equations (FBSDE). We present an algorithm for solving this system in Section (ref). This section presents also the results of a number of numerical experiments involving the three mitigation policies within various cost function regimes.

Acknowledgement. I would like to thank Nicholas Lesniewski for numerous discussions.

Stochastic SIR model

We consider a continuous time stochastic extension of the deterministic SIR model, see BBSG17 and references therein. Let $W_t$ denote the standard Brownian motion, and let $\dot W_t$ denote the white noise process. We assume that the infection rate $\beta$ is subject to random shocks, rather than being constant, namely

while the recovery rate $\gamma$ remains constant. Here $\sigma>0$ is a constant volatility parameter. This leads to the following system of stochastic differential equations (SDE), driven by $W_t$

with the initial conditions

The third component of the process, $X_3=R$, follows the dynamics

which implies that the conservation law

continues to hold in the stochastic model.

Notice that, under the stochastic SIR model (ref), an epidemic eventually comes to a natural end. More precisely, the solution $(X_{1,t},X_{2,t})=(0,0)$ to (ref) is stable in probability K12. In order to see it, we set

for $0\leq X_i\leq 1$, $i=1,2$, and fixed $0<\rho<\beta/\gamma$. Then $V_\rho(0,0)=0$, and $V_\rho(X_1,X_2)>0$ in a neighborhood of $(0.0)$. Furthermore, denoting by $\mathcal{L}$ the generator of the process(ref), we verify that

In other words, $V_\rho(X_1,X_2)$ is a Lyapunov function for (ref) and our claim follows from Theorem 5.3 in K12.

The model (ref) is a one factor model, driven by a single source of randomness. There is a natural two-factor version of the stochastic SIR model LL20, in which also the recovery rate $\gamma$ is allowed to be subject to white noise shocks. For simplicity, our analysis will focus on the one-factor model (ref).

HJB equation

We frame the problem of epidemic control as a stochastic control problem B57. We denote by $u=(u_1,u_2)\equiv(v,i)$ the vaccination and isolation controls, and we denote the controlled process by $X^u_t$. Generalizing the deterministic specification (ref) to the stochastic case, we assume that the dynamics of $X^u_t$ is given by:

subject to initial conditions (ref).

Two special cases of the controlled process are of interest. If a vaccine against the disease is unavailable, we set $u_1=0$ in the equation above, which yields the following controlled process:

We will refer to this policy as an isolation policy.

Similarly, we can consider a vaccination policy, for which $u_2=0$. In this case the controlled dynamics reads

We assume a finite time horizon $T<\infty$. The controlling agent's objective is to minimize a running cost function $c(X_t,u_t)$ and the terminal value function $G(X_T)$. In other words, we are seeking a policy $u^\ast$ such that

We consider the following cost function:

where

for $i=1,2$. In other words, the running cost of vaccination is assumed to be proportional to the number of susceptible individuals, while the cost of isolation is assumed to be proportional to the number of infected individuals. The coefficients $L_i>0, M_i, N_i$ are determined by the overall cost of following the mitigation policy. In particular, they should be selected so that the running cost functions are strictly positive.

As the terminal value function we take the transmission rate of the infection HD11, namely

We now invoke stochastic dynamic programming, see eg. FS92, P09. The key element of this approach is the value function $J(t,x,y)$. It is determined by two requirements:

Using Ito's lemma, we verify that these two conditions lead to the following nonlinear partial differential equation for the value function, namely the stochastic Hamilton-Jacobi-Bellman equation:

subject to the terminal condition

As the first step towards solving the HJB equation (ref), we let $u=u^\ast$ denote the minimizer of the expression inside the curly parentheses in (ref). In other words, $u^\ast$ satisfies

which leads to the following first order condition on $u^\ast$:

Substituting $u^*$ back to the HJB equation yields the equation

We do not believe that the solution to this equation can be explicitly represented in terms of standard functions. This is a three dimensional partial differential equation, and solving it numerically may pose challenges. Rather than following this path, we shall invoke the stochastic minimum principle and reformulate the problem as a system of FBSDEs. Among the advantages of this approach is that it might be amenable to a deep learning approach via the method advocated in HJE18.

Stochastic minimum principle

The stochastic minimum principle, see P09 and references therein, offers an alternative approach to stochastic optimal control. It is a stochastic version of Pontryagin's minimum principle introduced in the context of deterministic optimal control PBGM92. It also offers an effective numerical methodology to solving the HJB equation (ref) via Monte Carlo simulations. In this approach, the key object is the Hamiltonian function ${\mathcal{H}}={\mathcal{H}}(x,u,y,z)$ of four arguments and a system of stochastic differential equations, both forward and backward, determined by ${\mathcal{H}}$.

Specifically, for the case of the controlled stochastic SIR model (ref), we have $x,u,y,z\in\mathbb{R}^2$, and the Hamiltonian function reads:

We consider the following system of stochastic Hamilton's equations:

and

where

Equation (ref) is merely an alternative way of writing the underlying controlled diffusion process (ref), while equation (ref) reflects the dynamics of the control variables given the running cost function. Note that while the first of the equations (ref) is a standard (forward) SDE, the second one is a backward stochastic differential equation (BSDE), see e.g. EPQ97.

Explicitly, these four equations can be stated as

subject to the boundary conditions in (ref).

Let now choose $u^\ast$ to be the policy such that $X_t^\ast\equiv X_t^{u^\ast}$ is a solution to (ref) and

Then $u^\ast$ is an optimal control, i.e. it satisfies the optimality condition (ref).

In fact, there is a direct link between the Hamilton function ${\mathcal{H}}$ and the value function $J$ (also known in classical dynamics as the action function). Namely, we set P09, YZ99:

If $u^\ast$ is an optimal control and $X^\ast$ denotes the corresponding optimal diffusion process, then the pair

is the solution to the BSDE (ref).

Going back to the main line of reasoning, we find that $u^\ast$ has to satisfy

From equations (ref) we see that, up to a simple linear transformation, $Y_t$ is essentially the optimal policy. Furthermore, the process $Z_t$ can be thought of as the sensitivity of the optimal policy to the underlying process $X_t$ (multiplied by the instantaneous volatility of that process).

Substituting this expression into (ref), we find that the optimal process has to follow the dynamics:

subject to the boundary conditions in (ref) and (ref). In particular, the isolation only and vaccination only policies are given by the following systems of FBSDEs:

and

respectively.

Even though derived in the context of a meaningful underlying dynamics, there is no a priori reason why these equations should have solutions. The drivers of the backward equations above contain terms quadratic in $Y$, and so the standard existence theorems EPQ97 do not apply. We proceed in the following assuming that these systems do, in fact, have solutions.

Numerical solution

In this section we discuss a numerical algorithm for solving the system (ref). Applying this methodology in a number of numerical experiments, we present compelling evidence that solutions to (ref) exist and are meaningful over a wide range of model parameters.

Simulating the optimal FBSDE

We start by describing a method for discretization of the basic FBSDE system (ref). We notice first that the two dimensional system defined by the state variables $S,I$ is in fact Hamiltonian LR04. Namely, we define new (canonical) variables

and set\footnote{Note that the Lyapunov function $V$ of Section (ref) is the negative of a regularized version of $H$.}

It is easy to verify that

under the dynamics (ref). The system (ref) can be explicitly written in the Hamilton form

Notice also that the Hamiltonian function is separable, i.e. it is of the form $H(p,q)=T(p)+V(q)$\footnote{There is a practical dimension to these theoretical observations The fact that the SIR model is a Hamiltonian allows us to solve it numerically by means of the extremely efficient St\"ormer-Verlet method HLW03, LR04. We will not make use of this method, as it does not seem to directly extend to the stochastic case}.

A convenient discretization to (ref) can be formulated in terms of the canonical variables $(p_t,q_t)=(-\log X_{2,t},\,-\log X_{1,t})$ as follows LL20. Choose the number of steps $\frak{n}$, define the time step $\delta=T/\frak{n}$, set $t_n=n\delta$ for $n=0,1,\ldots,\frak{n}$, and use simplified $(p_n, q_n)\equiv(p_{t_n},q_{t_n})$. This yields the following Euler scheme:

for $n=0,1,\ldots,\frak{n}-1$. Here, the Brownian motion increments $\Delta W=W_{(n+1)\delta}-W_{n\delta}$ are independent variates drawn from the normal distribution $N(0,\delta)$. At each iteration step also have to to floor $q_{n+1}$ and $p_{n+1}$ at $0$, $q_{n+1}=\max(q_{n+1},0)$, $p_{n+1}=\max(p_{n+1},0)$, so that $e^{-q_n}, e^{-p_n}\leq 1$.

Approximating the backward equations of the system (ref) leads to the following backward Euler scheme:

where

denote the generators of the two BSDEs. Starting with the terminal condition

we will move backward in time with computing $Y_n$ and $Z_n$, for $n=\frak{n}-1,\ldots, 0$.

Notice that two difficulties arise while doing so: (i) the $Y_n$'s in (ref) are not necessarily adapted, and (ii) they depend on $Z_n$. These two problems can be solved by taking conditional expectations $\mathsf{E}_n(\,\cdot\,)=\mathsf{E}(\,\cdot\,|\,X_0,X_1,\ldots,X_n)$ on time $t_n$. This leads to the following condition:

This is an implicit scheme, which may slow down the computations. However, we can easily replace it with an explicit scheme of the same order:

In order to determine $Z_{i,n}$, $i=1,2$, we multiply (ref) by an increment $\Delta W_n$ and take conditional expectations. This yields

and hence we obtain the following expression for $Z_{i,n}$:

In summary, we are led to the following discrete time scheme for solving (ref):

for $n=\frak{n}-1,\ldots, 0$. Note that simulating this system requires numerical estimation of the conditional expected values $\mathsf{E}_n(\,\cdot \,)$ in the formulas above. We discuss this issue in the following section.

Computing the conditional expectations

A practical and powerful method of computing the conditional expected values in (ref) is the Longstaff-Schwartz regression method originally developed for modeling American options LS01. We use a variant of this method that involves the Hermite polynomials, and that was used for a similar purpose in LR16. This choice is natural, as conditional expectations of Hermite polynomials of a Gaussian random variable lead to simple closed form expressions.

Let $\mathrm{He}_k(x)$, $k=0,1,\ldots$, denote the $k$-th normalized Hermite polynomial corresponding to the standard Gaussian measure $d\mu(x)=(2\pi)^{-1/2}\,e^{-x^2/2} dx$. These functions form an orthonormal basis for the Hilbert space $L^2 \big(\mathbb{R}, d\mu\big)$.

The key property of $\mathrm{He}_k(x)$ is the following addition formula for $\chi \in [0,1]$ and $w,x \in \mathbb{R}$:

Consequently, integrating over $x$ with respect to to the measure $\mu$ yields the following conditioning rule:

Here, $w, x$ are independent standard normal random variables. We shall use this relation in order to estimate the conditional expected values in (ref).

We set $W_{t_i}=\sqrt{t_i}\,w_i$, for $i=1,\ldots, m$, where $w_i$ is an $n$-dimensional standard normal random variable. We notice that

where $\chi_i=t_i/t_{i+1}$, and where $X_i$ is standard normal and independent of $w_i$. In the following, we shall use this decomposition in conjunction with (ref).

Now, we assume the following linear architecture:

where $K$ is the cutoff value of the order of the Hermite polynomials. This is simply a truncated expansion of the random variable $Y_{i+1}$ in terms of the orthonormal basis $\mathrm{He}_k(w_{i+1})$. The values of the Fourier coefficients are estimated by means of ordinary least square regression. Then, as a consequence of the conditioning rule (ref),

In other words, conditioning $Y_{i+1}$ on $w_{i}$ is equivalent to multiplying its Fourier coefficients $g_k$ by the factor $\chi_i^{k/2}$. This allows us to calculate the first term on the right hand side of the third equation in (ref). In order to calculate the second term, we substitute the explicit formula for $Z_n$ into the generators (ref) and repeat the calculations in (ref) and (ref) with $Y_{n+1}$ replaced by $Y_{n+1} + f(X_n,Y_n,Z_n)\delta$.

Numerical experiments

For numerical experiments within the framework developed above, we assume a time horizon of 1 year (365 days), and choose the SIR model parameters as follows:

These parameters are purely hypothetical, and they do not arise from the calibration to any actually reported data. The corresponding value of the basic reproduction ratio $R_0=S_0\beta/\gamma$ is 3.30 and it indicates a highly infectious disease such as COVID-19. We choose the diffusion parameter

A typical scenario generated by this model is graphed in Figure (ref).

For solving the FBSDE system (ref) we generate $2,000$ scenarios (Monte Carlo paths) using the low variance spectral decomposition method. Each scenario is based on daily sampling, i.e. $\frak{n}=365$. For calculating the conditional expectations in the Longstaff-Schwartz regression we use the cutoff value of $K=4$. The expectation behind this choice is that we obtain the accuracy of four sigma.

We attempt to construct a numerical solution to (ref) through the following iterative procedure. We start by generating initial Monte Carlo paths $X_t^{(0)}$ using the (discretized version of the) raw process (ref). For the initial guess of $Y_t$ we take $Y_t^{(0)}=\nabla G(X_t^{(0)})$. Notice that this is not at solution to the backward equations in (ref), it merely satisfies the terminal condition. After this, we iterate

for $k=0,1,\ldots$, until the stopping criterion is satisfied or the maximum number of iterations is reached. For the stopping criterion we choose the condition that the average $L^2$-norm change of a Monte Carlo path falls below a tolerance threshold of $10^{-8}$.

We make various choices of the coefficients $L,M,N$ defining the running cost functions. Depending on the values of these coefficients, the iterative process described above converges to a meaningful solution or it diverges. At this point, it is unclear what choices of the coefficients lead to what outcomes, but it appears that there are well defined basins of convergence in the space of the parameters.

We first consider a high cost policy. Figures (ref) - (ref) show the graphs of the solutions assuming the following parameters of the quadratic polynomial in the running cost function $c_2(u_2,x_2)$:

Unlike Figure (ref), the curves in all graphs below show the averages of the corresponding quantities over 2,000 Monte Carlo paths.

Under this running cost function, the optimal policy is to implement a draconian isolation regime, which leads to a rapid drop in infections, while keeping the susceptible fraction of the population at a very high level.

On the other hand, Figures (ref) - (ref) show the plots of the solutions assuming a low cost policy, with the following parameters of the quadratic polynomial in the running cost function $c_2(u_2,x_2)$:

Under this running cost function, the optimal policy is a moderate isolation regime. Following this policy, the isolation rate is high, as the infections are low, and then it declines over time while the epidemic develops. Unlike the policy above, this leads to a gradual decline in both the infected and susceptible fractions of the population.

Consider now the case of an optimal vaccination strategy. Again, we make two choices of the running cost function: high cost and low cost.

For the high cost case we choose the parameters as follows:

The results of Monte Carlo simulations for this cost function are plotted in Figures (ref) - (ref). They parallel the results presented above in the case of high cost isolation mitigation. The optimal policy is a massive vaccination campaign that dramatically reduces the susceptible fraction of the population and leads to significantly lower infections.

The low cost running cost function is parameterized as follows:

The results of Monte Carlo simulations for this cost function are plotted in Figures (ref) - (ref).

Finally, we consider the case of the vaccination / isolation policy. An example of a cost running cost function is given by the following parameter values:

The results of Monte Carlo simulations are plotted in Figures (ref) - (ref).

A low cost running cost function is given by the following set of coefficients:

The results of Monte Carlo simulations are plotted in Figures (ref) - (ref).