非线性ODE系统用leastsq拟合实验数据时ODE求解报错如何解决

问题根源 #

你遇到的lsoda步长过小、积分步数超限错误，本质是优化迭代过程中生成了不符合物理意义的参数组合，导致ODE系统刚性极强甚至出现奇点，求解器被迫不断缩小步长直到数值精度溢出，单纯切换求解器或调整mxstep无法从根源解决问题。

具体解决步骤 #

1. 修复代码逻辑错误 #

• 删掉model3函数中V=x(5);后面多余的反引号，避免语法解析异常
• 补全model3的参数传递：F和Sf在Differences中定义但未传入model3，会导致求解时变量未定义返回异常值，修改两处代码：

  // model3函数定义修改为
  function dx=model3(t,x,F,Sf,Kap,Ksa,Ko,Ks,Kia,Kis,p_Amax,q_Amax,qm,q_Smax,Yas,Yoa,Yxa,Yem,Yos,Yxsof,H)
  

  // ode调用处的list修改为
  x_calc=ode('stiff',x0',t0,t,ODEOPTIONS,list(model3,F,Sf,Kap,Ksa,Ko,Ks,Kia,Kis,p_Amax,q_Amax,qm,q_Smax,Yas,Yoa,Yxa,Yem,Yos,Yxsof,H));
  

2. 优化参数约束范围 #

• 所有半饱和常数（Kap、Ksa、Ko、Ks、Kia、Kis）的下界从0改为1e-6，避免参数取0导致动力学项消失、系统进入刚性区域
• 给速率、产率参数加合理上下界，例如q_Smax上界设为2，所有产率系数上界不超过对应反应的理论最大值，避免优化器搜索到无物理意义的参数区间

3. 优化残差函数逻辑 #

• 给ODE求解加错误捕获，求解失败时返回大残差让优化器自动抛弃当前参数组合：

  function f=Differences(k)
      global MYDATA
      // ... 原有参数赋值代码 ...
      try
          x_calc=ode('stiff',x0',t0,t,ODEOPTIONS,list(model3,F,Sf,Kap,Ksa,Ko,Ks,Kia,Kis,p_Amax,q_Amax,qm,q_Smax,Yas,Yoa,Yxa,Yem,Yos,Yxsof,H));
          diffmat=x_calc'-x_exp;
          f=diffmat(:);
      catch
          // 求解失败返回极大残差
          f = 1e6 * ones(length(MYDATA.t)*5,1);
      end
      MYDATA.funeval=MYDATA.funeval+1;
  endfunction
  

• 给不同状态量的残差加权重，抵消量级差异：例如S和A的量级远小于V，可以给V的残差乘以0.01，避免拟合优先匹配大量级变量导致参数偏离合理范围

4. 调整ODE求解配置 #

• 优先选用刚性求解器stiff，手动设置求解精度和最大步数：

  ODEOPTIONS = odeset('mxstep',1e6,'rtol',1e-4,'atol',1e-6);
  

  优化中间步不需要过高的求解精度，适当调松rtol可以大幅提升求解速度，减少步长过小问题
• 初始阶段可以先固定部分确定性高的参数（例如产率系数），只优化动力学参数，收敛后再放开所有参数拟合，降低参数空间维度，避免优化器随机搜索到不合理参数

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