Advanced models in pypsa: capacity expansion planning with sector coupling

Advanced models in `pypsa`: capacity expansion planning with sector coupling#

Note

If you have not yet set up Python on your computer, you can execute this tutorial in your browser via Google Colab. Click on the rocket in the top right corner and launch “Colab”. If that doesn’t work download the .ipynb file and import it in Google Colab.

Then install the following packages by executing the following command in a Jupyter cell at the top of the notebook.

!pip install pypsa pandas numpy matplotlib plotly highspy

Problem description#

To explore capacity expansion problem with sector-coupling options, let’s model an “greenfield” energy system model with:

electricity demand (here: Germany 2015)
several technologies for electricity generation (wind, solar, and gas for peak load)
hydrogen consumption (e.g., an offtaker in the industrial sector)
hydrogen production from electrolysis
hydrogen storage
hydrogen fuel cell

import pypsa
import pandas as pd

pd.options.plotting.backend = "plotly"

Prepare technology data#

At TU Berlin, we maintain a database (PyPSA/technology-data) which collects assumptions and projections for energy system technologies (such as costs, efficiencies, lifetimes, etc.) for given years, which we use for our research.

Reading this data into a useable pandas.DataFrame requires some pre-processing (e.g. converting units, setting defaults, re-arranging dimensions):

YEAR = 2030
url = f"https://raw.githubusercontent.com/PyPSA/technology-data/master/outputs/costs_{YEAR}.csv"
costs = pd.read_csv(url, index_col=[0, 1])

costs.loc[costs.unit.str.contains("/kW"), "value"] *= 1e3
costs.unit = costs.unit.str.replace("/kW", "/MW")

defaults = {
    "FOM": 0,
    "VOM": 0,
    "efficiency": 1,
    "fuel": 0,
    "investment": 0,
    "lifetime": 25,
    "CO2 intensity": 0,
    "discount rate": 0.07,
}
costs = costs.value.unstack().fillna(defaults)

costs.at["OCGT", "fuel"] = costs.at["gas", "fuel"]
# costs.at["CCGT", "fuel"] = costs.at["gas", "fuel"]
costs.at["OCGT", "CO2 intensity"] = costs.at["gas", "CO2 intensity"]
# costs.at["CCGT", "CO2 intensity"] = costs.at["gas", "CO2 intensity"]

# costs

Let’s also write a small utility function that calculates the annuity to annualise investment costs. The formula is

\[ a(r, n) = \frac{r}{1-(1+r)^{-n}} \]

where \(r\) is the discount rate and \(n\) is the lifetime.

def annuity_factor(r, n):
    return r / (1.0 - 1.0 / (1.0 + r) ** n)

The resulting annuity_factor (often called capital recovery factor) represents the factor that, when multiplied by the initial capital cost, yields the equivalent annual cost accounting for the time value of money:

annuity_factor(0.07, 20)

0.09439292574325567

Based on this, we can calculate the short-term marginal generation costs (€/MWh)

costs["marginal_cost"] = costs["VOM"] + costs["fuel"] / costs["efficiency"]

and the annualised investment costs (capital_cost in PyPSA terms, €/MW/a):

annuity = costs.apply(lambda x: annuity_factor(0.07, x["lifetime"]), axis=1)

costs["capital_cost"] = (annuity + costs["FOM"] / 100) * costs["investment"]

Prepare load and renewable generation time-series#

We are also going to need some time series for wind, solar and load.

For now, we are going to recycle the time series we used in the pandas tutorial.

url = (
    "https://tubcloud.tu-berlin.de/s/pKttFadrbTKSJKF/download/time-series-lecture-2.csv"
)
ts = pd.read_csv(url, index_col=0, parse_dates=True)

# fallback if url is not available
# ts = pd.read_csv("../data/time-series-lecture-2.csv", index_col=0, parse_dates=True)

ts.head()

	load	onwind	offwind	prices
2015-01-01 00:00:00	41.151	0.1566	0.7030	NaN
2015-01-01 01:00:00	40.135	0.1659	0.6875	NaN
2015-01-01 02:00:00	39.106	0.1746	0.6535	NaN
2015-01-01 03:00:00	38.765	0.1745	0.6803	NaN
2015-01-01 04:00:00	38.941	0.1826	0.7272	NaN

Let’s also convert the load time series from GW to MW, the base unit of PyPSA:

ts.load *= 1e3

Optionally, we can downscale temporal resolution of the time series to save some computation time:

# here we sample only every other hour:
resolution = 4
ts = ts.resample(f"{resolution}h").first()

Build our energy system to simulate#

As always, let’s initialize an empty network:

n = pypsa.Network()

Then, we add a single electricity bus…

n.add("Bus", "electricity")

Index(['electricity'], dtype='object')

…and tell the pypsa.Network object n what the snapshots of the model will be using the utility function n.set_snapshots().

n.set_snapshots(ts.index)

n.snapshots

DatetimeIndex(['2015-01-01 00:00:00', '2015-01-01 04:00:00',
               '2015-01-01 08:00:00', '2015-01-01 12:00:00',
               '2015-01-01 16:00:00', '2015-01-01 20:00:00',
               '2015-01-02 00:00:00', '2015-01-02 04:00:00',
               '2015-01-02 08:00:00', '2015-01-02 12:00:00',
               ...
               '2015-12-30 08:00:00', '2015-12-30 12:00:00',
               '2015-12-30 16:00:00', '2015-12-30 20:00:00',
               '2015-12-31 00:00:00', '2015-12-31 04:00:00',
               '2015-12-31 08:00:00', '2015-12-31 12:00:00',
               '2015-12-31 16:00:00', '2015-12-31 20:00:00'],
              dtype='datetime64[ns]', name='snapshot', length=2190, freq='4h')

If we resampled the time series above, we need to adjust the weighting of the snapshots (i.e. how many hours they represent). We can do that with n.snapshot_weightings:

n.snapshot_weightings.head(3)

	objective	stores	generators
snapshot
2015-01-01 00:00:00	1.0	1.0	1.0
2015-01-01 04:00:00	1.0	1.0	1.0
2015-01-01 08:00:00	1.0	1.0	1.0

n.snapshot_weightings.loc[:, :] = resolution

n.snapshot_weightings.head(3)

	objective	stores	generators
snapshot
2015-01-01 00:00:00	4.0	4.0	4.0
2015-01-01 04:00:00	4.0	4.0	4.0
2015-01-01 08:00:00	4.0	4.0	4.0

Adding components for electricity part#

Then, we add all the technologies we are going to include as carriers.

carriers = [
    "onwind",
    "offwind",
    "solar",
    "OCGT",
    "hydrogen storage underground",
]

n.add(
    "Carrier",
    carriers,
    color=["dodgerblue", "aquamarine", "gold", "indianred", "magenta"],
    co2_emissions=[costs.at[c, "CO2 intensity"] for c in carriers],
)

Index(['onwind', 'offwind', 'solar', 'OCGT', 'hydrogen storage underground'], dtype='object')

Next, we add the demand time series to the model.

n.add(
    "Load",
    "demand",
    bus="electricity",
    p_set=ts.load,
)

Index(['demand'], dtype='object')

Let’s have a check whether the data was read-in correctly.

n.loads_t.p_set.plot()