python - GEKKO 和 Scipy.optimize 导致非线性参数估计的不同结果
问题描述
我正在学习如何使用 GEKKO 来解决参数估计问题，作为第一步，我正在开发我之前使用 Scipy 最小化例程实现的示例问题。这些是根据 APMonitor.com 中提供的信息和其中提供的课程完成的。当前的问题是从以下获得的甲醇制烃过程的间歇反应器模拟： http ://www.daetools.com/docs/tutorials-all.html#tutorial-che-opt-5
模型描述可以在下面进一步描述的代码中遵循，但考虑的基本步骤是：
A --> B
A + B --> C
C + B --> P
A --> C
A --> P
A + B --> P
其中实验数据可用于 A、C 和 P 的浓度随时间的变化。该模型的目标是估计六个基本反应 (k1-k6) 的速率常数。我现在遇到的困难是我的 GEKKO 模型和基于 Scipy.optimize 的模型导致不同的参数估计，尽管使用相同的实验数据和参数的初始猜测。我还将这个模型与使用 gPROMS 和 Athena Visual Studio 开发的模型进行了比较，scipy 模型与这些闭源程序获得的参数估计一致。每个程序的估计参数如下所示：
Scipy 模型（L-BFGS-B 优化器）：[k1 k2 k3 k4 k5 k6] = [2.779, 0., 0.197, 3.042, 2.148, 0.541]
GEKKO 模型（IPOPT 优化器）：[k1 k2 k3 k4 k5 k6] = [3.7766387559, 1.1826920269e-07, 0.21242442412, 4.130394645, 2.4232122905, 3.3140978171]
有趣的是，两种模型在优化结束时都导致相同的目标函数值 0.0123，并且在物种浓度与时间的图中看起来相似。我尝试更改 GEKKO 的优化器并将公差收紧到 1E-8 无济于事。我的猜测是我的 GEKKO 模型没有正确设置，但我找不到它的问题。如果能向我指出可能导致模型差异的可能问题，我们将不胜感激。我附上以下两个脚本：
Scipy模型
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import minimize
import matplotlib.pyplot as plt
#Experimental data
times = np.array([0.0, 0.071875, 0.143750, 0.215625, 0.287500, 0.359375, 0.431250,
0.503125, 0.575000, 0.646875, 0.718750, 0.790625, 0.862500,
0.934375, 1.006250, 1.078125, 1.150000])
A_obs = np.array([1.0, 0.552208, 0.300598, 0.196879, 0.101175, 0.065684, 0.045096,
0.028880, 0.018433, 0.011509, 0.006215, 0.004278, 0.002698,
0.001944, 0.001116, 0.000732, 0.000426])
C_obs = np.array([0.0, 0.187768, 0.262406, 0.350412, 0.325110, 0.367181, 0.348264,
0.325085, 0.355673, 0.361805, 0.363117, 0.327266, 0.330211,
0.385798, 0.358132, 0.380497, 0.383051])
P_obs = np.array([0.0, 0.117684, 0.175074, 0.236679, 0.234442, 0.270303, 0.272637,
0.274075, 0.278981, 0.297151, 0.297797, 0.298722, 0.326645,
0.303198, 0.277822, 0.284194, 0.301471])
def rxn(x, k): #rate equations in power law form r = k [A][B]
A = x[0]
B = x[1]
C = x[2]
P = x[3]
k1 = k[0]
k2 = k[1]
k3 = k[2]
k4 = k[3]
k5 = k[4]
k6 = k[5]
r1 = k1 * A
r2 = k2 * A * B
r3 = k3 * C * B
r4 = k4 * A
r5 = k5 * A
r6 = k6 * A * B
return [r1, r2, r3, r4, r5, r6] #returns reaction rate of each equation
#mass balance diff eqs, function calls rxn function
def mass_balances(t, x, *args):
k = args
r = rxn(x, k)
dAdt = - r[0] - r[1] - r[3] - r[4] - r[5]
dBdt = + r[0] - r[1] - r[2] - r[5]
dCdt = + r[1] - r[2] + r[3]
dPdt = + r[2] + r[4] + r[5]
return [dAdt, dBdt, dCdt, dPdt]
IC = [1.0, 0, 0, 0] #Initial conditions of species A, B, C, P
ki= [1, 1, 1, 1, 1, 1]
#Objective function definition
def obj_fun(k):
#solve initial value problem over time span of data
sol = solve_ivp(mass_balances,[min(times),max(times)],IC, args = (k), t_eval=(times))
y_model = np.vstack((sol.y[0],sol.y[2],sol.y[3])).T
obs = np.vstack((A_obs, C_obs, P_obs)).T
err = np.sum((y_model-obs)**2)
return err
bnds = ((0, None), (0, None),(0, None),(0, None),(0, None),(0, None))
model = minimize(obj_fun,ki, bounds=bnds, method = 'L-BFGS-B')
k_opt = model.x
print(k_opt.round(decimals = 3))
y_calc = solve_ivp(mass_balances,[min(times),max(times)],IC, args = (model.x), t_eval=(times))
plt.plot(y_calc.t, y_calc.y.T)
plt.plot(times,A_obs,'bo')
plt.plot(times,C_obs,'gx')
plt.plot(times,P_obs,'rs')
GEKKO模型
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
#Experimental data
times = np.array([0.0, 0.071875, 0.143750, 0.215625, 0.287500, 0.359375, 0.431250,
0.503125, 0.575000, 0.646875, 0.718750, 0.790625, 0.862500,
0.934375, 1.006250, 1.078125, 1.150000])
A_obs = np.array([1.0, 0.552208, 0.300598, 0.196879, 0.101175, 0.065684, 0.045096,
0.028880, 0.018433, 0.011509, 0.006215, 0.004278, 0.002698,
0.001944, 0.001116, 0.000732, 0.000426])
C_obs = np.array([0.0, 0.187768, 0.262406, 0.350412, 0.325110, 0.367181, 0.348264,
0.325085, 0.355673, 0.361805, 0.363117, 0.327266, 0.330211,
0.385798, 0.358132, 0.380497, 0.383051])
P_obs = np.array([0.0, 0.117684, 0.175074, 0.236679, 0.234442, 0.270303, 0.272637,
0.274075, 0.278981, 0.297151, 0.297797, 0.298722, 0.326645,
0.303198, 0.277822, 0.284194, 0.301471])
m = GEKKO(remote = False)
t = m.time = times
Am = m.CV(value=A_obs, lb = 0)
Cm = m.CV(value=C_obs, lb = 0)
Pm = m.CV(value=P_obs, lb = 0)
A = m.Var(1, lb = 0)
B = m.Var(0, lb = 0)
C = m.Var(0, lb = 0)
P = m.Var(0, lb = 0)
Am.FSTATUS = 1
Cm.FSTATUS = 1
Pm.FSTATUS = 1
k1 = m.FV(1, lb = 0)
k2 = m.FV(1, lb = 0)
k3 = m.FV(1, lb = 0)
k4 = m.FV(1, lb = 0)
k5 = m.FV(1, lb = 0)
k6 = m.FV(1, lb = 0)
k1.STATUS = 1
k2.STATUS = 1
k3.STATUS = 1
k4.STATUS = 1
k5.STATUS = 1
k6.STATUS = 1
r1 = m.Var(0, lb = 0)
r2 = m.Var(0, lb = 0)
r3 = m.Var(0, lb = 0)
r4 = m.Var(0, lb = 0)
r5 = m.Var(0, lb = 0)
r6 = m.Var(0, lb = 0)
m.Equation(r1 == k1 * A)
m.Equation(r2 == k2 * A * B)
m.Equation(r3 == k3 * C * B)
m.Equation(r4 == k4 * A)
m.Equation(r5 == k5 * A)
m.Equation(r6 == k6 * A * B)
#mass balance diff eqs, function calls rxn function
m.Equation(A.dt() == - r1 - r2 - r4 - r5 - r6)
m.Equation(B.dt() == r1 - r2 - r3 - r6)
m.Equation(C.dt() == r2 - r3 + r4)
m.Equation(P.dt() == r3 + r5 + r6)
m.Obj((A-Am)**2+(P-Pm)**2+(C-Cm)**2)
m.options.IMODE = 5
m.options.SOLVER = 3 #IPOPT optimizer
m.options.RTOL = 1E-8
m.options.OTOL = 1E-8
m.solve()
k_opt = [k1.value[0],k2.value[0], k3.value[0], k4.value[0], k5.value[0], k6.value[0]]
print(k_opt)
plt.plot(t,A)
plt.plot(t,C)
plt.plot(t,P)
plt.plot(t,B)
plt.plot(times,A_obs,'bo')
plt.plot(times,C_obs,'gx')
plt.plot(times,P_obs,'rs')
解决方案
这里有几个建议：
- 设置
m.options.NODES=3或更高设置为 6 以获得更好的积分精度。 - 将
Am,Cm,设置Pm为参数而不是变量。它们是固定输入。 - 尝试不同的初始条件。可能有多个局部最小值。
- 目标函数可以是平坦的，以便不同的参数值给出相同的目标函数值。您可以测试参数置信区间以查看数据是否给出了窄或宽的联合置信区域。
以下是修改后的结果：
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
#Experimental data
times = np.array([0.0, 0.071875, 0.143750, 0.215625, 0.287500, 0.359375, 0.431250,
0.503125, 0.575000, 0.646875, 0.718750, 0.790625, 0.862500,
0.934375, 1.006250, 1.078125, 1.150000])
A_obs = np.array([1.0, 0.552208, 0.300598, 0.196879, 0.101175, 0.065684, 0.045096,
0.028880, 0.018433, 0.011509, 0.006215, 0.004278, 0.002698,
0.001944, 0.001116, 0.000732, 0.000426])
C_obs = np.array([0.0, 0.187768, 0.262406, 0.350412, 0.325110, 0.367181, 0.348264,
0.325085, 0.355673, 0.361805, 0.363117, 0.327266, 0.330211,
0.385798, 0.358132, 0.380497, 0.383051])
P_obs = np.array([0.0, 0.117684, 0.175074, 0.236679, 0.234442, 0.270303, 0.272637,
0.274075, 0.278981, 0.297151, 0.297797, 0.298722, 0.326645,
0.303198, 0.277822, 0.284194, 0.301471])
m = GEKKO(remote=False)
t = m.time = times
Am = m.Param(value=A_obs, lb = 0)
Cm = m.Param(value=C_obs, lb = 0)
Pm = m.Param(value=P_obs, lb = 0)
A = m.Var(1, lb = 0)
B = m.Var(0, lb = 0)
C = m.Var(0, lb = 0)
P = m.Var(0, lb = 0)
k = m.Array(m.FV,6,value=1,lb=0)
for ki in k:
ki.STATUS = 1
k1,k2,k3,k4,k5,k6 = k
r1 = m.Var(0, lb = 0)
r2 = m.Var(0, lb = 0)
r3 = m.Var(0, lb = 0)
r4 = m.Var(0, lb = 0)
r5 = m.Var(0, lb = 0)
r6 = m.Var(0, lb = 0)
m.Equation(r1 == k1 * A)
m.Equation(r2 == k2 * A * B)
m.Equation(r3 == k3 * C * B)
m.Equation(r4 == k4 * A)
m.Equation(r5 == k5 * A)
m.Equation(r6 == k6 * A * B)
#mass balance diff eqs, function calls rxn function
m.Equation(A.dt() == - r1 - r2 - r4 - r5 - r6)
m.Equation(B.dt() == r1 - r2 - r3 - r6)
m.Equation(C.dt() == r2 - r3 + r4)
m.Equation(P.dt() == r3 + r5 + r6)
m.Minimize((A-Am)**2)
m.Minimize((P-Pm)**2)
m.Minimize((C-Cm)**2)
m.options.IMODE = 5
m.options.SOLVER = 3 #IPOPT optimizer
m.options.RTOL = 1E-8
m.options.OTOL = 1E-8
m.options.NODES = 5
m.solve()
k_opt = []
for ki in k:
k_opt.append(ki.value[0])
print(k_opt)
plt.plot(t,A)
plt.plot(t,C)
plt.plot(t,P)
plt.plot(t,B)
plt.plot(times,A_obs,'bo')
plt.plot(times,C_obs,'gx')
plt.plot(times,P_obs,'rs')
plt.show()
推荐阅读
- saml - Okta 作为 SP（入站 SAML）如何在成功登录后配置重定向 URL
- applescript - mac 上的 Spotify 流甲板插件缺少重复切换
- android - Unity - Android - BLE 蓝牙权限
- c# - 在 C# 中结合动态绑定和非托管约束
- javascript - 必须表达包含两个逻辑运算符（|| 和 &&）的复杂条件的 switch 语句看起来如何？
- python - 我想通过麦克风发送我正在计算机上收听的音频
- mysql - 如何将“当前”mysql时间戳值迁移到postgreSQL
- linux - GNU 并行 - 通过多个 ssh 跳转主机运行命令
- c - 当应用于无符号操作数和有符号操作数时，整数除法会产生不同的结果吗？
- c# - C# 使用 count 属性反序列化对象的最佳方法