Future Scenario: Electric Engagement (2050)#
This tutorial explores the Electric Engagement 2050 scenario from NESO’s FES - an ambitious pathway achieving Net Zero with extensive electrification across transport, heating, and industry. This represents the most aggressive decarbonization trajectory.
What You’ll Learn#
Modeling highly electrified future energy systems
Analyzing massive renewable capacity deployment
Understanding seasonal storage and flexibility requirements
Evaluating hydrogen’s role in long-duration energy storage
Identifying critical infrastructure and investment needs
Prerequisites#
Run the workflow to generate the solved network:
snakemake resources/network/EE50_clustered_solved.nc -j 4
Scenario Overview#
Parameter |
Value |
|---|---|
Modelled Year |
2050 |
FES Pathway |
Electric Engagement |
Network Model |
Clustered (~100 buses) |
Data Source |
FES 2024 projections |
Solve Period |
Representative week |
Key Features |
100+ GW wind/solar, large H2 system, high demand |
[1]:
# Import required libraries
import os
import sys
import pypsa
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import plotly.express as px
import plotly.graph_objects as go
import plotly.io as pio
from plotly.subplots import make_subplots
from pyproj import Transformer
from pathlib import Path
# Configure matplotlib and plotly
plt.rcParams.update({'font.size': 12})
plt.style.use('ggplot')
%matplotlib inline
# Configure Plotly to output HTML for Sphinx/nbsphinx compatibility
pio.renderers.default = 'notebook_connected'
# Suppress warnings
import warnings
warnings.filterwarnings('ignore')
print('✓ Libraries imported successfully')
✓ Libraries imported successfully
[2]:
# Load solved network (resolve path robustly from different working directories)
network_rel = Path('resources/network/EE50_clustered_solved.nc')
network_path = network_rel
# If the relative path doesn't exist from the current cwd, search upward for the repository root
if not network_path.exists():
for parent in [Path.cwd()] + list(Path.cwd().parents)[:6]:
candidate = parent / network_rel
if candidate.exists():
network_path = candidate
break
else:
# leave as relative; we'll raise a clear error below if not found
network_path = network_rel
print(f'Loading network from: {network_path}')
if not network_path.exists():
raise FileNotFoundError(f'Network file not found. Tried: {network_path}.\\nMake sure "resources/network/EE50_clustered_solved.nc" exists relative to the repository root.')
network = pypsa.Network(str(network_path))
# Derive a scenario id from the network filename when not provided
scenario_id = 'EE50_clustered'
print(f'✓ Network loaded successfully')
print(f' Buses: {len(network.buses)}')
print(f' Generators: {len(network.generators)}')
print(f' Lines: {len(network.lines)}')
print(f' Storage units: {len(network.storage_units)}')
print(f' Timesteps: {len(network.snapshots)}')
print(f' Period: {network.snapshots[0]} to {network.snapshots[-1]}')
Loading network from: c:\Users\alyden\OneDrive - University of Edinburgh\Python\PyPSA-GB v0.0.1\resources\network\EE50_clustered_solved.nc
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'EE50_clustered (Clustered)' has buses, carriers, generators, lines, links, loads, storage_units, stores, sub_networks
✓ Network loaded successfully
Buses: 110
Generators: 5147
Lines: 202
Storage units: 779
Timesteps: 168
Period: 2050-02-01 00:00:00 to 2050-02-07 23:00:00
Network Summary#
[3]:
# Calculate generation by carrier
p_by_carrier = network.generators_t.p.groupby(
network.generators.carrier, axis=1).sum()
# Add storage output
storage_by_carrier = network.storage_units_t.p.groupby(
network.storage_units.carrier, axis=1).sum()
storage_by_carrier[storage_by_carrier < 0] = 0
p_by_carrier = pd.concat([p_by_carrier, storage_by_carrier], axis=1)
# Add H2_turbine links (hydrogen power plants modeled as links, not generators)
if len(network.links) > 0:
h2_turbine_links = network.links[network.links['carrier'] == 'H2_turbine']
if len(h2_turbine_links) > 0 and len(network.links_t.p1) > 0:
# p1 is the power output at bus1 (GB side), which is negative
h2_output = network.links_t.p1[h2_turbine_links.index].abs()
h2_ts = pd.DataFrame({'H2_turbine': h2_output.sum(axis=1)})
p_by_carrier = pd.concat([p_by_carrier, h2_ts], axis=1)
print(f'Added H2_turbine generation: {h2_turbine_links.p_nom.sum():.0f} MW capacity')
# Separate interconnectors from internal power lines
# Interconnectors connect to external buses (e.g., 'HVDC_External_France')
if len(network.links) > 0:
# Identify interconnectors vs internal links by checking for 'external' in bus names
is_interconnector = (
network.links.bus0.str.lower().str.contains('external', na=False) |
network.links.bus1.str.lower().str.contains('external', na=False)
)
interconnector_links = network.links[is_interconnector].index
internal_links = network.links[~is_interconnector].index
# Get interconnector imports (positive p0 = flow into GB from external)
if len(interconnector_links) > 0:
ic_flows = network.links_t.p0[interconnector_links].copy()
ic_flows[ic_flows < 0] = 0
interconnector_import = pd.DataFrame({'Interconnectors Import': ic_flows.sum(axis=1)})
p_by_carrier = pd.concat([p_by_carrier, interconnector_import], axis=1)
# Report internal links separately (not included in generation mix)
if len(internal_links) > 0:
print(f'Note: {len(internal_links)} internal power transmission links (not shown in generation mix)')
else:
print('No links in network')
print('Generation Mix (MWh):')
print(p_by_carrier.sum().sort_values(ascending=False))
print()
print(f'Total demand satisfied: {network.loads_t.p_set.sum().sum():,.0f} MWh')
Added H2_turbine generation: 27520 MW capacity
Note: 14 internal power transmission links (not shown in generation mix)
Generation Mix (MWh):
nuclear 3.614886e+06
EU_import 3.217637e+06
wind_offshore 2.398144e+06
wind_onshore 2.215468e+06
solar_pv 1.923031e+06
Interconnectors Import 1.778327e+06
H2_turbine 3.639446e+05
marine 3.578051e+05
Battery 3.237777e+05
large_hydro 1.439953e+05
Domestic Battery 1.028002e+05
Pumped Storage Hydroelectricity 1.002960e+05
landfill_gas 7.865624e+04
LAES 7.170805e+04
load_shedding 1.486904e+04
waste_to_energy 4.791646e+03
sewage_gas 4.425736e+03
geothermal 1.848000e+03
biogas 7.626765e-02
biomass 3.222728e-02
advanced_biofuel 2.525980e-03
CCGT 1.070547e-03
CHP 8.378449e-04
OCGT 5.628356e-04
gas_engine 4.833162e-04
oil 4.893091e-05
dtype: float64
Total demand satisfied: 13,231,717 MWh
Interactive Generation Mix Visualization#
[18]:
# Interactive stacked area chart with Plotly
# Define distinct colors for different technologies (ordered: baseload -> renewables -> peakers)
color_map = {
# Baseload (bottom of stack)
'nuclear': '#9467BD', # Purple - distinctive for nuclear
# Low-carbon dispatchable
'waste_to_energy': '#8C564B', # Brown
'landfill_gas': '#D2691E', # Chocolate
'biogas': '#2E8B57', # Sea green
'sewage_gas': '#556B2F', # Dark olive
'advanced_biofuel': '#6B8E23', # Olive drab
'biomass': '#228B22', # Forest green
'Other Bioenergy': '#3CB371', # Medium sea green (grouped)
# Renewables (variable)
'wind_offshore': '#1E90FF', # Dodger blue
'wind_onshore': '#00CED1', # Dark turquoise
'solar_pv': '#FFD700', # Gold
'large_hydro': '#4169E1', # Royal blue
'small_hydro': '#87CEEB', # Sky blue
'tidal_stream': '#20B2AA', # Light sea green
'shoreline_wave': '#48D1CC', # Medium turquoise
'Other Renewables': '#5F9EA0', # Cadet blue (grouped)
# Storage (mid-merit)
'battery': '#32CD32', # Lime green
'Pumped Storage Hydroelectricity': '#00FA9A', # Medium spring green
'pumped_hydro': '#00FA9A',
# Gas (flexible)
'CCGT': '#FF6347', # Tomato red
'OCGT': '#FF4500', # Orange red (peaker)
# Coal and other fossil
'coal': '#696969', # Dim grey
# Imports
'Interconnectors Import': '#DDA0DD', # Plum
'EU_import': '#DA70D6', # Orchid
# Hydrogen system
'H2': '#00FFFF', # Cyan (legacy H2 generators)
'H2_turbine': '#00CED1', # Dark turquoise (H2 to power)
'electrolysis': '#40E0D0', # Turquoise (power to H2)
'H2_gas': '#7FFFD4', # Aquamarine (H2 storage)
# Emergency
'load_shedding': '#DC143C' # Crimson
}
# Define stacking order: baseload at bottom, peakers at top
stacking_order = [
'nuclear', # Baseload (bottom)
'waste_to_energy', 'landfill_gas', 'biogas', 'sewage_gas', 'advanced_biofuel', 'biomass', 'Other Bioenergy',
'wind_offshore', 'wind_onshore', 'solar_pv', # Renewables
'large_hydro', 'small_hydro', 'tidal_stream', 'shoreline_wave', 'Other Renewables',
'battery', 'Pumped Storage Hydroelectricity', 'pumped_hydro', # Storage
'H2', 'H2_turbine', # Hydrogen power
'CCGT', # Mid-merit gas
'Interconnectors Import', 'EU_import', # Imports
'coal', # Coal
'OCGT', # Peaker
'load_shedding' # Emergency (top)
]
# Group small contributors (<1% of total) into 'Other' categories
total_gen = p_by_carrier.sum().sum()
threshold = total_gen * 0.01 # 1% threshold
# Identify small bioenergy and renewable types to group
bioenergy_types = ['biogas', 'sewage_gas', 'advanced_biofuel', 'landfill_gas']
renewable_types = ['tidal_stream', 'shoreline_wave', 'small_hydro']
p_by_carrier_grouped = p_by_carrier.copy()
# Group small bioenergy
small_bio = [c for c in bioenergy_types if c in p_by_carrier.columns and p_by_carrier[c].sum() < threshold]
if len(small_bio) > 0:
p_by_carrier_grouped['Other Bioenergy'] = p_by_carrier_grouped[small_bio].sum(axis=1)
p_by_carrier_grouped = p_by_carrier_grouped.drop(columns=small_bio)
# Group small renewables
small_ren = [c for c in renewable_types if c in p_by_carrier.columns and p_by_carrier[c].sum() < threshold]
if len(small_ren) > 0:
p_by_carrier_grouped['Other Renewables'] = p_by_carrier_grouped[small_ren].sum(axis=1)
p_by_carrier_grouped = p_by_carrier_grouped.drop(columns=small_ren)
# Select columns with generation > 0
cols_with_gen = p_by_carrier_grouped.sum()[p_by_carrier_grouped.sum() > 0].index.tolist()
# Sort by stacking order
cols_to_plot = [c for c in stacking_order if c in cols_with_gen]
# Add any remaining columns not in stacking_order
cols_to_plot += [c for c in cols_with_gen if c not in cols_to_plot]
# Create interactive stacked area chart
fig = go.Figure()
for col in cols_to_plot: # Stack in order
fig.add_trace(go.Scatter(
x=p_by_carrier_grouped.index,
y=p_by_carrier_grouped[col] / 1e3, # Convert to GW
name=col.replace('_', ' ').title(),
mode='lines',
stackgroup='one',
fillcolor=color_map.get(col, '#CCCCCC'),
line=dict(width=0.5, color=color_map.get(col, '#CCCCCC')),
hovertemplate='<b>%{fullData.name}</b><br>' +
'Time: %{x}<br>' +
'Power: %{y:.2f} GW<br>' +
'<extra></extra>'
))
fig.update_layout(
title=dict(text='Generation Mix - EE50_clustered', x=0.5, xanchor='center'),
xaxis_title='Time',
yaxis_title='Generation (GW)',
hovermode='x unified',
height=600,
legend=dict(
orientation='v',
yanchor='top',
y=1,
xanchor='left',
x=1.02
),
template='plotly_white'
)
fig.show()
print('Stacking order: Baseload (bottom) -> Renewables -> Storage -> Gas -> Peakers (top)')
print('Tip: Click legend items to show/hide technologies. Double-click to isolate one.')
Stacking order: Baseload (bottom) -> Renewables -> Storage -> Gas -> Peakers (top)
Tip: Click legend items to show/hide technologies. Double-click to isolate one.
Interactive Network Topology#
Explore the network topology with interactive maps showing bus locations, transmission lines/links, and generators overlaid on a GB map.
Understanding Network Topology#
The maps below show the physical layout of Great Britain’s electricity infrastructure:
Buses: Connection points (substations, generation sites) shown as orange markers
Lines/Links: Transmission corridors carrying power between buses
Geography: Real GB coordinates using British National Grid projection
The network model captures how electricity flows are constrained by transmission capacity, creating locational price differences and potential congestion.
[19]:
# Interactive network map using PyPSA's explore() function
# Select a representative snapshot for visualization
snapshot_idx = len(network.snapshots) // 2
snapshot = network.snapshots[snapshot_idx]
print(f'Visualizing network at snapshot: {snapshot}')
# Detect network type based on bus count
is_zonal = len(network.buses) < 50 and len(network.lines) == 0 and len(network.links) > 0
is_reduced = len(network.buses) < 100 and len(network.buses) >= 30
network_type = 'Zonal' if is_zonal else ('Reduced' if is_reduced else 'ETYS')
print(f'Network type detected: {network_type} ({len(network.buses)} buses)')
# Network topology visualization (electricity only, no hydrogen)
try:
import copy
import folium
if 'x' in getattr(network.buses, 'columns', []) and 'y' in getattr(network.buses, 'columns', []):
# Filter out hydrogen buses and links
h2_carriers = ['electrolysis', 'H2_turbine', 'H2_power', 'H2', 'H2_gas', 'H2 pipeline']
h2_buses = set()
if 'carrier' in network.buses.columns:
h2_buses = set(network.buses[network.buses.carrier.isin(h2_carriers)].index)
# Create filtered network copy
network_plot = copy.deepcopy(network)
# Remove hydrogen buses if any found
if len(h2_buses) > 0:
for bus_id in h2_buses:
if bus_id in network_plot.buses.index:
network_plot.remove('Bus', bus_id)
# Also filter hydrogen links
if len(network_plot.links) > 0 and 'carrier' in network_plot.links.columns:
h2_links = network_plot.links[network_plot.links.carrier.isin(h2_carriers)].index
for link_id in h2_links:
if link_id in network_plot.links.index:
network_plot.remove('Link', link_id)
# Convert OSGB -> WGS84 lon/lat for mapping
t = Transformer.from_crs('EPSG:27700', 'EPSG:4326', always_xy=True)
lon, lat = t.transform(network_plot.buses['x'].to_numpy(), network_plot.buses['y'].to_numpy())
network_plot.buses['x'] = lon
network_plot.buses['y'] = lat
print(f'Electricity network: {len(network_plot.buses)} buses, {len(network_plot.lines)} lines, {len(network_plot.links)} links')
# Create interactive map
m = network_plot.plot.explore(map_style='light', tooltip=True)
display(m)
print('✓ Network topology displayed (electricity only)')
else:
print('⚠️ No bus coordinates (x, y) found. Network visualization requires coordinate data.')
except Exception as e:
print(f'⚠️ Network topology map unavailable: {e}')
import traceback
traceback.print_exc()
Visualizing network at snapshot: 2050-02-04 12:00:00
Network type detected: ETYS (110 buses)
Electricity network: 109 buses, 202 lines, 27 links
✓ Network topology displayed (electricity only)
[20]:
# Interactive network map showing transmission loading using PyPSA explore
try:
import copy
import folium
from folium import plugins
# Detect transmission element type (electricity only)
has_lines = len(network.lines) > 0 and 's_nom' in network.lines.columns
has_internal_links = False
if has_lines:
# ETYS/Reduced network - use lines
loading = (network.lines_t.p0.loc[snapshot].abs() / network.lines.s_nom).fillna(0)
component_type = 'lines'
elif len(network.links) > 0:
# Zonal network - use internal electricity links only
h2_carriers = ['electrolysis', 'H2_turbine', 'H2_power', 'H2', 'H2_gas', 'H2 pipeline']
is_external = (
network.links.bus0.str.lower().str.contains('external', na=False) |
network.links.bus1.str.lower().str.contains('external', na=False)
)
is_h2_link = network.links.carrier.isin(h2_carriers)
internal_links = network.links[~is_external & ~is_h2_link]
if len(internal_links) > 0 and snapshot in network.links_t.p0.index:
link_flows = network.links_t.p0.loc[snapshot, internal_links.index].abs()
loading = (link_flows / internal_links.p_nom).fillna(0)
component_type = 'links'
has_internal_links = True
else:
loading = pd.Series(dtype=float)
if len(loading) > 0:
# Convert bus coordinates
t = Transformer.from_crs('EPSG:27700', 'EPSG:4326', always_xy=True)
bus_lon, bus_lat = t.transform(network.buses['x'].to_numpy(), network.buses['y'].to_numpy())
bus_coords = pd.DataFrame({'lon': bus_lon, 'lat': bus_lat}, index=network.buses.index)
# Create folium map
center_lat = bus_coords['lat'].mean()
center_lon = bus_coords['lon'].mean()
m = folium.Map(location=[center_lat, center_lon], zoom_start=5, tiles='CartoDB positron')
# Color function
def get_color_by_loading(load_val):
if load_val >= 0.95: return '#DC143C' # Crimson (>95%)
elif load_val >= 0.80: return '#FF8C00' # Orange (80-95%)
elif load_val >= 0.60: return '#FFD700' # Gold (60-80%)
else: return '#228B22' # Green (<60%)
# Draw transmission elements with loading colors
if has_lines:
for line_id, load_val in loading.items():
line = network.lines.loc[line_id]
bus0_coords = [bus_coords.loc[line.bus0, 'lat'], bus_coords.loc[line.bus0, 'lon']]
bus1_coords = [bus_coords.loc[line.bus1, 'lat'], bus_coords.loc[line.bus1, 'lon']]
color = get_color_by_loading(load_val)
weight = 0.5 + (load_val * 4) # Width scales with loading
folium.PolyLine(
[bus0_coords, bus1_coords],
color=color,
weight=weight,
opacity=0.7,
tooltip=f'{line_id}: {load_val*100:.1f}% loaded'
).add_to(m)
elif has_internal_links:
for link_id, load_val in loading.items():
link = network.links.loc[link_id]
bus0_coords = [bus_coords.loc[link.bus0, 'lat'], bus_coords.loc[link.bus0, 'lon']]
bus1_coords = [bus_coords.loc[link.bus1, 'lat'], bus_coords.loc[link.bus1, 'lon']]
color = get_color_by_loading(load_val)
weight = 0.5 + (load_val * 4)
folium.PolyLine(
[bus0_coords, bus1_coords],
color=color,
weight=weight,
opacity=0.7,
tooltip=f'{link_id}: {load_val*100:.1f}% loaded'
).add_to(m)
# Add small bus markers
for bus_id, row in bus_coords.iterrows():
folium.CircleMarker(
location=[row['lat'], row['lon']],
radius=2,
color='orange',
fill=True,
fillOpacity=0.6,
tooltip=bus_id
).add_to(m)
display(m)
# Print statistics
print(f'\nTransmission Loading ({component_type}) at {snapshot}:')
print(f' Max: {loading.max()*100:.1f}% | Mean: {loading.mean()*100:.1f}%')
print(f' >95% loaded: {(loading >= 0.95).sum()} | >80% loaded: {(loading >= 0.80).sum()}')
print('Color legend: Green (<60%) → Gold (60-80%) → Orange (80-95%) → Red (>95%)')
else:
print(f'⚠️ No transmission loading data available at {snapshot}')
except Exception as e:
print(f'⚠️ Transmission loading map unavailable: {e}')
import traceback
traceback.print_exc()
Make this Notebook Trusted to load map: File -> Trust Notebook
Transmission Loading (lines) at 2050-02-04 12:00:00:
Max: 100.0% | Mean: 23.7%
>95% loaded: 13 | >80% loaded: 17
Color legend: Green (<60%) → Gold (60-80%) → Orange (80-95%) → Red (>95%)
[21]:
# Interactive network map showing generator dispatch (electricity only)
try:
import folium
# Get total generation per bus at snapshot (exclude hydrogen generators)
h2_carriers = ['electrolysis', 'H2_turbine', 'H2_power', 'H2']
elec_gens = network.generators[~network.generators.carrier.isin(h2_carriers)]
gen_per_bus = network.generators_t.p.loc[snapshot, elec_gens.index].groupby(elec_gens.bus).sum()
gen_per_bus = gen_per_bus[gen_per_bus > 0] # Only buses with generation
if len(gen_per_bus) > 0 and 'x' in getattr(network.buses, 'columns', []):
# Convert bus coordinates
t = Transformer.from_crs('EPSG:27700', 'EPSG:4326', always_xy=True)
bus_lon, bus_lat = t.transform(network.buses['x'].to_numpy(), network.buses['y'].to_numpy())
bus_coords = pd.DataFrame({'lon': bus_lon, 'lat': bus_lat}, index=network.buses.index)
# Create folium map
center_lat = bus_coords['lat'].mean()
center_lon = bus_coords['lon'].mean()
m = folium.Map(location=[center_lat, center_lon], zoom_start=5, tiles='CartoDB positron')
# Draw light gray transmission lines in background
if len(network.lines) > 0:
for line_id, line in network.lines.iterrows():
bus0_coords = [bus_coords.loc[line.bus0, 'lat'], bus_coords.loc[line.bus0, 'lon']]
bus1_coords = [bus_coords.loc[line.bus1, 'lat'], bus_coords.loc[line.bus1, 'lon']]
folium.PolyLine([bus0_coords, bus1_coords], color='lightgray', weight=0.5, opacity=0.3).add_to(m)
# Color and size encoding for generation
max_gen = gen_per_bus.max()
def color_by_generation(gen_val):
norm = gen_val / max_gen
if norm < 0.25: return '#FFE4E1' # Misty rose
elif norm < 0.50: return '#FF9999' # Light red
elif norm < 0.75: return '#FF6666' # Red
else: return '#CC0000' # Dark red
# Draw bus markers sized and colored by generation
for bus_id, gen_mw in gen_per_bus.items():
coords = bus_coords.loc[bus_id]
norm_gen = gen_mw / max_gen
radius = 3 + (norm_gen * 12) # Size 3-15 pixels
color = color_by_generation(gen_mw)
folium.CircleMarker(
location=[coords['lat'], coords['lon']],
radius=radius,
color=color,
fill=True,
fillColor=color,
fillOpacity=0.7,
tooltip=f'{bus_id}: {gen_mw:.1f} MW'
).add_to(m)
display(m)
# Statistics
print(f'\nGenerator Dispatch at {snapshot}:')
print(f' Total generation: {gen_per_bus.sum():,.0f} MW')
print(f' Active buses: {len(gen_per_bus[gen_per_bus > 0])}')
print(f' Max at single bus: {gen_per_bus.max():,.0f} MW')
print('Color legend: Light gray (no gen) → Light red → Dark red (peak generation)')
else:
print('⚠️ No generator dispatch data available at this snapshot')
except Exception as e:
print(f'⚠️ Generator dispatch map unavailable: {e}')
import traceback
traceback.print_exc()
Make this Notebook Trusted to load map: File -> Trust Notebook
Generator Dispatch at 2050-02-04 12:00:00:
Total generation: 116,924 MW
Active buses: 108
Max at single bus: 6,990 MW
Color legend: Light gray (no gen) → Light red → Dark red (peak generation)
[22]:
# Interactive installed capacity by carrier
generators_p_nom = network.generators.p_nom.groupby(
network.generators.carrier).sum().sort_values(ascending=True)
# Remove small contributors
generators_p_nom = generators_p_nom[generators_p_nom > 50]
# Create interactive bar chart
fig = go.Figure()
fig.add_trace(go.Bar(
y=generators_p_nom.index,
x=generators_p_nom.values,
orientation='h',
marker=dict(
color=[color_map.get(c, '#CCCCCC') for c in generators_p_nom.index],
),
text=[f'{val:,.0f} MW' for val in generators_p_nom.values],
textposition='auto',
hovertemplate='<b>%{y}</b><br>Capacity: %{x:,.0f} MW<extra></extra>'
))
fig.update_layout(
title=dict(text='Generator Capacity by Technology - EE50_clustered', x=0.5, xanchor='center'),
xaxis_title='Installed Capacity (MW)',
yaxis_title='',
height=500,
template='plotly_white',
showlegend=False
)
fig.show()
Storage Analysis#
[23]:
# Interactive storage dispatch and state of charge visualization
if len(network.storage_units) > 0:
p_storage = network.storage_units_t.p.sum(axis=1)
state_of_charge = network.storage_units_t.state_of_charge.sum(axis=1)
# Create subplot with two y-axes
fig = make_subplots(
rows=2, cols=1,
subplot_titles=('Storage Dispatch', 'State of Charge'),
vertical_spacing=0.12
)
# Storage dispatch
fig.add_trace(
go.Scatter(
x=p_storage.index,
y=p_storage.values,
name='Storage Dispatch',
line=dict(color='#32CD32', width=1.5),
fill='tozeroy',
fillcolor='rgba(50, 205, 50, 0.3)',
hovertemplate='Time: %{x}<br>Power: %{y:.0f} MW<extra></extra>'
),
row=1, col=1
)
# State of charge
fig.add_trace(
go.Scatter(
x=state_of_charge.index,
y=state_of_charge.values,
name='State of Charge',
line=dict(color='#00CED1', width=1.5),
fill='tozeroy',
fillcolor='rgba(0, 206, 209, 0.3)',
hovertemplate='Time: %{x}<br>Energy: %{y:.0f} MWh<extra></extra>'
),
row=2, col=1
)
fig.update_xaxes(title_text='Time', row=2, col=1)
fig.update_yaxes(title_text='Power (MW)', row=1, col=1)
fig.update_yaxes(title_text='Energy (MWh)', row=2, col=1)
fig.update_layout(
title=dict(text='Storage Analysis - EE50_clustered', x=0.5, xanchor='center'),
height=700,
hovermode='x unified',
template='plotly_white',
showlegend=True
)
fig.show()
# Storage statistics
print(f'\nStorage Statistics:')
print(f' Total storage capacity: {network.storage_units.p_nom.sum():,.0f} MW')
print(f' Total energy capacity: {network.storage_units.max_hours.sum() * network.storage_units.p_nom.sum():,.0f} MWh')
print(f' Peak discharge: {p_storage.max():,.0f} MW')
print(f' Peak charge: {p_storage.min():,.0f} MW')
print(f' Max state of charge: {state_of_charge.max():,.0f} MWh')
else:
print('No storage units in network')
Storage Statistics:
Total storage capacity: 46,099 MW
Total energy capacity: 121,609,119 MWh
Peak discharge: 22,143 MW
Peak charge: -36,731 MW
Max state of charge: 158,096 MWh
Hydrogen System Analysis#
The hydrogen system models the Power-to-Gas-to-Power pathway:
Electrolysis: Converts electricity to hydrogen (stored in central H2 bus)
H2 Storage: Central hydrogen storage (copper-plate model)
H2 Turbines: Converts hydrogen back to electricity when needed
About the Hydrogen System#
The hydrogen subsystem models the Power-to-Gas-to-Power pathway:
Electrolysis: Converts surplus renewable electricity to hydrogen (dashed turquoise links)
H2 Storage: Stores hydrogen for later use (purple markers)
H2 Turbines: Converts hydrogen back to electricity during high demand (red markers)
This provides seasonal energy storage, shifting summer renewable surplus to winter demand peaks.
[24]:
# Hydrogen System Analysis
# Check if hydrogen system exists
has_electrolysis = len(network.links) > 0 and 'electrolysis' in network.links.carrier.values
has_h2_turbines = len(network.links) > 0 and 'H2_turbine' in network.links.carrier.values
has_h2_stores = len(network.stores) > 0 and 'H2_gas' in network.stores.carrier.values
if has_electrolysis or has_h2_turbines:
print('=' * 80)
print('HYDROGEN SYSTEM SUMMARY')
print('=' * 80)
# Electrolysis
if has_electrolysis:
electrolysis = network.links[network.links.carrier == 'electrolysis']
elec_capacity = electrolysis.p_nom.sum()
elec_consumed = network.links_t.p0[electrolysis.index].sum().sum() / 1000 # GWh
elec_efficiency = electrolysis.efficiency.mean()
h2_produced = elec_consumed * elec_efficiency # GWh H2
print(f'\nELECTROLYSIS:')
print(f' Capacity: {elec_capacity:,.0f} MW ({len(electrolysis)} units)')
print(f' Efficiency: {elec_efficiency*100:.0f}%')
print(f' Electricity consumed: {elec_consumed:,.1f} GWh')
print(f' Hydrogen produced: {h2_produced:,.1f} GWh_H2')
# H2 Turbines
if has_h2_turbines:
h2_turbines = network.links[network.links.carrier == 'H2_turbine']
turbine_capacity = h2_turbines.p_nom.sum()
h2_consumed = network.links_t.p0[h2_turbines.index].sum().sum() / 1000 # GWh H2
turbine_efficiency = h2_turbines.efficiency.mean()
elec_generated = h2_consumed * turbine_efficiency # GWh electricity
print(f'\nH2 TURBINES:')
print(f' Capacity: {turbine_capacity:,.0f} MW ({len(h2_turbines)} units)')
print(f' Efficiency: {turbine_efficiency*100:.0f}%')
print(f' Hydrogen consumed: {h2_consumed:,.1f} GWh_H2')
print(f' Electricity generated: {elec_generated:,.1f} GWh')
# H2 Storage
if has_h2_stores:
h2_stores = network.stores[network.stores.carrier == 'H2_gas']
h2_storage_capacity = h2_stores.e_nom.sum()
print(f'\nH2 STORAGE:')
print(f' Energy capacity: {h2_storage_capacity:,.0f} MWh ({h2_storage_capacity/1000:.0f} GWh)')
if h2_stores.index[0] in network.stores_t.e.columns:
soc = network.stores_t.e[h2_stores.index[0]]
print(f' State of charge range: {soc.min():,.0f} - {soc.max():,.0f} MWh')
print(f' Final state of charge: {soc.iloc[-1]:,.0f} MWh')
# Round-trip efficiency
if has_electrolysis and has_h2_turbines:
roundtrip_eff = elec_efficiency * turbine_efficiency
print(f'\nSYSTEM METRICS:')
print(f' Round-trip efficiency: {roundtrip_eff*100:.1f}%')
if elec_consumed > 0:
actual_roundtrip = elec_generated / elec_consumed
print(f' Actual energy ratio (elec out / elec in): {actual_roundtrip*100:.1f}%')
else:
print('No hydrogen system components found in this network.')
print('(Hydrogen components are added for future scenarios with H2 power generation)')
================================================================================
HYDROGEN SYSTEM SUMMARY
================================================================================
ELECTROLYSIS:
Capacity: 28,946 MW (24 units)
Efficiency: 70%
Electricity consumed: 1,046.4 GWh
Hydrogen produced: 732.5 GWh_H2
H2 TURBINES:
Capacity: 27,520 MW (192 units)
Efficiency: 50%
Hydrogen consumed: 727.9 GWh_H2
Electricity generated: 363.9 GWh
H2 STORAGE:
Energy capacity: 4,623,382 MWh (4623 GWh)
State of charge range: 0 - 712,488 MWh
Final state of charge: 56,491 MWh
SYSTEM METRICS:
Round-trip efficiency: 35.0%
Actual energy ratio (elec out / elec in): 34.8%
[25]:
# Hydrogen System Dispatch Visualization
has_electrolysis = len(network.links) > 0 and 'electrolysis' in network.links.carrier.values
has_h2_turbines = len(network.links) > 0 and 'H2_turbine' in network.links.carrier.values
has_h2_stores = len(network.stores) > 0 and 'H2_gas' in network.stores.carrier.values
if has_electrolysis and has_h2_turbines:
electrolysis = network.links[network.links.carrier == 'electrolysis']
h2_turbines = network.links[network.links.carrier == 'H2_turbine']
# Electrolysis power consumption (positive = consuming electricity)
elec_power = network.links_t.p0[electrolysis.index].sum(axis=1) / 1000 # GW
# H2 turbine power generation (p1 is output, negative means generation)
turbine_power = -network.links_t.p1[h2_turbines.index].sum(axis=1) / 1000 # GW
# Net hydrogen flow (positive = production, negative = consumption)
h2_production = elec_power * electrolysis.efficiency.mean() # GW H2 produced
h2_consumption = network.links_t.p0[h2_turbines.index].sum(axis=1) / 1000 # GW H2 consumed
# Create visualization
fig = make_subplots(
rows=3, cols=1,
subplot_titles=(
'Electrolysis (Electricity Consumed)',
'H2 Turbines (Electricity Generated)',
'H2 Storage State of Charge'
),
vertical_spacing=0.1
)
# Electrolysis power
fig.add_trace(
go.Scatter(
x=elec_power.index,
y=elec_power.values,
name='Electrolysis',
line=dict(color='#00CED1', width=1.5),
fill='tozeroy',
fillcolor='rgba(0, 206, 209, 0.3)',
hovertemplate='Time: %{x}<br>Power: %{y:.2f} GW<extra></extra>'
),
row=1, col=1
)
# H2 turbine power
fig.add_trace(
go.Scatter(
x=turbine_power.index,
y=turbine_power.values,
name='H2 Turbines',
line=dict(color='#FF6347', width=1.5),
fill='tozeroy',
fillcolor='rgba(255, 99, 71, 0.3)',
hovertemplate='Time: %{x}<br>Power: %{y:.2f} GW<extra></extra>'
),
row=2, col=1
)
# H2 storage state of charge
if has_h2_stores:
h2_stores = network.stores[network.stores.carrier == 'H2_gas']
if h2_stores.index[0] in network.stores_t.e.columns:
h2_soc = network.stores_t.e[h2_stores.index[0]] / 1000 # GWh
fig.add_trace(
go.Scatter(
x=h2_soc.index,
y=h2_soc.values,
name='H2 Storage',
line=dict(color='#9467BD', width=1.5),
fill='tozeroy',
fillcolor='rgba(148, 103, 189, 0.3)',
hovertemplate='Time: %{x}<br>Energy: %{y:.1f} GWh<extra></extra>'
),
row=3, col=1
)
fig.update_yaxes(title_text='Power (GW)', row=1, col=1)
fig.update_yaxes(title_text='Power (GW)', row=2, col=1)
fig.update_yaxes(title_text='Energy (GWh)', row=3, col=1)
fig.update_xaxes(title_text='Time', row=3, col=1)
fig.update_layout(
title=dict(text='Hydrogen System Dispatch - EE50_clustered', x=0.5, xanchor='center'),
height=900,
hovermode='x unified',
template='plotly_white',
showlegend=True
)
fig.show()
# Correlation analysis
print('\nDispatch Patterns:')
print(f' Peak electrolysis: {elec_power.max():.2f} GW at {elec_power.idxmax()}')
print(f' Peak H2 turbine: {turbine_power.max():.2f} GW at {turbine_power.idxmax()}')
print(f' Capacity factor (electrolysis): {elec_power.mean() / (electrolysis.p_nom.sum()/1000) * 100:.1f}%')
print(f' Capacity factor (H2 turbines): {turbine_power.mean() / (h2_turbines.p_nom.sum()/1000) * 100:.1f}%')
else:
print('Hydrogen dispatch visualization not available (no hydrogen system components)')
Dispatch Patterns:
Peak electrolysis: 28.95 GW at 2050-02-01 13:00:00
Peak H2 turbine: 13.29 GW at 2050-02-02 22:00:00
Capacity factor (electrolysis): 21.5%
Capacity factor (H2 turbines): 7.9%
[26]:
# Interactive spatial map of hydrogen infrastructure
try:
import folium
has_h2_system = (
(len(network.links) > 0 and any(network.links.carrier.isin(['electrolysis', 'H2_turbine', 'H2 pipeline']))) or
(len(network.stores) > 0 and 'H2_gas' in network.stores.carrier.values)
)
if has_h2_system and 'x' in network.buses.columns:
# Convert all bus coordinates
t = Transformer.from_crs('EPSG:27700', 'EPSG:4326', always_xy=True)
bus_lon, bus_lat = t.transform(network.buses['x'].to_numpy(), network.buses['y'].to_numpy())
bus_coords = pd.DataFrame({'lon': bus_lon, 'lat': bus_lat}, index=network.buses.index)
# Create map centered on GB
center_lat = bus_coords['lat'].mean()
center_lon = bus_coords['lon'].mean()
m = folium.Map(location=[center_lat, center_lon], zoom_start=5, tiles='CartoDB positron')
# 1. Draw H2 pipelines (H2 pipeline links)
if len(network.links) > 0:
h2_pipelines = network.links[network.links.carrier == 'H2 pipeline']
if len(h2_pipelines) > 0:
for link_id, link in h2_pipelines.iterrows():
if link.bus0 in bus_coords.index and link.bus1 in bus_coords.index:
bus0_coords = [bus_coords.loc[link.bus0, 'lat'], bus_coords.loc[link.bus0, 'lon']]
bus1_coords = [bus_coords.loc[link.bus1, 'lat'], bus_coords.loc[link.bus1, 'lon']]
capacity_mw = link.p_nom
folium.PolyLine(
[bus0_coords, bus1_coords],
color='#00CED1', # Cyan for H2 pipelines
weight=3,
opacity=0.7,
tooltip=f'H2 Pipeline: {link_id}<br>Capacity: {capacity_mw:.0f} MW'
).add_to(m)
print(f'Added {len(h2_pipelines)} H2 pipelines to map')
# 2. Draw electrolysis links (electricity -> H2) and facilities
if len(network.links) > 0:
electrolysis = network.links[network.links.carrier == 'electrolysis']
if len(electrolysis) > 0:
# Draw electrolysis links connecting electricity bus to H2 storage
for link_id, link in electrolysis.iterrows():
# bus0 is electricity input, bus1 is H2 output/storage
if link.bus0 in bus_coords.index and link.bus1 in bus_coords.index:
bus0_coords = [bus_coords.loc[link.bus0, 'lat'], bus_coords.loc[link.bus0, 'lon']]
bus1_coords = [bus_coords.loc[link.bus1, 'lat'], bus_coords.loc[link.bus1, 'lon']]
folium.PolyLine(
[bus0_coords, bus1_coords],
color='#40E0D0', # Turquoise for electrolysis links
weight=2,
opacity=0.6,
dash_array='5, 5', # Dashed line
tooltip=f'Electrolysis Link: {link_id}<br>Capacity: {link.p_nom:.0f} MW'
).add_to(m)
# Draw electrolysis facility markers at electricity input bus
for link_id, link in electrolysis.iterrows():
if link.bus0 in bus_coords.index:
coords = bus_coords.loc[link.bus0]
folium.CircleMarker(
location=[coords['lat'], coords['lon']],
radius=6 + (link.p_nom / 500), # Size by capacity
color='#40E0D0', # Turquoise
fill=True,
fillColor='#40E0D0',
fillOpacity=0.7,
tooltip=f'Electrolysis: {link_id}<br>Capacity: {link.p_nom:.0f} MW<br>Electricity bus: {link.bus0}<br>H2 bus: {link.bus1}'
).add_to(m)
print(f'Added {len(electrolysis)} electrolysis units and links to map')
# 3. Draw H2 turbine locations (H2 -> electricity)
if len(network.links) > 0:
h2_turbines = network.links[network.links.carrier == 'H2_turbine']
if len(h2_turbines) > 0:
for link_id, link in h2_turbines.iterrows():
# bus1 is electricity output
if link.bus1 in bus_coords.index:
coords = bus_coords.loc[link.bus1]
folium.CircleMarker(
location=[coords['lat'], coords['lon']],
radius=6 + (link.p_nom / 500), # Size by capacity
color='#FF6347', # Tomato red
fill=True,
fillColor='#FF6347',
fillOpacity=0.7,
tooltip=f'H2 Turbine: {link_id}<br>Capacity: {link.p_nom:.0f} MW<br>Bus: {link.bus1}'
).add_to(m)
print(f'Added {len(h2_turbines)} H2 turbine units to map')
# 4. Draw H2 storage locations
if len(network.stores) > 0:
h2_stores = network.stores[network.stores.carrier == 'H2_gas']
if len(h2_stores) > 0:
for store_id, store in h2_stores.iterrows():
if store.bus in bus_coords.index:
coords = bus_coords.loc[store.bus]
energy_gwh = store.e_nom / 1000
folium.CircleMarker(
location=[coords['lat'], coords['lon']],
radius=8,
color='#9467BD', # Purple
fill=True,
fillColor='#9467BD',
fillOpacity=0.8,
tooltip=f'H2 Storage: {store_id}<br>Capacity: {energy_gwh:.1f} GWh<br>Bus: {store.bus}'
).add_to(m)
print(f'Added {len(h2_stores)} H2 storage facilities to map')
# Add legend
legend_html = '''<div style=\"position: fixed; bottom: 50px; right: 50px; width: 240px;
background-color: white; border:2px solid grey; z-index:9999; font-size:14px;
padding: 10px\">\n
<p style=\"margin:0; font-weight:bold;\">Hydrogen Infrastructure</p>\n
<p style=\"margin:5px 0;\"><span style=\"color:#00CED1; font-size:20px;\">━━</span> H2 Pipeline</p>\n
<p style=\"margin:5px 0;\"><span style=\"color:#40E0D0; font-size:12px;\">- - -</span> Electrolysis Link</p>\n
<p style=\"margin:5px 0;\"><span style=\"color:#40E0D0; font-size:20px;\">●</span> Electrolysis</p>\n
<p style=\"margin:5px 0;\"><span style=\"color:#FF6347; font-size:20px;\">●</span> H2 Turbine</p>\n
<p style=\"margin:5px 0;\"><span style=\"color:#9467BD; font-size:20px;\">●</span> H2 Storage</p>\n
</div>'''
m.get_root().html.add_child(folium.Element(legend_html))
display(m)
print('\\n✓ Hydrogen network map displayed')
print('Marker size indicates capacity (larger = higher capacity)')
else:
print('No hydrogen infrastructure found or missing bus coordinates')
except Exception as e:
print(f'⚠️ Hydrogen network map unavailable: {e}')
import traceback
traceback.print_exc()
Added 24 electrolysis units and links to map
Added 192 H2 turbine units to map
Added 1 H2 storage facilities to map
Make this Notebook Trusted to load map: File -> Trust Notebook
\n✓ Hydrogen network map displayed
Marker size indicates capacity (larger = higher capacity)
Transmission Analysis#
[27]:
# Interactive transmission loading analysis with congestion hotspots
# Handles both lines (ETYS/Reduced) and links (Zonal)
snapshot_idx = len(network.snapshots) // 2
snapshot = network.snapshots[snapshot_idx]
# Determine which transmission elements to analyze
has_lines = len(network.lines) > 0 and 's_nom' in network.lines.columns
has_internal_links = len(network.links) > 0
if has_lines:
# ETYS/Reduced network - use lines
loading = (network.lines_t.p0.loc[snapshot].abs() / network.lines.s_nom).fillna(0)
components = network.lines
element_type = 'Line'
elif has_internal_links:
# Zonal network - use internal links only (exclude external and hydrogen links)
h2_carriers = ['electrolysis', 'H2_turbine', 'H2_power', 'H2', 'H2_gas']
is_external = (
network.links.bus0.str.lower().str.contains('external', na=False) |
network.links.bus1.str.lower().str.contains('external', na=False)
)
is_h2_link = network.links.carrier.isin(h2_carriers)
internal_links = network.links[~is_external & ~is_h2_link]
if len(internal_links) > 0 and snapshot in network.links_t.p0.index:
link_flows = network.links_t.p0.loc[snapshot, internal_links.index].abs()
loading = (link_flows / internal_links.p_nom).fillna(0)
components = internal_links
element_type = 'Link'
else:
loading = pd.Series(dtype=float)
element_type = 'None'
else:
loading = pd.Series(dtype=float)
element_type = 'None'
if len(loading) > 0:
loading_sorted = loading.sort_values(ascending=False)
top_n = min(20, len(loading_sorted))
top_elements = loading_sorted.head(top_n)
# Create labels with bus names
element_labels = []
for elem_id in top_elements.index:
bus0 = components.loc[elem_id, 'bus0']
bus1 = components.loc[elem_id, 'bus1']
element_labels.append(f'{bus0} - {bus1}')
# Color code by congestion level
colors_list = [
'#DC143C' if x >= 0.95 else '#FF8C00' if x >= 0.8 else '#FFD700' if x >= 0.6 else '#228B22'
for x in top_elements.values
]
fig = go.Figure()
fig.add_trace(go.Bar(
y=element_labels,
x=top_elements.values * 100,
orientation='h',
marker=dict(color=colors_list),
text=[f'{val*100:.1f}%' for val in top_elements.values],
textposition='auto',
hovertemplate='<b>%{y}</b><br>Loading: %{x:.1f}%<extra></extra>'
))
# Add reference lines
fig.add_vline(x=95, line_dash='dash', line_color='red', annotation_text='95% (Congestion)')
fig.add_vline(x=100, line_dash='solid', line_color='darkred', annotation_text='100% (Limit)')
fig.update_layout(
title=dict(text=f'Top {top_n} Most Loaded {element_type}s - {scenario_id}', x=0.5, xanchor='center'),
xaxis_title=f'{element_type} Loading (%)',
yaxis_title='',
height=max(400, top_n * 25),
template='plotly_white',
showlegend=False
)
fig.show()
# Statistics
print(f'\n{element_type} Loading Statistics (snapshot: {snapshot})')
print(f' Total {element_type.lower()}s: {len(loading)}')
print(f' Max loading: {loading.max()*100:.1f}%')
print(f' Mean loading: {loading.mean()*100:.1f}%')
print(f' {element_type}s >95% loaded: {(loading >= 0.95).sum()}')
print(f' {element_type}s >80% loaded: {(loading >= 0.80).sum()}')
else:
print('No transmission elements found for loading analysis')
Line Loading Statistics (snapshot: 2050-02-04 12:00:00)
Total lines: 202
Max loading: 100.0%
Mean loading: 23.7%
Lines >95% loaded: 13
Lines >80% loaded: 17
Renewable Curtailment Analysis#
[28]:
# Interactive renewable curtailment analysis
renewable_carriers = ['wind_onshore', 'wind_offshore', 'solar_pv']
# Create subplots for each renewable carrier
fig = make_subplots(
rows=len(renewable_carriers), cols=1,
subplot_titles=[carrier.replace('_', ' ').title() for carrier in renewable_carriers],
vertical_spacing=0.08
)
curtailment_stats = {}
for idx, carrier in enumerate(renewable_carriers, 1):
# Get generators of this carrier
gens = network.generators[network.generators.carrier == carrier]
if len(gens) == 0:
continue
# Calculate available and dispatched
p_nom_sum = gens.p_nom.sum()
# Check if any of these generators have p_max_pu time series data
# p_max_pu columns are indexed by generator name, not carrier
gens_with_pmax = [g for g in gens.index if g in network.generators_t.p_max_pu.columns]
if len(gens_with_pmax) > 0:
# Calculate available capacity from p_max_pu (weather-dependent availability)
available = network.generators_t.p_max_pu[gens_with_pmax].multiply(
gens.loc[gens_with_pmax, 'p_nom'], axis=1).sum(axis=1)
# Add any generators without p_max_pu at their full capacity
gens_without_pmax = [g for g in gens.index if g not in network.generators_t.p_max_pu.columns]
if len(gens_without_pmax) > 0:
available = available + gens.loc[gens_without_pmax, 'p_nom'].sum()
else:
# No time-varying availability - use installed capacity
available = pd.Series(p_nom_sum, index=network.snapshots)
dispatched = network.generators_t.p[gens.index].sum(axis=1)
curtailed = available - dispatched
# Add available capacity trace
fig.add_trace(
go.Scatter(
x=available.index,
y=available / 1000,
name='Available',
line=dict(color='orange', width=1.5),
legendgroup=carrier,
showlegend=(idx == 1),
hovertemplate='Available: %{y:.2f} GW<extra></extra>'
),
row=idx, col=1
)
# Add available capacity trace (orange line showing weather-dependent potential)
fig.add_trace(
go.Scatter(
x=available.index,
y=available / 1000,
name='Available',
line=dict(color='#FF8C00', width=1.5),
mode='lines',
legendgroup=carrier,
showlegend=(idx == 1),
hovertemplate='Available: %{y:.2f} GW<extra></extra>'
),
row=idx, col=1
)
# Add dispatched trace (filled area from zero)
fig.add_trace(
go.Scatter(
x=dispatched.index,
y=dispatched / 1000,
name='Dispatched',
line=dict(color='#228B22', width=1.5),
fill='tozeroy',
fillcolor='rgba(34, 139, 34, 0.6)',
legendgroup=carrier,
showlegend=(idx == 1),
hovertemplate='Dispatched: %{y:.2f} GW<extra></extra>'
),
row=idx, col=1
)
# Add curtailed as separate trace showing actual curtailment
curtailed_positive = curtailed.clip(lower=0) # Only show positive curtailment
fig.add_trace(
go.Scatter(
x=curtailed_positive.index,
y=curtailed_positive / 1000,
name='Curtailed',
line=dict(color='#DC143C', width=1.5, dash='dot'),
mode='lines',
legendgroup=carrier,
showlegend=(idx == 1),
hovertemplate='Curtailed: %{y:.2f} GW<extra></extra>'
),
row=idx, col=1
)
# Update y-axis label
fig.update_yaxes(title_text='Power (GW)', row=idx, col=1)
# Calculate statistics
curtailment_pct = (curtailed.sum() / available.sum() * 100) if available.sum() > 0 else 0
capacity_factor = (dispatched.sum() / (p_nom_sum * len(network.snapshots)) * 100)
curtailment_stats[carrier] = {
'capacity_mw': p_nom_sum,
'available_mwh': available.sum(),
'dispatched_mwh': dispatched.sum(),
'curtailed_mwh': curtailed.sum(),
'curtailment_pct': curtailment_pct,
'capacity_factor': capacity_factor
}
# Update layout
fig.update_xaxes(title_text='Time', row=len(renewable_carriers), col=1)
fig.update_layout(
title=dict(text='Renewable Curtailment Analysis - EE50_clustered', x=0.5, xanchor='center'),
height=900,
hovermode='x unified',
template='plotly_white'
)
fig.show()
# Print statistics
print('\nRenewable Curtailment Statistics:')
print('=' * 80)
for carrier, stats in curtailment_stats.items():
print(f"\n{carrier.replace('_', ' ').title()}:")
print(f" Installed Capacity: {stats['capacity_mw']:,.0f} MW")
print(f" Available Energy: {stats['available_mwh']:,.0f} MWh")
print(f" Dispatched Energy: {stats['dispatched_mwh']:,.0f} MWh")
print(f" Curtailed Energy: {stats['curtailed_mwh']:,.0f} MWh ({stats['curtailment_pct']:.1f}%)")
print(f" Capacity Factor: {stats['capacity_factor']:.1f}%")
Renewable Curtailment Statistics:
================================================================================
Wind Onshore:
Installed Capacity: 50,518 MW
Available Energy: 2,216,344 MWh
Dispatched Energy: 2,215,468 MWh
Curtailed Energy: 876 MWh (0.0%)
Capacity Factor: 26.1%
Wind Offshore:
Installed Capacity: 96,365 MW
Available Energy: 4,939,029 MWh
Dispatched Energy: 2,398,144 MWh
Curtailed Energy: 2,540,886 MWh (51.4%)
Capacity Factor: 14.8%
Solar Pv:
Installed Capacity: 99,807 MW
Available Energy: 1,924,436 MWh
Dispatched Energy: 1,923,031 MWh
Curtailed Energy: 1,406 MWh (0.1%)
Capacity Factor: 11.5%
Advanced Analysis#
Temporal Patterns and System Dynamics#
[29]:
# Interactive heatmap of hourly generation patterns
# Aggregate generation by hour of day and day of year
if len(network.snapshots) > 24:
total_gen = network.generators_t.p.sum(axis=1)
total_demand = network.loads_t.p_set.sum(axis=1)
# Create DataFrame with datetime index
df_temporal = pd.DataFrame({
'generation': total_gen,
'demand': total_demand,
'hour': total_gen.index.hour,
'day': total_gen.index.dayofyear
})
# Create hourly average pattern
hourly_pattern = df_temporal.groupby('hour')[['generation', 'demand']].mean()
# Interactive line plot of hourly patterns
fig = go.Figure()
fig.add_trace(go.Scatter(
x=hourly_pattern.index,
y=hourly_pattern['demand'] / 1000,
name='Average Demand',
line=dict(color='blue', width=2),
hovertemplate='Hour: %{x}<br>Demand: %{y:.1f} GW<extra></extra>'
))
fig.add_trace(go.Scatter(
x=hourly_pattern.index,
y=hourly_pattern['generation'] / 1000,
name='Average Generation',
line=dict(color='red', width=2),
hovertemplate='Hour: %{x}<br>Generation: %{y:.1f} GW<extra></extra>'
))
fig.update_layout(
title='Average Hourly Pattern',
xaxis_title='Hour of Day',
yaxis_title='Power (GW)',
template='plotly_white',
hovermode='x unified',
height=400
)
fig.show()
print(f'Peak demand hour: {hourly_pattern["demand"].idxmax()}:00')
print(f'Minimum demand hour: {hourly_pattern["demand"].idxmin()}:00')
else:
print('Insufficient timesteps for temporal pattern analysis')
Peak demand hour: 17:00
Minimum demand hour: 4:00
Congestion Hotspots Over Time#
Identify when and where transmission congestion occurs.
[30]:
# Calculate transmission loading over time (handles both lines and links)
has_lines = len(network.lines) > 0 and 's_nom' in network.lines.columns
has_internal_links = len(network.links) > 0
if has_lines:
# ETYS/Reduced network - use lines
loading_pct = (network.lines_t.p0.abs() / network.lines.s_nom * 100).fillna(0)
components = network.lines
element_type = 'Line'
elif has_internal_links:
# Zonal network - use internal links only (exclude external and hydrogen links)
h2_carriers = ['electrolysis', 'H2_turbine', 'H2_power', 'H2', 'H2_gas']
is_external = (
network.links.bus0.str.lower().str.contains('external', na=False) |
network.links.bus1.str.lower().str.contains('external', na=False)
)
is_h2_link = network.links.carrier.isin(h2_carriers)
internal_links = network.links[~is_external & ~is_h2_link]
if len(internal_links) > 0:
link_flows = network.links_t.p0[internal_links.index].abs()
loading_pct = (link_flows / internal_links.p_nom * 100).fillna(0)
components = internal_links
element_type = 'Link'
else:
loading_pct = pd.DataFrame()
element_type = 'None'
else:
loading_pct = pd.DataFrame()
element_type = 'None'
if len(loading_pct.columns) > 0:
# Count congested elements at each timestep (>95% loading)
congested_count = (loading_pct >= 95).sum(axis=1)
# Create time series plot
fig = go.Figure()
fig.add_trace(go.Scatter(
x=congested_count.index,
y=congested_count.values,
name=f'Congested {element_type}s',
line=dict(color='#DC143C', width=1.5),
fill='tozeroy',
fillcolor='rgba(220, 20, 60, 0.2)',
hovertemplate='Time: %{x}<br>Congested: %{y}<extra></extra>'
))
fig.update_layout(
title=f'Transmission Congestion Over Time ({element_type}s >95% Loaded)',
xaxis_title='Time',
yaxis_title=f'Number of Congested {element_type}s',
template='plotly_white',
height=400
)
fig.show()
# Find most frequently congested elements
congestion_frequency = (loading_pct >= 95).sum(axis=0).sort_values(ascending=False)
top_congested = congestion_frequency.head(10)
if len(top_congested) > 0 and top_congested.iloc[0] > 0:
print(f'\nMost Frequently Congested {element_type}s:')
print(f'{element_type}'.ljust(50) + ' Congested Timesteps')
print('=' * 75)
for elem_id, count in top_congested.items():
if count > 0:
bus0 = components.loc[elem_id, 'bus0']
bus1 = components.loc[elem_id, 'bus1']
pct = count / len(network.snapshots) * 100
print(f'{bus0} - {bus1:<40} {count:>4} ({pct:.1f}%)')
else:
print(f'\n✓ No significant transmission congestion detected')
else:
print('No transmission elements found for congestion analysis')
Most Frequently Congested Lines:
Line Congested Timesteps
===========================================================================
cluster_33 - cluster_89 168 (100.0%)
cluster_49 - cluster_50 168 (100.0%)
cluster_25 - cluster_6 159 (94.6%)
cluster_42 - cluster_48 153 (91.1%)
cluster_62 - cluster_72 152 (90.5%)
cluster_32 - cluster_95 151 (89.9%)
cluster_56 - cluster_84 128 (76.2%)
cluster_49 - cluster_82 126 (75.0%)
cluster_52 - cluster_79 117 (69.6%)
cluster_2 - cluster_73 113 (67.3%)
Key Performance Indicators#
[31]:
# Calculate KPIs
print(f'SCENARIO: {scenario_id}')
print('=' * 80)
print()
total_demand = network.loads_t.p_set.sum().sum()
total_generation = network.generators_t.p.sum().sum()
# Load shedding
ls_gens = network.generators[network.generators.carrier == 'load_shedding']
load_shedding = network.generators_t.p[ls_gens.index].sum().sum() if len(ls_gens) > 0 else 0
print(f'DEMAND:')
print(f' Total demand: {total_demand:,.0f} MWh')
print(f' Load shedding: {load_shedding:,.0f} MWh ({load_shedding/total_demand*100:.2f}%)')
print(f' Demand satisfied: {total_demand - load_shedding:,.0f} MWh ({(total_demand - load_shedding)/total_demand*100:.2f}%)')
print()
# Renewable share
renewable_gens = network.generators[network.generators.carrier.isin(['wind_onshore', 'wind_offshore', 'solar_pv', 'large_hydro', 'small_hydro'])]
renewable_gen = network.generators_t.p[renewable_gens.index].sum().sum()
renewable_share = renewable_gen / (total_generation + load_shedding) * 100
print(f'GENERATION:')
print(f' Total generation: {total_generation:,.0f} MWh')
print(f' Renewable generation: {renewable_gen:,.0f} MWh')
print(f' Renewable share: {renewable_share:.1f}%')
print()
# System cost
if hasattr(network, 'objective'):
print(f'OPTIMIZATION:')
print(f' Total system cost: £{network.objective:,.2f}')
if total_generation + load_shedding > 0:
cost_per_mwh = network.objective / (total_generation + load_shedding)
print(f' Cost per MWh: £{cost_per_mwh:.2f}/MWh')
print()
# Line statistics
max_loading = abs(network.lines_t.p0 / network.lines.s_nom).max().max()
mean_loading = abs(network.lines_t.p0 / network.lines.s_nom).mean().mean()
n_congested = (abs(network.lines_t.p0 / network.lines.s_nom) >= 0.95).sum().sum()
print(f'TRANSMISSION:')
print(f' Total lines: {len(network.lines)}')
print(f' Max line loading: {max_loading*100:.1f}%')
print(f' Mean line loading: {mean_loading*100:.1f}%')
print(f' Timesteps with congestion (≥95%): {n_congested}')
print()
# Hydrogen system
has_electrolysis = len(network.links) > 0 and 'electrolysis' in network.links.carrier.values
has_h2_turbines = len(network.links) > 0 and 'H2_turbine' in network.links.carrier.values
if has_electrolysis or has_h2_turbines:
print(f'HYDROGEN SYSTEM:')
if has_electrolysis:
electrolysis = network.links[network.links.carrier == 'electrolysis']
elec_consumed = network.links_t.p0[electrolysis.index].sum().sum()
print(f' Electrolysis capacity: {electrolysis.p_nom.sum():,.0f} MW')
print(f' Electricity consumed: {elec_consumed:,.0f} MWh')
if has_h2_turbines:
h2_turbines = network.links[network.links.carrier == 'H2_turbine']
elec_generated = -network.links_t.p1[h2_turbines.index].sum().sum()
print(f' H2 turbine capacity: {h2_turbines.p_nom.sum():,.0f} MW')
print(f' Electricity generated: {elec_generated:,.0f} MWh')
if has_electrolysis and has_h2_turbines:
roundtrip = elec_generated / elec_consumed * 100 if elec_consumed > 0 else 0
print(f' Round-trip efficiency: {roundtrip:.1f}%')
print()
print('=' * 80)
SCENARIO: EE50_clustered
================================================================================
DEMAND:
Total demand: 13,231,717 MWh
Load shedding: 14,869 MWh (0.11%)
Demand satisfied: 13,216,848 MWh (99.89%)
GENERATION:
Total generation: 13,975,556 MWh
Renewable generation: 6,680,638 MWh
Renewable share: 47.8%
OPTIMIZATION:
Total system cost: £193,005,907.38
Cost per MWh: £13.80/MWh
TRANSMISSION:
Total lines: 202
Max line loading: 100.0%
Mean line loading: 22.9%
Timesteps with congestion (≥95%): 2054
HYDROGEN SYSTEM:
Electrolysis capacity: 28,946 MW
Electricity consumed: 1,046,447 MWh
H2 turbine capacity: 27,520 MW
Electricity generated: 363,945 MWh
Round-trip efficiency: 34.8%
================================================================================