Sector-Coupling with pypsa#
Note
Also in this tutorial, you might want to refer to the PyPSA documentation: https://pypsa.readthedocs.io.
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 highspy
import pypsa
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from urllib.request import urlretrieve
plt.style.use("bmh")
Show code cell content
fn = "network-cem.nc"
url = "https://tubcloud.tu-berlin.de/s/kpWaraGc9LeaxLK/download/" + fn
urlretrieve(url, fn)
('network-cem.nc', <http.client.HTTPMessage at 0x7fcb729e2250>)
Previous Capacity Expansion Model#
To explore sector-coupling options with PyPSA, let’s load the capacity expansion model we built for the electricity system and add sector-coupling technologies and demands on top.
This example has single node for Germany and 4-hourly temporal resolution for a year. It has wind and solar solar generation, an OCGT generator as well as battery and hydrogen storage to supply a fixed electricity demand.
Some sector-coupling technologies have multiple ouputs (e.g. CHP plants producing heat and power). PyPSA can automatically handle links have more than one input (bus0)
and/or output (i.e. bus1, bus2, bus3) with a given efficieny (efficiency, efficiency2, efficiency3).
n = pypsa.Network("network-cem.nc")
WARNING:pypsa.io:Importing network from PyPSA version v0.21.3 while current version is v0.27.1. Read the release notes at https://pypsa.readthedocs.io/en/latest/release_notes.html to prepare your network for import.
INFO:pypsa.io:Imported network network-cem.nc has buses, carriers, generators, global_constraints, loads, storage_units
n
PyPSA Network
Components:
- Bus: 1
- Carrier: 6
- Generator: 4
- GlobalConstraint: 1
- Load: 1
- StorageUnit: 2
Snapshots: 2190
n.buses.index
Index(['Germany'], dtype='object', name='Bus')
n.generators.index
Index(['OCGT', 'onwind', 'offwind', 'solar'], dtype='object', name='Generator')
n.storage_units.index
Index(['battery storage', 'hydrogen storage underground'], dtype='object', name='StorageUnit')
Hydrogen Production#
The following example shows how to model the components of hydrogen storage separately, i.e. electrolysis, fuel cell and storage.
First, let’s remove the simplified hydrogen storage representation:
n.remove("StorageUnit", "hydrogen storage underground")
Add a separate Bus for the hydrogen energy carrier:
n.add("Bus", "hydrogen")
Add a Link for the hydrogen electrolysis:
n.add(
"Link",
"electrolysis",
bus0="electricity",
bus1="hydrogen",
carrier="electrolysis",
p_nom_extendable=True,
efficiency=0.7,
capital_cost=50e3, # €/MW/a
)
WARNING:pypsa.components:The bus name `electricity` given for bus0 of Link `electrolysis` does not appear in network.buses
Add a Link for the fuel cell which reconverts hydrogen to electricity:
n.add(
"Link",
"fuel cell",
bus0="hydrogen",
bus1="electricity",
carrier="fuel cell",
p_nom_extendable=True,
efficiency=0.5,
capital_cost=120e3, # €/MW/a
)
WARNING:pypsa.components:The bus name `electricity` given for bus1 of Link `fuel cell` does not appear in network.buses
Add a Store for the hydrogen storage:
n.add(
"Store",
"hydrogen storage",
bus="hydrogen",
carrier="hydrogen storage",
capital_cost=140, # €/MWh/a
e_nom_extendable=True,
e_cyclic=True, # cyclic state of charge
)
We can also add a hydrogen demand to the hydrogen bus.
In the example below, we add a constant hydrogen demand the size of the electricity demand.
p_set = n.loads_t.p_set["demand"].mean()
p_set
54671.88812785388
n.add("Load", "hydrogen demand", bus="hydrogen", carrier="hydrogen", p_set=p_set) # MW
Heat Demand#
For the heat demand, we create another bus and connect a load with the heat demand time series to it:
n.add("Bus", "heat")
url = "https://tubcloud.tu-berlin.de/s/mSkHERH8fJCKNXx/download/heat-load-example.csv"
p_set = pd.read_csv(url, index_col=0, parse_dates=True).squeeze()
p_set.head()
snapshot
2015-01-01 00:00:00 61726.043437
2015-01-01 04:00:00 108787.133591
2015-01-01 08:00:00 101508.988082
2015-01-01 12:00:00 90475.260586
2015-01-01 16:00:00 96307.755312
Name: 0, dtype: float64
n.add("Load", "heat demand", carrier="heat", bus="heat", p_set=p_set)
n.loads_t.p_set.div(1e3).plot(figsize=(12, 4), ylabel="GW")
<Axes: xlabel='snapshot', ylabel='GW'>
Heat pumps#
To model heat pumps, first we have to calculate the coefficient of performance (COP) profile based on the temperature profile of the heat source.
In the example below, we calculate the COP for an air-sourced heat pump with a sink temperature of 55° C and a population-weighted ambient temperature profile for Germany.
The heat pump performance is given by the following function:
where \(\Delta T = T_{sink} - T_{source}\).
def cop(t_source, t_sink=55):
delta_t = t_sink - t_source
return 6.81 - 0.121 * delta_t + 0.000630 * delta_t**2
url = "https://tubcloud.tu-berlin.de/s/S4jRAQMP5Te96jW/download/ninja_weather_country_DE_merra-2_population_weighted.csv"
temp = pd.read_csv(url, skiprows=2, index_col=0, parse_dates=True).loc[
"2015", "temperature"
][::4]
cop(temp).plot(figsize=(10, 2), ylabel="COP");
plt.scatter(temp, cop(temp))
plt.xlabel("temperature [°C]")
plt.ylabel("COP [-]")
Text(0, 0.5, 'COP [-]')
Once we have calculated the heat pump coefficient of performance, we can add the heat pump to the network as a Link. We use the parameter efficiency to incorporate the COP.
n.add(
"Link",
"heat pump",
carrier="heat pump",
bus0="electricity",
bus1="heat",
efficiency=cop(temp),
p_nom_extendable=True,
capital_cost=3e5, # €/MWe/a
)
WARNING:pypsa.components:The bus name `electricity` given for bus0 of Link `heat pump` does not appear in network.buses
Let’s also add a resistive heater as backup technology:
n.add(
"Link",
"resistive heater",
carrier="resistive heater",
bus0="electricity",
bus1="heat",
efficiency=0.9,
capital_cost=1e4, # €/MWe/a
p_nom_extendable=True,
)
WARNING:pypsa.components:The bus name `electricity` given for bus0 of Link `resistive heater` does not appear in network.buses
Combined Heat-and-Power (CHP)#
In the following, we are going to add gas-fired combined heat-and-power plants (CHPs). Today, these would use fossil gas, but in the example below we assume green methane with relatively high marginal costs. Since we have no other net emission technology, we can remove the CO\(_2\) limit.
n.remove("GlobalConstraint", "CO2Limit")
Then, we explicitly represent the energy carrier gas:
n.add("Bus", "gas", carrier="gas")
And add a Store of gas, which can be depleted (up to 100 TWh) with fuel costs of 150 €/MWh.
n.add(
"Store",
"gas storage",
carrier="gas storage",
e_initial=100e6, # MWh
e_nom=100e6, # MWh
bus="gas",
marginal_cost=150, # €/MWh_th
)
When we do this, we have to model the OCGT power plant as link which converts gas to electricity, not as generator.
n.remove("Generator", "OCGT")
n.add(
"Link",
"OCGT",
bus0="gas",
bus1="electricity",
carrier="OCGT",
p_nom_extendable=True,
capital_cost=20000, # €/MW/a
efficiency=0.4,
)
WARNING:pypsa.components:The bus name `electricity` given for bus1 of Link `OCGT` does not appear in network.buses
Next, we are going to add a combined heat-and-power (CHP) plant with fixed heat-power ratio (i.e. backpressure operation). If you want to model flexible heat-power ratios, have a look at this example: https://pypsa.readthedocs.io/en/latest/examples/power-to-gas-boiler-chp.html
n.add(
"Link",
"CHP",
bus0="gas",
bus1="electricity",
bus2="heat",
carrier="CHP",
p_nom_extendable=True,
capital_cost=40000,
efficiency=0.4,
efficiency2=0.4,
)
WARNING:pypsa.components:The bus name `electricity` given for bus1 of Link `CHP` does not appear in network.buses
Electric Vehicles#
To model electric vehicles, we first create another bus for the electric vehicles.
n.add("Bus", "EV", carrier="EV")
Then, we can attach the electricity consumption of electric vehicles to this bus:
url = "https://tubcloud.tu-berlin.de/s/9r5bMSbzzQiqG7H/download/electric-vehicle-profile-example.csv"
p_set = pd.read_csv(url, index_col=0, parse_dates=True).squeeze()
p_set.loc["2015-01-01"].div(1e3).plot(figsize=(4, 4), ylabel="GW")
<Axes: xlabel='snapshot', ylabel='GW'>
n.add("Load", "EV demand", bus="EV", carrier="EV demand", p_set=p_set)
Let’s have a quick look at how the heat, electricity, constant hydrogen and electric vehicle demands relate to each other:
n.loads_t.p_set.div(1e3).plot(figsize=(10, 3), ylabel="GW")
plt.axhline(
n.loads.loc["hydrogen demand", "p_set"] / 1e3, label="hydrogen demand", color="m"
)
plt.legend()
<matplotlib.legend.Legend at 0x7fcb6e2e1650>
The electric vehicles can only be charged when they are plugged-in. Below we load an availability profile telling us what share of electric vehicles is plugged-in at home – we only assume home charging in this example.
url = "https://tubcloud.tu-berlin.de/s/E3PBWPfYaWwCq7a/download/electric-vehicle-availability-example.csv"
availability_profile = pd.read_csv(url, index_col=0, parse_dates=True).squeeze()
availability_profile.loc["2015-01-01"].plot(ylim=(0, 1))
<Axes: xlabel='snapshot'>
Then, we can add a link for the electric vehicle charger using assumption about the number of EVs and their charging rates.
number_cars = 40e6 # number of EV cars
bev_charger_rate = 0.011 # 3-phase EV charger with 11 kW
p_nom = number_cars * bev_charger_rate
n.add(
"Link",
"EV charger",
bus0="electricity",
bus1="EV",
p_nom=p_nom,
carrier="EV charger",
p_max_pu=availability_profile,
efficiency=0.9,
)
WARNING:pypsa.components:The bus name `electricity` given for bus0 of Link `EV charger` does not appear in network.buses
We can also allow vehicle-to-grid operation (i.e. electric vehicles inject power into the grid):
n.add(
"Link",
"V2G",
bus0="EV",
bus1="electricity",
p_nom=p_nom,
carrier="V2G",
p_max_pu=availability_profile,
efficiency=0.9,
)
WARNING:pypsa.components:The bus name `electricity` given for bus1 of Link `V2G` does not appear in network.buses
The demand-side management potential we model as a store. This is not unlike a battery storage, but we impose additional constraints on when the store needs to be charged to a certain level (e.g. 75% full every morning).
bev_energy = 0.05 # average battery size of EV in MWh
bev_dsm_participants = 0.5 # share of cars that do smart charging
e_nom = number_cars * bev_energy * bev_dsm_participants
url = "https://tubcloud.tu-berlin.de/s/K62yACBRTrxLTia/download/dsm-profile-example.csv"
dsm_profile = pd.read_csv(url, index_col=0, parse_dates=True).squeeze()
dsm_profile.loc["2015-01-01"].plot(figsize=(5, 2), ylim=(0, 1))
<Axes: xlabel='snapshot'>
n.add(
"Store",
"EV DSM",
bus="EV",
carrier="EV battery",
e_cyclic=True, # state of charge at beginning = state of charge at the end
e_nom=e_nom,
e_min_pu=dsm_profile,
)
Then, we can solve the fully sector-coupled model altogether including electricity, passenger transport, hydrogen and heating.
n.optimize(solver_name="highs")
Show code cell output
WARNING:pypsa.components:The following links have bus0 which are not defined:
Index(['electrolysis', 'heat pump', 'resistive heater', 'EV charger'], dtype='object', name='Link')
WARNING:pypsa.components:The following links have bus1 which are not defined:
Index(['fuel cell', 'OCGT', 'CHP', 'V2G'], dtype='object', name='Link')
WARNING:pypsa.components:The following links have bus2 which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'OCGT'], dtype='object', name='Link')
WARNING:pypsa.components:The following links have bus0 which are not defined:
Index(['electrolysis', 'heat pump', 'resistive heater', 'EV charger'], dtype='object', name='Link')
WARNING:pypsa.components:The following links have bus1 which are not defined:
Index(['fuel cell', 'OCGT', 'CHP', 'V2G'], dtype='object', name='Link')
WARNING:pypsa.components:The following links have bus2 which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'OCGT'], dtype='object', name='Link')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 0%| | 0/27 [00:00<?, ?it/s]
Writing constraints.: 48%|████▊ | 13/27 [00:00<00:00, 105.25it/s]
Writing constraints.: 89%|████████▉ | 24/27 [00:00<00:00, 84.98it/s]
Writing constraints.: 100%|██████████| 27/27 [00:00<00:00, 67.60it/s]
Writing continuous variables.: 0%| | 0/11 [00:00<?, ?it/s]
Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 146.13it/s]
INFO:linopy.io: Writing time: 0.49s
INFO:linopy.solvers:Log file at /tmp/highs.log.
Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]
Copyright (c) 2023 HiGHS under MIT licence terms
Presolving model
47129 rows, 42760 cols, 119487 nonzeros
42749 rows, 34000 cols, 106347 nonzeros
Presolve : Reductions: rows 42749(-51432); columns 34000(-9811); elements 106347(-61243)
Solving the presolved LP
Using EKK dual simplex solver - serial
Iteration Objective Infeasibilities num(sum)
0 -5.1304687500e+06 Ph1: 8758(9.8535e+06); Du: 2189(5130.47) 0s
22622 7.2640708420e+10 Pr: 6310(4.06246e+11); Du: 0(1.7237e-06) 5s
30752 1.2131406462e+11 7s
Using EKK dual simplex solver - serial
Iteration Objective Infeasibilities num(sum)
30752 1.2129041226e+11 Pr: 1(8.59499e-06) 7s
30856 1.2129041006e+11 Pr: 0(0) 7s
30856 1.2129041006e+11 Pr: 0(0) 7s
Solving the original LP from the solution after postsolve
Using EKK dual simplex solver - serial
Iteration Objective Infeasibilities num(sum)
30856 1.2129041006e+11 Pr: 1432(0.0134378); Du: 0(2.19504e-08) 8s
32170 1.2129041006e+11 Pr: 2024(0.154159); Du: 0(4.15029e-08) 13s
32469 1.2129041006e+11 Pr: 1722(0.056773); Du: 0(5.42683e-08) 18s
32650 1.2129041006e+11 Pr: 1520(0.0124723); Du: 0(5.76357e-08) 23s
32800 1.2129041006e+11 Pr: 1297(0.00293065); Du: 0(5.7636e-08) 28s
32857 1.2129041006e+11 Pr: 0(0) 30s
Model status : Optimal
Simplex iterations: 32857
Objective value : 1.2129041006e+11
HiGHS run time : 30.21
n.statistics()
n.statistics()["Capital Expenditure"].div(1e9).sort_values().dropna().plot.bar(
ylabel="bn€/a", cmap="tab20c", figsize=(7, 3)
)
Exercises#
Here are a few things you can do to explore this system:
exclude one of transport sector, heating sector, hydrogen demand and investigate what impact it has on system operation and investments
vary the charging rates and availability profiles of electric vehicles