使用FME::modFit/optim为参数子集添加拟合约束遇错求助

尝试通过最小二乘法拟合微分方程模型，待优化的4个参数中有3个（a、b、ks）不能小于0，rs无下限约束。使用FME::modFit和stats::optim时，无论设置部分下限还是全参数宽松上下限，都出现non-finite finite-difference value错误。

以下是数据、模型及辅助函数：

library(deSolve)
library(FME)

# Data
y.data <- data.frame(P = c(57537.4049311277, 58610.5565091595, 59380.7326528789, 59831.1412677501, 59956.9401326381, 59770.9438648549, 59339.7636243282, 58743.8966682831, 58062.6867570216, 57372.6802888231, 56729.6957034719, 56118.5748350006, 55507.8890166606, 54867.5694355373, 54168.9851576162, 53385.1758149806, 52500.082193378, 51534.2686505716, 50516.106686327), 
SSL = c(5918.49121715316, 6223.93819671156, 6522.34271426757, 6817.47760344158, 7115.49425740418, 7423.50644959141, 7748.62985898112, 8097.51802623853, 8472.32658437055, 8872.9327092582, 9294.62369710465, 9730.87922007949, 10175.6848086032, 10623.7661461289, 11075.7535333951, 11538.2076285642, 12018.4500914707, 12518.681172246, 13039.7328357312),
S = c(61.8032786885246, 69.8469945355191, 71.1475409836066, 75.0163934426229, 65.6393442622951, 63.7049180327869, 60.0437158469945, 67.9344262295082, 66.5573770491803, 67.0327868852459, 57.6010928961749, 43.7158469945355, 57.224043715847, 65.1366120218579, 49.2131147540984, 30.7158469945355, 34.3715846994535, 12.3661202185792, 27.4972677595628))

# Initial state conditions
y0 <- c(P=y.data[1,]$P, SSL=y.data[1,]$SSL, S=y.data[1,]$S)

# Initial free parameter values
a <- 0.00075
b <- 0.004464286
rs <- 0.01
ks <- 100
p0 <- c(a=a,b=b,rs=rs,ks=ks)

# Series along which to integrate
t <- 0:18
l <- length(t)

# Differential equation model
LVfn_dc <- function(time, y, param) {
  P = y[1]
  SSL = y[2]
  S = y[3]
  with(as.list(param), {
    dPdt <- ((-0.4101*4)*time^3) + ((18.058*3)*time^2) - ((297.9*2)*time) + 1511.2
    dSSLdt <- ((7.2468*2)*time) + 265.67      
    dSdt <- (rs*S)*((ks-S-(a*P)-(b*SSL))/ks)
    return(list(c(dPdt,dSSLdt,dSdt)))
  })
}

# Function to solve differential equation
sde <- function(param, time, f, y) {
  sol <- deSolve::ode(y = y, 
                      times = time, 
                      func = f, 
                      parms = param,
                      method = "rk4",
                      rtol = 1e-15, atol = 1e-15, maxsteps = 1e9)
  return(sol)
}

# Least squares function for optimization using FME::modFit
LS_modFit <- function(param, time, falala, y, y.values){
  # call sde function to compute y.theta(t) after initializing
  y.theta.all <- c()
  y.theta <- c()
  call.sde <- sde(param,time,falala,y)
  y.theta.all <- call.sde[,-1] # disregard t column
  # only keep y.theta(t.i) values for each data point t.i
  for(i in 1:length(time)){
    if(time[i]%%1 == 0){
      y.theta <- rbind(y.theta,y.theta.all[i,])
    }
  }
  error <- sum((y.values - y.theta)^2)
  return(error)
}

尝试的优化代码：

# Optimizing with limited lower bounds
fit <- FME::modFit(f = LS_modFit, p = p0,
                     time=t, falala=LVfn_dc, y=y0, y.values=y.data,
                     upper=c(Inf,Inf,Inf,Inf), lower=c(0,0,-Inf,1),
                     method = "L-BFGS-B")

设置全参数上下限后仍报错：

# Optimizing with upper and lower bounds for all parameters
fit <- FME::modFit(f = LS_modFit, p = p0,
                     time=t, falala=LVfn_dc, y=y0, y.values=y.data,
                     upper=c(0.02,0.02,1,400), lower=c(0,0,-2,65),
                     method = "L-BFGS-B")

1. 拟合目标与模型结构不匹配 #

你的微分方程中，dPdt和dSSLdt完全由time决定，和当前状态变量P、SSL无关。这意味着P和SSL的模型轨迹是固定的多项式曲线，参数变化无法影响这两个变量的拟合误差——但你却用这两个变量的观测值参与误差计算，导致优化算法无法通过调整参数来降低整体误差，有限差分计算时容易出现非有限值。

解决办法：
如果P和SSL是已知的时间函数（无需拟合），仅针对S计算拟合误差，修改LS_modFit的误差计算部分：

# 仅计算S的误差
error <- sum((y.values$S - y.theta[,3])^2)

2. S的微分方程存在数值不稳定风险 #

当参数变化时，dSdt的表达式(rs*S)*((ks-S-(a*P)-(b*SSL))/ks)可能出现以下问题：

• 若ks接近S + a*P + b*SSL，括号内的值会趋近于0，导致S的变化率异常；若ks小于该值，dSdt为负，S会持续衰减至负数，后续积分会产生数值问题；
• 过严的积分精度（rtol=1e-15、atol=1e-15）容易引发数值溢出。

解决办法：

• 给ks设置合理下限，确保ks始终大于S + a*P + b*SSL的最大值，避免dSdt为负：

  # 基于初始参数计算ks的最小合理值
  min_ks <- max(y.data$S + p0["a"]*y.data$P + p0["b"]*y.data$SSL)
  

• 在微分方程中添加约束，防止S变为负数：

  dSdt <- ifelse(S <= 0, 0, (rs*S)*((ks-S-(a*P)-(b*SSL))/ks))
  

• 放松积分精度，设置rtol=1e-8、atol=1e-8即可满足需求。

3. 冗余的时间点判断 #

你的time序列是0:18，和观测数据的时间点完全对应，无需通过time[i]%%1 ==0判断，直接提取call.sde[,-1]即可，冗余的循环判断可能引发不必要的数值问题。

简化后的LS_modFit函数：

LS_modFit <- function(param, time, falala, y, y.values){
  call.sde <- deSolve::ode(y = y, 
                           times = time, 
                           func = falala, 
                           parms = param,
                           method = "rk4",
                           rtol = 1e-8, atol = 1e-8)
  y.theta <- call.sde[,-1]
  # 仅计算S的拟合误差
  error <- sum((y.values$S - y.theta[,3])^2)
  return(error)
}

4. 调整参数约束与初始值 #

• ks初始值100过小，基于初始参数计算的S + a*P + b*SSL最大值约为75 + 37.5 + 26.76 = 139.26，建议将ks初始值设为150，下限设为140；
• 若rs允许为负，会导致S无限制下降，可根据实际场景考虑是否给rs添加非负约束。

# 计算ks的最小合理值
min_ks <- max(y.data$S + p0["a"]*y.data$P + p0["b"]*y.data$SSL)

# 重新设置初始参数
p0 <- c(a=0.00075, b=0.00446, rs=0.01, ks=150)

# 执行优化
fit <- FME::modFit(f = LS_modFit, p = p0,
                   time=t, falala=LVfn_dc, y=y0, y.values=y.data,
                   upper=c(Inf,Inf,Inf,Inf), lower=c(0,0,-Inf, min_ks),
                   method = "L-BFGS-B")

内容的提问来源于stack exchange，提问作者e-stred