python - 求解非线性方程: add constraints to Gibbs free energy problem

Question

我正在尝试求解一个非线性系统，该系统将使用拉格朗日方法和指数公式来最小化吉布斯自由能。 方程中已包含指数形式 Y1...Y6 的拉格朗日量，随后将其转换为化学物质 n1...n9 的摩尔数。

问题是 fsolve()给出的答案差异很大，即使使用相同的猜测重新运行问题，也会给出不同的值。 但主要问题是，我通过不同猜测得出的所有解决方案都没有物理意义，因为在将 Y 转换为 n 后，我得到的质量为负值。

因此，根据所涉及的物理原理，我可以确定所有 [n1...n9] >= 0。还可以确定[n1...n9]的所有最大值。

如何将其添加到代码中？

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]
# 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 F
#
zGuess = np.ones(9)
z = scipy.optimize.fsolve(GibbsEq, zGuess)
end_time = time.time()
time_solution = (end_time - start_time)
print('Solving time: {} s'.format(time_solution))
#
n1 = z[7] * B[0] * z[0]
n2 = z[6] * B[1] * z[0] * z[2]
n3 = z[6] * B[2] * z[0] * z[2]**2
n4 = z[6] * B[3] * z[1]**2
n5 = z[6] * B[4] * z[0] * z[1]**4
n6 = z[6] * B[5] * z[1]**2 * z[4]
n7 = z[6] * B[6] * z[3]**2
n8 = z[8] * B[7] * z[5]
n9 = z[6] * B[8] * z[1]**2 * z[4]
N_T = [n1, n2, n3, n4, n5, n6, n7, n8, n9]
print(z)
print(z[6],z[7],z[8])
print(N_T)
for n in N_T:
    if n < 0:
        print('Error: there is negative values for mass in the solution!')
        break

1. 如何在 fsolve 中添加约束？
2. 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采用 (min,max) 对的参数可用于对根施加正约束。

输出:

 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语句，并更改 fsolve至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]