如何通过拟合实验数据来估计一组复杂微分方程的参数？

Question

我试图通过在 python 中使用 gekko 来估计以下一组微分方程的参数，但出现了一些错误。这些是我正在使用的方程式：

使用的代码是：

from gekko import GEKKO

#experimental data
t_data=[0,2,4,6,8,10,12,14,16,18,20]
x_data=[0.0844,0.2068,0.8046,1.8019,2.4655,2.7467,2.7765,2.6966,2.6529,2.6605,2.5464]    
L_data=[0,18.194,36.389,540.069,958.987,1326.418,1069.154,1195.935,1256.966,1422.267,1267.442]
c_data=[9.845,9.4193,9.0340,7.6427,5.9962,5.2468,4.1849,4.4343,4.2462,3.8870,3.6511]
s_data=[5.0619,4.3798,2.6438,0.6220,0.557,0.492,0.4268,0.415,0.4017,0.395,0.3906]

m = GEKKO(remote=False)
m.time = t_data
x = m.CV(value=x_data); x.FSTATUS = 1  # fit to measurements
L = m.CV(value=L_data); L.FSTATUS = 1
c = m.CV(value=c_data); c.FSTATUS = 1
s = m.CV(value=s_data); s.FSTATUS = 1

umax = m.FV(); umax.STATUS = 1         # adjustable parameters
kc = m.FV(); kc.STATUS = 1
smin = m.FV(); smin.STATUS = 1              
ks = m.FV(); ks.STATUS = 1
kd = m.FV(); kd.STATUS = 1
a = m.FV(); a.STATUS = 1              
b = m.FV(); b.STATUS = 1
yxc = m.FV(); yxc.STATUS = 1
q = m.FV(); q.STATUS = 1              
yxs = m.FV(); yxs.STATUS = 1

# differential equations
m.Equation(x.dt() == umax*(c/(kc+c))*((s-smin)/(ks+s-smin))*x-kd*x )           
m.Equation(L.dt() == a*x.dt()+b*x)
m.Equation(c.dt() == (-1/yxc)*x.dt()-q*x*(c/(kc+c)))
m.Equation(s.dt() == (-1/yxs)*umax*(c/(kc+c))*((s-smin)/(ks+s-smin))*x)

#Global Options
m.options.IMODE   = 5
m.options.EV_TYPE = 2 
m.options.NODES   = 3 
m.options.SOLVER  = 3 

m.solve()

执行代码后出现如下错误：

**Error**: Exception: Access Violation
At line 2683 of file MUMPS/src/dmumps_part2.F
Traceback: not available, compile with -ftrace=frame or -ftrace=full

**Error**: 'results.json' not found. Check above for additional error details
Traceback (most recent call last):

  File "C:\Users\vcham\Desktop\2020\belgium\python\parameter estimation bg\test2.py", line 56, in <module>
    m.solve()

  File "C:\Users\vcham\anaconda3\lib\site-packages\gekko\gekko.py", line 2145, in solve
    self.load_JSON()

  File "C:\Users\vcham\anaconda3\lib\site-packages\gekko\gk_post_solve.py", line 13, in load_JSON
    f = open(os.path.join(self._path,'options.json'))

FileNotFoundError: [Errno 2] No such file or directory: 'C:\\Users\\vcham\\AppData\\Local\\Temp\\tmplna43pbhgk_model1\\options.json'

请建议我解决这些错误或估计参数的方法。提前谢谢！！！

Answer 1

带有线性求解器 MUMPS 的 IPOPT 求解器在尝试解决问题时崩溃。获得成功解决方案的一种方法是使用 m.options.SOLVER=1 切换到 APOPT 求解器，并在 FV 参数上设置上限和下限。

from gekko import GEKKO

#experimental data
t_data=[0,2,4,6,8,10,12,14,16,18,20]
x_data=[0.0844,0.2068,0.8046,1.8019,2.4655,2.7467,2.7765,2.6966,2.6529,2.6605,2.5464]    
L_data=[0,18.194,36.389,540.069,958.987,1326.418,1069.154,1195.935,1256.966,1422.267,1267.442]
c_data=[9.845,9.4193,9.0340,7.6427,5.9962,5.2468,4.1849,4.4343,4.2462,3.8870,3.6511]
s_data=[5.0619,4.3798,2.6438,0.6220,0.557,0.492,0.4268,0.415,0.4017,0.395,0.3906]

m = GEKKO(remote=False)
m.time = t_data
x = m.CV(value=x_data); x.FSTATUS = 1  # fit to measurements
L = m.CV(value=L_data); L.FSTATUS = 1
c = m.CV(value=c_data); c.FSTATUS = 1
s = m.CV(value=s_data); s.FSTATUS = 1

umax = m.FV(1,lb=0.1,ub=10); umax.STATUS = 1         # adjustable parameters
kc = m.FV(1,lb=0.1,ub=10); kc.STATUS = 1
smin = m.FV(1,lb=0.1,ub=10); smin.STATUS = 1              
ks = m.FV(1,lb=0.1,ub=10); ks.STATUS = 1
kd = m.FV(1,lb=0.1,ub=10); kd.STATUS = 1
a = m.FV(1,lb=0.1,ub=10); a.STATUS = 1              
b = m.FV(1,lb=0.1,ub=10); b.STATUS = 1
yxc = m.FV(1,lb=0.1,ub=10); yxc.STATUS = 1
q = m.FV(1,lb=0.1,ub=10); q.STATUS = 1              
yxs = m.FV(1,lb=0.1,ub=10); yxs.STATUS = 1

# differential equations
m.Equation(x.dt() == umax*(c/(kc+c))*((s-smin)/(ks+s-smin))*x-kd*x )           
m.Equation(L.dt() == a*x.dt()+b*x)
m.Equation(c.dt() == (-1/yxc)*x.dt()-q*x*(c/(kc+c)))
m.Equation(s.dt() == (-1/yxs)*umax*(c/(kc+c))*((s-smin)/(ks+s-smin))*x)

#Global Options
m.options.IMODE   = 5
m.options.EV_TYPE = 2 
m.options.NODES   = 3 
m.options.SOLVER  = 1 

m.solve()

print('umax: ' + str(umax.value[0]))
print('kc: ' + str(kc.value[0]))
print('smin: ' + str(smin.value[0]))
print('ks: ' + str(ks.value[0]))
print('kd: ' + str(kd.value[0]))
print('a: ' + str(a.value[0]))
print('b: ' + str(b.value[0]))
print('yxc: ' + str(yxc.value[0]))
print('q: ' + str(q.value[0]))
print('yxs: ' + str(yxs.value[0]))

import matplotlib.pyplot as plt
plt.subplot(4,1,1)
plt.plot(t_data,x.value)
plt.plot(t_data,x_data)
plt.subplot(4,1,2)
plt.plot(t_data,L.value)
plt.plot(t_data,L_data)
plt.subplot(4,1,3)
plt.plot(t_data,c.value)
plt.plot(t_data,c_data)
plt.subplot(4,1,4)
plt.plot(t_data,s.value)
plt.plot(t_data,s_data)
plt.show()

我对所有参数应用了0.1 的下限和10.0 的上限，但这些可能需要修改，因为一些解决方案在边界处完成。我建议您扩大 (*) 的界限，直到最佳解决方案没有在界限内结束。报告的解决方案是：

umax: 10.0*
kc: 0.1*
smin: 0.11236933654
ks: 6.2350087285
kd: 0.20089406528
a: 10.0*
b: 10.0*
yxc: 1.6785099655
q: 0.1*
yxs: 6.7395055839