使用nls.lm拟合PDE模型遇ftol错误，求排查模型/数据问题

我正尝试将实验数据拟合至如下3个PDE方程组：

\frac{dX}{dt} = \mu_{max} \cdot \frac{S}{K_s+S} \cdot (1-\frac{X}{X_m}) \cdot X
\frac{dS}{dt} = -\frac{1}{Y_{X/S}} \cdot \frac{dX}{dt} - m_s \cdot X - k_{LD} \cdot \frac{dL}{dt}
\frac{dL}{dt} = -k_l \cdot L

运行代码时出现两个错误：

Error in chol.default(object$hessian) :
the leading minor of order 6 is not positive definite

以及

reason terminated: Relative error in the sum of squares is at most `ftol'

无法确定问题出在模型还是数据上，相关代码如下：

library(deSolve)
library(ggplot2)
library(minpack.lm)
library(reshape2)
time <- c(0,2,4,5,6,7,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42)
X <- c(0.010044683, 0.011653262, 0.013306524, 0.014155496, 0.014959786, 
       0.015630027, 0.016434316, 0.018132261, 0.019651475, 0.021394102, 
       0.023002681, 0.024432529, 0.026085791, 0.027828418, 0.029347632, 
       0.031090259, 0.032698838, 0.034262735, 0.035915996, 0.037569258, 
       0.039311886, 0.040786416, 0.042350313, 0.04409294)
S <- c(0.585443495, 0.537584046, 0.505306743, 0.476368471, 0.464125356, 
       0.451882241, 0.441308642, 0.426283001, 0.403466287, 0.394005698, 
       0.37508452, 0.371188984, 0.353937322, 0.355050332, 0.33668566, 
       0.343920228, 0.328338082, 0.333346629, 0.319433998, 0.314425451, 
       0.309416904, 0.302182336, 0.294391263, 0.288269706)
L <- c(0.372443284, 0.370216418, 0.361865672, 0.362979104, 0.357968657, 
  0.359638806, 0.355741791, 0.345720896, 0.34961791, 0.334029851, 
  0.329019403, 0.334029851, 0.327349254, 0.322338806, 0.319555224, 
  0.315101493, 0.31398806, 0.305080597, 0.302297015, 0.297286567, 
  0.295059701, 0.289492537, 0.282255224, 0.27668806)
data <- data.frame(time,X,S,L)
attach(data)

Hybrid <- function(t,c,parms) {
  k1 <- parms$k1 # mumax
  k2 <- parms$k2 # Ks
  k3 <- parms$k3 # Xm
  k4 <- parms$k4 # Y_x/s
  k5 <- parms$k5 # m_s
  k6 <- parms$k6 # k_LD
  k7 <- parms$k7 # k_l
  r <- numeric(length(c)) 
  r[1] <-  k1 * (c["S"]  / ( k2 + c["S"] )) * c["X"] * (1-c["X"]/k3)  # r[1] = dX/dt = mumax.(S/(Ks+S)).(1-X/Xm).X
  r[2] <- -1/k4 * r[1] - k5 * c["X"] - k6 * r[3] # r[2] = dS/dt = -1/Y_x/s * dX/dt - m_s.X - k_LD * dL/dt
  r[3] <- -k7*c["L"] # dL/dt = -kl * L
  return(list(r))       
}

residuals <- function(parms){
  cinit <- c( X = data[1,2], S = data[1,3], L = data[1,4] )
  t <- sort(unique(c(seq(0, 42, 1), data$time)))
  k1 <- parms[1]
  k2 <- parms[2]
  k3 <- parms[3]
  k4 <- parms[4]
  k5 <- parms[5]
  k6 <- parms[6]
  k7 <- parms[7]
  out <- ode( y = cinit, 
              times = t, 
              func = Hybrid, 
              parms = list( k1 = k1, k2 = k2, k3 = k3, k4 = k4, k5 = k5, k6 = k6, k7 = k7) )
  out_data <- data.frame(out)
  out_data <- out_data[out_data$time %in% data$time,]
  pred_data <- melt(out_data,id.var="time",variable.name="Substance",value.name="Conc")
  exp_data <- melt(data,id.var="time",variable.name="Substance",value.name="Conc")
  residuals <- pred_data$Conc-exp_data$Conc
  return(residuals)
}

parms <- c( k1=0.8,  # mumax
            k2=1.2,  # Ks 
            k3=0.06, # Xm 
            k4=0.05, # Y_x/s
            k5=0.001,# m_s 
            k6=0.05, # k_LD
            k7=0.003)# kl

fitval <- nls.lm(par=parms,fn=residuals)
summary(fitval)
fitval

parest=as.list(coef(fitval))
cinit <- c( X=data[1,2], S=data[1,3], L=data[1,4] )
t<-seq(0, 42, 1)
parms<-as.list(parest)
out<-ode(y=cinit,times=t,func=Hybrid,parms=parms)
out_data<-data.frame(out)
names(out_data)<-c("time","X_pred","S_pred","L_pred")

tmppred<-melt(out_data,id.var=c("time"),variable.name="Substance",value.name="Concentration")
tmpexp<-melt(data,id.var=c("time"),variable.name="Substance",value.name="Concentration")
p<-ggplot(data=tmppred,aes(x=time,y=Concentration,color=Substance,linetype=Substance))+geom_line()
print(p)
p<-p+geom_line(data=tmpexp,aes(x=time,y=Concentration,color=Substance,linetype=Substance))
p<-p+geom_point(data=tmpexp,aes(x=time,y=Concentration,color=Substance))
p<-p+scale_linetype_manual(values=c(0,1,0,1,0,1))
p<-p+scale_color_manual(values=rep(c("red","blue","purple"),each=2))+theme_bw()
print(p)

错误原因分析 #

1. 海森矩阵非正定错误 #

该错误说明模型参数间存在强共线性，或部分参数对残差的影响极小，导致拟合算法无法区分这些参数的独立效应。从模型结构看，k6（k_LD）和k7（k_l）耦合紧密：dS/dt依赖dL/dt，而dL/dt仅由k7控制，两个参数难以同时被数据约束。

2. ftol终止提示 #

表示拟合算法认为当前参数下的残差平方和相对误差已达ftol阈值，但这大概率是参数陷入局部最优，或参数空间存在平坦区域，算法无法继续优化。

排查与修复步骤 #

步骤1：修正初始参数合理性 #

• k4（Y_{X/S}）初始值0.05偏离实际数据：X从0.01增长到0.04，S从0.58降到0.28，消耗0.3单位S，对应Y约为0.1，建议调整初始值为k4=0.1。
• k2（K_s）初始值1.2远大于S的初始值0.58，导致S/(K_s+S)始终处于低底物限制区，k1与k2难以区分，建议将初始k2调整为0.3（接近S的范围）。

步骤2：分步拟合简化模型 #

• 先固定k6=0，忽略L对S的影响，拟合X和S的子模型，得到k1,k2,k3,k4,k5的合理估计后，再加入L的拟合。
• 单独拟合L的衰减模型：对L(t) = L0 * exp(-k7*t)取对数，用线性回归先得到k7的估计值，再代入整体模型，减少参数自由度。

步骤3：优化残差函数效率 #

残差函数中无需生成额外时间点，直接使用实验时间点即可，减少计算量：

residuals <- function(parms){
  cinit <- c( X = data[1,2], S = data[1,3], L = data[1,4] )
  t <- data$time  # 直接使用实验时间点
  k1 <- parms[1]
  k2 <- parms[2]
  k3 <- parms[3]
  k4 <- parms[4]
  k5 <- parms[5]
  k6 <- parms[6]
  k7 <- parms[7]
  out <- ode( y = cinit, 
              times = t, 
              func = Hybrid, 
              parms = list( k1 = k1, k2 = k2, k3 = k3, k4 = k4, k5 = k5, k6 = k6, k7 = k7) )
  out_data <- data.frame(out)
  pred_data <- melt(out_data,id.var="time",variable.name="Substance",value.name="Conc")
  exp_data <- melt(data,id.var="time",variable.name="Substance",value.name="Conc")
  residuals <- pred_data$Conc - exp_data$Conc
  return(residuals)
}

步骤4：添加参数物理约束 #

使用nls.lm的lower和upper参数给参数设定正数边界（符合物理意义）：

fitval <- nls.lm(par=parms, fn=residuals, 
                 lower=c(0.01, 0.01, 0.04, 0.01, 0, 0, 1e-4),
                 upper=c(2, 1, 0.1, 0.5, 0.01, 1, 0.01))

（边界值可根据实验背景调整，比如Xm需大于X的最大值0.04）

步骤5：检查数据质量 #

L数据存在轻微回升（如t=5比t=4高），可能是实验噪声，可考虑平滑L数据，或在拟合时给L的残差赋予较低权重。

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