如何在lsoda中捕获状态变量满足条件的首个时间点？

针对你基于deSolve的模型需求，以下是两种可行方案，其中根查找方法为最优解：

方法一：使用deSolve内置根查找功能（推荐） #

deSolve的根查找机制能精准捕获状态变量达到阈值的时间点，无需事后插值，且保证模型在触发点的连续性。

1. 定义根检测函数 #

该函数返回I - 阈值的值，当I达到设定阈值时，此值会穿过零点，lsoda会自动检测这个触发点：

root_func <- function(t, x, parms) {
  # 从参数中读取I的阈值
  I_threshold <- parms["I_threshold"]
  # 返回I与阈值的差值，零点即为触发点
  return(x[3] - I_threshold)
}

2. 修改模型函数 #

调整Q的计算逻辑，引入I_change_time参数（初始为NA），仅当该时间点被捕获后，才触发Q的非零值：

model <- function(t, x, parms) {
  S <- x[1]
  E <- x[2]
  I <- x[3]
  R <- x[4]
  K <- x[5]
  V <- x[6]
  
  with(as.list(parms), {
    # 仅当I_change_time已确定，且t处于触发窗口内时，Q=0.04
    Q <- ifelse(!is.na(I_change_time) & t >= I_change_time + day0 & t <= I_change_time + day0 + duration, 0.04, 0)
    
    dS <- -B * S * I - Q * S
    dE <- B * S * I - r * E
    dI <- r * E - g * I
    dR <- g * I 
    dK <- r * E
    dV <- Q * S
    res <- c(dS, dE, dI, dR, dK, dV)
    list(res)
  })
}

3. 设置事件并运行模型 #

通过events参数定义根触发时的动作：将当前时间赋值给I_change_time，然后运行lsoda：

library(deSolve)

# 定义参数，加入I阈值和初始为NA的I_change_time
pp <- c(day0 = 30, duration = 5, B = 0.04, r = 1/7, g = 1/7, I_threshold = 5, I_change_time = NA)
init <- c(S = 99, E = 1, I = 0, R = 0, K = 0, V = 0)
mtime <- 120
step_size <- 0.2

# 根触发时执行的事件：更新I_change_time为当前时间
event_func <- function(t, x, parms) {
  parms["I_change_time"] <- t
  return(parms)
}

# 运行模型，启用根查找和事件触发
output <- as.data.frame(lsoda(
  y = init,
  times = seq(0, mtime, step_size),
  func = model,
  parms = pp,
  rootfun = root_func,
  events = list(func = event_func, root = TRUE),
  rtol = 1e-6, atol = 1e-6 # 提高精度确保捕获准确时间点
))

# 查看捕获的时间点
cat("I首次达到阈值的时间点：", pp["I_change_time"], "\n")

方法二：事后插值查找时间点（替代方案） #

如果不想使用根查找，可以先运行无Q触发的模型，再通过插值找到阈值时间点，最后重新运行完整模型：

1. 运行初始模型获取I的时间序列 #

# 先运行不带Q触发的基础模型
model_initial <- function(t, x, parms) {
  S <- x[1]
  E <- x[2]
  I <- x[3]
  R <- x[4]
  K <- x[5]
  V <- x[6]
  
  with(as.list(parms), {
    dS <- -B * S * I
    dE <- B * S * I - r * E
    dI <- r * E - g * I
    dR <- g * I 
    dK <- r * E
    dV <- 0
    res <- c(dS, dE, dI, dR, dK, dV)
    list(res)
  })
}

output_initial <- as.data.frame(lsoda(
  y = init,
  times = seq(0, mtime, step_size),
  func = model_initial,
  parms = pp
))

2. 线性插值计算I_change_time #

I_threshold <- 5
# 找到I首次超过阈值的前后索引
idx <- which(output_initial$I >= I_threshold)[1]

if (!is.na(idx)) {
  # 提取前后时间点和对应的I值
  t_prev <- output_initial$time[idx-1]
  t_curr <- output_initial$time[idx]
  I_prev <- output_initial$I[idx-1]
  I_curr <- output_initial$I[idx]
  # 线性插值计算精确时间点
  I_change_time <- t_prev + (I_threshold - I_prev) * (t_curr - t_prev) / (I_curr - I_prev)
} else {
  # 若整个模拟期I未达阈值，设为无穷大
  I_change_time <- Inf
}

3. 用捕获的时间点重新运行完整模型 #

# 更新参数加入I_change_time
pp_updated <- c(pp, I_change_time = I_change_time)

# 定义带Q触发逻辑的模型
model_updated <- function(t, x, parms) {
  S <- x[1]
  E <- x[2]
  I <- x[3]
  R <- x[4]
  K <- x[5]
  V <- x[6]
  
  with(as.list(parms), {
    Q <- ifelse(t >= I_change_time + day0 & t <= I_change_time + day0 + duration, 0.04, 0)
    
    dS <- -B * S * I - Q * S
    dE <- B * S * I - r * E
    dI <- r * E - g * I
    dR <- g * I 
    dK <- r * E
    dV <- Q * S
    res <- c(dS, dE, dI, dR, dK, dV)
    list(res)
  })
}

# 运行最终模型
output_final <- as.data.frame(lsoda(
  y = init,
  times = seq(0, mtime, step_size),
  func = model_updated,
  parms = pp_updated
))

方案对比 #

• 根查找方法：精度高，无需二次模拟，能实时响应阈值触发，是最优选择。
• 事后插值方法：实现简单，但依赖时间步长精度，步长越大误差越大，且需要两次模拟。

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