问题描述
我正在做天气数据模拟并且有四个变量的四个方程。还有很多其他参数，但提供了它们的值。我正在使用 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,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() 的壁虎模型。您可以为模型使用不同的名称，例如 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 的集成只需要返回第一个元素。默认情况下也会返回错误。
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,1)[0]
def f2(y):
return (y**(8/3*rd/cp))*((1-y)**(m-1))
int2 = quad(f2,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。这意味着至少有一个变量发散到无穷大。我建议您检查模型，可能包括上限和下限，或其他故障排除方法。