如何用GEKKO拟合两类种群至测量数据并优化初始点

模型定义 #

种群动态变化的微分方程模型如下：

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

优化需求 #

现有代码中种群初始值为固定值，需调整为引入取值范围在0到1之间的变量f，用其表示两类种群的初始比例：

S0 = x_meas[0] * f
R0 = x_meas[0] * (1-f)

修改后的完整代码 #

from gekko import GEKKO
m = GEKKO(remote=False)

# 时间与测量数据
m.time = [0,0.1,0.2,0.3,0.4,0.5,0.8,0.85,0.9,0.95,1]
x_meas = [1.1,1.2,1.56,1.77,1.89,2.1,2.3,2.5,2.7,2.2,2.1]

# 待优化参数a、b、c，取值范围0.5到10
p = m.Array(m.FV,3,lb=0.5,ub=10); a,b,c = p
# 开启参数优化状态
for fv in p:
   fv.STATUS=1

# 初始比例变量f，取值范围0到1，初始值设为0.01
f0 = m.FV(value=0.01, lb=0, ub=1)
f0.STATUS=1

# 定义状态变量S、R、X，下限为0
S,R,X = m.Array(m.Var,3,lb=0)

# 根据f设置初始条件
S0 = x_meas[0] * f0
R0 = x_meas[0] * (1 - f0)
# 绑定初始值到状态变量（确保求解器可优化f0）
m.fix(S, pos=0, val=S0)
m.fix(R, pos=0, val=R0)
X.value = x_meas[0]

# 测量值参数
Xm = m.Param(x_meas)

# 建立微分方程与总和方程
m.Equations([S.dt()==(a-b)*S,
             R.dt()==(a-b-c)*R,
             X==S+R])

# 目标：最小化拟合值与测量值的平方差
m.Minimize((X-Xm)**2)

# 设置求解器模式与配置
m.options.IMODE = 5   # 动态估计模式
m.options.NODES = 3   # 配置配点节点数
m.solve(disp=False)   # 求解，关闭输出

# 打印优化后的参数
print(f'a: {a.value[0]:.4f}, b: {b.value[0]:.4f}, c: {c.value[0]:.4f}, f: {f0.value[0]:.4f}')

# 绘图展示结果
import matplotlib.pyplot as plt
plt.figure(figsize=(7,3))
plt.plot(m.time,S.value,'b:',linewidth=3,label='S种群')
plt.plot(m.time,R.value,'r--',linewidth=2,label='R种群')
plt.plot(m.time,X.value,'k--',label='拟合总和X')
plt.plot(m.time,Xm.value,'kx',label='测量值Xm')
plt.legend(); plt.grid(); plt.xlabel('时间')
plt.tight_layout()
plt.savefig('fit.png',dpi=300)
plt.show()

关键修改说明 #

• 移除原代码中固定的初始值S0,R0 = [0.35,0.65]，替换为基于f0和测量初始值x_meas[0]计算的动态初始条件
• 使用m.fix()方法绑定状态变量S、R的初始值，确保求解器可通过优化f0调整初始比例（直接赋值S.value=S0会因S0为变量而非数值报错）
• 新增打印优化后的f值，便于查看最优初始比例

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