Gekko：使用 if3 函数后代码性能下降

如何解决Gekko：使用 if3 函数后代码性能下降

关于 Python 中的 Gekko 库，我有两个问题。

1. 有什么办法可以提高下面代码的性能（Input1，Input2）？代码求解很快，排除170~195行（注释）结果正确。但是，当我执行整个代码（包括第170-195行）时，性能急剧下降，等待30多分钟后没有得到结果。我推测是因为第 170 行和第 171 行中的 if3 函数，因为一旦我包含这些行，代码就无法解决。

   即使我在第 61 行到 1 行定义了 days_to_consider 变量，但我会在代码成功执行后将变量值增加到 365 天，这意味着解决这个代码性能问题对我来说至关重要.

import numpy as np
from gekko import GEKKO

np.set_printoptions(linewidth=2000)
np.set_printoptions(formatter={'float': '{: 0.2f}'.format})

# Given parameters (Please ignore,this is just given input parameters)
left_numeric = np.zeros((14,1))
tempo = [3738.0,8.5,1,2.25,3.28,0.6,0.5,4,285000,0]
for i in range(14):
    left_numeric[i,0] = tempo[i]

left_string = ['Urban / City','Light : 110,000 * Af','Yes','23. Variable air volume / Water or Water&Air / Air / Yes','1. Mechanical system only; no provision for natural ventilation','Slowly rotating or intermittent heat exchangers (0.7)','No exhaust air recirculation','Automatic control more than 50%','Taps More Than 3m from Heat Generation','Co-Generation (0.9)','20','','Electricity','Electricity']

weatherData = np.load("Atlanta_TMY3_climate.npy")
weatherData = np.insert(weatherData,np.zeros((1,15)),0)

Info_Array = np.load("Info_Array.npy")

vsite = 0.8
totalAppliance = 12.0
totalOccupants = 8.39748075577327
hs_unoccupied = 16.0
Cm = 30.55555555
Htr_ms = 22.75
Htr_w = 1.5824055209225467
Htr_is = 15.525
Htr_em = 1.0723552727080168
f_BAC_hc = 0.7
f_BAC_e = 0.87
fcntrl_heat = 0.5
fcntrl_cool = 0.5
dist_heat = 0.9259259
dist_cool = 1.00
totalArea = 877.4

total_DHW_system = 76.40691666666666
effi_gen_DHW = 0.9
occu_equi_hours = 3093.5

Prs = 0.40580206298113436
Prm = 0.5555555555555556
HR_efficiency = 0.7

T_supply_h = 28.0
T_supply_c = 17.0

ventSupply = 2.519244226731981
ventRecirculation = 1
ventilation_cooling_type = 1

dist_heat = 0.9259259
dist_cool = 1.00000


# Create GEKKO model
m = GEKKO(remote=False)
#m.open_folder()
print(m.path)

days_to_consider = 1
m.time = np.linspace(0,24*days_to_consider,24*days_to_consider+1)

# Define MVs       
fDim = m.MV(lb=1,ub=1,name="fDim") # dummy MV just for testing. In the future more MVs will be added and the interval will also be modified.
fDim.STATUS = 1

f_VT = 1 # this will be an MV,but for Now,it is just a constant

# 2D numppy array manipulation (add zero row in the first row of the "Info_Array" and "weatherData")
T_heating_set = m.Param(value = Info_Array[0:24*days_to_consider+1,9],name="T_heating_set")
T_cooling_set = m.Param(value = Info_Array[0:24*days_to_consider+1,10],name="T_cooling_set")

# Define parameters     
Te          = m.Param(value = weatherData[0:24*days_to_consider+1,0],name="Te")
Wind_speed  = m.Param(value = weatherData[0:24*days_to_consider+1,1],name="Wind_speed")
fapp        = m.Param(value = Info_Array[0:24*days_to_consider+1,6],name="fapp")
Qapp = totalAppliance 
focc        = m.Param(value = Info_Array[0:24*days_to_consider+1,5],name="focc")
Qocc = totalOccupants 
Qlight      = m.Param(value = Info_Array[0:24*days_to_consider+1,8],name="Qlight")
Qsol        = m.Param(value = Info_Array[0:24*days_to_consider+1,11],name="Qsol")

qv_infil_wind   = m.Param(value = 0.0769*left_numeric[8,0]*(0.75*vsite*weatherData[0:24*days_to_consider+1,1]**2)**0.667,name="qv_infil_wind")

# Define variables
qv_infil_stack    = m.Var(0,name="qv_infil_stack")

T_m0_t        = m.Var(name="T_m0_t")
T_m10_t       = m.Var(name="T_m10_t")
T_m0          = m.Var(name="T_m0")
T_m10         = m.Var(name="T_m10")
T_s0          = m.Var(name="T_s0")
T_s10         = m.Var(name="T_s10")
T_air0        = m.Var(value=hs_unoccupied,name="T_air0")
T_air10       = m.Var(name="T_air10")
T_m_ac_t      = m.Var(value=hs_unoccupied,name="T_m_ac_t")
T_m_ac        = m.Var(value=hs_unoccupied,name="T_m_ac")
T_s_ac        = m.Var(name="T_s_ac")
T_air_ac      = m.Var(name="T_air_ac")
T_air_set     = m.Var(name="T_air_set")

T_air_set_prev = m.Var(name="T_air_set_prev")
m.delay(T_m_ac_t,T_air_set_prev,1)

Q_HC_nd     = m.Var(value=0,name="Q_HC_nd")


# Define delivered energy
E_heating   = m.Var(name="E_heating")
E_cooling   = m.Var(name="E_cooling")
E_lighting  = m.Var(name="E_lighting")
E_fan       = m.Var(name="E_fan")
E_pump      = m.Var(name="E_pump")
E_equip     = m.Var(name="E_equip")
E_dhw       = m.Var(name="E_dhw")

V_heating_sup_air = m.Var(name="V_heating_sup_air")
V_cooling_sup_air = m.Var(name="V_cooling_sup_air")
q_v_sup           = m.Var(name="q_v_sup")

# Build building Equations
# 0) Internal heat gains
Qint     = m.Intermediate(fapP*Qapp + focc*Qocc + fDim*Qlight,name="Qint")
Qia      = m.Intermediate(0.5*Qint,name="Qia")
Qst      = m.Intermediate(Prs*(Qia + Qsol),name="Qst")
Qm       = m.Intermediate(Prm*(Qia + Qsol),name="Qm")

# 1) Airflow
qv_mech_sup    = m.Intermediate(ventSupply*ventRecirculation*focc,name="qv_mech_sup")

m.Equation(qv_infil_stack == 0.0146*left_numeric[8,0]*(0.7*left_numeric[1,0]*m.abs(Te-T_air_set))**0.667)
##qv_infil_stack = m.Intermediate(0.0146*left_numeric[8,0]*m.abs(Te-T_air_set))**0.667,name="qv_infil_stack")

qv_infil_sw    = m.Intermediate(m.max3(qv_infil_stack,qv_infil_wind) + 0.14*qv_infil_stack*qv_infil_wind/left_numeric[8,name="qv_infil_sw")

Hve      = m.Intermediate((qv_mech_suP*(1-HR_efficiency) + qv_infil_sw)*0.3333,name="Hve")
Htr_1    = m.Intermediate(Hve*Htr_is/(Hve + Htr_is),name="Htr_1")

# 2) Temperature
Tm_denom = m.Intermediate(Htr_ms + Htr_w + Htr_1,name="Tm_denom")
T_m_intermediate = m.Intermediate(0.5*((Htr_1+Htr_w)*Htr_ms/Tm_denom+Htr_em),name="T_m_intermediate")

m.Equations([\
T_m0_t  == (1/(Cm+T_m_intermediate))*((Cm-T_m_intermediate)*T_air_set_prev+(Htr_ms*(Htr_w+Htr_1)/Tm_denom+Htr_em)*Te+(Htr_ms/Tm_denom)*Qst+(Htr_ms*Htr_1)/Hve/Tm_denom*(Qia)+Qm),\
T_m10_t == (1/(Cm+T_m_intermediate))*((Cm-T_m_intermediate)*T_air_set_prev+(Htr_ms*(Htr_w+Htr_1)/Tm_denom+Htr_em)*Te+(Htr_ms/Tm_denom)*Qst+(Htr_ms*Htr_1)/Hve/Tm_denom*(Qia+10)+Qm),\
T_m0  == 0.5*(T_air_set_prev+T_m0_t),\
T_m10 == 0.5*(T_air_set_prev+T_m10_t),\
T_s0 == (Htr_ms*T_m0 +(Htr_w+Htr_1)*Te+Qst+Htr_1*Qia/Hve)/Tm_denom,\
T_s10 == (Htr_ms*T_m10+(Htr_w+Htr_1)*Te+Htr_1*10/Hve+Qst+Htr_1*Qia/Hve)/Tm_denom,\
T_air0 == (Htr_is*T_s0 +Hve*Te+Qia)/(Htr_is+Hve),\
T_air10 == (Htr_is*T_s10+Hve*Te+10+Qia)/(Htr_is+Hve)])
 
# 3. Heating/Cooling need
# Tair set & Q_HC_nd
m.Equations([T_air_set == m.if3(T_air0-T_heating_set,T_heating_set,m.if3(T_cooling_set-T_air0,T_cooling_set,T_air0)),\
             Q_HC_nd == 10*(T_air_set-T_air0)/(T_air10-T_air0)])

# 4. Tac & Q_nd
m.Equations([\
T_m_ac_t == (1/(Cm+T_m_intermediate))*((Cm-T_m_intermediate)*T_air_set_prev+(Htr_ms*(Htr_w+Htr_1)/Tm_denom+Htr_em)*Te+(Htr_ms/Tm_denom)*Qst+(Htr_ms*Htr_1)/Hve/Tm_denom*(Qia+Q_HC_nd)+Qm),\
T_m_ac == 0.5*(T_air_set_prev+T_m_ac_t),\
T_s_ac == (Htr_ms*T_m_ac+(Htr_w+Htr_1)*Te+Htr_1*Q_HC_nd/Hve+Qst+Htr_1*Qia/Hve)/Tm_denom,\
T_air_ac == (Htr_is*T_s_ac+Hve*Te+Q_HC_nd+Qia)/(Htr_is+Hve)])

m.Minimize(Q_HC_nd) # please use (uncomment) this line when the upper part of the code (line 1-163) is executed. 

#_________________________When below code is included,the solver Couldn't solve________________________#

##Q_heat_nd = m.Intermediate(m.if3(Q_HC_nd,Q_HC_nd),name="Q_heat_nd") # SUSPECT that this line slows down the code performance
##Q_cool_nd = m.Intermediate(m.if3(Q_HC_nd,-Q_HC_nd,0),name="Q_cool_nd")
##
### 5. System delivered energy use (Unit: kW)
##m.Equations([E_heating == Q_heat_nd*f_BAC_hc/(dist_heat*left_numeric[5,0])*totalArea/1000,\
##             E_cooling == f_VT*Q_cool_nd*f_BAC_hc/(dist_cool*left_numeric[6,\
##             E_lighting == fDim*Qlight*totalArea/1000,\
##             E_equip    == fapP*QapP*totalArea/1000,\
##             E_pump == 8*(fcntrl_heat+fcntrl_cool)/3.6/occu_equi_hours*totalArea,\
##             E_dhw == total_DHW_system/effi_gen_DHW/occu_equi_hours*totalArea/1000,\
##             V_heating_sup_air == Q_heat_nd*0.0036/(T_supply_h-T_air_ac)/0.001239913,\
##             V_cooling_sup_air == Q_cool_nd*0.0036/(T_air_ac-T_supply_c)/0.001239913])
##
##
##if ventilation_cooling_type == 3 and left_numeric[7,0] == 0:
##    m.Equation(E_fan == (m.max3(V_heating_sup_air,V_cooling_sup_air))*left_numeric[9,0]*left_numeric[10,0]*f_BAC_e*totalArea/1000)
##
##elif ventilation_cooling_type == 3 and left_numeric[7,0] != 0:
##    m.Equation(E_fan == (m.max3(V_heating_sup_air,V_cooling_sup_air)+left_numeric[7,0]*3.6/totalArea*focc)*left_numeric[9,0]*f_BAC_e)
##
##elif ventilation_cooling_type != 3 and left_numeric[7,0] == 0:
##    m.Equation(E_fan == (m.max3(m.max3(V_heating_sup_air,V_cooling_sup_air),qv_mech_suP*(1-HR_efficiency)))*left_numeric[9,0]*f_BAC_e*totalArea/1000)
##
##elif ventilation_cooling_type != 3 and left_numeric[7,0] != 0:
##     m.Equation(E_fan == (m.max3(m.max3(V_heating_sup_air,qv_mech_suP*(1-HR_efficiency))+left_numeric[7,0]*f_BAC_e)
##m.Minimize(E_heating + E_cooling + E_lighting + E_fan + E_equip + E_pump + E_dhw) # please use (uncomment) this line when the entire code is executed.

#_______________________________________________________________________________________________________#
m.options.IMODE = 6
m.options.DIAGLEVEL = 4
m.options.soLVER = 1
m.options.MAX_ITER = 1000000
m.solve(disp=True,GUI=False,debug=False)

2. 我计划进行动态最优控制问题，这将包含来自约 10 个建筑物的总共约 45 个 MV（代码如上所示，每个建筑物有 4-5 个 MV）和 1 个社区 PV（1 个 MV，该模型显示为 here) 8,760 个时间步长（= 365 天 x 24 小时）。您认为在普通的个人笔记本电脑上使用 IPOPT 求解器可以在几个小时内解决这个规模的问题吗（或者如果 IPOPT 求解器算法的复杂性可用，我将不胜感激）？我认为这很难解决，因为 MV 的数量大约为 400,000（=~45 MV x 8,760 小时），我想知道 IPOPT 求解器是否可以处理这种规模的优化问题。

解决方法

注释代码有 6 个 max3() 函数和 2 个 if3() 函数。那些为每个时间点添加一个二进制变量，因此总共 25 x 8 = 200 个二进制变量用于 1 天和 (24*365 + 1) x 8 = 70,088 个二进制变量用于 365 天。长达一年的解决方案可能需要太多时间。以下是一些改善求解时间的建议：

• 使用 m.options.DIAGLEVEL=0 而不是 4。额外的诊断需要更长的预处理时间。
• 如果可能，尝试使用 m.options.REDUCE=3 进行预处理以减少方程的数量。
• 用 m.options.SOLVER=3 (IPOPT) 初始化，然后切换到 m.options.SOLVER=1 (APOPT) with 以获得整数解。再次求解时使用 m.options.TIME_SHIFT=0 以防止初始条件在第二次求解时发生变化。
• 尝试将 if2() 和 max2() 作为 MPCC 形式，它们比 if3() 和 max3() 求解速度更快，尤其是对于大型问题。如果解决方案是在切换条件下，他们有时会遇到问题。
• 使用 .STATUS=0 关闭自由度，首先在没有决策变量的情况下求解以获得初始可行解。

未注释行的增强问题似乎不可行。 APOPT 和IPOPT 都报告问题不可行。

EXIT: Converged to a point of local infeasibility. Problem may be infeasible.

 An error occured.
 The error code is  2

 
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :  86.0874 sec
 Objective      :  165.28211287911208
 Unsuccessful with error code  0
 ---------------------------------------------------

文件 infeasibilities.txt 将包含其他见解。