含4个ODE的4参数动态参数估计求解求助

我参考动态估计相关示例，针对含4个微分方程的模型做4参数动态参数估计，仅其中1个方程有实验数据。调整代码后求解，所有参数返回0，求适配4个微分方程场景的解决办法。

from gekko import GEKKO
import matplotlib.pyplot as plt  # plot solution
import numpy as np

t_data = [1.22, 5.15, 23.67, 51.17, 74.58, 97.83, 118.97, 143.33, 166.73, 192.08, 222.83, 245.08, 266.67, 286.18, 309.90]
m4_data = [3634.8, 7035.9, 7797.8, 9351.4, 10041.0, 10674.6, 11339.5, 11115.7, 11225.1, 11465.4, 11383.2, 11456.5, 11506.8, 11683.6, 11588.2]


m = GEKKO()
m.time = t_data
m4 = m.CV(value=m4_data); m4.FSTATUS = 1  # fit to measurement
m1,m2,m3 = m.Array(m.Var,3,value=3)
a,b,c,d = m.Array(m.FV,4)
t = 50 
v = 7 
s = 10
#DiffEqs
#m1' = -aS*(b-m3/v)
#m2' = -cs*(b-m3/v)
#m3' = (-d*t/v)*m3 -m1'-m2'
#m4' = d*t*m3/v

m.Equations([m1.dt() == -a*s*(b-m3/v),m2.dt() == -c*s*(b-m3/v),m3.dt() ==(-d*t/v)*m3-m1.dt()-m2.dt() , m4.dt() == d*t*m3/v])

m.options.IMODE = 5
m.options.NODES = 5  
m.solve()   
print(a.value[0],b.value[0],c.value[0],d.value[0])

问题原因与解决步骤 #

参数全为0通常是因为自由变量（FV）未启用估计、初始值设置不合理，或者模型方程存在依赖关系导致参数不可识别。按以下步骤修改：

1. 启用参数估计功能
   默认情况下Gekko的FV变量不会参与优化，必须设置FSTATUS=0（允许调整），并给参数设置合理初始值（不能全为0）：

a,b,c,d = m.Array(m.FV,4, value=1)  # 设置初始值为1（根据实际物理意义调整）
for param in [a,b,c,d]:
    param.FSTATUS = 0  # 启用参数估计

2. 修正方程中的变量冲突
   代码里定义了t=50，但Gekko的时间变量是m.time，这里的t会被当成常数，而不是随时间变化的自变量。如果模型中的t是时间变量，应该用Gekko内置的m.t：

# 方程中把t换成m.t
m.Equations([m1.dt() == -a*s*(b-m3/v),
             m2.dt() == -c*s*(b-m3/v),
             m3.dt() == (-d*m.t/v)*m3 - m1.dt() - m2.dt(),
             m4.dt() == d*m.t*m3/v])

3. 设置CV变量的拟合权重（可选）
   可以给m4设置MEAS_GAP或WMEAS来调整拟合精度要求，比如：

m4.MEAS_GAP = 0.1  # 允许测量值与模型预测值有10%的偏差范围

4. 添加参数上下界（可选）
   如果参数有物理意义（比如正数），设置上下界帮助优化器收敛：

a.LOWER = 0; a.UPPER = 100
b.LOWER = 0; b.UPPER = 1000
c.LOWER = 0; c.UPPER = 100
d.LOWER = 0; d.UPPER = 10

5. 调整求解器选项（可选）
   若收敛困难，可尝试调整节点数或改用其他求解器：

m.options.NODES = 3  # 减少节点数降低计算量
m.options.SOLVER = 3  # 使用IPOPT求解器（默认是APOPT）

修改后的完整代码 #

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

t_data = [1.22, 5.15, 23.67, 51.17, 74.58, 97.83, 118.97, 143.33, 166.73, 192.08, 222.83, 245.08, 266.67, 286.18, 309.90]
m4_data = [3634.8, 7035.9, 7797.8, 9351.4, 10041.0, 10674.6, 11339.5, 11115.7, 11225.1, 11465.4, 11383.2, 11456.5, 11506.8, 11683.6, 11588.2]

m = GEKKO()
m.time = t_data

# 测量变量
m4 = m.CV(value=m4_data)
m4.FSTATUS = 1  # 使用测量数据
m4.MEAS_GAP = 0.1  # 拟合偏差范围

# 状态变量
m1, m2, m3 = m.Array(m.Var, 3, value=3)

# 待估计参数
a, b, c, d = m.Array(m.FV, 4, value=1)
for param in [a, b, c, d]:
    param.FSTATUS = 0  # 启用参数估计
# 设置参数上下界（根据物理意义调整）
a.LOWER = 0; a.UPPER = 100
b.LOWER = 0; b.UPPER = 1000
c.LOWER = 0; c.UPPER = 100
d.LOWER = 0; d.UPPER = 10

v = 7 
s = 10

# 微分方程（使用m.t作为时间变量）
m.Equations([
    m1.dt() == -a*s*(b - m3/v),
    m2.dt() == -c*s*(b - m3/v),
    m3.dt() == (-d*m.t/v)*m3 - m1.dt() - m2.dt(),
    m4.dt() == d*m.t*m3/v
])

m.options.IMODE = 5  # 动态参数估计模式
m.options.NODES = 3
m.options.SOLVER = 3  # 使用IPOPT求解器
m.solve(disp=True)  # 显示求解过程

print(f"参数估计结果：a={a.value[0]:.4f}, b={b.value[0]:.4f}, c={c.value[0]:.4f}, d={d.value[0]:.4f}")

# 绘制拟合曲线
plt.plot(m.time, m4_data, 'bo', label='实验数据')
plt.plot(m.time, m4.value, 'r-', label='模型预测')
plt.xlabel('时间')
plt.ylabel('m4')
plt.legend()
plt.show()

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