Storage inflow and outflow depending on SOC

Hello,
I want to ask if anyone knows how to implement in oemof.solph the following problem:
I have Energysystem which inculdes a storage.
I have limeted the maximum input and outpu flow of the storage.
Now i realised that the maximum input and output flow depends on the state of charge of the storage.
E.g.: Wenn the storage is empty I can load it with the full maximum input flow. But when the storage is full at 90% I am only able to load it with 0.2*maximum input flow.
Does any one know if this is possible to implement in oemof and if yes how?
THX for any help
Yours Thomas

Hi @TomK,

Doing this will need you to either change the class GenericStorage or to manually add a constraint after the model is build. If the relation between SOC and the flow limit is linear the second option is most probably the easier one. The constraint could look like oemof.solph.constraints.equate_variables (documentation of oemof.solph 0.5.0 at oemof-solph.readthedocs.io) but using “<=” instead of “==”. So, the question is whether something in the style of “Flow[grid, battery] <= (1 - SOC) * P_max” is good enough to describe your problem. If not, additional (integer) variables are needed.

Cheers,
Patrik

Hi @pschoen ,

THX for the information!
I will try it.

Do you know if it is possible to use a “if” in the constraints?
So I can say between a SOC of 0% and 80% this linear function should be used and between a SOC of 80% and 100% the other linear function should be used.

Yours
Thomas

Hi @TomK,

The description is purely based on declarative relations (such as equations), so unfortunately, no “if” is possible. However, there is a solution that uses a binary variable. End of last year, I commented on it in arXiv:2211.14080. Equations 3 formulates the “if” condition and Eq. 4 reformulates it using binary variables.

Cheers,
Patrik

Hi @pschoen,

Hmmm…
I try to clearfy my problem a little better.

I have at the moment my storage as followed:

el_storage = solph.components.GenericStorage(
    nominal_storage_capacity=3000,
    label="electrical_storage",
    inputs={bus_el: solph.Flow(nominal_value=2400)},
    outputs={bus_el: solph.Flow(nominal_value=2400)},
    loss_rate=0.0,
    initial_storage_level=0.5,
    inflow_conversion_factor=0.965,
    outflow_conversion_factor=0.965,
    )

So when i interpert this in the right way i will get as solution space for the allowed input flow the blue area in the image.
Now i want to inculde the linear function which represents the red line. So i want to limit the solution space to the red area.

So can i do this in they way you showed before with the constraints.equate_variables?

Yours
Thomas

Thank you for the visual representation. I got your problem before but having it graphic helped me to find another way that works without “if” (or binary variables). You can add a constraint that limits charging to 5 * nominal_value * (1 - soc) in the same way constraints.equate_variables is implemented. Below a SOC of 80 %, it will not have any effect as the regular limit is lower. Above, it will automatically be the only relevant constraint.

import pandas as pd
import numpy as np
import logging
import os
import pprint as pp
import matplotlib.pyplot as plt
from pyomo.environ import Constraint


from oemof import solph
from pyomo import environ as po

def equate_variables(model, var1, var2, var3, factor1=1, name=None):

    if name is None:
        name = "_".join(["equate", str(var1), str(var2),str(var3)])

    def equate_variables_rule(m):
        return var3 <= 5*var2*(1-var1)

    setattr(model, name, po.Constraint(rule=equate_variables_rule))


def run_oemof():
    my_index = pd.date_range(start='04/08/2019', periods=672, freq='15min')
    my_energysystem = solph.EnergySystem(timeindex=my_index, infer_last_interval=True)

    # csv Datei einlesen
    csvFile = pd.read_csv(r"..\Data\Input_Data_Dymola\Input_Dataframe_FMU.txt",
                          sep=";", decimal=".", thousands=",", header=None)

    # Definition der Buse des Energiesystems und hinzufügen zum Energiesystem
    bus_el = solph.Bus(label="electricity")
    bus_gas = solph.Bus(label="gas")
    bus_heat = solph.Bus(label="heat")
    bus_cool = solph.Bus(label="cool")

    my_energysystem.add(bus_el)
    my_energysystem.add(bus_gas)
    my_energysystem.add(bus_heat)
    my_energysystem.add(bus_cool)

    # Definition der Sources
    my_energysystem.add(solph.components.Source(label='el_Netz', outputs={bus_el: solph.Flow(
        variable_costs=csvFile[13]+0.04406,investment=solph.Investment(ep_costs=51.998))}))

    my_energysystem.add(solph.components.Source(label='el_PV', outputs={bus_el: solph.Flow(
        fix=csvFile[4], nominal_value=4)}))

    my_energysystem.add(solph.components.Source(label='gas_Netz', outputs={bus_gas: solph.Flow(
        variable_costs=0.027)}))

    # Definition der Sinks
    my_energysystem.add(solph.components.Sink(label='electricity_demand', inputs={bus_el: solph.Flow(
        fix=csvFile[1], nominal_value=4)}))

    my_energysystem.add(solph.components.Sink(label='heat_demand', inputs={bus_heat: solph.Flow(
        fix=csvFile[2], nominal_value=4)}))

    my_energysystem.add(solph.components.Sink(label='cool_demand', inputs={bus_cool: solph.Flow(
        fix=csvFile[3], nominal_value=4)}))

    # create excess components to allow overproduction
    my_energysystem.add(solph.components.Sink(label="excess_bus_el", inputs={bus_el: solph.Flow()}))
    my_energysystem.add(solph.components.Sink(label="excess_bus_heat", inputs={bus_heat: solph.Flow()}))

    # Definition der Transformers

    my_energysystem.add(solph.components.Transformer(
        label="Chiller",
        inputs={bus_el: solph.Flow()},
        outputs={bus_cool: solph.Flow(nominal_value=4000),
                 bus_heat: solph.Flow(nominal_value=4000)},
        conversion_factors={bus_cool: 2.8, bus_heat: 2.85}))

    my_energysystem.add(solph.components.Transformer(
        label="gas_boiler",
        inputs={bus_gas: solph.Flow()},
        outputs={bus_heat: solph.Flow(nominal_value=10e10)},
        conversion_factors={bus_heat: 0.9}))

    my_energysystem.add(solph.components.Transformer(
        label="heat_pump",
        inputs={bus_el: solph.Flow(nominal_value=40)},
        outputs={bus_heat: solph.Flow()},
        conversion_factors={bus_heat: 4.971538}))

    # Definition der Storages and add them to esys

    el_storage = solph.components.GenericStorage(
        nominal_storage_capacity=3000,
        label="electrical_storage",
        inputs={bus_el: solph.Flow(nominal_value=2400)},
        outputs={bus_el: solph.Flow(nominal_value=2400)},
        loss_rate=0.0,
        initial_storage_level=0.5,
        inflow_conversion_factor=0.965,
        outflow_conversion_factor=0.965,
        )

    heat_storage = solph.components.GenericStorage(
        nominal_storage_capacity=4000,
        label="heat_storage",
        inputs={bus_heat: solph.Flow(nominal_value=1396.67)},
        outputs={bus_heat: solph.Flow(nominal_value=1396.67)},
        loss_rate=0.0001,
        initial_storage_level=0.5,
        inflow_conversion_factor=0.99999,
        outflow_conversion_factor=0.99999,
        )

    my_energysystem.add(el_storage)
    my_energysystem.add(heat_storage)

    om = solph.Model(my_energysystem)

#------------Adaption------------------------------------------
    soc = 1010 # state of charge
    nominal_value = 2400
    flow = 1010 # electrical_storage input flow from el_bus
    solph.constraints.equate_variables(
    om,
    soc,
    nominal_value,
    flow)
#---------------------------------------------------------------

    om.solve(solver='gurobi', solve_kwargs={'tee': True})


    # Ausgabe für Berechnung
    results = solph.processing.results(om)
    oemof_result=pd.DataFrame()

    for ob in results.keys():
        new_df = results[ob]['sequences']
        label = 'VON: ' + ob[0].label+' [kWh]'
        faktor=1
        if hasattr(ob[1], 'label'):
            label = 'VON: ' + ob[0].label + '   ' + 'ZU: ' + ob[1].label+' [kWh]'
            faktor=4
        new_df = new_df.rename(columns={new_df.columns[0]: label})/faktor
        oemof_result = pd.concat([oemof_result, new_df], axis=1)

    return oemof_result, my_energysystem

result_tuple=run_oemof()

Hi @pschoen ,

sorry for the late answer.

I tried now to implement in my code what we have discussed before.
I copied the code into here as you can see above.

I am stil struggeling with some points:

  1. I don’t know how to adrers the state of charge. It is changed during the optimisation and I don’t know how to get it before the optimisation starts.

  2. I am also not shure how to select the input flow from the electrical storage.

  3. I hope the formular is working in that way

Do you have any suggestions how to solve it?
The changed parts are in “adaption” and the “def equate_variables”.

THX
Thomas

Hi @TomK,

The additional constraints can be added after the model is created. You already did that. The variable for the storage content actually has two dimensions, one index is for the storage component (the GenericStorage object), one is for the time step t in m.TIMESTEPS. The equation in the constraint should look something like this in your case:

1/ full_charging_limit * nominal_value * (
      1 - m.GenericStorageBlock.storage_content[storage_component, t + 1] / storage_capacity
) >= m.flow[i, o, t]

The parameter full_charging_limit gives the share that cannot be charged at the full rate. Note that I put “+ 1” to the time, so that charging slows down if the limit is passed within the time step.

Cheers,
Patrik

Hi @pschoen,

I tried to insert your formular. I got the following Error:

ERROR: Rule failed when generating expression for Constraint equate with index
0: TypeError: equate_variables2..equate_variables_rule() missing 1
required positional argument: ‘t’

I am still not shure how the program knows that model.flow[i,o,t] is the input flow from the electrical storage. Because “model” includes the whole energy model. As “storage_component” I insert the name of the electrical storage. Here I am not shure if I have to use “el_storage” or the lable “electrical_storage”.

THX
Thomas

Here is the updated code:

import pandas as pd
import numpy as np
import logging
import os
import pprint as pp
import matplotlib.pyplot as plt
from pyomo.environ import Constraint


from oemof import solph
from pyomo import environ as po

def equate_variables2(model, name=None):

    if name is None:
        name = "_".join(["equate"])

    full_charging_limit = 0.8
    nominal_value = 2400
    storage_capacity = 3000


    def equate_variables_rule(i,o,t):
        return 1/ full_charging_limit * nominal_value * (1 - model.GenericStorageBlock.storage_content["el_storage", t + 1] / storage_capacity) >= model.flow[i, o, t]

    setattr(model, name, po.Constraint(model.TIMESTEPS, rule=equate_variables_rule))


def run_oemof():
    my_index = pd.date_range(start='04/08/2019', periods=672, freq='15min')
    my_energysystem = solph.EnergySystem(timeindex=my_index, infer_last_interval=True)

    # csv Datei einlesen
    csvFile = pd.read_csv(r"..\Data\Input_Data_Dymola\Input_Dataframe_FMU.txt",
                          sep=";", decimal=".", thousands=",", header=None)

    # Definition der Buse des Energiesystems und hinzufügen zum Energiesystem
    bus_el = solph.Bus(label="electricity")
    bus_gas = solph.Bus(label="gas")
    bus_heat = solph.Bus(label="heat")
    bus_cool = solph.Bus(label="cool")

    my_energysystem.add(bus_el)
    my_energysystem.add(bus_gas)
    my_energysystem.add(bus_heat)
    my_energysystem.add(bus_cool)

    # Definition der Sources
    my_energysystem.add(solph.components.Source(label='el_Netz', outputs={bus_el: solph.Flow(
        variable_costs=csvFile[13]+0.04406,investment=solph.Investment(ep_costs=51.998))}))

    my_energysystem.add(solph.components.Source(label='el_PV', outputs={bus_el: solph.Flow(
        fix=csvFile[4], nominal_value=4)}))

    my_energysystem.add(solph.components.Source(label='gas_Netz', outputs={bus_gas: solph.Flow(
        variable_costs=0.027)}))

    # Definition der Sinks
    my_energysystem.add(solph.components.Sink(label='electricity_demand', inputs={bus_el: solph.Flow(
        fix=csvFile[1], nominal_value=4)}))

    my_energysystem.add(solph.components.Sink(label='heat_demand', inputs={bus_heat: solph.Flow(
        fix=csvFile[2], nominal_value=4)}))

    my_energysystem.add(solph.components.Sink(label='cool_demand', inputs={bus_cool: solph.Flow(
        fix=csvFile[3], nominal_value=4)}))

    # create excess components to allow overproduction
    my_energysystem.add(solph.components.Sink(label="excess_bus_el", inputs={bus_el: solph.Flow()}))
    my_energysystem.add(solph.components.Sink(label="excess_bus_heat", inputs={bus_heat: solph.Flow()}))

    # Definition der Transformers

    my_energysystem.add(solph.components.Transformer(
        label="Chiller",
        inputs={bus_el: solph.Flow()},
        outputs={bus_cool: solph.Flow(nominal_value=4000),
                 bus_heat: solph.Flow(nominal_value=4000)},
        conversion_factors={bus_cool: 2.8, bus_heat: 2.85}))

    my_energysystem.add(solph.components.Transformer(
        label="gas_boiler",
        inputs={bus_gas: solph.Flow()},
        outputs={bus_heat: solph.Flow(nominal_value=10e10)},
        conversion_factors={bus_heat: 0.9}))

    my_energysystem.add(solph.components.Transformer(
        label="heat_pump",
        inputs={bus_el: solph.Flow(nominal_value=40)},
        outputs={bus_heat: solph.Flow()},
        conversion_factors={bus_heat: 4.971538}))

    # Definition der Storages and add them to esys

    el_storage = solph.components.GenericStorage(
        nominal_storage_capacity=3000,
        label="electrical_storage",
        inputs={bus_el: solph.Flow(nominal_value=2400)},
        outputs={bus_el: solph.Flow(nominal_value=2400)},
        loss_rate=0.0,
        initial_storage_level=0.5,
        inflow_conversion_factor=0.965,
        outflow_conversion_factor=0.965,
        )

    heat_storage = solph.components.GenericStorage(
        nominal_storage_capacity=4000,
        label="heat_storage",
        inputs={bus_heat: solph.Flow(nominal_value=1396.67)},
        outputs={bus_heat: solph.Flow(nominal_value=1396.67)},
        loss_rate=0.0001,
        initial_storage_level=0.5,
        inflow_conversion_factor=0.99999,
        outflow_conversion_factor=0.99999,
        )

    my_energysystem.add(el_storage)
    my_energysystem.add(heat_storage)

    om = solph.Model(my_energysystem)

#------------Adaption------------------------------------------
    equate_variables2(
    om)
#---------------------------------------------------------------

    om.solve(solver='gurobi', solve_kwargs={'tee': True})


    # Ausgabe für Berechnung
    results = solph.processing.results(om)
    oemof_result=pd.DataFrame()

    for ob in results.keys():
        new_df = results[ob]['sequences']
        label = 'VON: ' + ob[0].label+' [kWh]'
        faktor=1
        if hasattr(ob[1], 'label'):
            label = 'VON: ' + ob[0].label + '   ' + 'ZU: ' + ob[1].label+' [kWh]'
            faktor=4
        new_df = new_df.rename(columns={new_df.columns[0]: label})/faktor
        oemof_result = pd.concat([oemof_result, new_df], axis=1)

    return oemof_result, my_energysystem

result_tuple=run_oemof()

I was a bit sloppy here, in particular because I also don’t always do everything correct instantly and often have to try myself.

m.flow[i, o, t] was the standard formulation, denoting the flow from i to o at point t. In your case, it would translate to:

    def equate_variables_rule(t):
        return 1/ full_charging_limit * nominal_value * (1 - model.GenericStorageBlock.storage_content[el_storage, t + 1] / storage_capacity) >= model.flow[bus_el, el_storage, t]

For the index, it should be the instance of the object and not the string representation of the label. (Both might work, but the string is definitely not the commonly used option.)

Doesn’t matter. I am happy that someone can help me.

I tried it know in this way:

return 1 / full_charging_limit * nominal_value * (
                    1 - model.GenericStorageBlock.storage_content['el_storage', t + 1] / storage_capacity) >= \
               model.flow['bus_el', 'el_storage', t]

But I get the error message:

File "C:/Users/TKurz/Documents/DISS/GitHub/Dymola_Neural_Network/Scripts/oemof_model.py", line 31, in equate_variables_rule
    1 - model.GenericStorageBlock.storage_content['el_storage', t + 1] / storage_capacity) >= \
TypeError: unsupported operand type(s) for +: 'Model' and 'int'

Therefore i also tried this:

        return 1 / full_charging_limit * nominal_value * (
                1 - model.GenericStorageBlock.storage_content['electrical_storage', t + 1] / storage_capacity) >= \
           model.flow['electricity', 'electrical_storage', t]

But I got the same error message.

I also tryed to insert instead of ‘el_storage’ something like model.el_storage or model.xxx.storage but I was only able to receive the error that this component dosen’t exist.

Sorry that I have so many questions :grimacing:

THX
Thomas

For some reason, t seems to be your model. I think, your code should contain something like the following:

my_energysystem.add(el_storage, bus_el)

model = solph.Model(my_energysystem)

def soc_limit_rule(t):
    return (
        1/ full_charging_limit * nominal_value * (
            1 - model.GenericStorageBlock.storage_content[el_storage, t + 1]
            / storage_capacity)
        >= model.flow[bus_el, el_storage, t]
    )
setattr(
    model,
    "soc_limit",
    po.Constraint(
        model.TIMESTEPS,
        noruleinit=True,
    ),
)
setattr(
    model,
    "soc_limit_build",
    po.BuildAction(rule=soc_limit_rule),
)

(Did not test it, but maybe it helps you.)

PS: I just added a pretty similar constraint to solph, that allows for stepwise limits: Add storage_level_constraint by p-snft · Pull Request #936 · oemof/oemof-solph (github.com). Maybe, also that one can serve as a basis for your own constraint.