Market Dispatch#

This tutorial explains PyPSA-GB’s two-stage market workflow using the 2020 validation notebook as a worked example. It follows the model from the copperplate wholesale schedule through constrained balancing redispatch, then compares the result with ELEXON, NESO, and ESPENI evidence.

The notebook is intentionally analysis-first. It does not build the networks itself; it reads the market outputs produced by Snakemake and shows how to interpret them.

What This Tutorial Covers#

The market workflow answers a different question from the standard network-constrained OPF tutorials. The standard solve asks for the least-cost physical dispatch subject to network limits. The market workflow first asks what the system would schedule in an unconstrained GB wholesale market, then asks what redispatch is required once the real network constraints are restored.

By the end of the tutorial you should be able to:

  • identify the files produced by the wholesale and balancing stages;

  • explain why the wholesale price is read from the exported demand-bus price series;

  • compare wholesale dispatch with final physical dispatch;

  • interpret BM redispatch volumes, constraint costs, and congestion diagnostics;

  • use historical ELEXON and NESO validation outputs to judge whether a market scenario is behaving plausibly.

For configuration details, see the user guide page on market dispatch.

Before Running The Notebook#

This notebook expects the monthly 2020 validation market scenarios to have already been solved. The setup cell checks for the required outputs and raises a clear FileNotFoundError if any are missing.

The key output groups are:

File pattern

Produced by

Meaning

resources/market/{scenario}_wholesale.nc

solve_wholesale_market

Copperplate wholesale network

resources/market/{scenario}_wholesale_price.csv

solve_wholesale_market

Demand-bus wholesale SMP series

resources/market/{scenario}_balancing.nc

solve_balancing_mechanism

Physical constrained network

resources/market/{scenario}_redispatch_summary.csv

solve_balancing_mechanism

Asset-level BM increase/decrease volumes and costs

resources/market/{scenario}_bm_validation.csv

validate_bm_results

ELEXON BM validation metrics

resources/market/{scenario}_neso_validation.csv

validate_neso_constraints

NESO constraint validation metrics

Typical targets are:

snakemake resources/analysis/Validation_Jan2020_UniformNetworkBaseline_bm_validation.html -j 4
snakemake resources/analysis/Validation_Jan2020_UniformNetworkBaseline_neso_validation.html -j 4

The notebook aggregates all twelve monthly validation scenarios listed in the setup cell. If you only have one month, edit MONTH_SCENARIOS and MONTH_ORDER before running the rest of the notebook.

[1]:

from functools import lru_cache
from pathlib import Path
import json

import numpy as np
import pandas as pd
import plotly.graph_objects as go
import plotly.express as px
import plotly.io as pio
import pypsa
import yaml
from IPython.display import HTML, display
from plotly.subplots import make_subplots

ROOT = Path.cwd()
if not (ROOT / "resources" / "market").exists():
    for parent in ROOT.parents:
        if (parent / "resources" / "market").exists():
            ROOT = parent
            break

MARKET = ROOT / "resources" / "market"
DATA = ROOT / "data"

MONTH_SCENARIOS = [
    "Validation_Jan2020_UniformNetworkBaseline",
    "Validation_Feb2020_UniformNetworkBaseline",
    "Validation_Mar2020_UniformNetworkBaseline",
    "Validation_Apr2020_UniformNetworkBaseline",
    "Validation_May2020_UniformNetworkBaseline",
    "Validation_Jun2020_UniformNetworkBaseline",
    "Validation_Jul2020_UniformNetworkBaseline",
    "Validation_Aug2020_UniformNetworkBaseline",
    "Validation_Sep2020_UniformNetworkBaseline",
    "Validation_Oct2020_UniformNetworkBaseline",
    "Validation_Nov2020_UniformNetworkBaseline",
    "Validation_Dec2020_UniformNetworkBaseline",
]
MONTH_ORDER = [
    "Jan 2020", "Feb 2020", "Mar 2020", "Apr 2020", "May 2020", "Jun 2020",
    "Jul 2020", "Aug 2020", "Sep 2020", "Oct 2020", "Nov 2020", "Dec 2020",
]
SCENARIO_TO_PERIOD = dict(zip(MONTH_SCENARIOS, MONTH_ORDER))
SCENARIO_EXAMPLE = "Validation_Dec2020_UniformNetworkBaseline"
BOUNDARY_CASES = ["SCOTEX", "SSE-SP", "SSHARN::SSHARN_SPARSE_GATE_FIT", "ESTEX", "SEIMP", "SWALEX"]

pd.options.display.float_format = "{:,.3f}".format
pio.renderers.default = "notebook_connected"
plot_template = "plotly_white"

CARRIER_COLORS = {
    "CCGT": "#4C78A8", "OCGT": "#72B7B2", "coal": "#3B3B3B",
    "nuclear": "#F58518", "biomass": "#54A24B", "waste_to_energy": "#8CD17D",
    "wind": "#59A14F", "solar": "#EDC948", "hydro_non_ps": "#1F77B4",
    "pumped storage": "#9467BD", "Battery": "#B07AA1", "load_shedding": "#D62728",
    "other": "#C7C7C7", "unknown": "#C7C7C7",
}

def gbp_m(value):
    return f"GBP {value / 1e6:,.1f}m"

def ratio(value):
    return f"{value:,.2f}x" if pd.notna(value) else "n/a"

def style_plot(fig, height=430):
    fig.update_layout(
        template=plot_template,
        height=height,
        margin=dict(l=60, r=30, t=70, b=55),
        font=dict(size=12),
        legend_title_text="",
    )
    return fig

def _scenario_paths(scenario):
    return {
        "wholesale": MARKET / f"{scenario}_wholesale.nc",
        "balancing": MARKET / f"{scenario}_balancing.nc",
        "wholesale_price": MARKET / f"{scenario}_wholesale_price.csv",
        "boav_by_carrier": MARKET / f"{scenario}_boav_by_carrier.csv",
    }

missing = [
    str(path.relative_to(ROOT))
    for scenario in MONTH_SCENARIOS
    for key, path in _scenario_paths(scenario).items()
    if key in {"wholesale", "balancing", "wholesale_price"} and not path.exists()
]
if missing:
    raise FileNotFoundError("Missing solved outputs:\n" + "\n".join(missing))

@lru_cache(maxsize=4)
def load_network(path):
    return pypsa.Network(path)

def wholesale_network(scenario):
    return load_network(str(_scenario_paths(scenario)["wholesale"]))

def balancing_network(scenario):
    return load_network(str(_scenario_paths(scenario)["balancing"]))

def _snapshot_hours(index):
    index = pd.DatetimeIndex(index)
    if len(index) <= 1:
        return pd.Series(1.0, index=index, dtype=float)
    dt = pd.Series(index=index, data=index.to_series().diff().dt.total_seconds() / 3600.0)
    inferred = dt.dropna().median()
    if not np.isfinite(inferred) or inferred <= 0:
        inferred = 1.0
    dt.iloc[0] = inferred
    return dt.fillna(inferred).clip(lower=0).astype(float)

def _carrier_compare(carrier):
    carrier = str(carrier)
    if carrier in {"wind_onshore", "wind_offshore", "embedded_wind"}:
        return "wind"
    if carrier in {"solar_pv", "embedded_solar"}:
        return "solar"
    if carrier in {"large_hydro", "small_hydro"}:
        return "hydro_non_ps"
    if carrier in {"Pumped Storage Hydroelectricity", "pumped_hydro"}:
        return "pumped storage"
    return carrier

def _load_wholesale_price(scenario):
    path = _scenario_paths(scenario)["wholesale_price"]
    df = pd.read_csv(path, index_col=0, parse_dates=True)
    if "wholesale_price" not in df.columns:
        raise ValueError(f"{path} does not contain a wholesale_price column")
    return df["wholesale_price"].astype(float).rename("Model SMP")

def _load_mid_prices():
    path = MARKET / "elexon" / "mid_prices_2020.csv"
    if not path.exists():
        return pd.Series(dtype=float, name="ELEXON MID")
    df = pd.read_csv(path, parse_dates=["datetime"])
    price_col = "mid_price" if "mid_price" in df.columns else df.columns[-1]
    return df.set_index("datetime")[price_col].astype(float).rename("ELEXON MID")

def _load_neso_constraint_monthly():
    frames = []
    csv_path = DATA / "validation" / "constraint_breakdown_2019-2020.csv"
    if csv_path.exists():
        frames.append(pd.read_csv(csv_path, parse_dates=["Date"]))
    json_path = DATA / "validation" / "constraint_breakdown_2020_2021.json"
    if json_path.exists():
        payload = json.loads(json_path.read_text(encoding="utf-8-sig"))
        records = payload.get("result", {}).get("records", [])
        if records:
            frames.append(pd.DataFrame(records).assign(Date=lambda d: pd.to_datetime(d["Date"])))
    if not frames:
        return pd.DataFrame(columns=["period", "neso_thermal_cost_gbp", "neso_thermal_volume_abs_mwh", "neso_thermal_volume_signed_mwh"])
    df = pd.concat(frames, ignore_index=True)
    df = df[(df["Date"] >= "2020-01-01") & (df["Date"] < "2021-01-01")].copy()
    df["period"] = pd.Categorical(df["Date"].dt.strftime("%b %Y"), MONTH_ORDER, ordered=True)
    df["Thermal constraints cost"] = pd.to_numeric(df["Thermal constraints cost"], errors="coerce").fillna(0.0)
    df["Thermal constraints volume"] = pd.to_numeric(df["Thermal constraints volume"], errors="coerce").fillna(0.0)
    monthly = (
        df.groupby("period", observed=True)
        .agg(
            neso_thermal_cost_gbp=("Thermal constraints cost", "sum"),
            neso_thermal_volume_signed_mwh=("Thermal constraints volume", "sum"),
        )
        .reset_index()
    )
    abs_monthly = df.assign(abs_volume=df["Thermal constraints volume"].abs()).groupby("period", observed=True)["abs_volume"].sum()
    monthly["neso_thermal_volume_abs_mwh"] = monthly["period"].map(abs_monthly).astype(float)
    return monthly

def _load_boav_by_carrier():
    frames = []
    for scenario in MONTH_SCENARIOS:
        path = _scenario_paths(scenario)["boav_by_carrier"]
        if not path.exists():
            continue
        df = pd.read_csv(path)
        df["scenario"] = scenario
        df["period"] = SCENARIO_TO_PERIOD[scenario]
        df["carrier_compare"] = df["carrier"].map(_carrier_compare)
        for col in ["increase_mwh", "decrease_mwh"]:
            df[col] = pd.to_numeric(df[col], errors="coerce").fillna(0.0)
        frames.append(df)
    if not frames:
        return pd.DataFrame(columns=["scenario", "period", "carrier", "carrier_compare", "increase_mwh", "decrease_mwh"])
    out = pd.concat(frames, ignore_index=True)
    out["period"] = pd.Categorical(out["period"], MONTH_ORDER, ordered=True)
    return out

def _load_boundary_mapping():
    path = DATA / "network" / "neso_boundary_mapping_ssesp_strict.yaml"
    with open(path, "r", encoding="utf-8") as fh:
        raw = yaml.safe_load(fh)
    flat = {}
    for name, info in raw.get("boundaries", {}).items():
        item = dict(info)
        item.setdefault("neso_boundary", name)
        flat[name] = item
        for sub_name, sub_info in (info.get("subconstraints") or {}).items():
            sub = dict(sub_info)
            sub["parent_boundary"] = name
            sub["neso_boundary"] = name
            flat[f"{name}::{sub_name}"] = sub
    return {"boundaries": flat}

def _boundary_model_flow(network, mapping, boundary):
    spec = mapping["boundaries"][boundary]
    groups = spec.get("flow_groups", {}) or {}
    pos = groups.get("positive", spec.get("lines", [])) or []
    neg = groups.get("negative", []) or []
    flow = pd.Series(0.0, index=network.snapshots)
    for line in pos:
        if line in network.lines_t.p0.columns:
            flow = flow.add(network.lines_t.p0[line], fill_value=0.0)
    for line in neg:
        if line in network.lines_t.p0.columns:
            flow = flow.sub(network.lines_t.p0[line], fill_value=0.0)
    for link in spec.get("links", []) or []:
        name = link.get("name") if isinstance(link, dict) else link
        sign = float(link.get("sign", 1.0)) if isinstance(link, dict) else 1.0
        if name in network.links_t.p0.columns:
            flow = flow.add(sign * network.links_t.p0[name], fill_value=0.0)
    return flow.abs().rename(boundary)

def _load_neso_boundary_raw():
    path = DATA / "validation" / "day_ahead_constraint_flows_limits.csv"
    if not path.exists():
        return pd.DataFrame()
    df = pd.read_csv(path, low_memory=False)
    df["datetime"] = pd.to_datetime(df["Date (GMT/BST)"], errors="coerce")
    df = df[(df["datetime"] >= "2020-01-01") & (df["datetime"] < "2021-01-01")].copy()
    df["Limit (MW)"] = pd.to_numeric(df["Limit (MW)"], errors="coerce").abs()
    df["Flow (MW)"] = pd.to_numeric(df["Flow (MW)"], errors="coerce").abs()
    df["hour"] = df["datetime"].dt.floor("h")
    return df

def _model_dispatch_by_carrier(network, scenario):
    dt = _snapshot_hours(network.snapshots)
    records = []
    gen_p = network.generators_t.p
    common = gen_p.columns.intersection(network.generators.index)
    if len(common):
        energy = gen_p.loc[:, common].clip(lower=0).mul(dt, axis=0).sum()
        carriers = network.generators.loc[common, "carrier"].map(_carrier_compare)
        for carrier, value in energy.groupby(carriers).sum().items():
            records.append({"scenario": scenario, "period": SCENARIO_TO_PERIOD[scenario], "comparison": carrier, "model_mwh": value})
    su_p = network.storage_units_t.p
    common_su = su_p.columns.intersection(network.storage_units.index)
    if len(common_su):
        energy = su_p.loc[:, common_su].clip(lower=0).mul(dt, axis=0).sum()
        carriers = network.storage_units.loc[common_su, "carrier"].map(_carrier_compare)
        for carrier, value in energy.groupby(carriers).sum().items():
            records.append({"scenario": scenario, "period": SCENARIO_TO_PERIOD[scenario], "comparison": carrier, "model_mwh": value})
    demand_mwh = network.loads_t.p_set.mul(dt, axis=0).sum().sum()
    records.append({"scenario": scenario, "period": SCENARIO_TO_PERIOD[scenario], "comparison": "model_demand_profile", "model_mwh": demand_mwh})
    out = pd.DataFrame(records)
    out["demand_mwh"] = demand_mwh
    out["share_of_demand"] = out["model_mwh"] / demand_mwh if demand_mwh else np.nan
    return out

def _redispatch_summary(wholesale, balancing, scenario):
    dt = _snapshot_hours(wholesale.snapshots)
    records = []
    def add_component(wholesale_p, physical_p, components, component_type):
        common = wholesale_p.columns.intersection(physical_p.columns).intersection(components.index)
        for name in common:
            diff = physical_p[name].reindex(dt.index).fillna(0.0) - wholesale_p[name].reindex(dt.index).fillna(0.0)
            inc = diff.clip(lower=0).mul(dt).sum()
            dec = (-diff).clip(lower=0).mul(dt).sum()
            mc = pd.to_numeric(components.loc[name].get("marginal_cost", 0.0), errors="coerce")
            mc = 0.0 if pd.isna(mc) else float(mc)
            carrier = str(components.loc[name].get("carrier", "unknown"))
            records.append({
                "scenario": scenario, "period": SCENARIO_TO_PERIOD[scenario],
                "component": name, "type": component_type, "carrier": carrier,
                "carrier_compare": _carrier_compare(carrier),
                "increase_mwh": float(inc), "decrease_mwh": float(dec),
                "gross_mwh": float(inc + dec), "proxy_cost_gbp": float((inc + dec) * abs(mc)),
            })
    add_component(wholesale.generators_t.p, balancing.generators_t.p, balancing.generators, "generator")
    if not wholesale.storage_units_t.p.empty and not balancing.storage_units_t.p.empty:
        add_component(wholesale.storage_units_t.p, balancing.storage_units_t.p, balancing.storage_units, "storage_unit")
    return pd.DataFrame(records)

def _load_shedding_mwh(network):
    load_shedding = network.generators.index[network.generators["carrier"].astype(str).eq("load_shedding")]
    if len(load_shedding) == 0:
        return 0.0
    cols = network.generators_t.p.columns.intersection(load_shedding)
    if len(cols) == 0:
        return 0.0
    return float(network.generators_t.p.loc[:, cols].clip(lower=0).mul(_snapshot_hours(network.snapshots), axis=0).sum().sum())

def _load_espeni_dispatch():
    path = DATA / "demand" / "espeni.csv"
    if not path.exists():
        return pd.DataFrame(columns=["period", "comparison", "espeni_mwh"])
    header = pd.read_csv(path, nrows=0).columns
    date_col = "ELEC_elex_startTime[utc](datetime)"
    candidates = {
        "CCGT": ["ELEC_POWER_ELEX_CCGT[MW](float32)"],
        "coal": ["ELEC_POWER_ELEX_COAL[MW](float32)"],
        "nuclear": ["ELEC_POWER_ELEX_NUCLEAR[MW](float32)"],
        "OCGT": ["ELEC_POWER_ELEX_OCGT[MW](float32)"],
        "oil": ["ELEC_POWER_ELEX_OIL[MW](float32)"],
        "biomass": ["ELEC_POWER_ELEX_BIOMASS_POSTCALC[MW](float32)", "ELEC_POWER_ELEX_BIOMASS_PRECALC[MW](float32)"],
        "hydro_non_ps": ["ELEC_POWER_ELEX_NPSHYD[MW](float32)"],
        "pumped storage": ["ELEC_POWER_ELEX_PS_DISCHARGING[MW](float32)"],
        "wind": ["ELEC_POWER_TOTAL_WIND[MW](float32)"],
        "solar": ["ELEC_POWER_NGEM_EMBEDDED_SOLAR_GENERATION[MW](float32)"],
        "model_demand_profile": ["ELEC_POWER_TOTAL_ESPENI[MW](float32)"],
    }
    col_for = {name: next((c for c in cols if c in header), None) for name, cols in candidates.items()}
    usecols = [date_col] + [c for c in col_for.values() if c]
    df = pd.read_csv(path, usecols=usecols, parse_dates=[date_col])
    idx = pd.to_datetime(df[date_col], utc=True).dt.tz_convert(None)
    df = df.set_index(idx).drop(columns=[date_col])
    df = df[(df.index >= "2020-01-01") & (df.index < "2021-01-01")]
    dt = _snapshot_hours(df.index)
    records = []
    for comparison, col in col_for.items():
        if not col:
            continue
        series = pd.to_numeric(df[col], errors="coerce").fillna(0.0).clip(lower=0.0)
        monthly = series.mul(dt, axis=0).groupby(series.index.strftime("%b %Y")).sum()
        for period, value in monthly.items():
            records.append({"period": period, "comparison": comparison, "espeni_mwh": value})
    out = pd.DataFrame(records)
    out["period"] = pd.Categorical(out["period"], MONTH_ORDER, ordered=True)
    return out

mid_prices = _load_mid_prices()
neso_constraints = _load_neso_constraint_monthly()
thermal_costs = neso_constraints.set_index("period")["neso_thermal_cost_gbp"]
thermal_volumes = neso_constraints.set_index("period")["neso_thermal_volume_abs_mwh"]
boav_by_carrier = _load_boav_by_carrier()
boundary_mapping = _load_boundary_mapping()
neso_boundary_raw = _load_neso_boundary_raw()
espeni_dispatch = _load_espeni_dispatch()

score_records = []
dispatch_parts = []
redispatch_parts = []
wholesale_price_parts = []
model_boundary_parts = {boundary: [] for boundary in BOUNDARY_CASES}

for scenario in MONTH_SCENARIOS:
    wholesale = wholesale_network(scenario)
    balancing = balancing_network(scenario)
    period = SCENARIO_TO_PERIOD[scenario]
    price = _load_wholesale_price(scenario)
    wholesale_price_parts.append(price)
    mid_aligned = mid_prices.reindex(price.index).interpolate(limit_direction="both") if len(mid_prices) else pd.Series(dtype=float)
    redispatch_month = _redispatch_summary(wholesale, balancing, scenario)
    redispatch_parts.append(redispatch_month)
    dispatch_parts.append(_model_dispatch_by_carrier(balancing, scenario))
    for boundary in BOUNDARY_CASES:
        model_boundary_parts[boundary].append(_boundary_model_flow(balancing, boundary_mapping, boundary))
    score_records.append({
        "scenario": scenario, "period": period,
        "bm_cost_gbp": redispatch_month["proxy_cost_gbp"].sum(),
        "gross_redispatch_mwh": redispatch_month["gross_mwh"].sum(),
        "increase_mwh": redispatch_month["increase_mwh"].sum(),
        "decrease_mwh": redispatch_month["decrease_mwh"].sum(),
        "model_one_sided_redispatch_mwh": max(redispatch_month["increase_mwh"].sum(), redispatch_month["decrease_mwh"].sum()),
        "load_shedding_mwh": _load_shedding_mwh(balancing),
        "mean_smp_gbp_mwh": price.mean(),
        "mean_mid_gbp_mwh": mid_aligned.reindex(price.index).mean() if len(mid_aligned) else np.nan,
        "neso_thermal_cost_gbp": thermal_costs.get(period, np.nan),
        "neso_thermal_volume_abs_mwh": thermal_volumes.get(period, np.nan),
    })

score = pd.DataFrame(score_records)
score["period"] = pd.Categorical(score["period"], MONTH_ORDER, ordered=True)
score = score.sort_values("period").reset_index(drop=True)
score["bm_cost_gbp_m"] = score["bm_cost_gbp"] / 1e6
score["neso_thermal_cost_gbp_m"] = score["neso_thermal_cost_gbp"] / 1e6
score["bm_vs_neso_thermal_ratio"] = score["bm_cost_gbp"] / score["neso_thermal_cost_gbp"].replace(0, np.nan)
score["smp_mid_ratio"] = score["mean_smp_gbp_mwh"] / score["mean_mid_gbp_mwh"].replace(0, np.nan)

redispatch = pd.concat(redispatch_parts, ignore_index=True)
dispatch = pd.concat(dispatch_parts, ignore_index=True)
wholesale_price_series = pd.concat(wholesale_price_parts).sort_index()
model_boundary_flows = {boundary: pd.concat(parts).sort_index() for boundary, parts in model_boundary_parts.items() if parts}
boav_monthly = (
    boav_by_carrier.groupby("period", as_index=False, observed=True)[["increase_mwh", "decrease_mwh"]].sum()
    if not boav_by_carrier.empty
    else pd.DataFrame(columns=["period", "increase_mwh", "decrease_mwh"])
)

boundary_records = []
for scenario in MONTH_SCENARIOS:
    period = SCENARIO_TO_PERIOD[scenario]
    for boundary in BOUNDARY_CASES:
        spec = boundary_mapping["boundaries"][boundary]
        neso_boundary = spec.get("neso_boundary", boundary)
        model = model_boundary_flows[boundary]
        model_month = model.loc[model.index.to_period("M") == pd.Period(period)]
        nb = neso_boundary_raw[neso_boundary_raw["Constraint Group"] == neso_boundary]
        nb = nb[nb["hour"].dt.to_period("M") == pd.Period(period)]
        hourly = nb.groupby("hour")[["Limit (MW)", "Flow (MW)"]].mean() if not nb.empty else pd.DataFrame()
        if hourly.empty or model_month.empty:
            boundary_records.append({"scenario": scenario, "period": period, "boundary": boundary, "neso_boundary": neso_boundary})
            continue
        model_aligned = model_month.reindex(hourly.index).interpolate(limit_direction="both")
        limits = hourly["Limit (MW)"].where(lambda s: s < 50000).dropna()
        model_limited = model_aligned.reindex(limits.index).dropna()
        limit_aligned = limits.reindex(model_limited.index)
        loading = model_limited / limit_aligned.replace(0, np.nan)
        exceed = (model_limited - limit_aligned).clip(lower=0)
        boundary_records.append({
            "scenario": scenario, "period": period, "boundary": boundary, "neso_boundary": neso_boundary,
            "model_mean_flow_mw": model_month.mean(), "model_max_flow_mw": model_month.max(),
            "neso_mean_flow_mw": hourly["Flow (MW)"].mean(), "neso_max_flow_mw": hourly["Flow (MW)"].max(),
            "neso_mean_limit_mw": limits.mean(), "model_boundary_mean_loading": loading.mean(),
            "model_boundary_hours_above_90": int((loading >= 0.9).sum()),
            "model_boundary_limit_exceed_hours": int((exceed > 1e-6).sum()),
            "model_boundary_max_limit_exceed_mw": exceed.max() if len(exceed) else 0.0,
            "flow_ratio": model_month.mean() / hourly["Flow (MW)"].mean() if hourly["Flow (MW)"].mean() else np.nan,
        })

boundaries = pd.DataFrame(boundary_records)
boundaries["period"] = pd.Categorical(boundaries["period"], MONTH_ORDER, ordered=True)
boundaries = boundaries.sort_values(["period", "boundary"]).reset_index(drop=True)

dispatch_compare = dispatch.merge(espeni_dispatch, on=["period", "comparison"], how="left")
dispatch_compare["model_espeni_ratio"] = dispatch_compare["model_mwh"] / dispatch_compare["espeni_mwh"].replace(0, np.nan)

INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Jan2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Jan2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Feb2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Feb2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Mar2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Mar2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Apr2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Apr2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_May2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_May2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Jun2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Jun2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Jul2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Jul2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Aug2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Aug2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Sep2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Sep2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Oct2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Oct2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Nov2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Nov2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Dec2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers
INFO:pypsa.network.io:New version 1.2.1 available! (Current: 1.0.7)
INFO:pypsa.network.io:Imported network 'Validation_Dec2020_UniformNetworkBaseline (Full)' has buses, carriers, generators, lines, links, loads, storage_units, transformers

1. Method and Raw Data Sources#

Two-Step Dispatch Method#

The validation uses a two-stage market approximation.

  1. Wholesale day-ahead dispatch solves GB as a copperplate system. Transmission limits are relaxed, historical interconnector positions are fixed, and the model finds the least-cost dispatch position for generators and storage. This stage produces the wholesale schedule and the reported system marginal price (SMP).

  2. Balancing mechanism redispatch restores the transmission network and NESO boundary constraints. The balancing solve is anchored to the wholesale schedule, so assets can move up or down from their wholesale positions to manage network constraints. The difference between balancing dispatch and wholesale dispatch is reported as model BM redispatch volume.

This is intended to test whether the model produces credible market and network outcomes, not to reproduce the full GB settlement process.

Raw Data Sources#

Area

Raw source

Use in the model or validation

Transmission network

NESO Electricity Ten Year Statement (ETYS) network data

GB transmission topology, substations, lines, transformers, voltage levels, and network areas.

Network constraints

NESO day-ahead boundary flow and limit publications

Validation benchmark for boundary utilisation and limit exceedance, including SSE-SP, SCOTEX, SSHARN, ESTEX, SEIMP, and SWALEX.

Thermal fleet

DESNZ Digest of UK Energy Statistics (DUKES), NESO TEC register, and BM Unit reference mapping where available

Historical thermal station capacities, technologies, locations, and mapping to the transmission network.

Renewable and storage assets

DESNZ Renewable Energy Planning Database (REPD)

Site-level renewable and storage capacity, status, location, technology classification, and spatial allocation.

Renewable availability

ERA5 weather data processed through the PyPSA/atlite renewable-profile workflow

Weather-driven wind and solar availability profiles, with monthly performance factors used for the validation runs.

Demand and physical generation

ESPENI

Historical demand, embedded wind/solar, interconnector flows, and observed technology-level physical dispatch used for end-of-notebook validation.

Interconnector capacities and locations

DESNZ DUKES and NESO Interconnector Register

Built interconnector capacities, connection points, commissioning status, and mapping to GB buses.

Wholesale price benchmark

ELEXON Market Index Data (MID)

Observed wholesale price benchmark for model SMP.

Balancing action benchmark

ELEXON BOAV accepted bid/offer volumes

Observed accepted BM turn-up and turn-down volumes, used by carrier where mappings are available.

Constraint cost and volume benchmark

NESO thermal constraint cost and volume publications

System-level benchmark for total thermal constraint cost and volume; this is not carrier-resolved.

Benchmark Interpretation#

The validation uses several public datasets that do not measure exactly the same operational concept.

Evidence stream

What it checks

Main caveat

SMP vs ELEXON MID

Wholesale price level and time-series shape

MID is an observed traded-price benchmark, while SMP is a model marginal price from the copperplate wholesale solve.

Model redispatch vs BOAV

BM movement scale and carrier composition

BOAV is accepted bid/offer action volume; mapping to model carriers is only as good as the BMU-to-asset mapping.

Model redispatch vs NESO thermal volume

System-level constraint-management scale

NESO thermal volume is not carrier-resolved, so it is not used in the carrier plot.

Boundary flows vs NESO limits

Spatial plausibility of network congestion

Boundary limits are operational limits, not static asset ratings; exceedances flag periods needing closer review.

Physical dispatch vs ESPENI

Final technology-level generation plausibility

ESPENI includes embedded wind and solar, so this is a dispatch comparison rather than a renewable-availability calibration target.

The BM cost shown here is a redispatch cost indicator reconstructed from exported network dispatch and marginal-cost data. It is useful for directional comparison, but it is not a settlement-grade reconstruction of the full ELEXON bid/offer ladder objective.

How To Read The Results#

Treat the plots as a consistency check across three levels of evidence:

  1. The wholesale section checks whether the unconstrained schedule and SMP are credible before physical network limits are applied.

  2. The BM section checks how much redispatch the model needs once network constraints are restored, and which carriers provide that movement.

  3. The boundary and ESPENI sections check whether the physical dispatch is plausible against independent operational data.

A mismatch in one chart is not automatically a model failure. For example, high redispatch can be consistent with a copperplate schedule that overuses generation behind a real transmission boundary. The useful signal is whether the wholesale schedule, congestion pattern, redispatch mix, and final physical dispatch tell a coherent story.

2. High-Level Validation Results#

This section is the executive view. The KPI cards summarise the model against observed 2020 benchmarks, and the overview chart keeps the market and network evidence together before the detailed sections.

For review, the most important checks are directional: whether the model has the right order of magnitude for redispatch, whether the price comparison is distorted by internal network artefacts, and whether boundary congestion appears in plausible locations and seasons.

[2]:

boav_total_inc = boav_monthly["increase_mwh"].sum() if len(boav_monthly) else np.nan
boav_total_dec = boav_monthly["decrease_mwh"].sum() if len(boav_monthly) else np.nan
kpis = {
    "Months covered": f"{len(score):.0f}",
    "Model one-sided redispatch": f"{score['model_one_sided_redispatch_mwh'].sum() / 1000:,.0f} GWh",
    "BOAV increase volume": f"{boav_total_inc / 1000:,.0f} GWh",
    "BOAV decrease volume": f"{boav_total_dec / 1000:,.0f} GWh",
    "NESO thermal volume": f"{score['neso_thermal_volume_abs_mwh'].sum() / 1000:,.0f} GWh",
    "Load shedding": f"{score['load_shedding_mwh'].sum():,.1f} MWh",
    "Mean SMP / MID": ratio(score["smp_mid_ratio"].mean()),
}
cards = "".join(
    f"""
    <div style="border:1px solid #d9dee7;border-radius:6px;padding:14px 16px;background:#fbfcfe">
      <div style="font-size:12px;color:#5d6778;text-transform:uppercase;letter-spacing:.04em">{name}</div>
      <div style="font-size:24px;font-weight:650;color:#172033;margin-top:4px">{value}</div>
    </div>
    """
    for name, value in kpis.items()
)
display(HTML(f'<div style="display:grid;grid-template-columns:repeat(3,minmax(0,1fr));gap:12px">{cards}</div>'))

fig = make_subplots(
    rows=2,
    cols=1,
    shared_xaxes=True,
    vertical_spacing=0.12,
    subplot_titles=(
        "System-level BM volume benchmarks",
        "Monthly wholesale price fit",
    ),
)
fig.add_bar(
    row=1,
    col=1,
    x=score["period"].astype(str),
    y=score["model_one_sided_redispatch_mwh"] / 1000,
    name="Model one-sided redispatch",
    marker_color="#345995",
)
fig.add_bar(
    row=1,
    col=1,
    x=score["period"].astype(str),
    y=score["neso_thermal_volume_abs_mwh"] / 1000,
    name="NESO thermal constraint volume",
    marker_color="#E26D5A",
)
fig.add_scatter(
    row=2,
    col=1,
    x=score["period"].astype(str),
    y=score["mean_smp_gbp_mwh"],
    mode="lines+markers",
    name="Model SMP",
    line=dict(color="#345995", width=3),
)
fig.add_scatter(
    row=2,
    col=1,
    x=score["period"].astype(str),
    y=score["mean_mid_gbp_mwh"],
    mode="lines+markers",
    name="ELEXON MID",
    line=dict(color="#E26D5A", width=3),
)
fig.update_layout(title="Validation Overview", barmode="group")
fig.update_yaxes(row=1, col=1, title_text="GWh")
fig.update_yaxes(row=2, col=1, title_text="GBP/MWh")
style_plot(fig, height=720).show()
Months covered
12
Model one-sided redispatch
5,516 GWh
BOAV increase volume
16,414 GWh
BOAV decrease volume
12,122 GWh
NESO thermal volume
5,609 GWh
Load shedding
0.1 MWh
Mean SMP / MID
1.00x

3. Wholesale Market Evidence#

The wholesale section checks the copperplate dispatch before network constraints are restored. The merit-order chart is a stack and marginal-cost diagnostic only; it deliberately does not include a mean SMP reference line. SMP is compared to ELEXON MID in the price-ratio and time-series charts.

The SMP series comes from the exported wholesale-price output rather than a raw average of all bus marginal prices, because internal/link buses can carry artefactual prices that are not representative of the GB demand-bus wholesale price.

[3]:

def plot_merit_order(scenario=SCENARIO_EXAMPLE):
    network = wholesale_network(scenario)
    gen = network.generators.copy()
    gen = gen[(gen["p_nom"] > 0) & (gen["carrier"] != "load_shedding")].copy()
    snapshots = network.snapshots

    static_pmax = gen["p_max_pu"].reindex(gen.index).fillna(1.0).astype(float)
    pmax = pd.DataFrame(
        np.tile(static_pmax.values, (len(snapshots), 1)),
        index=snapshots,
        columns=gen.index,
    )
    if not network.generators_t.p_max_pu.empty:
        tv_cols = [c for c in gen.index if c in network.generators_t.p_max_pu.columns]
        pmax.loc[:, tv_cols] = network.generators_t.p_max_pu.reindex(
            snapshots,
            columns=tv_cols,
        ).fillna(1.0)

    gen["available_mw"] = gen["p_nom"] * pmax.mean(axis=0).reindex(gen.index).fillna(1.0)
    gen = gen[gen["available_mw"] > 1.0].copy()
    gen["marginal_cost"] = pd.to_numeric(gen["marginal_cost"], errors="coerce").fillna(0.0)
    gen["carrier_plot"] = gen["carrier"].map(_carrier_compare).where(
        gen["carrier"].map(_carrier_compare).isin(CARRIER_COLORS),
        "other",
    )
    gen["mc_bin"] = gen["marginal_cost"].round(1)

    blocks = (
        gen.groupby(["carrier_plot", "mc_bin"], as_index=False)["available_mw"]
        .sum()
        .rename(columns={"carrier_plot": "carrier", "mc_bin": "marginal_cost"})
        .sort_values(["marginal_cost", "carrier"])
    )
    blocks["left_gw"] = blocks["available_mw"].cumsum().shift(fill_value=0) / 1000.0
    blocks["width_gw"] = blocks["available_mw"] / 1000.0
    blocks["center_gw"] = blocks["left_gw"] + blocks["width_gw"] / 2.0
    blocks["base"] = np.minimum(blocks["marginal_cost"], 0.0)
    blocks["height"] = blocks["marginal_cost"].abs().clip(lower=1.0)

    served = network.loads_t.p_set.sum(axis=1) / 1000.0
    fig = go.Figure()
    for carrier, group in blocks.groupby("carrier", sort=False):
        fig.add_bar(
            x=group["center_gw"],
            y=group["height"],
            width=group["width_gw"],
            base=group["base"],
            name=carrier,
            marker=dict(color=CARRIER_COLORS.get(carrier, CARRIER_COLORS["other"]), line=dict(width=0)),
            customdata=group["marginal_cost"],
            hovertemplate=(
                "Carrier: %{fullData.name}<br>"
                "Capacity block: %{width:.2f} GW<br>"
                "Marginal cost: %{customdata:.1f} GBP/MWh<extra></extra>"
            ),
        )

    fig.add_vline(
        x=float(served.mean()),
        line_color="#D62728",
        annotation_text=f"Mean load: {served.mean():.1f} GW",
        annotation_position="top left",
    )
    fig.add_vline(
        x=float(served.max()),
        line_dash="dot",
        line_color="#D62728",
        annotation_text=f"Peak load: {served.max():.1f} GW",
        annotation_position="top right",
    )
    fig.update_layout(
        title=f"Wholesale Merit Order - {scenario}",
        xaxis_title="Cumulative average available capacity (GW)",
        yaxis_title="Marginal cost (GBP/MWh)",
        barmode="overlay",
        bargap=0,
        legend=dict(orientation="h", y=-0.22, x=0),
        yaxis=dict(range=[-120, 220], zeroline=True, zerolinewidth=1, zerolinecolor="#333333"),
    )
    return style_plot(fig, height=650)


plot_merit_order().show()

fig = go.Figure()
fig.add_scatter(
    x=score["period"].astype(str),
    y=score["smp_mid_ratio"],
    mode="lines+markers+text",
    text=score["smp_mid_ratio"].map(lambda x: f"{x:.2f}x"),
    textposition="top center",
    name="SMP / MID",
    line=dict(color="#4C9F70", width=3),
)
fig.add_hline(y=1.0, line_dash="dash", line_color="#545b66", annotation_text="Parity with ELEXON MID")
fig.update_layout(title="Wholesale Price Ratio by Month")
fig.update_yaxes(title_text="Model SMP / ELEXON MID", range=[0, max(1.25, score["smp_mid_ratio"].max() * 1.15)])
style_plot(fig, height=430).show()

def plot_wholesale_price_trace():
    model_price = wholesale_price_series.rename("Model SMP")
    elexon_price = mid_prices.reindex(model_price.index).interpolate(limit_direction="both").rename("ELEXON MID")
    comparison = pd.concat([model_price, elexon_price], axis=1).dropna().loc["2020-01-01":"2020-12-31 23:30"]
    comparison["Delta"] = comparison["Model SMP"] - comparison["ELEXON MID"]
    corr = comparison["Model SMP"].corr(comparison["ELEXON MID"])
    mae = comparison["Delta"].abs().mean()
    bias = comparison["Delta"].mean()

    fig = make_subplots(rows=2, cols=1, shared_xaxes=True, row_heights=[0.72, 0.28], vertical_spacing=0.08)
    fig.add_scatter(row=1, col=1, x=comparison.index, y=comparison["ELEXON MID"], mode="lines", name=f"ELEXON MID mean: GBP {comparison['ELEXON MID'].mean():.1f}/MWh", line=dict(color="#F58518", width=1.5))
    fig.add_scatter(row=1, col=1, x=comparison.index, y=comparison["Model SMP"], mode="lines", name=f"Model SMP mean: GBP {comparison['Model SMP'].mean():.1f}/MWh", line=dict(color="#1F4E79", width=1.8))
    fig.add_scatter(row=2, col=1, x=comparison.index, y=comparison["Delta"].where(comparison["Delta"] >= 0), mode="lines", name="Model above ELEXON", line=dict(color="#4C78A8", width=0), fill="tozeroy", fillcolor="rgba(76,120,168,0.55)", connectgaps=False)
    fig.add_scatter(row=2, col=1, x=comparison.index, y=comparison["Delta"].where(comparison["Delta"] < 0), mode="lines", name="Model below ELEXON", line=dict(color="#E45756", width=0), fill="tozeroy", fillcolor="rgba(228,87,86,0.55)", connectgaps=False)
    fig.add_hline(row=2, col=1, y=0, line_color="#333333", line_width=1)
    fig.add_annotation(row=1, col=1, text=f"r = {corr:.2f}<br>MAE = GBP {mae:.1f}/MWh<br>Bias = GBP {bias:+.1f}/MWh", xref="paper", yref="paper", x=0.99, y=0.94, showarrow=False, align="right", bgcolor="rgba(255,255,255,0.85)", bordercolor="#DDDDDD", borderwidth=1)
    fig.update_layout(title="Wholesale Price: Model SMP vs ELEXON MID - 2020")
    fig.update_yaxes(row=1, col=1, title_text="GBP/MWh")
    fig.update_yaxes(row=2, col=1, title_text="Delta<br>GBP/MWh")
    fig.update_xaxes(row=2, col=1, title_text="2020")
    fig.update_xaxes(range=[pd.Timestamp("2020-01-01"), pd.Timestamp("2021-01-01")])
    return style_plot(fig, height=650)


plot_wholesale_price_trace().show()

4. Balancing Mechanism Evidence#

All BM evidence is grouped here so the model redispatch, BOAV actions, and NESO thermal constraint volumes can be read together. Model redispatch volume is computed from the change between wholesale and final balancing dispatch.

BOAV is used for accepted bid/offer action volume and can be compared by carrier after mapping. NESO thermal constraint volume is system-level, so it is included in the monthly volume comparison but not allocated across carriers.

[4]:

monthly = score[["period", "model_one_sided_redispatch_mwh", "neso_thermal_volume_abs_mwh"]].merge(boav_monthly, on="period", how="left")
monthly_plot = monthly.rename(columns={
    "model_one_sided_redispatch_mwh": "Model one-sided redispatch",
    "increase_mwh": "BOAV increase",
    "decrease_mwh": "BOAV decrease",
    "neso_thermal_volume_abs_mwh": "NESO thermal constraint volume",
}).melt("period", var_name="series", value_name="MWh")
monthly_plot["GWh"] = monthly_plot["MWh"] / 1000
fig = px.bar(
    monthly_plot,
    x=monthly_plot["period"].astype(str),
    y="GWh",
    color="series",
    barmode="group",
    labels={"x": "", "GWh": "GWh", "series": ""},
    title="BM Redispatch Volume vs BOAV and NESO Thermal Constraint Volume",
    template=plot_template,
    height=560,
)
fig.show()

carrier_order = ["CCGT", "coal", "pumped storage", "hydro_non_ps", "biomass", "wind", "OCGT", "nuclear", "unknown", "load_shedding"]
model_carrier = (
    redispatch.groupby("carrier_compare", as_index=False)[["increase_mwh", "decrease_mwh"]]
    .sum()
    .rename(columns={"carrier_compare": "carrier", "increase_mwh": "model_increase_mwh", "decrease_mwh": "model_decrease_mwh"})
)
boav_carrier = (
    boav_by_carrier.groupby("carrier_compare", as_index=False, observed=True)[["increase_mwh", "decrease_mwh"]]
    .sum()
    .rename(columns={"carrier_compare": "carrier", "increase_mwh": "boav_increase_mwh", "decrease_mwh": "boav_decrease_mwh"})
)
carrier_vol = model_carrier.merge(boav_carrier, on="carrier", how="outer").fillna(0.0)
carrier_vol["carrier"] = pd.Categorical(carrier_vol["carrier"], carrier_order, ordered=True)
carrier_vol = carrier_vol.sort_values("carrier")
for col in ["model_increase_mwh", "model_decrease_mwh", "boav_increase_mwh", "boav_decrease_mwh"]:
    carrier_vol[col.replace("_mwh", "_gwh")] = carrier_vol[col] / 1000

fig = make_subplots(rows=1, cols=2, shared_yaxes=True, horizontal_spacing=0.1, subplot_titles=("Increase / turn-up", "Decrease / turn-down"))
fig.add_bar(row=1, col=1, y=carrier_vol["carrier"].astype(str), x=carrier_vol["model_increase_gwh"], orientation="h", name="Model", marker_color="#4C78A8")
fig.add_bar(row=1, col=1, y=carrier_vol["carrier"].astype(str), x=carrier_vol["boav_increase_gwh"], orientation="h", name="BOAV", marker_color="#F58518")
fig.add_bar(row=1, col=2, y=carrier_vol["carrier"].astype(str), x=carrier_vol["model_decrease_gwh"], orientation="h", name="Model ", marker_color="#4C78A8", showlegend=False)
fig.add_bar(row=1, col=2, y=carrier_vol["carrier"].astype(str), x=carrier_vol["boav_decrease_gwh"], orientation="h", name="BOAV ", marker_color="#F58518", showlegend=False)
fig.update_layout(
    title=f"BM Volumes by Carrier: Model vs BOAV (NESO thermal volume is system-level: {score['neso_thermal_volume_abs_mwh'].sum() / 1000:,.0f} GWh)",
    barmode="group",
)
fig.update_xaxes(title_text="Volume (GWh)", row=1, col=1)
fig.update_xaxes(title_text="Volume (GWh)", row=1, col=2)
fig.update_yaxes(autorange="reversed", row=1, col=1)
style_plot(fig, height=650).show()

fig = go.Figure()
fig.add_bar(x=score["period"].astype(str), y=score["bm_cost_gbp_m"], name="Model redispatch cost indicator", marker_color="#345995")
fig.add_bar(x=score["period"].astype(str), y=score["neso_thermal_cost_gbp_m"], name="NESO thermal constraint cost", marker_color="#E26D5A")
fig.update_layout(title="Monthly Redispatch Cost Indicator Compared With NESO Real Thermal Constraint Cost", barmode="group")
fig.update_yaxes(title_text="GBP million")
style_plot(fig, height=430).show()

5. Network Boundary Evidence#

This is the physical network check. Boundary flows are recomputed from the monthly balancing networks using the strict NESO boundary mapping used for the uniform-baseline validation case, with particular attention to SSE-SP.

The plots are intended to show both whether constraints bind and whether any periods exceed the published boundary limits. A limit exceedance is a review flag: it can indicate a genuine model issue, a boundary-mapping issue, or a difference between static validation limits and operational limit setting.

[5]:

def plot_neso_boundary_case_study(selected=("SCOTEX", "SSE-SP", "SSHARN::SSHARN_SPARSE_GATE_FIT")):
    fig = make_subplots(rows=len(selected), cols=1, shared_xaxes=True, vertical_spacing=0.05, subplot_titles=[b.replace("::SSHARN_SPARSE_GATE_FIT", " sparse gate") for b in selected])
    for row, boundary in enumerate(selected, start=1):
        spec = boundary_mapping["boundaries"][boundary]
        neso_boundary = spec.get("neso_boundary", boundary)
        nb = neso_boundary_raw[neso_boundary_raw["Constraint Group"] == neso_boundary]
        hourly = nb.groupby("hour")[["Limit (MW)", "Flow (MW)"]].mean()
        model_flow = model_boundary_flows.get(boundary, pd.Series(dtype=float)).reindex(hourly.index).interpolate(limit_direction="both")
        limit_gw = (hourly["Limit (MW)"] / 1000.0).where(lambda s: s < 20.0)
        fig.add_scatter(row=row, col=1, x=hourly.index, y=limit_gw, mode="lines", name="NESO DA limit", line=dict(color="#9EC5F8", width=0), fill="tozeroy", fillcolor="rgba(158,197,248,0.45)", showlegend=row == 1)
        fig.add_scatter(row=row, col=1, x=hourly.index, y=hourly["Flow (MW)"] / 1000.0, mode="lines", name="NESO DA flow", line=dict(color="#F58518", width=1.4), showlegend=row == 1)
        fig.add_scatter(row=row, col=1, x=model_flow.index, y=model_flow / 1000.0, mode="lines", name="Model BM flow", line=dict(color="#1F4E79", width=1.8), showlegend=row == 1)
        constrained = hourly["Limit (MW)"] < 50000
        util = (hourly.loc[constrained, "Flow (MW)"] / hourly.loc[constrained, "Limit (MW)"]).replace([np.inf, -np.inf], np.nan)
        model_util = (model_flow.loc[hourly.index[constrained]] / hourly.loc[constrained, "Limit (MW)"]).replace([np.inf, -np.inf], np.nan)
        fig.add_annotation(row=row, col=1, text=f"NESO >=90%: {(util >= 0.9).mean() * 100:.0f}%<br>Model exceed hrs: {(model_util > 1.0).sum():.0f}", xref=f"x{row if row > 1 else ''} domain", yref=f"y{row if row > 1 else ''} domain", x=0.99, y=0.88, showarrow=False, align="right", font=dict(size=10, color="#333333"))
        fig.update_yaxes(row=row, col=1, title_text=f"{boundary.split('::')[0]}<br>GW")
    fig.update_layout(title="NESO Real Boundary Limits and Flows vs Model Boundary Flows - 2020")
    fig.update_xaxes(row=len(selected), col=1, title_text="2020")
    fig.update_xaxes(range=[pd.Timestamp("2020-01-01"), pd.Timestamp("2021-01-01")])
    return style_plot(fig, height=760)

plot_neso_boundary_case_study().show()

top_boundaries = boundaries.sort_values("model_boundary_hours_above_90", ascending=False).head(18).copy()
top_boundaries["label"] = top_boundaries["period"].astype(str) + " - " + top_boundaries["boundary"].str.replace("::SSHARN_SPARSE_GATE_FIT", " sparse", regex=False)
fig = px.bar(
    top_boundaries.sort_values("model_boundary_hours_above_90"),
    x="model_boundary_hours_above_90",
    y="label",
    color="model_boundary_limit_exceed_hours",
    orientation="h",
    color_continuous_scale="Reds",
    labels={"model_boundary_hours_above_90": "Model hours above 90% of NESO limit", "label": "", "model_boundary_limit_exceed_hours": "Exceed hours"},
    title="Highest Model Boundary Utilisation Cases Derived From Balancing Networks",
    template=plot_template,
    height=560,
)
fig.show()

6. ESPENI Dispatch Comparison#

ESPENI is used as an independent observed physical-dispatch benchmark after the market and boundary checks. This gives a final view of whether the solved balancing dispatch produces credible annual technology-level output.

Wind and solar include ESPENI embedded generation, so this is a dispatch comparison, not an availability calibration target. Large differences should be read alongside the renewable capacity, embedded-generation treatment, and curtailment evidence in the earlier sections.

[6]:

focus = ["CCGT", "coal", "nuclear", "biomass", "hydro_non_ps", "wind", "solar", "pumped storage", "model_demand_profile"]
heat = (
    dispatch_compare[dispatch_compare["comparison"].isin(focus)]
    .pivot_table(index="comparison", columns="period", values="model_espeni_ratio", aggfunc="sum", observed=True)
    .reindex(focus)
)
heat = heat[MONTH_ORDER]
fig = px.imshow(
    heat,
    color_continuous_scale="RdBu",
    zmin=0,
    zmax=2,
    aspect="auto",
    labels={"x": "", "y": "", "color": "Model / ESPENI"},
    title="Model Final Physical Dispatch vs ESPENI",
    template=plot_template,
    height=520,
)
fig.update_traces(text=np.round(heat.values, 2), texttemplate="%{text}")
fig.show()

annual_espeni = (
    dispatch_compare[dispatch_compare["comparison"].isin(focus)]
    .groupby("comparison", as_index=False)[["model_mwh", "espeni_mwh"]]
    .sum()
)
annual_espeni["model_espeni_ratio"] = annual_espeni["model_mwh"] / annual_espeni["espeni_mwh"].replace(0, np.nan)
display(annual_espeni.set_index("comparison").reindex(focus).round(2))
model_mwh espeni_mwh model_espeni_ratio
comparison
CCGT 88,403,686.200 94,990,067.500 0.930
coal 226,359.000 4,381,873.500 0.050
nuclear 47,138,094.630 47,407,780.000 0.990
biomass 19,669,898.210 17,997,149.500 1.090
hydro_non_ps 2,368,261.240 4,314,369.500 0.550
wind 76,562,155.890 72,943,684.000 1.050
solar 12,335,410.050 12,062,276.000 1.020
pumped storage 474,501.240 1,493,022.000 0.320
model_demand_profile 280,043,421.000 275,586,455.500 1.020

7. Appendix Tables#

The appendix keeps the monthly audit table visible for reviewers who want to trace the plotted KPIs back to the underlying month-level values.

[7]:

display_cols = [
    "period", "bm_cost_gbp_m", "neso_thermal_cost_gbp_m", "bm_vs_neso_thermal_ratio",
    "model_one_sided_redispatch_mwh", "neso_thermal_volume_abs_mwh",
    "load_shedding_mwh", "mean_smp_gbp_mwh", "mean_mid_gbp_mwh", "smp_mid_ratio",
]
table = score[display_cols].rename(columns={
    "period": "Month",
    "bm_cost_gbp_m": "Model redispatch cost indicator (GBPm)",
    "neso_thermal_cost_gbp_m": "NESO thermal cost (GBPm)",
    "bm_vs_neso_thermal_ratio": "Cost indicator / NESO cost",
    "model_one_sided_redispatch_mwh": "Model one-sided redispatch (MWh)",
    "neso_thermal_volume_abs_mwh": "NESO thermal volume (MWh)",
    "load_shedding_mwh": "Load shedding (MWh)",
    "mean_smp_gbp_mwh": "Model SMP (GBP/MWh)",
    "mean_mid_gbp_mwh": "ELEXON MID (GBP/MWh)",
    "smp_mid_ratio": "SMP / MID",
})
display(table.round(2))

display(boundaries[[
    "period", "boundary", "model_max_flow_mw", "neso_mean_limit_mw",
    "model_boundary_hours_above_90", "model_boundary_limit_exceed_hours",
    "model_boundary_max_limit_exceed_mw",
]].round(2))
Month Model redispatch cost indicator (GBPm) NESO thermal cost (GBPm) Cost indicator / NESO cost Model one-sided redispatch (MWh) NESO thermal volume (MWh) Load shedding (MWh) Model SMP (GBP/MWh) ELEXON MID (GBP/MWh) SMP / MID
0 Jan 2020 71.530 72.600 0.990 1,040,076.870 736,380.000 0.010 36.530 34.940 1.050
1 Feb 2020 58.540 76.150 0.770 844,840.450 779,093.000 0.000 34.040 29.490 1.150
2 Mar 2020 45.630 40.510 1.130 633,025.290 460,646.000 0.000 29.400 28.490 1.030
3 Apr 2020 18.600 10.700 1.740 234,057.150 122,679.000 0.000 24.220 23.050 1.050
4 May 2020 22.160 34.830 0.640 277,844.880 327,806.000 0.000 21.270 22.030 0.970
5 Jun 2020 29.640 33.380 0.890 385,869.270 473,286.000 0.010 28.340 25.520 1.110
6 Jul 2020 29.740 40.340 0.740 380,802.050 452,897.000 0.000 28.250 27.670 1.020
7 Aug 2020 21.210 17.390 1.220 268,024.690 178,959.000 0.000 33.050 35.230 0.940
8 Sep 2020 12.220 40.080 0.300 153,472.200 319,840.000 0.000 38.610 43.160 0.890
9 Oct 2020 22.670 54.220 0.420 288,996.150 452,850.000 0.000 39.800 41.470 0.960
10 Nov 2020 41.520 111.080 0.370 556,707.240 886,526.000 0.010 36.350 36.120 1.010
11 Dec 2020 34.270 64.770 0.530 452,166.550 417,607.000 0.010 46.720 54.050 0.860
period boundary model_max_flow_mw neso_mean_limit_mw model_boundary_hours_above_90 model_boundary_limit_exceed_hours model_boundary_max_limit_exceed_mw
0 Jan 2020 ESTEX 7,033.870 8,488.000 11 0 0.000
1 Jan 2020 SCOTEX 4,636.080 4,441.370 208 0 0.000
2 Jan 2020 SEIMP 7,141.510 7,427.860 3 0 0.000
3 Jan 2020 SSE-SP 3,036.900 2,557.430 60 0 0.000
4 Jan 2020 SSHARN::SSHARN_SPARSE_GATE_FIT 8,750.000 7,433.990 306 4 425.000
... ... ... ... ... ... ... ...
67 Dec 2020 SCOTEX 4,642.380 5,425.000 8 0 0.000
68 Dec 2020 SEIMP 7,081.180 6,021.670 136 23 2,281.180
69 Dec 2020 SSE-SP 2,468.710 2,333.330 12 0 0.000
70 Dec 2020 SSHARN::SSHARN_SPARSE_GATE_FIT 7,450.680 7,178.660 112 0 0.000
71 Dec 2020 SWALEX 2,263.120 NaN 0 0 0.000

72 rows × 7 columns