I have asked about a year ago about a dispatch scheduler I’ve written using Gekko. Over the previous months I have modified the original code several times to fit what I want and needed adding new constraints and more variables. The code arrives at a solution that is acceptable. However, the run time takes about 10-12 mins, sometimes more. I would like to ask if there is any part in my code that I can change to significantly reduce run time.
I have tried simplifying everything but I still cannot get it to arrive at a solution faster. I understand it could be that the problem is too large and that the long run time is expected. However, as I am not a Gekko expert and this is the first problem I’ve really used this on, there might be workarounds that I am not aware of. I have tried using sum instead of m.sum as Ive read it as a solution in another question, but it didnt really reduce the runtime by much. I have posted the code below.
from gekko import GEKKO
import numpy as np
import pandas as pd
import itertools
pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', None)
m = GEKKO(remote=False) # Initialize gekko
#Sample Demand Without RE
# grid_load = [94.17, 83.35, 79.16, 78.43, 91.74, 75.26, 93.83, 70.09, 81.09, 85.11, 99.24, 97.26, 94.66, 97.86, 96.67, 95.06, 82.66, 84.26, 110.97, 104.45, 102.45, 108.59, 97.7, 90.86]
# grid_load = [34.17, 32.35, 31.16, 30.43, 29.74, 30.26, 31.83, 34.09, 35.09, 36.11, 37.24, 38.26, 37.66, 36.86, 37.67, 37.06, 36.66, 38.26, 45.97, 46.45, 45.45, 43.59, 38.7, 35.86]
grid_load = [57.34, 54.16, 48.35, 45.55, 44.04, 44.33, 45.18, 45.43, 52.61, 56.56, 60.64, 63.12, 67.55, 68.33, 69.09, 68.58, 67.15, 63.55, 63.14, 71.09, 71.67, 68.24, 66.14, 60.45,]
renewables = [
[6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 3.7, 3.7, 3.7, 3.7], #HYDRO1
[5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7], #HYDRO2
[10, 9, 8, 7, 6.5, 6.5, 7.5, 9, 9, 9, 9.5, 9.5, 9.5, 9.5, 9.5, 9, 9, 9.5, 12, 12, 12, 12, 12, 10.5] #WIND
]
demand = grid_load - np.sum(renewables,axis=0)
# print(demand)
n_units_per_plant = [4,5,5,5,7,4,5,4,4,5,4,7,6]
num_plants = len(n_units_per_plant)
pmax = [3.75,3.75,3.75,3.75, #DMCI 4 units
1.33,1.33,1.33,1.33,1.33, #ORMIN 5 units
2.8,3.5,1.8,1.8,1.8, #PONE 5 units
2,1.4,1.4,1.4,1.4, #MHEC 5 units
0.7,1,0.9,1,1,0.35,0.8, #TTR 7 units
1,1,1,1, #MSS 4 units
0.8,0.8,0.8,0.8,1.1, #PSI69 5 units
0.8,0.8,0.8,0.5, #PSI138 4 units
1,1,1,1, #PP 4 units
0.9,0.9,0.9,0.9,0.9, #UPRC 5 units
0.85,0.85,0.85,0.85, #UPRV 4 units
0.35,0.35,0.6,0.4,0.85,0.85,0.85, #TCGR 7 units
0.85,0.85,0.85,0.85,0.75,1, #TTC 6 units
]
pmin = [2.5,2.5,2.5,2.5, #DMCI 4 units (15 MW)
1.33,1.33,1.33,1.33,1.33, #ORMIN 5 units (9 MW)
2.8,3.5,1.8,1.8,1.8, #PONE 5 units
2,1.4,1.4,1.4,1.4, #MHEC 5 units
0.7,1,0.9,1,1,0.35,0.8, #TTR 7 units
1,1,1,1, #MSS 4 units
0.8,0.8,0.8,0.8,1.1, #PSI69 5 units
0.8,0.8,0.8,0.5, #PSI138 4 units
1,1,1,1, #PP 4 units
0.9,0.9,0.9,0.9,0.9, #UPRC 5 units
0.85,0.85,0.85,0.85, #UPRV 4 units
0.35,0.35,0.6,0.4,0.85,0.85,0.85, #TCGR 7 units
0.85,0.85,0.85,0.85,0.75,1, #TTC 6 units
]
v_cost = [14.52,14.52,14.52,14.52,
16.545,16.545,16.545,16.545,16.545,
18.000178,18.000178,18.000178,18.000178,18.000178,
17.953,17.953,17.953,17.953,17.953,
18.531,18.531,18.531,18.531,18.531,18.531,18.531, #TTR
18.182,18.182,18.182,18.182, #MSS
18.12,18.12,18.12,18.12,18.12, #PSI69
18.12,18.12,18.12,18.12, #PSI138
18.23,18.23,18.23,18.23,
22.23,22.23,22.23,22.23,22.23, #UPRC
22.23,22.23,22.23,22.23, #UPRV
18.531,18.531,18.531,18.531,18.531,18.531,18.531, #TCGR
18.531,18.531,18.531,18.531,18.531,18.531, #TTC
]
startup_cost = [200] * len(pmax)
meots = [180,160,221,0,59.5,59.5,59.5,72,51,59.5,59.5,76.5,43.35]
min_ol = [2,1,1,1,0,0,0,0,0,0,0,0,0]
availability = [[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], #DMCI
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], #OPI
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], #PONE
[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], #MHEC
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], #TTR
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], #MSS
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], #PSI69
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], #PSI138
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], #PP
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], #UPRC
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], #UPRV
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], #TCGR
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], #TTC
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
]
availability = np.array(availability).T
# print(availability.T)
n_plants = len(pmax)
hours_tot = len(demand)
unit_indexes = []
start_index = 0
start_index = 0
for num_units in n_units_per_plant:
end_index = start_index + num_units
unit_indexes.append([start_index, end_index])
start_index = end_index
unit_indexes = np.array(list(map(list, zip(*unit_indexes)))) # Transpose the list
# print(unit_indexes[0,1])
print('Total Dependable Demand', sum(pmax))
def solver(slice,demandido,meoty):
hours = len(demandido)
print(meoty)
m.options.SOLVER=1 # APOPT is an MINLP solver
m.options.IMODE= 3
# Initialize Variables
prod = m.Array(m.Var,(hours,n_plants),lb=0, integer=False)
switch_on = m.Array(m.Var,(hours,n_plants),lb=0,ub=1,integer=True)
plant_stat = m.Array(m.Var,(hours,n_plants),lb=0,ub=1,integer=True)
# meot_stat = m.Array(m.Var,(1,num_plants),lb=0,ub=1,integer=True)
for i in range(hours):
#DEMAND CONSTRAINT
# print('DEMAND:',demandido[i])
m.Equation(m.sum(prod[i,:]) >= demandido[i])
#RESERVE CONSTRAINT: Current Dependable Capacity Requires a 2MW reserve on the Regulating Plant
m.Equation(m.sum(np.multiply(pmax[0:4]-prod[i,0:4],plant_stat[i,0:4])) >= 2)
#MINIMUM ONLINE UNITS CONSTRAINT: ORMIN AND POWER ONE
for g in range(num_plants):
if(min_ol[g]!=0):
m.Equation(m.sum(np.multiply(plant_stat[i,unit_indexes[0,g]:unit_indexes[1,g]],availability[i,unit_indexes[0,g]:unit_indexes[1,g]])) >= min_ol[g])
for j in range(n_plants):
#CHECK IF THE UNITS ARE AVAILABLE. IF A AVAILABLE SET THE CONSTRAINTS. ELSE MAKE IT ZERO
m.Equation(prod[i,j] <= pmax[j]*availability[i,j]*plant_stat[i,j]*1.02) #2% tolerance
m.Equation(prod[i,j] >= pmin[j]*availability[i,j]*plant_stat[i,j]*0.98) #2% tolerance
#Plant State constraint
q = m.if2(prod[i,j],0,1)
m.Equation(plant_stat[i,j] == q)
#SWITCHING CONSTRAINT
if i == 0:
m.Equation(switch_on[i,j] == plant_stat[i,j])
else:
m.Equation(switch_on[i,j] >= plant_stat[i,j] - plant_stat[i-1,j])
m.Equation(switch_on[i,j] <= 1 - plant_stat[i-1,j])
for g in range(num_plants):
if(min_ol[g]!=0):
m.Equation(m.sum(np.multiply(plant_stat[i,unit_indexes[0,g]:unit_indexes[1,g]],availability[i,unit_indexes[0,g]:unit_indexes[1,g]])) >= min_ol[g])
if(meots[g]>0):
start = unit_indexes[0,g]
end = unit_indexes[1,g]
#DAILY MEOT CONSTRAINT
if g == 2: #PONE
pia = (np.sum(prod[0:hours, unit_indexes[0, g]:unit_indexes[1, g]]) + np.sum(prod[0:hours, unit_indexes[0, g+1]:unit_indexes[1, g+1]]))
rig = sum(sum(np.multiply(availability[0:hours,unit_indexes[0,g]:unit_indexes[1,g]], pmax[unit_indexes[0,g]:unit_indexes[1,g]]))) +
sum(sum(np.multiply(availability[0:hours,unit_indexes[0,g+1]:unit_indexes[1,g+1]], pmax[unit_indexes[0,g+1]:unit_indexes[1,g+1]])))
else:
pia = np.sum(prod[0:hours,start:end]) #24 hours production
rig = sum(sum(np.multiply(availability[0:hours,start:end], pmax[start:end]))) #maximum allowable dispatch as per unit pmax and availability
mik = meoty[g]/(hours_tot/hours) #plant meot per demand slice
# print("START, END, RIG, MIK:", start, end, rig, mik)
#Sanity Checker for MEOT
#Check if its possible to reach meot constraint. Waive if not possible
if rig >= mik:
m.Equation(pia >= mik)
else:
m.Equation(pia >= rig)
#OBJECTIVE FUNCTION
abc = [np.sum(prod[:, unit_indexes[0,g]:unit_indexes[1,g]]) for g in range(num_plants)]
meot_stat = [m.if2(meots[i]-abc[i], 0, 1) for i in range(num_plants)] #if prod of plant is greater than their meot its zero
# start_up_cost = m.Intermediate(m.sum([switch_on[t,i]*0.1*v_cost[i]*pmax[i] for i in range(n_plants) for t in range(hours)])) #starup cost is 10% of the fullload cost
start_up_cost = m.Intermediate(m.sum([switch_on[t,i]*startup_cost[i] for i in range(n_plants) for t in range(hours)])) #starup cost is 10% of the fullload cost
prod_cost = m.Intermediate(m.sum([prod[t,i]*v_cost[i] for i in range(n_plants) for t in range(hours)]))
meot_cost = m.Intermediate(m.sum([(meots[i] - abc[i])*meot_stat[i]*v_cost[unit_indexes[0,i]] for i in range(num_plants)]))
# reserve - m.sum(reserve)
total_cost = start_up_cost +
prod_cost
m.Minimize(total_cost)
m.Minimize(meot_cost)
m.solve(disp=True) # Solve
print("Total Cost", m.options.OBJFCNVAL-start_up_cost.value[0])
# meot_cost_value = meot_cost.value[0]
#OUTPUTS
prody = np.round(np.array([prod[i][j].value[0] for j in range(n_plants) for i in range(hours)]).reshape(n_plants,hours,),4)
plant_stat = np.array([plant_stat[i][j].value[0] for j in range(n_plants) for i in range(hours)]).reshape(n_plants,hours,)
switch_on = np.array([switch_on[i][j].value[0] for j in range(n_plants) for i in range(hours)]).reshape(n_plants,hours,)
return(plant_stat,switch_on,prody)
slice = len(demand)
# # print(hours_tot)
reshaped_prody_list = []
updated_meots = meots
x, y, prody = solver(slice, demand, updated_meots)
# # Convert dispatch_array to a DataFrame
dispatch_schedule = pd.DataFrame(prody)