python - 最小化温度向量中每个 T 值的多变量函数

Question

我已尝试根据社区中找到的答案使用 GEKKO 构建代码，但无法解决我的问题。它是一个函数 G(T)，应该求解向量 T 的每个分量。是不是m.Obj部分的错误？

错误及代码如下:

      File "C:/Users/Roberto/PycharmProjects/Testes iniciais/Dry reforming - isothermic.py", line 119, in <module>
        m.solve()
      File "C:\Users\Roberto\PycharmProjects\Testes iniciais\venv\lib\site-packages\gekko\gekko.py", line 1882, in solve
        raise Exception(response)
    Exception:  @error: Equation Definition
     Warning: there is insufficient data in CSV file 1 
   51.21.105.203_gk_model0.csv
     @error: Equation Definition
     Equation without an equality (=) or inequality (>,<)
     true
     Equation without an equality (=) or inequality (>,<)
     STOPPING...
     true
     STOPPING...

进程已完成，退出代码为 1。

我只需要知道错误消息的来源。

# Equilibrium calculations - Isothermal conditions

# library with mathematical package

import math
import numpy as np

# ############   Substances ID:       
##########################################
# A=methane           B=oxygen
# C=hydrogen          D=carbon monoxide
# E=carbon dioxide    F=water
# G=carbon            H=Nitrogen
# I=argon

# ############ Thermodynamic data             
##########################################
T0 = 298.15 #[K]
R = 8.134 #[J/mol.K]
p = 1 #[atm]

# A) Standard Enthalpy (hi0) [J/mol] and Entropy (si0) [J/mol] at T0=298.15K
# Source: [Perry's Chemical, Engineers' Handbook seventh edition, Mc Graw-Hill, 1997]

hA0 = -74520;      sA0 = 1.8627*10**2
hB0 = 0;           sB0 = 2.05149*10**2
hC0 = 0;           sC0 = 1.3057*10**2
hD0 = -110530;     sD0 = 1.9756*10**2
hE0 = -393510;     sE0 = 2.1368*10**2
hF0 = -241814;     sF0 = 1.8872*10**2
hG0 = 0;           sG0 = 5.734
hH0 = 0;           sH0 = 1.91610*10**2
#hI0 = 0;           sI0 = 1.5484*10**2

# B) Ideal gas heat capacity polynomial: Coefficients for Cp in [J/mol.K]
# Source: [Chemical thermodynamics for process simulation, Appendix G, pg:711]

C1A = 19.238;   C2A = 52.090*10**-3;   C3A = 11.966*10**-6;   C4A = -11.309*10**-9
C1B = 28.087;   C2B = -0.004*10**-3;   C3B = 17.447*10**-6;   C4B = -10.644*10**-9
C1C = 27.124;   C2C = 9.2670*10**-3;   C3C = -13.799*10**-6;  C4C = 7.64000*10**-9
C1D = 30.848;   C2D = -12.84*10**-3;   C3D = 27.870*10**-6;   C4D = -12.710*10**-9
C1E = 19.780;   C2E = 73.390*10**-3;   C3E = -55.98*10**-6;   C4E = 17.1400*10**-9
C1F = 32.220;   C2F = 1.9225*10**-3;   C3F = 10.548*10**-6;   C4F = -3.5940*10**-9
C1G = -5.416;   C2G = 58.981*10**-3;   C3G = -43.559*10**-6;  C4G = 11.604*10**-9
C1H = 31.128;   C2H = -13.556*10**-3;  C3H = 26.777*10**-6;   C4H = -11.673*10**-9

# ####### Thermodynamic Calculations     #####################################
# A) Standard Enthalpy (hi0), Entropy (si0) and Gibbs free molar energy (mi) at T0 [J/mol]

T = np.arange(50, 1000, 50)
for T in T:
hA = hA0 + C1A*(T - T0) + (C2A/2) * (T**2 - T0**2) + (C3A/3) * (T**3 - T0**3) + (C4A/4) * (T**4 - T0**4)
hB = hB0 + C1B*(T - T0) + (C2B/2) * (T**2 - T0**2) + (C3B/3) * (T**3 - T0**3) + (C4B/4) * (T**4 - T0**4)
hC = hC0 + C1C*(T - T0) + (C2C/2) * (T**2 - T0**2) + (C3C/3) * (T**3 - T0**3) + (C4C/4) * (T**4 - T0**4)
hD = hD0 + C1D*(T - T0) + (C2D/2) * (T**2 - T0**2) + (C3D/3) * (T**3 - T0**3) + (C4D/4) * (T**4 - T0**4)
hE = hE0 + C1E*(T - T0) + (C2E/2) * (T**2 - T0**2) + (C3E/3) * (T**3 - T0**3) + (C4E/4) * (T**4 - T0**4)
hF = hF0 + C1F*(T - T0) + (C2F/2) * (T**2 - T0**2) + (C3F/3) * (T**3 - T0**3) + (C4F/4) * (T**4 - T0**4)
hG = hG0 + C1G*(T - T0) + (C2G/2) * (T**2 - T0**2) + (C3G/3) * (T**3 - T0**3) + (C4G/4) * (T**4 - T0**4)
 hH = hH0 + C1H*(T - T0) + (C2H/2) * (T**2 - T0**2) + (C3H/3) * (T**3 - T0**3) + (C4H/4) * (T**4 - T0**4)

sA = sA0 + 1/6 * (6 * C2A * T + 3 * C3A * T**2 + 2 * C4A * T**3 - 6 * C2A * T0 - 3 * C3A * T0**2 - 2 * C4A * T0**3 +6 * C1A * (np.log(T) - np.log(T0)))
sB = sB0 + 1/6 * (6 * C2B * T + 3 * C3B * T**2 + 2 * C4B * T**3 - 6 * C2B * T0 - 3 * C3B * T0**2 - 2 * C4B * T0**3 +6 * C1B * (np.log(T) - np.log(T0)))
sC = sC0 + 1/6 * (6 * C2C * T + 3 * C3C * T**2 + 2 * C4C * T**3 - 6 * C2C * T0 - 3 * C3C * T0**2 - 2 * C4C * T0**3 +6 * C1C * (np.log(T) - np.log(T0)))
sD = sD0 + 1/6 * (6 * C2D * T + 3 * C3D * T**2 + 2 * C4D * T**3 - 6 * C2D * T0 - 3 * C3D * T0**2 - 2 * C4D * T0**3 +6 * C1D * (np.log(T) - np.log(T0)))
sE = sE0 + 1/6 * (6 * C2E * T + 3 * C3E * T**2 + 2 * C4E * T**3 - 6 * C2E * T0 - 3 * C3E * T0**2 - 2 * C4E * T0**3 +6 * C1E * (np.log(T) - np.log(T0)))
sF = sF0 + 1/6 * (6 * C2F * T + 3 * C3F * T**2 + 2 * C4F * T**3 - 6 * C2F * T0 - 3 * C3F * T0**2 - 2 * C4F * T0**3 +6 * C1F * (np.log(T) - np.log(T0)))
sG = sG0 + 1/6 * (6 * C2G * T + 3 * C3G * T**2 + 2 * C4G * T**3 - 6 * C2G * T0 - 3 * C3G * T0**2 - 2 * C4G * T0**3 +6 * C1G * (np.log(T) - np.log(T0)))
sH = sH0 + 1/6 * (6 * C2H * T + 3 * C3H * T**2 + 2 * C4H * T**3 - 6 * C2H * T0 - 3 * C3H * T0**2 - 2 * C4H * T0**3 +6 * C1H * (np.log(T) - np.log(T0)))

mA = hA-T*sA
mB = hB-T*sB
mC = hC-T*sC
mD = hD-T*sD
mE = hE-T*sE
mF = hF-T*sF
mG = hG-T*sG
mH = hH-T*sH

# ################################# Gibbs free energy minimization routine ################################################################

# Inform from a to b the respective molar inlet quantity:
a = 1; # CH4
b = 1; # O2
c = 0; # H2
d = 0; # CO
e = 0; # CO2
f = 0; # H2O
g = 0; # C
h = 1-(a+b+c+d+e+f+g); # N2
#i = 0; # Ar

from gekko import GEKKO
m = GEKKO()

# Variables to be minimized:
nA, nB, nC, nD, nE, nF, nG, nH = [m.Var() for i in range(8)]

# Initial values:
nA = 0.2; nB = 0.2; nC = 0.2; nD = 0.1; nE = 0.1; nF = 0.1; nG = 0.05; nH = 2
nt = nA + nB + nC + nD + nE + nF + nG + nH

# Boundary conditions (in this case, boundary conditions are the atomic balances for H, O and C)

m.Equation(nA>=0); m.Equation(nB>=0); m.Equation(nC>=0); m.Equation(nD>=0);         m.Equation(nE>=0); m.Equation(nF>=0); m.Equation(nG>=0); m.Equation(nH>=0)
m.Equation(nA+nD+nE+nG==a)
m.Equation(2*nB+nD+2*nE+nF==2*b)
m.Equation(4*nA+2*nC+2*nF==4*a)
m.Equation(nA + nB + nC + nD + nE + nF + nG + nH == nt)

#Gibbs free energy function to be minimized (here for gekko it should be the "Objective")

# Objective:
m.Obj(nA*mA + nB*mB + nC*mC + nD*mD + nE*mE + nF*mF + nG*mG + nH*mH + \
  R*T*(nA*np.log(((nA*p)/nt)) + nB*np.log(((nB*p)/nt))+     nC*np.log(((nC*p)/nt)) + \
       nD*np.log(((nD*p)/nt)) + nE*np.log(((nE*p)/nt)) +     nF*np.log(((nF*p)/nt)) + nH*np.log(((nH*p)/nt))))
# Set global options
m.options.IMODE = 3

# Solve minimization
m.solve()

# Results
print('')
print('Results')
print('nA: ' + str(nA.value))

我该如何解决这个问题？

Answer 1

最小化吉布斯自由能的好应用!一些需要修复的东西:

• 在定义 Gekko 方程时使用 m.log 而不是 np.log。这允许自动微分为求解器 (IPOPT) 提供精确的一阶和二阶导数。
• 不要用数值覆盖 nA-nH。使用 nA.value 属性。
• 您可以添加诸如 nA>=0 之类的方程式，或者通过设置 nA.lower=0 来添加下限。对变量使用上限和下限比添加不等式约束更有效。
• 您发布的代码中的缩进似乎不正确。我假设 for T in T: 语句之后的所有内容都应该作为该循环的一部分缩进。请检查此项。
• 我为您的变量添加了 5.0 的上限。否则，解是无界的。对于吉布斯自由能最小化问题，您通常不应该在边界处设置变量，因此也请检查一下。我还将下限设置为 0.01，这样 m.log 项就不会为零并导致计算错误。

# Equilibrium calculations - Isothermal conditions

# library with mathematical package

import math
import numpy as np

# ############   Substances ID:       
##########################################
# A=methane           B=oxygen
# C=hydrogen          D=carbon monoxide
# E=carbon dioxide    F=water
# G=carbon            H=Nitrogen
# I=argon

# ############ Thermodynamic data             
##########################################
T0 = 298.15 #[K]
R = 8.134 #[J/mol.K]
p = 1 #[atm]

# A) Standard Enthalpy (hi0) [J/mol] and Entropy (si0) [J/mol] at T0=298.15K
# Source: [Perry's Chemical, Engineers' Handbook seventh edition, Mc Graw-Hill, 1997]

hA0 = -74520;      sA0 = 1.8627*10**2
hB0 = 0;           sB0 = 2.05149*10**2
hC0 = 0;           sC0 = 1.3057*10**2
hD0 = -110530;     sD0 = 1.9756*10**2
hE0 = -393510;     sE0 = 2.1368*10**2
hF0 = -241814;     sF0 = 1.8872*10**2
hG0 = 0;           sG0 = 5.734
hH0 = 0;           sH0 = 1.91610*10**2
#hI0 = 0;           sI0 = 1.5484*10**2

# B) Ideal gas heat capacity polynomial: Coefficients for Cp in [J/mol.K]
# Source: [Chemical thermodynamics for process simulation, Appendix G, pg:711]

C1A = 19.238;   C2A = 52.090*10**-3;   C3A = 11.966*10**-6;   C4A = -11.309*10**-9
C1B = 28.087;   C2B = -0.004*10**-3;   C3B = 17.447*10**-6;   C4B = -10.644*10**-9
C1C = 27.124;   C2C = 9.2670*10**-3;   C3C = -13.799*10**-6;  C4C = 7.64000*10**-9
C1D = 30.848;   C2D = -12.84*10**-3;   C3D = 27.870*10**-6;   C4D = -12.710*10**-9
C1E = 19.780;   C2E = 73.390*10**-3;   C3E = -55.98*10**-6;   C4E = 17.1400*10**-9
C1F = 32.220;   C2F = 1.9225*10**-3;   C3F = 10.548*10**-6;   C4F = -3.5940*10**-9
C1G = -5.416;   C2G = 58.981*10**-3;   C3G = -43.559*10**-6;  C4G = 11.604*10**-9
C1H = 31.128;   C2H = -13.556*10**-3;  C3H = 26.777*10**-6;   C4H = -11.673*10**-9

# ####### Thermodynamic Calculations     #####################################
# A) Standard Enthalpy (hi0), Entropy (si0) and Gibbs free molar energy (mi) at T0 [J/mol]

T = np.arange(50, 1000, 50)
for T in T:
   hA = hA0 + C1A*(T - T0) + (C2A/2) * (T**2 - T0**2) + (C3A/3) * (T**3 - T0**3) + (C4A/4) * (T**4 - T0**4)
   hB = hB0 + C1B*(T - T0) + (C2B/2) * (T**2 - T0**2) + (C3B/3) * (T**3 - T0**3) + (C4B/4) * (T**4 - T0**4)
   hC = hC0 + C1C*(T - T0) + (C2C/2) * (T**2 - T0**2) + (C3C/3) * (T**3 - T0**3) + (C4C/4) * (T**4 - T0**4)
   hD = hD0 + C1D*(T - T0) + (C2D/2) * (T**2 - T0**2) + (C3D/3) * (T**3 - T0**3) + (C4D/4) * (T**4 - T0**4)
   hE = hE0 + C1E*(T - T0) + (C2E/2) * (T**2 - T0**2) + (C3E/3) * (T**3 - T0**3) + (C4E/4) * (T**4 - T0**4)
   hF = hF0 + C1F*(T - T0) + (C2F/2) * (T**2 - T0**2) + (C3F/3) * (T**3 - T0**3) + (C4F/4) * (T**4 - T0**4)
   hG = hG0 + C1G*(T - T0) + (C2G/2) * (T**2 - T0**2) + (C3G/3) * (T**3 - T0**3) + (C4G/4) * (T**4 - T0**4)
   hH = hH0 + C1H*(T - T0) + (C2H/2) * (T**2 - T0**2) + (C3H/3) * (T**3 - T0**3) + (C4H/4) * (T**4 - T0**4)

   sA = sA0 + 1/6 * (6 * C2A * T + 3 * C3A * T**2 + 2 * C4A * T**3 - 6 * C2A * T0 - 3 * C3A * T0**2 - 2 * C4A * T0**3 +6 * C1A * (np.log(T) - np.log(T0)))
   sB = sB0 + 1/6 * (6 * C2B * T + 3 * C3B * T**2 + 2 * C4B * T**3 - 6 * C2B * T0 - 3 * C3B * T0**2 - 2 * C4B * T0**3 +6 * C1B * (np.log(T) - np.log(T0)))
   sC = sC0 + 1/6 * (6 * C2C * T + 3 * C3C * T**2 + 2 * C4C * T**3 - 6 * C2C * T0 - 3 * C3C * T0**2 - 2 * C4C * T0**3 +6 * C1C * (np.log(T) - np.log(T0)))
   sD = sD0 + 1/6 * (6 * C2D * T + 3 * C3D * T**2 + 2 * C4D * T**3 - 6 * C2D * T0 - 3 * C3D * T0**2 - 2 * C4D * T0**3 +6 * C1D * (np.log(T) - np.log(T0)))
   sE = sE0 + 1/6 * (6 * C2E * T + 3 * C3E * T**2 + 2 * C4E * T**3 - 6 * C2E * T0 - 3 * C3E * T0**2 - 2 * C4E * T0**3 +6 * C1E * (np.log(T) - np.log(T0)))
   sF = sF0 + 1/6 * (6 * C2F * T + 3 * C3F * T**2 + 2 * C4F * T**3 - 6 * C2F * T0 - 3 * C3F * T0**2 - 2 * C4F * T0**3 +6 * C1F * (np.log(T) - np.log(T0)))
   sG = sG0 + 1/6 * (6 * C2G * T + 3 * C3G * T**2 + 2 * C4G * T**3 - 6 * C2G * T0 - 3 * C3G * T0**2 - 2 * C4G * T0**3 +6 * C1G * (np.log(T) - np.log(T0)))
   sH = sH0 + 1/6 * (6 * C2H * T + 3 * C3H * T**2 + 2 * C4H * T**3 - 6 * C2H * T0 - 3 * C3H * T0**2 - 2 * C4H * T0**3 +6 * C1H * (np.log(T) - np.log(T0)))

   mA = hA-T*sA
   mB = hB-T*sB
   mC = hC-T*sC
   mD = hD-T*sD
   mE = hE-T*sE
   mF = hF-T*sF
   mG = hG-T*sG
   mH = hH-T*sH

   # ################################# Gibbs free energy minimization routine ################################################################

   # Inform from a to b the respective molar inlet quantity:
   a = 1; # CH4
   b = 1; # O2
   c = 0; # H2
   d = 0; # CO
   e = 0; # CO2
   f = 0; # H2O
   g = 0; # C
   h = 1-(a+b+c+d+e+f+g); # N2
   #i = 0; # Ar

   from gekko import GEKKO
   m = GEKKO()

   # Variables to be minimized:
   nA, nB, nC, nD, nE, nF, nG, nH = [m.Var() for i in range(8)]
   var = [nA, nB, nC, nD, nE, nF, nG, nH]

   # Initial values:
   n0 = [0.2,0.2,0.2,0.1,0.1,0.1,0.05,2]
   nL = np.ones(len(n0))*0.01
   nU = np.ones(len(n0))*5.0
   for i,x in enumerate(var):
      x.value = n0[i]
      x.lower = nL[i]
      x.upper = nU[i]

   nt = m.Intermediate(nA + nB + nC + nD + nE + nF + nG + nH)

   # Boundary conditions (in this case, boundary conditions are the atomic balances for H, O and C)

   #m.Equation(nA>=0); m.Equation(nB>=0); m.Equation(nC>=0); m.Equation(nD>=0)
   #m.Equation(nE>=0); m.Equation(nF>=0); m.Equation(nG>=0); m.Equation(nH>=0)
   m.Equation(nA+nD+nE+nG==a)
   m.Equation(2*nB+nD+2*nE+nF==2*b)
   m.Equation(4*nA+2*nC+2*nF==4*a)
   m.Equation(nA + nB + nC + nD + nE + nF + nG + nH == nt)

   #Gibbs free energy function to be minimized (here for gekko it should be the "Objective")

   # Objective:
   m.Obj(nA*mA + nB*mB + nC*mC + nD*mD + nE*mE + nF*mF + nG*mG + nH*mH + \
     R*T*(nA*m.log(((nA*p)/nt)) + nB*m.log(((nB*p)/nt))+ nC*m.log(((nC*p)/nt)) + \
          nD*m.log(((nD*p)/nt)) + nE*m.log(((nE*p)/nt)) + nF*m.log(((nF*p)/nt)) + \
          nH*m.log(((nH*p)/nt))))
   # Set global options
   m.options.IMODE = 3

   # Solve minimization
   m.solve()

   # Results
   print('')
   print('Results')
   print('nA: ' + str(nA.value))