求解非线性方程：为吉布斯自由能问题添加约束

Question

我正在尝试使用拉格朗日方法和指数公式求解非线性系统，该系统将最小化吉布斯自由能。这些方程已经具有指数形式的拉格朗日函数Y1...Y6，稍后将其转换为化学物质的摩尔数n1...n9。问题是fsolve()给出的答案差异很大，即使用相同的猜测重新运行问题，它也会给出不同的值。但主要问题是我用不同的猜测得出的所有解决方案都没有物理意义，因为在将Ys转换为ns 后，我得到了负质量值。因此，通过所涉及的物理学，我可以确定所有[n1...n9] >= 0. 还可以确定 的所有最大值[n1...n9]。如何将其添加到代码中？import numpy as npimport scipyfrom scipy.optimize import fsolveimport time## "B" is the energy potentials of the species [C_gr , CO , CO2 , H2 , CH4 , H2O , N2* , SiO2* , H2S]B = [-11.0, -309.3632404425132, -613.3667287153355, -135.61840658777166, -269.52018727412405, -434.67499662354476, -193.0773646004259, -980.0, -230.02942769438977]# "a_atoms" is the number of atoms in the reactants [C, H, O, N*, S, SiO2*]  # * Elements that doesn't react. 'a_atoms = [4.27311296e-02, 8.10688756e-02, 6.17738749e-02, 1.32864225e-01, 3.18931655e-05, 3.74477901e-04]P_zero = 100.0 # Standard energy pressureP_eq = 95.0 # Reaction pressure# Standard temperature 298.15K, reaction temperature 940K.#start_time = time.time()def GibbsEq(z):# Lambda's exponentials:    Y1 = z[0]    Y2 = z[1]     Y3 = z[2]    Y4 = z[3]     Y5 = z[4]     Y6 = z[5]# Number of moles in each phase:    N1 = z[6]    N2 = z[7]    N3 = z[8]# Equations of energy conservation and mass conservation:    F = np.zeros(9)     F[0] = (P_zero/P_eq) * N1 * ((B[1] * (Y1 * Y3) + B[2] * (Y1 * Y3**2) + B[4] * (Y1 * Y2**2)) + N2 * (B[0] * Y1)) - a_atoms[0]    F[1] = (P_zero/P_eq) * N1 * (2 * B[3] * Y2**2 + 4 * B[4] * (Y1 * Y2**4) + 2 * B[5] * ((Y2**2) * Y3) + 2 * B[8] * ((Y2**2) * Y5)) - a_atoms[1]    F[2] = (P_zero/P_eq) * N1 * (B[1] * (Y1 * Y3) + 2 * B[2] * (Y1 * Y3**2) + B[5] * ((Y2**2) * Y3)) - a_atoms[2]    F[3] = (P_zero/P_eq) * N1 * (2 * B[6]**2) - a_atoms[3]    F[4] = (P_zero/P_eq) * N1 * (B[8] * ((Y2**2) * Y5)) - a_atoms[4]    F[5] = N3 * (B[7] * Y5)  - a_atoms[5]# 如何在 中添加约束fsolve？python 中是否有其他求解器具有更多约束选项以获得稳定性和初始猜测的更多独立性？谢谢！

Answer 1

两个问题的一个答案。

fsolve不支持约束。您可以将初始估计值提供为正值，但这并不能保证正根。但是，您可以将问题重新表述为优化问题，并使用任何优化函数（例如 ）最小化施加约束的成本函数 scipy.optimize.minimize。

作为一个最小的例子，如果你想找到方程 x*x -4 的正根，你可以做如下。

scipy.optimize.minimize(lambda x:(x*x-4)**2,x0= [5], bounds =((0,None),))

在bounds这需要（最小，最大）对参数可用于施于根正约束。

输出 ：

 fun: array([1.66882981e-17])

 hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>

      jac: array([1.27318954e-07])

  message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'

     nfev: 20

      nit: 9

   status: 0

  success: True

        x: array([2.])

按照这个，你的代码可以修改如下。只需添加边界，更改函数return语句，然后将fsolveto更改scipy.optimize.minimize为bounds。

import numpy as np

import scipy

from scipy.optimize import fsolve

import time

#

# "B" is the energy potentials of the species [C_gr , CO , CO2 , H2 , CH4 , H2O , N2* , SiO2* , H2S]

B = [-11.0, -309.3632404425132, -613.3667287153355, -135.61840658777166, -269.52018727412405, -434.67499662354476, -193.0773646004259, -980.0, -230.02942769438977]

# "a_atoms" is the number of atoms in the reactants [C, H, O, N*, S, SiO2*]  

# * Elements that doesn't react. '

a_atoms = [4.27311296e-02, 8.10688756e-02, 6.17738749e-02, 1.32864225e-01, 3.18931655e-05, 3.74477901e-04]

P_zero = 100.0 # Standard energy pressure

P_eq = 95.0 # Reaction pressure

# Standard temperature 298.15K, reaction temperature 940K.

#

start_time = time.time()

def GibbsEq(z):

# Lambda's exponentials:

    Y1 = z[0]

    Y2 = z[1] 

    Y3 = z[2]

    Y4 = z[3] 

    Y5 = z[4] 

    Y6 = z[5]

# Number of moles in each phase:

    N1 = z[6]

    N2 = z[7]

    N3 = z[8]

    bounds =((0,None),)*9

# Equations of energy conservation and mass conservation:

    F = np.zeros(9) 

    F[0] = (P_zero/P_eq) * N1 * ((B[1] * (Y1 * Y3) + B[2] * (Y1 * Y3**2) + B[4] * (Y1 * Y2**2)) + N2 * (B[0] * Y1)) - a_atoms[0]

    F[1] = (P_zero/P_eq) * N1 * (2 * B[3] * Y2**2 + 4 * B[4] * (Y1 * Y2**4) + 2 * B[5] * ((Y2**2) * Y3) + 2 * B[8] * ((Y2**2) * Y5)) - a_atoms[1]

    F[2] = (P_zero/P_eq) * N1 * (B[1] * (Y1 * Y3) + 2 * B[2] * (Y1 * Y3**2) + B[5] * ((Y2**2) * Y3)) - a_atoms[2]

    F[3] = (P_zero/P_eq) * N1 * (2 * B[6]**2) - a_atoms[3]

    F[4] = (P_zero/P_eq) * N1 * (B[8] * ((Y2**2) * Y5)) - a_atoms[4]

    F[5] = N3 * (B[7] * Y5)  - a_atoms[5]

# 

    F[6] = (P_zero/P_eq) * (B[1] * (Y1 * Y3) + B[2] * (Y1 * Y3**2) + B[3] * Y2**2 + B[4] * (Y1 * Y2**4) + B[5] * ((Y2**2) * Y3) + B[6] * Y4 + B[8] * Y5) - 1 

    F[7] = B[0] * Y1 - 1 

    F[8] = B[7] * Y6 - 1

    return (np.sum(F)**2)

#

zGuess = np.ones(9)

z = scipy.optimize.minimize(GibbsEq, zGuess , bounds=bounds)

end_time = time.time()

time_solution = (end_time - start_time)

print('Solving time: {} s'.format(time_solution))

#

print(z.x)

print(N_T)

for n in N_T:

    if n < 0:

        print('Error: there is negative values for mass in the solution!')

        break 

输出：

Solving time: 0.012451648712158203 s

[1.47559173 2.09905553 1.71722403 1.01828262 1.17529548 1.08815712

 1.00294916 1.00104157 1.08815763]