gpt4 book ai didi

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

转载 作者:行者123 更新时间:2023-12-01 08:04:31 29 4
gpt4 key购买 nike

我正在尝试求解一个非线性系统,该系统将使用拉格朗日方法和指数公式来最小化吉布斯自由能。方程中已包含指数形式 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 中是否还有其他求解器具有更多约束选项,以获得稳定性和初始猜测的更多独立性?

谢谢!

最佳答案

两个问题都有一个答案。

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语句,并更改 fsolvescipy.optimize.minimizebounds

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]

关于python - 求解非线性方程: add constraints to Gibbs free energy problem,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/55599588/

29 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com