警告：lmdif：info = 0。nls.lm() 函数的输入参数不正确

Question

我有以下解决方案。

$$ rac{{dx}}{{dt}} = k_+(2d11 + d12 - 2x - z) - k_{-}x$$

$$ rac{{dy}}{{dt}} = k_+(2d22 + d12 - 2y - z) - k_{-}y$$

$$ rac{{dz}}{{dt}} = 2k_+(2d11 + d12 - 2x - z)(2d22 + d12 - 2y - z) - k_{-}z$$

我想用 R 编写一个程序来求解模型方程，然后将文章中表格中的数据拟合到方程的解中。我想根据 Levenberg-Marquardt 最小二乘误差算法找到最适合的 k+ 和 k−。

以下代码是我尝试这样做的。

library(deSolve)
library(minpack.lm)

model = function(t, state, parms) {
  with(as.list(c(state, parms)), {
    dxdt = k_plus*(2*d11 + d12 - 2*x - z) - k_minus*x
    dydt = k_plus*(2*d22 + d12 - 2*y - z)-k_minus * y
    dzdt = 2*k_plus*(2*d11 + d12 - 2*x - z)*(2*d22 + d12 - 2*y - z) - k_minus*z
    return(list(c(dxdt, dydt, dzdt)))
  })
}

# Defining an objective function
objective = function(par){
  k_plus = par[1]
  k_minus = par[2]
  x0 = 1.71
  y0 = 1.75
  z0 = 0.62
  d11 = 1.71
  d12 = 1.75
  d22 = 0.62
  
# Vector of times
times_vec = c(1/30, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 16, 18, 20)

initial_state = c(x = x0, y = y0, z = z0)
  
  # Solving the ODE
  out = lsode(y = initial_state, 
              times = times_vec,
              func = model,
              parms = c(k_plus = k_plus, k_minus = k_minus, d11 = d11, d12 = d12, d22 = d22),
              maxsteps = 100000)
  
  error = sum((out[,"x"] - xdata)^2 + (out[,"y"] - ydata)^2 + (out[,"z"] - zdata)^2)
  return(error)
}



# The table from the article
xdata = c(1.71, 1.56, 1.40, 1.28, 1.14, 1.17, 1.14, 1.06, 1.10, 1.10, 1.07, 1.09, 1.06, 1.03, 1.06, 1.06, 1.05)
ydata = c(1.75, 1.59, 1.43, 1.31, 1.17, 1.20, 1.17, 1.09, 1.13, 1.13, 1.10, 1.13, 1.09, 1.06, 1.09, 1.09, 1.08)
zdata = c(0.62, 0.93, 1.26, 1.50, 1.78, 1.71, 1.78, 1.93, 1.86, 1.86, 1.91, 1.87, 1.93, 1.99, 1.94, 1.94, 1.95)


# 
fit = nls.lm(par = c(k_plus = 200,
                     k_minus = 2),
             fn = objective)

当我运行代码时，除了使用 nls.lm 函数之外，一切正常。我收到以下错误

警告：lmdif：info = 0。输入参数不正确。

我尝试过不同的参数。我查找了该函数以确保我的参数具有正确的格式，但我似乎找不到问题所在。

Answer 1

为了调试问题，我建议在尝试拟合模型之前比较模型和数据。我们可以清楚地看到，提供数值的模型结果与数据相差甚远。

library(deSolve)
library(minpack.lm)

model = function(t, state, parms) {
  with(as.list(c(state, parms)), {
    dxdt = k_plus*(2*d11 + d12 - 2*x - z) - k_minus*x
    dydt = k_plus*(2*d22 + d12 - 2*y - z) - k_minus*y
    dzdt = 2*k_plus*(2*d11 + d12 - 2*x - z)*(2*d22 + d12 - 2*y - z) - k_minus*z
    return(list(c(dxdt, dydt, dzdt)))
  })
}

times_vec = c(1/30, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 16, 18, 20)
k_plus = 200; k_minus = 2
x0 = 1.71; y0 = 1.75; z0 = 0.62
d11 = 1.71; d12 = 1.75; d22 = 0.62
initial_state = c(x = x0, y = y0, z = z0)

out = lsode(y = initial_state, 
            times = times_vec,
            func = model,
            parms = c(k_plus = k_plus, k_minus = k_minus, d11 = d11, d12 = d12, d22 = d22),
            maxsteps = 100000)


# The table from the article
data = data.frame(
  time = times_vec,
  x = c(1.71, 1.56, 1.40, 1.28, 1.14, 1.17, 1.14, 1.06, 1.10, 1.10, 1.07, 1.09, 1.06, 1.03, 1.06, 1.06, 1.05),
  y = c(1.75, 1.59, 1.43, 1.31, 1.17, 1.20, 1.17, 1.09, 1.13, 1.13, 1.10, 1.13, 1.09, 1.06, 1.09, 1.09, 1.08),
  z = c(0.62, 0.93, 1.26, 1.50, 1.78, 1.71, 1.78, 1.93, 1.86, 1.86, 1.91, 1.87, 1.93, 1.99, 1.94, 1.94, 1.95)
)

plot(out, obs=data, mfrow=c(1, 3))

因此，我建议检查：

• 模型方程和
• 数值，特别是初始状态显然与 d..-参数相同。

此外，我建议使用

lsoda

求解器而不是

lsode

。人们还可以考虑另一种优化包，例如FME，请参阅 https://doi.org/10.18637/jss.v033.i03 该包是多种优化算法（包括

nls.lm

）的包装器，但对拟合 ODE 模型有更多支持。