带时间依赖传播率的SEIARWF传染病动力学模型参数估计优化时lsoda非法输入错误的解决咨询

阿华AIGC实验室

2026-4-28

让我来帮你一步步排查并解决这个lsoda报错的问题——这个错误通常是因为微分方程求解时遇到了非法输入（比如负值状态变量、数值溢出或者状态变量维度不匹配），我们来逐个修复：

1. 移除状态变量中的`N`，修正维度不匹配问题

你把N作为初始状态变量传入了ode函数，但在微分方程里N = S+E+I+A+R是派生变量，不需要单独作为状态变量。lsoda要求每个状态变量都有对应的导数，而你之前返回的导数列表里包含了N的导数，这会导致求解器混淆状态变量的数量，直接触发非法输入错误。

修改点：

初始状态移除N：

initials <- c(S = S0, E = E0, I = I0, A =A0, R = R0, W = W0, FD = FD0)

微分方程函数的返回值去掉N的导数：

list(c(dS, dE, dI, dA, dR, dW, dFD))

同时记得在函数开头先计算N，避免后续引用出错：

SEIARWF <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    N <- S + E + I + A + R # 先计算N
    # 后续微分方程代码不变
  })
}

2. 修正`dW`和`dFD`的系数，避免数值溢出

你代码里dW和dFD的系数是3014000000000和1507000000000，这会导致W和FD的数值瞬间膨胀到天文数字（比如初始I=1429，乘以3e12直接得到4e15），超出R的数值处理范围，引发lsoda的非法输入错误。这大概率是单位换算错误，比如把“每小时排放”写成了“每月排放”，或者系数多写了几个零。

修改点：
根据文献中的实际单位调整系数，比如先缩小到合理量级测试：

dW=3014*I+1507*A-1.3333*W
dFD=3014*I+1507*A-1.3333*FD

之后再根据文献校准到正确数值。

3. 添加参数约束，确保传播率始终为正

传播率betaH=aH*(1+bH cos(2πt/12))必须大于0，否则会出现负的感染率，导致E等状态变量的导数为负，进而变成负值触发报错。虽然你给bH设了-1到1的范围，但当cos(2πt/12)=-1时，1+bH*(-1)可能等于0，导致传播率为0甚至负数。

修改点：
在RSS函数中添加参数检查，一旦传播率出现非正情况，返回一个极大的RSS值，让优化算法自动避开这些参数组合：

RSS <- function(parameters) {
  names(parameters) <- c("aH", "bH", "aW", "bW", "aF", "bF")
  # 确保传播率始终为正
  betaH_min <- parameters["aH"]*(1 - abs(parameters["bH"]))
  betaW_min <- parameters["aW"]*(1 - abs(parameters["bW"]))
  betaF_min <- parameters["aF"]*(1 - abs(parameters["bF"]))
  if(betaH_min <= 1e-8 || betaW_min <= 1e-8 || betaF_min <= 1e-8){
    return(1e18) # 返回极大值，让优化算法放弃该参数组合
  }
  out <- ode(y = initials, times = t, func = SEIARWF, parms = parameters)
  fit <- out[ , 4] # 第4列对应I变量，和Infected数据匹配
  sum((Infected - fit)^2)
}

同时把bH等参数的上下限调整为-0.99到0.99，彻底避免传播率为0的情况：

Opt <- optim(par=c(0.5, 0.3, 0.1, 0.3, 0.2, 0.3), fn = RSS, method = "L-BFGS-B",
             lower = c(1e-6, -0.99, 1e-6, -0.99, 1e-6, -0.99), 
             upper = c(5, 0.99, 5, 0.99, 5, 0.99),
             control = list(maxit=1000), hessian = TRUE)

4. 调整初始参数，降低优化收敛难度

你原来的初始参数aH=0.0005太小，可能导致初始感染率远低于实际数据，让优化算法一开始就进入不合理的参数空间。建议把初始参数调整到更合理的范围，比如aH=0.5、aW=0.1、aF=0.2，帮助优化算法更快找到可行区域。

完整修改后的代码

整合所有修改后的代码如下，你可以直接测试：

library(deSolve)
# SEIARWF微分方程 - 修改后
SEIARWF <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    N <- S + E + I + A + R
    dS <- 0.001048*N + 0.00595*R - 
      (aH*(1+bH*cos(2*pi*t/12)))*S*(I+A)/N -
      (aW*(1+bW*cos(2*pi*t/12)))*S*W/N -
      (aF*(1+bF*cos(2*pi*t/12)))*S*FD/N - 
      0.000582*S
    dE <- (aH*(1+bH*cos(2*pi*t/12)))*S*(I+A)/N +
      (aW*(1+bW*cos(2*pi*t/12)))*S*W/N +
      (aF*(1+bF*cos(2*pi*t/12)))*S*FD/N -
      0.75*E - 0.000582*E
    dI <- 0.1*0.75*E - 1*I - 0.000582*I
    dA <- 0.9*0.75*E - 1*A - 0.000582*A
    dR <- 1*I + 1*A - 0.00595*R - 0.000582*R
    # 修正dW和dFD的系数（请根据文献校准到正确值）
    dW <- 3014*I + 1507*A - 1.3333*W
    dFD <- 3014*I + 1507*A - 1.3333*FD
    
    list(c(dS, dE, dI, dA, dR, dW, dFD))
  })
}
# 初始状态 - 移除N
N0 <- 1292270000
Infected <- c(1429,1813,2227,1752,1541,1467,1249,1188,1022,973,1052,731,1539,1604,
              2190,1838,1327,1183,1038,1052,942,997,1045,786,1704,2219,2768,2116,
              1818,1490,1228,1207,1042,1179,1213,1023,2183,2039,2616,2176,1842,1643,
              1571,1612,1334,1306,1292,963,1843,1871,2726,2112,1786,1425,1305,1226,
              1061,1096,1097,977,1930,2444,2730,2048,1735,1487,1338,1424,1157,1195,
              1237,1550,2288,2304,3076,2567,2180,1719,1584,1461,1421,1468,1684,1930,
              2342,2930,4178,3270,2459,2118,2010,2116,1871,1872,2062,1974,2338,3180,
              3557,2662,2354,1889,1843,1906,1737,1763,1959,2083,3076,3145,3628,2772,
              2241,1794,1874,1899,1723,1764,1960,2026,2655,3018,3140,2510,2052,1836,
              1958,2054,1908,1849,1883,2125,2893,2301,3032,2668,2165,1956,1991,1931,
              1792,1982,2131,2327,2497,2468,3318,2676,2410,2143,2094,2105,1966,1885,
              2065,2295,2306,2752,3141,2568,2559,2397,2390,2439,2049,2006,2191,2216,
              2811,2441,3163,2651,2489,2160,2324,2249,1928,1863,2304,2220)
I0 <- Infected[1]
A0 <- Infected[1]*9
E0 <- 19054
R0 <- 303166542
S0 = N0 - E0 - I0 - A0 - R0
W0 <- 5169000000000000
FD0 <- 11080000000000000
initials <- c(S = S0, E = E0, I = I0, A =A0, R = R0, W = W0, FD = FD0)
# 时间序列
t <- 0:(length(Infected)-1)
# 基于RSS的参数优化 - 添加参数检查
RSS <- function(parameters) {
  names(parameters) <- c("aH", "bH", "aW", "bW", "aF", "bF")
  # 确保传播率始终为正
  betaH_min <- parameters["aH"]*(1 - abs(parameters["bH"]))
  betaW_min <- parameters["aW"]*(1 - abs(parameters["bW"]))
  betaF_min <- parameters["aF"]*(1 - abs(parameters["bF"]))
  if(betaH_min <= 1e-8 || betaW_min <= 1e-8 || betaF_min <= 1e-8){
    return(1e18)
  }
  out <- ode(y = initials, times = t, func = SEIARWF, parms = parameters)
  fit <- out[ , 4] # 第4列是I，对应Infected数据
  sum((Infected - fit)^2)
}
# 调整初始参数和上下限
Opt <- optim(par=c(0.5, 0.3, 0.1, 0.3, 0.2, 0.3), fn = RSS, method = "L-BFGS-B",
             lower = c(1e-6, -0.99, 1e-6, -0.99, 1e-6, -0.99), 
             upper = c(5, 0.99, 5, 0.99, 5, 0.99),
             control = list(maxit=1000), hessian = TRUE)
Opt
Opt_par <- setNames(Opt$par, c("aH", "bH", "aW", "bW", "aF", "bF"))
Opt_par

额外测试建议

先单独运行ode函数测试，不进行优化：比如用初始参数运行ode(y=initials, times=t, func=SEIARWF, parms=c(aH=0.5,bH=0.3,aW=0.1,bW=0.3,aF=0.2,bF=0.3))，查看输出是否有负值或异常大的数值，逐步排查问题。
如果优化仍不收敛，可以尝试更换优化方法（比如method="Nelder-Mead"），或者增大maxit参数的迭代次数。

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