Home PV installation with battery storage

To get in touch with basic functionalities in oemof.solph, we will use this tutorial to model the electricity supply of a single family household that intents to install a PV plant in combination with a battery storage. Therefore several steps are considered:

  • Step 1: Definition of the scenario

  • Step 2: Adding fixed size PV

  • Step 3: PV investment optimisation

  • Step 4: PV investment optimisation with existing battery

  • Step 5: Full investment optimisation

  • Step 6: Autarky of the system

Step 1: Definition of the scenario

Imagine you want to set up a PV plant on top of a single family house. How would you find out which system fits best and why? Here are some possible points to think about:

  • Which area (e.g. on the roof) can be covered and how large is it?

  • How much energy will you consume or be able to feed into the grid?

  • How much will you have to pay for the system?

Especially for single-family homes without energy storage, the load can be very volatile: Switching on and off e.g. the oven can change the demand from several 10 Watts to several kW. This is particularly important when alignment of supply and demand. For this tutorial, we use hourly aggregated data from real-world measurements (pv_example_data.csv). This can be seen as a compromise: A finer resolution increases computatonal time but we want fast results in tutorial. The full dataset (which has 10-minute resolution) is described in the paper Dataset on electrical single-family house and heat pump load profiles in Germany. The PV time series has been created using PVGIS.

First, let us import all needed packages: os is handy to for file handling, matplotlib is for plotting, networkx as well as oemof.network.graph can be used to visualise the energy system graph, and pandas is for (input and output) data handling.

import os

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd

from oemof import solph

Now, let us load the input data and have an initial look at it. If the data provides time information, it is a good idea to directly consider it. This also allows to work with time stamps when handling the data. For illustation, we plot the first week of March 2020.


file_path = os.path.dirname(__file__)
filename = os.path.join(file_path, "pv_example_data.csv")
input_data = pd.read_csv(
    filename, index_col="timestep", parse_dates=["timestep"]
)

input_data.plot(drawstyle="steps")
plt.xlim(pd.Timestamp("2020-03-01 00:00"), pd.Timestamp("2020-03-07 00:00"))
plt.show()

The style has been chosen to reflect that the data is defined on the time intervals maked by the beginning time step. This means that 1.5 kW is probably not the peak load within the displayed week but only the peek hourly average. In contrast to the demand, the PV gain is given in normalised units and has to be multiplied by the installed PV capacity to yield the electricity production.

Input data
Input data

The model differentiates between time points and time steps. As solph supports mixing differnt length of time steps in one model, we need N+1 time points to define N time steps. If the index does not contain the last time point, the last interval can be infered based on index.freq.

Note

You need N+1 point in time to define N time steps.

The EnergySystem acts as a container. You will need to add everything that is part of the energy system.


energy_system = solph.EnergySystem(
    timeindex=input_data.index,
    infer_last_interval=True,
)

Now, to the actual energy system. It is modelled in the form of a mathematical graph. Energy flows along the edges (Flow) from node to node. First of all, where is an electricity Bus, which is used as a central point for the rest of the network to connect. Secondly, there is a Sink to model the enery demand. It takes an input from el_bus. Note that we set nominal_capacity=1 as the time series is already in kW. While this can be considered as a hack, it is common practice to do so: In the model, it is no issue to run in overload. As the last component, we add another Bus wich represents the electricity grid. Often, a Source would be used instead. However, to emphasise that the same grid is used for purchasing electricity and for feeding in, we use an nbalanced bus, also referred to as a “slack bus”.


el_bus = solph.Bus(label="electricity")

demand = solph.components.Sink(
    label="demand",
    inputs={
        el_bus: solph.Flow(
            nominal_capacity=1,
            fix=input_data["electricity demand (kW)"],
        )
    },
)

energy_system.add(el_bus, demand)

grid = solph.Bus(
    label="grid",
    outputs={el_bus: solph.Flow(variable_costs=0.3)},
    balanced=False,
)

energy_system.add(grid)

To obtain results, we have to create an optimisation model from the EnergySystem and solve it. To see diagnostic output, solve_kwargs={"tee": True} can be set, by default (False) it will be silent.

model = solph.Model(energy_system)

results = model.solve(solver="cbc", solve_kwargs={"tee": True})

If the model cannot find an optimal solution, model.solve() will fail. You can try this, e.g. by setting nominal_capacity=0.5 to the Flow between grid and el_bus (the one that has the variable costs defined). When the optimisation is finished, we can proceed with evaluating the results.


tce = results["objective"]
print(f"The total annual costs are {tce:.2f} €.")
el_costs = 0.3 * results["flow"][(grid, el_bus)].sum()
print(f"The annual costs for grid electricity are {el_costs:.2f} €.")

flows = results["flow"]

baseline = np.zeros(len(flows))

mode = "light"
# mode = "dark"
if mode == "dark":
    plt.style.use("dark_background")

plt.fill_between(
    flows.index,
    baseline,
    flows[("grid", "electricity")],
    step="pre",
    label="Grid supply",
)

plt.step(
    flows.index,
    flows[("electricity", "demand")],
    "-",
    color="darkgrey",
    label="Electricity demand",
)
plt.legend()
plt.ylabel("Power (kW)")
plt.xlim(pd.Timestamp("2020-01-01 00:00"), pd.Timestamp("2020-01-07 00:00"))
plt.gcf().autofmt_xdate()

plt.savefig(f"home_pv_result-1_{mode}.svg")

plt.show()

At first, we look at the total costs. As we used these for our objective value, it is directly accassible in the meta_results. Because we want to evaluate differnt contibutions later, we also manually calculate the costs manually. To no surprise the numbers are the same:

Table 4 Result overview

Quantity

Unit

Value

Annual costs for grid electricity

€

629.90

Total annual costs

€

629.90

To have a look at the flows from and to the “electricity” bus. You can use DataFrame.xs(...). For example flows_from = flows.xs("electricity", axis=1, level=0, drop_level=False), and with level=1 for flows entering the bus.

Result time-series
Result time-series

As grid supply is the only option, both lines of course overlap perfectly.

Learning

The model balances supply and demand along flows in a graph based model.

You can get the complete (uncommented) code for this step: home_pv_1.py

Click to display the code
# -*- coding: utf-8 -*-

"""
SPDX-FileCopyrightText: Patrik Schönfeldt
SPDX-FileCopyrightText: Daniel Niederhöfer
SPDX-FileCopyrightText: DLR e.V.

SPDX-License-Identifier: MIT
"""

# %%[imports]
import os

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd

from oemof import solph

# %%[input_data]

file_path = os.path.dirname(__file__)
filename = os.path.join(file_path, "pv_example_data.csv")
input_data = pd.read_csv(
    filename, index_col="timestep", parse_dates=["timestep"]
)

input_data.plot(drawstyle="steps")
plt.xlim(pd.Timestamp("2020-03-01 00:00"), pd.Timestamp("2020-03-07 00:00"))
plt.show()

# %%[energy_system]

energy_system = solph.EnergySystem(
    timeindex=input_data.index,
    infer_last_interval=True,
)

# %%[dispatch_model]

el_bus = solph.Bus(label="electricity")

demand = solph.components.Sink(
    label="demand",
    inputs={
        el_bus: solph.Flow(
            nominal_capacity=1,
            fix=input_data["electricity demand (kW)"],
        )
    },
)

energy_system.add(el_bus, demand)

grid = solph.Bus(
    label="grid",
    outputs={el_bus: solph.Flow(variable_costs=0.3)},
    balanced=False,
)

energy_system.add(grid)
# %%[graph_plotting]
plt.figure()
graph = energy_system.to_networkx()
nx.draw(graph, with_labels=True, font_size=8)
plt.show()
# %%[model_optimisation]
model = solph.Model(energy_system)

results = model.solve(solver="cbc", solve_kwargs={"tee": True})


# %%[results]

tce = results["objective"]
print(f"The total annual costs are {tce:.2f} €.")
el_costs = 0.3 * results["flow"][(grid, el_bus)].sum()
print(f"The annual costs for grid electricity are {el_costs:.2f} €.")

flows = results["flow"]

baseline = np.zeros(len(flows))

mode = "light"
# mode = "dark"
if mode == "dark":
    plt.style.use("dark_background")

plt.fill_between(
    flows.index,
    baseline,
    flows[("grid", "electricity")],
    step="pre",
    label="Grid supply",
)

plt.step(
    flows.index,
    flows[("electricity", "demand")],
    "-",
    color="darkgrey",
    label="Electricity demand",
)
plt.legend()
plt.ylabel("Power (kW)")
plt.xlim(pd.Timestamp("2020-01-01 00:00"), pd.Timestamp("2020-01-07 00:00"))
plt.gcf().autofmt_xdate()

plt.savefig(f"home_pv_result-1_{mode}.svg")

plt.show()

Step 2: Adding fixed size PV

So far so good. But what would happen if we had an existing PV plant? Let us simply add one to the system by creating a suitable source. We might have five kW installed and therefore enter a nominal capacity of five.


pv_size = 5  # kW
pv_specific_costs = 1500  # €/kW
pv_lifetime = 20  # years

pv_system = solph.components.Source(
    label="PV",
    outputs={
        el_bus: solph.Flow(
            nominal_capacity=pv_size,
            maximum=input_data["pv yield (kW/kW)"],
        )
    },
)

energy_system.add(pv_system)

Note that we set the PV yield time series using the max parameter: It is always possible to not take the energy. In our case, however, we want to allow feeding into the grid. Thus, the corresponding Bus needs to accept an incoming Flow. It is possible, to add Flow s to existing nodes:


grid = solph.Bus(
    label="grid",
    outputs={el_bus: solph.Flow(variable_costs=0.3)},
    balanced=False,
)

grid.inputs[el_bus] = solph.Flow(variable_costs=-0.06)

However, the Flow is typically directly added when defining the grid. So, alternatively, you can just change the previos definition to look like:


grid = solph.Bus(
    label="grid",
    inputs={el_bus: solph.Flow(variable_costs=-0.06)},
    outputs={el_bus: solph.Flow(variable_costs=0.3)},
    balanced=False,
)

Note that feeding into the grid will make sense if there is compensation (or if the PV output is set to have a fix production instead of a maximum one). For our example, we take 6 ct/kWh as compensation for feed-in.

As building the PV system is free for the model (it is already present), we manually add an annuity to have the PV system included in the total costs.


pv_annuity = pv_size * pv_specific_costs / pv_lifetime
annual_grid_supply = results["flow"][(grid, el_bus)].sum()
el_costs = 0.3 * annual_grid_supply

el_revenue = 0.1 * results["flow"][(el_bus, grid)].sum()

tce = pv_annuity + results["objective"]
Result time-series
Result time-series

In this scenario, the objective value is negative, so a revenue is achieved. This is costly because of feeding in, demand and supply do hardly match in this scenario. However, correcting this by adding the annuity of the PV system leaves us with (positive) costs.

Table 5 Result overview

Quantity

Unit

Value (no PV)

Value (5 kW PV)

Annual costs for grid electricity

€

629.90

365.32

Annual revenue from feed-in

€

451.61

Annuity for the PV system

€

375.00

Total annual costs

€

629.90

479.03

Autarky

%

0.00

42.00

Learning

You can give negative costs to model a revenue.

You can get the complete (uncommented) code for this step: home_pv_2.py

Click to display the code
# -*- coding: utf-8 -*-

"""
SPDX-FileCopyrightText: Patrik Schönfeldt
SPDX-FileCopyrightText: Daniel Niederhöfer
SPDX-FileCopyrightText: DLR e.V.

SPDX-License-Identifier: MIT
"""

# %%[imports]
import os

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd

from oemof import solph

# %%[input_data]

file_path = os.path.dirname(__file__)
filename = os.path.join(file_path, "pv_example_data.csv")
input_data = pd.read_csv(
    filename, index_col="timestep", parse_dates=["timestep"]
)

# %%[energy_system]

energy_system = solph.EnergySystem(
    timeindex=input_data.index,
    infer_last_interval=True,
)

# %%[dispatch_model]

el_bus = solph.Bus(label="electricity")

demand = solph.components.Sink(
    label="demand",
    inputs={
        el_bus: solph.Flow(
            nominal_capacity=1,
            fix=input_data["electricity demand (kW)"],
        )
    },
)

energy_system.add(el_bus, demand)

# %%[grid_feedin]

grid = solph.Bus(
    label="grid",
    outputs={el_bus: solph.Flow(variable_costs=0.3)},
    balanced=False,
)

grid.inputs[el_bus] = solph.Flow(variable_costs=-0.06)

# %%[add_grid]

energy_system.add(grid)

# %%[pv_system]

pv_size = 5  # kW
pv_specific_costs = 1500  # €/kW
pv_lifetime = 20  # years

pv_system = solph.components.Source(
    label="PV",
    outputs={
        el_bus: solph.Flow(
            nominal_capacity=pv_size,
            maximum=input_data["pv yield (kW/kW)"],
        )
    },
)

energy_system.add(pv_system)
# %%[graph_plotting]
plt.figure()
graph = energy_system.to_networkx()
nx.draw(graph, with_labels=True, font_size=8)
# %%[model_optimisation]
model = solph.Model(energy_system)

results = model.solve(solver="cbc", solve_kwargs={"tee": True})

# %%[results]

pv_annuity = pv_size * pv_specific_costs / pv_lifetime
annual_grid_supply = results["flow"][(grid, el_bus)].sum()
el_costs = 0.3 * annual_grid_supply

el_revenue = 0.1 * results["flow"][(el_bus, grid)].sum()

tce = pv_annuity + results["objective"]
# %%[result_plotting]

print(f"The annual costs for grid electricity are {el_costs:.2f} €.")
print(f"The annual revenue from feed-in is {el_revenue:.2f} €.")
print(f"The annuity for the PV system is {pv_annuity:.2f} €.")
print(f"The total annual costs are {tce:.2f} €.")

annual_demand = input_data["electricity demand (kW)"].sum()

print(
    f"Autarky is 1 - {annual_grid_supply:.2f} kWh / {annual_demand:.2f} kWh"
    + f" = {100 - 100 * annual_grid_supply / annual_demand:.2f} %."
)

flows = results["flow"]
baseline = np.zeros(len(flows))

plt.figure()

mode = "light"
# mode = "dark"
if mode == "dark":
    plt.style.use("dark_background")

plt.fill_between(
    flows.index,
    baseline,
    baseline + flows[("grid", "electricity")],
    step="pre",
    label="Grid supply",
)

baseline += flows[("grid", "electricity")]

plt.fill_between(
    flows.index,
    baseline,
    baseline + flows[("PV", "electricity")],
    step="pre",
    label="PV supply",
)

plt.step(
    flows.index,
    flows[("electricity", "demand")],
    "-",
    color="darkgrey",
    label="Electricity demand",
)

plt.step(
    flows.index,
    flows[("electricity", "demand")] + flows[("electricity", "grid")],
    ":",
    color="darkgrey",
    label="Feed-In",
)

plt.legend()
plt.ylabel("Power (kW)")
plt.xlim(pd.Timestamp("2020-02-21 00:00"), pd.Timestamp("2020-02-28 00:00"))
plt.gcf().autofmt_xdate()

plt.savefig(f"home_pv_result-2_{mode}.svg")

plt.show()

Step 3: PV investment optimisation

Now, compared to dispatch with a fixed PV plant, almost everything else stays the same, except for one thing: the nominal capacity of the PV plant. That’s because we want to find out which peak power the system should have to minimise our costs. Therefore an investment object with a periodical cost of 75 Euros per kW is assigned. These costs represent an assumed investment cost of 1500 Euros per kW divided by an estimated life time of 20 years as a straight-line deprecation. To make sure the model converges, we set a maximum capacity: If building the PV system was profitable, the optimiser would try to build an infinite size PV system, which would effectively prevent the model from converging.


pv_specific_costs = 1500  # €/kW
pv_lifetime = 20  # years
pv_epc = pv_specific_costs / pv_lifetime

pv_system = solph.components.Source(
    label="PV",
    outputs={
        el_bus: solph.Flow(
            nominal_capacity=solph.Investment(ep_costs=pv_epc, maximum=10),
            maximum=input_data["pv yield (kW/kW)"],
        )
    },
)

energy_system.add(pv_system)

However, in this case we did not need to set the limit. The assumed annuity is so high that the costs of the produced electricity are above 6 ct/kWh. The resulting optimal PV size


# Potentially, there are more points in time to invest.
# For upfront invest, we need to select the initial.
pv_size = results["invest"][(pv_system, el_bus)][0]

is 4.13 kW, which is even smaller than what we tried before.

Table 6 Result overview

Quantity

Unit

Value

Optimal PV size

kW

4.13

Annual costs for grid electricity

€

376.59

Annual revenue from feed-in

€

347.64

Annuity for the PV system

€

309.40

Total annual costs

€

477.41

Autarky

%

40.21

You can get the complete (uncommented) code for this step: home_pv_3.py

Click to display the code
# -*- coding: utf-8 -*-

"""
SPDX-FileCopyrightText: Patrik Schönfeldt
SPDX-FileCopyrightText: Daniel Niederhöfer
SPDX-FileCopyrightText: DLR e.V.

SPDX-License-Identifier: MIT
"""

# %%[imports]
import os

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd

from oemof import solph

# %%[input_data]

file_path = os.path.dirname(__file__)
filename = os.path.join(file_path, "pv_example_data.csv")
input_data = pd.read_csv(
    filename, index_col="timestep", parse_dates=["timestep"]
)

# %%[energy_system]

energy_system = solph.EnergySystem(
    timeindex=input_data.index,
    infer_last_interval=True,
)

# %%[dispatch_model]

el_bus = solph.Bus(label="electricity")

demand = solph.components.Sink(
    label="demand",
    inputs={
        el_bus: solph.Flow(
            nominal_capacity=1,
            fix=input_data["electricity demand (kW)"],
        )
    },
)

energy_system.add(el_bus, demand)


# %%[grid_conection]

grid = solph.Bus(
    label="grid",
    inputs={el_bus: solph.Flow(variable_costs=-0.06)},
    outputs={el_bus: solph.Flow(variable_costs=0.3)},
    balanced=False,
)

# %%[add_grid]

energy_system.add(grid)

# %%[pv_system]

pv_specific_costs = 1500  # €/kW
pv_lifetime = 20  # years
pv_epc = pv_specific_costs / pv_lifetime

pv_system = solph.components.Source(
    label="PV",
    outputs={
        el_bus: solph.Flow(
            nominal_capacity=solph.Investment(ep_costs=pv_epc, maximum=10),
            maximum=input_data["pv yield (kW/kW)"],
        )
    },
)

energy_system.add(pv_system)

# %%[graph_plotting]
plt.figure()
graph = energy_system.to_networkx()
nx.draw(graph, with_labels=True, font_size=8)
# %%[model_optimisation]
model = solph.Model(energy_system)

results = model.solve(solver="cbc", solve_kwargs={"tee": True})

# %%[result_pv]

# Potentially, there are more points in time to invest.
# For upfront invest, we need to select the initial.
pv_size = results["invest"][(pv_system, el_bus)][0]

# %%[results]

pv_annuity = pv_epc * pv_size
annual_grid_supply = results["flow"][(grid, el_bus)].sum()
el_costs = 0.3 * annual_grid_supply
el_revenue = 0.1 * results["flow"][(el_bus, grid)].sum()

tce = results["objective"]

# %%[result_plotting]

print(f"The optimal PV size is {pv_size:.2f} kW.")

print(f"The annual costs for grid electricity are {el_costs:.2f} €.")
print(f"The annual revenue from feed-in is {el_revenue:.2f} €.")
print(f"The annuity for the PV system is {pv_annuity:.2f} €.")
print(f"The total annual costs are {tce:.2f} €.")

annual_demand = input_data["electricity demand (kW)"].sum()

print(
    f"Autarky is 1 - {annual_grid_supply:.2f} kWh / {annual_demand:.2f} kWh"
    + f" = {100 - 100 * annual_grid_supply / annual_demand:.2f} %."
)

flows = results["flow"]

baseline = np.zeros(len(flows))

plt.figure()

mode = "light"
# mode = "dark"
if mode == "dark":
    plt.style.use("dark_background")

plt.fill_between(
    flows.index,
    baseline,
    baseline + flows[("grid", "electricity")],
    step="pre",
    label="Grid supply",
)

baseline += flows[("grid", "electricity")]

plt.fill_between(
    flows.index,
    baseline,
    baseline + flows[("PV", "electricity")],
    step="pre",
    label="PV supply",
)

plt.step(
    flows.index,
    flows[("electricity", "demand")],
    "-",
    color="darkgrey",
    label="Electricity demand",
)

plt.step(
    flows.index,
    flows[("electricity", "demand")] + flows[("electricity", "grid")],
    ":",
    color="darkgrey",
    label="Feed-In",
)

plt.legend()
plt.ylabel("Power (kW)")
plt.xlim(pd.Timestamp("2020-02-21 00:00"), pd.Timestamp("2020-02-28 00:00"))
plt.gcf().autofmt_xdate()

plt.savefig(f"home_pv_result-3_{mode}.svg")

plt.show()

Learning

You can use an Investment object for the nominal_capacity to make the capacity an optimisation variable.

Step 4: PV investment optimisation with existing battery

In the previous step we learned that in our assumed scenario, compensation for feeding into the grid will not pay the PV system. So, what if we had a big battery? Let us just add one:


battery_specific_costs = 1000  # €/kW
battery_lifetime = 10  # years
battery_epc = battery_specific_costs / battery_lifetime
battery_size = 10  # kWh

battery = solph.components.GenericStorage(
    label="Battery",
    nominal_capacity=battery_size,
    inputs={el_bus: solph.Flow()},
    outputs={el_bus: solph.Flow()},
    inflow_conversion_factor=0.9,
    loss_rate=0.01,
)

energy_system.add(battery)

As before with the PV system, we manually add the deprecation of the battery to the objective value to get the total costs of the system.

Result time-series
Result time-series

As can be seen in the time series, the Battery allows to use much of the produced PV electricity. Thus, total PV size is not increased too much. With a battery that large, the it now dominates the total costs.

Table 7 Result overview

Quantity

Unit

Value

Optimal PV size

kW

5.67

Annual costs for grid electricity

€

72.56

Annual revenue from feed-in

€

383.39

Annuity for the PV system

€

425.12

Annuity for the battery

€

1000.00

Total annual costs

€

1267.65

Autarky

%

88.48

Learning

Energy storage is modelled using the GenericStorage class.

You can get the complete (uncommented) code for this step: home_pv_4.py

Click to display the code
# -*- coding: utf-8 -*-

"""
SPDX-FileCopyrightText: Patrik Schönfeldt
SPDX-FileCopyrightText: Daniel Niederhöfer
SPDX-FileCopyrightText: DLR e.V.

SPDX-License-Identifier: MIT
"""

# %%[imports]
import os

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd

from oemof import solph

# %%[input_data]

file_path = os.path.dirname(__file__)
filename = os.path.join(file_path, "pv_example_data.csv")
input_data = pd.read_csv(
    filename, index_col="timestep", parse_dates=["timestep"]
)

# %%[energy_system]

energy_system = solph.EnergySystem(
    timeindex=input_data.index,
    infer_last_interval=True,
)

# %%[dispatch_model]

el_bus = solph.Bus(label="electricity")

demand = solph.components.Sink(
    label="demand",
    inputs={
        el_bus: solph.Flow(
            nominal_capacity=1,
            fix=input_data["electricity demand (kW)"],
        )
    },
)

energy_system.add(el_bus, demand)

grid = solph.Bus(
    label="grid",
    inputs={el_bus: solph.Flow(variable_costs=-0.06)},
    outputs={el_bus: solph.Flow(variable_costs=0.3)},
    balanced=False,
)

energy_system.add(grid)

pv_specific_costs = 1500  # €/kW
pv_lifetime = 20  # years
pv_epc = pv_specific_costs / pv_lifetime

pv_system = solph.components.Source(
    label="PV",
    outputs={
        el_bus: solph.Flow(
            nominal_capacity=solph.Investment(ep_costs=pv_epc, maximum=10),
            maximum=input_data["pv yield (kW/kW)"],
        )
    },
)

energy_system.add(pv_system)

# %%[battery]

battery_specific_costs = 1000  # €/kW
battery_lifetime = 10  # years
battery_epc = battery_specific_costs / battery_lifetime
battery_size = 10  # kWh

battery = solph.components.GenericStorage(
    label="Battery",
    nominal_capacity=battery_size,
    inputs={el_bus: solph.Flow()},
    outputs={el_bus: solph.Flow()},
    inflow_conversion_factor=0.9,
    loss_rate=0.01,
)

energy_system.add(battery)

# %%[graph_plotting]
plt.figure()
graph = energy_system.to_networkx()
nx.draw(graph, with_labels=True, font_size=8)
# %%[model_optimisation]
model = solph.Model(energy_system)

results = model.solve(solver="cbc", solve_kwargs={"tee": True})

# %%[results]

# Potentially, there are more points in time to invest.
# For upfront invest, we need to select the initial.
pv_size = results["invest"][(pv_system, el_bus)][0]

battery_annuity = battery_epc * battery_size
pv_annuity = pv_epc * results["invest"][(pv_system, el_bus)][0]
annual_grid_supply = results["flow"][(grid, el_bus)].sum()
el_costs = 0.3 * annual_grid_supply
el_revenue = 0.1 * results["flow"][(el_bus, grid)].sum()

tce = results["objective"] + battery_annuity

print(f"The optimal PV size is {pv_size:.2f} kW.")

print(f"The annual costs for grid electricity are {el_costs:.2f} €.")
print(f"The annual revenue from feed-in is {el_revenue:.2f} €.")
print(f"The annuity for the PV system is {pv_annuity:.2f} €.")
print(f"The annuity for the battery is {battery_annuity:.2f} €.")
print(f"The total annual costs are {tce:.2f} €.")

annual_demand = input_data["electricity demand (kW)"].sum()

print(
    f"Autarky is 1 - {annual_grid_supply:.2f} kWh / {annual_demand:.2f} kWh"
    + f" = {100 - 100 * annual_grid_supply / annual_demand:.2f} %."
)

flows = results["flow"]

baseline = np.zeros(len(flows))

plt.figure()

mode = "light"
# mode = "dark"
if mode == "dark":
    plt.style.use("dark_background")

plt.fill_between(
    flows.index,
    baseline,
    baseline + flows[("grid", "electricity")],
    step="pre",
    label="Grid supply",
)

baseline += flows[("grid", "electricity")]

plt.fill_between(
    flows.index,
    baseline,
    baseline + flows[("PV", "electricity")],
    step="pre",
    label="PV supply",
)

baseline += flows[("PV", "electricity")]

plt.fill_between(
    flows.index,
    baseline,
    baseline + flows[("Battery", "electricity")],
    step="pre",
    label="Battery supply",
)

plt.step(
    flows.index,
    flows[("electricity", "demand")],
    "-",
    color="darkgrey",
    label="Electricity demand",
)

plt.step(
    flows.index,
    flows[("electricity", "demand")] + flows[("electricity", "Battery")],
    "--",
    color="darkgrey",
    label="Battery charging",
)

plt.step(
    flows.index,
    flows[("electricity", "demand")]
    + flows[("electricity", "Battery")]
    + flows[("electricity", "grid")],
    ":",
    color="darkgrey",
    label="Feed-In",
)

plt.legend()
plt.ylabel("Power (kW)")
plt.xlim(pd.Timestamp("2020-02-21 00:00"), pd.Timestamp("2020-02-28 00:00"))
plt.gcf().autofmt_xdate()

plt.savefig(f"home_pv_result-4_{mode}.svg")

plt.show()

Step 5: Full investment optimisation

For a full investment optimisation, of course, we include optimising the battery size. First, however, we will separate the PV system into PV panels and inverter. For clearity, we rename the former el_bus to ac_bus and create a separate dc_bus. As the numbers given in the input data already consider inverter losses, we compensate for that by allowing higer DC gains from the panels, the inverter then has a conversion_factor so that the Flow to the AC bus is reduced by 5 %.

dc_bus = solph.Bus(label="DC")

pv_specific_costs = 1200  # €/kW
pv_lifetime = 20  # years
pv_epc = pv_specific_costs / pv_lifetime

pv_panels = solph.components.Source(
    label="PV",
    outputs={
        dc_bus: solph.Flow(
            nominal_capacity=solph.Investment(ep_costs=pv_epc, maximum=10),
            maximum=input_data["pv yield (kW/kW)"] / 0.95,
        )
    },
)

inverter_specific_costs = 300  # €/kW
inverter_lifetime = 20  # years
inverter_epc = inverter_specific_costs / inverter_lifetime

inverter = solph.components.Converter(
    label="inverter",
    inputs={
        dc_bus: solph.Flow(
            nominal_capacity=solph.Investment(ep_costs=inverter_epc)
        )
    },
    outputs={ac_bus: solph.Flow()},
    conversion_factors={ac_bus: 0.95},
)

energy_system.add(dc_bus, pv_panels, inverter)

Now, the optimiser can choose on the inverter and number of PV panels separately. Note that the costs for inverter and pv panels add up to the value we used in the previous step. For the sake of simplicity in this tutorial, the battery will remain coupled to the AC system. We only let the optimisation decide on the size.

battery_specific_costs = 750  # €/kW
battery_lifetime = 20  # years
battery_epc = battery_specific_costs / battery_lifetime
battery_size = solph.Investment(ep_costs=battery_epc)

battery = solph.components.GenericStorage(
    label="Battery",
    nominal_capacity=battery_size,
    inputs={ac_bus: solph.Flow()},
    outputs={ac_bus: solph.Flow()},
    inflow_conversion_factor=0.9,
    loss_rate=0.01,
)

energy_system.add(battery)

We find that it is advisible to have an inverter that can only output two thirds of the nominal power of the PV panels. While the inverter is even smaller than in the previous steps, the total power of the panels is increased.

Table 8 Result overview

Quantity

Unit

Value

Optimal PV size

kW

6.04

Optimal inverter size

kW

4.10

Optimal battery size

kWh

2.39

Annual costs for grid electricity

€

156.35

Annual revenue from feed-in

€

451.61

Annuity for the PV system

€

362.27

Annuity for the battery

€

89.53

Total annual costs

€

398.66

Autarky

%

75.18

Learning

Converter s are used to model (possibly lossy) conversion.

You can get the complete (uncommented) code for this step: home_pv_5.py

Click to display the code
# -*- coding: utf-8 -*-

"""
SPDX-FileCopyrightText: Patrik Schönfeldt
SPDX-FileCopyrightText: Daniel Niederhöfer
SPDX-FileCopyrightText: DLR e.V.

SPDX-License-Identifier: MIT
"""

# %%[imports]
import os

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd

from oemof import solph

# %%[input_data]

file_path = os.path.dirname(__file__)
filename = os.path.join(file_path, "pv_example_data.csv")
input_data = pd.read_csv(
    filename, index_col="timestep", parse_dates=["timestep"]
)

# %%[energy_system]

energy_system = solph.EnergySystem(
    timeindex=input_data.index,
    infer_last_interval=True,
)

# %%[dispatch_model]

ac_bus = solph.Bus(label="electricity")

demand = solph.components.Sink(
    label="demand",
    inputs={
        ac_bus: solph.Flow(
            nominal_capacity=1,
            fix=input_data["electricity demand (kW)"],
        )
    },
)

energy_system.add(ac_bus, demand)

grid = solph.Bus(
    label="grid",
    inputs={ac_bus: solph.Flow(variable_costs=-0.06)},
    outputs={ac_bus: solph.Flow(variable_costs=0.3)},
    balanced=False,
)

energy_system.add(grid)

# %%[pv_system]
dc_bus = solph.Bus(label="DC")

pv_specific_costs = 1200  # €/kW
pv_lifetime = 20  # years
pv_epc = pv_specific_costs / pv_lifetime

pv_panels = solph.components.Source(
    label="PV",
    outputs={
        dc_bus: solph.Flow(
            nominal_capacity=solph.Investment(ep_costs=pv_epc, maximum=10),
            maximum=input_data["pv yield (kW/kW)"] / 0.95,
        )
    },
)

inverter_specific_costs = 300  # €/kW
inverter_lifetime = 20  # years
inverter_epc = inverter_specific_costs / inverter_lifetime

inverter = solph.components.Converter(
    label="inverter",
    inputs={
        dc_bus: solph.Flow(
            nominal_capacity=solph.Investment(ep_costs=inverter_epc)
        )
    },
    outputs={ac_bus: solph.Flow()},
    conversion_factors={ac_bus: 0.95},
)

energy_system.add(dc_bus, pv_panels, inverter)

# %%[battery]
battery_specific_costs = 750  # €/kW
battery_lifetime = 20  # years
battery_epc = battery_specific_costs / battery_lifetime
battery_size = solph.Investment(ep_costs=battery_epc)

battery = solph.components.GenericStorage(
    label="Battery",
    nominal_capacity=battery_size,
    inputs={ac_bus: solph.Flow()},
    outputs={ac_bus: solph.Flow()},
    inflow_conversion_factor=0.9,
    loss_rate=0.01,
)

energy_system.add(battery)

# %%[graph_plotting]
plt.figure()
graph = energy_system.to_networkx()
nx.draw(graph, with_labels=True, font_size=8)
# %%[model_optimisation]
model = solph.Model(energy_system)

results = model.solve(solver="cbc", solve_kwargs={"tee": True})

# %%[results]

# Potentially, there are more points in time to invest.
# For upfront invest, we need to select the initial.
pv_size = results["invest"][(pv_panels, dc_bus)][0]
pv_annuity = pv_epc * pv_size

inverter_size = results["invest"][(dc_bus, inverter)][0]
inverter_annuity = inverter_epc * inverter_size

battery_size = results["invest"][battery][0]
battery_annuity = battery_epc * battery_size

annual_grid_supply = results["flow"][(grid, ac_bus)].sum()
el_costs = 0.3 * annual_grid_supply
el_revenue = 0.1 * results["flow"][(ac_bus, grid)].sum()

tce = results["objective"]

print(f"The optimal PV size is {pv_size:.2f} kW.")
print(f"The optimal inverter size is {inverter_size:.2f} kW.")
print(f"The optimal battery size is {battery_size:.2f} kWh.")

print(f"The annual costs for grid electricity are {el_costs:.2f} €.")
print(f"The annual revenue from feed-in is {el_revenue:.2f} €.")
print(f"The annuity for the PV system is {pv_annuity:.2f} €.")
print(f"The annuity for the battery is {battery_annuity:.2f} €.")
print(f"The total annual costs are {tce:.2f} €.")

annual_demand = input_data["electricity demand (kW)"].sum()

print(
    f"Autarky is 1 - {annual_grid_supply:.2f} kWh / {annual_demand:.2f} kWh"
    + f" = {100 - 100 * annual_grid_supply / annual_demand:.2f} %."
)

flows = results["flow"]

baseline = np.zeros(len(flows))

plt.figure()

mode = "light"
# mode = "dark"
if mode == "dark":
    plt.style.use("dark_background")

plt.fill_between(
    flows.index,
    baseline,
    baseline + flows[("grid", "electricity")],
    step="pre",
    label="Grid supply",
)

baseline += flows[("grid", "electricity")]

plt.fill_between(
    flows.index,
    baseline,
    baseline + flows[("inverter", "electricity")],
    step="pre",
    label="PV supply",
)

baseline += flows[("inverter", "electricity")]

plt.fill_between(
    flows.index,
    baseline,
    baseline + flows[("Battery", "electricity")],
    step="pre",
    label="Battery supply",
)

plt.step(
    flows.index,
    flows[("electricity", "demand")],
    "-",
    color="darkgrey",
    label="Electricity demand",
)

plt.step(
    flows.index,
    flows[("electricity", "demand")] + flows[("electricity", "Battery")],
    "--",
    color="darkgrey",
    label="Battery charging",
)

plt.step(
    flows.index,
    flows[("electricity", "demand")]
    + flows[("electricity", "Battery")]
    + flows[("electricity", "grid")],
    ":",
    color="darkgrey",
    label="Feed-In",
)

plt.legend()
plt.ylabel("Power (kW)")
plt.xlim(pd.Timestamp("2020-02-21 00:00"), pd.Timestamp("2020-02-28 00:00"))
plt.gcf().autofmt_xdate()

plt.savefig(f"home_pv_result-5_{mode}.svg")

plt.show()

Step 6: Autarky of the system

The final step of the tutorial shows a way how to limit the energy supply from the grid in order to guarantee a minimum autarky grade of our system consisting of a combined PV and battery storage investment.

Since our total energy demand is almost 2100 kWh, we could accomplish an autarky grade of 90 percent if we limited the sum of thegrid energy supply to 210 kWh. Choosing 42 kW as a reasonable dimensioned power, multiplying it by five full load hours does the trick.


grid = solph.Bus(
    label="grid",
    inputs={ac_bus: solph.Flow(variable_costs=-0.06)},
    outputs={
        ac_bus: solph.Flow(
            nominal_capacity=42,
            full_load_time_max=5,
            variable_costs=0.3,
        )
    },
    balanced=False,
)

energy_system.add(grid)

So, finally let us compare the results from all different steps.

Table 9 Result overview

Quantity

Unit

Value (no PV)

Value (5 kW PV)

Value (opt. PV)

Value (big bat.)

Value (all opt.)

Value (Autarky)

Optimal PV size

kW

4.13

5.67

6.04

10

Optimal inverter size

kW

4.1

6.79

Optimal battery size

kWh

2.39

5.01

Annual costs for grid electricity

€

629.9

365.32

376.59

72.56

156.35

63

Annual revenue from feed-in

€

451.61

347.64

383.39

451.61

821.78

Annuity for the PV system

€

375

309.4

425.12

362.27

600

Annuity for the battery

€

1000.00

89.53

188.02

Total annual costs

€

629.9

479.03

477.41

1267.65

398.66

459.76

Autarky

%

0

42

40.21

88.48

75.18

90

Learning

You can define the full load time of a Flow to limit its energy transfer. Here, this has been used to accieve a minimum autarky.

You can get the complete (uncommented) code for this step: home_pv_6.py

Click to display the code
# -*- coding: utf-8 -*-

"""
SPDX-FileCopyrightText: Patrik Schönfeldt
SPDX-FileCopyrightText: Daniel Niederhöfer
SPDX-FileCopyrightText: DLR e.V.

SPDX-License-Identifier: MIT
"""

# %%[imports]
import os

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd

from oemof import solph

# %%[input_data]

file_path = os.path.dirname(__file__)
filename = os.path.join(file_path, "pv_example_data.csv")
input_data = pd.read_csv(
    filename, index_col="timestep", parse_dates=["timestep"]
)

# %%[energy_system]

energy_system = solph.EnergySystem(
    timeindex=input_data.index,
    infer_last_interval=True,
)

# %%[dispatch_model]

ac_bus = solph.Bus(label="electricity")

demand = solph.components.Sink(
    label="demand",
    inputs={
        ac_bus: solph.Flow(
            nominal_capacity=1,
            fix=input_data["electricity demand (kW)"],
        )
    },
)

energy_system.add(ac_bus, demand)

# %%[grid]

grid = solph.Bus(
    label="grid",
    inputs={ac_bus: solph.Flow(variable_costs=-0.06)},
    outputs={
        ac_bus: solph.Flow(
            nominal_capacity=42,
            full_load_time_max=5,
            variable_costs=0.3,
        )
    },
    balanced=False,
)

energy_system.add(grid)

# %%[pv_system]
dc_bus = solph.Bus(label="DC")

pv_specific_costs = 1200  # €/kW
pv_lifetime = 20  # years
pv_epc = pv_specific_costs / pv_lifetime

pv_panels = solph.components.Source(
    label="PV",
    outputs={
        dc_bus: solph.Flow(
            nominal_capacity=solph.Investment(ep_costs=pv_epc, maximum=10),
            maximum=input_data["pv yield (kW/kW)"] / 0.95,
        )
    },
)

inverter_specific_costs = 300  # €/kW
inverter_lifetime = 20  # years
inverter_epc = inverter_specific_costs / inverter_lifetime

inverter = solph.components.Converter(
    label="inverter",
    inputs={
        dc_bus: solph.Flow(
            nominal_capacity=solph.Investment(ep_costs=inverter_epc)
        )
    },
    outputs={ac_bus: solph.Flow()},
    conversion_factors={ac_bus: 0.95},
)

energy_system.add(dc_bus, pv_panels, inverter)

# %%[battery]
battery_specific_costs = 750  # €/kW
battery_lifetime = 20  # years
battery_epc = battery_specific_costs / battery_lifetime
battery_size = solph.Investment(ep_costs=battery_epc)

battery = solph.components.GenericStorage(
    label="Battery",
    nominal_capacity=battery_size,
    inputs={ac_bus: solph.Flow()},
    outputs={ac_bus: solph.Flow()},
    inflow_conversion_factor=0.9,
    loss_rate=0.01,
)

energy_system.add(battery)

# %%[graph_plotting]
plt.figure()
graph = energy_system.to_networkx()
nx.draw(graph, with_labels=True, font_size=8)
# %%[model_optimisation]
model = solph.Model(energy_system)

results = model.solve(solver="cbc", solve_kwargs={"tee": True})

# %%[results]

pv_size = results["invest"][(pv_panels, dc_bus)][0]
pv_annuity = pv_epc * pv_size

inverter_size = results["invest"][(dc_bus, inverter)][0]
inverter_annuity = inverter_epc * inverter_size

battery_size = results["invest"][battery][0]
battery_annuity = battery_epc * battery_size

annual_grid_supply = results["flow"][(grid, ac_bus)].sum()
el_costs = 0.3 * annual_grid_supply
el_revenue = 0.1 * results["flow"][(ac_bus, grid)].sum()

tce = results["objective"]

print(f"The optimal PV size is {pv_size:.2f} kW.")
print(f"The optimal inverter size is {inverter_size:.2f} kW.")
print(f"The optimal battery size is {battery_size:.2f} kWh.")

print(f"The annual costs for grid electricity are {el_costs:.2f} €.")
print(f"The annual revenue from feed-in is {el_revenue:.2f} €.")
print(f"The annuity for the PV system is {pv_annuity:.2f} €.")
print(f"The annuity for the battery is {battery_annuity:.2f} €.")
print(f"The total annual costs are {tce:.2f} €.")

annual_demand = input_data["electricity demand (kW)"].sum()

print(
    f"Autarky is 1 - {annual_grid_supply:.2f} kWh / {annual_demand:.2f} kWh"
    + f" = {100 - 100 * annual_grid_supply / annual_demand:.2f} %."
)


flows = results["flow"]

baseline = np.zeros(len(flows))

plt.figure()

mode = "light"
# mode = "dark"
if mode == "dark":
    plt.style.use("dark_background")

plt.fill_between(
    flows.index,
    baseline,
    baseline + flows[("grid", "electricity")],
    step="pre",
    label="Grid supply",
)

baseline += flows[("grid", "electricity")]

plt.fill_between(
    flows.index,
    baseline,
    baseline + flows[("inverter", "electricity")],
    step="pre",
    label="PV supply",
)

baseline += flows[("inverter", "electricity")]

plt.fill_between(
    flows.index,
    baseline,
    baseline + flows[("Battery", "electricity")],
    step="pre",
    label="Battery supply",
)

plt.step(
    flows.index,
    flows[("electricity", "demand")],
    "-",
    color="darkgrey",
    label="Electricity demand",
)

plt.step(
    flows.index,
    flows[("electricity", "demand")] + flows[("electricity", "Battery")],
    "--",
    color="darkgrey",
    label="Battery charging",
)

plt.step(
    flows.index,
    flows[("electricity", "demand")]
    + flows[("electricity", "Battery")]
    + flows[("electricity", "grid")],
    ":",
    color="darkgrey",
    label="Feed-In",
)

plt.legend()
plt.ylabel("Power (kW)")
plt.xlim(pd.Timestamp("2020-02-21 00:00"), pd.Timestamp("2020-02-28 00:00"))
plt.gcf().autofmt_xdate()

plt.savefig(f"home_pv_result-6_{mode}.svg")

plt.show()