主页 / coding / 问题 / 79221744
Accepted
Ronny M
Asked: 2024-11-25 13:35:44 +0800 CST

Gekko 的 IPOPT 解决方案与 Excel-Solver 中使用的 GRG 算法的比较

  • 772

目的是根据乙烷蒸汽气化反应产物和反应物（蒸汽+C2H6 的摩尔比为 4:1）的吉布斯生成能，计算 1000K 时混合物的热力学平衡组成，如下所示：

from gekko import GEKKO     
m = GEKKO(remote=True)
x = m.Array(m.Var,13,value=1,lb=1e-27,ub=20.0)
H2,H2O,CO,O2,CO2,CH4,C2H6,C2H4,C2H2,lamda1,lamda2,lamda3,summe= x
H2.value = 3.0
H2O.value = 1.0
CO.value = 0.5
O2.value = 0.001
CO2.value = 0.5
CH4.value = 0.1
C2H6.value = 0.000215
C2H4.value = 0.00440125
C2H2.value = 0.0041294
summe.value = 8.0
lamda1.value = 1.0
lamda2.value = 1.0
lamda3.value = 1.0
eq1 = m.Param(value=14)
eq2 = m.Param(value=4)
eq3 = m.Param(value=2)
summe = m.Var(H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2)
lamda2 = m.Var((-1)*m.log(H2 / summe) / 2)
lamda1 = m.Var(46.03 / 1.9872 - m.log(H2O / summe) + 2 * lamda2)
lamda3 = m.Var(47.942 / 1.9872 - m.log(CO / summe) + lamda1)
m.Equation(m.exp(-4.61 / 1.9872 - 4 * lamda2 - lamda3) * summe == CH4)
m.Equation(m.exp(-28.249 / 1.9872 - 4 * lamda2 - 2 * lamda3) * summe == C2H4)
m.Equation(m.exp(-40.604 / 1.9872 - 2 * lamda2 - 2 * lamda3) * summe == C2H2)
m.Equation(m.exp(-26.13 / 1.9872 - 6 * lamda2 - 2 * lamda3) * summe == C2H6)
m.Equation(m.exp(94.61 / 1.9872 - 2 * lamda1 - lamda3) * summe == CO2)
m.Equation(m.exp(-2 * lamda1) * summe == O2)
m.Equation(2 * CO2 + CO + 2 * O2 + H2O  == eq2)
m.Equation(4 * CH4 + 4 * C2H4 + 2 * C2H2 + 2 * H2 + 2 * H2O + 6 * C2H6  == eq1)
m.Equation(CH4 + 2 * C2H4 + 2 * C2H2 + CO2 + CO + 2 * C2H6  == eq3)
m.Minimize((summe-(H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2))**2)
m.options.IMODE = 3 #IPOPT
m.options.MAX_ITER = 1000
m.options.OTOL = 1e-10
m.options.RTOL = 1e-10
m.solve()
print('x: ', x)
print('Objective: ',m.options.OBJFCNVAL)

EXIT: Optimal Solution Found. 
The solution was found. 
The final value of the objective function is    81.0000000000000     

---------------------------------------------------
Solver         :  IPOPT (v3.12)
Solution time  :   1.169999998819549E-002 sec
Objective      :    81.0000000000000     
Successful solution
---------------------------------------------------

x:  [[5.797458326] [1.202541674] [1.202541674] [1.6791020526e-21]
[0.79745832603] [1.6503662452e-22] [3.2694282596e-27] [1.1255712085e-27]
[1e-27] [2.185089969] [2.185089969] [2.185089969] [9.5605307091]]

然后使用 GRG 作为优化器：

Excel 工作表输出

两种解决方案仍然不同。有趣的是，（约束）根查找器的输出仍然不同，但同时收敛：

Total number of equations:  13  
Number of implicit equations:  4  
Number of explicit equations:  9  
Solution method  CONSTRAINED  
Convergence tolerance:  1e-07  
# of iterations used:  8  

CO  1.685236
H2  5.118105
H2O  1.387357
SUM  8.899115  
C2H2  0.0041294  
C2H4  0.0044224  
C2H6  0.0092403  
CH4  0.2269213  
CO2  0.0522585  
lamda1  1.537016  
lamda2  0.2765838  
lamda3  2.53957  
O2  0.4114448

以下是各种解决方案的比较：

解向量 x 壁虎 EXCEL表格 寻根器 GEKKO-最终版
H2：1 5.797458326 5.344360712 5.118105 5.219
水：2 1.202541674 1.521995663 1.387357 1.681
二氧化碳：3 1.202541674 1.388351672 1.685236 1.581
氧气：4 1.6791020526e-21 5.82E-21 0.4114448 1e-5
二氧化碳：5 0.79745832603 0.544826 0.0522585 0.3693
通道4: 6 1.6503662452e-22 0.0668213 0.2269213 0.05
乙烷：7 3.2694282596e-27 1.70E-07 0.0092403 1e-5
乙烷：8 1.1255712085e-27 9.74E-08 0.0044224 1e-5
乙烷：9 1e-27 3.25E-10 0.0041294 1e-5
lamda1：10 2.185089969 24.3878385 1.537016
lamda2：11 2.185089969 0.253110984 0.2765838
lamda3：12 2.185089969 1.558979624 2.53957
总数：13 9.5605307091 8.866356107 8.899115 8.90
python

1 个回答

  1. Best Answer

    变量lamda1、lamda2、lamda3、summe被定义两次，与这些变量相关的方程仅用于初始化变量的第二次定义。

    H2,H2O,CO,O2,CO2,CH4,C2H6,C2H4,C2H2,lamda1,lamda2,lamda3,summe= x
summe = m.Var(H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2)
lamda2 = m.Var((-1)*m.log(H2 / summe) / 2)
lamda1 = m.Var(46.03 / 1.9872 - m.log(H2O / summe) + 2 * lamda2)
lamda3 = m.Var(47.942 / 1.9872 - m.log(CO / summe) + lamda1)

    将它们切换到Intermediate()定义是解决问题的简单方法。

    summe = m.Intermediate(H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2)
lamda2 = m.Intermediate((-1)*m.log(H2 / summe) / 2)
lamda1 = m.Intermediate(46.03 / 1.9872 - m.log(H2O / summe) + 2 * lamda2)
lamda3 = m.Intermediate(47.942 / 1.9872 - m.log(CO / summe) + lamda1)

    完整脚本如下：

    from gekko import GEKKO   
import numpy as np  
m = GEKKO(remote=True)
x = m.Array(m.Var,9,value=1,lb=1e-27,ub=20.0)
H2,H2O,CO,O2,CO2,CH4,C2H6,C2H4,C2H2=x
H2.value = 3.0
H2O.value = 1.0
CO.value = 0.5
O2.value = 0.001
CO2.value = 0.5
CH4.value = 0.1
C2H6.value = 0.000215
C2H4.value = 0.00440125
C2H2.value = 0.0041294
eq1 = m.Param(value=14)
eq2 = m.Param(value=4)
eq3 = m.Param(value=2)
summe = m.Intermediate(H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2)
lamda2 = m.Intermediate((-1)*m.log(H2 / summe) / 2)
lamda1 = m.Intermediate(46.03 / 1.9872 - m.log(H2O / summe) + 2 * lamda2)
lamda3 = m.Intermediate(47.942 / 1.9872 - m.log(CO / summe) + lamda1)
m.Equation(m.exp(-4.61 / 1.9872 - 4 * lamda2 - lamda3) * summe == CH4)
m.Equation(m.exp(-28.249 / 1.9872 - 4 * lamda2 - 2 * lamda3) * summe == C2H4)
m.Equation(m.exp(-40.604 / 1.9872 - 2 * lamda2 - 2 * lamda3) * summe == C2H2)
m.Equation(m.exp(-26.13 / 1.9872 - 6 * lamda2 - 2 * lamda3) * summe == C2H6)
m.Equation(m.exp(94.61 / 1.9872 - 2 * lamda1 - lamda3) * summe == CO2)
m.Equation(m.exp(-2 * lamda1) * summe == O2)
m.Equation(2 * CO2 + CO + 2 * O2 + H2O  == eq2)
m.Equation(4 * CH4 + 4 * C2H4 + 2 * C2H2 + 2 * H2 + 2 * H2O + 6 * C2H6  == eq1)
m.Equation(CH4 + 2 * C2H4 + 2 * C2H2 + CO2 + CO + 2 * C2H6  == eq3)
m.Minimize((summe-(H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2))**2)
m.options.IMODE = 3 #IPOPT
m.options.MAX_ITER = 1000
m.options.OTOL = 1e-10
m.options.RTOL = 1e-10
m.solve()
print('x: ', x)
print('Objective: ',m.options.OBJFCNVAL)
print(f'H2: {H2.value[0]}')
print(f'H2O: {H2O.value[0]}')
print(f'CO: {CO.value[0]}')
print(f'O2: {O2.value[0]}')
print(f'CO2: {CO2.value[0]}')
print(f'CH4: {CH4.value[0]}')
print(f'C2H6: {C2H6.value[0]}')
print(f'C2H4: {C2H4.value[0]}')
print(f'C2H2: {C2H2.value[0]}')

    现在，目标值0具有0求根所需的自由度。解决方案是：

     The solution was found.
 
 The final value of the objective function is   0.000000000000000E+000
 
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :   1.099999999860302E-002 sec
 Objective      :   0.000000000000000E+000
 Successful solution
 ---------------------------------------------------
 
x:  [[5.0] [2.0] [2.0] [1.0421441742e-21] [3.9704669403e-23]
 [2.1920286233e-23] [1e-27] [1e-27] [1e-27]]
Objective:  0.0
H2: 5.0
H2O: 2.0
CO: 2.0
O2: 1.0421441742e-21
CO2: 3.9704669403e-23
CH4: 2.1920286233e-23
C2H6: 1e-27
C2H4: 1e-27
C2H2: 1e-27

    对于水煤气变换反应来说，这个解决方案看起来不太正确。我不得不回去重新回忆一下 1000 K 下含有乙烷 (C₂H₆) 和蒸汽 (H₂O) 的反应混合物的平衡组成。

    1. 蒸汽重整反应

    氢气产生的主要反应：

    C2H6 + 2H2O → 2CO + 5H2
    • 用途：在蒸汽存在下将碳氢化合物转化为一氧化碳（CO）和氢气（H₂）。
    • 热力学效应：高度吸热，在高温下更有利。

    2. 水煤气变换反应

    一氧化碳与蒸汽的二次反应：

    CO + H2O ↔ CO2 + H2
    • 目的：将 CO 转化为 CO₂ 并产生额外的 H₂。
    • 热力学效应：放热，在较低温度下有利，但即使在高温下也会达到平衡。

    3. 其他碳氢化合物反应

    该系统还包括碳氢化合物的分解和部分氧化：

    • 乙烯形成（C₂H₄）： 
      C2H6 → C2H4 + H2
    • 乙炔形成（C₂H₂）： 
      C2H4 → C2H2 + H2

    4. 积碳

    在某些条件下会发生碳沉积：

    CO → C (solid) + CO2

    或者

    CH4 → C (solid) + 2H2
    • 目的：表示碳以固体形式沉积的不良副反应。
    • 热力学效应：根据系统条件，可以在较低的氢碳比或氧碳比下发生。

    热力学建模

    通过最小化系统的吉布斯自由能来计算平衡组成：

    G = Σ (ni * (gi° + RT * ln(ni / ntotal)))

    在哪里：

    • ni： 物质的摩尔量i。
    • gi°：物种的标准吉布斯生成自由能i。
    • ntotal：所有物种的总摩尔数。

    限制：

    • 碳、氢和氧的质量平衡。
    • 对摩尔数量进行限制以避免非物理解决方案。

    以下是基于这些方程的完整脚本：

    from gekko import GEKKO

# Create GEKKO model
m = GEKKO(remote=True)

# Define variables for species molar amounts
n_C2H6 = m.Var(value=1.0, lb=0, ub=1)  # Ethane
n_H2O = m.Var(value=4.0, lb=0, ub=1)   # Water
n_CO = m.Var(value=0.0, lb=0, ub=1)    # Carbon monoxide
n_CO2 = m.Var(value=0.0, lb=0, ub=1)   # Carbon dioxide
n_H2 = m.Var(value=0.0, lb=0, ub=1)    # Hydrogen
n_CH4 = m.Var(value=0.0, lb=0, ub=1)   # Methane
n_C = m.Var(value=0.0, lb=0, ub=1)     # Solid carbon
n_C2H4 = m.Var(value=0.0, lb=0, ub=1)  # Ethylene
n_C2H2 = m.Var(value=0.0, lb=0, ub=1)  # Acetylene

# Universal gas constant
R = 1.9872 # cal/(K mol)

# Gibbs free energies of formation at 1000K
gibbs_energies = {
    "C2H6": -26.13,    # Ethane
    "H2O": 46.03,      # Water (vapor)
    "CO": 47.942,      # Carbon monoxide
    "CO2": 94.61,      # Carbon dioxide
    "H2": 0.0,         # Hydrogen gas
    "CH4": -4.61,      # Methane
    "C": 0.0,          # Solid carbon
    "C2H4": -28.249,   # Ethylene
    "C2H2": -40.604    # Acetylene
}


# Define total moles
total_moles = (
    n_C2H6 + n_H2O + n_CO + n_CO2 + n_H2 + n_CH4 + n_C + n_C2H4 + n_C2H2
)

# Define the objective function (Gibbs free energy minimization)
m.Obj(
    n_C2H6 * (gibbs_energies["C2H6"]/R + m.log(n_C2H6 / total_moles + 1e-10)) +
    n_H2O * (gibbs_energies["H2O"]/R + m.log(n_H2O / total_moles + 1e-10)) +
    n_CO * (gibbs_energies["CO"]/R + m.log(n_CO / total_moles + 1e-10)) +
    n_CO2 * (gibbs_energies["CO2"]/R + m.log(n_CO2 / total_moles + 1e-10)) +
    n_H2 * (gibbs_energies["H2"]/R + m.log(n_H2 / total_moles + 1e-10)) +
    n_CH4 * (gibbs_energies["CH4"]/R + m.log(n_CH4 / total_moles + 1e-10)) +
    n_C * (gibbs_energies["C"]/R + m.log(n_C / total_moles + 1e-10)) +
    n_C2H4 * (gibbs_energies["C2H4"]/R + m.log(n_C2H4 / total_moles + 1e-10)) +
    n_C2H2 * (gibbs_energies["C2H2"]/R + m.log(n_C2H2 / total_moles + 1e-10))
)

# Mass balance constraints
# Carbon balance: 4(C2H6) = CO + CO2 + CH4 + C + C2H4 + C2H2
m.Equation(4*n_C2H6 == n_CO + n_CO2 + n_CH4 + n_C + n_C2H4 + n_C2H2)

# Oxygen balance: 3*H2O = CO + CO2
m.Equation(3*n_H2O == n_CO + n_CO2)

# Hydrogen balance with 4:1 ratio: C2H6 + 4*H2O = 2*H2 + CH4 + C2H4 + C2H2
m.Equation(4*n_H2O + n_C2H6 == 2*n_H2 + n_CH4 + n_C2H4 + n_C2H2)

# Total mole constraint for mole fraction result
m.Equation(total_moles == 1)

# Solve the model
m.solve(disp=True)

# Extract and display results
results = {
    "C2H6 (mol)": n_C2H6.value[0],
    "H2O (mol)": n_H2O.value[0],
    "CO (mol)": n_CO.value[0],
    "CO2 (mol)": n_CO2.value[0],
    "H2 (mol)": n_H2.value[0],
    "CH4 (mol)": n_CH4.value[0],
    "C (mol)": n_C.value[0],
    "C2H4 (mol)": n_C2H4.value[0],
    "C2H2 (mol)": n_C2H2.value[0]
}
[print(f'{r}: {results[r]:.4}') for r in results]

    结果：

    C2H6 (mol): 0.1993
H2O (mol): 0.003599
CO (mol): 0.0108
CO2 (mol): 0.0
H2 (mol): 0.0
CH4 (mol): 2.915e-09
C (mol): 0.5726
C2H4 (mol): 0.0004254
C2H2 (mol): 0.2133
    • 1

相关问题