求助：用Python least_squares求解带约束热力学非线性方程组失败

我需要求解一个热力学方程组，但作为Python初学者遇到了问题：使用scipy.optimize.least_squares求解带约束的非线性方程组时，结果始终是初始向量。以下是我的实现代码：

from scipy.optimize import least_squares
import numpy as np

def equations(x):
    eq1 = x[7] - 1 / (np.exp((5.24677 - 1598.673 / (60.0 - -46.424)) * np.log(10.0)))
    eq2 = x[4] - 1 / (np.exp((5.08354 - 1663.125 / (60.0 - -45.622)) * np.log(10.0)))
    eq3 = 0.0 - (0.5 * 0.4) + x[2] * x[9] - x[14] * (0.02 * 0.1)
    eq4 = 0.0 - (0.5 * 0.6) + x[2] * x[5] - x[12] * (0.02 * 0.1)
    eq5 = 0.0 - (0.4 * 0.6) + x[3] * x[0] + x[14] * (0.02 * 0.1)
    eq6 = 0.0 - (0.4 * 0.4) + x[3] * x[15] + x[12] * (0.02 * 0.1)
    eq7 = 0.0 - (x[10] - x[7] * x[13])
    eq8 = 0.0 - (x[8] - x[4] * x[11])
    eq9 = x[14] - (1.0 * x[6] * (x[13] - x[9]) + x[9] * (x[14] + x[12]))
    eq10 = x[12] - (1.0 * x[6] * (x[11] - x[5]) + x[5] * (x[14] + x[12]))
    eq11 = x[14] - (0.01 * x[1] * (x[0] - x[10]) + x[0] * (x[14] + x[12]))
    eq12 = x[12] - (0.01 * x[1] * (x[15] - x[8]) + x[15] * (x[14] + x[12]))
    eq13 = 1.0 - (x[9] + x[5])
    eq14 = 1.0 - (x[0] + x[15])
    eq15 = 1.0 - (x[13] + x[11])
    eq16 = 1.0 - (x[10] + x[8])
    return [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13, eq14, eq15, eq16]

# 约束条件
bounds = (np.zeros(16), np.array([1, np.inf, np.inf, np.inf, 1, np.inf, np.inf, 1, 1, 1, 1, 1, np.inf, 1, np.inf, 1]))

# 初始向量
x0 = np.zeros(16)
# x0 = np.array([1, 1e-4, 1, 1, 1, 1, 1e-4, 1, 1, 1, 1, 1, 1, 1, 1, 1])

res = least_squares(equations, x0, bounds=bounds)

print(res.x)

问题排查与解决建议 #

• 初始值优化：全零初始值x0 = np.zeros(16)大概率落在函数的平坦区域或约束边界，导致优化器无法找到有效下降方向。直接启用你注释掉的非零初始值，非线性优化对初始值敏感度极高，符合物理意义的初始值能大幅提升收敛概率。
• 简化固定变量：eq1和eq2的右侧是固定常量，可提前计算并直接赋值给x[4]和x[7]，减少优化变量数量，降低问题复杂度：

  # 提前计算固定值
  val7 = 1 / (np.exp((5.24677 - 1598.673 / (60.0 - -46.424)) * np.log(10.0)))
  val4 = 1 / (np.exp((5.08354 - 1663.125 / (60.0 - -45.622)) * np.log(10.0)))
  # 初始化时直接设置固定变量
  x0 = np.array([1, 1e-4, 1, 1, val4, 1, 1e-4, val7, 1, 1, 1, 1, 1, 1, 1, 1])
  

• 约束合理性检查：检查bounds中变量的上下限是否符合热力学物理意义，避免变量被卡在边界无法更新。
• 调整优化器参数：尝试切换优化方法（如method='dogbox'）或提高精度阈值，让优化器更敏感：

  res = least_squares(equations, x0, bounds=bounds, method='dogbox', ftol=1e-8, xtol=1e-8, gtol=1e-8)
  

• 残差验证：先打印初始值对应的残差print(equations(x0))，确认方程逻辑是否正确，针对性调整初始值或修正方程。

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