simple_dispatch

Complexity

Total Complexity

4

Size/Duplication

Total Lines	278
Duplicated Lines	0 %

Importance

Changes

0

1 Function

Rating	Name	Duplication	Size	Complexity
B	main()	0	224	4

# -*- coding: utf-8 -*-

"""
General description
-------------------
This example shows how to create an energysystem with oemof objects and
solve it with the solph module. Results are plotted with solph.

Dispatch modelling is a typical thing to do with solph. However cost does not
have to be monetary but can be emissions etc. In this example a least cost
dispatch of different generators that meet an inelastic demand is undertaken.
Some of the generators are renewable energies with marginal costs of zero.
Additionally, it shows how combined heat and power units may be easily modelled
as well.

Data
----
input_data.csv

Installation requirements
-------------------------

This example requires oemof.solph (v0.5.x), install by:

    pip install oemof.solph[examples]


License
-------
`MIT license <https://github.com/oemof/oemof-solph/blob/dev/LICENSE>`_
"""

import os
import warnings

import matplotlib.pyplot as plt
import pandas as pd

from oemof.solph import Bus
from oemof.solph import EnergySystem
from oemof.solph import Flow
from oemof.solph import Model
from oemof.solph import create_time_index
from oemof.solph import views
from oemof.solph.components import Sink
from oemof.solph.components import Source
from oemof.solph.components import Transformer


def main():
    # Read data file
    filename = os.path.join(os.getcwd(), "input_data.csv")
    try:
        data = pd.read_csv(filename)
    except FileNotFoundError:
        msg = "Data file not found: {0}. Only one value used!"
        warnings.warn(msg.format(filename), UserWarning)
        data = pd.DataFrame(
            {
                "pv": [0.3, 0.7],
                "wind": [0.6, 0.5],
                "demand_el": [500, 400],
                "demand_th": [400, 300],
            }
        )

    solver = "cbc"

    # Create an energy system and optimize the dispatch at least costs.
    # ####################### initialize and provide data #####################
    datetimeindex = create_time_index(2016, number=len(data))
    energysystem = EnergySystem(
        timeindex=datetimeindex, infer_last_interval=False
    )

    # ######################### create energysystem components ################

    # resource buses
    bcoal = Bus(label="coal", balanced=False)
    bgas = Bus(label="gas", balanced=False)
    boil = Bus(label="oil", balanced=False)
    blig = Bus(label="lignite", balanced=False)

    # electricity and heat
    bel = Bus(label="bel")
    bth = Bus(label="bth")

    energysystem.add(bcoal, bgas, boil, blig, bel, bth)

    # an excess and a shortage variable can help to avoid infeasible problems
    energysystem.add(Sink(label="excess_el", inputs={bel: Flow()}))
    # shortage_el = Source(label='shortage_el',
    #                      outputs={bel: Flow(variable_costs=200)})

    # sources
    energysystem.add(
        Source(
            label="wind",
            outputs={bel: Flow(fix=data["wind"], nominal_value=66.3)},
        )
    )

    energysystem.add(
        Source(
            label="pv", outputs={bel: Flow(fix=data["pv"], nominal_value=65.3)}
        )
    )

    # demands (electricity/heat)
    energysystem.add(
        Sink(
            label="demand_el",
            inputs={bel: Flow(nominal_value=85, fix=data["demand_el"])},
        )
    )

    energysystem.add(
        Sink(
            label="demand_th",
            inputs={bth: Flow(nominal_value=40, fix=data["demand_th"])},
        )
    )

    # power plants
    energysystem.add(
        Transformer(
            label="pp_coal",
            inputs={bcoal: Flow()},
            outputs={bel: Flow(nominal_value=20.2, variable_costs=25)},
            conversion_factors={bel: 0.39},
        )
    )

    energysystem.add(
        Transformer(
            label="pp_lig",
            inputs={blig: Flow()},
            outputs={bel: Flow(nominal_value=11.8, variable_costs=19)},
            conversion_factors={bel: 0.41},
        )
    )

    energysystem.add(
        Transformer(
            label="pp_gas",
            inputs={bgas: Flow()},
            outputs={bel: Flow(nominal_value=41, variable_costs=40)},
            conversion_factors={bel: 0.50},
        )
    )

    energysystem.add(
        Transformer(
            label="pp_oil",
            inputs={boil: Flow()},
            outputs={bel: Flow(nominal_value=5, variable_costs=50)},
            conversion_factors={bel: 0.28},
        )
    )

    # combined heat and power plant (chp)
    energysystem.add(
        Transformer(
            label="pp_chp",
            inputs={bgas: Flow()},
            outputs={
                bel: Flow(nominal_value=30, variable_costs=42),
                bth: Flow(nominal_value=40),
            },
            conversion_factors={bel: 0.3, bth: 0.4},
        )
    )

    # heat pump with a coefficient of performance (COP) of 3
    b_heat_source = Bus(label="b_heat_source")
    energysystem.add(b_heat_source)

    energysystem.add(
        Source(label="heat_source", outputs={b_heat_source: Flow()})
    )

    cop = 3
    energysystem.add(
        Transformer(
            label="heat_pump",
            inputs={bel: Flow(), b_heat_source: Flow()},
            outputs={bth: Flow(nominal_value=10)},
            conversion_factors={bel: 1 / 3, b_heat_source: (cop - 1) / cop},
        )
    )

    # ################################ optimization ###########################

    # create optimization model based on energy_system
    optimization_model = Model(energysystem=energysystem)

    # 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 = views.node(optimization_model.results(), "bel")
    data["sequences"].info()
    print("Optimization successful. Showing some results:")

    # see: https://pandas.pydata.org/pandas-docs/stable/visualization.html
    node_results_bel = views.node(optimization_model.results(), "bel")
    node_results_flows = node_results_bel["sequences"]

    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)

    dates = node_results_flows.index
    tick_distance = int(len(dates) / 7) - 1
    ax.set_xticks(range(0, len(dates), tick_distance), minor=False)
    if tick_distance > 0:
        ax.set_xticklabels(
            [
                item.strftime("%d-%m-%Y")
                for item in dates.tolist()[0::tick_distance]
            ],
            rotation=90,
            minor=False,
        )

    plt.show()

    node_results_bth = views.node(optimization_model.results(), "bth")
    node_results_flows = node_results_bth["sequences"]

    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)

    dates = node_results_flows.index
    tick_distance = int(len(dates) / 7) - 1
    ax.set_xticks(range(0, len(dates), tick_distance), minor=False)
    if tick_distance > 0:
        ax.set_xticklabels(
            [
                item.strftime("%d-%m-%Y")
                for item in dates.tolist()[0::tick_distance]
            ],
            rotation=90,
            minor=False,
        )

    plt.show()


if __name__ == "__main__":
    main()