解决CO₂-H₂反应动力学NLP模拟中非物理结果及速率常数异常问题的技术问询

阿华AIGC实验室

2026-4-14

我最近在研究CO₂和H₂反应体系的动力学问题，模拟温度升高时反应混合物的变化。我用Gekko工具来求解这个非线性规划（NLP）模型，但遇到了两个棘手的问题：一是模拟结果出现了非物理的情况，特别是CO的浓度变成了负值；二是当我调整了反应方程的形式后，拟合效果有所提升，但得到的动力学速率常数k却都是接近零的噪声值。有没有大佬能给我一些改进代码的建议呀？万分感谢！

当前模拟结果展示

反应体系浓度模拟与测量值对比（出现CO浓度负值）：
物种相对浓度模拟结果：
调整方程后速率常数结果（k值接近零）：

方程调整后的情况

当我改用如下形式的单一反应方程时，拟合效果好了很多，但得到的动力学速率常数k都是接近零的噪声值：

m.Equation(0 == -(k1 * S * CO2 * H2 * (1 - ((CO * H2O) / (CO2 * H2)) / K1eq_param)))

完整模拟代码

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
import csv

CH4_meas = np.array([0.1361,0.1885,0.2481,0.3341,0.4266]) # and further values...
CO_meas = np.array([0.0281,0.0376,0.0476,0.0628,0.0817])
H2O_meas = np.array([0.1201,0.234,0.3821,0.5821,0.8274])
H2_meas = np.array([62.6176,62.3346,62.0223,61.7904,61.5901])
CO2_meas = np.array([30.4669,30.3408,30.3541,30.3541,30.2984])


K1eq = np.array([0.01690,0.01881,0.02088,0.02312,0.02555])
K2eq = np.array([5.7990E+08,3.3957E+08,2.0132E+08,1.2079E+08,7.3315E+07])
K3eq = np.array([9.7992E+06,6.3860E+06,4.2031E+06,2.7930E+06,1.8732E+06])


CH4_values = []
CO_values = []
H2O_values = []
H2_values = []
CO2_values = []
k1_values = []
k2_values = []
k3_values = []
S=25

for i in range(len(CH4_meas)):
    m = GEKKO(remote=False)

    # Initial guess for the rate constants, but when I set lb=0, no solution can be found
    k1 = m.FV(value=0.0025, lb=None, ub=None)
    k2 = m.FV(value=0.0025, lb=None, ub=None)
    k3 = m.FV(value=0.2, lb=None, ub=None)
    k1.STATUS = 1
    k2.STATUS = 1
    k3.STATUS = 1

    # Controlled Variables (CV) with measurement normalization
    CH4 = m.CV(CH4_meas[i] / (H2_meas[0] + CO2_meas[0]))
    CO = m.CV(CO_meas[i] / (H2_meas[0] + CO2_meas[0]))
    H2O = m.CV(H2O_meas[i] / (H2_meas[0] + CO2_meas[0]))
    H2 = m.CV(H2_meas[i] / (H2_meas[0] + CO2_meas[0]))
    CO2 = m.CV(CO2_meas[i] / (H2_meas[0] + CO2_meas[0]))    
    CH4.FSTATUS = 1
    CO.FSTATUS = 1
    H2O.FSTATUS = 1
    H2.FSTATUS = 1
    CO2.FSTATUS = 1

    # Parameters
    K1eq_param = m.Param(value=K1eq[i])
    K2eq_param = m.Param(value=K2eq[i])
    K3eq_param = m.Param(value=K3eq[i])

    # Steady-state equations
    m.Equation(CO2 == -(k1 * S * CO2 * H2 * (1 - ((CO * H2O) / (CO2 * H2)) / K1eq_param)) - k3 * S * CO2 * (1 - ((CH4 * H2O**2) / (CO2 * H2**4)) / K3eq_param))
    m.Equation(H2 == -(k1 * S * CO2 * H2 * (1 - ((CO * H2O) / (CO2 * H2)) / K1eq_param) - 3 * k2 * S * CO * (1 - ((CH4 * H2O) / (CO * H2**3)) / K2eq_param) - 4 * k3 * S * CO2 * (1 - ((CH4 * H2O**2) / (CO2 * H2**4)) / K3eq_param)))
    m.Equation(CO == (k1 * S * CO2 * H2 * (1 - ((CO * H2O) / (CO2 * H2)) / K1eq_param) - k2 * S * CO * (1 - ((CH4 * H2O) / (CO * H2**3)) / K2eq_param)))
    m.Equation(H2O == (k1 * S * CO2 * H2 * (1 - ((CO * H2O) / (CO2 * H2)) / K1eq_param) + k2 * S * CO * (1 - ((CH4 * H2O) / (CO * H2**3)) / K2eq_param) + k3 * S * CO2 * (1 - ((CH4 * H2O**2) / (CO2 * H2**4)) / K3eq_param)))
    m.Equation(CH4 == (k2 * S * CO * (1 - ((CH4 * H2O) / (CO * H2**3)) / K2eq_param) + k3 * S * CO2 * (1 - ((CH4 * H2O**2) / (CO2 * H2**4)) / K3eq_param)))

    # Objective function for regression fitting
    m.Obj((CH4.value - CH4_meas[i])**2 + (CO.value - CO_meas[i])**2 + (H2O.value - H2O_meas[i])**2 +
           (H2.value - H2_meas[i])**2 + (CO2.value - CO2_meas[i])**2)
  

    # Set the mode for steady-state optimization
    m.options.IMODE = 3  # Steady-state optimization
    m.options.EV_TYPE = 1 #2 = Sum of squared error, 1 = Sum of absoluted error

    m.options.MAX_ITER = 300
    m.options.OTOL = 1e-6
    m.options.RTOL = 1e-6

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

    # Store results
    CH4_values.append(CH4.value[0])
    CO_values.append(CO.value[0])
    H2O_values.append(H2O.value[0])
    H2_values.append(H2.value[0])
    CO2_values.append(CO2.value[0])
    k1_values.append(k1.value[0])
    k2_values.append(k2.value[0])
    k3_values.append(k3.value[0])
    
    print("k1 = " + str(k1.value[0]))
    print("k2 = " + str(k2.value[0]))
    print("k3 = " + str(k3.value[0]))


with open('OptimGekko_Methane_RWGS.csv', 'w', newline='') as csvfile:
    csvwriter = csv.writer(csvfile)
    csvwriter.writerow(['T', 'k1', 'k2', 'k3','CO2','H2','CO','H2O','CH4'])
    for i in range(len(k1_values)):
        csvwriter.writerow([time[i], k1_values[i], k2_values[i], k3_values[i], CO2_values[i], H2_values[i],CO_values[i], H2O_values[i], CH4_values[i]])


plt.figure(figsize=(10, 12))


plt.subplot(5, 1, 1)
plt.plot(time[:len(CH4_values)], CH4_values, label='CH4_simul')
plt.plot(time[:len(CH4_meas)], CH4_meas/(30.47+62.62), 'ro', label='CH4_meas')
plt.ylabel('c/c0 (-)')
plt.xlabel('T (°C)')
plt.legend()
plt.grid()
plt.title('CH4')
plt.subplot(5, 1, 2)
plt.plot(time[:len(CO_values)], CO_values, label='CO_simul')
plt.plot(time[:len(CO_meas)], CO_meas/(30.47+62.62), 'ro', label='CO_meas')
plt.ylabel('c/c0 (-)')
plt.xlabel('T (°C)')
plt.legend()
plt.grid()
plt.title('CO')
plt.subplot(5, 1, 3)
plt.plot(time[:len(H2O_values)], H2O_values, label='H2O_simul')
plt.plot(time[:len(H2O_meas)], H2O_meas/(30.47+62.62), 'ro', label='H2O_meas')
plt.ylabel('c/c0 (-)')
plt.xlabel('T (°C)')
plt.legend()
plt.grid()
plt.title('H2O')
plt.subplot(5, 1, 4)
plt.plot(time[:len(H2_values)], H2_values, label='H2_simul')
plt.plot(time[:len(H2_meas)], H2_meas/(30.47+62.62), 'ro', label='H2_meas')
plt.ylabel('c/c0 (-)')
plt.xlabel('T (°C)')
plt.legend()
plt.grid()
plt.title('H2')
plt.subplot(5, 1, 5)
plt.plot(time[:len(CO2_values)], CO2_values, label='CO2_simul')
plt.plot(time[:len(CO2_meas)], CO2_meas/(30.47+62.62), 'ro', label='CO2_meas')
plt.ylabel('c/c0 (-)')
plt.xlabel('T (°C)')
plt.legend()
plt.grid()
plt.title('CO2')

plt.tight_layout()
plt.savefig('Optim_RWGS_CH4.png')
plt.show()

备注：内容来源于stack exchange，提问作者Ronny M