У меня есть два вопроса о библиотеке Gekko на Python.
Есть ли способ повысить производительность приведенного ниже кода (Input1, Input2)? Код был решен быстро, и результаты были правильными, когда строки 170 ~ 195 были исключены (прокомментированы). Однако, когда я выполнил весь код (включая строки 170–195), производительность резко снизилась, и я не получил результатов после ожидания более 30 минут. Я предположил, что это было из-за функций
if3в строках 170 и 171, потому что, как только я включил эти строки, код не мог быть решен.Несмотря на то, что я определил переменную
days_to_considerв строках с 61 до 1, я увеличу значение переменной до 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, 1, 1, 2.25, 3.28, 0, 0.6, 1, 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%', 'Automatic control more than 50%', 'Taps More Than 3m from Heat Generation', 'Co-Generation (0.9)', '20', '', 'Electricity', 'Electricity', 'Electricity']
weatherData = np.load("Atlanta_TMY3_climate.npy")
weatherData = np.insert(weatherData, 0, 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]*(0.7*left_numeric[1,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,0], 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, 0, 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,0])*totalArea/1000,\
## 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]*left_numeric[10,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]*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(m.max3(V_heating_sup_air, V_cooling_sup_air),qv_mech_sup*(1-HR_efficiency))+left_numeric[7,0]*3.6/totalArea*focc)*left_numeric[9,0]*left_numeric[10,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)
- Я планирую провести динамические задачи оптимального управления, которые будут состоять из всего ~ 45 MV от ~ 10 зданий (код показан выше, каждое здание имеет 4-5 MV) и 1 PV сообщества (1 MV, модель показано здесь) для 8760 временных шагов (= 365 дней x 24 часа) . Считаете ли вы, что проблема такого размера решается с помощью решателя IPOPT в течение нескольких часов на обычном персональном ноутбуке (или, если доступна сложность алгоритма решателя IPOPT, я был бы признателен)? Я думаю, что это будет сложно решить, потому что количество MV становится примерно ~ 400000 (= ~ 45 MV x 8760 часов), и мне интересно, сможет ли решатель IPOPT справиться с задачей оптимизации такого размера.
