Feature/constraint flow idle

docs/reference/oemof.solph.rst

@@ -113,7 +113,7 @@ oemof.solph.constraints module
 ------------------------------
 
 .. automodule:: oemof.solph.constraints
-    :members: equate_variables, 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, limit_active_flow_count, limit_active_flow_count_by_keyword, emission_limit, generic_integral_limit, additional_investment_flow_limit, investment_limit, shared_limit, set_idle_time
     :undoc-members:
     :show-inheritance:
 

examples/idle_time_example/idle_time_example.py

@@ -0,0 +1,119 @@
+# -*- coding: utf-8 -*-
+
+"""
+
+"""
+
+__copyright__ = "oemof developer group"
+__license__ = "GPLv3"
+
+import os
+import pandas as pd
+import oemof.solph as solph
+
+import matplotlib.pyplot as plt
+
+
+solver = "cbc"
+
+# Create an energy system and optimize the dispatch at least costs.
+# ####################### initialize and provide data #####################
+
+datetimeindex = pd.date_range("1/1/2016", periods=10, freq="H")
+energysystem = solph.EnergySystem(timeindex=datetimeindex)
+
+# ######################### create energysystem components ################
+
+bus = solph.buses.Bus(label="Bus")
+sink = solph.components.Sink(
+    label="Sink",
+    inputs={
+        bus: solph.flows.Flow(
+            nominal_value=1,
+            variable_costs=-1,
+            min=0.2,
+            nonconvex=solph.NonConvex(),
+        )
+    },
+)
+source = solph.components.Source(
+    label="Source",
+    outputs={
+        bus: solph.flows.Flow(
+            nominal_value=1,
+            min=0.2,
+            nonconvex=solph.NonConvex(),
+        )
+    },
+)
+storage = solph.components.GenericStorage(
+    nominal_storage_capacity=10,
+    label="Storage",
+    inputs={bus: solph.Flow(nominal_value=10 / 6)},
+    outputs={bus: solph.Flow(nominal_value=10 / 6, variable_costs=0.001)},
+    initial_storage_level=0,
+    loss_rate=0,
+    inflow_conversion_factor=1,
+    outflow_conversion_factor=1,
+)
+
+energysystem.add(bus, sink, source, storage)
+
+# ################################ optimization ###########################
+
+# create optimization model based on energy_system
+optimization_model = solph.Model(energysystem=energysystem)
+
+solph.constraints.set_idle_time(
+    optimization_model, (source, bus), (bus, sink), n=2
+)
+
+optimization_model.write(
+    "/home/jann/Desktop/lp-idle.lp",
+    io_options={"symbolic_solver_labels": True},
+)
+
+# solve problem
+optimization_model.solve(
+    solver=solver, solve_kwargs={"tee": True, "keepfiles": False}
+)
+
+# write back results from optimization object to energysystem
+optimization_model.results()
+
+# ################################ results ################################
+
+# subset of results that includes all flows into and from electrical bus
+# sequences are stored within a pandas.DataFrames and scalars e.g.
+# investment values within a pandas.Series object.
+# in this case the entry data['scalars'] does not exist since no investment
+# variables are used
+data = solph.views.node(optimization_model.results(), "Bus")
+
+print("Optimization successful. Showing some results:")
+
+# see: https://pandas.pydata.org/pandas-docs/stable/visualization.html
+node_results_bel = solph.views.node(optimization_model.results(), "Bus")
+node_results_flows = node_results_bel["sequences"]
+node_results_flows = node_results_flows[
+    [
+        (("Bus", "Sink"), "flow"),
+        (("Source", "Bus"), "flow"),
+        (("Bus", "Storage"), "flow"),
+        (("Storage", "Bus"), "flow"),
+    ]
+]
+node_results_flows.loc[
+    :, [(("Bus", "Sink"), "flow"), (("Bus", "Storage"), "flow")]
+] *= -1
+fig, ax = plt.subplots(figsize=(10, 5))
+node_results_flows.plot(ax=ax, kind="bar", stacked=True, linewidth=0, width=1)
+ax.set_title("Sums for optimization period")
+ax.legend(loc="upper right", bbox_to_anchor=(1, 1))
+ax.set_xlabel("Energy (MWh)")
+ax.set_ylabel("Flow")
+plt.legend(loc="center left", prop={"size": 8}, bbox_to_anchor=(1, 0.5))
+
+fig.subplots_adjust(right=0.8)
+
+plt.show()

src/oemof/solph/constraints/__init__.py

@@ -11,6 +11,8 @@
 from .investment_limit import additional_investment_flow_limit
 from .investment_limit import investment_limit
 from .shared_limit import shared_limit
+from .set_idle_time import set_idle_time
+
 
 __all__ = [
     "equate_variables",
@@ -21,4 +23,5 @@
     "additional_investment_flow_limit",
     "investment_limit",
     "shared_limit",
+    "set_idle_time",
 ]

src/oemof/solph/constraints/set_idle_time.py

@@ -0,0 +1,90 @@
+# -*- coding: utf-8 -*-
+
+"""Constraints to relate variables in an existing model.
+
+SPDX-FileCopyrightText: Jann Launer
+
+SPDX-License-Identifier: MIT
+"""
+from pyomo import environ as po
+
+
+def set_idle_time(model, f1, f2, n, name_constraint="constraint_idle_time"):
+    r"""
+    Adds a constraint to the given model that enforces f1 to be inactive
+    for n timesteps before f2 can be active.
+
+    For each timestep status of f2 can only be active if f1 has been inactive
+    the previous n timesteps.
+
+    **Constraints:**
+
+    .. math:: X_1(s) + X_2(t) <= 1 \forall 0 < t < n, 0 \le s \le t
+    .. math:: X_1(s) + X_2(t) <= 1 \forall t \ge n, t-n \le s \le t
+
+    Parameters
+    ----------
+    model : oemof.solph.Model
+        Model to which the constraint is added.
+    f1 : tuple
+        First flow tuple.
+    f2 : tuple
+        Second flow tuple. Has to be inactive for a defined number of
+        timesteps after first flow was active.
+    n : int
+        Number of timesteps f2 has to be inactive after f1 has been active.
+    name_constraint : str, default='constraint_idle_time'
+        Name for the equation e.g. in the LP file.
+
+    Returns
+    -------
+    the updated model.
+    """
+    # make sure that idle time is not longer than number of timesteps
+    n_timesteps = len(model.TIMESTEPS)
+    assert n_timesteps > n, (
+        f"Selected idle time {n}"
+        f"is longer than total number of timesteps {n_timesteps}"
+    )
+
+    # Create an index for the idle time
+    model.idle_time = po.RangeSet(0, n)
+
+    def _idle_rule(m):
+        # In the first n steps, the status of f1 has to be inactive
+        # for f2 to be active. For all following timesteps, f1 has to be
+        # inactive in the preceding window of n timesteps for f2 to be active.
+        for ts in list(m.TIMESTEPS):
+            for delay in model.idle_time:
+
+                if ts - delay < 0:
+                    continue
+
+                expr = (
+                    m.NonConvexFlowBlock.status[f2[0], f2[1], ts]
+                    + m.NonConvexFlowBlock.status[f1[0], f1[1], ts - delay]
+                    <= 1
+                )
+                if expr is not True:
+                    getattr(m, name_constraint).add((ts, ts - delay), expr)
+
+    setattr(
+        model,
+        name_constraint,
+        po.Constraint(
+            [
+                (t, t - delay)
+                for t in model.TIMESTEPS
+                for delay in model.idle_time
+            ],
+            noruleinit=True,
+        ),
+    )
+
+    setattr(
+        model,
+        name_constraint + "_build",
+        po.BuildAction(rule=_idle_rule),
+    )
+
+    return model

tests/constraint_tests.py

@@ -886,6 +886,35 @@ def test_equate_variables_constraint(self):
 
         self.compare_lp_files("connect_investment.lp", my_om=om)
 
+    def test_flow_idle(self):
+        bus = solph.buses.Bus(label="Bus")
+        sink = solph.components.Sink(
+            label="Sink",
+            inputs={
+                bus: solph.flows.Flow(
+                    nominal_value=1,
+                    min=0.2,
+                    nonconvex=solph.NonConvex(),
+                )
+            },
+        )
+        source = solph.components.Source(
+            label="Source",
+            outputs={
+                bus: solph.flows.Flow(
+                    nominal_value=1,
+                    min=0.2,
+                    nonconvex=solph.NonConvex(),
+                )
+            },
+        )
+
+        om = self.get_om()
+        solph.constraints.set_idle_time(om, (source, bus), (bus, sink), n=1)
+
+        self.compare_lp_files("set_idle_time.lp", my_om=om)
+
+
     def test_gradient(self):
         """Testing gradient constraints and costs."""
         bel = solph.buses.Bus(label="electricityBus")