GEKKO optimizer @TypeError: 'GK_Operators' object is not callable
GEKKO optimizer @TypeError: 'GK_Operators' object is not callable
我正在做天气数据模拟，有四个方程和四个变量。还有很多其他参数，但提供了它们的值。我正在使用 GEKKO 优化器并收到此错误。我是编程新手。请帮忙。我已尽力提供尽可能多的信息。
# all libraries
from scipy.integrate import quad
import numpy as np
# all input parameters are as below
h = 1000 # height of air slab in meters
ps = 100000 # surface presuure in pascals
ph = 88000 # pressure at slab top in pascals
atop = 0.2 # entrainment parameter
b = 0.3 # moistening parameter
zh = 0.2 # hydrologically active soil depth in meters
zt = 0.4 # thermally active soil depth in m
n = 0.25 # soil porosity
vmin = 0.5 # Mineral volume fraction of soil
vorg = 0.25 # Organic volume fraction of soil
c = 1 # Exponent in evaporation efficiency
eta = 1 # Coefficient in runoff ratio
r = 2 # Exponent in runoff ratio
uz = 4 # Mixed layer wind speed in m/s
l = 500000 # Length scale of region
qin = 8/1000 #Effective humidity of incoming air
g = 9.81 # gravitational accelaration in m/s**2
rd = 287.058 # gas constant for dry air in SI uits
cp = 1005 # dry air specific heat at constant pressure in SI units
sigma = 5.670374419*10**(-8)
A = 0.75 # integration term
m = 1/7 # integration term
c1 = 0.001 # coefficient for forced velocity
c2 = 0.00025 # coefficient for buoyancy velocity
LHV = 22.5*10**5 # latent heat of vaporization of water in SI units
rv = 461 # gas consatant for water vapor in SI units
rho = 1.225 # air density in SI units
RS = 1000 # incoming solar radiation in SI units
yh = (ph/ps)**(3/2)
# calculation of incoming water vapor flux
Qin = ((ps-ph)/g)/l*uz*qin
def f1(y):
return (y**(8/3*rd/cp))*((y-yh)**(m-1))
int1 = quad(f1, yh, 1)
def f2(y):
return (y**(8/3*rd/cp))*((1-y)**(m-1))
int2 = quad(f2, yh, 1)
from gekko import GEKKO
m = GEKKO()
s = m.Var(value = 0.5)
qm = m.Var(value = 1)
tg = m.Var(value = 250)
thetam = m.Var(value = 350)
m.Equation(((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))-(((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))*eta*(s**r))-((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))==0)
m.Equation(Qin-(((ps-ph)/g)/l*uz*qm)+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))-((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))==0)
m.Equation(RS*(1-(0.17-0.085*s))+(1.24*((qm*ph)/(0.622*thetam*(ph/ps)**(rd/cp)))**(1/7)*sigma*(thetam*(ph/ps)**(rd/cp))**(4))(1-(0.75*((2/3*qm*ps/g*(1-(ph/ps)**(3/2)))**(1/7))))+(sigma*(thetam**(4))*m*A*((2/3*qm*ps/g)**(m))*int2)-(sigma*(tg**(4)))-((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*(tg-thetam)*rho*cp)-LHV*((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))==0)
m.Equation(((1.24*((qm*ph)/(0.622*thetam*(ph/ps)**(rd/cp)))**(1/7)*sigma*(thetam*(ph/ps)**(rd/cp))**(4))+(sigma*(tg**(4))))*0.75*((2/3*qm*ps/g*(1-(ph/ps)**(3/2)))**(1/7))-sigma*(thetam**(4))*m*A*((2/3*qm*ps/g)**(m))*int2-sigma*thetam**(4)*m*A*(2/3*qm*ps/g)**(m)*int1+1.2*(c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*(tg-thetam)*rho*cp==0)
m.solve()
错误指向模型中的特定方程。
m.Equation(RS*(1-(0.17-0.085*s))+(1.24*((qm*ph)/(0.622*thetam*\
(ph/ps)**(rd/cp)))**(1/7)*sigma*\
(thetam*(ph/ps)**(rd/cp))**(4))(1-(0.75*((2/3*qm*ps/g*\
(1-(ph/ps)**(3/2)))**(1/7))))+(sigma*(thetam**(4))*m*A*\
((2/3*qm*ps/g)**(m))*int2)-(sigma*(tg**(4)))-((c1*uz\
+ c2*((g/thetam*h*(thetam-tg))**0.5))*(tg-thetam)*rho*cp)\
-LHV*((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))\
*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))==0)
Traceback (most recent call last):
File "C:\Users\johnh\Desktop\test.py", line 59, in <module>
m.Equation(RS*(1-(0.17-0.085*s))+(1.24*((qm*ph)/(0.622*thetam*\
TypeError: 'GK_Operators' object is not callable
方程式至少缺少 1 个运算符，例如该方程式的第 3 行应该在中间 )( 括号之间包含一个缺失的运算符。
(thetam*(ph/ps)**(rd/cp))**(4))*(1-(0.75*((2/3*qm*ps/g*\
另一个问题是变量 m=1/7 被重新定义为具有 m=GEKKO() 的 gekko 模型。您可以为模型使用不同的名称，例如 gk=GEKKO() 以避免冲突。
排除模型故障并使其更具可读性的一个好方法是将模型的各个部分分解为子方程，例如使用带 gk.Intermediate() 的中间变量。另一种选择是在方程定义之外定义部分。我在下面展示了一个例子：
m.Equation(a+b*c+d**e)
简化版：
p1 = m.Intermediate(b*c) # define sub-equation with Intermediate
p2 = d**e # define sub-equation without Intermediate
m.Equation(a+p1+p2)
使用 m.Intermediate() 的优势在于，对于较大的模型，模型通常求解得更快。使用这种方法揭示了另一个问题。与 quad 的集成只需要 return 第一个元素。默认情况下，错误也是 returned.
int2 = quad(f2, yh, 1)[0]
完成所有这些修复后，模型现已完成：
# all libraries
from scipy.integrate import quad
import numpy as np
# all input parameters are as below
h = 1000 # height of air slab in meters
ps = 100000 # surface presuure in pascals
ph = 88000 # pressure at slab top in pascals
atop = 0.2 # entrainment parameter
b = 0.3 # moistening parameter
zh = 0.2 # hydrologically active soil depth in meters
zt = 0.4 # thermally active soil depth in m
n = 0.25 # soil porosity
vmin = 0.5 # Mineral volume fraction of soil
vorg = 0.25 # Organic volume fraction of soil
c = 1 # Exponent in evaporation efficiency
eta = 1 # Coefficient in runoff ratio
r = 2 # Exponent in runoff ratio
uz = 4 # Mixed layer wind speed in m/s
l = 500000 # Length scale of region
qin = 8/1000 #Effective humidity of incoming air
g = 9.81 # gravitational accelaration in m/s**2
rd = 287.058 # gas constant for dry air in SI uits
cp = 1005 # dry air specific heat at constant pressure in SI units
sigma = 5.670374419*10**(-8)
A = 0.75 # integration term
m = 1/7 # integration term
c1 = 0.001 # coefficient for forced velocity
c2 = 0.00025 # coefficient for buoyancy velocity
LHV = 22.5*10**5 # latent heat of vaporization of water in SI units
rv = 461 # gas consatant for water vapor in SI units
rho = 1.225 # air density in SI units
RS = 1000 # incoming solar radiation in SI units
yh = (ph/ps)**(3/2)
# calculation of incoming water vapor flux
Qin = ((ps-ph)/g)/l*uz*qin
def f1(y):
return (y**(8/3*rd/cp))*((y-yh)**(m-1))
int1 = quad(f1, yh, 1)[0]
def f2(y):
return (y**(8/3*rd/cp))*((1-y)**(m-1))
int2 = quad(f2, yh, 1)[0]
from gekko import GEKKO
gk = GEKKO(remote=False)
s = gk.Var(value = 0.5)
qm = gk.Var(value = 1)
tg = gk.Var(value = 250)
thetam = gk.Var(value = 350)
gk.Equation(((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))\
*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))\
-(((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))\
*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))*\
eta*(s**r))-((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))\
*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))==0)
gk.Equation(Qin-(((ps-ph)/g)/l*uz*qm)+((s**c)*((c1*uz + c2*((g/thetam*h*\
(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))\
-qm)*rho))-((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*\
(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)\
*LHV/rv)))-qm)*rho))))==0)
p1 = gk.Intermediate(1-(0.17-0.085*s))
p2 = gk.Intermediate(((qm*ph)/(0.622*thetam*(ph/ps)**(rd/cp))))
p3 = gk.Intermediate((1.24*p2**(1/7)*sigma*(thetam*(ph/ps)**(rd/cp))**(4)))
p4 = gk.Intermediate((1-(0.75*((2/3*qm*ps/g*(1-(ph/ps)**(3/2)))**(1/7)))))
p5 = gk.Intermediate(sigma*(thetam**(4))*m*A*((2/3*qm*ps/g)**(m))*int2)
p6 = gk.Intermediate((sigma*(tg**(4))))
p7 = gk.Intermediate((c1*uz+ c2*((g/thetam*h*(thetam-tg))**0.5)))
p8 = gk.Intermediate((p7*(tg-thetam)*rho*cp))
p9 = gk.Intermediate((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5)))
p10 = gk.Intermediate(((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm))
p11 = gk.Intermediate(LHV*((s**c)*(p9*p10*rho)))
gk.Equation(RS*p1+p3*p4+p5-p6-p8-p11==0)
gk.Equation(((1.24*((qm*ph)/(0.622*thetam*(ph/ps)**(rd/cp)))**(1/7)*\
sigma*(thetam*(ph/ps)**(rd/cp))**(4))+(sigma*(tg**(4))))*\
0.75*((2/3*qm*ps/g*(1-(ph/ps)**(3/2)))**(1/7))-sigma*(thetam**(4))\
*m*A*((2/3*qm*ps/g)**(m))*int2-sigma*thetam**(4)*m*A*(2/3*qm*ps/g)**(m)\
*int1+1.2*(c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*(tg-thetam)*rho*cp==0)
gk.options.SOLVER = 1
gk.open_folder()
gk.solve()
但是，求解器无法找到问题的解决方案：
----------------------------------------------
Steady State Optimization with APOPT Solver
----------------------------------------------
Iter Objective Convergence
0 2.99819E-17 3.88868E-01
1 4.00000E-20 1.35463E-01
2 4.00000E-20 7.81250E-03
3 4.00000E-20 0.00000E+00
No feasible solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 0.0291 sec
Objective : 0.
Unsuccessful with error code 0
---------------------------------------------------
Creating file: infeasibilities.txt
Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
@error: Solution Not Found
运行 文件夹以 gk.open_folder() 打开。 infeasibilities.txt 文件显示所有方程式的计算结果为 NaN。这意味着至少有一个变量发散到无穷大。我建议您检查模型，可能包括上限和下限，或其他故障排除方法。
我正在做天气数据模拟，有四个方程和四个变量。还有很多其他参数，但提供了它们的值。我正在使用 GEKKO 优化器并收到此错误。我是编程新手。请帮忙。我已尽力提供尽可能多的信息。
# all libraries
from scipy.integrate import quad
import numpy as np
# all input parameters are as below
h = 1000 # height of air slab in meters
ps = 100000 # surface presuure in pascals
ph = 88000 # pressure at slab top in pascals
atop = 0.2 # entrainment parameter
b = 0.3 # moistening parameter
zh = 0.2 # hydrologically active soil depth in meters
zt = 0.4 # thermally active soil depth in m
n = 0.25 # soil porosity
vmin = 0.5 # Mineral volume fraction of soil
vorg = 0.25 # Organic volume fraction of soil
c = 1 # Exponent in evaporation efficiency
eta = 1 # Coefficient in runoff ratio
r = 2 # Exponent in runoff ratio
uz = 4 # Mixed layer wind speed in m/s
l = 500000 # Length scale of region
qin = 8/1000 #Effective humidity of incoming air
g = 9.81 # gravitational accelaration in m/s**2
rd = 287.058 # gas constant for dry air in SI uits
cp = 1005 # dry air specific heat at constant pressure in SI units
sigma = 5.670374419*10**(-8)
A = 0.75 # integration term
m = 1/7 # integration term
c1 = 0.001 # coefficient for forced velocity
c2 = 0.00025 # coefficient for buoyancy velocity
LHV = 22.5*10**5 # latent heat of vaporization of water in SI units
rv = 461 # gas consatant for water vapor in SI units
rho = 1.225 # air density in SI units
RS = 1000 # incoming solar radiation in SI units
yh = (ph/ps)**(3/2)
# calculation of incoming water vapor flux
Qin = ((ps-ph)/g)/l*uz*qin
def f1(y):
return (y**(8/3*rd/cp))*((y-yh)**(m-1))
int1 = quad(f1, yh, 1)
def f2(y):
return (y**(8/3*rd/cp))*((1-y)**(m-1))
int2 = quad(f2, yh, 1)
from gekko import GEKKO
m = GEKKO()
s = m.Var(value = 0.5)
qm = m.Var(value = 1)
tg = m.Var(value = 250)
thetam = m.Var(value = 350)
m.Equation(((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))-(((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))*eta*(s**r))-((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))==0)
m.Equation(Qin-(((ps-ph)/g)/l*uz*qm)+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))-((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))==0)
m.Equation(RS*(1-(0.17-0.085*s))+(1.24*((qm*ph)/(0.622*thetam*(ph/ps)**(rd/cp)))**(1/7)*sigma*(thetam*(ph/ps)**(rd/cp))**(4))(1-(0.75*((2/3*qm*ps/g*(1-(ph/ps)**(3/2)))**(1/7))))+(sigma*(thetam**(4))*m*A*((2/3*qm*ps/g)**(m))*int2)-(sigma*(tg**(4)))-((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*(tg-thetam)*rho*cp)-LHV*((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))==0)
m.Equation(((1.24*((qm*ph)/(0.622*thetam*(ph/ps)**(rd/cp)))**(1/7)*sigma*(thetam*(ph/ps)**(rd/cp))**(4))+(sigma*(tg**(4))))*0.75*((2/3*qm*ps/g*(1-(ph/ps)**(3/2)))**(1/7))-sigma*(thetam**(4))*m*A*((2/3*qm*ps/g)**(m))*int2-sigma*thetam**(4)*m*A*(2/3*qm*ps/g)**(m)*int1+1.2*(c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*(tg-thetam)*rho*cp==0)
m.solve()
错误指向模型中的特定方程。
m.Equation(RS*(1-(0.17-0.085*s))+(1.24*((qm*ph)/(0.622*thetam*\
(ph/ps)**(rd/cp)))**(1/7)*sigma*\
(thetam*(ph/ps)**(rd/cp))**(4))(1-(0.75*((2/3*qm*ps/g*\
(1-(ph/ps)**(3/2)))**(1/7))))+(sigma*(thetam**(4))*m*A*\
((2/3*qm*ps/g)**(m))*int2)-(sigma*(tg**(4)))-((c1*uz\
+ c2*((g/thetam*h*(thetam-tg))**0.5))*(tg-thetam)*rho*cp)\
-LHV*((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))\
*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))==0)
Traceback (most recent call last):
File "C:\Users\johnh\Desktop\test.py", line 59, in <module>
m.Equation(RS*(1-(0.17-0.085*s))+(1.24*((qm*ph)/(0.622*thetam*\
TypeError: 'GK_Operators' object is not callable
方程式至少缺少 1 个运算符，例如该方程式的第 3 行应该在中间 )( 括号之间包含一个缺失的运算符。
(thetam*(ph/ps)**(rd/cp))**(4))*(1-(0.75*((2/3*qm*ps/g*\
另一个问题是变量 m=1/7 被重新定义为具有 m=GEKKO() 的 gekko 模型。您可以为模型使用不同的名称，例如 gk=GEKKO() 以避免冲突。
排除模型故障并使其更具可读性的一个好方法是将模型的各个部分分解为子方程，例如使用带 gk.Intermediate() 的中间变量。另一种选择是在方程定义之外定义部分。我在下面展示了一个例子：
m.Equation(a+b*c+d**e)
简化版：
p1 = m.Intermediate(b*c) # define sub-equation with Intermediate
p2 = d**e # define sub-equation without Intermediate
m.Equation(a+p1+p2)
使用 m.Intermediate() 的优势在于，对于较大的模型，模型通常求解得更快。使用这种方法揭示了另一个问题。与 quad 的集成只需要 return 第一个元素。默认情况下，错误也是 returned.
int2 = quad(f2, yh, 1)[0]
完成所有这些修复后，模型现已完成：
# all libraries
from scipy.integrate import quad
import numpy as np
# all input parameters are as below
h = 1000 # height of air slab in meters
ps = 100000 # surface presuure in pascals
ph = 88000 # pressure at slab top in pascals
atop = 0.2 # entrainment parameter
b = 0.3 # moistening parameter
zh = 0.2 # hydrologically active soil depth in meters
zt = 0.4 # thermally active soil depth in m
n = 0.25 # soil porosity
vmin = 0.5 # Mineral volume fraction of soil
vorg = 0.25 # Organic volume fraction of soil
c = 1 # Exponent in evaporation efficiency
eta = 1 # Coefficient in runoff ratio
r = 2 # Exponent in runoff ratio
uz = 4 # Mixed layer wind speed in m/s
l = 500000 # Length scale of region
qin = 8/1000 #Effective humidity of incoming air
g = 9.81 # gravitational accelaration in m/s**2
rd = 287.058 # gas constant for dry air in SI uits
cp = 1005 # dry air specific heat at constant pressure in SI units
sigma = 5.670374419*10**(-8)
A = 0.75 # integration term
m = 1/7 # integration term
c1 = 0.001 # coefficient for forced velocity
c2 = 0.00025 # coefficient for buoyancy velocity
LHV = 22.5*10**5 # latent heat of vaporization of water in SI units
rv = 461 # gas consatant for water vapor in SI units
rho = 1.225 # air density in SI units
RS = 1000 # incoming solar radiation in SI units
yh = (ph/ps)**(3/2)
# calculation of incoming water vapor flux
Qin = ((ps-ph)/g)/l*uz*qin
def f1(y):
return (y**(8/3*rd/cp))*((y-yh)**(m-1))
int1 = quad(f1, yh, 1)[0]
def f2(y):
return (y**(8/3*rd/cp))*((1-y)**(m-1))
int2 = quad(f2, yh, 1)[0]
from gekko import GEKKO
gk = GEKKO(remote=False)
s = gk.Var(value = 0.5)
qm = gk.Var(value = 1)
tg = gk.Var(value = 250)
thetam = gk.Var(value = 350)
gk.Equation(((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))\
*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))\
-(((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))\
*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))*\
eta*(s**r))-((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))\
*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))==0)
gk.Equation(Qin-(((ps-ph)/g)/l*uz*qm)+((s**c)*((c1*uz + c2*((g/thetam*h*\
(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))\
-qm)*rho))-((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*\
(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)\
*LHV/rv)))-qm)*rho))))==0)
p1 = gk.Intermediate(1-(0.17-0.085*s))
p2 = gk.Intermediate(((qm*ph)/(0.622*thetam*(ph/ps)**(rd/cp))))
p3 = gk.Intermediate((1.24*p2**(1/7)*sigma*(thetam*(ph/ps)**(rd/cp))**(4)))
p4 = gk.Intermediate((1-(0.75*((2/3*qm*ps/g*(1-(ph/ps)**(3/2)))**(1/7)))))
p5 = gk.Intermediate(sigma*(thetam**(4))*m*A*((2/3*qm*ps/g)**(m))*int2)
p6 = gk.Intermediate((sigma*(tg**(4))))
p7 = gk.Intermediate((c1*uz+ c2*((g/thetam*h*(thetam-tg))**0.5)))
p8 = gk.Intermediate((p7*(tg-thetam)*rho*cp))
p9 = gk.Intermediate((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5)))
p10 = gk.Intermediate(((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm))
p11 = gk.Intermediate(LHV*((s**c)*(p9*p10*rho)))
gk.Equation(RS*p1+p3*p4+p5-p6-p8-p11==0)
gk.Equation(((1.24*((qm*ph)/(0.622*thetam*(ph/ps)**(rd/cp)))**(1/7)*\
sigma*(thetam*(ph/ps)**(rd/cp))**(4))+(sigma*(tg**(4))))*\
0.75*((2/3*qm*ps/g*(1-(ph/ps)**(3/2)))**(1/7))-sigma*(thetam**(4))\
*m*A*((2/3*qm*ps/g)**(m))*int2-sigma*thetam**(4)*m*A*(2/3*qm*ps/g)**(m)\
*int1+1.2*(c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*(tg-thetam)*rho*cp==0)
gk.options.SOLVER = 1
gk.open_folder()
gk.solve()
但是，求解器无法找到问题的解决方案：
----------------------------------------------
Steady State Optimization with APOPT Solver
----------------------------------------------
Iter Objective Convergence
0 2.99819E-17 3.88868E-01
1 4.00000E-20 1.35463E-01
2 4.00000E-20 7.81250E-03
3 4.00000E-20 0.00000E+00
No feasible solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 0.0291 sec
Objective : 0.
Unsuccessful with error code 0
---------------------------------------------------
Creating file: infeasibilities.txt
Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
@error: Solution Not Found
运行 文件夹以 gk.open_folder() 打开。 infeasibilities.txt 文件显示所有方程式的计算结果为 NaN。这意味着至少有一个变量发散到无穷大。我建议您检查模型，可能包括上限和下限，或其他故障排除方法。