使用least_squares求解镍-氨体系化学形态时结果不准确的问题

方程组定义 #

import numpy as np
from scipy.optimize import least_squares

def equations(x):
    Ni0, Ni1, Ni2, Ni3, Ni4, Ni5, Ni6, NH3 = x
    eq1 = np.log10(Ni1)-np.log10(Ni0)-np.log10(NH3)-K1
    eq2 = np.log10(Ni2)-np.log10(Ni1)-np.log10(NH3)-K2
    eq3 = np.log10(Ni3)-np.log10(Ni2)-np.log10(NH3)-K3
    eq4 = np.log10(Ni4)-np.log10(Ni3)-np.log10(NH3)-K4
    eq5 = np.log10(Ni5)-np.log10(Ni4)-np.log10(NH3)-K5
    eq6 = np.log10(Ni6)-np.log10(Ni5)-np.log10(NH3)-K6
    eq7 = Ni0 + Ni1 + Ni2 + Ni3 + Ni4 + Ni5 + Ni6 - Ni_tot
    eq8 = Ni1 + 2*Ni2 + 3*Ni3 + 4*Ni4 + 5*Ni5 + 6*Ni6 - NH3_tot[i]
    
    return np.array([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8])

求解代码 #

通过循环计算不同NH3_tot值下的参数，观察各物种浓度变化，求解代码如下：

# Arrays for data storage
NH3_tot = np.arange(0, 6, 0.1)
Ni0_values = []
Ni1_values = []
Ni2_values = []
Ni3_values = []
Ni4_values = []
Ni5_values = []
Ni6_values = []
Solutions = []

# Speciation determination with least_square at various NH3_tot
for i in range(len(NH3_tot)):
    x0 = [1,0.1,0.1,0.1,0.1,0.1,0.1,0.1]
    result = least_squares(equations, x0, bounds=([0,0,0,0,0,0,0,0],[Ni_tot,Ni_tot,Ni_tot,Ni_tot,Ni_tot,Ni_tot,Ni_tot,Ni_tot]))
    
    Solution = result.x
    Solutions.append(Solution)
    Ni0_value,Ni1_value,Ni2_value,Ni3_value,Ni4_value,Ni5_value,Ni6_value, NH3 = Solution
    
    # Extract concentrations of species
    Ni0_values.append(Ni0_value)
    Ni1_values.append(Ni1_value)
    Ni2_values.append(Ni2_value)
    Ni3_values.append(Ni3_value)
    Ni4_values.append(Ni4_value)
    Ni5_values.append(Ni5_value)
    Ni6_values.append(Ni6_value)

遇到的问题 #

部分NH3_tot值（通常是低和高值）的结果完全错误，尤其不满足eq.7，物种总和超过Ni_tot。尝试过调整边界、更改方程写法、切换求解方法，但均无明显改善，尽管函数返回成功终止的信息。

后续尝试 #

借助@jlandercy的修正，使用各配合物比例的数学表达式（alpha_i）降低了曲线形状的误差，但配体浓度轴仍存在异常，无法得到符合预期的结果（预期为Ni总浓度1mol/L、NH3总浓度0-12mol/L时的物种分布）。

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