Adding availability constraint with limit_active_flow_count

Q&A
1

Hello oemof community,

I try to model the dispatch of various heat generators in a district heating system. In the following minimal example, I only have a coal-fired CHP and gas boilers. I already successfully implemented constraints such as startup costs and minimum downtime restrictions in order to have a more realistic dispatch. I struggle with implementing an availability constraint and think that the limit_active_flow_count method is the thing I’m looking for.

However, it does not work as I want. In the minimal example with 100 timesteps I want to restrict the coal CHP to be available in only 80 hours. (sum of status variables <= 80) In the output graphic can be seen that the coal CHP is running all the time and that the restriction is not working properly. (sum of status variables = 100) I assume the error is in line 95, but I can’t help myself anymore. Any help appreciated :slight_smile: Thanks

import pandas as pd
from oemof import solph
from oemof.solph import constraints
from oemof.network.network import Node

try:
    import matplotlib.pyplot as plt
except ImportError:
    plt = None

d = {'heatload': [569.894450, 554.2135219,576.837488,592.7529451,618.0502412,632.558263,678.3786056,726.6801345,806.329521,848.8491301,843.8206541,813.0734607,807.3272522,778.3711548,778.4368546,784.8020551,779.5431271,761.8519632,777.5371374,761.8624275,726.3735687,	675.0647848,	618.6561332,	554.6097775,	513.7508669,	539.4350185,	535.7783665,	529.6785499,	549.5875261,	574.5496884,	600.239436,	637.075229,	628.5732578,	601.4101048,	541.8661765,	588.8839689,	598.4106573,	606.8333384,	617.5686688,	639.4357284,	632.9272146,	655.3578554,	702.1020292,	706.0926676,	652.98615,	640.1256493,	582.4460837,	545.8409157,	506.8142917,	472.7983128,	476.7137177,	510.0998035,	550.6710718,	581.9813453,	634.0033481,	686.6767652,	727.9725173,	757.7538129,	767.6656191,	736.8802585,	719.4589034,	715.998717,	707.87435,	727.4863255,	737.984244,	760.8643837,	787.3455427,	791.1650559,	765.8598904,	705.039758,	667.8975652,	609.6473977,	581.9407515,	566.2968668,	556.7446213,	568.0255219,	592.5596927,	641.3621953,	696.5671288,	731.2254129,	765.7689387,	791.0106244,	780.5576145,	738.871454	,703.7412025,	694.3974001	,691.4191767	,731.6054528,	756.4210631	,780.9160513,	787.4793288	,761.8624275,	751.0210896	,703.7914218,	686.5752325	,629.104551	,569.4304585	,555.1185683,	564.725173,	589.2807187], 
     'elprice': [23.86,	22.39,	20.59,	16.81,	17.41,	17.02,	15.86,	18.16,	17.73,	19.77,	23.75,	26.03,	27.06,	26.59, 25,	24.43, 28.87, 37.44, 37.41,	35.34, 33.07, 29.52, 30.1, 24.57, 22.2, 16.57, 15.35, 12.77, 11.27, 11.91, 12.62, 13.83, 16.12,	18.12,	19.59,	21.07,	22.54,	19.79,	16.9,	16.97,	19.25,	28.01	,28.36	,26.56	,17.38	,15.83,	16.97,	15.31,	6.99,	-0.01,	0.23,	4.51,	4.12,	6.76,	0.75,	4.27,	9.11,	15.92,	17.06,	19.28,	20.25	,17.66	,18.49	,16.51,	18.76	,22.21	,26.76	,27.08,	25.33,	22.11,	20.91	,14.43,	13.78,	12.77	,10.56,	3.87,	3.2,	8.67,	18.01,	28.52	,34.74	,33.46	,33.24,	35.09	,35,	35	,34.94,	34.92,	38	,41.21,	42.95,	41.96,	34.94,	30.66	,30,	23.9,	22.34,	17.93,	15.17	,16.38]}
data2 = pd.DataFrame(data=d)

# select periods
periods = len(data2)

# create an energy system
idx = pd.date_range("1/1/2015", periods=periods, freq="H")
es = solph.EnergySystem(timeindex=idx)
Node.registry = es

#fuel cost
ccoal= 12.5
cgas = 20

#Emission pricing
cprice = 8
efgas = 0.2
efcoal = 0.337

# resources gas, coal, waste
bgas = solph.Bus(label="natural gas")
rgas = solph.Source(label="natural gas market", outputs={bgas: solph.Flow(variable_costs = cgas + cprice*efgas)})

bcoal = solph.Bus(label = "hard coal")
rcoal = solph.Source(label = "hard coal market", outputs = {bcoal: solph.Flow(variable_costs = ccoal + cprice*efcoal)})

# heat bus
bheat = solph.Bus(label="heat")

#heat load profile that needs to be served
demand_th = solph.Sink(
    label="demand_th",
    inputs={bheat: solph.Flow(fix=data2['heatload'], nominal_value = 1)},
)

# power bus
bel = solph.Bus(label="electricity")

demand_el = solph.Sink(
    label="demand_el",
    inputs={bel: solph.Flow(variable_costs=-1*data2['elprice'])},
)

# coal CHP
coal_chp = solph.components.GenericCHP(
    label="coal_CHP",
    fuel_input={
        bcoal: solph.Flow(H_L_FG_share_max=[0.1 for p in range(0, periods)],
                          nonconvex=solph.NonConvex(startup_costs=60*31, minimum_downtime=4, initial_status=1),
                          nominal_value=194/0.39,
                          min= 60/194*0.39/0.35,
                          max= 1
                          )
    },
    electrical_output={
        bel: solph.Flow(
            P_max_woDH=[194 for p in range(0, periods)],
            P_min_woDH=[60 for p in range(0, periods)],
            Eta_el_max_woDH=[0.39 for p in range(0, periods)],
            Eta_el_min_woDH=[0.35 for p in range(0, periods)],
        )
    },
    heat_output={bheat: solph.Flow(Q_CW_min=[15 for p in range(0, periods)])},
    Beta=[0.16 for p in range(0, periods)],
    back_pressure=False,
)

# gas-fired heat-only plants
gas_boilers = solph.Transformer(
    label ="Gas boilers",
    inputs = {bgas: solph.Flow(nonconvex=solph.NonConvex(startup_costs=10*70),
                          nominal_value=700/0.88,
                          min= 0.01,
                          max= 1)},
    outputs = {bheat: solph.Flow(max = [700 for p in range(0, periods)], nominal_value=1)},
    conversion_factors = {bheat: 0.88}
)

# create an optimization problem and solve it
om = solph.Model(es)

om = constraints.limit_active_flow_count(om, 'availability', [(bcoal, coal_chp)], lower_limit=0, upper_limit=80)

# solve model
om.solve(solver="gurobi", solve_kwargs={"tee": True})

# create result object
results = solph.processing.results(om)

# plot data
if plt is not None:

    # plot thermal bus
    data = solph.views.node(results, "heat")["sequences"]
    data.columns = ['gas boilers', 'coal CHP', 'Demand']
    ax = data.plot(kind="line", drawstyle="steps-post", grid=True)
    ax.set_xlabel("Time (h)")
    ax.set_ylabel("Q (MW)")
    plt.show()
2

No one got any clue how to achieve this? :confused: I somehow successfully created an availability constraint with the generic_integral_limit but this is not useful when there are more plants in the model. I have no idea, why the limit_active_flow_count is not working and would still be happy about any advice.

3

The limit_active_flow_count function is not right for you. You can pass two flows and limit the number of active flows to one. This means that only one of the flows can be active.

If you want to reduce the usage of the coal plant you can use the annual_limit parameter to reduce the overall resource usage or you can limit the emissions of the coal plant.

If you want to model an unplanned shut down or a revision you can set the capacity to zero. If you want to have an unplanned shut down you can use a random function.

First use a variable within the CHP component e.g. pmax:

- electrical_output={
-     bel: solph.Flow(
-         P_max_woDH=[194 for p in range(0, periods)],
+ electrical_output={
+     bel: solph.Flow(
+         P_max_woDH=pmax,

Then you can manipulate the variable beforehand:

pmax = numpy.empty(periods)
pmax.fill(194)
pmax[30:40].fill(0)

As I said you can use a random function to create shutdown blocks.

I hope this helps.

4

Thank you for your advice. I think the summed_max property of the Flow object is rather what I want. Your proposal also looks fine. I should have looked through all the properties of Flow before as this is quite powerful.