问题描述
我在 Gekko 中编写了一个模型并遇到了错误。于是,我开始逐行检查,发现了一个我无法理解的问题。下面的代码(显示部分模型代码,输入文件link)可以执行,结果也正确。但是,当我取消注释第 31 行时,发生错误。是不是因为我按顺序使用了逻辑条件 (ab3
,max3
)?有什么办法可以解决这个问题吗?
import numpy as np
from gekko import GEKKO
# Read weather data
weatherData = np.load("Atlanta_TMY3_climate.npy")
# create GEKKO model
m = GEKKO(remote=False)
m.open_folder()
print(m.path)
weatherData = np.insert(weatherData,np.zeros((1,15)),0)
days_to_consider = 1
m.time = np.linspace(0,24*days_to_consider,24*days_to_consider+1)
# Declare parameters
focc = m.Param(value = [0,0.05,0.1,\
0.1,0.1],name="focc")
setpoint = [16.00,16.00,\
16.00,16.38,17.43,17.88,18.08,18.35,\
17.85,17.13,16.27,16.00]
T_air_set = m.Param(value=setpoint,name="T_air_set")
Te = m.Param(value = weatherData[0:24*days_to_consider+1,0],name="Te")
qv_infil_wind = m.Param(value = 0.0769*0.6*(0.75*0.8* \
weatherData[0:24*days_to_consider+1,1]**2)**0.667,\
name="qv_infil_wind")
# 1) Airflow
qv_mech_sup = m.Intermediate(2.519244226731981*focc*1,name="qv_mech_sup")
qv_infil_stack = m.Intermediate(0.0146*0.6*(0.7*8.5* \
m.abs3(Te-T_air_set))**0.667,\
name="qv_infil_stack")
# When below line code is uncommented,an error occurs.
##qv_infil_sw = m.Intermediate(m.max3(qv_infil_stack,qv_infil_wind) + \
## 0.14*qv_infil_stack*qv_infil_wind/0.6,\
## name="qv_infil_sw")
m.options.IMODE = 6
m.options.DIAGLEVEL = 4
m.options.soLVER = 3
m.options.MAX_ITER = 1000
m.solve(disp=True,GUI=False)
解决方法
我改用 m.options.SOLVER=1
是因为 m.abs3()
和 m.max3()
函数需要混合整数非线性规划 (MINLP) 求解器。然而,这也没有解决问题。似乎每个二进制变量都有难以解决的组合。我不确定为什么,但这是一个求解器问题,而不是 Gekko 公式问题。这是一个简单的问题,它表明 m.abs3()
和 m.max3()
在同一模型中并作为依赖项不是问题。
from gekko import GEKKO
m = GEKKO()
x = m.Param(2)
y = m.Param(3)
z = m.abs3(x-y)
w = m.max3(z,x)
m.solve()
print(z.value[0],w.value[0])
切换到 m.abs()
有一个成功的解决方案。这避免了求解器收敛问题。我通常不建议使用 m.abs()
,但所有因变量都是参数,因此非连续梯度不存在问题。
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
# Read weather data
weatherData = np.load("Atlanta_TMY3_climate.npy")
# create GEKKO model
m = GEKKO(remote=False)
print(m.path)
weatherData = np.insert(weatherData,np.zeros((1,15)),0)
days_to_consider = 1
m.time = np.linspace(0,24*days_to_consider,24*days_to_consider+1)
# Declare parameters
focc = m.Param(value = [0,0.05,0.1,\
0.1,0.1],name="focc")
setpoint = [16.00,16.00,\
16.00,16.38,17.43,17.88,18.08,18.35,\
17.85,17.13,16.27,16.00]
T_air_set = m.Param(value=setpoint,name="T_air_set")
Te = m.Param(value = weatherData[0:24*days_to_consider+1,0],name="Te")
qv_infil_wind = m.Param(value = 0.0769*0.6*(0.75*0.8* \
weatherData[0:24*days_to_consider+1,1]**2)**0.667,\
name="qv_infil_wind")
# 1) Airflow
qv_mech_sup = m.Intermediate(2.519244226731981*focc*1,name="qv_mech_sup")
qv_infil_stack = m.Intermediate(0.0146*0.6*(0.7*8.5* \
m.abs(Te-T_air_set))**0.667,\
name="qv_infil_stack")
mx = m.max3(qv_infil_stack,qv_infil_wind)
qv_infil_sw = m.Intermediate(mx+0.14*qv_infil_stack*qv_infil_wind/0.6,\
name="qv_infil_sw")
m.options.IMODE = 6
m.options.DIAGLEVEL = 0
m.options.SOLVER = 1
m.options.MAX_ITER = 1000
m.solve(disp=True,GUI=False)
plt.plot(m.time,mx.value,color='orange',lw=2)
plt.plot(m.time,qv_infil_stack.value,'r:')
plt.plot(m.time,qv_infil_wind.value,'k:')
plt.show()
我知道这是问题的简化版本,但如果所有这些变量都基于常量 m.Param()
类型,您可以将所有这些变量预先计算为 Numpy 数组。这将避免隐式解决方案中额外的不等式约束和二元变量。