Skip to contents

Epidemic control based-on deaths

Sandor Beregi

2026-04-11

control_est_D.Rmd

Overview

This vignette demonstrates model predictive optimal control of deaths for COVID-19 based on death data only

Load Libraries and Setup

Load the required libraries and initialize the simulation environment.

library(EpiControl)
library(VGAM)
library(parallel)
library(pbapply)
library(zoo)  # For rolling sum operations
library(dplyr)
library(ggplot2)

Initialize the cluster for parallel computing:

cores <- detectCores() - 1
cl <- makeCluster(cores)

clusterSetRNGStream(cl, iseed = 20250501)

Epidemiological and Noise Parameters

Define parameters for pathogens and noise: Can include different pathogens, in this case we consider COVID-19 and Ebola virus.

Epi_pars <- data.frame(
  Pathogen = c("COVID-19", "Ebola"),
  R0 = c(2.2, 2.5),
  gen_time = c(6.5, 15.0),
  gen_time_var = c(2.1, 2.1),
  CFR = c(0.0132, 0.5),
  mortality_mean = c(10.0, 10.0),
  mortality_var = c(1.1, 1.1)
)

Noise_pars <- data.frame (
  repd_mean = 14.5, #Reporting delay mean 10.5
  del_disp = 1.1, #Reporting delay variance 5.0
  ur_mean = 0.3, #Under reporting, mean/variance
  ur_beta_a = 50.0
)

Action Space for Policy Interventions

Define non-pharmaceutical interventions: We have 2 actions, introducing a “lockdown” and no interventions.

Action_space <- data.frame (
  NPI = c("No restrictions", "Social distancing", "Lockdown"),
  R_coeff = c(1.0, 0.5, 0.2), #R0_act = R0 * ctrl_states
  R_beta_a = c(0.0, 5.0, 5.0), #R0_act uncertainty
  cost_of_NPI = c(0.0, 0.01, 0.15)
)

Simulation Setup

Set up key parameters and initialize simulation data:

C_target <- 200
D_target <- 100
r_trans_len <- 7

sim_settings <- list(
  ndays = 61 * 7, #simulation length
  start_day = 1,
  N = 1e7, # population size
  I0 = 10, # initial infections
  C_target = C_target, #target cases
  C_target_pen = C_target*1.5, #overshoot penalty threshold
  R_target = 1.0,
  D_target = D_target, #one way to get peaks at 400 is to increase this to 15
  D_target_pen = 300, #max death
  #alpha = 0*1.3/C_target, #~proportional gain (regulates error in cases) covid
  alpha = 0*3.25/C_target, #~proportional gain (regulates error in cases) ebola
  alpha_d = 3.25/D_target,
  ovp = 0*5.0, #overshoot penalty
  dovp = 5.0, #death overshoot penalty
  gamma = 0.95, #discounting factor
  n_ens = 10L, #MC assembly size for 4
  sim_ens = 20L, #assembly size for full simulation
  rf = 21L, #days 14
  R_est_wind = 5L, #rf-2 #window for R estimation
  pred_days = 28L,
  r_trans_steep = 1.5,  # Growth rate
  r_trans_len = r_trans_len,  # Number of days for the transition
  t0 = r_trans_len  / 2, # Midpoint of the transition
  pathogen = 2,
  susceptibles = 0,
  delay = 0,
  ur = 0,
  r_dir = 2
)

column_names <- c("days", "sim_id", "I", "Lambda", "C", "Lambda_C", "S", "Deaths", "Re", "Rew", "Rest", "R0est", "policy", "R_coeff")
episim_data <- data.frame(matrix(0, nrow = (sim_settings$ndays), ncol = length(column_names)))
colnames(episim_data) <- column_names

episim_data$policy <- rep(1, sim_settings$ndays)
episim_data$days <- 1:(sim_settings$ndays)
episim_data[1, ] <- c(1, 1, sim_settings$I0, sim_settings$I0, Noise_pars$ur_mean * sim_settings$I0, Noise_pars$ur_mean * sim_settings$I0, sim_settings$N - sim_settings$I0, 0, Epi_pars[1, "R0"], Epi_pars[1, "R0"], 1, 1, 1, 1)

episim_data_ens <- replicate(sim_settings$sim_ens, episim_data, simplify = FALSE)

for (ii in 1:sim_settings$sim_ens) {
  episim_data_ens[[ii]]$sim_id <- rep(ii, sim_settings$ndays)
}

Running the Simulation

Run simulations in parallel using the pblapply function:


episettings <- list(
  sim_function = Epi_MPC_run_wd,
  reward_function = reward_fun_wd,
  R_estimator = R_estim,
  noise_par = Noise_pars,
  epi_par = Epi_pars,
  actions = Action_space,
  sim_settings = sim_settings,
  parallel = TRUE
)


clusterExport(cl, ls())

parallel::clusterEvalQ(cl, {
  library(pbapply)
  library(VGAM)
  library(EpiControl)
})
#> [[1]]
#>  [1] "EpiControl" "VGAM"       "splines"    "stats4"     "pbapply"   
#>  [6] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
#> [11] "methods"    "base"      
#> 
#> [[2]]
#>  [1] "EpiControl" "VGAM"       "splines"    "stats4"     "pbapply"   
#>  [6] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
#> [11] "methods"    "base"      
#> 
#> [[3]]
#>  [1] "EpiControl" "VGAM"       "splines"    "stats4"     "pbapply"   
#>  [6] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
#> [11] "methods"    "base"      
#> 
#> [[4]]
#>  [1] "EpiControl" "VGAM"       "splines"    "stats4"     "pbapply"   
#>  [6] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
#> [11] "methods"    "base"      
#> 
#> [[5]]
#>  [1] "EpiControl" "VGAM"       "splines"    "stats4"     "pbapply"   
#>  [6] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
#> [11] "methods"    "base"      
#> 
#> [[6]]
#>  [1] "EpiControl" "VGAM"       "splines"    "stats4"     "pbapply"   
#>  [6] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
#> [11] "methods"    "base"      
#> 
#> [[7]]
#>  [1] "EpiControl" "VGAM"       "splines"    "stats4"     "pbapply"   
#>  [6] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
#> [11] "methods"    "base"      
#> 
#> [[8]]
#>  [1] "EpiControl" "VGAM"       "splines"    "stats4"     "pbapply"   
#>  [6] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
#> [11] "methods"    "base"      
#> 
#> [[9]]
#>  [1] "EpiControl" "VGAM"       "splines"    "stats4"     "pbapply"   
#>  [6] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
#> [11] "methods"    "base"      
#> 
#> [[10]]
#>  [1] "EpiControl" "VGAM"       "splines"    "stats4"     "pbapply"   
#>  [6] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
#> [11] "methods"    "base"      
#> 
#> [[11]]
#>  [1] "EpiControl" "VGAM"       "splines"    "stats4"     "pbapply"   
#>  [6] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
#> [11] "methods"    "base"

# Add the cluster object to your settings
episettings$cl <- cl

results <- epicontrol(episim_data_ens, episettings)
#> [1] 0
#>   |                                                          |                                                  |   0%

episim_data_ens <- results
stopCluster(cl)

#for (jj in 1:sim_ens) {
#  episim_data_ens[[jj]] <- Epi_MPC_run_wd(episim_data_ens[[jj]], Epi_pars, Noise_pars, Action_space, pred_days = pred_days, n_ens = n_ens, ndays = nrow(episim_data), R_est_wind = R_est_wind, pathogen = 1, susceptibles = 0, delay = 0, ur = 0, r_dir = 2, N = N)
#  setTxtProgressBar(pb,jj)
#}
#close(pb)

for (jj in 1:sim_settings$sim_ens) {
  episim_data_ens[[jj]] <- head(episim_data_ens[[jj]], -sim_settings$pred_days)
  episim_data_ens[[jj]]["D_roll"] <- rollsum(episim_data_ens[[jj]]["Deaths"], 7, fill = NA)
  episim_data_ens[[jj]]["I_roll"] <- rollsum(episim_data_ens[[jj]]["I"], 7, fill = NA)
  episim_data_ens[[jj]]["D_cum"] <- cumsum(episim_data_ens[[jj]]["Deaths"])
  episim_data_ens[[jj]]["I_cum"] <- cumsum(episim_data_ens[[jj]]["I"])
}

# Combine Simulation Results
combined_data <- do.call(rbind, episim_data_ens)

Plotting results 


combined_data <- combined_data %>%
  group_by(sim_id, policy) %>%
  arrange(days) %>%
  mutate(group = cumsum(c(1, diff(days) != 1))) %>%
  ungroup()

combined_data <- combined_data %>%
  arrange(sim_id, days, policy)

policy_labels <- c("1" = "No intervention", "2" = "Lockdown")

cls <- rep("grey", sim_settings$sim_ens)
cls[1] <- "red"

ggplot() +
  geom_line(data = subset(combined_data, sim_id != 1), aes(x = days, y = Deaths, color = as.factor(sim_id)), alpha = 0.3) +
  geom_line(data = subset(combined_data, sim_id == 1), aes(x = days, y = Deaths), color = "darkred", alpha = 1.0) +
  labs(x = "Days", y = "Deaths") +
  scale_color_manual(values = cls) +
  geom_hline(yintercept = D_target, linetype = "dashed", color = "purple", size=0.25)+
  guides(color = FALSE)