Creating energy systems from spreadsheet¶

[1]:

import os
import pkg_resources as pkg
import pandas as pd

from pyomo.opt import SolverFactory
from oemof.solph import EnergySystem, Model, Bus
from oemof.tools.economics import annuity as annuity
from oemof.solph import constraints
import oemof.tabular.tools.postprocessing as pp
import oemof.tabular.facades as fc

Creating and Setting the Datapaths¶

Setting the datapath for raw-data and results. Data handling looks more complex than it is. You can easily adapt this to a simple pd.read_excel(filepath,...) in the next block if your file is located somewhere else. Otherwise we will use data from the oemof tabular repository.

In addition a results directory will be created in home/user/oemof-results/dispatch/output.

[2]:

scenario_name = "base-scenario"

# datapath for input data from the oemof tabular pacakge
datapath = os.path.join(
    pkg.resource_filename("oemof.tabular", ""),
    "data/data.xls",
)

# results path
results_path = os.path.join(
    os.path.expanduser("~"), "oemof-results"
)

scenario_path = os.path.join(
    results_path, scenario_name, "output"
)

if not os.path.exists(scenario_path):
    os.makedirs(scenario_path)
print(scenario_path)

/home/docs/oemof-results/base-scenario/output

Next we will read the required input data. The profiles index will be used for the EnergySystem object below. All generator data etc. will also be loaded.

[3]:

file = pd.ExcelFile(datapath)

sheet_names = [typ for typ in file.sheet_names if typ in fc.TYPEMAP.keys()]

data = {}

for sheet in sheet_names:
    data[sheet]= pd.read_excel(datapath, sheet_name=sheet, index_col=0)

# profiles and tech data not be part of TYPEMAP and need to be read seperately
profiles = pd.read_excel(
    datapath,
    sheet_name="profiles",
    index_col=[0],
    parse_dates=True,
)
profiles.index.freq = "1H"

technology = pd.read_excel(
    datapath, sheet_name="technology-data", index_col=[0, 1]
)
carrier = pd.read_excel(
    datapath, sheet_name="carrier", index_col=[0]
)

[4]:

all_components = pd.concat([v for k, v in data.items() if k != "bus"], sort=False)
# Only be used for Latex export of tables
#columns = ['profile', 'capacity_potential']
#print(all_components.to_latex(columns=columns, na_rep="-"))

Creating the EnergySystem and its Nodes¶

We are starting by creating a EnergySystem object which will hold all information (nodes, etc.) of hour energy system that we will add below. This is just the standard way of using the oemof.solph library for your modelling.

[5]:

es = EnergySystem(timeindex=profiles.index)

Add Bus¶

Before we add any components we will create all bus objects for our model and add it to the energy system object.

[6]:

buses = {
    name: Bus(label=name, balanced=bool(arg.balanced))
    for name, arg in data["bus"].iterrows()
}
es.add(*buses.values())

Bus Constraints¶

With the set of all Buses \(B\) all inputs \(x^{flow}_{i(b),b}\) to a bus \(b\) must equal all its outputs \(x^{flow}_{b,o(b)}\)

\[\sum_i x^{flow}_{i(b), b}(t) - \sum_o x^{flow}_{b, o(b)}(t) = 0 \qquad \forall t \in T, \forall b \in B\]

This equation will be build once the complete energy system is setup with its component. Every time a Component is created, the connected bus inputs/outputs will be updated. By this update every bus has all required information of its inputs and outputs available to construct the constraints.

Add Load¶

[7]:

for name, l in data["load"].iterrows():
    es.add(
        fc.Load(
            label=name,
            bus=buses[
                l.bus
            ],  # reference the bus in the buses dictionary
            amount=l.amount,  # amount column
            profile=profiles[l.profile],
        )
    )

Load Constraint¶

For the set of all Load denoted with \(l \in L\) the load \(x_l\) at timestep t equals the exogenously defined profile value \(c^{profile}_l\) multiplied by the amount of this load \(c^{amount}_l\)

\[x^{flow}_{l}(t) = c^{profile}_{l}(t) \cdot c^{amount}_{l} \qquad \forall t \in T, \forall l \in L\]

Add Generators¶

[8]:

for name, g in data["dispatchable"].iterrows():
    es.add(
        fc.Dispatchable(
            label=name,
            bus=buses[g.bus],
            carrier=g.carrier,
            tech=g.tech,
            marginal_cost=(
                carrier.at[g.carrier, "cost"]
                / technology.at[
                    (g.carrier, g.tech), "efficiency"
                ]
            ),
            # efficiency=technology.at[(g.carrier, g.tech), 'efficiency'],
            expandable=g.expandable,
            capacity=g.capacity,
            capacity_potential=g.capacity_potential,
            capacity_cost=annuity(
                technology.at[
                    (g.carrier, g.tech), "capex"
                ],  # to $/MW
                technology.at[
                    (g.carrier, g.tech), "lifetime"
                ],
                0.07,
            )
            * 1000,
            output_parameters={
                "emission_factor": (
                    carrier.at[g.carrier, "emission_factor"]
                    / technology.at[
                        (g.carrier, g.tech), "efficiency"
                    ]
                )
            },
        )
    )

Dispatchable Generator Constraint¶

A Generator component can be used to model all types of dispatchble units in a energy system. This can include diesel generators oder coal fired power plants but also hot water boilers for heat. Every generator must be connected to an Bus object.

This basic mathematical model for the component with the set of all dispatchable generators being \(d \in D\) looks as follows:

\[x^{flow}_{d}(t) \leq x^{capacity}_{d} \qquad \forall t \in T, \forall d \in D\]

Meaning, the production of the generator \(x^{flow}\) must be less than its maximum capacity \(c^{max}\) in every timestep.

[9]:

for name, v in data["volatile"].iterrows():
    es.add(
        fc.Volatile(
            label=name,
            bus=buses[v.bus],
            carrier=v.carrier,
            tech=v.tech,
            expandable=v.expandable,
            capacity=v.capacity,
            capacity_potential=v.capacity_potential,
            capacity_cost=annuity(
                technology.at[(v.carrier, v.tech), "capex"],
                technology.at[
                    (v.carrier, v.tech), "lifetime"
                ],
                0.07,
            )
            * 1000,  # $/kW -> $/MW
            profile=profiles[v.profile],
        )
    )

Volatile Generator Constraint¶

Using the Generator component with output_parameters={"fixed": True} is very similar to the Dispatchable component. However, in this case the flow of the volatile components denoted with \(v \in V\) will be fixed to a specific value.

\[x^{flow}_{v}(t) = c^{profile}_{v}(t) \cdot x^{capacity}_{v} \qquad \forall t \in T, \forall v \in V\]

Alternatively you can use the Volatile component which automatically enforced the fixed behaviour.

Add Storage¶

[10]:

for name, s in data["storage"].iterrows():
    es.add(
        fc.Storage(
            label=name,
            bus=buses[s.bus],
            carrier=s.carrier,
            tech=s.tech,
            marginal_cost=s.marginal_cost,
            capacity=s.capacity,
            storage_capacity=s.storage_capacity,
            expandable=s.expandable,
            efficiency=technology.at[
                (s.carrier, s.tech), "efficiency"
            ],
            loss_rate=s.loss_rate,
            storage_capacity_cost=annuity(
                technology.at[
                    (s.carrier, s.tech), "storage_capex"
                ],
                technology.at[
                    (s.carrier, s.tech), "lifetime"
                ],
                0.07,
            )
            * 1000,  # $/kW -> $/MW
            capacity_cost=annuity(
                technology.at[(s.carrier, s.tech), "capex"],
                technology.at[
                    (s.carrier, s.tech), "lifetime"
                ],
                0.07,
            )
            * 1000,  # $/kW -> $/MW
        )
    )

Storage Constraints¶

The mathematical representation of the storage for all storages \(s \in S\) will include the flow into the storage, out of the storage and a storage level. The defaul efficiency for input/output is 1. Note that this is is included during charge and discharge. If you want to set the round trip efficiency you need to do for example: \(\eta = \sqrt{\eta^{roundtrip}}\)

Intertemporal energy balance of the storage:

\[x^{level}_{s}(t) = \eta^{loss} x^{level}_{s}(t) + \eta^{in} x^{flow}_{s, in} - \frac{x^{flow}_{s, out}(t)}{\eta^{out}} \qquad \forall t \in T, \forall s \in S\]

Bounds of the storage level variable \(x^{level}_s(t)\):

\[x^{level}_s(t) \leq c_s^{max,level} \qquad \forall t \in T, \forall s \in S\]

\[x^{level}_s(1) = x_s^{level}(t_{e}) = 0.5 \cdot c_s^{max,level} \qquad \forall t \in T, \forall s \in S\]

Of course, in addition the inflow/outflow of the storage also needs to be within the limit of the minimum and maximum power.

\[-c_s^{capacity} \leq x^{flow}_s(t) \leq c_s^{capacity} \qquad \forall t \in T, \forall s \in S\]

Add Conversion¶

A conversion unit will take from a bus and feed into another:

\[x^{flow}_{c, to}(t) = c^{efficiencty}_{c} \cdot x^{flow}_{c, from}(t), \qquad \forall c \in C, \forall t \in T\]

[11]:

for name, c in data["conversion"].iterrows():
    es.add(
        fc.Conversion(
            label=name,
            from_bus=buses[c.from_bus],
            to_bus=buses[c.to_bus],
            carrier=c.carrier,
            tech=c.tech,
            efficiency=technology.at[
                (c.carrier, c.tech), "efficiency"
            ],
            marginal_cost=(
                carrier.at[c.carrier, "cost"]
                / technology.at[
                    (c.carrier, c.tech), "efficiency"
                ]
            ),
            expandable=c.expandable,
            capacity=c.capacity,
            capacity_potential=c.capacity_potential,
            capacity_cost=annuity(
                technology.at[(c.carrier, c.tech), "capex"],
                technology.at[
                    (c.carrier, c.tech), "lifetime"
                ],
                0.07,
            )
            * 1000,  # $/kW -> $/MW
            output_parameters={
                "emission_factor": (
                    carrier.at[c.carrier, "emission_factor"]
                    / technology.at[
                        (c.carrier, c.tech), "efficiency"
                    ]
                )
            },
        )
    )

Add Commodity¶

[12]:

for name, c in data["commodity"].iterrows():
    es.add(
        fc.Commodity(
            label=name,
            bus=buses[c.bus],
            carrier=c.carrier,
            tech=c.tech,
            amount=c.amount,
        )
    )

Add Link¶

[13]:

for name, c in data["link"].iterrows():
    es.add(
        fc.Link(
            label=name,
            from_bus=buses[c.from_bus],
            to_bus=buses[c.to_bus],
            capacity=c.capacity,
            expandable=True,
            capacity_cost=99999,
            loss=c.loss
        )
    )

Objective Function¶

The objective function is created from all instantiated objects. It will use all operating costs (i.e. marginal_cost argument) and if set all investment costs (i.e. capacity_cost argument)

\[\begin{split} \text{min:} \sum_g \sum_t \overbrace{c^{marginal\_cost}_g \cdot x^{flow}_{g}(t)}^{\text{operating cost}} \\ \sum_g \sum_t \overbrace{c^{capacity\_cost}_g \cdot x^{capacity}_{g}(t)}^{\text{investment cost}}\end{split}\]

Add Shortage/Excess Slack Components¶

[14]:

for name, e in data["excess"].iterrows():
    es.add(fc.Excess(label=name, bus=buses[e.bus]))

for name, s in data["shortage"].iterrows():
    es.add(
        fc.Shortage(
            label=name,
            carrier="electricity",
            tech="shortage",
            bus=buses[s.bus],
            marginal_cost=s.marginal_cost,
        )
    )

Creating the Mathematical Model¶

[15]:

# create model based on energy system and its components
m = Model(es)

# inspect objective function
# m.objective.pprint()

m.receive_duals()

WARNING: Reassigning the non-component attribute dual on block (model).Model
    with a new Component with type <class 'pyomo.core.base.suffix.Suffix'>.
    This is usually indicative of a modelling error. To avoid this warning,
    explicitly delete the attribute:
        del Model.dual
WARNING: Reassigning the non-component attribute rc on block (model).Model
    with a new Component with type <class 'pyomo.core.base.suffix.Suffix'>.
    This is usually indicative of a modelling error. To avoid this warning,
    explicitly delete the attribute:
        del Model.rc

Add CO2 Constraint¶

To add a CO2-constraint we will use the oemof.solph.constraints module which allows to add such a constraint in a easy way.

\[\sum_t \sum_f x^{flow}_f(t) \cdot c^{emission\_factor}_f \leq \overline{L_{CO_2}}\]

The constraint will sum all flows for the complete time horzion that have an attribute emission_factor and multiple the flow value with this factor.

[16]:

#m = constraints.emission_limit(m, limit=0)
m.write(io_options={"symbolic_solver_labels":True})

[16]:

('Model.lp', 139886760483792)

Solving the Model and Writing Results¶

[17]:

# check if cbc solver library is available
cbc = SolverFactory('cbc').available()

if cbc:
    #  solve the model using cbc solver
    m.solve("cbc")

    # write results back to the model object
    m.results = m.results()

    # writing results with the standard oemof-tabular output formatt
    pp.write_results(m, m.results, scenario_path)

    print(
        "Optimization done. Results are in {}.".format(
            results_path
        )
    )

    # write the lp-file
    # m.write(io_options={'symbolic_solver_labels': True})

WARNING: Could not locate the 'cbc' executable, which is required for solver
    cbc

Result Analysis¶

Plotting the Results¶

[18]:

if cbc:
    import os
    from plotly import offline
    from oemof.tabular.tools.plots import (
        hourly_plot,
        stacked_plot,
    )

    offline.init_notebook_mode(connected=True)


    offline.iplot(
        hourly_plot(
            scenario_name,
            "north-bus",
            os.path.join(
                os.path.expanduser("~"), "oemof-results"
            ),
            plot_filling_levels=False,
        ),
        filename=os.path.join(
            scenario_path, "hourly-plot.html"
        ),
    )