Single Node Sector Coupling

Note

You can download this example as a Jupyter notebook or start it in interactive mode.

Single Node Sector Coupling#

[1]:

from urllib.request import urlretrieve

import matplotlib.pyplot as plt
import pandas as pd

import pypsa

plt.style.use("bmh")

[2]:

fn = "network-cem.nc"
url = "https://tubcloud.tu-berlin.de/s/kpWaraGc9LeaxLK/download/" + fn
urlretrieve(url, fn)

[2]:

('network-cem.nc', <http.client.HTTPMessage at 0x7fb2958c6a50>)

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).

[3]:

n = pypsa.Network("network-cem.nc")

WARNING:pypsa.io:Importing network from PyPSA version v0.21.3 while current version is v0.31.0. 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

[4]:

n

[4]:

PyPSA Network
Components:
 - Bus: 1
 - Carrier: 6
 - Generator: 4
 - GlobalConstraint: 1
 - Load: 1
 - StorageUnit: 2
Snapshots: 2190

[5]:

n.buses.index

[5]:

Index(['Germany'], dtype='object', name='Bus')

[6]:

n.generators.index

[6]:

Index(['OCGT', 'onwind', 'offwind', 'solar'], dtype='object', name='Generator')

[7]:

n.storage_units.index

[7]:

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:

[8]:

n.remove("StorageUnit", "hydrogen storage underground")

Add a separate Bus for the hydrogen energy carrier:

[9]:

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

[9]:

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

Add a Link for the hydrogen electrolysis:

[10]:

n.add(
    "Link",
    "electrolysis",
    bus0="electricity",
    bus1="hydrogen",
    carrier="electrolysis",
    p_nom_extendable=True,
    efficiency=0.7,
    capital_cost=50e3,  # €/MW/a
)

WARNING:pypsa.io:The following Link have buses which are not defined:
Index(['electrolysis'], dtype='object')

[10]:

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

Add a Link for the fuel cell which reconverts hydrogen to electricity:

[11]:

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.io:The following Link have buses which are not defined:
Index(['fuel cell'], dtype='object')

[11]:

Index(['fuel cell'], dtype='object')

Add a Store for the hydrogen storage:

[12]:

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
)

[12]:

Index(['hydrogen storage'], dtype='object')

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.

[13]:

p_set = n.loads_t.p_set["demand"].mean()

[14]:

p_set

[14]:

54671.88812785388

[15]:

n.add("Load", "hydrogen demand", bus="hydrogen", carrier="hydrogen", p_set=p_set)  # MW

[15]:

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

Heat Demand#

For the heat demand, we create another bus and connect a load with the heat demand time series to it:

[16]:

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

[16]:

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

[17]:

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()

[18]:

p_set.head()

[18]:

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

[19]:

n.add("Load", "heat demand", carrier="heat", bus="heat", p_set=p_set)

[19]:

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

[20]:

n.loads_t.p_set.div(1e3).plot(figsize=(12, 4), ylabel="GW")

[20]:

<Axes: xlabel='snapshot', ylabel='GW'>
../_images/examples_sector-coupling-single-node_31_1.png

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:

\[COP(\Delta T) = 6.81 - 0.121 \Delta T + 0.00063^\Delta T^2\]

where \(\Delta T = T_{sink} - T_{source}\).

[21]:

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

[22]:

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]

[23]:

cop(temp).plot(figsize=(10, 2), ylabel="COP");
../_images/examples_sector-coupling-single-node_36_0.png 
[24]:

plt.scatter(temp, cop(temp))
plt.xlabel("temperature [°C]")
plt.ylabel("COP [-]")

[24]:

Text(0, 0.5, 'COP [-]')
../_images/examples_sector-coupling-single-node_37_1.png

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.

[25]:

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.io:The following Link have buses which are not defined:
Index(['heat pump'], dtype='object')

[25]:

Index(['heat pump'], dtype='object')

Let’s also add a resistive heater as backup technology:

[26]:

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.io:The following Link have buses which are not defined:
Index(['resistive heater'], dtype='object')

[26]:

Index(['resistive heater'], dtype='object')

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.

[27]:

n.remove("GlobalConstraint", "CO2Limit")

Then, we explicitly represent the energy carrier gas:

[28]:

n.add("Bus", "gas", carrier="gas")

[28]:

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

And add a Store of gas, which can be depleted (up to 100 TWh) with fuel costs of 150 €/MWh.

[29]:

n.add(
    "Store",
    "gas storage",
    carrier="gas storage",
    e_initial=100e6,  # MWh
    e_nom=100e6,  # MWh
    bus="gas",
    marginal_cost=150,  # €/MWh_th
)

[29]:

Index(['gas storage'], dtype='object')

When we do this, we have to model the OCGT power plant as link which converts gas to electricity, not as generator.

[30]:

n.remove("Generator", "OCGT")

[31]:

n.add(
    "Link",
    "OCGT",
    bus0="gas",
    bus1="electricity",
    carrier="OCGT",
    p_nom_extendable=True,
    capital_cost=20000,  # €/MW/a
    efficiency=0.4,
)

WARNING:pypsa.io:The following Link have buses which are not defined:
Index(['OCGT'], dtype='object')

[31]:

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

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

[32]:

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.io:The following Link have buses which are not defined:
Index(['CHP'], dtype='object')

[32]:

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

Electric Vehicles#

To model electric vehicles, we first create another bus for the electric vehicles.

[33]:

n.add("Bus", "EV", carrier="EV")

[33]:

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

Then, we can attach the electricity consumption of electric vehicles to this bus:

[34]:

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()

[35]:

p_set.loc["2015-01-01"].div(1e3).plot(figsize=(4, 4), ylabel="GW")

[35]:

<Axes: xlabel='snapshot', ylabel='GW'>
../_images/examples_sector-coupling-single-node_59_1.png 
[36]:

n.add("Load", "EV demand", bus="EV", carrier="EV demand", p_set=p_set)

[36]:

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

Let’s have a quick look at how the heat, electricity, constant hydrogen and electric vehicle demands relate to each other:

[37]:

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()

[37]:

<matplotlib.legend.Legend at 0x7fb28ff1bd10>
../_images/examples_sector-coupling-single-node_62_1.png

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.

[38]:

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()

[39]:

availability_profile.loc["2015-01-01"].plot(ylim=(0, 1))

[39]:

<Axes: xlabel='snapshot'>
../_images/examples_sector-coupling-single-node_65_1.png

Then, we can add a link for the electric vehicle charger using assumption about the number of EVs and their charging rates.

[40]:

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

[41]:

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.io:The following Link have buses which are not defined:
Index(['EV charger'], dtype='object')

[41]:

Index(['EV charger'], dtype='object')

We can also allow vehicle-to-grid operation (i.e. electric vehicles inject power into the grid):

[42]:

n.add(
    "Link",
    "V2G",
    bus0="EV",
    bus1="electricity",
    p_nom=p_nom,
    carrier="V2G",
    p_max_pu=availability_profile,
    efficiency=0.9,
)

WARNING:pypsa.io:The following Link have buses which are not defined:
Index(['V2G'], dtype='object')

[42]:

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

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).

[43]:

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

[44]:

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()

[45]:

dsm_profile.loc["2015-01-01"].plot(figsize=(5, 2), ylim=(0, 1))

[45]:

<Axes: xlabel='snapshot'>
../_images/examples_sector-coupling-single-node_74_1.png 
[46]:

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,
)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[46], line 1
----> 1 n.add(
      2     "Store",
      3     "EV DSM",
      4     bus="EV",
      5     carrier="EV battery",
      6     e_cyclic=True,  # state of charge at beginning = state of charge at the end
      7     e_nom=e_nom,
      8     e_min_pu=dsm_profile,
      9 )

File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.12/site-packages/pypsa/components.py:986, in Network.add(self, class_name, name, suffix, overwrite, **kwargs)
    984 if isinstance(v, pd.Series) and single_component:
    985     if not v.index.equals(self.snapshots):
--> 986         raise ValueError(msg.format(f"Series {k}", "network snapshots"))
    987 elif isinstance(v, pd.Series):
    988     # Cast names index to string + suffix
    989     v = v.rename(
    990         index=lambda s: str(s) if str(s).endswith(suffix) else s + suffix
    991     )

ValueError: Series e_min_pu has an index which does not align with the passed network snapshots.

Then, we can solve the fully sector-coupled model altogether including electricity, passenger transport, hydrogen and heating.

[47]:

n.optimize(solver_name="highs")

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['Germany', 'hydrogen', 'heat', 'gas', 'EV'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand', 'heat demand', 'EV demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following links have buses which are not defined:
Index(['electrolysis', 'heat pump', 'resistive heater', 'EV charger'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following links have buses which are not defined:
Index(['fuel cell', 'OCGT', 'CHP', 'V2G'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'CHP',
       'EV charger', 'V2G'],
      dtype='object', name='Link')
WARNING:pypsa.consistency:Encountered nan's in static data for columns ['efficiency2'] of component 'Link'.
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage', 'gas storage'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['Germany', 'hydrogen', 'heat', 'gas', 'EV'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand', 'heat demand', 'EV demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following links have buses which are not defined:
Index(['electrolysis', 'heat pump', 'resistive heater', 'EV charger'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following links have buses which are not defined:
Index(['fuel cell', 'OCGT', 'CHP', 'V2G'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'CHP',
       'EV charger', 'V2G'],
      dtype='object', name='Link')
WARNING:pypsa.consistency:Encountered nan's in static data for columns ['efficiency2'] of component 'Link'.
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage', 'gas storage'], dtype='object', name='Store')
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.12/site-packages/linopy/common.py:147: UserWarning: coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
  warn(
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 29/29 [00:00<00:00, 49.57it/s]
Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 106.56it/s]
INFO:linopy.io: Writing time: 0.72s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 39431 primals, 87614 duals
Objective: 1.21e+11
Solver model: available
Solver message: optimal


Coefficient ranges:
  Matrix [1e-04, 6e+00]
  Cost   [4e-02, 3e+05]
  Bound  [0e+00, 0e+00]
  RHS    [2e+03, 1e+08]
Presolving model
42749 rows, 34000 cols, 106347 nonzeros  0s
38370 rows, 29621 cols, 97589 nonzeros  0s
Presolve : Reductions: rows 38370(-49244); columns 29621(-9810); elements 97589(-63431)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Ph1: 0(0) 0s
      18887     5.4939739200e+10 Pr: 6549(2.18291e+13); Du: 0(3.03056e-06) 5s
      30128     1.2129041006e+11 Pr: 0(0) 10s
      30128     1.2129041006e+11 Pr: 0(0) 10s
Solving the original LP from the solution after postsolve
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
      30128     1.2129041006e+11 Pr: 1997(0.169348); Du: 0(1.51501e-08) 12s
      34004     1.2129041006e+11 Pr: 2184(0.0318987); Du: 0(4.19575e-08) 18s
      38830     1.2129041006e+11 Pr: 2183(0.0429197); Du: 0(4.24779e-08) 23s
      43234     1.2129041006e+11 Pr: 2196(0.0363101); Du: 0(5.14587e-08) 28s
      46944     1.2129041006e+11 Pr: 2183(0.0271657); Du: 0(6.55028e-08) 34s
      48743     1.2129041006e+11 Pr: 0(0) 36s
Model   status      : Optimal
Simplex   iterations: 48743
Objective value     :  1.2129041006e+11
HiGHS run time      :         36.34
Writing the solution to /tmp/linopy-solve-qvbjs5kh.sol

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, Link-fix-p-lower, Link-fix-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-fix-e-lower, Store-fix-e-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance, Store-energy_balance were not assigned to the network.

[47]:

('ok', 'optimal')

[48]:

n.statistics()

[48]:
Optimal Capacity Installed Capacity Supply Withdrawal Energy Balance Transmission Capacity Factor Curtailment Capital Expenditure Operational Expenditure Revenue Market Value
Generator offwind 2.281124e+04 0.0 3.196011e+07 0.000000e+00 3.196011e+07 0.0 0.159939 4.053516e+07 3.606999e+09 6.392023e+05 3.607639e+09 112.879405
onwind 3.437762e+05 0.0 3.866098e+07 0.000000e+00 3.866098e+07 0.0 0.012838 5.809284e+08 3.303204e+10 5.219232e+07 3.308423e+10 855.752513
solar 7.594874e+05 0.0 4.245393e+08 0.000000e+00 4.245393e+08 0.0 0.063811 4.016830e+08 3.532292e+10 4.245393e+06 3.532716e+10 83.212930
Link EV charger 4.400000e+05 440000.0 1.242776e+08 1.380862e+08 -1.380862e+07 0.0 0.035826 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 NaN
V2G 4.400000e+05 440000.0 0.000000e+00 0.000000e+00 0.000000e+00 0.0 0.000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000
electrolysis 7.810270e+04 0.0 4.789257e+08 6.841796e+08 -2.052539e+08 0.0 1.000000 0.000000e+00 3.905135e+09 0.000000e+00 3.905135e+09 8.153947
resistive heater 2.277568e+05 0.0 5.936692e+08 6.596325e+08 -6.596325e+07 0.0 0.330618 0.000000e+00 2.277568e+09 0.000000e+00 2.277568e+09 3.836426
Load - 0.000000e+00 0.0 0.000000e+00 4.789257e+08 -4.789257e+08 0.0 NaN 0.000000e+00 0.000000e+00 0.000000e+00 -1.151077e+11 NaN
EV demand 0.000000e+00 0.0 0.000000e+00 1.242776e+08 -1.242776e+08 0.0 NaN 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000
heat 0.000000e+00 0.0 0.000000e+00 5.936692e+08 -5.936692e+08 0.0 NaN 0.000000e+00 0.000000e+00 0.000000e+00 -2.277568e+09 NaN
hydrogen 0.000000e+00 0.0 0.000000e+00 4.789257e+08 -4.789257e+08 0.0 NaN 0.000000e+00 0.000000e+00 0.000000e+00 -3.905135e+09 NaN
StorageUnit battery storage 4.468159e+05 0.0 1.908404e+08 2.070751e+08 -1.623469e+07 0.0 0.101662 3.930342e+09 4.308868e+10 0.000000e+00 4.308868e+10 225.783806
Store gas storage 1.000000e+08 100000000.0 0.000000e+00 0.000000e+00 0.000000e+00 0.0 1.000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000
hydrogen storage 0.000000e+00 0.0 0.000000e+00 0.000000e+00 0.000000e+00 0.0 0.000000 0.000000e+00 9.000000e-05 0.000000e+00 0.000000e+00 0.000000 
[49]:

n.statistics()["Capital Expenditure"].div(1e9).sort_values().dropna().plot.bar(
    ylabel="bn€/a", cmap="tab20c", figsize=(7, 3)
)

[49]:

<Axes: ylabel='bn€/a'>
../_images/examples_sector-coupling-single-node_79_1.png 
[ ]:
Contents