使用 Gekko 和 Python 来拟合数据的数值 ODE 解
Use Gekko and Python to fit a numerical ODE solution to data
使用 Gekko 对数据拟合数值 ODE 解。
大家好!
我想知道是否可以使用 GEKKO 来拟合 ODE 的系数。
我尝试复制 example given here.
失败
这是我想出的(但有缺陷——我或许应该提一下,不幸的是我的数学技能很差):
import numpy as np
from gekko import GEKKO
tspan = [0, 0.1, 0.2, 0.4, 0.8, 1]
Ca_data = [2.0081, 1.5512, 1.1903, 0.7160, 0.2562, 0.1495]
m = GEKKO(remote=False)
t = m.Param(value=tspan)
m.time = t
Ca_m = m.Param(value=Ca_data)
Ca = m.Var()
k = m.FV(value=1.3)
k.STATUS = 1
m.Equation( Ca.dt() == -k * Ca)
m.Obj( ((Ca-Ca_m)**2)/Ca_m )
m.options.IMODE = 2
m.solve(disp=True)
print(k.value[0]) #2.58893455 is the solution
有人可以帮我吗?
非常感谢,
马丁
(这是我第一次 post 在这里——如果我做了不恰当的事情,请保持温和。)
您的解决方案很接近,但您需要:
- 更多节点(默认值=2)以提高准确性。 Gekko 只添加您定义的点。参见 additional information on collocation。
- 将
Ca
定义为 m.CV()
以使用内置错误模型,而不是节点 >=3 的 m.Var()
和 m.Obj
。否则,每个配置区间的内部节点也与测量值匹配,这给出了一个稍微错误的答案。
- 设置
EV_TYPE=2
使用平方误差。绝对值 objective EV_TYPE=1
(默认值)给出了正确但略有不同的答案。
import numpy as np
from gekko import GEKKO
m = GEKKO(remote=False)
m.time = [0, 0.1, 0.2, 0.4, 0.8, 1]
Ca_data = [2.0081, 1.5512, 1.1903, 0.7160, 0.2562, 0.1495]
Ca = m.CV(value=Ca_data); Ca.FSTATUS = 1 # fit to measurement
k = m.FV(value=1.3); k.STATUS = 1 # adjustable parameter
m.Equation(Ca.dt()== -k * Ca) # differential equation
m.options.IMODE = 5 # dynamic estimation
m.options.NODES = 5 # collocation nodes
m.options.EV_TYPE = 2 # squared error
m.solve(disp=True) # display solver output
print(k.value[0]) # 2.58893455 is the curve_fit solution
解决方案是k=2.5889717102
。绘图显示与测量值的匹配。
import matplotlib.pyplot as plt # plot solution
plt.plot(m.time,Ca_data,'ro')
plt.plot(m.time,Ca.value,'bx')
plt.show()
关于微分和代数方程模型的参数估计有additional tutorials and course material。
使用 Gekko 对数据拟合数值 ODE 解。
大家好! 我想知道是否可以使用 GEKKO 来拟合 ODE 的系数。 我尝试复制 example given here.
失败这是我想出的(但有缺陷——我或许应该提一下,不幸的是我的数学技能很差):
import numpy as np
from gekko import GEKKO
tspan = [0, 0.1, 0.2, 0.4, 0.8, 1]
Ca_data = [2.0081, 1.5512, 1.1903, 0.7160, 0.2562, 0.1495]
m = GEKKO(remote=False)
t = m.Param(value=tspan)
m.time = t
Ca_m = m.Param(value=Ca_data)
Ca = m.Var()
k = m.FV(value=1.3)
k.STATUS = 1
m.Equation( Ca.dt() == -k * Ca)
m.Obj( ((Ca-Ca_m)**2)/Ca_m )
m.options.IMODE = 2
m.solve(disp=True)
print(k.value[0]) #2.58893455 is the solution
有人可以帮我吗? 非常感谢, 马丁
(这是我第一次 post 在这里——如果我做了不恰当的事情,请保持温和。)
您的解决方案很接近,但您需要:
- 更多节点(默认值=2)以提高准确性。 Gekko 只添加您定义的点。参见 additional information on collocation。
- 将
Ca
定义为m.CV()
以使用内置错误模型,而不是节点 >=3 的m.Var()
和m.Obj
。否则,每个配置区间的内部节点也与测量值匹配,这给出了一个稍微错误的答案。 - 设置
EV_TYPE=2
使用平方误差。绝对值 objectiveEV_TYPE=1
(默认值)给出了正确但略有不同的答案。
import numpy as np
from gekko import GEKKO
m = GEKKO(remote=False)
m.time = [0, 0.1, 0.2, 0.4, 0.8, 1]
Ca_data = [2.0081, 1.5512, 1.1903, 0.7160, 0.2562, 0.1495]
Ca = m.CV(value=Ca_data); Ca.FSTATUS = 1 # fit to measurement
k = m.FV(value=1.3); k.STATUS = 1 # adjustable parameter
m.Equation(Ca.dt()== -k * Ca) # differential equation
m.options.IMODE = 5 # dynamic estimation
m.options.NODES = 5 # collocation nodes
m.options.EV_TYPE = 2 # squared error
m.solve(disp=True) # display solver output
print(k.value[0]) # 2.58893455 is the curve_fit solution
解决方案是k=2.5889717102
。绘图显示与测量值的匹配。
import matplotlib.pyplot as plt # plot solution
plt.plot(m.time,Ca_data,'ro')
plt.plot(m.time,Ca.value,'bx')
plt.show()
关于微分和代数方程模型的参数估计有additional tutorials and course material。