R语言nls拟合光合作用非矩形双曲线报错，如何确定最优参数起始值

当前报错的核心诱因 #

你遇到的singular gradient、NaNs produced、步长小于最小因子等报错，本质是原代码的初值设置脱离了非矩形双曲线模型的物理参数范围，尤其是Theta = 1的设置已经触发了模型退化：

• 非矩形双曲线的凸度参数Theta物理范围为0 < Theta < 1，取Theta=1时，模型表达式会完全退化为直角双曲线的分段形式，求导时会出现梯度为0的奇异点，直接触发nls的梯度报错。
• 原设置的AQY=0.05与你当前数据集的低光段实测斜率偏差较大，优化过程中会频繁出现根号内计算项为负的情况，生成NaN值，导致优化步长不断缩减直到低于默认的minFactor阈值。

可直接落地的初值确定方法 #

所有方法都适配你当前的数据集，不需要额外安装拓展包：

1. 基于模型物理意义手动锚定初值 #

直接从观测数据拆解每个参数的合理初值，避免拍脑袋设数：

• Am：为光饱和阶段的总同化速率，取最高光强（1500 μmol·m⁻²·s⁻¹）下的平均净光合+Rd初值，你这个数据集可直接设为21，远好于原来max(Cultivar)-min(Cultivar)的计算方式。
• AQY：为低光段（PPFD<200）净光合随光强变化的线性斜率，取观测值计算为(平均200 PPFD下的净光合 - 平均0 PPFD下的净光合)/200，你这个数据集算出来约为0.02，可直接设为0.022。
• Rd：为PPFD=0时的净光合绝对值，取0光强下的三次重复平均，你这个数据集约为0.7，设为0.7即可。
• Theta：基于作物观测的常规范围取0.8，完全避开0或1的参数边界，不会出现梯度奇异。

2. 网格搜索自动找最优初值 #

如果手动锚定的初值还是跑不通，直接用全网格搜索找残差最小的参数组合当start，完全规避猜初值的问题，直接运行以下代码即可：

# 划定参数的合理搜索范围
param_grid <- expand.grid(
  Am = seq(18, 25, length.out = 12),
  AQY = seq(0.01, 0.05, length.out = 12),
  Rd = seq(0.3, 1.2, length.out = 12),
  Theta = seq(0.2, 0.9, length.out = 12)
)
# 计算每组参数的残差平方和
sse <- apply(param_grid, 1, function(p) {
  pred <- (1/(2*p["Theta"])) * (p["AQY"]*NLRC$PPFD + p["Am"] - 
           sqrt((p["AQY"]*NLRC$PPFD + p["Am"])^2 - 4*p["AQY"]*p["Theta"]*p["Am"]*NLRC$PPFD)) - p["Rd"]
  sum((NLRC$Cultivar - pred)^2, na.rm = TRUE)
})
# 选残差最小的组作为nls的起始值
best_start <- as.list(param_grid[which.min(sse), ])

用这个方法选的初值和最优解的距离足够近，nls迭代时几乎不会触发步长缩减、梯度奇异的报错。

3. 调整nls的内置容错参数 #

就算初值有小的偏差，也可以通过修改nls的控制参数降低报错概率：

• 把默认的minFactor从1e-3调低到1e-6，给优化器留出足够的步长调整空间；
• 把最大迭代次数从默认的50提到1000，避免迭代提前终止。
  修改后的可运行代码片段：

reg_nrh <- nls(
  Cultivar ~ (1/(2*Theta))*(AQY*PPFD+Am-sqrt((AQY*PPFD+Am)^2-4*AQY*Theta*Am*PPFD))-Rd,
  data = NLRC,
  start = best_start,
  control = nls.control(minFactor = 1e-6, maxiter = 1000, warnOnly = TRUE)
)

拟合后的绘图代码修正 #

原代码里参数提取时存在顺序依赖问题、y轴范围设置过宽的问题，修正后的参数提取与绘图部分可参考：

# 按参数名提取拟合值，避免顺序变动导致的赋值错误
Amnrh <- coef(reg_nrh)["Am"]
AQYnrh <- coef(reg_nrh)["AQY"]
Rdnrh <- coef(reg_nrh)["Rd"]
Thetanrh <- coef(reg_nrh)["Theta"]

# 生成拟合曲线
ppfd = seq(from = 0, to = 1500)
pnnrh <- (1/(2*Thetanrh))*(AQYnrh*ppfd+Amnrh-sqrt((AQYnrh*ppfd+Amnrh)^2-4*AQYnrh*Thetanrh*Amnrh*ppfd)) - Rdnrh

# 绘图
plot(NLRC$PPFD, NLRC$Cultivar, main = "Cultivar", xlab="", ylab="", 
     ylim=c(-2, 25), cex.lab=1.2, cex.axis=1.5, cex=2)
mtext(expression("PPFD ("*mu*"mol photons "*m^-2*s^-1*")"),side=1,line=3.3,cex=1.5)
mtext(expression(P[net]*" ("*mu*"mol "*CO[2]*" "*m^-2*s^-1*")"),side=2,line=2.5,cex=1.5)
lines(ppfd, pnnrh, col="Green", lwd=2)

注：拟合过程中如果出现sqrt内为负的警告，不需要额外处理，只要最终收敛的参数值在合理物理范围内，就说明结果是可靠的。

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