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

编程问答 2022-09-30

问题描述

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


    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] #Initial conditions of species A,B,C,P
    ki= [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))
    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,args = (model.x),t_eval=(times)) 
    
    plt.plot(y_calc.t,y_calc.y.T)
    plt.plot(times,A_obs,'bo')
    plt.plot(times,'gx')
    plt.plot(times,P_obs,'rs')

    GEKKO Model

    import numpy as np
    import matplotlib.pyplot as plt
    from gekko import GEKKO
    
    #Experimental data
    times  = np.array([0.0,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,'rs')

解决方法

以下是一些建议：

  • 设置m.options.NODES=3或更高（最多6）以获得更好的集成精度。
  • 将Am，Cm，Pm设置为参数而不是变量。它们是固定输入。
  • 尝试不同的初始条件。可能有多个局部最小值。
  • 目标函数可能是平坦的，因此不同的参数值会给出相同的目标函数值。您可以test the parameter confidence intervals来查看数据是给出窄的还是宽的联合置信区间。

以下是修改后的结果：

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()
gekkopython

相关问答

matplotlib报错：AttributeError: module 'backend_interagg' has no attribute 'FigureCanvas'. Did you mean: 'FigureCanvasAgg'?
使用本地python环境可以成功执行 import pandas as pd impor...
gitlab登录失败，报错：This challenge page was accidentally cached by an intermediary and is no longer available.
设置时间 控制面板
后端开发常见错误
错误1：Request method ‘DELETE‘ not supported 错误还原：...
docker常见错误
错误1：启动docker镜像时报错：Error response from daemon:...
idea常见错误
错误1：private field ‘xxx‘ is never assigned 按Alt...
pip安装依赖失败
报错如下，通过源不能下载，最后警告pip需升级版本 Requirem...