Min-cost Vaccination

Lecture 1 of IRC, LIAM, Fields-CQAM MfPH Training Series

2020-03-09

@York University

Model without vaccination

\begin{align*} S^{\prime}(t) &=& -\beta I S\\ I^{\prime}(t) &=& \beta I S - \gamma I \tag{2.1}\\ R^{\prime}(t) &=& \gamma I\\ S(0) &=& S_0\\ I(0) &=& I_0\\ R(0) &=& 0\\ \end{align*}

Relative removal rate

We divide the second equation of (2.1) by the first to get

\begin{align*} \frac{dI}{dS} &=& -1 + \frac{\rho}{S} \tag{2.2}\\ \end{align*}

where

\(\rho = \frac{\gamma}{\beta}\) is the relative removal rate.

Solution

(2.2) has the following solution form:

\begin{equation*} S + I - \rho\ln(S) = C \end{equation*}

where the constant \(C=1 - \rho\ln(S_0)\) since \(S_0+I_0=1\), hence

\begin{equation*} I = 1 - S + \rho\ln(S/S_0) \tag{2.3} \end{equation*}

SI-plane

Peak infection

From (2.3) it is clear at \(S=\rho\) the infection achieves its peak:

\begin{equation*} I_{max} = 1 - \rho + \rho\ln(\rho/S_0) \end{equation*}

Asymptotics

\(S(\infty)\) exists since \(S(t)\) is monotone and bounded.
\(I(\infty)=0\)
hence (2.3) gives \(0 = 1 - S(\infty) + \rho\ln(S(\infty)/S_0)\) which can be solved for \(S(\infty)\)
\(R(\infty)\) alone, the intensity, could indicate the severity of an epidemic, since \(I(\infty)=0\) and \(S(\infty)=1-R(\infty)\)
For \(S(\infty)\), \(R(\infty)\), only \(\rho\) matters.

Isolation effect

if we increase the removal rate \(\gamma\) by isolating some infected
the relative removal rate \(\rho\) increases
the peak infection \(I_{max}\) decreases since \(\rho<S_0\)
\(S(\infty)\) increases
\(R(\infty)\) decreases since \(R(\infty)+S(\infty)=1\)

Vaccination effect

Qualitatively vaccination is equivalent to jump from one solution path in Fig.1 to a lower solution path.

Temporal views

\begin{align*} \frac{dS}{dt} &=& -\beta S(t) [1-S(t)+\rho \ln (S(t)/S_0) ]\\ I(t) &=& 1 - S(t) + \rho \ln (S(t)/S_0)\\ R(t) &=& -\rho \ln [S(t)/S_0]\\ \end{align*}

I(t)

R(t)

S(t)

Model with vaccination

\begin{align*} S^{\prime}(t) &=& -\beta I S - \alpha\\ I^{\prime}(t) &=& \beta I S - \gamma I \tag{3.1}\\ R^{\prime}(t) &=& \gamma I\\ V^{\prime}(t) &=& \alpha \\ S(0) &=& S_0\\ I(0) &=& I_0\\ R(0) &=& 0\\ V(0) &=& 0\\ \end{align*}

V(t)

Vaccination function

Consider a finite control period \(T=15\Delta\) instead of \(\infty\), and a stepwise vaccination function instead of a continuous one,

\begin{align*} \alpha(t)= M_j d \;\; \forall t \in \small{\left[ j\Delta, (j+1)\Delta \right) } \\ M_j \in \small{\{0,\cdots,9\}}, j \in \small{\{0,\cdots,14\}} \tag{5.1} \end{align*}

Control requirements

\begin{align*} \forall t \in [0,T] \\ R(t) + I(t) & \le & A\\ I(t) & \le & B\\ \end{align*}

Optimization model

\begin{align*} \min_{\substack{M_j \\ j=1,\cdots,14}} & \sum_{j=1}^{14} (M_j d \Delta)^2\\ & R(j\Delta) + I(j\Delta)\le A \;\;\forall j\\ & I(j\Delta) \le B \;\; \forall j\\ & M_j \in \small{\{0,\cdots,9\}} \end{align*}

Typical Solution

\begin{equation*} [3,2,2,2,2,1,2,1,0,0,0,0,0,0] \end{equation*}

Table 1

Dynamic programming

State transition \(\Gamma(\cdot)\) is given by the first two equations and initial conditions in (3.1)

\begin{equation*} (S_{j+1}, I_{j+1}) = \Gamma (S_j,I_j,M) \end{equation*}

\(\Gamma(\cdot)\) shall be rounded to integer in DP after Runte-Kutta.

Terminal reward

Terminal reward function involves only one forward step,

\begin{equation*} f_1(S,I) = \min_{M \in \small{\{0,\cdots,9\}}} C(M) \tag{A.1} \end{equation*}

where \(C(M) = (M d \Delta)^2\) if feasible else \(\infty\)

Bellman equation

\begin{equation*} f_{k+1}(S,I) = \min_{M \in \small{\{0,\cdots,9\}}} \left( C(M) + f_k(\Gamma(S,I,M)) \right) \tag{A.2} \end{equation*}

Memoization technique

#setup, see Table 1
T   = 14
rho = 0.1

@functools.lru_cache()
def Gamma(S,I,M): #state transition
    pass          #Runte-Kutta Eq 3.1

@functools.lru_cache()
def f(k,S,I):     #cost table
    if k==1:
        pass      # (A.1)
    else:
        pass      # (A.2)

I0 = 0.02
f(14, 1 - I0, S0) # seek optimal schedul at I0

I0 = 0.04
f(14, 1 - I0, S0) # seek optimal schedule at I0

Optimal solution

\begin{equation*} f_{14}(S_0, I_0) \end{equation*}

Cost

Table 1

Discretization error

With the found control, we may

re-calculate the dynamic system to see state difference,
check violation of the path and point constraints.

One might need to adaptively refine the grid and start over.

Questions

Homework

implement the DP algorithm using memoization and lazy evaluation techniques;
reproduce the Table 1;
reproduce Fig 3-8.
due March 9 in class.

Next lecture

We will solve the problem as a continuous optimal control problem by collocation, single shooting and multiple shooting algorithms.

Bibliography

HERBERT W. HETHCOTE AND PAUL WALTMAN, Optimal Vaccination Schedules in a Deterministic Epidemic Model, MATHEMATICAL BIOSCIENCES 18, 365-381 (1973)