双位点Langmuir-Hinshelwood动力学模型拟合参数优选方法咨询

云原生阿岚

2026-7-17

双位点Langmuir-Hinshelwood模型参数优选问题

我正在将自行推导的双位点Langmuir-Hinshelwood模型（反应速率为反应物与产物分压的函数）拟合至实验数据。模型能够收敛，但初始参数的变化会导致拟合结果差异显著，而这些参数对我的分析至关重要。由于需要对比不同模型，模型与实验数据的拟合度依赖于K1、K2等参数。

我希望了解除拟合度外，是否存在其他方法判断某组参数更优；此前采用@jlandercy建议的MSE landscape计算分析参数间的凸性关系，但无法从中获取足够信息来确定K参数。

*更新：根据@jlandercy的建议修改了代码，目前代码可求解下方微分方程以优化产物气体流量，样本数据已更新。

拟合代码

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import least_squares
from scipy.optimize import fsolve
from scipy.integrate import odeint

# Step 1: Read experimental data
data1 = pd.read_excel("Local address of file")

# Extract the partial pressures (A, B, C) and the reaction rate
pp = data1.iloc[:49, 2:5].values
TT = data1.iloc[:49, 16].values
rr = data1.iloc[:49, 17].values
F0 = data1.iloc[:49, 5:10].values
W = data1.iloc[:49, 15].values
FF = data1.iloc[:49, 10:15].values

data = np.hstack((pp, F0, FF, W.reshape(-1,1), TT.reshape(-1,1)))
#-----------------------------------------------------------------------------

class Kinetic:
    
    """Langmuir-Hinshelwood Two-Site Kinetic Model Solver"""
    
    R = 8.314              # J/mol.K
    Ea = 45000  # J/mol  
    theta_g = [0.9, 0.1]    # 1
    
    w = np.linspace(0,1,100)
    
    @staticmethod
    def k3(flow, data, K1, K2, k4, K5, k0):
        """Arrhenius estimation for k3"""
        pA, pB, pC, F0_a, F0_b, F0_c, F0_ar, F0_h, FF_a, FF_b, FF_c, FF_ar, FF_h, mg, T = data
        return k0 * np.exp(-Kinetic.Ea / (Kinetic.R * T))

    @staticmethod
    def s(flow, data, K1, K2, k4, K5, k0):
        """Partial kinetic rate"""
        pA, pB, pC, F0_a, F0_b, F0_c, F0_ar, F0_h, FF_a, FF_b, FF_c, FF_ar, FF_h, mg, T = data
        
        #Calculating gas pressures at every position of PFR reactor
        pA_t = flow[0] / np.sum(flow)
        pB_t = flow[1] / np.sum(flow)
        pC_t = flow[2] / np.sum(flow)
        
        k3 = Kinetic.k3(flow, data, K1, K2, k4, K5, k0)
        return (K1 * K2 * k3 / k4) * pA_t * pB_t

    @staticmethod
    def equation(x, y, C, s):
        return 1/(1 + C + np.sqrt(s*y/x)) - x
    
    @staticmethod
    def system(theta, flow, data, K1, K2, k4, K5, k0):
        """Isothermal coverages system"""
        pA, pB, pC, F0_a, F0_b, F0_c, F0_ar, F0_h, FF_a, FF_b, FF_c, FF_ar, FF_h, mg, T = data
        t0, t1 = theta
        s = Kinetic.s(flow, data, K1, K2, k4, K5, k0)
        
        pA_t = flow[0] / np.sum(flow)
        pB_t = flow[1] / np.sum(flow)
        pC_t = flow[2] / np.sum(flow)
        
        C0 = K2 * pB_t
        C1 = K1 * pA_t + pC_t / K5
        return np.array([
            Kinetic.equation(t0, t1, C0, s),
            Kinetic.equation(t1, t0, C1, s),
        ])

    @staticmethod
    def _theta(flow, data, K1, K2, k4, K5, k0):
        """Single Isothermal coverages"""
        solution = fsolve(
            Kinetic.system,
            Kinetic.theta_g,
            args=(flow, data, K1, K2, k4, K5, k0),
            full_output=False
        )
        if not(all((solution >= 0) & (solution <= 1.0))):
            raise ValueError("Theta are not constrained to domain: %s" % solution)
        return solution

    @staticmethod
    def theta(flow, data, k1, k2, k4, k5, k0):
        """Isothermal coverages"""
        return np.apply_along_axis(Kinetic._theta, 0, flow, data, k1, k2, k4, k5, k0)

    @staticmethod
    def r1(flow, w, data, K1, K2, k4, K5, k0):
        """Global kinetic rate"""
        pA, pB, pC, F0_a, F0_b, F0_c, F0_ar, F0_h, FF_a, FF_b, FF_c, FF_ar, FF_h, mg, T = data
         
        pA_t = flow[0] / np.sum(flow)
        pB_t = flow[1] / np.sum(flow)
        pC_t = flow[2] / np.sum(flow)

        #print(pA)
        
        k3 = Kinetic.k3(flow, data, K1, K2, k4, K5, k0)
        t0, t1 = Kinetic.theta(flow, data, K1, K2, k4, K5, k0)
        v=np.array([-1,-1, 1, 0, 0])
        #v = v[:, np.newaxis]  # Convert v to a 2D array with shape (5, 1)
        return ((K1 * K2 * k3) * pA_t * pB_t * t0 * t1) * v
        

    @staticmethod
    def dFdm(data, K1, K2, k4, K5, k0):
        """Global kinetic rate"""
        residual =[]
        sq=[]
        for i in data:
            pA, pB, pC, F0_a, F0_b, F0_c, F0_ar, F0_h, FF_a, FF_b, FF_c, FF_ar, FF_h, mg, T = i
            F0 = np.array([F0_a, F0_b, F0_c, F0_ar, F0_h])
            FF = np.array([FF_a, FF_b, FF_c, FF_ar, FF_h]) # Final Experimental flows
            
            flow = odeint(Kinetic.r1, F0, Kinetic.w, args = (i, K1, K2, k4, K5, k0))
            flow_pred = flow[-1,2]  # Predicted final model flows
            res = (flow_pred - FF_c)*1000
            residual.append(res)
        return residual

bounds=((0, 0, 0, 0, 0), (np.inf, np.inf, np.inf, np.inf, np.inf))
initial_guesses = [8, 8, 50, 16, 6]
result = least_squares(lambda k: Kinetic.dFdm(data, *k), x0=initial_guesses, bounds=bounds, ftol=1e-14)

K1_fit, K2_fit, k4_fit, K5_fit, k0_fit = result.x
print(K1_fit, K2_fit, k4_fit, K5_fit, k0_fit)

内容的提问来源于stack exchange，提问作者Suyash Sachin Damir

双位点Langmuir-Hinshelwood动力学模型拟合参数优选方法咨询

相关资源 #

拟合代码 #

相关资源

拟合代码