Yofre H. Garcia, UNACH

Saúl Diaz-Infante Velasco

Jesús Adolfo Minjárez Sosa, UNISON

saul.diazinfante@unison.mx

https://slides.com/sauldiazinfantevelasco/sym_psp_xv

Argument. When there is a shortage of vaccines, sometimes the best response is to refrain from vaccination, at least for a while.

Hypothesis. Under these conditions, inventory management suffers significant random fluctuations

Objective. Optimize the management of vaccine inventory and its effect on a vaccination campaign

On October 13 2020, the Mexican government announced a vaccine delivery plan from Pfizer-BioNTech and other companies as part of the COVID-19 vaccination campaign.

Methods. Given a vaccine shipping schedule, we describe stock management with a backup protocol and quantify the random fluctuations due to a program under high uncertainty.

Then, we incorporate this dynamic into a system of ODE that describes the disease and evaluate its response.

G_t = R_{t+1} + R_{t+2} + R_{t+3} + \dots R_{T}

G_t = R_{t+1} + \gamma R_{t+2} + \gamma^2 R_{t+3} + \dots

\begin{aligned} %C_{t+1} &= C(x_{t}, a_{t}) \\ %\Phi_{t+1}^{h}(x_t,a_t) &= x_t + F(x_t, a_t, \xi_t) \end{aligned}

Agent

R_{t+1}

x_{t+1}

x_0,a_0,R_1,

a_t\in \mathcal{A}(x_t)

action

state

x_{t}

reward

R_{t}

x_0, a_0, R_1, x_1, a_1, R_2, \cdots, x_t, a_t, R_{t+1} \cdots, x_{T-1}, a_{T-1}, R_T, x_T

\begin{equation*} \begin{aligned} S'(t) &= -\beta IS \\ I'(t) &= \beta IS - \gamma I \\ R'(t) & = \gamma I \\ & S(0) = S_0, I(0)=I_0, R(0)=0 \\ & S(t) + I(t) + R(t )= 1 \end{aligned} \end{equation*}

\lambda_V:= \underbrace{ \textcolor{orange}{\Psi_V}}_{cte.} \cdot \ S(t)

\begin{equation*} \begin{aligned} S'(t) &= -\beta IS - \textcolor{red}{\lambda_V(t)} \\ I'(t) &= \beta IS - \gamma I \\ R'(t) & = \gamma I \\ V'(t) & = \textcolor{red}{\lambda_V(t)} \\ & S(0) = S_0, I(0)=I_0, \\ &R(0)=0, V(0) = 0 \\ & S(t) + I(t) + R(t) + V(t)= 1 \end{aligned} \end{equation*}

\begin{equation*} \begin{aligned} S'(t) &= -\beta IS - \textcolor{red}{\lambda_V(x, t)} \\ I'(t) &= \beta IS - \gamma I \\ R'(t) & = \gamma I \\ V'(t) & = \textcolor{red}{\lambda_V(x,t)} \end{aligned} \end{equation*}

\lambda_V:= \underbrace{ \textcolor{orange}{\Psi_V}}_{cte.} \cdot \ S(t)

\begin{aligned} S'(t) &= \cancel{-\beta IS} - \underbrace{\textcolor{red}{\lambda_V(x, t)}}_{=\Psi_V S(t)} \\ I'(t) &= \cancel{\beta IS} - \cancel{\gamma I} \\ {R'(t)} & = \cancel{\gamma I} \\ V'(t) & = \underbrace{\textcolor{red}{\lambda_V(x, t)}}_{=\Psi_V S(t)} \\ & S(0) \approx 1, I(0) \approx 0, \\ &R(0)\approx 0, V(0) = 0 \\ & S(t) + \cancel{I(t)} + \cancel{R(t)} +\cancel{ V(t)}= 1 \end{aligned}

\begin{aligned} S'(t) &= - \Psi_V S(t) \\ V'(t) & = \Psi_V S(t) \\ & S(0) \approx 1, V(0) = 0 \\ & S(t) + V(t) \approx 1 \end{aligned}

The effort invested in preventing or mitigating an epidemic through vaccination is proportional to the vaccination rate

\textcolor{orange}{\Psi_V}

Let us assume at the beginning of the outbreak:

\textcolor{orange}{S(0)\approx 1}

\begin{aligned} S'(t) &= - \Psi_V S(t) \\ V'(t) & = \Psi_V S(t) \\ & S(0) \approx 1, V(0) = 0 \\ & S(t) + V(t) \approx 1 \end{aligned}

S(t) \approx N \exp(-\Psi_V t)

\begin{aligned} V(t) \approx \cancel{V(0)} + \int_0^t \Psi_V S(\tau) d\tau \end{aligned}

\begin{aligned} V(t) &\approx \cancel{V(0)} + \int_0^t \Psi_V S(\tau) d\tau \\ & \approx \int_0^t N \Psi_V \exp(-\Psi_V \tau) d\tau \\ & \approx N \exp(-\Psi_V \tau) \mid_{\tau=t}^{\tau=0} \\ &= \cancel{N} (1 - \exp(-\Psi_V t)) \end{aligned}

Then we estimate the number of vaccines with

X_{vac}(t) := \int_{0}^t \Psi_V \Bigg( \underbrace{S(\tau)}_{ \substack{ \text{target} \\ \text{pop.} } }\Bigg )

Then, for a vaccination campaign, let:

\begin{aligned} \textcolor{orange}{T}: &\text{ time horizon} \\ \textcolor{green}{ X_{cov}:= X_{vac}(} \textcolor{orange}{T} \textcolor{green}{)}: &\text{ covering at time $T$} \end{aligned}

\begin{aligned} X_{cov} = &X_{vac}(T) \\ \approx & 1 - \exp(-\textcolor{teal}{\Psi_V} T). \\ \therefore & \textcolor{teal}{\Psi_V} \approx -\frac{1}{T} \ln(1 - X_{cov}) \end{aligned}

Then we estimate the number of vaccines with

X_{vac}(t) := \int_{0}^t \Psi_v \Bigg( \underbrace{S(\tau)}_{ \substack{ \text{target} \\ \text{pop} } }\Bigg )

Then, for a vaccination campaign, let:

\begin{aligned} \textcolor{orange}{T}: &\text{ horizon time} \\ \textcolor{green}{ X_{cov}:= X_{vac}(} \textcolor{orange}{T} \textcolor{green}{)}: &\text{ Coverage at time $T$} \end{aligned}

\begin{aligned} X_{cov} = &X_{vac}(T) \\ \approx & 1 - \exp(-\textcolor{teal}{\Psi_V} T). \\ \therefore & \textcolor{teal}{\Psi_V} = -\frac{1}{T} \ln(1 - X_{cov}) \end{aligned}

Estimated population of Hermosillo, Sonora in 2024 is 930,000.

So to vaccinate 70% of this population in one year:

\begin{aligned} 930,000 &\times 0.00329 \\ \approx & 3,059.7 \ \text{jabs/day} \end{aligned}

Base Model

\begin{equation*} \begin{aligned} \frac{dS}{dt} &= \mu N - (\mu + \psi_V) S + \omega_V V + \delta_R R, \\ \frac{dE}{dt} &= \lambda_f S + (1 - \varepsilon)\lambda_f V - (\mu + \delta_E) E, \\ \frac{dI_S}{dt} &= \rho \delta_E E - (\mu + \alpha_S) I_S, \\ \frac{dI_A}{dt} &= (1 - \rho)\delta_E E - (\mu + \alpha_A) I_A, \\ \frac{dR}{dt} &= (1 - \theta)\alpha_S I_S + \alpha_A I_A - (\mu + \delta_R) R, \\ \frac{dD}{dt} &= \theta \alpha_S I_S, \\ \frac{dV}{dt} &= \psi_V S - \big[(1 - \varepsilon)\lambda_f + \mu + \omega_V\big] V, \\ X_{\text{vac}} &= \psi_V (S + E + I_A + R), \\ \hat{N} &= S + E + I_S + I_A + R, \quad \hat{N} + D = 1, \\ \lambda_f &= \frac{\beta_S I_S + \beta_A I_A}{\hat{N}}. \end{aligned} \end{equation*}

After discretizise with a NSFDS

\begin{equation*} \begin{aligned} S^{t_{n+1}^{(k-1)}} &= \frac{ (1 - \varphi(h_k)\mu)\, S^{t_{n}^{(k-1)}} + \varphi(h_k)\!\left[ \mu \hat{N}^{t_{n}^{(k-1)}} + \omega_V V^{t_{n}^{(k-1)}} + \delta_R R^{t_{n}^{(k-1)}} \right] }{ \bigl( 1 + \varphi(h_k) \bigr)\, \bigl( \lambda_f^{t_{n}^{(k-1)}} + a_{k-1} \bigr) }, \\ E^{t_{n+1}^{(k-1)}} &= \frac{ \bigl(1 - \mu\,\varphi(h_k)\bigr) E^{t_{n}^{(k-1)}} + \varphi(h_k)\,\lambda_f^{t_{n}^{(k-1)}} \!\left[ S^{t_{n+1}^{(k-1)}} + (1 - \epsilon)\, V^{t_{n}^{(k-1)}} \right] }{ 1 + \varphi(h_k)\,\delta_E }, \\ I_S^{t_{n+1}^{(k-1)}} &= \frac{ \bigl(1 - \varphi(h_k)\mu\bigr)\, I_S^{t_{n}^{(k-1)}} + \varphi(h_k)\, p\,\delta_E\, E^{t_{n+1}^{(k-1)}} }{ 1 + \varphi(h_k)\,\alpha_S }, \\ I_A^{t_{n+1}^{(k-1)}} &= \frac{ \bigl(1 - \varphi(h_k)\mu\bigr)\, I_A^{t_{n}^{(k-1)}} + \varphi(h_k)\, (1 - p)\,\delta_E\, E^{t_{n+1}^{(k-1)}} }{ 1 + \varphi(h_k)\,\alpha_A }, \\ R^{t_{n+1}^{(k-1)}} &= \bigl[ 1 - \varphi(h_k)\,(\mu + \delta_R) \bigr] R^{t_{n}^{(k-1)}} + \varphi(h_k)\!\left[ (1 - \theta)\,\alpha_S\, I_S^{t_{n+1}^{(k-1)}} + \alpha_A\, I_A^{t_{n+1}^{(k-1)}} \right], \\ D^{t_{n+1}^{(k-1)}} &= D^{t_{n}^{(k-1)}} + \varphi(h_k)\, \theta\, \alpha_S\, I_S^{t_{n+1}^{(k-1)}}, \\ V^{t_{n+1}^{(k-1)}} &= a_{k-1}\, \varphi(h_k)\, S^{t_{n+1}^{(k-1)}} + V^{t_{n}^{(k-1)}} \!\left( 1 - \varphi(h_k)\!\left[ (1 - \epsilon)\, \lambda_f^{t_{n}^{(k-1)}} + \mu + \omega_V \right] \right), \\ X_{\text{vac}}^{t_{n+1}^{(k-1)}} &= \varphi(h_k)\, \psi_V^{(k-1)}\!\left( S^{t_{n}^{(k-1)}} + E^{t_{n}^{(k-1)}} + I_A^{t_{n}^{(k-1)}} + R^{t_{n}^{(k-1)}} \right) + X_{\text{vac}}^{t_{n}^{(k-1)}}, \\ \lambda_{f}^{t_n^{(k-1)}} &:= \frac{ \beta_S I_{S}^{t_n^{(k-1)}} + \beta_A I_{A}^{t_n^{(k-1)}} }{\widehat{N}^{t_n^{(k-1)}}}, \\ \widehat{N}^{t_n} &= S^{t_n^{(k-1)}} + E^{t_n^{(k-1)}} + I_S^{t_n^{(k-1)}} + I_A^{t_n^{(k-1)}} + R^{t_n^{(k-1)}} + V^{t_n^{(k-1)}}, \\ & a_{k-1} = p_i \cdot \psi_V^{(k-1)}, \quad p_i \in [0,1], \qquad i = 1, \dots, l, \quad k = 1, \ldots, M. \\ \\ &X_{\bullet}^{t_0^{(0)}} :=( S^{t_0^{(0)}}, E^{t_0^{(0)}}, I_S^{t_0^{(0)}}, I_A^{t_0^{(0)}}, R^{t_0^{(0)}}, V^{t_0^{(0)}}, K_{\text{vac}}^{t_0^{(0)}}, T^{t_0^{(0)}}, L^{t_0^{(0)}} )^{\top} \\ &X_{\bullet}^{t_0^{(k-1)}} = X_{\bullet}^{t_{N_{k-1}}^{(k-2)}}, \qquad k=2, \dots, M-1. \end{aligned} \end{equation*}

\begin{equation*} \begin{aligned} X_{vac}^{t_{n+1}^{(k)}} &= \varphi(h_k)\psi_V^{(k)} \left( S^{t_{n}^{(k)}} + E^{t_{n}^{(k)}} + I_A^{t_{n}^{(k)}} + R^{t_{n}^{(k)}} \right) + X_{vac}^{t_{n}^{(k)}}, \\ \widetilde{K}^{t_{{n+1}}^{(k-1)}} = & K^{t_n^{(k-1)}} - \underbrace{ \left( X_{\text{vac}}^{t_{n+1}^{(k-1)}} - X^{t_0^{(k-1)}}_{\text{vac}} \right) }_{ \substack{ \text{Coverage } \\ \text{during} \\ \text{period } k } } \\ \\ K^{t_{n+1}^{(k-1)}} &= \max \left\{ 0, \widetilde{K}^{t_{{n+1}}^{(k-1)}} \right\} \end{aligned} \end{equation*}

\begin{align*} L^{t_{n+1}^{(k)}} &= K^{t_{n}^{(k)}} \cdot \left( 1 - e^{ - \textcolor{orange}{\lambda(T^{t_{n}^{(k)}})} \cdot t_{n}^{(k)} } \right), \\ \lambda(T^{t_{n+1}^{(k)}}) &= \kappa \cdot \max\{0, \ \textcolor{orange}{T^{t_{n}^{(k)}}} - \mu_{T}\}, \end{align*}

Stock degradation due to Temperature

Sources of uncertainty

Random fluctuations in the stock replenishment
Stochastic degradation due to temperature

Replenishment random perturbation

order size
time-delivery

Stock degradation due to Temperature

d\textcolor{orange}{T_t} = \theta_T (\mu_T - T_t) dt + \sigma_T dW_t

DALY (Disability-Adjusted Life Years)

For a given cause (c), sex (s), age (a), and year (t).

🩸 YLL — Years of Life Lost

Measures years of life lost due to premature death.

N(c, s, a, t) — number of deaths due to cause c
L(s, a) — standard life expectancy at age a for sex s

♿ YLD — Years Lived with Disability

Measures years lived with disease or disability.

I(c, s, a, t) — incident cases for cause c
DW(c, s, a) — disability weight for cause c
L(c, s, a, t) — average duration of the case

YLL(c, s, a, t) = N(c, s, a, t) \times L(s, a)

YLD(c, s, a, t) = I(c, s, a, t) \times DW(c, s, a) \times L(c, s, a, t)

DALY measures the total burden of disease, combining:

🕯️ Mortality (YLL) — early death
🧩 Morbidity (YLD) — time lived with illness or disability

DALY(c, s, a, t) = YLL(c, s, a, t) + YLD(c, s, a, t)

\begin{aligned} x_{t_{n+1}} & = x_{t_n} + F(x_{t_n}, \theta, a_0) \cdot h, \quad x_{t_0} = x(0), \\ \text{where: }& \\ t_n &:= n \cdot h, \quad n = 0, \cdots, N, \quad t_{N} = T. \end{aligned}

\begin{aligned} &\min_{a_{0}^{(k)} \in \mathcal{A}_0} c(x_{\cdot}, a_0):= c_1(a_0^{(k)})\cdot T^{(k)} + \sum_{n=0}^{N-1} c_0(t_n, x_{t_n}) \cdot h \\ \text{ s.t.} & \\ & x_{t_{n+1}} = x_{t_n} + F(x_{t_n}, \theta, a_0) \cdot h, \quad x_{t_0} = x(0), \\ \text{where: }& \\ t_n &:= n \cdot h, \quad n = 0, \cdots, N, \quad t_{N} = T. \end{aligned}

\begin{aligned} C(x_{t^{(k+1)}}, & a_{t^{(k+1)}}) = \\ & C_{YLL}(x_{t^{(k+1)}},a_{t^{(k+1)}}) \\ +& C_{YLD}(x_{t^{(k+1)}},a_{t^{(k+1)}}) \\ +& C_{stock}(x_{t^{(k+1)}},a_{t^{(k+1)}}) \\ +& C_{campaign}(x_{t^{(k+1)}},a_{t^{(k+1)}}) \end{aligned}

J(x,\pi) = E\left[ \sum_{k=0}^M C(x_{t^{(k)}},a_{t^{(k)}}) | x_{t^{(0)}} = x , \pi \right]

\begin{aligned} C_{YLL}(x_{t^{(k+1)}},a_{t^{(k+1)}}) &= \int_{t^{(k)}}^{t^{(k+1)}} YLL dt, \\ C_{YLD}(x_{t^{(k+1)}},a_{t^{(k+1)}}) &= \int_{t^{(k)}}^{t^{(k+1)}} YLD(x_t, a_t) dt \\ YLL(x_t, a_t) &:= m_1 p \delta_E (E(t) - E^{t^{(k)}} ), \\ YLD(x_t, a_t) &:= m_2 \theta \alpha_S(E(t) - E^{t^{(k)}}), \\ t &\in [t^{(k)},t^{(k + 1)}] \end{aligned}

\begin{aligned} C_{stock}(x_{t^{(k+1)}},a_{t^{(k+1)}}) & = \int_{t^{(k)}}^{t^{(k+1)}} m_3(K_{Vac}(t) - K_{Vac}^{t^{(k)}}) dt \\ C_{campaign}(x_{t^{(k+1)}},a_{t^{(k+1)}}) &=\int_{t^{(k)}}^{t^{(k+1)}} m_4(X_{vac}(t) - X_{vac}^{t^{(k)}}) dt \end{aligned}

Q-learning for multiple stages and uncertainty in policies and other terms related to epidemics.

Move to Reinforcement Learning

Onésimo Hernández-Lerma, Leonardo R. Laura-Guarachi, Saul Mendoza-alacios, David González-Sánchez

An Introduction to Optimal Control Theory

The Dynamic Programming Approach

978-3-031-21138-6

Reinforcement learning and optimal control

Bertsekas, Dimitri P.

Athena Sci. Optim. Comput. Ser.

Athena Scientific, Belmont, MA, 2019, xiv+373 pp.

ISBN: 978-1-886529-39-7

Powell, Warren B.

Reinforcement Learning and Stochastic Optimization: A Unified Framework for Sequential Decisions. United Kingdom: Wiley, 2022.

ISBN: 978-1-119-81505-1

Reinforcement Learning: An Introduction
Richard S. Sutton and Andrew G. Barto
Publisher: MIT Press
Year: 2018 (Second Edition)
ISBN: 978-0-262-03924-6

https://github.com/SaulDiazInfante/rl_vac.jl

https://slides.com/sauldiazinfantevelasco/sym_psp_xv

GRACIAS!!