0

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

fabmid
  • 11
  • 2
  • the question needs sufficient code for a minimal reproducible example: https://stackoverflow.com/help/minimal-reproducible-example – D.L Feb 28 '23 at 12:23
  • Using blocks will incur an overhead over a "flat" model: both in memory (more modeling objects, and a larger number of "non slotized" objects), and in processing time (the hierarchy has to be traversed looking for objectives, constraints, and variables). That said, if the overhead is "significant", the development team would be interested / appreciate you sharing the complete model to use as a test case to aid in performance profiling. – jsiirola Feb 28 '23 at 16:20
  • @jsiirola The total model is still in development, so I can not yet share it with you. But it will be open source so you may use it in the future. The performance difference I noticed so far with my current model was factor 20 for the model building step. – fabmid Mar 01 '23 at 14:39
  • Interesting. I'm curious why this construct would be suggested in "the book." It would not occur to me to make blocks in a time-stepped problem that can be simply & elegantly expressed with basic time indexing. Of course if there is any overhead at all with the Vars and other elements (likely a lot), the blocked model is going to be massive compared to flat. In the flat, you have 5 variables. In the blocked, you have ~54,000. – AirSquid Mar 01 '23 at 16:10
  • 1
    There are several cases where blocks can be useful. The most common modeling case is when developing frameworks or a hierarchy of models: e.g., I have a model for a process, then I want to make it a multiperiod process. I can rewrite the model and add time indices to all the Vars/Params/Constraints, or I can reuse the original model verbatim and place copies of it on time-indexed blocks. There is overhead (mostly memory) because while `VarData` objects are slotized (making them pretty compact), general containers (e.g., `ScalarVar` and `IndexedVar`) are not (and consume more memory) – jsiirola Mar 01 '23 at 18:18
  • @jsiirola I like the encapsulation concept ^^^. That makes good sense. Maybe I should get the book. :/ – AirSquid Mar 01 '23 at 20:02

0 Answers0