如何将插值法引入常微分方程组参数拟合优化？

核心思路 #

将原ODE拟合中x0的微分方程替换为带可调参数的插值函数，把插值的控制点（如分段插值的节点取值）作为优化变量，与原ODE参数（k0、k1、k2、p1、p2）一起纳入lmfit的优化框架。通过求解x1、x2的ODE（基于插值得到的x0(t)），结合x0的实验数据残差、x1/x2的模拟与实验残差构建目标函数，完成整体优化。

实现步骤与代码示例 #

1. 依赖库导入与实验数据准备 #

先导入所需工具库，同时模拟实验数据（实际使用时替换为真实数据）：

import numpy as np
from scipy.integrate import solve_ivp
from scipy.interpolate import interp1d
import lmfit
import matplotlib.pyplot as plt

# 生成模拟实验数据（替换为真实数据）
t_exp = np.linspace(0, 10, 20)
# 真实ODE用于生成带噪声的实验数据
def true_ode(t, y, k0, k1, k2, p1, p2):
    x0, x1, x2 = y
    dx0dt = -k0*x0 + p1
    dx1dt = k0*x0 - k1*x1 + p2
    dx2dt = k1*x1 - k2*x2
    return [dx0dt, dx1dt, dx2dt]

true_params = {'k0':0.5, 'k1':0.3, 'k2':0.2, 'p1':0.1, 'p2':0.05}
sol_true = solve_ivp(true_ode, [t_exp[0], t_exp[-1]], [1.0, 0.0, 0.0], 
                     args=(true_params['k0'], true_params['k1'], true_params['k2'], 
                           true_params['p1'], true_params['p2']), t_eval=t_exp)
# 加噪声模拟实验数据
x0_exp = sol_true.y[0] + np.random.normal(0, 0.05, len(t_exp))
x1_exp = sol_true.y[1] + np.random.normal(0, 0.05, len(t_exp))
x2_exp = sol_true.y[2] + np.random.normal(0, 0.05, len(t_exp))

2. 定义插值节点与优化参数 #

选择插值节点（建议从实验时间点中选取子集，平衡自由度与过拟合风险），并构建lmfit参数组：

# 选取插值节点（每隔4个实验点取1个）
t_nodes = t_exp[::4]
n_nodes = len(t_nodes)

# 初始化lmfit参数
params = lmfit.Parameters()
# 添加原ODE参数，设置合理初始值与边界
params.add('k0', value=0.4, min=0)
params.add('k1', value=0.2, min=0)
params.add('k2', value=0.1, min=0)
params.add('p1', value=0.08, min=0)
params.add('p2', value=0.04, min=0)
# 添加x0插值节点的取值作为优化参数
for i in range(n_nodes):
    params.add(f'x0_node_{i}', value=x0_exp[::4][i], min=0)

3. 自定义目标函数 #

在目标函数中完成插值构建、ODE求解与残差计算：

def objective(params, t_exp, x0_exp, x1_exp, x2_exp, t_nodes):
    # 提取所有参数
    k0 = params['k0'].value
    k1 = params['k1'].value
    k2 = params['k2'].value
    p1 = params['p1'].value
    p2 = params['p2'].value
    x0_nodes = np.array([params[f'x0_node_{i}'].value for i in range(n_nodes)])
    
    # 构建x0的分段线性插值函数（可改为kind='cubic'使用三次样条）
    x0_interp = interp1d(t_nodes, x0_nodes, kind='linear', fill_value='extrapolate')
    
    # 定义仅包含x1、x2的ODE（依赖插值得到的x0(t)）
    def ode(t, y):
        x1, x2 = y
        dx1dt = k0 * x0_interp(t) - k1*x1 + p2
        dx2dt = k1*x1 - k2*x2
        return [dx1dt, dx2dt]
    
    # 求解x1、x2的ODE（初始条件根据实际场景调整）
    sol = solve_ivp(ode, [t_exp[0], t_exp[-1]], [0.0, 0.0], t_eval=t_exp)
    x1_sim = sol.y[0]
    x2_sim = sol.y[1]
    
    # 计算残差：x0插值与实验值的差 + x1、x2模拟与实验值的差
    res_x0 = x0_interp(t_exp) - x0_exp
    res_x1 = x1_sim - x1_exp
    res_x2 = x2_sim - x2_exp
    
    # 返回一维残差数组供lmfit优化
    return np.concatenate([res_x0, res_x1, res_x2])

4. 执行拟合与结果可视化 #

运行优化并输出、可视化结果：

# 执行拟合
result = lmfit.minimize(objective, params, args=(t_exp, x0_exp, x1_exp, x2_exp, t_nodes))

# 输出拟合报告
lmfit.report_fit(result)

# 绘制拟合结果
t_plot = np.linspace(0, 10, 100)
# 获取拟合后的x0插值函数
x0_nodes_fit = np.array([result.params[f'x0_node_{i}'].value for i in range(n_nodes)])
x0_fit_interp = interp1d(t_nodes, x0_nodes_fit, kind='linear', fill_value='extrapolate')
x0_fit = x0_fit_interp(t_plot)

# 获取拟合后的x1、x2曲线
k0_fit = result.params['k0'].value
k1_fit = result.params['k1'].value
k2_fit = result.params['k2'].value
p2_fit = result.params['p2'].value

def ode_fit(t, y):
    x1, x2 = y
    dx1dt = k0_fit * x0_fit_interp(t) - k1_fit*x1 + p2_fit
    dx2dt = k1_fit*x1 - k2_fit*x2
    return [dx1dt, dx2dt]

sol_fit = solve_ivp(ode_fit, [t_plot[0], t_plot[-1]], [0.0, 0.0], t_eval=t_plot)
x1_fit = sol_fit.y[0]
x2_fit = sol_fit.y[1]

# 绘图
plt.figure(figsize=(12,8))
plt.subplot(311)
plt.scatter(t_exp, x0_exp, label='x0 实验数据')
plt.plot(t_plot, x0_fit, label='x0 拟合插值曲线', color='red')
plt.scatter(t_nodes, x0_nodes_fit, color='red', marker='s', s=100, label='x0 插值节点')
plt.legend()
plt.title('x0: 实验数据 vs 拟合插值')

plt.subplot(312)
plt.scatter(t_exp, x1_exp, label='x1 实验数据')
plt.plot(t_plot, x1_fit, label='x1 拟合曲线', color='green')
plt.legend()
plt.title('x1: 实验数据 vs 拟合结果')

plt.subplot(313)
plt.scatter(t_exp, x2_exp, label='x2 实验数据')
plt.plot(t_plot, x2_fit, label='x2 拟合曲线', color='blue')
plt.legend()
plt.title('x2: 实验数据 vs 拟合结果')
plt.tight_layout()
plt.show()

关键注意事项 #

• 插值方式选择：示例用分段线性插值，若需更平滑曲线可改用三次样条（kind='cubic'）；若尝试UnivariateSpline，可将样条的节点系数作为优化参数，但需额外处理节点单调性约束。
• 节点数量控制：插值节点过多易过拟合，过少则自由度不足，建议从少量节点开始调试，逐步调整。
• 残差权重调整：若x0实验数据误差较大，可降低其残差权重（如res_x0 * 0.5），反之则提高权重，以平衡各组分的拟合精度。
• 初始值设置：x0节点的初始值建议用对应实验点数据，ODE参数初始值可复用原ODE拟合的结果，提升优化收敛效率。

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