diff --git a/docs/reference/oemof.solph.rst b/docs/reference/oemof.solph.rst index c7bbf159a..c35f4261b 100644 --- a/docs/reference/oemof.solph.rst +++ b/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: diff --git a/examples/idle_time_example/idle_time_example.py b/examples/idle_time_example/idle_time_example.py new file mode 100644 index 000000000..85af3f2a4 --- /dev/null +++ b/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() diff --git a/src/oemof/solph/constraints/__init__.py b/src/oemof/solph/constraints/__init__.py index 828aaf70f..965492fd2 100644 --- a/src/oemof/solph/constraints/__init__.py +++ b/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", ] diff --git a/src/oemof/solph/constraints/set_idle_time.py b/src/oemof/solph/constraints/set_idle_time.py new file mode 100644 index 000000000..c9d29a3ca --- /dev/null +++ b/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 diff --git a/tests/constraint_tests.py b/tests/constraint_tests.py index ef231b38e..0c688526d 100644 --- a/tests/constraint_tests.py +++ b/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") diff --git a/tests/lp_files/set_idle_time.lp b/tests/lp_files/set_idle_time.lp new file mode 100644 index 000000000..cc49bb67f --- /dev/null +++ b/tests/lp_files/set_idle_time.lp @@ -0,0 +1,132 @@ +\* Source Pyomo model name=Model *\ + +min +objective: ++0 ONE_VAR_CONSTANT + +s.t. + +c_u_constraint_idle_time(0_0)_: ++1 NonConvexFlowBlock_status(Bus_Sink_0) ++1 NonConvexFlowBlock_status(Source_Bus_0) +<= 1 + +c_u_constraint_idle_time(1_0)_: ++1 NonConvexFlowBlock_status(Bus_Sink_1) ++1 NonConvexFlowBlock_status(Source_Bus_0) +<= 1 + +c_u_constraint_idle_time(1_1)_: ++1 NonConvexFlowBlock_status(Bus_Sink_1) ++1 NonConvexFlowBlock_status(Source_Bus_1) +<= 1 + +c_u_constraint_idle_time(2_1)_: ++1 NonConvexFlowBlock_status(Bus_Sink_2) ++1 NonConvexFlowBlock_status(Source_Bus_1) +<= 1 + +c_u_constraint_idle_time(2_2)_: ++1 NonConvexFlowBlock_status(Bus_Sink_2) ++1 NonConvexFlowBlock_status(Source_Bus_2) +<= 1 + +c_e_BusBlock_balance(Bus_0)_: +-1 flow(Bus_Sink_0) ++1 flow(Source_Bus_0) += 0 + +c_e_BusBlock_balance(Bus_1)_: +-1 flow(Bus_Sink_1) ++1 flow(Source_Bus_1) += 0 + +c_e_BusBlock_balance(Bus_2)_: +-1 flow(Bus_Sink_2) ++1 flow(Source_Bus_2) += 0 + +c_u_NonConvexFlowBlock_min(Bus_Sink_0)_: ++0.20000000000000001 NonConvexFlowBlock_status(Bus_Sink_0) +-1 flow(Bus_Sink_0) +<= 0 + +c_u_NonConvexFlowBlock_min(Bus_Sink_1)_: ++0.20000000000000001 NonConvexFlowBlock_status(Bus_Sink_1) +-1 flow(Bus_Sink_1) +<= 0 + +c_u_NonConvexFlowBlock_min(Bus_Sink_2)_: ++0.20000000000000001 NonConvexFlowBlock_status(Bus_Sink_2) +-1 flow(Bus_Sink_2) +<= 0 + +c_u_NonConvexFlowBlock_min(Source_Bus_0)_: ++0.20000000000000001 NonConvexFlowBlock_status(Source_Bus_0) +-1 flow(Source_Bus_0) +<= 0 + +c_u_NonConvexFlowBlock_min(Source_Bus_1)_: ++0.20000000000000001 NonConvexFlowBlock_status(Source_Bus_1) +-1 flow(Source_Bus_1) +<= 0 + +c_u_NonConvexFlowBlock_min(Source_Bus_2)_: ++0.20000000000000001 NonConvexFlowBlock_status(Source_Bus_2) +-1 flow(Source_Bus_2) +<= 0 + +c_u_NonConvexFlowBlock_max(Bus_Sink_0)_: +-1 NonConvexFlowBlock_status(Bus_Sink_0) ++1 flow(Bus_Sink_0) +<= 0 + +c_u_NonConvexFlowBlock_max(Bus_Sink_1)_: +-1 NonConvexFlowBlock_status(Bus_Sink_1) ++1 flow(Bus_Sink_1) +<= 0 + +c_u_NonConvexFlowBlock_max(Bus_Sink_2)_: +-1 NonConvexFlowBlock_status(Bus_Sink_2) ++1 flow(Bus_Sink_2) +<= 0 + +c_u_NonConvexFlowBlock_max(Source_Bus_0)_: +-1 NonConvexFlowBlock_status(Source_Bus_0) ++1 flow(Source_Bus_0) +<= 0 + +c_u_NonConvexFlowBlock_max(Source_Bus_1)_: +-1 NonConvexFlowBlock_status(Source_Bus_1) ++1 flow(Source_Bus_1) +<= 0 + +c_u_NonConvexFlowBlock_max(Source_Bus_2)_: +-1 NonConvexFlowBlock_status(Source_Bus_2) ++1 flow(Source_Bus_2) +<= 0 + +c_e_ONE_VAR_CONSTANT: +ONE_VAR_CONSTANT = 1.0 + +bounds + 0 <= flow(Bus_Sink_0) <= 1 + 0 <= flow(Bus_Sink_1) <= 1 + 0 <= flow(Bus_Sink_2) <= 1 + 0 <= flow(Source_Bus_0) <= 1 + 0 <= flow(Source_Bus_1) <= 1 + 0 <= flow(Source_Bus_2) <= 1 + 0 <= NonConvexFlowBlock_status(Bus_Sink_0) <= 1 + 0 <= NonConvexFlowBlock_status(Bus_Sink_1) <= 1 + 0 <= NonConvexFlowBlock_status(Bus_Sink_2) <= 1 + 0 <= NonConvexFlowBlock_status(Source_Bus_0) <= 1 + 0 <= NonConvexFlowBlock_status(Source_Bus_1) <= 1 + 0 <= NonConvexFlowBlock_status(Source_Bus_2) <= 1 +binary + NonConvexFlowBlock_status(Bus_Sink_0) + NonConvexFlowBlock_status(Bus_Sink_1) + NonConvexFlowBlock_status(Bus_Sink_2) + NonConvexFlowBlock_status(Source_Bus_0) + NonConvexFlowBlock_status(Source_Bus_1) + NonConvexFlowBlock_status(Source_Bus_2) +end