Demand#

This tutorial covers electricity demand modeling in PyPSA-GB, including historical profiles, future projections, and spatial disaggregation.

What You’ll Learn#

  • Demand data sources (ESPENI, eLOAD, DESSTINEE)

  • Temporal profiles and patterns

  • Spatial demand distribution

  • Future demand scaling with FES

  • Demand-side flexibility

1. Setup#

[1]:
import pypsa
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import warnings
from pyproj import Transformer
from _map_utils import prepare_map_network, explore_network_map

warnings.filterwarnings('ignore')
plt.style.use('seaborn-v0_8-whitegrid')
plt.rcParams['figure.figsize'] = [12, 6]
plt.rcParams['figure.dpi'] = 100

print(f"PyPSA version: {pypsa.__version__}")
PyPSA version: 1.0.7

2. Demand Data Sources#

PyPSA-GB supports multiple demand data sources:

Source

Type

Resolution

Use Case

ESPENI

Historical

Half-hourly, national

Historical scenarios

eLOAD

Synthetic

Hourly, sectoral

Future projections

DESSTINEE

Synthetic

Hourly

Alternative future profiles

FES

Projections

Annual totals

Capacity scaling

3. Load a Network with Demand#

[2]:
# Load a historical network
n = pypsa.Network("../../../resources/network/Historical_2015_reduced_solved.nc")

print("Network loaded")
print(f"  Loads: {len(n.loads)}")
print(f"  Snapshots: {len(n.snapshots)}")
print(f"  Period: {n.snapshots[0]} to {n.snapshots[-1]}")
INFO:pypsa.network.io:Imported network 'Historical_2015_reduced (Full)' has buses, carriers, generators, lines, links, loads, storage_units, sub_networks
Network loaded
  Loads: 32
  Snapshots: 168
  Period: 2015-01-01 00:00:00 to 2015-01-07 23:00:00

4. Demand Structure#

4.1 Loads in PyPSA#

[3]:
# Load components
print("Load DataFrame:")
n.loads.head()
Load DataFrame:
[3]:
bus carrier type p_set q_set sign active
name
load_Beauly Beauly 0.0 0.0 -1.0 True
load_Peterhead Peterhead 0.0 0.0 -1.0 True
load_Errochty Errochty 0.0 0.0 -1.0 True
load_Denny/Bonnybridge Denny/Bonnybridge 0.0 0.0 -1.0 True
load_Neilston Neilston 0.0 0.0 -1.0 True
[4]:
# Demand time series
print("Demand time series shape:", n.loads_t.p_set.shape)
print(f"\nSample demand values (MW):")
n.loads_t.p_set.head()
Demand time series shape: (168, 31)

Sample demand values (MW):
[4]:
name load_Beauly load_Peterhead load_Errochty load_Denny/Bonnybridge load_Neilston load_Strathaven load_Torness load_Eccles load_Harker load_Stella West ... load_Sundon/East Claydon load_Melksham load_Bramley load_London load_Kemsley load_Sellindge load_Lovedean load_S.W.Penisula EU_demand_IFA EU_demand_East West Interconnector
snapshot
2015-01-01 00:00:00 125.398540 115.993649 34.484598 250.797080 250.797080 75.239124 106.588759 18.809781 172.422992 1686.610361 ... 1319.819632 2752.497950 1225.770727 4790.224222 1147.396640 87.778978 742.986349 2580.074957 0.0 120.0
2015-01-01 01:00:00 125.460546 116.051005 34.501650 250.921092 250.921092 75.276328 106.641464 18.819082 172.508251 1687.444344 ... 1320.472247 2753.858986 1226.376838 4792.592859 1147.963996 87.822382 743.353735 2581.350735 0.0 120.0
2015-01-01 02:00:00 120.162016 111.149865 33.044554 240.324032 240.324032 72.097210 102.137714 18.024302 165.222772 1616.179118 ... 1264.705221 2637.556256 1174.583708 4590.189019 1099.482448 84.113411 711.959946 2472.333483 0.0 8.0
2015-01-01 03:00:00 114.773477 106.165467 31.562706 229.546955 229.546955 68.864086 97.557456 17.216022 157.813531 1543.703270 ... 1207.990849 2519.277828 1121.910741 4384.346835 1050.177318 80.341434 680.032853 2361.464296 0.0 0.0
2015-01-01 04:00:00 108.474847 100.339234 29.830583 216.949695 216.949695 65.084908 92.203620 16.271227 149.152915 1458.986699 ... 1141.697770 2381.022902 1060.341634 4143.739174 992.544854 75.932393 642.713471 2231.869987 1022.0 0.0

5 rows × 31 columns

4.2 System-Wide Demand#

[5]:
# Total system demand
total_demand = n.loads_t.p_set.sum(axis=1) / 1000  # GW

print("System Demand Statistics (GW):")
print(f"  Mean: {total_demand.mean():.2f}")
print(f"  Peak: {total_demand.max():.2f}")
print(f"  Minimum: {total_demand.min():.2f}")
print(f"  Standard Dev: {total_demand.std():.2f}")

# Total energy
hours = len(n.snapshots)
energy_twh = total_demand.sum() / 1000  # TWh
print(f"\nTotal Energy: {energy_twh:.2f} TWh over {hours} hours")
System Demand Statistics (GW):
  Mean: 38.33
  Peak: 52.21
  Minimum: 27.11
  Standard Dev: 7.17

Total Energy: 6.44 TWh over 168 hours
[6]:
# Plot demand profile
fig, ax = plt.subplots(figsize=(14, 5))

ax.plot(total_demand.index, total_demand.values, linewidth=1.5, color='blue')
ax.fill_between(total_demand.index, total_demand.values, alpha=0.3)

ax.axhline(y=total_demand.mean(), color='red', linestyle='--', label=f'Mean: {total_demand.mean():.1f} GW')

ax.set_ylabel('Demand (GW)')
ax.set_xlabel('Time')
ax.set_title('System Electricity Demand')
ax.legend()

plt.tight_layout()
plt.show()
../_images/tutorials_6-demand_11_0.png

5. Temporal Patterns#

5.1 Daily Pattern#

[7]:
# Hourly pattern
demand_df = total_demand.to_frame('demand')
demand_df['hour'] = demand_df.index.hour

hourly_pattern = demand_df.groupby('hour')['demand'].mean()

fig, ax = plt.subplots(figsize=(10, 5))

ax.bar(hourly_pattern.index, hourly_pattern.values, color='steelblue', edgecolor='black')
ax.set_xlabel('Hour of Day')
ax.set_ylabel('Average Demand (GW)')
ax.set_title('Average Daily Demand Pattern')
ax.set_xticks(range(0, 24))

plt.tight_layout()
plt.show()

print(f"Morning peak: {hourly_pattern.loc[7:9].idxmax()}:00 ({hourly_pattern.loc[7:9].max():.2f} GW)")
print(f"Evening peak: {hourly_pattern.loc[17:21].idxmax()}:00 ({hourly_pattern.loc[17:21].max():.2f} GW)")
print(f"Night trough: {hourly_pattern.idxmin()}:00 ({hourly_pattern.min():.2f} GW)")
../_images/tutorials_6-demand_13_0.png
Morning peak: 9:00 (39.86 GW)
Evening peak: 17:00 (48.25 GW)
Night trough: 4:00 (29.30 GW)

5.2 Load Duration Curve#

[8]:
# Load duration curve
sorted_demand = total_demand.sort_values(ascending=False).values
hours = np.arange(1, len(sorted_demand) + 1)

fig, ax = plt.subplots(figsize=(12, 6))

ax.plot(hours, sorted_demand, linewidth=2, color='blue')
ax.fill_between(hours, sorted_demand, alpha=0.3)

# Mark key points
ax.axhline(y=total_demand.mean(), color='red', linestyle='--', label=f'Mean: {total_demand.mean():.1f} GW')
ax.axhline(y=total_demand.min(), color='green', linestyle=':', label=f'Min: {total_demand.min():.1f} GW')

ax.set_xlabel('Hours')
ax.set_ylabel('Demand (GW)')
ax.set_title('Load Duration Curve')
ax.legend()

plt.tight_layout()
plt.show()

# Load factor
load_factor = total_demand.mean() / total_demand.max() * 100
print(f"Load Factor: {load_factor:.1f}%")
../_images/tutorials_6-demand_15_0.png
Load Factor: 73.4%

6. Spatial Distribution#

6.1 Demand by Bus#

[9]:
# Average demand at each bus
bus_demand = n.loads_t.p_set.mean() / 1000  # GW

print("Top 10 Buses by Demand (GW):")
print(bus_demand.sort_values(ascending=False).head(10).round(3).to_string())
Top 10 Buses by Demand (GW):
name
load_London                    5.820
load_Melksham                  3.344
load_S.W.Penisula              3.135
load_Daines                    3.013
load_Feckenham                 2.727
load_Deeside                   2.643
load_Stella West               2.049
load_Th. Marsh/Stocksbridge    2.007
load_Sundon/East Claydon       1.603
load_Ratcliffe                 1.554
[10]:
# Interactive map showing demand distribution
if len(bus_demand) > 0:
    map_n = prepare_map_network(n)

    # Vectorised calculation of average demand per load (MW)
    mean_per_load = n.loads_t.p_set.mean()

    # Keep only loads that actually have a timeseries
    loads_with_ts = n.loads.loc[n.loads.index.intersection(mean_per_load.index)]

    # Join load-to-bus mapping with mean demand and aggregate per bus
    df = loads_with_ts[['bus']].join(mean_per_load.rename('mean_MW'))
    bus_demand_by_bus = df.groupby('bus')['mean_MW'].sum()

    # Reindex to include all buses and fill missing values with 0
    bus_demand_by_bus = bus_demand_by_bus.reindex(n.buses.index).fillna(0.0)

    # Add demand to buses for tooltip display (units: MW)
    map_n.buses["demand_MW"] = bus_demand_by_bus

    # Report any loads that lacked a timeseries (examples shown)
    missing_loads = set(n.loads.index) - set(mean_per_load.index)
    if len(missing_loads) > 0:
        print(f"Warning: {len(missing_loads)} loads have no timeseries and were ignored (examples: {list(missing_loads)[:5]})")

    print("Spatial Demand Distribution (Bus size = average demand)")
    print("lon range:", float(map_n.buses.x.min()), float(map_n.buses.x.max()))
    print("lat range:", float(map_n.buses.y.min()), float(map_n.buses.y.max()))

    # Interactive network map with demand-based bus sizing
    m = map_n.plot.explore(
        map_style="light",
        tooltip=True,
        bus_size=map_n.buses["demand_MW"].clip(lower=10),
        bus_size_factor=0.05,
        branch_width_factor=1.0,
        bus_columns=["demand_MW", "v_nom"],
    )
    display(m)
Warning: 1 loads have no timeseries and were ignored (examples: ['EU_demand_Britned'])
Spatial Demand Distribution (Bus size = average demand)
lon range: -6.09249939157784 4.058099999999997
lat range: 50.912041859956986 57.484467069924335

7. Future Demand Projections#

7.1 Compare Historical vs Future#

[11]:
# Try to load a future scenario for comparison
try:
    n_future = pypsa.Network("../../../resources/network/HT35_clustered_solved.nc")
    future_demand = n_future.loads_t.p_set.sum(axis=1) / 1000

    print("Demand Comparison:")
    print(f"\n  Historical 2015:")
    print(f"    Mean: {total_demand.mean():.2f} GW")
    print(f"    Peak: {total_demand.max():.2f} GW")

    print(f"\n  Future 2035 HT:")
    print(f"    Mean: {future_demand.mean():.2f} GW")
    print(f"    Peak: {future_demand.max():.2f} GW")

    print(f"\n  Change:")
    print(f"    Mean: {(future_demand.mean()/total_demand.mean() - 1)*100:+.1f}%")
    print(f"    Peak: {(future_demand.max()/total_demand.max() - 1)*100:+.1f}%")
except FileNotFoundError:
    print("Future network not found")
    n_future = None
INFO:pypsa.network.io:Imported network 'HT35_clustered (Clustered)' has buses, carriers, generators, lines, links, loads, storage_units, stores, sub_networks
Demand Comparison:

  Historical 2015:
    Mean: 38.33 GW
    Peak: 52.21 GW

  Future 2035 HT:
    Mean: 38.95 GW
    Peak: 48.43 GW

  Change:
    Mean: +1.6%
    Peak: -7.2%
[12]:
# Compare load duration curves
if n_future is not None:
    fig, ax = plt.subplots(figsize=(12, 6))

    # Historical
    sorted_hist = total_demand.sort_values(ascending=False).values
    hours_hist = np.arange(1, len(sorted_hist) + 1) / len(sorted_hist) * 100

    # Future
    sorted_future = future_demand.sort_values(ascending=False).values
    hours_future = np.arange(1, len(sorted_future) + 1) / len(sorted_future) * 100

    ax.plot(hours_hist, sorted_hist, linewidth=2, label='2015 Historical', color='blue')
    ax.plot(hours_future, sorted_future, linewidth=2, label='2035 HT', color='green')

    ax.set_xlabel('% of Time')
    ax.set_ylabel('Demand (GW)')
    ax.set_title('Load Duration Curve Comparison')
    ax.legend()

    plt.tight_layout()
    plt.show()
../_images/tutorials_6-demand_21_0.png

8. Demand vs Generation Balance#

[13]:
# Generation vs demand
generation = n.generators_t.p.sum(axis=1) / 1000  # GW

fig, ax = plt.subplots(figsize=(14, 5))

ax.plot(total_demand.index, total_demand.values, label='Demand', linewidth=2, color='blue')
ax.plot(generation.index, generation.values, label='Generation', linewidth=2, color='green', alpha=0.7)

ax.set_ylabel('Power (GW)')
ax.set_xlabel('Time')
ax.set_title('Generation vs Demand Balance')
ax.legend()

plt.tight_layout()
plt.show()

# Balance check
imbalance = (generation - total_demand).abs()
print(f"Max imbalance: {imbalance.max()*1000:.1f} MW")
print(f"Average imbalance: {imbalance.mean()*1000:.1f} MW")
../_images/tutorials_6-demand_23_0.png
Max imbalance: 300.0 MW
Average imbalance: 78.2 MW

9. Demand Time Series Configuration#

Demand time series can be configured in scenarios:

[14]:
demand_config = """
# ESPENI (default for historical)
demand_timeseries: "ESPENI"    # Half-hourly national demand

# eLOAD (synthetic profiles)
demand_timeseries: "eload"
profile_year: 2050             # 2010 or 2050 profile year

# DESSTINEE (alternative synthetic)
demand_timeseries: "desstinee"

# Demand disaggregation (future)
demand_disaggregation:
  enabled: true
  components: ["heat_pumps", "EVs"]
"""
print("Demand Configuration Options:")
print(demand_config)
Demand Configuration Options:

# ESPENI (default for historical)
demand_timeseries: "ESPENI"    # Half-hourly national demand

# eLOAD (synthetic profiles)
demand_timeseries: "eload"
profile_year: 2050             # 2010 or 2050 profile year

# DESSTINEE (alternative synthetic)
demand_timeseries: "desstinee"

# Demand disaggregation (future)
demand_disaggregation:
  enabled: true
  components: ["heat_pumps", "EVs"]