R语言deSolve包SIR模型实现7天周期触发疫苗接种逻辑求助

核心问题说明 #

原有代码的问题在于deSolve的ODE求解函数默认是无状态的，无法记录当前是否处于接种周期、已接种天数这类过程信息，因此需要新增辅助状态变量跟踪接种状态，同时使用离散欧拉法匹配按天触发的规则。

修改后完整可运行代码 #

library(deSolve)

sirV <- function(time, y, params) {
  # 状态变量：S易感, I感染, R康复, V接种, vac_active是否在接种周期, vac_days_passed当前周期已接种天数
  S <- y[1]
  I <- y[2]
  R <- y[3]
  V <- y[4]
  vac_active <- y[5]
  vac_days_passed <- y[6]
  
  with(as.list(params), {
    # 接种逻辑判断
    vac_helper <- 0
    new_vac_active <- vac_active
    new_vac_days <- vac_days_passed
    
    if (vac_active == 0) {
      # 非接种周期，校验触发条件
      if (I > 100 & S >= 50) {
        new_vac_active <- 1
        new_vac_days <- 1
        vac_helper <- 50
      }
    } else {
      # 接种周期内
      if (S < 50) {
        # 提前终止条件
        new_vac_active <- 0
        new_vac_days <- 0
        vac_helper <- 0
      } else {
        if (vac_days_passed < 7) {
          # 未完成7天接种，继续执行
          new_vac_days <- vac_days_passed + 1
          vac_helper <- 50
        } else {
          # 完成7天接种，关闭当前周期，下一轮重新判断触发条件
          new_vac_active <- 0
          new_vac_days <- 0
          vac_helper <- 0
        }
      }
    }
    
    # 微分方程计算
    N <- S + I + R + V
    dS <- -S * beta * I / N - vac_helper
    dI <- S * beta * I / N - gamma * I
    dR <- gamma * I
    dV <- vac_helper
    # 辅助状态变量为离散更新，导数设为0，实际值通过事件更新
    d_vac_active <- 0
    d_vac_days <- 0
    
    return(list(
      c(dS, dI, dR, dV, d_vac_active, d_vac_days),
      new_vac_active = new_vac_active,
      new_vac_days = new_vac_days
    ))
  })
}

# 事件函数：每时间步更新接种状态变量
update_vac_state <- function(time, y, params) {
  out <- sirV(time, y, params)
  y[5] <- out$new_vac_active
  y[6] <- out$new_vac_days
  return(y)
}

# 参数与初始值设置
myparameters <- c(gamma = 1/10, beta = 0.2)
times <- seq(0, 300, by = 1)
# 初始状态新增两个辅助变量，初始值均为0
init_state <- c(S = 100000, I = 10, R = 0, V = 0, vac_active = 0, vac_days_passed = 0)

# 求解ODE，使用离散欧拉法匹配按天计算的规则
my_ode <- as.data.frame(ode(
  y = init_state,
  times = times,
  func = sirV,
  parms = myparameters,
  method = "euler",
  events = list(func = update_vac_state, time = times)
))

# 查看结果时可剔除辅助变量
head(my_ode[,1:5])

逻辑匹配说明 #

• 完全符合触发条件：仅当感染人数I>100且易感人数S≥50时启动新的接种周期
• 接种周期内自动连续执行7天，中间I低于阈值不会中断，仅S<50时提前终止
• 每轮7天接种完成后自动重新校验触发条件，符合要求则开启下一轮接种

内容的提问来源于stack exchange，提问作者Parinn