R：根据 ODE 中的状态改变参数

Question

请原谅我，因为这是我第一次在 SO 上发帖。如果您认为我需要格式化我的更改，请告诉我:)

我目前正在尝试运行一个引入疫苗接种的 SIR 模型，这种疫苗接种将在达到阈值（例如，50% 的人口）后停止，而感染仍在继续。

在我的隔室模型中，我有以下状态：

1. S
2. 我
3. R
4. P
5. 锡
6. 在
7. Rn

来自S的人将以α率进行疫苗接种，疫苗接种要么成功，要么失败（基于某种疫苗有效率）。

我知道 SO 上有解决方案，但我能找到的大多数是在某个时间停止参数，因此是与时间相关的参数。但是，我对状态相关参数感兴趣。

我目前的解决方案是使用疫苗接种率，α，是一个状态。一旦达到阈值（应该被根函数捕获），事件函数会将 α 设置为 0。

我尝试过使用 rootfunc 和 eventfunc，但老实说我不确定我是否正确执行此操作，因此出现错误。

在这里，我创建了根和事件函数。

root_func <- function(t, states, parameters) {
    with(as.list(c(states, parameters)),
         {
             return(P + Rn - 0.5 * (S + I + R + P + Sn + In + Rn))  
         })
    
}

event_func <- function(t, states, parameters) {
    with(as.list(c(states, parameters)),
         {
             alpha <- 0
             return(list(c(S, I, R, P, Sn, In, Rn, alpha)))
         })
}

接下来，我创建了我的 SIR 模型。

model_1_take <- function(t, states, parameters) {
    with(
        as.list(c(states, parameters)),
        {
            N <- S + I + R + 
                P + 
                Sn + In + Rn
            
            alpha <- ifelse(t >= 100, 0, alpha)
            
            dS <- (mu * N) - 
                (beta * S * (I + In) / N) - 
                (alpha * epsilon * S) - 
                (alpha * (1 - epsilon) * S) - 
                (nu * S)
            dI <- (beta * S * (I + In) / N) - 
                (gamma * I) - 
                (nu * I)
            dR <- (gamma * I) - 
                (alpha * R) - 
                (nu * R)
            dP <- (alpha * epsilon * S) - 
                (nu * P)
            dSn <- (alpha * (1 - epsilon) * S) - 
                (beta * Sn * (I + In) / N) - 
                (nu * Sn)
            dIn <- (beta * Sn * (I + In) / N) - 
                (gamma * In) - 
                (nu * In)
            dRn <- (alpha * R) +
                (gamma * In) -
                (nu * Rn)
            
            return(list(c(dS, dI, dR, dP, dSn, dIn, dRn)))
        }
    )
}

在这里，我正在初始化我的状态。请注意，α 位于状态中，尽管从技术上讲它应该是一个参数。

initN <- 5e6
initRn <- 0
initR <- 0
initIn <- 0
initI <- 1
initP <- 0
initSn <- 0
initS <- initN - initRn - initR - initI - initIn - initP - initSn

states_1_take <- c(
  S = initS, 
  I = initI, 
  R = initR,
  P = initP,
  Sn = initSn,
  In = initIn,
  Rn = initRn,
  alpha = 0.02
)

最后，我初始化我的参数。

R_0 <- 2
gamma <- 1 / 5
epsilon <- 0.75
rho <- 0.2
mu <- 1 / (365 * 70)
nu <- 1 / (365 * 70)
beta <- R_0 * (gamma + nu)

parameters_1_take <- c(
    beta = beta, 
    gamma = gamma,
    epsilon = epsilon,
    rho = rho, 
    mu = mu,
    nu = nu
)

当我想运行模型时，我使用以下代码。

out <- ode(
        y = states_1_take,
        times = times, 
        func = model_1_take,
        parms = parameters_1_take,
        method = "lsodes",
    )

我得到的错误是：

checkEventFunc(Eventfunc, times, y, rho) 中的错误： events$func() 返回的值的数量 (1) 必须等于初始条件向量的长度 (8)

非常感谢任何帮助！

Answer 1

要调试问题，只需使用默认求解器调用模型函数，无需事件来重现错误，例如：

times <- 0:100
out <- ode(
        y = states_1_take,
        times = times, 
        func = model_1_take,
        parms = parameters_1_take,
        method = "lsoda",
    )

它表示模型返回 7 个导数，而初始向量包含 8 个状态。

也可以单独调用模型函数：

model_1_take(0, states_1_take, parameters_1_take)

错误的原因是，

alpha

显然没有状态变量。解决这个问题的一种方法是在模型函数中添加

alpha

作为通常具有导数

dAlpha = 0

的状态：

model_1_take <- function(t, states, parameters) {
  # ...
  dAlpha <- 0
            
  return(list(c(dS, dI, dR, dP, dSn, dIn, dRn, dAlpha)))
  # ...
}

然后仅在事件函数中更改。

现在我们可以先在没有事件的情况下运行模型，然后在有事件的情况下运行模型。这里我们得到另一个错误，因为事件函数应该只返回状态向量，所以它应该是：

event_func <- function(t, states, parameters) {
  with(as.list(c(states, parameters)), {
    alpha <- 0
    c(S, I, R, P, Sn, In, Rn, alpha)
  })
}

times <- 0:100
out <- ode(
         y = states_1_take,
         times = times, 
         func = model_1_take,
         parms = parameters_1_take,
         method = "lsoda", 
         events = list(func=event_func, root = TRUE), 
         rootfun = root_func)

plot(out)

最后两个注意事项：

• 可以在此处使用默认的

  lsoda

  解算器，它会在刚性和非刚性解算器之间自动切换，还具有求根和事件功能。仅当需要指定稀疏雅可比矩阵的结构时才需要

  lsodes

  求解器。
• 状态变量有不同的数量级。这有时会导致稳定性和准确性的数值问题。在这里，明智的做法是将状态重新调整到 0.001 到 1000 之间的“公共范围”，或者为状态分配单独的绝对容差：

  atol = c(.....)