程序员经验分享-hwhale

GEKKO 和 Scipy.optimize 导致非线性参数估计结果不同

2023-12-05

我正在学习如何使用 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')

壁虎模型

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')

这里有一些建议：

  • Set m.options.NODES=3或更高至 6，以获得更好的积分精度。
  • Set Am, Cm, Pm作为参数而不是变量。它们是固定输入。
  • 尝试不同的初始条件。可能存在多个局部最小值。
  • 目标函数可以是平坦的，使得不同的参数值给出相同的目标函数值。你可以测试参数置信区间查看数据是否给出窄或宽的联合置信区域。

以下是修改后的结果：

parameter estimation

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()
本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)

python

gekko

GEKKO 和 Scipy.optimize 导致非线性参数估计结果不同 的相关文章

随机推荐

  • 如何使用ivy发布原生库？

    对于 Java 库项目 要发布的工件非常简单 因为输出是单个 jar 文件 但是 我还有一个包含以下内容的项目要发布 MyLib jar armeabi libStuff so armeabi v7a libStuff so mips li

  • 在ajax调用中访问函数外部的变量时出现问题

    getJSON http 192 168 1 9 8983 solr db select wt json start 0 rows 100 q query json wrf function result each result respo

  • 在java中保存对话框的数据

    if e getActionCommand equals save to file System out println save is pressed StringBuffer fileContent new StringBuffer T

  • 为什么导出/导入默认 ES 模块属性比命名模块属性更快？

    我正在阅读 Material UI 文档 它指出 请注意 在上面的示例中 我们使用了 import RaisedButton from material ui RaisedButton 代替import RaisedButton from

  • jQuery $.post() + IE8

  • 如何卸载RVM？ [复制]

    这个问题在这里已经有答案了 可能的重复 如何从我的系统中删除rvm ruby版本管理器 如何在 Ubuntu 9 10 上卸载 或重新安装 RVM 我搞乱了当前的安装 这很容易 只需执行以下操作 rvm implode or rm rf r

  • Pandas Dataframe，当列不相等时连接两个dt

    假设我有两个形状相同的数据表 即 N 行和 2 列 它们具有相同的列名称 一 二 将第一个表称为 左 然后将第二个表称为 右 如何返回新的数据表FROM表 左 当 一 列中两个表的值分别为不等于 EX Table Left One Two

  • HTML 标签和元素有什么区别？

    我注意到大多数人都使用这些词HTML 标签 and HTML 元素可以互换 但它们之间有什么区别呢 我的看法是 标签位于源代码中 元素是 DOM 中 由浏览器 处理的标签 我错了吗 HTML 标签只是打开或关闭实体 例如 p and p 称

  • Firebase addValueEventListener 未触发

    上周左右 我在从 Firebase 数据库检索数据时遇到了麻烦 我在 StackOverflow 和 google 上搜索了这个问题的答案 要么我不理解所提供的解决方案 要么它们根本不适合我 一切都被成功调用 直到我到达 valueEven

  • 从 HTTP 响应中获取 JSON 对象

    我想得到一个JSON从 Http 对象获取响应 这是我当前的 Http get 代码 protected String doInBackground String params HttpClient client new DefaultHt

  • 我们如何使用 asp.net、webservice 和 sql 数据库集成 jQuery 自动完成？

    我正在尝试实现 jQuery Autocomplete 和 ASP NET 给出的代码 但无法集成它 因为您正在使用亚音速来查询数据库 那么你能告诉我如何使用C 查询sql数据库并将查询结果从asp net中的Web服务绑定到插件吗 这是一

  • 动态选择 TableLayout 中的 tableRow

    我正在动态创建一个包含许多 TableRows 的 Tablelayout 例如 for int i 0 i

  • Spring 4 中的 DeferredResult 支持 Servlet 3.1 (Read|Write)Listener 吗？

    我正在读杰威文章关于 Spring 对 Servlet 的异步支持 有趣的部分是 如果您的服务预计会接收大量请求或响应主体 特别是如果客户端写入或读取速度较慢 那么您将受益于使用 Servlet 3 1 中引入的非阻塞 IO 功能 如前所述

  • 删除 R 矩阵中所有数据均为 NA 的行[重复]

    这个问题在这里已经有答案了 可能的重复 在 R 中删除数据文件的空行 如何从矩阵或数据框中删除行all该行中的元素是否为 NA 所以要从中得到 1 2 3 1 1 6 11 2 NA NA NA 3 3 8 13 4 4 NA NA 5 5

  • 为什么我不能使用 HttpContext 或 HttpCookie？ （ASP.NET核心1.0）

    为什么我不能使用HttpContext or HttpCookie 有什么特殊用途吗 我的实际使用情况 using System using System Collections Generic using System Linq usin

  • 元描述/标签不起作用

    我有一个非常奇怪的问题 元标记之前在我的主页上工作 但现在由于某种原因不再工作了 它在谷歌中没有正确显示标题 也没有我正在寻找的描述 我感觉元标记以某种方式被阻止了 我所做的唯一一件大事就是对网站进行 gzip 压缩 但我不确定这就是问题所

  • 用python将多页pdf文件分割成多个pdf文件？

    我想获取一个多页 pdf 文件并为每页创建单独的 pdf 文件 我已经下载了报告实验室并浏览了文档 但它似乎旨在生成 pdf 我还没有看到任何有关处理 PDF 文件本身的内容 有没有一种简单的方法可以在 python 中做到这一点 from

  • 如何找出从共享对象导出的所有符号？

    我有一个共享对象 dll 我如何找出从中导出的所有符号 您是否有 共享对象 通常是 AIX 上的共享库 UNIX 共享库或 Windows DLL 这些都是不同的事情 你的问题将它们全部混为一谈 对于 AIX 共享对象 请使用dump Tv

  • R bookdown 中标题前的封面页和版权声明？

    早在三月份 我就在 R bookdown 渲染的 pdf 文档中包含封面页提出了一个问题并得到了答案 R bookdown 封面页和附录 我尝试了该解决方案并得出以下结果 在index rmd yaml中使用 output pdf docu

  • GEKKO 和 Scipy.optimize 导致非线性参数估计结果不同

    我正在学习如何使用 GEKKO 来解决参数估计问题以及 作为第一步 我正在开发我遇到的示例问题 之前使用 Scipy 最小化例程实现 这些有 已按照 APMonitor com 中提供的信息和 内提供的课程 目前的问题是间歇式反应器 甲醇转

热门标签