Merge 354aa0d into 1687fa9

docs/reference/oemof.solph.constraints.rst

@@ -2,6 +2,6 @@ oemof.solph.constraints
 -----------------------
 
 .. automodule:: oemof.solph.constraints
-    :members: equate_variables, equate_flows, limit_active_flow_count, limit_active_flow_count_by_keyword, emission_limit, generic_integral_limit, additional_investment_flow_limit, investment_limit, shared_limit
+    :members: equate_variables, equate_flows, limit_active_flow_count, limit_active_flow_count_by_keyword, emission_limit, generic_integral_limit, additional_investment_flow_limit, investment_limit, shared_limit, storage_level_constraint
     :undoc-members:
     :show-inheritance:

docs/whatsnew/v0-5-1.rst

@@ -25,6 +25,7 @@ New features
       `old_end` holding model-endogenously installed capacity to be decommissioned after its lifetime.
     * Include discounting and calculating annuities in the objective function terms. Introduce attribute `discount_rate`
       of :class:`oemof.solph.models.Model` and `interest_rate` for individual investment objects (options.Investment).
+* Add `storage_level_constraint` that allows to set flows from/to storage (in)active based on storage content.
 
 Documentation
 #############
@@ -52,4 +53,5 @@ Contributors
 
 * Johannes Kochems
 * Tobi Rohrer
+* Patrik Schönfeldt
 

examples/storage_level_constraint/storage_level_constraint.py

@@ -0,0 +1,96 @@
+import pandas as pd
+from oemof.solph import Bus, EnergySystem, Flow, Model
+from oemof.solph.components import GenericStorage, Source, Sink
+from oemof.solph.processing import results
+
+import matplotlib.pyplot as plt
+
+from oemof.solph.constraints import storage_level_constraint
+
+
+es = EnergySystem(
+    timeindex=pd.date_range("2022-01-01", freq="1H", periods=24),
+    infer_last_interval=True,
+)
+
+multiplexer = Bus(
+    label="multiplexer",
+)
+
+storage = GenericStorage(
+    label="storage",
+    nominal_storage_capacity=3,
+    initial_storage_level=1,
+    balanced=True,
+    loss_rate=0.05,
+    inputs={multiplexer: Flow()},
+    outputs={multiplexer: Flow()},
+)
+
+es.add(multiplexer, storage)
+
+in_0 = Source(
+    label="in_0",
+    outputs={multiplexer: Flow(nominal_value=0.5, variable_costs=0.15)},
+)
+es.add(in_0)
+
+in_1 = Source(label="in_1", outputs={multiplexer: Flow(nominal_value=0.1)})
+es.add(in_1)
+
+out_0 = Sink(
+    label="out_0",
+    inputs={multiplexer: Flow(nominal_value=0.25, variable_costs=-0.1)},
+)
+es.add(out_0)
+
+out_1 = Sink(
+    label="out_1",
+    inputs={multiplexer: Flow(nominal_value=0.15, variable_costs=-0.1)},
+)
+es.add(out_1)
+
+
+model = Model(es)
+
+storage_level_constraint(
+    model=model,
+    name="multiplexer",
+    storage_component=storage,
+    multiplexer_bus=multiplexer,
+    input_levels={in_1: 1 / 3},  # in_0 is always active
+    output_levels={out_0: 1 / 6, out_1: 1 / 2},
+)
+model.solve()
+
+my_results = results(model)
+
+df = pd.DataFrame(my_results[(storage, None)]["sequences"])
+df["in1_status"] = my_results[(in_1, None)]["sequences"]
+df["out1_status"] = my_results[(out_1, None)]["sequences"]
+df["out0_status"] = my_results[(out_0, None)]["sequences"]
+
+df["in1"] = my_results[(in_1, multiplexer)]["sequences"]
+df["in0"] = my_results[(in_0, multiplexer)]["sequences"]
+df["out0"] = my_results[(multiplexer, out_0)]["sequences"]
+df["out1"] = my_results[(multiplexer, out_1)]["sequences"]
+
+plt.step(df.index, df["in0"], where="post", label="inflow (<= 1)")
+plt.step(df.index, df["in1"], where="post", label="inflow (< 1/3)")
+plt.step(df.index, df["out0"], where="post", label="outflow (> 1/6)")
+plt.step(df.index, df["out1"], where="post", label="outflow (> 1/2)")
+
+plt.grid()
+plt.legend()
+plt.ylabel("Flow Power (arb. units)")
+
+plt.twinx()
+
+plt.plot(df.index, df["storage_content"], "r-", label="storage content")
+plt.ylim(0, 3)
+plt.legend(loc="center right")
+plt.ylabel("Stored Energy (arb. units)")
+
+print(df)
+
+plt.show()

src/oemof/solph/constraints/__init__.py

@@ -16,6 +16,7 @@
 from .investment_limit import investment_limit
 from .investment_limit import investment_limit_per_period
 from .shared_limit import shared_limit
+from .storage_level import storage_level_constraint
 
 __all__ = [
     "equate_flows",
@@ -31,4 +32,5 @@
     "investment_limit",
     "investment_limit_per_period",
     "shared_limit",
+    "storage_level_constraint",
 ]

src/oemof/solph/constraints/storage_level.py

@@ -0,0 +1,177 @@
+"""A constraint to allow flows from to a storage based on the storage level.
+
+SPDX-FileCopyrightText: Patrik Schönfeldt
+SPDX-FileCopyrightText: Johannes Kochems
+SPDX-FileCopyrightText: Deutsches Zentrum für Luft- und Raumfahrt (DLR)
+
+SPDX-License-Identifier: MIT
+"""
+
+from pyomo import environ as po
+
+
+def storage_level_constraint(
+    model,
+    name,
+    storage_component,
+    multiplexer_bus,
+    input_levels,
+    output_levels,
+):
+    r"""
+    Add constraints to implement storage content based access.
+
+    As a GenericStorage just allows exactly one input and one output,
+    an additional bus, the multiplexer_bus, is used for the connections.
+    Note that all Flow objects connected to the multiplexer_bus have to have
+    a nominal_value.
+
+    Parameters
+    ----------
+    model : oemof.solph.Model
+        Model to which the constraint is added.
+    name : string
+        Name of the multiplexer.
+    storage_component : oemof.solph.components.GenericStorage
+        Storage component whose content should mandate
+        the possible inputs and outputs.
+    multiplexer_bus : oemof.solph.Bus
+        Bus which connects the input and output levels to the storage.
+    input_levels : dictionary with oemof.solph.Bus as keys and float as values
+        Dictionary of buses which act as inputs and corresponding levels
+    output_levels : dictionary with oemof.solph.Bus as keys and float as values
+        Dictionary of buses which act as outputs and corresponding level
+
+    Verbose description can be found in https://arxiv.org/abs/2211.14080
+    """
+
+    def _outputs():
+        OUTPUTS = po.Set(initialize=output_levels.keys())
+        setattr(model, f"{name}_OUTPUTS", OUTPUTS)
+
+        active_output = po.Var(
+            OUTPUTS, model.TIMESTEPS, domain=po.Binary, bounds=(0, 1)
+        )
+        setattr(model, f"{name}_active_output", active_output)
+
+        constraint_name = f"{name}_output_active_constraint"
+
+        def _output_active_rule(m):
+            r"""
+            .. math::
+                y_n \le E(t) / E_n
+            """
+            for t in m.TIMESTEPS:
+                for o in output_levels:
+                    getattr(m, constraint_name).add(
+                        (o, t),
+                        m.GenericStorageBlock.storage_content[
+                            storage_component, t + 1
+                        ]
+                        / storage_component.nominal_storage_capacity
+                        >= active_output[o, t] * output_levels[o],
+                    )
+
+        setattr(
+            model,
+            constraint_name,
+            po.Constraint(
+                OUTPUTS,
+                model.TIMESTEPS,
+                noruleinit=True,
+            ),
+        )
+        setattr(
+            model,
+            constraint_name + "build",
+            po.BuildAction(rule=_output_active_rule),
+        )
+
+        # Define constraints on the output flows
+        def _constraint_output_rule(m, o, p, t):
+            return (
+                m.flow[multiplexer_bus, o, p, t]
+                / m.flows[multiplexer_bus, o].nominal_value
+                <= active_output[o, t]
+            )
+
+        setattr(
+            model,
+            f"{name}_output_constraint",
+            po.Constraint(
+                OUTPUTS,
+                model.PERIODS,
+                model.TIMESTEPS,
+                rule=_constraint_output_rule,
+            ),
+        )
+
+    _outputs()
+
+    def _inputs():
+        INPUTS = po.Set(initialize=input_levels.keys())
+        setattr(model, f"{name}_INPUTS", INPUTS)
+
+        inactive_input = po.Var(
+            INPUTS, model.TIMESTEPS, domain=po.Binary, bounds=(0, 1)
+        )
+        setattr(model, f"{name}_active_input", inactive_input)
+
+        constraint_name = f"{name}_input_active_constraint"
+
+        def _input_active_rule(m):
+            r"""
+            .. math::
+                \hat{y}_n \ge (E(t) - E_n) - E_{max}
+            """
+            for t in m.TIMESTEPS:
+                for o in input_levels:
+                    getattr(m, constraint_name).add(
+                        (o, t),
+                        (
+                            m.GenericStorageBlock.storage_content[
+                                storage_component, t
+                            ]
+                            - input_levels[o]
+                        )
+                        / storage_component.nominal_storage_capacity
+                        <= inactive_input[o, t],
+                    )
+
+        setattr(
+            model,
+            constraint_name,
+            po.Constraint(
+                INPUTS,
+                model.TIMESTEPS,
+                noruleinit=True,
+            ),
+        )
+        setattr(
+            model,
+            constraint_name + "build",
+            po.BuildAction(rule=_input_active_rule),
+        )
+
+        # Define constraints on the input flows
+        def _constraint_input_rule(m, i, p, t):
+            return (
+                m.flow[i, multiplexer_bus, p, t]
+                / m.flows[i, multiplexer_bus].nominal_value
+                <= 1 - inactive_input[i, t]
+            )
+
+        setattr(
+            model,
+            f"{name}_input_constraint",
+            po.Constraint(
+                INPUTS,
+                model.PERIODS,
+                model.TIMESTEPS,
+                rule=_constraint_input_rule,
+            ),
+        )
+
+    _inputs()
+
+    return

tests/constraint_tests.py

@@ -1909,3 +1909,74 @@ def test_nonequidistant_storage(self):
         es.add(b_gas, b_th, boiler, storage)
         om = solph.Model(es)
         self.compare_lp_files("nonequidistant_timeindex.lp", my_om=om)
+
+    def test_storage_level_constraint(self):
+        """Constraint test of an energy system
+        with storage_level_constraint
+        """
+        es = solph.EnergySystem(
+            timeindex=pd.date_range("2022-01-01", freq="1H", periods=2),
+            infer_last_interval=True,
+        )
+
+        multiplexer = solph.Bus(
+            label="multiplexer",
+        )
+
+        storage = solph.components.GenericStorage(
+            label="storage",
+            nominal_storage_capacity=4,
+            initial_storage_level=1,
+            balanced=True,
+            loss_rate=0.25,
+            inputs={multiplexer: solph.Flow()},
+            outputs={multiplexer: solph.Flow()},
+        )
+
+        es.add(multiplexer, storage)
+
+        in_0 = solph.components.Source(
+            label="in_0",
+            outputs={
+                multiplexer: solph.Flow(nominal_value=0.5, variable_costs=0.25)
+            },
+        )
+        es.add(in_0)
+
+        in_1 = solph.components.Source(
+            label="in_1",
+            outputs={multiplexer: solph.Flow(nominal_value=0.125)},
+        )
+        es.add(in_1)
+
+        out_0 = solph.components.Sink(
+            label="out_0",
+            inputs={
+                multiplexer: solph.Flow(
+                    nominal_value=0.25, variable_costs=-0.125
+                )
+            },
+        )
+        es.add(out_0)
+
+        out_1 = solph.components.Sink(
+            label="out_1",
+            inputs={
+                multiplexer: solph.Flow(
+                    nominal_value=0.125, variable_costs=-0.125
+                )
+            },
+        )
+        es.add(out_1)
+
+        om = solph.Model(es)
+
+        solph.constraints.storage_level_constraint(
+            model=om,
+            name="multiplexer",
+            storage_component=storage,
+            multiplexer_bus=multiplexer,
+            input_levels={in_1: 1 / 4},  # in_0 is always active
+            output_levels={out_0: 1 / 8, out_1: 1 / 2},
+        )
+        self.compare_lp_files("storage_level_constraint.lp", my_om=om)