Решение нелинейных уравнений: добавьте ограничения к проблеме свободной энергии Гиббса - PullRequest
Новая
       40

Решение нелинейных уравнений: добавьте ограничения к проблеме свободной энергии Гиббса

1 голос
/ 09 апреля 2019

Я пытаюсь решить нелинейную систему, которая минимизирует свободную энергию Гиббса, используя метод Лагранжа, с экспоненциальной формулировкой.В уравнениях уже есть лагранжианы в экспоненциальной форме Y1...Y6, которые впоследствии преобразуются в число молей химических частиц n1...n9.

Проблема в том, что fsolve() дает ответыкоторые сильно различаются, даже если повторять проблему с одним и тем же предположением, это дает разные значения.Но главная проблема заключается в том, что все решения, которые я получаю с разными догадками, не имеют физического смысла, потому что после преобразования Y s в n s я получаю отрицательные значения для массы.

Итак, по физикеЯ могу определить, что все [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 другие решатели, которые имеют больше опций ограничений для получения стабильности и большей независимости от начального предположения?

Спасибо!

1 Ответ

1 голос
/ 10 апреля 2019

Один ответ на оба вопроса.

fsolve не поддерживает ограничения.Вы можете предоставить начальную оценку как положительные значения, но это не гарантирует положительные корни.Однако вы можете переформулировать вашу проблему как задачу оптимизации и минимизировать функцию затрат, накладывающую ограничения, используя любую функцию оптимизации, такую ​​как scipy.optimize.minimize.

В качестве минимального примера, если вы хотите найти положительный корень уравнения x* х-4, вы можете сделать, как показано ниже.

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

Параметр bounds, который принимает (мин, макс) пару, можно использовать для наложения положительного ограничения на корень.

output: 

 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]
...