常微分方程组参数优化：改造Lotka-Volterra模型拟合异常排查

核心错误：残差函数逻辑问题 #

你当前的残差函数返回了平方后的误差，这是拟合失效的最主要原因：
lmfit的minimize接口（包括leastsq等最小二乘实现）要求残差函数返回未平方的误差序列，算法会自动计算平方和作为损失函数。你提前对误差做了平方操作，相当于把损失函数从标准的最小二乘sum((y_pred-y_obs)^2)变成了sum((y_pred-y_obs)^4)，会大幅提升大误差点的权重，且让优化器极易陷入局部极小值，完全无法得到合理拟合结果。

修正后的残差函数 #

def residual(paras, t, data):
    arg0 = paras['T0'].value, paras['L0'].value
    model = sol(t, arg0, paras)
    x2_model = model[:, 0]
    # 去掉平方操作，直接返回误差序列
    return (x2_model - data).ravel()

其他影响拟合效果的问题 #

• L0不应固定为常数：你没有L(t)的观测数据，L的初始值对系统演化的影响极大，固定L0=1相当于人为限制了系统的可行动力学路径，大概率和真实数据的演化规律不匹配。建议将L0设为可变参数，设置合理的取值边界，示例如下：

params.add('L0', value=L0, min=1e-3, max=1e5, vary=True)

• 时间单位核对：代码中b = 60*24的赋值逻辑是将速率单位从分钟转换为天，你需要确认t_measured的单位是否为天，如果时间单位是小时，b的取值会偏差24倍，直接导致动力学速率和数据完全不匹配。
• 优化方法选择：修正残差后如果仍然存在局部极小问题，可将优化方法换成全局优化算法differential_evolution，对参数少的ODE拟合问题适配性更好：

result = minimize(residual, params, args=(t_measured, T_measured), method='differential_evolution')

• 数值稳定性验证：你设置的G=1.7e9远大于初始T和L的取值，初始阶段分母T+L+G近似等于G，TL项会非常小，T的初始增长近似为纯指数增长dT/dt ≈ aT，你可以先手动验证a的初始值对应的指数增长速率是否和数据的增长趋势匹配，避免参数边界设置的隐含问题。

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