I am setting up a big time-dependent energy optimization problem over many timesteps with pyomo.The problem is formulated using indexed blocks, as proposed in the Pyomo Book in chapter 8.6.2., where the model object contains a block for each timestep.
I notice now that the model creation using indexed blocks takes much longer than the creation of a model with a simple formulation without blocks.
I could not find any remarks on the performance of indexed blocks vs. standard formulation in the pyomo book or other posts.
Do you agree on these findings that standard formulation is faster then indexed blocks in case of biggish optimization problems with many timesteps?
Here is my code example for a PV-Battery system with cost minimzation. The standard formualation follows the indexed Block formulation.
import pyomo.environ as pyo
from pyomo.common.timing import TicTocTimer
import random
##% Generate some random data for PV and Load
pv = [random.randint(0, 5) for _ in range(8760)]
pv_dict = (dict(enumerate(pv,1)))
load_el = [random.randint(0, 8) for _ in range(8760)]
load_el_dict = (dict(enumerate(load_el,1)))
#%%
# Create standard formulation pyomo model
timer = TicTocTimer()
timer.tic('Start Sandard formulation')
## Define model
model = pyo.ConcreteModel()
# Define timeperiod set
model.T = pyo.RangeSet(len(pv_dict))
# Define model parameters
model.pv = pyo.Param(model.T, initialize=pv_dict)
model.load_el= pyo.Param(model.T, initialize=load_el_dict)
model.grid_cost_buy = pyo.Param(model.T, initialize=0.4)
model.battery_eoCH = pyo.Param(initialize=1.0)
model.battery_eoDCH = pyo.Param(initialize=0.1)
model.battery_capacity = pyo.Param(initialize=5)
# Define the variables
model.battery_power_CH = pyo.Var(model.T, domain=pyo.NonNegativeReals) # battery charging power
model.battery_power_DCH = pyo.Var(model.T, domain=pyo.NonNegativeReals) # battery discharging power
model.battery_soc = pyo.Var(model.T, bounds=(model.battery_eoDCH, model.battery_eoCH)) # battery soct with end of ch/DCH levels
model.grid_power_import = pyo.Var(model.T, domain=pyo.NonNegativeReals) # grid import power
model.grid_power_export = pyo.Var(model.T, domain=pyo.NonNegativeReals) # grid export power
## Define battery constraints
# Battery end of CH/ constraints
def battery_end_of_CH_rule(m, t):
return m.battery_soc[t] <= model.battery_eoCH
model.battery_eoCH_c = pyo.Constraint(model.T, rule=battery_end_of_CH_rule)
# Battery end of DCH constraints
def battery_end_of_DCH_rule(m, t):
return m.battery_soc[t] >= model.battery_eoDCH
model.battery_eoDCH_c = pyo.Constraint(model.T, rule=battery_end_of_DCH_rule)
# Define battery SoC constraint
def battery_soc_rule(m, t):
if t == m.T.first():
return m.battery_soc[t] == ((m.battery_power_CH[t] - m.battery_power_DCH[t]) / model.battery_capacity)
return m.battery_soc[t] == m.battery_soc[t-1] + ((m.battery_power_CH[t] - m.battery_power_DCH[t]) / model.battery_capacity)
model.battery_soc_c = pyo.Constraint(model.T, rule=battery_soc_rule)
# Define balanced electricity bus rule
def balanced_bus_rule(m, t):
return (0 == (m.pv[t] - m.load_el[t]
+ m.battery_power_DCH[t] - m.battery_power_CH[t]
+ m.grid_power_import[t] - m.grid_power_export[t]))
model.bus_c = pyo.Constraint(model.T, rule=balanced_bus_rule)
## Define the cost function
def obj_rule(m):
return sum(m.grid_power_import[t]*m.grid_cost_buy[t] for t in m.T)
model.obj = pyo.Objective(rule=obj_rule, sense=1)
timer.toc('Built model')
## Solve the problem
solver = pyo.SolverFactory('gurobi')
results = solver.solve(model)#, report_timing=True)#, tee=True)
timer.toc('Wrote LP file and solved')
print('Total operation costs:',pyo.value(model.obj))
#%%
# Create time indexed Block formulation
timer.tic('Start Indexed Block formulation')
# Define model
model = pyo.ConcreteModel()
# Define timeperiods
model.T = pyo.RangeSet(len(pv_dict))
def block_rule(b, t):
# Define model parameters
b.pv = pyo.Param(initialize=pv_dict[t])
b.load_el = pyo.Param(initialize=load_el_dict[t])
b.grid_cost_buy = pyo.Param(initialize=0.4)
b.battery_eoCH = pyo.Param(initialize=1.0)
b.battery_eoDCH = pyo.Param(initialize=0.1)
b.battery_capacity = pyo.Param(initialize=5)
# define the variables
b.battery_power_CH = pyo.Var(domain=pyo.NonNegativeReals) # battery charging power
b.battery_power_DCH = pyo.Var(domain=pyo.NonNegativeReals) # battery discharging power
b.battery_soc = pyo.Var(bounds=(b.battery_eoDCH, b.battery_eoCH)) # battery soct with end of ch/DCH levels
b.battery_soc_initial = pyo.Var() # Initial battery soc
b.grid_power_import = pyo.Var(domain=pyo.NonNegativeReals) # grid import power
b.grid_power_export = pyo.Var(domain=pyo.NonNegativeReals) # grid export power
# define the constraints
# Balanced electricity bus rule
def balanced_bus_rule(_b):
return (0 == (_b.pv - _b.load_el
+ _b.battery_power_DCH - _b.battery_power_CH
+ _b.grid_power_import - _b.grid_power_export))
b.bus_c = pyo.Constraint(rule=balanced_bus_rule)
# Battery end of CH/ constraints
def battery_end_of_CH_rule(_b):
return (_b.battery_eoCH >= _b.battery_soc)
b.battery_eoCH_c = pyo.Constraint(rule=battery_end_of_CH_rule)
# Battery end of DCH constraints
def battery_end_of_DCH_rule(_b):
return (_b.battery_eoDCH <= _b.battery_soc)
b.battery_eoDCH_c = pyo.Constraint(rule=battery_end_of_DCH_rule)
# Battery SoC constraint
def battery_soc_rule(_b):
return (_b.battery_soc == _b.battery_soc_initial + ((_b.battery_power_CH - _b.battery_power_DCH) / _b.battery_capacity))
b.battery_soc_c = pyo.Constraint(rule=battery_soc_rule)
# Initialize Blocks for each timestep defined in T
model.pvbatb = pyo.Block(model.T, rule=block_rule)
# link the battery SoC variables between blocks
# setting the initial SoC of one block equal to the final SoC of the previous block
def battery_soc_linking_rule(m, t):
if t == m.T.first():
return m.pvbatb[t].battery_soc_initial == 0
return m.pvbatb[t].battery_soc_initial == m.pvbatb[t-1].battery_soc
model.battery_soc_linking = pyo.Constraint(model.T, rule=battery_soc_linking_rule)
## define the cost function
def obj_rule(m):
return sum(m.pvbatb[t].grid_power_import * m.pvbatb[t].grid_cost_buy for t in m.T)
model.obj = pyo.Objective(rule=obj_rule, sense=1)
timer.toc('Built model')
## solve the problem
solver = pyo.SolverFactory('gurobi')
results = solver.solve(model)#, report_timing=True)
timer.toc('Wrote LP file and solved')
print('Total operation costs:',pyo.value(model.obj))
This is the output, which shows that model creation takes about 4x more time with indexed blocks then standard formulation. This further creases with model complexity.
[ 0.00] Start Sandard formulation
[+ 0.92] Built model
[+ 3.37] Wrote LP file and solved
Total operation costs: 5706.199999999935
[ 4.31] Start Indexed Block formulation
[+ 12.50] Built model
[+ 4.62] Wrote LP file and solved
Total operation costs: 5706.199999999934
Thank you so much for your support! Sorry for the not very compact example.
Fabian