Python方程组拟合失效：参数出现负值问题排查

尝试将γ与Cₛ的实验数据拟合到Frumkin等温相关方程组中，拟合参数包括gamma₀、Gamma_inf、a、K，所有参数均不能为负值。代码运行无报错，但拟合失效，部分参数出现负值，无法定位问题。

公式1：Frumkin平衡方程（x与Cₛ的关系）
公式2：表面张力γ与x的关系

原代码：

import numpy as np
from scipy.optimize import curve_fit, fsolve
import matplotlib.pyplot as plt
import pandas as pd

# Prompt the user for the CSV file location
csv_file = input("Enter the path to your CSV file: ")

# Read the CSV file into a DataFrame
df = pd.read_csv(csv_file,header=None)

# Alternatively, select each column by index
# Remember that indices start from 0
Cs_data = df.iloc[:, 0]
gamma_data = df.iloc[:, 1]

# Define the Frumkin equilibrium function to solve for x
def frumkin_x(C_s, a, K, x_guess=0.5):
    # Function to find the root of
    func = lambda x: x - C_s / (C_s + a * np.exp(K * x))
    x_solution, = fsolve(func, x_guess)
    return x_solution

# Define the equilibrium isotherm function
def equilibrium_isotherm(C_s, gamma_0, Gamma_inf, a, K):
    x = np.array([frumkin_x(cs, a, K) for cs in C_s])
    return gamma_0 + Gamma_inf * R * T * (np.log(1 - x) - 0.5 * K * x**2)

# Constants
R = 8.314  # Universal gas constant, J/(mol*K)
T = 298    # Temperature, K

# Initial parameter guesses for curve_fit
# gamma_0, Gamma_inf, a, K
initial_guesses = [72, 0.0000000004, 0.00000007, 0]

params_opt, params_cov = curve_fit(lambda Cs, gamma_0, Gamma_inf, a, K: equilibrium_isotherm(Cs, gamma_0, Gamma_inf, a, K), Cs_data, gamma_data, p0=initial_guesses)

# Print optimized parameters
print("Optimized Parameters:")
print("gamma_0:", params_opt[0])
print("Gamma_inf:", params_opt[1])
print("a (alpha_0/beta_0):", params_opt[2])
print("K:", params_opt[3])

# Plot the original data and fitted curve for visualization
Cs_fit = np.linspace(min(Cs_data), max(Cs_data), 100)
gamma_fit = equilibrium_isotherm(Cs_fit, *params_opt)

plt.scatter(Cs_data, gamma_data, label='Data')
plt.plot(Cs_fit, gamma_fit, label='Fitted Curve', color='red')
plt.xlabel('Sublayer Concentration (Cs)')
plt.ylabel('Surface Tension (gamma)')
plt.legend()
plt.show()

核心问题原因 #

• 无参数非负约束：curve_fit默认允许参数取任意值，优化过程中容易出现不符合物理意义的负值。
• 初始猜测不合理：部分参数初始值数量级偏离实际范围，K设为0可能让优化陷入局部最优。
• Frumkin方程求解不稳定：固定初始猜测值0.5，当参数或Cₛ变化时，可能无法收敛到正确的覆盖度x。

具体修正方案 #

1. 添加参数非负边界约束 #

在curve_fit中使用bounds参数强制所有参数≥0，同时设置合理上限：

# 设置参数边界：下限全为0，上限根据物理意义设定
bounds = (0, [100, 1e-8, 1e-6, 10])
params_opt, params_cov = curve_fit(
    lambda Cs, gamma_0, Gamma_inf, a, K: equilibrium_isotherm(Cs, gamma_0, Gamma_inf, a, K),
    Cs_data, gamma_data,
    p0=initial_guesses,
    bounds=bounds
)

2. 优化初始猜测值 #

调整初始值至更符合物理意义的范围：

# gamma0为纯溶剂表面张力，Gamma_inf为饱和吸附量，a为平衡常数相关参数，K为相互作用参数
initial_guesses = [72, 1e-9, 1e-7, 1]

3. 提升Frumkin方程求解稳定性 #

根据Cₛ动态调整x的初始猜测，并限制解的物理范围：

def frumkin_x(C_s, a, K):
    # 低浓度时x接近0，高浓度时x接近1，动态设置初始猜测
    x_guess = 0.1 if C_s < a else 0.9
    func = lambda x: x - C_s / (C_s + a * np.exp(K * x))
    x_solution, = fsolve(func, x_guess)
    # 强制x在(0,1)范围内，避免数值异常
    return np.clip(x_solution, 1e-6, 1 - 1e-6)

4. 避免对数数值异常 #

在计算np.log(1-x)时添加小量，防止出现log(0)错误：

def equilibrium_isotherm(C_s, gamma_0, Gamma_inf, a, K):
    x = np.array([frumkin_x(cs, a, K) for cs in C_s])
    log_term = np.log(1 - x + 1e-12)  # 添加1e-12避免log(0)
    return gamma_0 + Gamma_inf * R * T * (log_term - 0.5 * K * x**2)

完整修正代码 #

import numpy as np
from scipy.optimize import curve_fit, fsolve
import matplotlib.pyplot as plt
import pandas as pd

# 读取CSV数据
csv_file = input("请输入CSV文件路径: ")
df = pd.read_csv(csv_file, header=None)
Cs_data = df.iloc[:, 0]
gamma_data = df.iloc[:, 1]

# 定义Frumkin平衡方程求解覆盖度x
def frumkin_x(C_s, a, K):
    x_guess = 0.1 if C_s < a else 0.9
    func = lambda x: x - C_s / (C_s + a * np.exp(K * x))
    x_solution, = fsolve(func, x_guess)
    return np.clip(x_solution, 1e-6, 1 - 1e-6)

# 定义表面张力拟合函数
def equilibrium_isotherm(C_s, gamma_0, Gamma_inf, a, K):
    x = np.array([frumkin_x(cs, a, K) for cs in C_s])
    log_term = np.log(1 - x + 1e-12)
    return gamma_0 + Gamma_inf * R * T * (log_term - 0.5 * K * x**2)

# 常数定义
R = 8.314  # J/(mol*K)
T = 298    # K

# 优化初始猜测值
initial_guesses = [72, 1e-9, 1e-7, 1]

# 设置参数非负边界
bounds = (0, [100, 1e-8, 1e-6, 10])

# 执行拟合
params_opt, params_cov = curve_fit(
    lambda Cs, gamma_0, Gamma_inf, a, K: equilibrium_isotherm(Cs, gamma_0, Gamma_inf, a, K),
    Cs_data, gamma_data,
    p0=initial_guesses,
    bounds=bounds
)

# 输出结果
print("优化后参数:")
print("gamma_0:", params_opt[0])
print("Gamma_inf:", params_opt[1])
print("a (alpha_0/beta_0):", params_opt[2])
print("K:", params_opt[3])

# 绘图可视化
Cs_fit = np.linspace(min(Cs_data), max(Cs_data), 100)
gamma_fit = equilibrium_isotherm(Cs_fit, *params_opt)

plt.scatter(Cs_data, gamma_data, label='实验数据')
plt.plot(Cs_fit, gamma_fit, label='拟合曲线', color='red')
plt.xlabel('亚层浓度 (Cs)')
plt.ylabel('表面张力 (gamma)')
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者Woojin Jung