使用GEKKO将双种群微分模型拟合至测量数据的技术问询

模型定义 #

需要拟合的双种群微分方程及输出如下：

dS/dt = (a - b) * S
dR/dt = (a - b - c)* R
X(t) = S(t) + R(t)

目标是通过(time, measurement)格式的测量数据，估计参数 a、b、c，让模型输出X(t)与测量数据的误差最小。

完整实现代码 #

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

# 替换为你的真实测量数据
time_data = np.array([0, 1, 2, 3, 4, 5])
meas_data = np.array([2.0, 3.1, 4.5, 6.2, 8.1, 10.5])

# 初始化GEKKO模型
m = GEKKO(remote=False)
m.time = time_data

# 定义待估计参数，设置初始值与合理范围
a = m.FV(value=0.5, lb=0, ub=2)
b = m.FV(value=0.1, lb=0, ub=1)
c = m.FV(value=0.05, lb=0, ub=0.5)
# 标记参数参与优化
a.STATUS = 1
b.STATUS = 1
c.STATUS = 1

# 定义状态变量，初始值可根据实际调整
S = m.Var(value=meas_data[0]/2)
R = m.Var(value=meas_data[0]/2)

# 写入微分方程
m.Equation(S.dt() == (a - b) * S)
m.Equation(R.dt() == (a - b - c) * R)

# 定义模型输出X(t)
X = m.Intermediate(S + R)

# 定义目标函数：最小化拟合值与测量值的平方误差和
meas = m.Param(value=meas_data)
m.Obj((X - meas)**2)

# 配置求解器并运行
m.options.IMODE = 5  # 动态数据拟合模式
m.options.NODES = 3   # 每个时间点的节点数，提升拟合精度
m.solve(disp=True)

# 输出估计参数
print("估计参数值：")
print(f"a = {a.value[0]:.4f}")
print(f"b = {b.value[0]:.4f}")
print(f"c = {c.value[0]:.4f}")

# 绘制拟合结果对比图
plt.figure(figsize=(8,5))
plt.plot(time_data, meas_data, 'bo', label='测量数据')
plt.plot(m.time, X.value, 'r-', label='拟合曲线')
plt.xlabel('时间')
plt.ylabel('X(t) = S(t) + R(t)')
plt.legend()
plt.grid(True)
plt.show()

关键细节说明 #

• 参数约束：通过lb/ub给参数设置合理范围，能避免求解过程中出现无意义的参数值，提升收敛稳定性。
• 状态变量初始值：如果对初始种群分布不确定，也可以将S、R的初始值设为待估计变量（把m.Var()改为m.FV()并设置STATUS=1）。
• 求解模式：IMODE=5是GEKKO专门针对动态模型拟合时间序列数据的模式，会自动处理微分方程的数值积分与参数优化。

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