ode求解器报错：`hmax'必须为非负值，求R语言水蚤模型调试帮助

问题描述 #

运行下方水蚤模型脚本时，持续收到报错：

Error in checkInput(y, times, func, rtol, atol, jacfunc, tcrit, hmin, : `hmax' must be a non-negative value

已尝试调整时间序列、初始状态、参数及ODE求解器方法，问题仍未解决。

报错根源分析 #

1. 语法错误：代码中使用了中文上标符号ˆ代替R的幂运算符^，导致Reproduction计算时语法异常，破坏了ODE函数的正常返回值。
2. 触发时机逻辑错误：蜕皮和换水的判断条件使用了循环起始时间Time而非结束时间TimeOut，导致状态更新时机延迟，可能引发时间序列异常。
3. 全局变量依赖：Moulting函数依赖全局state变量，可能导致状态读取不一致，引发计算错误。
4. 求解器步长异常：当模型出现非连续状态修改（如蜕皮）或导数异常时，默认求解器可能自动设置负的最大步长hmax。

修正后的完整代码 #

# 模型参数定义
neonateWeight <- 1.1
reproductiveWeight <- 7.5
maximumWeight <- 60.0
instarDuration <- 3.0

ksFood <- 85.0
IngestWeight <- 132.0
maxIngest <- 1.05
assimilEff <- 0.8
maxReproduction <- 0.8
respirationRate <- 0.25

transferTime <- 2
foodInMedium <- 509
numberIndividuals <- 32

# 初始状态
state <- c(INDWEIGHT = neonateWeight,
           EGGWEIGHT = 0,
           FOOD = foodInMedium)

# ODE模型函数：修正幂运算符为英文^
model <- function(t, state, parameters) {
  with(as.list(c(state)), {
    WeightFactor <- (IngestWeight - INDWEIGHT) / (IngestWeight - neonateWeight)
    MaxIngestion <- maxIngest * WeightFactor
    Ingestion <- MaxIngestion * INDWEIGHT * FOOD / (FOOD + ksFood)
    Respiration <- respirationRate * INDWEIGHT
    Growth <- Ingestion * assimilEff - Respiration
    
    # 修正幂运算符为^
    if (Growth <= 0. | INDWEIGHT < reproductiveWeight) {
      Reproduction <- 0.
    } else {
      WeightRatio <- reproductiveWeight / INDWEIGHT
      Reproduction <- maxReproduction * (1 - WeightRatio^2)
    }
    
    dINDWEIGHT <- (1 - Reproduction) * Growth
    dEGGWEIGHT <- Reproduction * Growth
    # 限制FOOD不低于0，避免负食物导致计算异常
    dFOOD <- ifelse(FOOD <= 0, 0, -Ingestion * numberIndividuals)
    
    list(c(dINDWEIGHT, dEGGWEIGHT, dFOOD),
         c(Ingestion = Ingestion,
           Respiration = Respiration,
           Reproduction = Reproduction))
  }) 
}

# 蜕皮函数：改为传入state参数，避免全局变量依赖
Moulting <- function(current_state) {
  with(as.list(current_state), {
    refLoss <- 0.24
    cLoss <- 3.1
    INDLength <- (INDWEIGHT / 3.0)^(1/2.6) # Mariam公式
    WeightLoss <- refLoss * INDLength^cLoss
    # 确保蜕皮后体重为正
    max(INDWEIGHT - WeightLoss, neonateWeight * 0.8)
  })  
}

# 主循环：修正触发时机判断逻辑，显式设置hmax
TimeFrom <- 0
TimeEnd <- 40
TimeMoult <- TimeFrom + instarDuration
TimeTransfer <- TimeFrom + transferTime
Time <- TimeFrom
Outdt <- 0.1
out <- NULL

while (Time < TimeEnd) {
  TimeOut <- min(TimeMoult, TimeTransfer, TimeEnd)
  times <- seq(Time, TimeOut, by = Outdt)
  
  # 显式设置hmax为正数，指定求解器（可选）
  out1 <- as.data.frame(ode(state, times, model, parms = 0, hmax = 1, method = "ode45"))
  out <- rbind(out, out1)
  
  lout <- nrow(out1)
  state <- c(
    INDWEIGHT = out1[lout, "INDWEIGHT"],
    EGGWEIGHT = out1[lout, "EGGWEIGHT"],
    FOOD = out1[lout, "FOOD"]
  )
  
  # 用TimeOut判断是否到达蜕皮/换水时间，修正触发时机
  if (TimeOut >= TimeMoult) {
    state[1] <- Moulting(state)
    state[2] <- 0 # 蜕皮后卵重清零
    TimeMoult <- TimeOut + instarDuration # 更新下一次蜕皮时间
  }
  if (TimeOut >= TimeTransfer) {
    state[3] <- foodInMedium # 换水重置食物量
    TimeTransfer <- TimeOut + transferTime # 更新下一次换水时间
  }
  
  Time <- TimeOut
}

# 绘图代码
par(mfrow = c(2,2)) # 调整子图布局
plot(out$time, out$FOOD, type = "l", main = "Food",
     xlab = "time, days", ylab = "μgC/l")
plot(out$time, out$INDWEIGHT, type = "l", main = "Individual weight",
     xlab = "time, days", ylab = "μgC")
plot(out$time, out$EGGWEIGHT, type = "l", main = "Egg weight",
     xlab = "time, days", ylab = "μgC")
plot(out$time, out$Ingestion, type = "l", main = "Ingestion",
     xlab = "time, days", ylab = "μgC/ind/day")
mtext(outer = TRUE, side = 3, "DAPHNIA model", cex = 1.5)
par(mfrow = c(1,1)) # 恢复默认布局

关键修改说明 #

• 语法修正：将WeightRatioˆ2替换为WeightRatio^2，解决幂运算语法错误。
• 状态更新时机：将判断条件从Time >= TimeMoult改为TimeOut >= TimeMoult，确保到达指定时间时立即触发蜕皮和换水。
• 全局变量消除：Moulting函数改为接收current_state参数，避免依赖全局状态导致的计算不一致。
• 求解器参数设置：在ode函数中显式设置hmax=1，强制求解器使用非负最大步长；可选指定method="ode45"提升稳定性。
• 边界保护：添加FOOD下限判断（不低于0），蜕皮后体重限制为不低于幼体体重的80%，避免极端值导致模型崩溃。

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