Nonconvex CHP plant vs. running time

Hey community,

I am struggling with the running time of my simulation. If I only consider one day than everything is fine and the simulation reaches the end. But if I consider one week or even longer than the simulation is limitless. Since I can simulate one day, my code should be fine or? The main reason for the long simulation time might be the usage of the chp plant (mixed integer), but nevertheless I was expecting any results after some time, which was not possible. The running time of one night for simulating one week was not enough.

Thanks in advance for your help.

Br,
Alex

import os
import pandas as pd
import pprint as pp
from oemof import solph
from oemof.network.network import Node

solver = "cbc"
number_of_time_steps = 24 * 7

date_time_index = pd.date_range("1/1/2012", periods=number_of_time_steps,
                                freq="H")

# add es
energysystem = solph.EnergySystem(timeindex=date_time_index)
#Node.registry = energysystem

# add buses
brapeseedoil = solph.Bus(label="rapeseedoil")
bel = solph.Bus(label="electricity")
bheat = solph.Bus(label="heat")
energysystem.add(brapeseedoil, bel, bheat)

# add gas
energysystem.add(
    solph.Source(
        label="r_rapeseedoil",
        outputs={brapeseedoil: solph.Flow()},
    )
)

heat_extra_source = energysystem.add(solph.Source(label="heat_extra_source", outputs={
    bheat: solph.Flow(variable_costs=1000000)
    }))

el_extra_source = energysystem.add(solph.Source(label="el_extra_source", outputs={
    bel: solph.Flow(variable_costs=1000000)
    }))

# add heat demand
hd = [2, 2, 2, 2, 2, 2, 3, 3, 4, 4, 5, 7, 7, 5, 4, 4, 3, 3, 2, 2, 2, 2, 2, 2] * 7

#add el demand
eld = [1, 1, 1, 1, 1, 1, 4, 4, 8, 8, 10, 11, 11, 10, 8, 8, 4, 4, 1, 1, 1, 1, 1, 1] * 7

energysystem.add(
    solph.Sink(
        label="demand_heat",
        inputs={bheat: solph.Flow(fix=hd, nominal_value=1)},
    )
)

energysystem.add(
    solph.Sink(
        label="excess_heat",
        inputs={bheat: solph.Flow()},
    )
)

energysystem.add(
    solph.Sink(
        label="demand_el",
        inputs={bel: solph.Flow(fix=eld, nominal_value=1)},
    )
)

energysystem.add(
    solph.Sink(
        label="excess_el",
        inputs={bel: solph.Flow()},
    )
)

# add bhkw
bhkw = energysystem.add(solph.components.GenericCHP(
    label="pp_chp",
    fuel_input={
        brapeseedoil: solph.Flow(
            H_L_FG_share_max=[0.14 for p in range(0, number_of_time_steps)],
            H_L_FG_share_min=[0.14 for p in range(0, number_of_time_steps)],
            variable_costs=0.085,
        )
    },
    electrical_output={
        bel: solph.Flow(
            P_max_woDH=[11 for p in range(0, number_of_time_steps)],
            P_min_woDH=[6 for p in range(0, number_of_time_steps)],
            Eta_el_max_woDH=[0.30 for p in range(0, number_of_time_steps)],
            Eta_el_min_woDH=[0.30 for p in range(0, number_of_time_steps)],
            variable_costs=0.03,
        )
    },
    heat_output={bheat: solph.Flow(
        Q_CW_min=[0 for p in range(0, number_of_time_steps)])},
        Beta=[0 for p in range(0, number_of_time_steps)],
        fixed_costs=0,
        back_pressure=False,
))

"""
pp_chp = energysystem.add(
        solph.Transformer(
        label="pp_chp",
        inputs={brapeseedoil: solph.Flow(variable_costs=0.085)},
        outputs={
            bel: solph.Flow(variable_costs=0.03),
            bheat: solph.Flow(nominal_value=11, min=0.25, nonconvex=solph.NonConvex())
            },
        conversion_factors={bel: 0.30, bheat: 0.56},
   )
)
"""

electric_heater = energysystem.add(
    solph.Transformer(
    label="pp_electric_heater",
    inputs={bel: solph.Flow(nominal_value=7)},
    outputs={bheat: solph.Flow()},
    conversion_factors={bheat: 0.90},
   )
)

storage_el_lithium = energysystem.add(solph.components.GenericStorage(
    nominal_storage_capacity=20,
    label="storage_lithium",
    inputs={bel: solph.Flow()},
    outputs={bel: solph.Flow()},
    loss_rate=0.00001,
    initial_storage_level=0,
    min_storage_level=0,
    max_storage_level=1,
    nominal_output_capacity_ratio=1,
    nominal_input_capactiy_ratio=1,
    inflow_conversion_factor=0.95,
    outflow_conversion_factor=0.95,
))

storage_heat_puffer = energysystem.add(solph.components.GenericStorage(
    nominal_storage_capacity=200,
    label="storage_puffer",
    inputs={bheat: solph.Flow()},
    outputs={
        bheat: solph.Flow()
    },
    loss_rate=0.004,
    initial_storage_level=None,
    inflow_conversion_factor=1,
    outflow_conversion_factor=1,
    nominal_input_capacity_ratio=1,
    nominal_output_capacity_ratio=1,
    min_storage_level=0,
    max_storage_level=1,
))

model = solph.Model(energysystem)
model.solve(solver=solver, solve_kwargs={"tee": True})

results = solph.processing.results(model)

electricity_bus = solph.views.node(results, "electricity")
heat_bus = solph.views.node(results, "heat")

print(electricity_bus["sequences"].sum(axis=0))
print(heat_bus["sequences"].sum(axis=0))

Hey Alexander,

I never used the chp plant, but struggled with long running times for MILP’s before. Generally, I found out that the commercial solver do perform much better when solving MILP. So I tried the Gurobi Solver with your code and found the optimal solution in less than 1s. For educational purpose you can get a free license here:
https://www.gurobi.com/academia/academic-program-and-licenses/
I also tried the free GLPK solver, which atleast found a feasible solution quickly, but taking longer to find global optimum.

Hope this helps you,
Best Henry

Hey Henry,

thanks a lot for your answer! I checked out the Gurobi Solver and the simulation is running! Unbelievable! Since I am new, I thought of having some mistakes in my code. Therefore, I was not expected the problem in the usage of the solver.
Thanks!

Br,
Alex

Dear Alexander,

I am happy that your problem has been solved, in my experience working with MILP can be tricky. I find this paper a good reading material to understand the messages from of the solver and to have some idea on how to deal with the problems. The study case is cplex but gurobi has a similar log structure.