From 0f5598212422521f918c1ba6207a011f59208360 Mon Sep 17 00:00:00 2001 From: Kamran Date: Mon, 31 Aug 2026 19:11:35 -0700 Subject: [PATCH 01/13] instrument: add benchmark directives to the 11 rules missing them Ten rules in build_electricity.smk and add_sectors in build_sector.smk had no benchmark: directive, so snakemake recorded no wall time or max RSS for them. That included cluster_resources (the largest walltime budget in the config) and add_extra_components (the terminal rule for network-only builds), which are exactly the rules needed to tune Slurm resource requests. Benchmark paths mirror each rule's log path under BENCHMARKS, and carry every wildcard the rule takes so concurrent jobs of the same rule cannot collide. Co-Authored-By: Claude Opus 5 --- workflow/rules/build_electricity.smk | 21 +++++++++++++++++++++ workflow/rules/build_sector.smk | 3 +++ 2 files changed, 24 insertions(+) diff --git a/workflow/rules/build_electricity.smk b/workflow/rules/build_electricity.smk index 9274be168..cf5422079 100644 --- a/workflow/rules/build_electricity.smk +++ b/workflow/rules/build_electricity.smk @@ -35,6 +35,8 @@ rule build_shapes: county_shapes=GEOSPATIAL + "{interconnect}/county_shapes.geojson", log: "logs/build_shapes/{interconnect}.log", + benchmark: + BENCHMARKS + "{interconnect}/build_shapes" threads: 1 resources: walltime=config_provider("walltime", "build_shapes", default="00:30:00"), @@ -83,6 +85,8 @@ rule build_base_network: network=NETWORKS + "{interconnect}/elec_base_network.nc", log: "logs/create_network/{interconnect}.log", + benchmark: + BENCHMARKS + "{interconnect}/build_base_network" threads: 1 resources: mem_mb=5000, @@ -114,6 +118,8 @@ rule build_bus_regions: regions_offshore=GEOSPATIAL + "{interconnect}/regions_offshore.geojson", log: "logs/build_bus_regions/{interconnect}.log", + benchmark: + BENCHMARKS + "{interconnect}/build_bus_regions" threads: 1 resources: mem_mb=3000, @@ -137,6 +143,8 @@ rule build_cost_data: sector_costs=COSTS + "sector_costs_{year}.csv", log: LOGS + "costs_{year}.log", + benchmark: + BENCHMARKS + "build_cost_data_{year}" threads: 1 resources: mem_mb=5000, @@ -754,6 +762,8 @@ rule build_powerplants: powerplants="resources/powerplants/powerplants.csv", log: "logs/build_powerplants.log", + benchmark: + BENCHMARKS + "build_powerplants" resources: mem_mb=30000, walltime=config_provider("walltime", "build_powerplants", default="00:30:00"), @@ -889,6 +899,8 @@ rule aggregate_to_substations: busmap=BUSMAPS + "{interconnect}/busmap_b.csv", log: "logs/aggregate_to_substations/{interconnect}.log", + benchmark: + BENCHMARKS + "{interconnect}/aggregate_to_substations" threads: 1 resources: mem_mb=lambda wildcards, input, attempt: (input.size // 150000) * attempt * 1.5, @@ -938,6 +950,8 @@ rule cluster_resources: busmap=BUSMAPS + "{interconnect}/busmap_s{simpl}.csv", log: "logs/cluster_resources/{interconnect}/elec_s{simpl}.log", + benchmark: + BENCHMARKS + "{interconnect}/cluster_resources_elec_s{simpl}" threads: 1 resources: mem_mb=lambda wildcards, input, attempt: (input.size // 150000) * attempt * 1.5, @@ -956,6 +970,8 @@ rule build_servm_load_weights: weights=DEMAND + "{interconnect}/servm_load_weights_s{simpl}.csv", log: LOGS + "{interconnect}/build_servm_load_weights_s{simpl}.log", + benchmark: + BENCHMARKS + "{interconnect}/build_servm_load_weights_s{simpl}" threads: 1 resources: mem_mb=lambda wildcards, input, attempt: (input.size // 200000) * attempt * 2, @@ -1071,6 +1087,8 @@ rule add_extra_components: NETWORKS + "{interconnect}/elec_s{simpl}_c{clusters}_ec.nc", log: "logs/add_extra_components/{interconnect}/elec_s{simpl}_c{clusters}_ec.log", + benchmark: + BENCHMARKS + "{interconnect}/add_extra_components_elec_s{simpl}_c{clusters}_ec" threads: 1 resources: mem_mb=lambda wildcards, input, attempt: (input.size // 100000) * attempt * 2, @@ -1109,6 +1127,9 @@ rule prepare_network: NETWORKS + "{interconnect}/elec_s{simpl}_c{clusters}_ec_l{ll}_{opts}.nc", log: "logs/prepare_network/{interconnect}/elec_s{simpl}_c{clusters}_ec_l{ll}_{opts}.log", + benchmark: + BENCHMARKS + +"{interconnect}/prepare_network_elec_s{simpl}_c{clusters}_ec_l{ll}_{opts}" threads: 1 resources: walltime=config_provider("walltime", "prepare_network", default="00:30:00"), diff --git a/workflow/rules/build_sector.smk b/workflow/rules/build_sector.smk index 02b3a69bd..e93b11cce 100644 --- a/workflow/rules/build_sector.smk +++ b/workflow/rules/build_sector.smk @@ -57,6 +57,9 @@ rule add_sectors: + "{interconnect}/elec_s{simpl}_c{clusters}_ec_l{ll}_{opts}_{sector}.nc", log: "logs/add_sectors/{interconnect}/elec_s{simpl}_c{clusters}_ec_l{ll}_{opts}_{sector}.log", + benchmark: + BENCHMARKS + +"{interconnect}/add_sectors_elec_s{simpl}_c{clusters}_ec_l{ll}_{opts}_{sector}" group: "prepare" threads: 1 From 0dbc17a0ae3c2614efa4ef2b880462cd848d8cb6 Mon Sep 17 00:00:00 2001 From: Kamran Date: Mon, 31 Aug 2026 19:11:42 -0700 Subject: [PATCH 02/13] feat: multi-config Slurm runner for scenario ensembles run_slurm.sh took a single hard-coded configfile and passed -A {cluster.account}, which Sherlock rejects (it has no accounts). It now takes a list of overlay configs, runs each as its own `snakemake --cluster` invocation with bounded concurrency, and reads each build target from config/weather_years/manifest.tsv. slurm_submit.sh wraps sbatch because several rules compute mem_mb as a float (e.g. (input.size // 150000) * attempt * 1.5) and `sbatch --mem 4500.0` is rejected; it rounds up and floors at 2000 MB. Note: the target must precede --configfile, which is nargs='+' and would otherwise swallow it as a third config file. Co-Authored-By: Claude Opus 5 --- workflow/run_slurm.sh | 99 ++++++++++++++++++++++++++++++++++++++-- workflow/slurm_submit.sh | 34 ++++++++++++++ 2 files changed, 130 insertions(+), 3 deletions(-) create mode 100644 workflow/slurm_submit.sh diff --git a/workflow/run_slurm.sh b/workflow/run_slurm.sh index 56bdd4776..b7a1e595b 100644 --- a/workflow/run_slurm.sh +++ b/workflow/run_slurm.sh @@ -1,3 +1,96 @@ -# SLURM specifications live in config/config.slurm.yaml & the individual rules -# GRB_LICENSE_FILE=/share/software/user/restricted/gurobi/11.0.2/licenses/gurobi.lic⁠ -snakemake --cluster "sbatch -A {cluster.account} --mail-type ALL --mail-user {cluster.email} -p {cluster.partition} -o {cluster.output} -e {cluster.error} -c {threads} --mem {resources.mem_mb} --time {resources.walltime}" --cluster-config config/config.slurm.yaml --jobs 20 --latency-wait 60 --rerun-incomplete --configfile config/CH1/config.tamu.single_horizon.bau.yaml +#!/usr/bin/env bash +# Run one or more PyPSA-USA scenario overlays, fanning each rule out to Slurm. +# +# bash run_slurm.sh ca2040_wy2019_z4.yaml ca2040_wy2019_county.yaml +# bash run_slurm.sh $(cd config/weather_years && ls ca2040_wy*_z4.yaml) +# +# Each overlay becomes its own `snakemake --cluster` invocation. Overlays run +# CONCURRENCY-at-a-time; within each, up to JOBS rule-jobs are in flight. +# +# Config layering (left to right, later wins): +# repo_data/config/config.{slurm,common,plotting,api,sector,default}.yaml auto-loaded by the Snakefile +# repo_data/config/config.california.yaml passed here +# config/weather_years/.yaml passed here +# +# The build target comes from config/weather_years/manifest.tsv and is the +# add_extra_components output, so ll/opts/sector never enter the DAG. +# +# Sherlock notes: no --account (there are no accounts); this script itself must +# run inside a job or an sh_dev shell, never on the login node. +set -uo pipefail + +cd "$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" + +PARTITION="${PARTITION:-serc}" +EMAIL="${EMAIL:-ctehran@stanford.edu}" +JOBS="${JOBS:-20}" # concurrent Slurm rule-jobs per overlay +CONCURRENCY="${CONCURRENCY:-4}" # concurrent snakemake drivers +RESTART_TIMES="${RESTART_TIMES:-0}" # 0 while benchmarking: a retry re-runs with + # attempt*mem and would pollute MaxRSS data +# Use the venv binary directly rather than `uv run`, which would re-resolve the +# environment on every invocation. See CLAUDE.md for the uv setup. +SNAKEMAKE="${SNAKEMAKE:-$(cd .. && pwd)/.venv/bin/snakemake}" + +CA_CONFIG=repo_data/config/config.california.yaml +OVERLAY_DIR=config/weather_years +MANIFEST="$OVERLAY_DIR/manifest.tsv" + +[ $# -ge 1 ] || { echo "usage: $0 [overlay.yaml ...]" >&2; exit 2; } +[ -f "$MANIFEST" ] || { echo "missing $MANIFEST -- run $OVERLAY_DIR/generate_overlays.sh" >&2; exit 2; } + +mkdir -p logs/slurm logs/drivers + +export PARTITION EMAIL +# slurm_submit.sh rounds the float mem_mb some rules produce; see that file. +SBATCH_CMD="bash slurm_submit.sh {rule} {threads} {resources.mem_mb} {resources.walltime}" + +run_one() { + local overlay="$1" row run_name target + row=$(awk -F'\t' -v o="$overlay" '$1==o{print; exit}' "$MANIFEST") + if [ -z "$row" ]; then + echo "[SKIP] $overlay not in $MANIFEST" >&2 + return 1 + fi + run_name=$(printf '%s' "$row" | cut -f2) + target=$(printf '%s' "$row" | cut -f3) + + echo "[START] $run_name -> $target" + # shellcheck disable=SC2086 + # NB: the target MUST precede --configfile. `--configfile` is nargs='+' and + # would otherwise swallow the target as a third config file. + $SNAKEMAKE \ + "$target" \ + --cluster "$SBATCH_CMD" \ + --configfile "$CA_CONFIG" "$OVERLAY_DIR/$overlay" \ + --jobs "$JOBS" \ + --latency-wait 60 \ + --rerun-incomplete \ + --restart-times "$RESTART_TIMES" \ + --default-resources "mem_mb=8000" "walltime='02:00:00'" \ + --printshellcmds \ + > "logs/drivers/${run_name}.log" 2>&1 + local rc=$? + if [ $rc -eq 0 ]; then + echo "[DONE ] $run_name" + else + echo "[FAIL ] $run_name (rc=$rc) -- see logs/drivers/${run_name}.log" >&2 + fi + return $rc +} + +fails=0 +for overlay in "$@"; do + while [ "$(jobs -rp | wc -l)" -ge "$CONCURRENCY" ]; do wait -n; done + run_one "$overlay" & +done +wait + +for overlay in "$@"; do + row=$(awk -F'\t' -v o="$overlay" '$1==o{print; exit}' "$MANIFEST") + [ -n "$row" ] || continue + t=$(printf '%s' "$row" | cut -f3) + [ -f "$t" ] || { echo "[MISSING] $t"; fails=$((fails + 1)); } +done + +echo "=== ${#@} overlays attempted, ${fails} missing target(s) ===" +exit $(( fails > 0 ? 1 : 0 )) diff --git a/workflow/slurm_submit.sh b/workflow/slurm_submit.sh new file mode 100644 index 000000000..f49950b0d --- /dev/null +++ b/workflow/slurm_submit.sh @@ -0,0 +1,34 @@ +#!/usr/bin/env bash +# Snakemake --cluster submit wrapper. +# +# slurm_submit.sh +# +# Exists because several rules compute mem_mb as a float +# (e.g. `(input.size // 150000) * attempt * 1.5` in build_electricity.smk), and +# `sbatch --mem 4500.0` is rejected. This rounds up to an integer MB and floors +# it at a usable minimum. Prints the job id (--parsable) for snakemake to track. +set -euo pipefail + +rule="$1"; threads="$2"; mem_raw="$3"; walltime="$4"; jobscript="$5" + +# Round up to whole MB; floor at 2000 MB so a tiny input.size can't request ~0. +mem_mb=$(awk -v m="$mem_raw" 'BEGIN{ v=int(m); if (v Date: Mon, 31 Aug 2026 19:20:55 -0700 Subject: [PATCH 03/13] fix: point Python at certifi so HTTPS retrieve rules can verify certs The uv-managed CPython ships no CA bundle -- ssl.get_default_verify_paths() returns cafile=None -- so retrieve_cpuc_servm_load died with CERTIFICATE_VERIFY_FAILED. A curl reachability test does not catch this, because curl reads the system trust store that Python never consults. run_slurm.sh now derives SSL_CERT_FILE/REQUESTS_CA_BUNDLE from certifi and relies on sbatch --export=ALL to propagate them to every rule job. Also adds probe_2019.sbatch, the driver job for the resource-tuning probe. Co-Authored-By: Claude Opus 5 --- workflow/probe_2019.sbatch | 33 +++++++++++++++++++++++++++++++++ workflow/run_slurm.sh | 11 +++++++++++ 2 files changed, 44 insertions(+) create mode 100644 workflow/probe_2019.sbatch diff --git a/workflow/probe_2019.sbatch b/workflow/probe_2019.sbatch new file mode 100644 index 000000000..00bc580c2 --- /dev/null +++ b/workflow/probe_2019.sbatch @@ -0,0 +1,33 @@ +#!/bin/bash +# Probe run: weather year 2019, both resolutions, purely to collect per-rule +# resource data before launching the remaining 45 overlays. +# +# sbatch probe_2019.sbatch +# +# This is the DRIVER job. It holds two `snakemake --cluster` processes, which +# submit each rule as its own Slurm job. Its own CPU/memory needs are tiny; its +# walltime must cover the whole build span including queue wait for rule jobs. +#SBATCH -p serc +#SBATCH -J ca2040_probe +#SBATCH --time=06:00:00 +#SBATCH --cpus-per-task=2 +#SBATCH --mem=8G +#SBATCH -o logs/drivers/probe-%j.out +#SBATCH -e logs/drivers/probe-%j.err +#SBATCH --mail-type=FAIL +#SBATCH --mail-user=ctehran@stanford.edu + +set -uo pipefail +cd /oak/stanford/groups/iazevedo/kamran/CH3/pypsa-usa/workflow + +export PARTITION=serc +export EMAIL=ctehran@stanford.edu +export JOBS=20 # concurrent Slurm rule-jobs per overlay +export CONCURRENCY=2 # both overlays at once +export RESTART_TIMES=0 # a retry would re-run with attempt*mem and skew MaxRSS + +echo "=== probe start $(date -Is) on $(hostname) ===" +bash run_slurm.sh ca2040_wy2019_z4.yaml ca2040_wy2019_county.yaml +rc=$? +echo "=== probe end $(date -Is) rc=$rc ===" +exit $rc diff --git a/workflow/run_slurm.sh b/workflow/run_slurm.sh index b7a1e595b..941d20aab 100644 --- a/workflow/run_slurm.sh +++ b/workflow/run_slurm.sh @@ -40,6 +40,17 @@ MANIFEST="$OVERLAY_DIR/manifest.tsv" mkdir -p logs/slurm logs/drivers +# The uv-managed CPython ships no CA bundle -- ssl.get_default_verify_paths() +# returns cafile=None -- so every retrieve_* rule that downloads over HTTPS dies +# with CERTIFICATE_VERIFY_FAILED. (curl works, because it reads the system +# store, which is why a curl reachability test does not catch this.) Point +# Python at certifi. sbatch --export=ALL propagates these to the rule jobs. +if [ -z "${SSL_CERT_FILE:-}" ]; then + _certifi=$("$(dirname "$SNAKEMAKE")/python" -c 'import certifi; print(certifi.where())' 2>/dev/null || true) + [ -n "$_certifi" ] && export SSL_CERT_FILE="$_certifi" REQUESTS_CA_BUNDLE="$_certifi" +fi +[ -n "${SSL_CERT_FILE:-}" ] || echo "WARNING: SSL_CERT_FILE unset; HTTPS retrieve rules will likely fail" >&2 + export PARTITION EMAIL # slurm_submit.sh rounds the float mem_mb some rules produce; see that file. SBATCH_CMD="bash slurm_submit.sh {rule} {threads} {resources.mem_mb} {resources.walltime}" From 8b64feb70644f4041539b0acca1b78b1e8dab681 Mon Sep 17 00:00:00 2001 From: Kamran Date: Mon, 31 Aug 2026 19:26:31 -0700 Subject: [PATCH 04/13] tooling: per-rule Slurm resource collector for tuning sacct splits what we need across two rows: the allocation row has JobName and ReqMem, the .batch row has MaxRSS. `sacct -X` shows MaxRSS blank, which is a quiet trap. This joins them on the base job id and aggregates per rule. Co-Authored-By: Claude Opus 5 --- workflow/collect_benchmarks.sh | 49 ++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) create mode 100644 workflow/collect_benchmarks.sh diff --git a/workflow/collect_benchmarks.sh b/workflow/collect_benchmarks.sh new file mode 100644 index 000000000..cbb2039f2 --- /dev/null +++ b/workflow/collect_benchmarks.sh @@ -0,0 +1,49 @@ +#!/usr/bin/env bash +# Per-rule resource report for tuning Slurm requests. +# +# bash collect_benchmarks.sh [SINCE] # SINCE defaults to now-6hours +# +# sacct splits the data we need across two rows per job: the allocation row +# carries JobName/ReqMem, the .batch row carries MaxRSS. This joins them on the +# base job id and aggregates per rule, so `sacct -X` (which shows MaxRSS blank) +# is never the thing you look at. +set -uo pipefail + +SINCE="${1:-now-6hours}" + +sacct -S "$SINCE" -P --format=JobID,JobName,State,Elapsed,MaxRSS,ReqMem 2>/dev/null | +awk -F'|' ' + NR == 1 { next } + { + split($1, a, ".") + id = a[1] + if ($1 !~ /\./) { name[id] = $2; state[id] = $3; elapsed[id] = $4; req[id] = $6 } + else if ($1 ~ /\.batch$/) { gsub(/K$/, "", $5); rss[id] = $5 + 0 } + } + END { + for (id in name) { + n = name[id] + if (n == "batch" || n == "extern" || n == "" ) continue + if (n ~ /^(bash|uv|python|snakemake|ca2040_probe)$/) continue + cnt[n]++ + sum_rss[n] += rss[id] + if (rss[id] > max_rss[n]) max_rss[n] = rss[id] + split(elapsed[id], t, ":") + secs = t[1] * 3600 + t[2] * 60 + t[3] + if (secs > max_secs[n]) max_secs[n] = secs + reqm[n] = req[id] + if (state[id] !~ /COMPLETED/) bad[n] = bad[n] " " state[id] + } + printf "%-34s %5s %11s %11s %10s %8s %s\n", "RULE", "N", "MAXRSS_MB", "PEAK_ELAPSED", "REQ_MEM", "RATIO", "NONCOMPLETE" + for (n in cnt) { + mb = max_rss[n] / 1024 + r = reqm[n]; gsub(/[MGn]$/, "", r) + if (reqm[n] ~ /G/) r = r * 1024 + ratio = (mb > 0 && r > 0) ? sprintf("%.1fx", r / mb) : "-" + printf "%-34s %5d %11.0f %11s %10s %8s %s\n", \ + n, cnt[n], mb, \ + sprintf("%d:%02d:%02d", max_secs[n]/3600, (max_secs[n]%3600)/60, max_secs[n]%60), \ + reqm[n], ratio, bad[n] + } + } +' | { read -r hdr; echo "$hdr"; sort -k3 -nr; } From a128c2e9149714c7ae1a059e5aad292037df5f0e Mon Sep 17 00:00:00 2001 From: Kamran Date: Mon, 31 Aug 2026 22:01:48 -0700 Subject: [PATCH 05/13] chore: campaign scheduler grouping and CA clustering aggregation strategies --- .../repo_data/config/config.california.yaml | 11 +++++ workflow/run_slurm.sh | 48 ++++++++++++++++++- 2 files changed, 57 insertions(+), 2 deletions(-) diff --git a/workflow/repo_data/config/config.california.yaml b/workflow/repo_data/config/config.california.yaml index 05bd6d255..4c0eb1262 100644 --- a/workflow/repo_data/config/config.california.yaml +++ b/workflow/repo_data/config/config.california.yaml @@ -69,6 +69,17 @@ clustering: aggregation_strategies: generators: p_min_pu: 'capacity_weighted_average' # keep the UC feasibility invariant through clustering + # CPUC out-of-state units (REMOTE_PREFIX "R ") are attached to CA buses with + # p_nom_extendable=False, while in-state units of the same carrier are extendable + # because their carrier is in extendable_carriers.Generator. Clustering groups + # one-ports by (bus, carrier), so those two land in one group and PyPSA's default + # `consense` raises. Same situation `committable: any` already handles in the base. + p_nom_extendable: any + # Hydro units disagree on up_time_before ({0, 1}). `max` keeps the integer + # dtype (capacity_weighted_average would yield a float) and is the + # UC-feasibility-safe direction: the cluster inherits the longest online + # history, so it is never spuriously forced to stay committed. + up_time_before: max solving: mem: 30000 diff --git a/workflow/run_slurm.sh b/workflow/run_slurm.sh index 941d20aab..0f09e229b 100644 --- a/workflow/run_slurm.sh +++ b/workflow/run_slurm.sh @@ -76,6 +76,7 @@ run_one() { --jobs "$JOBS" \ --latency-wait 60 \ --rerun-incomplete \ + --rerun-triggers $RERUN_TRIGGERS \ --restart-times "$RESTART_TIMES" \ --default-resources "mem_mb=8000" "walltime='02:00:00'" \ --printshellcmds \ @@ -89,10 +90,53 @@ run_one() { return $rc } -fails=0 +# Scheduling. +# +# Snakemake locks a run's input/output FILE SETS, not the working directory, so +# two runs collide only when both want to CREATE the same file. Almost all +# outputs live under resources//, but two tiers do not: +# +# global resources/powerplants/powerplants.csv, data/cpuc/..., data/caiso/..., +# data/nrel/... -- shared by every overlay +# per-year data/godeeep/historical/*_{year}_*.nc, cutouts/*_{year}.nc +# -- shared by the z4 and county overlays of the SAME weather year +# +# Hence: bootstrap one overlay alone so the global tier exists, then serialise +# within a weather year and parallelise across years. +BOOTSTRAP="${BOOTSTRAP:-1}" + +# Only rebuild on mtime. Without this a code or config edit mid-campaign pulls +# the shared retrieve rules back into every DAG at once and they all collide. +RERUN_TRIGGERS="${RERUN_TRIGGERS:-mtime}" + +year_of() { printf '%s' "$1" | sed -n 's/.*_wy\([0-9]\{4\}\)_.*/\1/p'; } + +declare -A GROUP +YEARS=() for overlay in "$@"; do + y=$(year_of "$overlay") + [ -n "$y" ] || y="ungrouped_$overlay" + [ -n "${GROUP[$y]:-}" ] || YEARS+=("$y") + GROUP[$y]="${GROUP[$y]:-} $overlay" +done + +fails=0 + +if [ "$BOOTSTRAP" = "1" ] && [ $# -gt 1 ]; then + echo "=== bootstrap: building globally-shared artifacts via $1 ===" + run_one "$1" || fails=$((fails + 1)) +fi + +run_group() { + local rc=0 o + for o in "$@"; do run_one "$o" || rc=1; done + return $rc +} + +for y in "${YEARS[@]}"; do while [ "$(jobs -rp | wc -l)" -ge "$CONCURRENCY" ]; do wait -n; done - run_one "$overlay" & + # shellcheck disable=SC2086 + run_group ${GROUP[$y]} & done wait From bcaeac96d108ba85da2eaa15572b2b0540fe8495 Mon Sep 17 00:00:00 2001 From: Kamran Date: Mon, 31 Aug 2026 22:11:08 -0700 Subject: [PATCH 06/13] tooling: benchmark TSV aggregator with per-rule mem/walltime recommendations Complements collect_benchmarks.sh: reads snakemake benchmark TSVs, resolves rule names against rules/*.smk, and prints recommended mem_mb (1.5x peak RSS) and walltime (3x peak runtime) per rule. Offered upstream in issue #808. Co-Authored-By: Claude Opus 5 --- workflow/report_benchmarks.py | 141 ++++++++++++++++++++++++++++++++++ 1 file changed, 141 insertions(+) create mode 100644 workflow/report_benchmarks.py diff --git a/workflow/report_benchmarks.py b/workflow/report_benchmarks.py new file mode 100644 index 000000000..6a3e8f1e9 --- /dev/null +++ b/workflow/report_benchmarks.py @@ -0,0 +1,141 @@ +#!/usr/bin/env python3 +"""Turn snakemake benchmark TSVs into per-rule Slurm resource recommendations. + + ../.venv/bin/python report_benchmarks.py [benchmarks_dir] + +Snakemake writes one TSV per job under benchmarks/, with columns +s / h:m:s / max_rss / max_vms / max_uss / max_pss / io_in / io_out / +mean_load / cpu_time. Memory is in MB, s in seconds. + +The rule name is recovered from the path, which mirrors each rule's `log:` +layout, so it is stripped of run name, interconnect, and wildcard-bearing +prefixes/suffixes. + +Recommendations use headroom on the observed peak: + mem = max(2 GB, ceil(1.5 x peak max_rss / 500 MB) * 500 MB) + time = max(10 min, ceil(3 x peak runtime / 5 min) * 5 min) +The multipliers differ because memory overrun is a hard kill while a generous +walltime only costs queue priority; 3x also absorbs a slower/contended node. +""" + +import re +import sys +from pathlib import Path + +import pandas as pd + +root = Path(sys.argv[1] if len(sys.argv) > 1 else "benchmarks") + +# Wildcard debris to strip from the file stem to recover a rule name. The simpl +# wildcard is a number or a word like "county"/"all", so anchor on the known +# forms -- a bare `_s\w+$` would eat `_shapes` and `_substations`. +SIMPL = r"(?:\d+|county|all)" +STEM_SUBS = [ + (rf"_?elec_s{SIMPL}(?:_c\d+[a-z]?)?", ""), # elec_s75_c4 / elec_scounty + (rf"_s{SIMPL}$", ""), # trailing _s75 / _scounty + (r"_c?\d{2,}[a-z]?$", ""), # trailing cluster counts / horizons + (r"^elec_", ""), +] + + +# Benchmark paths are not uniform: most are benchmarks///, +# but some rules (cluster_network) write benchmarks///. +# So resolve against the real rule list rather than inferring from position. +KNOWN_RULES = set() +for smk in Path("rules").glob("*.smk"): + KNOWN_RULES.update(re.findall(r"^rule\s+(\w+)\s*:", smk.read_text(), re.M)) + + +def rule_from(path: Path) -> str: + # A rule name appearing as a directory component wins outright. + for part in path.parts: + if part in KNOWN_RULES: + return part + stem = path.stem + # Longest known rule that prefixes the stem, so build_renewable_profiles_onwind + # keeps its tech suffix but build_shapes is not truncated. + hits = [r for r in KNOWN_RULES if stem.startswith(r)] + if hits: + best = max(hits, key=len) + return stem if stem.startswith(f"{best}_") and best != stem else best + for pat, rep in STEM_SUBS: + stem = re.sub(pat, rep, stem) + stem = stem.strip("_") + return stem or path.parent.name + + +rows = [] +for f in sorted(root.rglob("*")): + if not f.is_file(): + continue + try: + df = pd.read_csv(f, sep="\t") + except Exception: + continue + if df.empty or "s" not in df.columns: + continue + parts = f.relative_to(root).parts + run = parts[0] if len(parts) > 1 else "-" + r = df.iloc[0] + rows.append( + { + "rule": rule_from(f), + "run": run, + "secs": float(r["s"]), + "max_rss": float(r.get("max_rss", 0) or 0), + "max_pss": float(r.get("max_pss", 0) or 0), + "cpu_time": float(r.get("cpu_time", 0) or 0), + "mean_load": float(r.get("mean_load", 0) or 0), + } + ) + +if not rows: + sys.exit(f"no benchmark files under {root}") + +d = pd.DataFrame(rows) + + +def fmt_hms(secs: float) -> str: + s = int(round(secs)) + return f"{s // 3600:d}:{(s % 3600) // 60:02d}:{s % 60:02d}" + + +def rec_mem(peak_mb: float) -> int: + import math + + return max(2000, int(math.ceil(1.5 * peak_mb / 500.0) * 500)) + + +def rec_time(peak_s: float) -> str: + import math + + mins = max(10, int(math.ceil(3 * peak_s / 60.0 / 5.0) * 5)) + return f"{mins // 60:02d}:{mins % 60:02d}:00" + + +agg = ( + d.groupby("rule") + .agg( + n=("secs", "size"), + peak_s=("secs", "max"), + peak_rss=("max_rss", "max"), + mean_load=("mean_load", "max"), + ) + .sort_values("peak_rss", ascending=False) +) + +agg["rec_mem_mb"] = agg.peak_rss.map(rec_mem) +agg["rec_walltime"] = agg.peak_s.map(rec_time) +agg["peak_time"] = agg.peak_s.map(fmt_hms) + +print(f"{'RULE':<44} {'N':>2} {'PEAK_RSS_MB':>11} {'PEAK_TIME':>10} {'LOAD%':>6} {'REC_MEM_MB':>10} {'REC_WALLTIME':>12}") +print("-" * 108) +for name, r in agg.iterrows(): + print( + f"{name:<44} {int(r.n):>2} {r.peak_rss:>11.0f} {r.peak_time:>10} " + f"{r.mean_load:>6.0f} {int(r.rec_mem_mb):>10} {r.rec_walltime:>12}" + ) + +print() +print(f"total jobs benchmarked: {len(d)} distinct rules: {len(agg)}") +print(f"sum of peak runtimes (serial lower bound): {fmt_hms(agg.peak_s.sum())}") From 53a46759443190c72f481e7803c75d5faffb5155 Mon Sep 17 00:00:00 2001 From: Kamran Date: Mon, 31 Aug 2026 22:30:03 -0700 Subject: [PATCH 07/13] feat: external-region representations for out-of-footprint supply MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds `electricity.imports.representation` (store | generator) and moves all import/export construction into a new self-contained `external_regions` module. * `store` (default) is the existing behaviour, byte-for-byte: a bottomless Store per external zone and PRICED interface links. * `generator` puts a priced, non-extendable import generator (carrier `unspecified_imports`, carrying `imports.co2_emissions`) on the external bus, leaves the interface link unpriced, and moves the CPUC contracted out-of-state units BEHIND the boundary. Their deliveries then traverse a carrier-`imports` link, so `opts/interchange.py` meters them against the import volume cap and the interface capacity bounds them — which is the point. `add_electricity` now splits the contracted-unit derivation from its attachment: in `generator` mode it serializes the fully-derived bundle (unit table, borrowed VRE profiles, cost table) to a new `remote_units_s{simpl}.pkl` output that `add_extra_components` attaches at the external buses. Each unit is placed behind the boundary zone whose state matches its plants, preferring the zone with the largest inbound interface capacity, falling back (with a warning) to the largest boundary zone overall. Also deletes the dead `trim_network` path: the schema typed `model_topology.trim` as a boolean while the code indexed `trim_topology["zone"]`, so it could never run. The California config switches to `generator` and drops the now-obsolete `p_nom_extendable: any` clustering workaround — remote units no longer exist inside the footprint at clustering time. Co-Authored-By: Claude Opus 5 --- docs/source/configtables/electricity.csv | 1 + .../repo_data/config/config.california.yaml | 7 +- workflow/repo_data/config/config.default.yaml | 1 + workflow/rules/build_electricity.smk | 11 +- workflow/schemas/config.schema.yaml | 2 +- workflow/scripts/add_electricity.py | 198 ++++- workflow/scripts/add_extra_components.py | 442 +---------- workflow/scripts/external_regions.py | 691 ++++++++++++++++++ workflow/scripts/opts/interfaces.py | 2 +- workflow/scripts/plot_statistics.py | 4 +- 10 files changed, 923 insertions(+), 436 deletions(-) create mode 100644 workflow/scripts/external_regions.py diff --git a/docs/source/configtables/electricity.csv b/docs/source/configtables/electricity.csv index 7741f971d..20b79ec66 100644 --- a/docs/source/configtables/electricity.csv +++ b/docs/source/configtables/electricity.csv @@ -52,6 +52,7 @@ demand_response:,,,Settings to activate and configure demand response ,,, imports:,,,Configure electric imports from regions outside of model scope -- enable,,``true`` or ``false``,Enable electric imports +-- representation,,``store`` or ``generator``,"How out-of-footprint supply is represented. ``store`` (default) puts a bottomless ``Store`` behind each external zone and PRICES the interface ``Link``. ``generator`` puts a priced import ``Generator`` (carrier ``unspecified_imports``, carrying ``co2_emissions``) on the external bus, leaves the interface link unpriced, and moves the CPUC contracted out-of-state units behind the boundary so their deliveries are metered against ``volume_limit`` and the interface capacity." -- costs,$/MWh,``wholesale`` or ```` or ``float``,Cost of electric imports from regions outside of model scope. ``Wholesale`` will use in monthly wholesales electric prices. ``float`` will assign a user specified value. ``carrier`` will take average marginal cost of the carrier. -- co2_emissions,CO2/MWh,``float``,CO2 emissions of electric imports from regions outside of model scope. -- capacity_limit,,``true`` or ``false``,Enable capacity limit for electric imports from regions outside of model scope diff --git a/workflow/repo_data/config/config.california.yaml b/workflow/repo_data/config/config.california.yaml index 4c0eb1262..f3faa1d60 100644 --- a/workflow/repo_data/config/config.california.yaml +++ b/workflow/repo_data/config/config.california.yaml @@ -55,6 +55,7 @@ electricity: profile: servm imports: enable: true + representation: generator # CPUC out-of-state contracts sit BEHIND the WECC boundary, so their deliveries are metered against the import cap volume_limit: 25 # % of annual CA load balancing_period: year exports: @@ -69,12 +70,6 @@ clustering: aggregation_strategies: generators: p_min_pu: 'capacity_weighted_average' # keep the UC feasibility invariant through clustering - # CPUC out-of-state units (REMOTE_PREFIX "R ") are attached to CA buses with - # p_nom_extendable=False, while in-state units of the same carrier are extendable - # because their carrier is in extendable_carriers.Generator. Clustering groups - # one-ports by (bus, carrier), so those two land in one group and PyPSA's default - # `consense` raises. Same situation `committable: any` already handles in the base. - p_nom_extendable: any # Hydro units disagree on up_time_before ({0, 1}). `max` keeps the integer # dtype (capacity_weighted_average would yield a float) and is the # UC-feasibility-safe direction: the cluster inherits the longest online diff --git a/workflow/repo_data/config/config.default.yaml b/workflow/repo_data/config/config.default.yaml index c23dd3123..5214d8aaa 100644 --- a/workflow/repo_data/config/config.default.yaml +++ b/workflow/repo_data/config/config.default.yaml @@ -228,6 +228,7 @@ electricity: # neighbors using historical EIA interchange data. imports: enable: false + representation: store # store = priced links off a bottomless Store; generator = priced generator (+ CPUC contracted units) behind an unpriced interface link costs: wholesale # wholesale | carrier | float — how imported energy is priced co2_emissions: 0.428 # tCO2/MWh assigned to imported energy capacity_limit: true # cap import power to historical maxima diff --git a/workflow/rules/build_electricity.smk b/workflow/rules/build_electricity.smk index cf5422079..ebda8c080 100644 --- a/workflow/rules/build_electricity.smk +++ b/workflow/rules/build_electricity.smk @@ -804,6 +804,9 @@ rule add_electricity: planning_horizons=config["scenario"]["planning_horizons"], eia_api=config["api"]["eia"], remote_contracted=config["electricity"].get("remote_contracted_resources", {}), + imports_representation=config_provider( + "electricity", "imports", "representation", default="store" + ), input: unpack(dynamic_fuel_price_files), unpack(remote_contracted_resource_files), @@ -870,7 +873,11 @@ rule add_electricity: else [] ), output: - NETWORKS + "{interconnect}/elec_s{simpl}_l_pp.pkl", + network=NETWORKS + "{interconnect}/elec_s{simpl}_l_pp.pkl", + # CPUC contracted out-of-state units, serialized for attachment behind + # the model boundary when electricity.imports.representation == generator. + # Always declared; a sentinel is written otherwise. + remote_units=NETWORKS + "{interconnect}/remote_units_s{simpl}.pkl", log: LOGS + "{interconnect}/elec_s{simpl}_add_electricity.log", benchmark: @@ -1064,6 +1071,7 @@ rule add_extra_components: + "{interconnect}/regions_onshore_s{simpl}_{clusters}.geojson", flowgates=flowgates_for_extra_components, reeds_memberships="repo_data/ReEDS_Constraints/membership.csv", + remote_units=NETWORKS + "{interconnect}/remote_units_s{simpl}.pkl", co2_storage=( CO2 + "{interconnect}/co2_storage_s{simpl}_{clusters}.csv" if config["scenario"]["sector"] == "" and config["co2"]["storage"] is True @@ -1073,7 +1081,6 @@ rule add_extra_components: params: retirement=config["electricity"].get("retirement", "technical"), demand_response=config["electricity"].get("demand_response", {}), - trim_network=config_provider("model_topology", "trim", default=False), imports=config_provider("electricity", "imports", default={}), exports=config_provider("electricity", "exports", default={}), weather_year=config_provider("renewable_weather_years"), diff --git a/workflow/schemas/config.schema.yaml b/workflow/schemas/config.schema.yaml index faabe09af..0368f2221 100644 --- a/workflow/schemas/config.schema.yaml +++ b/workflow/schemas/config.schema.yaml @@ -90,7 +90,6 @@ properties: transmission_network: {enum: [reeds, tamu]} topological_boundaries: {enum: [county, reeds_zone, state]} interface_transmission_limits: {type: boolean} - trim: {type: boolean} include: description: > Zone subset for footprint-scoped runs. null means "no filter"; see the @@ -285,6 +284,7 @@ properties: additionalProperties: false properties: enable: {type: boolean} + representation: {enum: [store, generator]} costs: {enum: [wholesale, carrier, float]} co2_emissions: {type: number} capacity_limit: {type: boolean} diff --git a/workflow/scripts/add_electricity.py b/workflow/scripts/add_electricity.py index c0c45b88f..e57947f57 100755 --- a/workflow/scripts/add_electricity.py +++ b/workflow/scripts/add_electricity.py @@ -1301,6 +1301,40 @@ def attach_remote_contracted_resources( One row per ledger row with its disposition (``status``, ``carrier``, ``bus``, ``p_nom``), for logging and testing. """ + unit_df, summary = build_remote_contracted_units(n, plants_prefilter, remote_df, weights, costs, tech_map) + if unit_df.empty: + return summary + + dropped_vre = attach_remote_units(n, unit_df, costs, conventional_carriers, unit_commitment) + return _finalize_remote_summary(summary, dropped_vre) + + +def build_remote_contracted_units( + n: pypsa.Network, + plants_prefilter: pd.DataFrame, + remote_df: pd.DataFrame, + weights: pd.DataFrame, + costs: pd.DataFrame, + tech_map: pd.DataFrame | None = None, +) -> tuple[pd.DataFrame, pd.DataFrame]: + """Derive the CPUC contracted-unit table without touching the network. + + This is the half of :func:`attach_remote_contracted_resources` that only + needs data: it resolves each ledger row onto its live ``powerplants.csv`` + rows, derives capacity-weighted techno-economics, and picks the California + attachment bus (the region's max-LAF bus). In ``imports.representation: + generator`` the units are attached at an EXTERNAL bus instead, but the + California bus is still resolved here because remote VRE units borrow that + bus's capacity-factor profile. + + Returns + ------- + (unit_df, summary) + ``unit_df`` is indexed by ``REMOTE_PREFIX + cpuc_unit_name`` and carries + every attribute the ``n.add`` calls need, plus the physical ``state`` of + the constituent plants (used to place the unit behind the right external + zone). ``summary`` is the per-ledger-row disposition table. + """ region_bus = _servm_region_buses(weights, n) carrier_to_category, category_to_carrier = _carrier_category_maps(tech_map) @@ -1349,6 +1383,7 @@ def attach_remote_contracted_resources( bus = region_bus[row["servm_region"]] p_nom = float(min(record["capmax_mw"], constituents["p_nom"].sum())) attrs = _remote_unit_attributes(constituents, p_nom, carrier, costs) + state = _remote_unit_state(constituents) record.update( { @@ -1359,7 +1394,9 @@ def attach_remote_contracted_resources( }, ) records.append(record) - units.append({"name": REMOTE_PREFIX + name, "carrier": carrier, "bus": bus, "p_nom": p_nom, **attrs}) + units.append( + {"name": REMOTE_PREFIX + name, "carrier": carrier, "bus": bus, "p_nom": p_nom, "state": state, **attrs}, + ) logger.info( f"Remote contracted unit '{name}' -> bus '{bus}' ({row['servm_region']} max-LAF bus), carrier " @@ -1369,7 +1406,7 @@ def attach_remote_contracted_resources( summary = pd.DataFrame(records) if not units: logger.warning("Remote contracted resources enabled but no ledger row could be attached.") - return summary + return pd.DataFrame(), summary unit_df = pd.DataFrame(units).set_index("name") if unit_df.index.has_duplicates: @@ -1377,6 +1414,40 @@ def attach_remote_contracted_resources( f"Duplicate cpuc_unit_name(s) in the remote contracted-resource file: " f"{sorted(unit_df.index[unit_df.index.duplicated()])}", ) + return unit_df, summary + + +def _remote_unit_state(constituents: pd.DataFrame) -> str: + """Physical state of a contract's constituent plants (dominant by capacity). + + Hoover is the reason this is capacity-weighted rather than a single lookup: + the contract spans EIA 154 (NV) and 8902 (AZ), two halves of the same dam. + """ + if "state" not in constituents.columns: + return "" + states = constituents.dropna(subset=["state"]) + if states.empty: + return "" + return str(states.groupby("state")["p_nom"].sum().idxmax()) + + +def attach_remote_units( + n: pypsa.Network, + unit_df: pd.DataFrame, + costs: pd.DataFrame, + conventional_carriers: list, + unit_commitment: bool, + profiles: pd.DataFrame | None = None, +) -> list[str]: + """Add a derived contracted-unit table to the network at ``unit_df['bus']``. + + Shared by both import representations: in ``store`` mode the buses are + California buses and the profiles are derived here; in ``generator`` mode + ``external_regions`` rewrites ``bus`` to an external bus and supplies the + profiles that were borrowed at the pre-clustering stage. + + Returns the ledger names dropped for want of a VRE profile. + """ add_missing_carriers(n, sorted(set(unit_df["carrier"]))) batteries = unit_df[unit_df["carrier"] == "battery"] @@ -1384,9 +1455,48 @@ def attach_remote_contracted_resources( firm = unit_df.drop(index=batteries.index.union(vre.index)) _attach_remote_firm(n, firm, costs, conventional_carriers, unit_commitment) - dropped_vre = _attach_remote_vre(n, vre, costs) + dropped_vre = _attach_remote_vre(n, vre, costs, profiles=profiles) _attach_remote_batteries(n, batteries, costs) + return dropped_vre + +def build_remote_unit_bundle( + n: pypsa.Network, + unit_df: pd.DataFrame, + summary: pd.DataFrame, + costs: pd.DataFrame, + conventional_carriers: list, + unit_commitment: bool, +) -> dict: + """Serialize the contracted units instead of attaching them. + + Used by ``imports.representation: generator``, where the units belong at an + external bus that only exists later, in ``add_extra_components``. Everything + that stage cannot recompute is carried across: the derived unit table, the + CF profiles borrowed from the (pre-clustering) California buses, and the + cost table the ``n.add`` calls resolve ``capital_cost``/``lifetime`` from. + """ + vre = unit_df[unit_df["carrier"].isin(VRE_PROFILE_CARRIERS)] + profiles, dropped_vre = collect_remote_vre_profiles(n, vre) + unit_df = unit_df.drop(index=[REMOTE_PREFIX + name for name in dropped_vre], errors="ignore") + + bundle = { + "units": unit_df, + "vre_profiles": pd.DataFrame(profiles) if profiles else pd.DataFrame(index=n.snapshots), + "costs": costs, + "conventional_carriers": list(conventional_carriers), + "unit_commitment": bool(unit_commitment), + "summary": _finalize_remote_summary(summary, dropped_vre), + } + logger.info( + f"Serialized {len(unit_df)} remote contracted unit(s) ({unit_df['p_nom'].sum():.1f} MW) for attachment " + "behind the model boundary (imports.representation: generator).", + ) + return bundle + + +def _finalize_remote_summary(summary: pd.DataFrame, dropped_vre: list[str]) -> pd.DataFrame: + """Fold the dropped-VRE outcome into the summary and log the totals.""" if dropped_vre: summary.loc[summary["cpuc_unit_name"].isin(dropped_vre), "status"] = "skipped_no_profile" summary.loc[summary["cpuc_unit_name"].isin(dropped_vre), "p_nom"] = 0.0 @@ -1470,16 +1580,14 @@ def _attach_remote_firm( _apply_remote_seasonal_derates(n, firm[["summer_derate", "winter_derate"]]) -def _attach_remote_vre(n: pypsa.Network, vre: pd.DataFrame, costs: pd.DataFrame) -> list[str]: - """Attach remote wind/solar contracts, borrowing the attachment bus's CF profile. +def collect_remote_vre_profiles(n: pypsa.Network, vre: pd.DataFrame) -> tuple[dict[str, pd.Series], list[str]]: + """Borrow a CF profile per remote VRE unit from its attachment bus. - Returns the ledger names that had to be dropped for want of any profile. + Returns ``(profiles, dropped)``; ``dropped`` holds the ledger names for which + the network carries no profile of that carrier at all. """ + profiles: dict[str, pd.Series] = {} dropped: list[str] = [] - if vre.empty: - return dropped - - profiles = {} for name, unit in vre.iterrows(): profile = _remote_vre_profile(n, unit["bus"], unit["carrier"]) if profile is None: @@ -1490,12 +1598,39 @@ def _attach_remote_vre(n: pypsa.Network, vre: pd.DataFrame, costs: pd.DataFrame) dropped.append(name[len(REMOTE_PREFIX) :]) continue profiles[name] = profile + return profiles, dropped + + +def _attach_remote_vre( + n: pypsa.Network, + vre: pd.DataFrame, + costs: pd.DataFrame, + profiles: pd.DataFrame | None = None, +) -> list[str]: + """Attach remote wind/solar contracts, borrowing the attachment bus's CF profile. + + ``profiles`` short-circuits the borrowing: in ``imports.representation: + generator`` the units land on an external bus that has no profile of its + own, so the profiles borrowed before clustering are carried in on the + serialized bundle instead. + + Returns the ledger names that had to be dropped for want of any profile. + """ + dropped: list[str] = [] + if vre.empty: + return dropped + + if profiles is None: + profiles, dropped = collect_remote_vre_profiles(n, vre) + else: + profiles = {name: profiles[name] for name in vre.index if name in profiles.columns} vre = vre.loc[list(profiles)] if vre.empty: return dropped p_max_pu = pd.DataFrame(profiles).reindex(columns=vre.index) + p_max_pu = p_max_pu.set_axis(n.snapshots) n.add( "Generator", vre.index, @@ -1825,6 +1960,7 @@ def main(snakemake): # (small) slice of the fleet the CPUC ledger points at BEFORE that filter runs; # attach_remote_contracted_resources adds them back at CA buses later. remote_contracted = dict(getattr(params, "remote_contracted", None) or {}) + imports_representation = getattr(params, "imports_representation", "store") or "store" remote_df = None plants_prefilter = None if remote_contracted.get("enable", False): @@ -1924,18 +2060,42 @@ def main(snakemake): # Runs last among the attach_* steps: remote VRE contracts copy the # capacity-factor profile of their California attachment bus, which only # exists once attach_wind_and_solar has run. + # + # `imports.representation: generator` moves these units BEHIND the model + # boundary, so they must not be added here — add_extra_components attaches + # them at the external import buses instead. The bundle output is declared + # unconditionally by the rule, so a sentinel is written when there is + # nothing to hand over. + remote_bundle = None if remote_df is not None: - attach_remote_contracted_resources( + unit_df, summary = build_remote_contracted_units( n, plants_prefilter, remote_df, pd.read_csv(snakemake.input["servm_load_weights"]), costs, - conventional_carriers, - extendable_carriers, tech_map=pd.read_csv(snakemake.input["servm_tech_map"]), - unit_commitment=params.conventional["unit_commitment"], ) + if unit_df.empty: + pass + elif imports_representation == "generator": + remote_bundle = build_remote_unit_bundle( + n, + unit_df, + summary, + costs, + conventional_carriers, + params.conventional["unit_commitment"], + ) + else: + dropped_vre = attach_remote_units( + n, + unit_df, + costs, + conventional_carriers, + params.conventional["unit_commitment"], + ) + _finalize_remote_summary(summary, dropped_vre) update_p_nom_max(n) @@ -2010,7 +2170,7 @@ def main(snakemake): axis=1, ) - output_folder = os.path.dirname(snakemake.output[0]) + "/base_network" + output_folder = os.path.dirname(snakemake.output.network) + "/base_network" export_network_for_gis_mapping(n, output_folder) clean_bus_data(n) @@ -2018,8 +2178,12 @@ def main(snakemake): n.meta = snakemake.config log_network_schema(n, stage="exit", baseline=schema_entry) - # n.export_to_netcdf(snakemake.output[0]) - pickle.dump(n, open(snakemake.output[0], "wb")) + # n.export_to_netcdf(snakemake.output.network) + pickle.dump(n, open(snakemake.output.network, "wb")) + + # Always written, even when it holds nothing: snakemake requires every + # declared output to exist. + pickle.dump(remote_bundle, open(snakemake.output.remote_units, "wb")) if __name__ == "__main__": diff --git a/workflow/scripts/add_extra_components.py b/workflow/scripts/add_extra_components.py index 061ddc93c..b167d5f2f 100644 --- a/workflow/scripts/add_extra_components.py +++ b/workflow/scripts/add_extra_components.py @@ -9,8 +9,13 @@ from _helpers import calculate_annuity, configure_logging, load_costs, log_network_schema from add_electricity import add_missing_carriers from constants import HOURS_PER_YEAR -from eia import FuelCosts -from opts._helpers import get_region_buses +from external_regions import ( + add_external_regions, + convert_flowgates_to_state, + format_flowgates_for_imports_exports, + load_remote_unit_bundle, + resolve_trade_costs, +) from shapely.geometry import Point idx = pd.IndexSlice @@ -696,370 +701,6 @@ def add_demand_response( ) -def trim_network(n, trim_topology): - """ - Trim_network splits the network into two parts: - - The internal network, which is the network within the specified zones. - - The external network, which is the network outside the specified zones. - - The internal network is retained and unchanged. While the external network components are removed. The external buses which are directly connected to the internal network are aggregated to the `nerc_reg` value of their buses. - The only generators kept are the OCGTs at the external buses, which are set to non-extendable. - - The external OCGT generators are set to the carrier name `imports` and retain the same emissions intensity. - - """ - retain_zones = trim_topology["zone"] - internal_buses = get_region_buses(n, retain_zones) - if internal_buses.empty: - logger.warning("No internal buses found, skipping trim_network") - return None - - # Get all lines and links connected to internal buses - retain_lines = n.lines[n.lines.bus0.isin(internal_buses.index) | n.lines.bus1.isin(internal_buses.index)] - retain_links = n.links[n.links.bus0.isin(internal_buses.index) | n.links.bus1.isin(internal_buses.index)] - - # Find buses to remove (those not connected to internal network) - buses_to_remove = n.buses[ - ~n.buses.index.isin(retain_lines.bus0) - & ~n.buses.index.isin(retain_lines.bus1) - & ~n.buses.index.isin(retain_links.bus0) - & ~n.buses.index.isin(retain_links.bus1) - ] - - # Find external buses to keep (connected to internal network but not internal) - external_buses_to_keep = n.buses.loc[ - ~n.buses.index.isin(buses_to_remove.index) & ~n.buses.index.isin(internal_buses.index) - ] - - # Remove components at buses that are being removed - for c in n.one_port_components: - component = n.components[c].static - rm = component[component.bus.isin(buses_to_remove.index)] - if not rm.empty: - n.remove(c, rm.index) - - # Remove lines and links at buses being removed - for c in ["Line", "Link"]: - component = n.components[c].static - rm = component[~component.bus0.isin(internal_buses.index) & ~component.bus1.isin(internal_buses.index)] - if not rm.empty: - n.remove(c, rm.index) - - # Remove the buses - n.remove("Bus", buses_to_remove.index) - - # Get OCGT generators and calculate average marginal cost - ocgt_gens = n.generators[n.generators.carrier == "OCGT"] - avg_marginal_cost = n.get_switchable_as_dense("Generator", "marginal_cost").loc[:, ocgt_gens.index].mean().mean() - n.add("Carrier", "imports", co2_emissions=0.428, nice_name="imports") - - # remove existing oneport components at bus - for c in n.one_port_components: - component = n.components[c].static - rm = component[component.bus.isin(external_buses_to_keep.index)] - if not rm.empty: - logger.info(f"Removing {c} at external buses {external_buses_to_keep.index} with components {rm.index}") - n.remove(c, rm.index) - - # Handle external buses and their generators - for bus in external_buses_to_keep.index: - # Create new import generator - bus_name = n.buses.loc[bus].name - n.add( - "Generator", - f"import_{bus_name}", - bus=bus, - carrier="imports", - p_nom=1e4, - p_nom_extendable=False, - marginal_cost=avg_marginal_cost, - efficiency=1, - build_year=n.investment_periods[0], - lifetime=100, - ) - - # Change location names of external buses, append imports to the ['reeds_state', 'reeds_zone', 'reeds_ba', 'interconnect', 'trans_reg', 'trans_grp'] - n.buses.loc[bus, "reeds_state"] = f"imports_{n.buses.loc[bus, 'reeds_state']}" - n.buses.loc[bus, "reeds_zone"] = f"imports_{n.buses.loc[bus, 'reeds_zone']}" - n.buses.loc[bus, "reeds_ba"] = f"imports_{n.buses.loc[bus, 'reeds_ba']}" - n.buses.loc[bus, "interconnect"] = f"imports_{n.buses.loc[bus, 'interconnect']}" - n.buses.loc[bus, "trans_reg"] = f"imports_{n.buses.loc[bus, 'trans_reg']}" - n.buses.loc[bus, "trans_grp"] = f"imports_{n.buses.loc[bus, 'trans_grp']}" - - # Set all links and lines connected to the bus as non-extendable - for c in ["Line", "Link"]: - attr_name = "p_nom_extendable" if c == "Link" else "s_nom_extendable" - component = n.components[c].static - mask = (component.bus0 == bus) | (component.bus1 == bus) - if mask.any(): - component.loc[mask, attr_name] = False - n.components[c].static.update(component) - - # Remove the links which have "exp" in the name and are connected to the external buses - links_to_remove = n.links[ - n.links.index.str.contains("exp") - & (n.links.bus0.isin(external_buses_to_keep.index) | n.links.bus1.isin(external_buses_to_keep.index)) - ] - n.remove("Link", links_to_remove.index) - - # Update network topology - n.determine_network_topology() - - -def calc_import_export_costs(n: pypsa.Network, carrier: str) -> float: - """Calculates the average marginal cost for a given carrier.""" - gens = n.generators[n.generators.carrier == carrier] - component = "Generator" - if gens.empty: - gens = n.links[n.links.carrier == carrier] - component = "Link" - if gens.empty: - raise ValueError(f"No generators or links found for carrier to calculate imports/exports costs: {carrier}") - costs = n.get_switchable_as_dense(component, "marginal_cost").loc[:, gens.index].mean().mean() - if costs <= 0.01: - raise ValueError( - f"Average marginal cost for {carrier} is less than or equal to 0.01. Check the fuel costs configuration.", - ) - return costs - - -def load_import_export_costs(eia_api: str, year: int) -> pd.DataFrame: - """Loads fuel costs from EIA.""" - return FuelCosts(fuel="electricity", year=year, api=eia_api).get_data() - - -def format_import_export_costs(n: pypsa.Network, fuel_costs: pd.DataFrame) -> pd.DataFrame: - """Formats fuel costs for BA mappings.""" - df = fuel_costs.copy() - data = [] - - buses = n.buses.copy() - - region_mapping = buses.set_index("country")["reeds_state"].to_dict() - for region, state in region_mapping.items(): - for period in df.index.unique(): - temp = df[(df.index == period) & (df.state == state)] - value = temp.value.mean() - data.append([period, region, value, "usd/mwh"]) - formatted = pd.DataFrame(data, columns=["period", "zone", "value", "units"]).set_index("period") - return formatted[~formatted.value.isna()] # regions outside of model scope - - -def format_flowgates_for_imports_exports(n: pypsa.Network, flowgates: pd.DataFrame, zone_col: str) -> pd.DataFrame: - """Formats flowgates for zone mappings.""" - zones_in_model = n.buses[zone_col].unique() - df = flowgates.copy() - - # only keep flowgates that connect inside to outside model scope - df = df[df.r.isin(zones_in_model) ^ df.rr.isin(zones_in_model)] - - # reformat to sinlge value column for easier addition to network - data = [] - for _, row in df.iterrows(): - if row.MW_f0 > 0: - data.append([row.r, row.rr, row.MW_f0]) - if row.MW_r0 > 0: - data.append([row.rr, row.r, row.MW_r0]) - - return pd.DataFrame(data, columns=["r", "rr", "value"]) - - -def convert_flowgates_to_state(flowgates: pd.DataFrame, membership: pd.DataFrame) -> pd.DataFrame: - """Converts flowgates to state level.""" - mbshp = membership.set_index("ba") - df = flowgates.copy() - - df["s"] = df.r.map(mbshp["st"]) - df["ss"] = df.rr.map(mbshp["st"]) - df = df.drop(columns=["r", "rr"]) - df = df.rename(columns={"s": "r", "ss": "rr"}) - return df - - -def add_elec_imports_exports( - n: pypsa.Network, - direction: str, - flowgates: pd.DataFrame, - fuel_costs: pd.DataFrame | float, - co2_emissions: float = 0, - zone_col: str = "reeds_zone", -): - """Add electricity imports and exports to the network. - - These are capacity constrianed links to/from states outside the model spatial scope. - """ - - def _get_regions_2_add(n: pypsa.Network, flowgates: pd.DataFrame, zone_col: str) -> list[str]: - """Gets regions to add import and export buses to.""" - unique_regions = set(flowgates.r.unique()) | set(flowgates.rr.unique()) - return [x for x in unique_regions if x not in n.buses[zone_col].unique()] - - def _add_import_export_carriers(n: pypsa.Network, direction: str, co2_emissions: float | None = None) -> None: - """Adds import and export carriers to the network.""" - if direction == "imports": - co2_emissions = 0 if not co2_emissions else co2_emissions - n.add("Carrier", "imports", co2_emissions=co2_emissions, nice_name="Imports") - elif direction == "exports": - n.add("Carrier", "exports", co2_emissions=0, nice_name="Exports") - else: - raise ValueError(f"direction must be either imports or exports; received: {direction}") - - def _add_import_export_buses(n: pypsa.Network, regions_2_add: list[str], direction: str) -> None: - """Adds import and export buses to the network.""" - if direction == "imports": - suffix = "_imports" - carrier = "imports" - elif direction == "exports": - suffix = "_exports" - carrier = "exports" - else: - raise ValueError(f"direction must be either imports or exports; received: {direction}") - - # cant add in the reeds_state, reeds_zone, reeds_ba, interconnect, trans_reg, trans_grp - # because this information has already been filtered out of the network - - n.add( - "Bus", - regions_2_add, - suffix=suffix, - carrier=carrier, - country=regions_2_add, - ) - - def _add_import_export_stores(n: pypsa.Network, regions_2_add: list[str], direction: str) -> None: - """Adds import and export stores to the network.""" - if direction == "imports": - n.add( - "Store", - regions_2_add, - bus=[f"{x}_imports" for x in regions_2_add], - suffix="_imports", - carrier="imports", - e_nom=0, - e_nom_extendable=True, - capital_cost=0, - e_nom_min=0, - e_nom_max=1e9, - e_min_pu=-1, - e_max_pu=0, - e_cyclic_per_period=False, - marginal_cost=0, - ) - elif direction == "exports": - n.add( - "Store", - regions_2_add, - bus=[f"{x}_exports" for x in regions_2_add], - suffix="_exports", - carrier="exports", - e_nom_extendable=True, - marginal_cost=0, - e_nom=0, - e_nom_max=1e9, - e_min=0, - e_min_pu=0, - e_max_pu=1, - ) - else: - raise ValueError(f"direction must be either imports or exports; received: {direction}") - - def _build_cost_timeseries(n: pypsa.Network, costs: pd.DataFrame, zone: str) -> pd.Series: - """Builds a cost timeseries for a given state.""" - timesteps = n.snapshots.get_level_values("timestep") - years = n.investment_periods - cost_by_zone = costs[costs.zone == zone].drop(columns=["zone", "units"]) - dfs = [] - for year in years: - df = cost_by_zone.copy() - df.index = pd.to_datetime(df.index).map(lambda x: x.replace(year=year)) - df = df.resample("h").ffill().reindex(timesteps).ffill() - df["year"] = year - df = df.set_index(["year", df.index]) # df.index is timestep - dfs.append(df) - df = pd.concat(dfs) - return df.reindex(n.snapshots) - - def _add_import_export_links( - n: pypsa.Network, - flowgates: pd.DataFrame, - fuel_costs: pd.DataFrame | float | str, - direction: str, - zone_col: str = "reeds_zone", - ) -> None: - """Adds import and export links to the network.""" - costs = {} - zones_in_model = n.buses[zone_col].dropna().unique() - - for _, row in flowgates.iterrows(): - zone_inside = row.r if row.r in zones_in_model else row.rr - zone_outside = row.r if row.r not in zones_in_model else row.rr - - # extremely crude caching for generating cost timeseries :| - # keyed by the INSIDE zone — checking the outside zone here skipped - # the write whenever an earlier row's inside zone happened to match, - # leaving costs[zone_inside] unset (KeyError on county networks). - if zone_inside not in costs: - if isinstance(fuel_costs, float | int): - costs[zone_inside] = fuel_costs - elif isinstance(fuel_costs, pd.DataFrame): - costs[zone_inside] = _build_cost_timeseries(n, fuel_costs, zone_inside) - else: - costs[zone_inside] = 0 - - marginal_cost = costs[zone_inside] - - capacity = row.value - - """Structre of flowgates is given by: - - r rr value - 0 p6 p8 488.117 - 1 p8 p6 378.458 - 2 p6 p9 4800.000 - ... - """ - - if direction == "imports": - if row.r == zone_inside: # originating at r is exports (ie r -> rr) - continue - name = f"{zone_inside}_{zone_outside}_imports" - bus0 = f"{zone_outside}_imports" - bus1 = zone_inside - carrier = "imports" - else: - if row.r == zone_outside: # originating at rr is exports (ie rr -> r) - continue - name = f"{zone_inside}_{zone_outside}_exports" - bus0 = zone_inside - bus1 = f"{zone_outside}_exports" - carrier = "exports" - if isinstance(marginal_cost, pd.Series): - marginal_cost = marginal_cost.mul(-1) # constraint will limit exports - - mc = marginal_cost.value if isinstance(marginal_cost, pd.DataFrame) else marginal_cost - - n.add( - "Link", - name, - bus0=bus0, - bus1=bus1, - carrier=carrier, - p_nom_extendable=False, - p_min_pu=0, - p_max_pu=1, - marginal_cost=mc, - p_nom=capacity, - ) - - assert direction in ["imports", "exports"], f"direction must be either imports or exports; received: {direction}" - - regions_2_add = _get_regions_2_add(n, flowgates, zone_col) - _add_import_export_carriers(n, direction, co2_emissions) - _add_import_export_buses(n, regions_2_add, direction) - _add_import_export_stores(n, regions_2_add, direction) - _add_import_export_links(n, flowgates, fuel_costs, direction, zone_col) - - def add_co2_storage(n: pypsa.Network, config: dict, co2_storage_csv: str, costs: pd.DataFrame, sector: bool): """Adds node level CO2 (underground) storage.""" # get node level CO2 (underground) storage potential and cost from CSV file @@ -1556,23 +1197,16 @@ def main(snakemake) -> None: if dr_config: add_demand_response(n, dr_config) - trim_network_config = snakemake.params.trim_network imports_config = snakemake.params.imports exports_config = snakemake.params.exports - - assert not ( - snakemake.params.trim_network and (imports_config.get("enable", False) or exports_config.get("enable", False)) - ), "trim_network and imports/exports cannot be used together" - - if snakemake.params.trim_network: - trim_network(n, trim_network_config) + representation = imports_config.get("representation", "store") if snakemake.params.transmission_network == "reeds": # flowgates to limit the capacity (removed later if configured capacity limit is inf) flowgates = pd.read_csv(snakemake.input.flowgates) + membership = pd.read_csv(snakemake.input.reeds_memberships) if snakemake.params.topological_boundaries == "state": zone_col = "reeds_state" - membership = pd.read_csv(snakemake.input.reeds_memberships) flowgates = convert_flowgates_to_state(flowgates, membership) flowgates = format_flowgates_for_imports_exports(n, flowgates, zone_col) flowgates = flowgates.groupby(["r", "rr"], as_index=False).sum() @@ -1597,22 +1231,24 @@ def main(snakemake) -> None: if not imports_config.get("capacity_limit", True): import_flowgates["value"] = np.inf - import_costs = imports_config.get("costs", False) + fuel_costs = resolve_trade_costs(n, imports_config, "imports", snakemake.params.eia_api, year) - if isinstance(import_costs, float | int): # user defined value - fuel_costs = import_costs - elif isinstance(import_costs, str): # 'wholesale' or name of carrier - if import_costs == "wholesale": - fuel_costs = load_import_export_costs(snakemake.params.eia_api, year) - fuel_costs = format_import_export_costs(n, fuel_costs) - else: - fuel_costs = calc_import_export_costs(n, import_costs) - else: - raise ValueError( - f"'imports.costs' must be 'wholesale', name of a carrier, or a float/int. Received: {import_costs}", - ) + # Only `generator` mode places the CPUC contracted units behind the + # boundary; in `store` mode add_electricity has already attached them + # and the bundle is a sentinel. + remote_bundle = load_remote_unit_bundle(snakemake.input.remote_units) if representation == "generator" else None - add_elec_imports_exports(n, "imports", import_flowgates, fuel_costs, co2_emissions, zone_col) + add_external_regions( + n, + "imports", + representation, + import_flowgates, + fuel_costs, + co2_emissions, + zone_col, + remote_bundle=remote_bundle, + membership=membership, + ) # Electricity exports configuration if exports_config.get("enable", False) and snakemake.params.transmission_network == "reeds": @@ -1627,25 +1263,17 @@ def main(snakemake) -> None: if not exports_config.get("capacity_limit", True): export_flowgates["value"] = np.inf - export_costs = exports_config.get("costs", False) + fuel_costs = resolve_trade_costs(n, exports_config, "exports", snakemake.params.eia_api, year) - if isinstance(export_costs, float | int): # user defined value - fuel_costs = export_costs - fuel_costs *= -1 # make money by exporting - elif isinstance(export_costs, str): # 'wholesale' or name of carrier - if export_costs == "wholesale": - fuel_costs = load_import_export_costs(snakemake.params.eia_api, year) - fuel_costs = format_import_export_costs(n, fuel_costs) - fuel_costs["value"] = fuel_costs.value.mul(-1) # make money by exporting - else: - fuel_costs = calc_import_export_costs(n, export_costs) - fuel_costs *= -1 # make money by exporting - else: - raise ValueError( - f"'exports.costs' must be 'wholesale', name of a carrier, or a float/int. Received: {export_costs}", - ) - - add_elec_imports_exports(n, "exports", export_flowgates, fuel_costs, co2_emissions, zone_col) + add_external_regions( + n, + "exports", + representation, + export_flowgates, + fuel_costs, + co2_emissions, + zone_col, + ) if snakemake.config["scenario"]["sector"] == "E": co2_storage = snakemake.config.get("co2", {}).get("storage", False) diff --git a/workflow/scripts/external_regions.py b/workflow/scripts/external_regions.py new file mode 100644 index 000000000..f05c3f11d --- /dev/null +++ b/workflow/scripts/external_regions.py @@ -0,0 +1,691 @@ +"""External (out-of-footprint) regions: imports, exports and contracted units. + +A regionally scoped run (California, say) cuts the synchronous grid at a +political boundary. Everything behind that cut still serves load inside the +footprint, and this module is the single place where it is represented. Two +representations are available, selected by ``electricity.imports.representation``: + +``store`` (default, unchanged behaviour) + Per external flowgate zone, a ``{zone}_imports`` bus carrying a bottomless + ``Store`` (``e_nom_max`` 1e9, ``e_min_pu`` -1) and one-way ``Link`` s into + the internal zone buses. Each link is rated at the NARIS flowgate capacity + and PRICED at the import price (a ``wholesale`` EIA timeseries, a carrier's + average marginal cost, or a flat float). Emissions are carried by the + ``imports`` carrier that the Store belongs to. Exports mirror this with an + absorbing Store and negatively-priced links. + +``generator`` + The external zone is modelled as a place with generation rather than as an + infinite energy tank: + + * the ``{zone}_imports`` bus carries a generic import ``Generator`` + (carrier ``unspecified_imports``) sized at the zone's total inbound + interface capacity and priced with the same machinery as ``store`` mode; + * California's CPUC-contracted out-of-state units (Palo Verde, Intermountain, + Hoover, Apex, the AZ/NV solar and battery contracts) are attached at that + same external bus instead of at a California bus, i.e. BEHIND the boundary; + * the ``Link`` s into the footprint are UNPRICED and rated at the flowgate + capacity, so the price sits on the generator and the link is pure transfer + capacity. + + The point of the second mode is accounting: because deliveries now traverse + a carrier-``imports`` link, they are counted by + :func:`opts.interchange.add_interchange_constraints` against the import + volume cap, and bounded by the interface capacity, exactly like generic + imports. In ``store`` mode the contracted units sit inside the footprint and + look like in-state generation, bypassing both limits. + +Both modes keep the link carriers ``imports`` / ``exports`` untouched, which is +what ``opts/interchange.py`` and ``opts/interfaces.py`` key off. + +External-bus naming +------------------- +Both modes use SEPARATE ``{zone}_imports`` and ``{zone}_exports`` buses rather +than one shared ``{zone}_external`` bus. Beyond keeping ``store`` mode +byte-identical, the separation is load-bearing in ``generator`` mode: a shared +bus would let the priced import generator (and the contracted units) sell +straight into the export sink, collecting the export price for free whenever the +export price exceeds the import cost — a pure arbitrage loop with no physical +meaning. + +Export pricing +-------------- +Export pricing is IDENTICAL in both modes: the negative price stays on the +export ``Link`` and the absorbing ``Store`` behind it is free. On the export +side there is no double-charge to worry about, because the only way energy can +reach the export Store is through exactly one export link — nothing else injects +into a ``{zone}_exports`` bus — so the negative price is earned exactly once per +exported MWh. That is what makes it safe to leave the export half of the +construction untouched by the representation switch, and it is only safe because +the import generator lives on a different bus (see above). + +CO2 +--- +In ``store`` mode the ``imports`` carrier carries ``imports.co2_emissions`` and +the Store's withdrawal is what the global CO2 constraint sees. In ``generator`` +mode there is no import Store, so the emission factor moves onto the +``unspecified_imports`` carrier of the generic import generator and ``imports`` +is set to zero. PyPSA attributes primary-energy emissions to generators and +stores, never to links, so the carrier-``imports`` links never double-count. +""" + +import logging + +import dill +import pandas as pd +import pypsa +from add_electricity import add_missing_carriers, attach_remote_units +from eia import FuelCosts + +logger = logging.getLogger(__name__) + +REPRESENTATIONS = ("store", "generator") + +#: Carrier of the generic import generator in ``generator`` mode. Deliberately +#: NOT ``imports``: that carrier is reserved for the transfer links which +#: ``opts/interchange.py`` sums over. +GENERIC_IMPORT_CARRIER = "unspecified_imports" + + +# --------------------------------------------------------------------------- +# Flowgate formatting +# --------------------------------------------------------------------------- + + +def format_flowgates_for_imports_exports(n: pypsa.Network, flowgates: pd.DataFrame, zone_col: str) -> pd.DataFrame: + """Formats flowgates for zone mappings.""" + zones_in_model = n.buses[zone_col].unique() + df = flowgates.copy() + + # only keep flowgates that connect inside to outside model scope + df = df[df.r.isin(zones_in_model) ^ df.rr.isin(zones_in_model)] + + # reformat to sinlge value column for easier addition to network + data = [] + for _, row in df.iterrows(): + if row.MW_f0 > 0: + data.append([row.r, row.rr, row.MW_f0]) + if row.MW_r0 > 0: + data.append([row.rr, row.r, row.MW_r0]) + + return pd.DataFrame(data, columns=["r", "rr", "value"]) + + +def convert_flowgates_to_state(flowgates: pd.DataFrame, membership: pd.DataFrame) -> pd.DataFrame: + """Converts flowgates to state level.""" + mbshp = membership.set_index("ba") + df = flowgates.copy() + + df["s"] = df.r.map(mbshp["st"]) + df["ss"] = df.rr.map(mbshp["st"]) + df = df.drop(columns=["r", "rr"]) + df = df.rename(columns={"s": "r", "ss": "rr"}) + return df + + +# --------------------------------------------------------------------------- +# Import / export pricing +# --------------------------------------------------------------------------- + + +def calc_import_export_costs(n: pypsa.Network, carrier: str) -> float: + """Calculates the average marginal cost for a given carrier.""" + gens = n.generators[n.generators.carrier == carrier] + component = "Generator" + if gens.empty: + gens = n.links[n.links.carrier == carrier] + component = "Link" + if gens.empty: + raise ValueError(f"No generators or links found for carrier to calculate imports/exports costs: {carrier}") + costs = n.get_switchable_as_dense(component, "marginal_cost").loc[:, gens.index].mean().mean() + if costs <= 0.01: + raise ValueError( + f"Average marginal cost for {carrier} is less than or equal to 0.01. Check the fuel costs configuration.", + ) + return costs + + +def load_import_export_costs(eia_api: str, year: int) -> pd.DataFrame: + """Loads fuel costs from EIA.""" + return FuelCosts(fuel="electricity", year=year, api=eia_api).get_data() + + +def format_import_export_costs(n: pypsa.Network, fuel_costs: pd.DataFrame) -> pd.DataFrame: + """Formats fuel costs for BA mappings.""" + df = fuel_costs.copy() + data = [] + + buses = n.buses.copy() + + region_mapping = buses.set_index("country")["reeds_state"].to_dict() + for region, state in region_mapping.items(): + for period in df.index.unique(): + temp = df[(df.index == period) & (df.state == state)] + value = temp.value.mean() + data.append([period, region, value, "usd/mwh"]) + formatted = pd.DataFrame(data, columns=["period", "zone", "value", "units"]).set_index("period") + return formatted[~formatted.value.isna()] # regions outside of model scope + + +def resolve_trade_costs( + n: pypsa.Network, + trade_config: dict, + direction: str, + eia_api: str, + year: int, +) -> pd.DataFrame | float: + """Resolve ``imports.costs`` / ``exports.costs`` into a price the network can use. + + ``wholesale`` pulls the monthly EIA electricity price per state, a carrier + name averages that carrier's marginal cost, and a float is taken at face + value. Export prices are negated: exporting earns money. + """ + if direction not in ("imports", "exports"): + raise ValueError(f"direction must be either imports or exports; received: {direction}") + + sign = 1 if direction == "imports" else -1 + costs = trade_config.get("costs", False) + + if isinstance(costs, float | int): # user defined value + return costs * sign + if isinstance(costs, str): # 'wholesale' or name of carrier + if costs == "wholesale": + fuel_costs = load_import_export_costs(eia_api, year) + fuel_costs = format_import_export_costs(n, fuel_costs) + if sign < 0: + fuel_costs["value"] = fuel_costs.value.mul(sign) # make money by exporting + return fuel_costs + return calc_import_export_costs(n, costs) * sign + raise ValueError( + f"'{direction}.costs' must be 'wholesale', name of a carrier, or a float/int. Received: {costs}", + ) + + +def _build_cost_timeseries(n: pypsa.Network, costs: pd.DataFrame, zone: str) -> pd.Series: + """Builds a cost timeseries for a given state.""" + timesteps = n.snapshots.get_level_values("timestep") + years = n.investment_periods + cost_by_zone = costs[costs.zone == zone].drop(columns=["zone", "units"]) + dfs = [] + for year in years: + df = cost_by_zone.copy() + df.index = pd.to_datetime(df.index).map(lambda x: x.replace(year=year)) + df = df.resample("h").ffill().reindex(timesteps).ffill() + df["year"] = year + df = df.set_index(["year", df.index]) # df.index is timestep + dfs.append(df) + df = pd.concat(dfs) + return df.reindex(n.snapshots) + + +# --------------------------------------------------------------------------- +# Shared construction helpers +# --------------------------------------------------------------------------- + + +def _get_regions_2_add(n: pypsa.Network, flowgates: pd.DataFrame, zone_col: str) -> list[str]: + """Gets regions to add import and export buses to.""" + unique_regions = set(flowgates.r.unique()) | set(flowgates.rr.unique()) + return [x for x in unique_regions if x not in n.buses[zone_col].unique()] + + +def _add_import_export_carriers(n: pypsa.Network, direction: str, co2_emissions: float | None = None) -> None: + """Adds import and export carriers to the network.""" + if direction == "imports": + co2_emissions = 0 if not co2_emissions else co2_emissions + n.add("Carrier", "imports", co2_emissions=co2_emissions, nice_name="Imports") + elif direction == "exports": + n.add("Carrier", "exports", co2_emissions=0, nice_name="Exports") + else: + raise ValueError(f"direction must be either imports or exports; received: {direction}") + + +def _add_import_export_buses(n: pypsa.Network, regions_2_add: list[str], direction: str) -> None: + """Adds import and export buses to the network.""" + if direction == "imports": + suffix = "_imports" + carrier = "imports" + elif direction == "exports": + suffix = "_exports" + carrier = "exports" + else: + raise ValueError(f"direction must be either imports or exports; received: {direction}") + + # cant add in the reeds_state, reeds_zone, reeds_ba, interconnect, trans_reg, trans_grp + # because this information has already been filtered out of the network + + n.add( + "Bus", + regions_2_add, + suffix=suffix, + carrier=carrier, + country=regions_2_add, + ) + + +def _add_import_export_stores(n: pypsa.Network, regions_2_add: list[str], direction: str) -> None: + """Adds import and export stores to the network.""" + if direction == "imports": + n.add( + "Store", + regions_2_add, + bus=[f"{x}_imports" for x in regions_2_add], + suffix="_imports", + carrier="imports", + e_nom=0, + e_nom_extendable=True, + capital_cost=0, + e_nom_min=0, + e_nom_max=1e9, + e_min_pu=-1, + e_max_pu=0, + e_cyclic_per_period=False, + marginal_cost=0, + ) + elif direction == "exports": + n.add( + "Store", + regions_2_add, + bus=[f"{x}_exports" for x in regions_2_add], + suffix="_exports", + carrier="exports", + e_nom_extendable=True, + marginal_cost=0, + e_nom=0, + e_nom_max=1e9, + e_min=0, + e_min_pu=0, + e_max_pu=1, + ) + else: + raise ValueError(f"direction must be either imports or exports; received: {direction}") + + +def _add_import_export_links( + n: pypsa.Network, + flowgates: pd.DataFrame, + fuel_costs: pd.DataFrame | float | str, + direction: str, + zone_col: str = "reeds_zone", + priced: bool = True, +) -> None: + """Adds import and export links to the network. + + ``priced`` is what distinguishes the two representations on the import side: + in ``store`` mode the link carries the import price, in ``generator`` mode + the price sits on the external generator instead and the link is a pure + transfer capacity. + """ + costs = {} + zones_in_model = n.buses[zone_col].dropna().unique() + + for _, row in flowgates.iterrows(): + zone_inside = row.r if row.r in zones_in_model else row.rr + zone_outside = row.r if row.r not in zones_in_model else row.rr + + # extremely crude caching for generating cost timeseries :| + # keyed by the INSIDE zone — checking the outside zone here skipped + # the write whenever an earlier row's inside zone happened to match, + # leaving costs[zone_inside] unset (KeyError on county networks). + if zone_inside not in costs: + if isinstance(fuel_costs, float | int): + costs[zone_inside] = fuel_costs + elif isinstance(fuel_costs, pd.DataFrame): + costs[zone_inside] = _build_cost_timeseries(n, fuel_costs, zone_inside) + else: + costs[zone_inside] = 0 + + marginal_cost = costs[zone_inside] + + capacity = row.value + + """Structre of flowgates is given by: + + r rr value + 0 p6 p8 488.117 + 1 p8 p6 378.458 + 2 p6 p9 4800.000 + ... + """ + + if direction == "imports": + if row.r == zone_inside: # originating at r is exports (ie r -> rr) + continue + name = f"{zone_inside}_{zone_outside}_imports" + bus0 = f"{zone_outside}_imports" + bus1 = zone_inside + carrier = "imports" + if not priced: + marginal_cost = 0 + else: + if row.r == zone_outside: # originating at rr is exports (ie rr -> r) + continue + name = f"{zone_inside}_{zone_outside}_exports" + bus0 = zone_inside + bus1 = f"{zone_outside}_exports" + carrier = "exports" + if isinstance(marginal_cost, pd.Series): + marginal_cost = marginal_cost.mul(-1) # constraint will limit exports + + mc = marginal_cost.value if isinstance(marginal_cost, pd.DataFrame) else marginal_cost + + n.add( + "Link", + name, + bus0=bus0, + bus1=bus1, + carrier=carrier, + p_nom_extendable=False, + p_min_pu=0, + p_max_pu=1, + marginal_cost=mc, + p_nom=capacity, + ) + + +# --------------------------------------------------------------------------- +# `store` representation +# --------------------------------------------------------------------------- + + +def add_elec_imports_exports( + n: pypsa.Network, + direction: str, + flowgates: pd.DataFrame, + fuel_costs: pd.DataFrame | float, + co2_emissions: float = 0, + zone_col: str = "reeds_zone", +): + """Add electricity imports and exports to the network. + + These are capacity constrianed links to/from states outside the model spatial scope. + """ + assert direction in ["imports", "exports"], f"direction must be either imports or exports; received: {direction}" + + regions_2_add = _get_regions_2_add(n, flowgates, zone_col) + _add_import_export_carriers(n, direction, co2_emissions) + _add_import_export_buses(n, regions_2_add, direction) + _add_import_export_stores(n, regions_2_add, direction) + _add_import_export_links(n, flowgates, fuel_costs, direction, zone_col) + + +# --------------------------------------------------------------------------- +# `generator` representation +# --------------------------------------------------------------------------- + + +def inbound_capacity_by_zone( + n: pypsa.Network, + flowgates: pd.DataFrame, + zone_col: str = "reeds_zone", +) -> pd.Series: + """Total interface capacity flowing INTO the footprint, per external zone.""" + zones_in_model = n.buses[zone_col].dropna().unique() + inbound = flowgates[~flowgates.r.isin(zones_in_model) & flowgates.rr.isin(zones_in_model)] + return inbound.groupby("r")["value"].sum() + + +def _internal_zones_served( + n: pypsa.Network, + flowgates: pd.DataFrame, + zone_col: str, +) -> dict[str, list[str]]: + """Map each external zone onto the internal zones it can deliver into.""" + zones_in_model = n.buses[zone_col].dropna().unique() + inbound = flowgates[~flowgates.r.isin(zones_in_model) & flowgates.rr.isin(zones_in_model)] + return {zone: sorted(set(rows.rr)) for zone, rows in inbound.groupby("r")} + + +def _external_generator_cost( + n: pypsa.Network, + fuel_costs: pd.DataFrame | float, + internal_zones: list[str], +) -> pd.Series | float: + """Price for the generic import generator of one external zone. + + The wholesale price table is keyed by the zones INSIDE the model (it is built + from bus attributes), so an external zone has no price of its own. Its + generator is priced at the mean of the internal zones it can reach, which is + the same set of prices the ``store``-mode links would have carried. + """ + if isinstance(fuel_costs, float | int): + return fuel_costs + if not isinstance(fuel_costs, pd.DataFrame): + return 0 + series = [_build_cost_timeseries(n, fuel_costs, zone)["value"] for zone in internal_zones] + series = [s for s in series if not s.isna().all()] + if not series: + return 0 + return pd.concat(series, axis=1).mean(axis=1) + + +def _add_generic_import_generators( + n: pypsa.Network, + regions_2_add: list[str], + flowgates: pd.DataFrame, + fuel_costs: pd.DataFrame | float, + co2_emissions: float, + zone_col: str, +) -> None: + """One priced, non-extendable import generator per external bus.""" + n.add("Carrier", GENERIC_IMPORT_CARRIER, co2_emissions=co2_emissions, nice_name="Unspecified Imports") + + capacity = inbound_capacity_by_zone(n, flowgates, zone_col) + served = _internal_zones_served(n, flowgates, zone_col) + + for zone in regions_2_add: + p_nom = float(capacity.get(zone, 0.0)) + if p_nom <= 0: + logger.info(f"External zone '{zone}' has no inbound interface capacity; no import generator added.") + continue + marginal_cost = _external_generator_cost(n, fuel_costs, served.get(zone, [])) + n.add( + "Generator", + f"{zone}_imports {GENERIC_IMPORT_CARRIER}", + bus=f"{zone}_imports", + carrier=GENERIC_IMPORT_CARRIER, + p_nom=p_nom, + p_nom_extendable=False, + efficiency=1, + marginal_cost=marginal_cost, + ) + + +def map_remote_units_to_zones( + unit_states: pd.Series, + external_zones: list[str], + inbound_capacity: pd.Series, + zone_col: str, + membership: pd.DataFrame | None = None, +) -> pd.Series: + """Assign each contracted remote unit to the external zone it sits behind. + + Candidate zones are the boundary zones whose state matches the unit's + physical state (from ``powerplants.csv``); the candidate with the largest + inbound interface capacity wins. A unit whose state has no direct interface + with the footprint (Utah's Intermountain, say) falls back to the boundary + zone with the largest inbound capacity overall, with a warning. + """ + ranked = [z for z in external_zones if inbound_capacity.get(z, 0.0) > 0] + ranked = sorted(ranked, key=lambda z: (-float(inbound_capacity.get(z, 0.0)), z)) + if not ranked: + raise ValueError("No external zone with inbound interface capacity; cannot place remote contracted units.") + + zone_state = _external_zone_states(ranked, zone_col, membership) + fallback = ranked[0] + + assignment = {} + for unit, state in unit_states.items(): + candidates = [z for z in ranked if zone_state.get(z) == state] + if candidates: + assignment[unit] = candidates[0] + else: + assignment[unit] = fallback + logger.warning( + f"Remote contracted unit '{unit}' sits in state '{state}', which has no direct interface with the " + f"model footprint; placing it behind external zone '{fallback}' (largest inbound capacity).", + ) + return pd.Series(assignment, dtype=object) + + +def _external_zone_states( + external_zones: list[str], + zone_col: str, + membership: pd.DataFrame | None, +) -> dict[str, str]: + """State code of each external zone. + + With ``topological_boundaries: state`` the zone IS the state code; otherwise + the ReEDS membership table maps the balancing area onto its state. + """ + if zone_col == "reeds_state": + return {zone: zone for zone in external_zones} + if membership is None: + logger.warning("No ReEDS membership table supplied; cannot map external zones onto states.") + return {} + mbshp = membership.set_index("ba")["st"] + return {zone: mbshp.get(zone) for zone in external_zones} + + +def attach_remote_contracted_units_externally( + n: pypsa.Network, + bundle: dict, + flowgates: pd.DataFrame, + zone_col: str, + membership: pd.DataFrame | None, +) -> pd.Series: + """Attach the serialized CPUC contracted units at their external buses. + + ``bundle`` is the ``add_electricity`` output described in + :func:`add_electricity.build_remote_unit_bundle`: the fully-derived unit + table (still keyed to the California bus whose VRE profile it borrowed) plus + the borrowed profiles and the cost table needed to resolve capital costs. + Only the ``bus`` column is rewritten here. + """ + units = bundle["units"].copy() + if units.empty: + return pd.Series(dtype=object) + + external_zones = [b[: -len("_imports")] for b in n.buses.index[n.buses.index.str.endswith("_imports")]] + inbound = inbound_capacity_by_zone(n, flowgates, zone_col) + zones = map_remote_units_to_zones(units["state"], external_zones, inbound, zone_col, membership) + + units["zone"] = zones + units["bus"] = zones.map(lambda z: f"{z}_imports") + + for name, unit in units.iterrows(): + logger.info( + f"Remote contracted unit '{name}' ({unit['carrier']}, {unit['p_nom']:.1f} MW, state {unit['state']}) " + f"attached behind external zone '{unit['zone']}' at bus '{unit['bus']}'.", + ) + + add_missing_carriers(n, sorted(set(units["carrier"]))) + attach_remote_units( + n, + units.drop(columns=["zone"]), + bundle["costs"], + bundle["conventional_carriers"], + bundle["unit_commitment"], + profiles=bundle.get("vre_profiles"), + ) + return units["zone"] + + +def _add_generator_representation( + n: pypsa.Network, + flowgates: pd.DataFrame, + fuel_costs: pd.DataFrame | float, + co2_emissions: float, + zone_col: str, + remote_bundle: dict | None, + membership: pd.DataFrame | None, +) -> None: + """Import side of the ``generator`` representation (exports are unchanged).""" + regions_2_add = _get_regions_2_add(n, flowgates, zone_col) + + # emissions move onto the generic import generator's carrier; the `imports` + # carrier now only labels links, which PyPSA never charges emissions to. + _add_import_export_carriers(n, "imports", 0) + _add_import_export_buses(n, regions_2_add, "imports") + _add_generic_import_generators(n, regions_2_add, flowgates, fuel_costs, co2_emissions, zone_col) + _add_import_export_links(n, flowgates, fuel_costs, "imports", zone_col, priced=False) + + if remote_bundle: + attach_remote_contracted_units_externally(n, remote_bundle, flowgates, zone_col, membership) + + +# --------------------------------------------------------------------------- +# Public entry point +# --------------------------------------------------------------------------- + + +def add_external_regions( + n: pypsa.Network, + direction: str, + representation: str, + flowgates: pd.DataFrame, + fuel_costs: pd.DataFrame | float, + co2_emissions: float = 0, + zone_col: str = "reeds_zone", + remote_bundle: dict | None = None, + membership: pd.DataFrame | None = None, +) -> None: + """Add the external-region representation for one direction of trade. + + Parameters + ---------- + direction + ``imports`` or ``exports``. + representation + ``store`` or ``generator`` — see the module docstring. + flowgates + Formatted NARIS flowgates (columns ``r``, ``rr``, ``value``), already + restricted to interfaces that cross the model boundary. + fuel_costs + Output of :func:`resolve_trade_costs`. + remote_bundle, membership + Only used for ``direction='imports'`` in ``generator`` mode: the CPUC + contracted-unit bundle written by ``add_electricity`` and the ReEDS + ``membership.csv`` used to map external zones onto states. + """ + if direction not in ("imports", "exports"): + raise ValueError(f"direction must be either imports or exports; received: {direction}") + if representation not in REPRESENTATIONS: + raise ValueError(f"representation must be one of {REPRESENTATIONS}; received: {representation}") + + if representation == "store" or direction == "exports": + add_elec_imports_exports(n, direction, flowgates, fuel_costs, co2_emissions, zone_col) + return + + _add_generator_representation(n, flowgates, fuel_costs, co2_emissions, zone_col, remote_bundle, membership) + + +def load_remote_unit_bundle(path: str | None) -> dict | None: + """Read the contracted-unit bundle written by ``add_electricity``. + + The rule always declares the file, so a run with contracted resources + disabled (or in ``store`` mode) writes a sentinel ``None``; that is not an + error, it just means there is nothing to place behind the boundary. + """ + if not path: + return None + with open(path, "rb") as f: + bundle = dill.load(f) + if not bundle or bundle.get("units") is None or bundle["units"].empty: + return None + return bundle + + +__all__ = [ + "add_elec_imports_exports", + "add_external_regions", + "calc_import_export_costs", + "convert_flowgates_to_state", + "format_flowgates_for_imports_exports", + "format_import_export_costs", + "inbound_capacity_by_zone", + "load_import_export_costs", + "load_remote_unit_bundle", + "map_remote_units_to_zones", + "resolve_trade_costs", +] diff --git a/workflow/scripts/opts/interfaces.py b/workflow/scripts/opts/interfaces.py index 4dd59c2fe..15fcacc1a 100644 --- a/workflow/scripts/opts/interfaces.py +++ b/workflow/scripts/opts/interfaces.py @@ -10,7 +10,7 @@ CAISO_Imports,"p9, p10, p11","p2, p5, p6, ...",9728,10208,RESOLVE The caps are applied to the import/export ``Link`` components created by -``add_extra_components.add_elec_imports_exports`` and are therefore a no-op +``external_regions.add_external_regions`` and are therefore a no-op when ``electricity.imports``/``electricity.exports`` are disabled. """ diff --git a/workflow/scripts/plot_statistics.py b/workflow/scripts/plot_statistics.py index 8041386f5..54314a243 100644 --- a/workflow/scripts/plot_statistics.py +++ b/workflow/scripts/plot_statistics.py @@ -229,8 +229,8 @@ def plot_capacity_additions_bar( optimal_capacity = optimal_capacity.fillna(0) # Drop the synthetic "imports" carrier — it represents external power - # injection (from trim_network or the imports/exports config), not built - # capacity, and would otherwise dominate the bar by orders of magnitude. + # injection (from the imports/exports config), not built capacity, and would + # otherwise dominate the bar by orders of magnitude. hidden_carriers = {"imports", "Imports"} optimal_capacity = optimal_capacity.drop( index=[c for c in hidden_carriers if c in optimal_capacity.index], From b2e630e18e8d09eeeccce9905159a93881cdeee9 Mon Sep 17 00:00:00 2001 From: Kamran Date: Mon, 31 Aug 2026 22:31:45 -0700 Subject: [PATCH 08/13] test: cover both external-region representations Extends test_interfaces.py with a toy one-zone footprint and a NARIS-shaped flowgate frame (two NV zones, one AZ zone, one outbound-only interface): * generator mode creates the external buses, a non-extendable import generator sized at the zone's inbound interface capacity, and UNPRICED interface links; store mode keeps the priced links and the bottomless Store; * the CO2 factor moves onto `unspecified_imports` in generator mode and stays on `imports` in store mode; * exports are byte-identical across representations, and nothing but export links injects into an export bus; * a toy contracted-unit bundle attaches at the state-matched external bus (AZ -> the AZ boundary zone), falls back with a warning when the unit's state has no direct interface (UT -> the largest inbound zone), and carries its borrowed VRE profile across. Co-Authored-By: Claude Opus 5 --- workflow/scripts/test/test_interfaces.py | 346 ++++++++++++++++++++++- 1 file changed, 343 insertions(+), 3 deletions(-) diff --git a/workflow/scripts/test/test_interfaces.py b/workflow/scripts/test/test_interfaces.py index 91b0f9c56..e72655d42 100644 --- a/workflow/scripts/test/test_interfaces.py +++ b/workflow/scripts/test/test_interfaces.py @@ -1,13 +1,16 @@ """ -Test the aggregate transmission interface limits. +Test the aggregate transmission interface limits and the external-region build. This module contains tests for the RESOLVE/NARIS style interface constraints -applied to the electricity import/export links in PyPSA-USA. +applied to the electricity import/export links in PyPSA-USA, and for the two +representations of out-of-footprint supply built by ``external_regions``. """ +import logging import os import sys +import numpy as np import pandas as pd import pypsa import pytest @@ -15,6 +18,12 @@ sys.path.append(os.path.join(os.path.dirname(__file__), "..")) from _helpers import get_multiindex_snapshots +from external_regions import ( + GENERIC_IMPORT_CARRIER, + add_external_regions, + inbound_capacity_by_zone, + map_remote_units_to_zones, +) from opts.interfaces import ( _boundary_links, _parse_regions, @@ -36,7 +45,7 @@ def interface_network(): """ Build a small network with electricity import and export links. - Mirrors the conventions of ``add_extra_components.add_elec_imports_exports``: + Mirrors the conventions of ``external_regions.add_elec_imports_exports``: external buses are named ``{zone}_imports`` / ``{zone}_exports`` with the matching carrier, import links run from the external bus into the model and export links run the other way. @@ -285,3 +294,334 @@ def extra_functionality(n, sns): n.optimize(solver_name="glpk", multi_investment_periods=True, extra_functionality=extra_functionality) assert not [c for c in n.model.constraints if c.startswith("interface_limit-")] + + +# --------------------------------------------------------------------------- +# external_regions: the two import representations +# --------------------------------------------------------------------------- + +MEMBERSHIP = pd.DataFrame( + [ + {"ba": "p9", "st": "CA"}, + {"ba": "p12", "st": "NV"}, + {"ba": "p13", "st": "NV"}, + {"ba": "p28", "st": "AZ"}, + ], +) + +# Mirrors the CA boundary: two NV zones (p13 much larger than p12), one AZ zone, +# and one outbound-only interface so the inbound/outbound split is exercised. +FLOWGATES = pd.DataFrame( + [ + {"r": "p13", "rr": "p9", "value": 2000.0}, + {"r": "p12", "rr": "p9", "value": 500.0}, + {"r": "p28", "rr": "p9", "value": 1000.0}, + {"r": "p9", "rr": "p13", "value": 1500.0}, + ], +) + + +@pytest.fixture +def footprint_network(): + """A one-zone footprint (California's p9) with nothing external attached yet.""" + n = pypsa.Network() + n.snapshots = get_multiindex_snapshots( + sns_config={"start": "2030-01-01 00:00", "end": "2030-01-01 03:00", "inclusive": "both"}, + invest_periods=[2030], + ) + n.set_investment_periods(periods=[2030]) + + n.add("Carrier", "AC", co2_emissions=0) + n.add("Carrier", "solar", co2_emissions=0) + n.add( + "Bus", + "p9", + carrier="AC", + country="p9", + interconnect="western", + reeds_state="CA", + reeds_zone="p9", + ) + n.add( + "Generator", + "p9 solar", + bus="p9", + carrier="solar", + p_nom=100, + p_max_pu=pd.Series(0.5, index=n.snapshots), + ) + return n + + +def wholesale_costs(): + """A `format_import_export_costs`-shaped price table for the internal zone.""" + return pd.DataFrame( + {"zone": ["p9"], "value": [42.0], "units": ["usd/mwh"]}, + index=pd.to_datetime(["2030-01-01"]), + ) + + +def import_links(n): + return n.links[n.links.carrier == "imports"] + + +def test_inbound_capacity_by_zone_ignores_outbound_rows(footprint_network): + inbound = inbound_capacity_by_zone(footprint_network, FLOWGATES, "reeds_zone") + assert inbound.to_dict() == {"p12": 500.0, "p13": 2000.0, "p28": 1000.0} + + +def test_generator_mode_builds_external_buses_generators_and_links(footprint_network): + n = footprint_network + add_external_regions(n, "imports", "generator", FLOWGATES, 30.0, co2_emissions=0.428, zone_col="reeds_zone") + + for zone in ("p12", "p13", "p28"): + assert f"{zone}_imports" in n.buses.index + assert n.buses.at[f"{zone}_imports", "carrier"] == "imports" + + gens = n.generators[n.generators.carrier == GENERIC_IMPORT_CARRIER] + assert sorted(gens.bus) == ["p12_imports", "p13_imports", "p28_imports"] + assert not gens.p_nom_extendable.any() + # p_nom is the zone's total INBOUND interface capacity + assert gens.set_index("bus")["p_nom"].to_dict() == { + "p12_imports": 500.0, + "p13_imports": 2000.0, + "p28_imports": 1000.0, + } + assert (gens.marginal_cost == 30.0).all() + + links = import_links(n) + assert sorted(links.index) == ["p9_p12_imports", "p9_p13_imports", "p9_p28_imports"] + assert links.set_index("bus0")["p_nom"].to_dict() == { + "p12_imports": 500.0, + "p13_imports": 2000.0, + "p28_imports": 1000.0, + } + assert (links.bus1 == "p9").all() + assert not links.p_nom_extendable.any() + + # generator mode prices the GENERATOR, so the interface link is free + assert (links.marginal_cost == 0).all() + # and there is no bottomless import Store + assert n.stores[n.stores.carrier == "imports"].empty + + +def test_store_mode_prices_the_links_and_keeps_the_store(footprint_network): + n = footprint_network + add_external_regions(n, "imports", "store", FLOWGATES, 30.0, co2_emissions=0.428, zone_col="reeds_zone") + + links = import_links(n) + assert (links.marginal_cost == 30.0).all() + assert sorted(n.stores[n.stores.carrier == "imports"].index) == [ + "p12_imports", + "p13_imports", + "p28_imports", + ] + assert n.generators[n.generators.carrier == GENERIC_IMPORT_CARRIER].empty + + +def test_representations_agree_on_interface_capacity(footprint_network): + """The transfer capacity a zone can deliver is the same in both modes.""" + store = footprint_network.copy() + add_external_regions(store, "imports", "store", FLOWGATES, 30.0, zone_col="reeds_zone") + generator = footprint_network.copy() + add_external_regions(generator, "imports", "generator", FLOWGATES, 30.0, zone_col="reeds_zone") + + pd.testing.assert_series_equal( + import_links(store).set_index("bus0")["p_nom"].sort_index(), + import_links(generator).set_index("bus0")["p_nom"].sort_index(), + ) + + +def test_generator_mode_moves_emissions_onto_the_generator_carrier(footprint_network): + n = footprint_network + add_external_regions(n, "imports", "generator", FLOWGATES, 30.0, co2_emissions=0.428, zone_col="reeds_zone") + + # links never carry emissions in PyPSA, so the factor must sit on the generator + assert n.carriers.at[GENERIC_IMPORT_CARRIER, "co2_emissions"] == pytest.approx(0.428) + assert n.carriers.at["imports", "co2_emissions"] == 0 + + +def test_store_mode_keeps_emissions_on_the_imports_carrier(footprint_network): + n = footprint_network + add_external_regions(n, "imports", "store", FLOWGATES, 30.0, co2_emissions=0.428, zone_col="reeds_zone") + assert n.carriers.at["imports", "co2_emissions"] == pytest.approx(0.428) + + +def test_generator_price_comes_from_the_wholesale_table(footprint_network): + n = footprint_network + add_external_regions(n, "imports", "generator", FLOWGATES, wholesale_costs(), zone_col="reeds_zone") + + name = f"p13_imports {GENERIC_IMPORT_CARRIER}" + assert name in n.generators_t.marginal_cost.columns + assert n.generators_t.marginal_cost[name].to_numpy() == pytest.approx(42.0) + + +@pytest.mark.parametrize("representation", ["store", "generator"]) +def test_exports_are_identical_across_representations(footprint_network, representation): + """Export construction is deliberately untouched by the representation switch.""" + n = footprint_network + add_external_regions(n, "exports", representation, FLOWGATES, -30.0, zone_col="reeds_zone") + + links = n.links[n.links.carrier == "exports"] + assert sorted(links.index) == ["p9_p13_exports"] + assert links.at["p9_p13_exports", "bus0"] == "p9" + assert links.at["p9_p13_exports", "bus1"] == "p13_exports" + assert links.at["p9_p13_exports", "p_nom"] == 1500.0 + # the export revenue stays on the link in both modes + assert links.at["p9_p13_exports", "marginal_cost"] == -30.0 + assert not n.stores[n.stores.carrier == "exports"].empty + # nothing but export links may inject into the export bus + assert n.generators[n.generators.bus == "p13_exports"].empty + + +# --------------------------------------------------------------------------- +# external_regions: contracted units behind the boundary +# --------------------------------------------------------------------------- + + +def test_map_remote_units_to_zones_prefers_the_largest_matching_zone(): + inbound = pd.Series({"p12": 500.0, "p13": 2000.0, "p28": 1000.0}) + zones = map_remote_units_to_zones( + pd.Series({"R Hoover": "NV", "R Palo Verde": "AZ"}), + ["p12", "p13", "p28"], + inbound, + "reeds_zone", + MEMBERSHIP, + ) + assert zones["R Hoover"] == "p13" # NV, and p13 > p12 + assert zones["R Palo Verde"] == "p28" + + +def test_map_remote_units_to_zones_falls_back_and_warns(caplog): + inbound = pd.Series({"p12": 500.0, "p13": 2000.0, "p28": 1000.0}) + with caplog.at_level(logging.WARNING): + zones = map_remote_units_to_zones( + pd.Series({"R Intermountain": "UT"}), + ["p12", "p13", "p28"], + inbound, + "reeds_zone", + MEMBERSHIP, + ) + assert zones["R Intermountain"] == "p13" # largest inbound capacity overall + assert "R Intermountain" in caplog.text + assert "UT" in caplog.text + + +def test_map_remote_units_to_zones_state_boundaries_need_no_membership(): + inbound = pd.Series({"NV": 2000.0, "AZ": 1000.0}) + zones = map_remote_units_to_zones( + pd.Series({"R Apex": "NV"}), + ["NV", "AZ"], + inbound, + "reeds_state", + None, + ) + assert zones["R Apex"] == "NV" + + +def remote_bundle(n): + """A toy `build_remote_unit_bundle` output: one AZ firm unit, one UT VRE unit.""" + units = pd.DataFrame( + [ + { + "name": "R Palo Verde", + "carrier": "nuclear", + "bus": "p9", + "p_nom": 600.0, + "state": "AZ", + "efficiency": 0.33, + "marginal_cost": 8.0, + "heat_rate": 10.0, + "summer_derate": 1.0, + "winter_derate": 1.0, + "ramp_limit_up": 1.0, + "ramp_limit_down": 1.0, + "min_up_time": 0, + "min_down_time": 0, + "start_up_cost": 0.0, + "fuel_cost": 1.0, + "min_load_pu": 0.0, + "build_year": 1988, + "duration": 4.0, + }, + { + "name": "R Cape Solar", + "carrier": "solar", + "bus": "p9", + "p_nom": 50.0, + "state": "UT", + "efficiency": 1.0, + "marginal_cost": 0.0, + "heat_rate": 0.0, + "summer_derate": 1.0, + "winter_derate": 1.0, + "ramp_limit_up": 1.0, + "ramp_limit_down": 1.0, + "min_up_time": 0, + "min_down_time": 0, + "start_up_cost": 0.0, + "fuel_cost": 0.0, + "min_load_pu": 0.0, + "build_year": 2024, + "duration": 4.0, + }, + ], + ).set_index("name") + return { + "units": units, + "vre_profiles": pd.DataFrame({"R Cape Solar": np.full(len(n.snapshots), 0.4)}, index=n.snapshots), + "costs": pd.DataFrame(), + "conventional_carriers": ["nuclear"], + "unit_commitment": False, + } + + +def test_remote_bundle_attaches_behind_the_matching_external_bus(footprint_network, caplog): + n = footprint_network + with caplog.at_level(logging.WARNING): + add_external_regions( + n, + "imports", + "generator", + FLOWGATES, + 30.0, + co2_emissions=0.428, + zone_col="reeds_zone", + remote_bundle=remote_bundle(n), + membership=MEMBERSHIP, + ) + + # AZ unit lands behind the AZ boundary zone + assert n.generators.at["R Palo Verde", "bus"] == "p28_imports" + assert n.generators.at["R Palo Verde", "p_nom"] == 600.0 + assert not n.generators.at["R Palo Verde", "p_nom_extendable"] + + # UT has no direct CA interface -> fallback to the largest inbound zone + assert n.generators.at["R Cape Solar", "bus"] == "p13_imports" + assert "R Cape Solar" in caplog.text + + # the borrowed profile travels with the bundle + assert n.generators_t.p_max_pu["R Cape Solar"].to_numpy() == pytest.approx(0.4) + + # and none of them sit inside the footprint any more + assert "p9" not in set(n.generators.loc[["R Palo Verde", "R Cape Solar"], "bus"]) + + +def test_remote_bundle_is_ignored_in_store_mode(footprint_network): + n = footprint_network + add_external_regions( + n, + "imports", + "store", + FLOWGATES, + 30.0, + zone_col="reeds_zone", + remote_bundle=remote_bundle(n), + membership=MEMBERSHIP, + ) + assert "R Palo Verde" not in n.generators.index + + +def test_unknown_representation_raises(footprint_network): + with pytest.raises(ValueError, match="representation"): + add_external_regions(footprint_network, "imports", "banana", FLOWGATES, 30.0, zone_col="reeds_zone") From 7d7f0d86862007477e52d4e100e1c0d59c83c2b1 Mon Sep 17 00:00:00 2001 From: Kamran Date: Tue, 1 Sep 2026 08:39:34 -0700 Subject: [PATCH 09/13] fix: resolve county-FIPS external zones to states for remote-unit placement At county resolution the boundary zones are 'p'+FIPS ids, absent from the ReEDS BA membership table, so every contracted unit fell through to the largest-inbound fallback (all 49 CPUC units behind Clark County). Derive the state from the FIPS prefix instead; BA zones still use membership. Co-Authored-By: Claude Fable 5 --- workflow/scripts/external_regions.py | 79 ++++++++++++++++++++++-- workflow/scripts/test/test_interfaces.py | 15 +++++ 2 files changed, 88 insertions(+), 6 deletions(-) diff --git a/workflow/scripts/external_regions.py b/workflow/scripts/external_regions.py index f05c3f11d..b4afa8aeb 100644 --- a/workflow/scripts/external_regions.py +++ b/workflow/scripts/external_regions.py @@ -70,6 +70,7 @@ """ import logging +import re import dill import pandas as pd @@ -528,6 +529,64 @@ def map_remote_units_to_zones( return pd.Series(assignment, dtype=object) +# County zones are "p" + 5-digit county FIPS; the first two digits are the state. +_STATE_BY_FIPS = { + "01": "AL", + "02": "AK", + "04": "AZ", + "05": "AR", + "06": "CA", + "08": "CO", + "09": "CT", + "10": "DE", + "11": "DC", + "12": "FL", + "13": "GA", + "15": "HI", + "16": "ID", + "17": "IL", + "18": "IN", + "19": "IA", + "20": "KS", + "21": "KY", + "22": "LA", + "23": "ME", + "24": "MD", + "25": "MA", + "26": "MI", + "27": "MN", + "28": "MS", + "29": "MO", + "30": "MT", + "31": "NE", + "32": "NV", + "33": "NH", + "34": "NJ", + "35": "NM", + "36": "NY", + "37": "NC", + "38": "ND", + "39": "OH", + "40": "OK", + "41": "OR", + "42": "PA", + "44": "RI", + "45": "SC", + "46": "SD", + "47": "TN", + "48": "TX", + "49": "UT", + "50": "VT", + "51": "VA", + "53": "WA", + "54": "WV", + "55": "WI", + "56": "WY", +} + +_COUNTY_ZONE_RE = re.compile(r"^p?(\d{2})\d{3}$") + + def _external_zone_states( external_zones: list[str], zone_col: str, @@ -535,16 +594,24 @@ def _external_zone_states( ) -> dict[str, str]: """State code of each external zone. - With ``topological_boundaries: state`` the zone IS the state code; otherwise - the ReEDS membership table maps the balancing area onto its state. + With ``topological_boundaries: state`` the zone IS the state code; county + zones ("p" + county FIPS) resolve through the state FIPS prefix; ReEDS + balancing areas resolve through the membership table. """ if zone_col == "reeds_state": return {zone: zone for zone in external_zones} - if membership is None: + mbshp = membership.set_index("ba")["st"] if membership is not None else pd.Series(dtype=object) + + def state_of(zone: str) -> str | None: + m = _COUNTY_ZONE_RE.match(zone) + if m and zone not in mbshp.index: + return _STATE_BY_FIPS.get(m.group(1)) + return mbshp.get(zone) + + states = {zone: state_of(zone) for zone in external_zones} + if membership is None and not all(states.values()): logger.warning("No ReEDS membership table supplied; cannot map external zones onto states.") - return {} - mbshp = membership.set_index("ba")["st"] - return {zone: mbshp.get(zone) for zone in external_zones} + return states def attach_remote_contracted_units_externally( diff --git a/workflow/scripts/test/test_interfaces.py b/workflow/scripts/test/test_interfaces.py index e72655d42..6035afbd7 100644 --- a/workflow/scripts/test/test_interfaces.py +++ b/workflow/scripts/test/test_interfaces.py @@ -507,6 +507,21 @@ def test_map_remote_units_to_zones_falls_back_and_warns(caplog): assert "UT" in caplog.text +def test_map_remote_units_to_zones_resolves_county_fips_states(): + # County zones are "p" + county FIPS and are absent from the BA membership + # table; the state must come from the FIPS prefix (32=NV, 04=AZ). + inbound = pd.Series({"p32003": 3000.0, "p04012": 1000.0, "p04027": 2000.0}) + zones = map_remote_units_to_zones( + pd.Series({"R Hoover": "NV", "R Palo Verde": "AZ"}), + ["p32003", "p04012", "p04027"], + inbound, + "county", + MEMBERSHIP, + ) + assert zones["R Hoover"] == "p32003" + assert zones["R Palo Verde"] == "p04027" # AZ, and p04027 > p04012 + + def test_map_remote_units_to_zones_state_boundaries_need_no_membership(): inbound = pd.Series({"NV": 2000.0, "AZ": 1000.0}) zones = map_remote_units_to_zones( From 6983cf1bcb579446c49c4faaced174bfcf7c3fe1 Mon Sep 17 00:00:00 2001 From: Kamran Date: Tue, 1 Sep 2026 09:46:49 -0700 Subject: [PATCH 10/13] feat: allow flat-number and carrier-name import/export costs in the schema; CA imports at 60 USD/MWh The costs code path has always accepted a float or a carrier name, but the schema enum only admitted the literal strings 'wholesale'/'carrier'/'float', so neither actually validated. California switches off the mislabeled 'wholesale' (retail-priced, #807) EIA series to a flat 60 USD/MWh. Co-Authored-By: Claude Fable 5 --- workflow/repo_data/config/config.california.yaml | 1 + workflow/schemas/config.schema.yaml | 12 ++++++++++-- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/workflow/repo_data/config/config.california.yaml b/workflow/repo_data/config/config.california.yaml index f3faa1d60..594e6a2bb 100644 --- a/workflow/repo_data/config/config.california.yaml +++ b/workflow/repo_data/config/config.california.yaml @@ -56,6 +56,7 @@ electricity: imports: enable: true representation: generator # CPUC out-of-state contracts sit BEHIND the WECC boundary, so their deliveries are metered against the import cap + costs: 60 # flat $/MWh for unspecified imports; the 'wholesale' EIA series is retail-priced (see issue #807) volume_limit: 25 # % of annual CA load balancing_period: year exports: diff --git a/workflow/schemas/config.schema.yaml b/workflow/schemas/config.schema.yaml index 0368f2221..0729b3871 100644 --- a/workflow/schemas/config.schema.yaml +++ b/workflow/schemas/config.schema.yaml @@ -285,7 +285,11 @@ properties: properties: enable: {type: boolean} representation: {enum: [store, generator]} - costs: {enum: [wholesale, carrier, float]} + costs: + description: > + 'wholesale' (EIA price series), a carrier name (that fleet's + average marginal cost), or a flat number in $/MWh. + type: [string, number] co2_emissions: {type: number} capacity_limit: {type: boolean} volume_limit: {} @@ -296,7 +300,11 @@ properties: additionalProperties: false properties: enable: {type: boolean} - costs: {enum: [wholesale, carrier, float]} + costs: + description: > + 'wholesale' (EIA price series), a carrier name (that fleet's + average marginal cost), or a flat number in $/MWh. + type: [string, number] capacity_limit: {type: boolean} volume_limit: {} balancing_period: {enum: [day, week, month, year]} From 8d02897c5ba9043d4c175191b711443ada1b22b6 Mon Sep 17 00:00:00 2001 From: Kamran Date: Tue, 1 Sep 2026 14:48:20 -0700 Subject: [PATCH 11/13] fix: clamp EIA coal price query to >=2008 for early weather years EIA coal/shipments/receipts has no data before 2008 (2007 returns only spillover from the inclusive end bound). Snapshot years <2008 crashed build_fuel_prices with AttributeError on the empty payload; 2007 would have silently produced all-NaN coal prices. Fetch 2008 prices instead and relabel the month-start index onto the snapshot year. Co-Authored-By: Claude Fable 5 --- workflow/scripts/build_fuel_prices.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/workflow/scripts/build_fuel_prices.py b/workflow/scripts/build_fuel_prices.py index dd29aa58a..5a3d0adc2 100644 --- a/workflow/scripts/build_fuel_prices.py +++ b/workflow/scripts/build_fuel_prices.py @@ -73,12 +73,24 @@ def get_state_ng_power_prices(sns: pd.date_range, eia_api: str) -> pd.DataFrame: def get_state_coal_power_prices(sns: pd.date_range, eia_api: str) -> pd.DataFrame: + year = sns.year[0] + # EIA coal shipment receipts (coal/shipments/receipts) have no data before + # 2008; earlier years return an empty payload. Use 2008 prices for earlier + # snapshot years, relabelled onto the snapshot year (index is month starts, + # so no Feb-29 hazard). + data_year = max(year, 2008) + if data_year != year: + logger.warning( + f"No EIA coal prices before 2008; using {data_year} prices for snapshot year {year}", + ) eia_coal = ( - eia.FuelCosts("coal", sns.year[0], eia_api, industry="power").get_data( + eia.FuelCosts("coal", data_year, eia_api, industry="power").get_data( pivot=True, ) * const.COAL_dol_ton_2_MWHthermal ) + if data_year != year: + eia_coal.index = eia_coal.index.map(lambda ts: ts.replace(year=year)) return make_hourly(eia_coal) From c86683bce370f1aa483c5314b6bf47443540cf75 Mon Sep 17 00:00:00 2001 From: Kamran Date: Tue, 1 Sep 2026 14:59:48 -0700 Subject: [PATCH 12/13] fix: clamp EIA gas price query to >=2002, refactor shared clamp helper Same failure class as the coal clamp: EIA electric-power gas prices begin 2002-01, so weather years 2000-2001 got an empty (or spillover-only) payload -- wy2000 crashed in format_data and wy2001 would have KeyError'd at the snapshot filter. Shift years rather than replace to keep the inclusive-end spillover January collision-free. Co-Authored-By: Claude Fable 5 --- workflow/scripts/build_fuel_prices.py | 57 +++++++++++++++------------ 1 file changed, 32 insertions(+), 25 deletions(-) diff --git a/workflow/scripts/build_fuel_prices.py b/workflow/scripts/build_fuel_prices.py index 5a3d0adc2..482ffe094 100644 --- a/workflow/scripts/build_fuel_prices.py +++ b/workflow/scripts/build_fuel_prices.py @@ -61,36 +61,43 @@ def make_hourly(df: pd.DataFrame) -> pd.DataFrame: ### -def get_state_ng_power_prices(sns: pd.date_range, eia_api: str) -> pd.DataFrame: - df = ( - eia.FuelCosts("gas", sns.year[0], eia_api, industry="power").get_data( - pivot=True, - ) - * 1000 - / const.NG_MWH_2_MMCF - ) # $/MCF -> $/MWh - return make_hourly(df) - - -def get_state_coal_power_prices(sns: pd.date_range, eia_api: str) -> pd.DataFrame: +def _clamped_fuel_costs( + fuel: str, + sns: pd.date_range, + eia_api: str, + floor: int, +) -> pd.DataFrame: + """ + Fetch EIA fuel costs, clamping the query year to the earliest with data. + + EIA electric-power gas prices begin 2002-01 and coal shipment receipts + begin 2008; earlier years return empty payloads that crash format_data. + For earlier snapshot years fetch the floor year and shift the month-start + index back onto the snapshot year (a shift, not a replace, because the + inclusive API end bound leaves a spillover January of the next year). + """ year = sns.year[0] - # EIA coal shipment receipts (coal/shipments/receipts) have no data before - # 2008; earlier years return an empty payload. Use 2008 prices for earlier - # snapshot years, relabelled onto the snapshot year (index is month starts, - # so no Feb-29 hazard). - data_year = max(year, 2008) + data_year = max(year, floor) + df = eia.FuelCosts(fuel, data_year, eia_api, industry="power").get_data( + pivot=True, + ) if data_year != year: logger.warning( - f"No EIA coal prices before 2008; using {data_year} prices for snapshot year {year}", + f"No EIA {fuel} power prices before {floor}; using {data_year} prices for snapshot year {year}", ) - eia_coal = ( - eia.FuelCosts("coal", data_year, eia_api, industry="power").get_data( - pivot=True, + df.index = df.index.map( + lambda ts: ts.replace(year=ts.year - (data_year - year)), ) - * const.COAL_dol_ton_2_MWHthermal - ) - if data_year != year: - eia_coal.index = eia_coal.index.map(lambda ts: ts.replace(year=year)) + return df + + +def get_state_ng_power_prices(sns: pd.date_range, eia_api: str) -> pd.DataFrame: + df = _clamped_fuel_costs("gas", sns, eia_api, floor=2002) * 1000 / const.NG_MWH_2_MMCF # $/MCF -> $/MWh + return make_hourly(df) + + +def get_state_coal_power_prices(sns: pd.date_range, eia_api: str) -> pd.DataFrame: + eia_coal = _clamped_fuel_costs("coal", sns, eia_api, floor=2008) * const.COAL_dol_ton_2_MWHthermal return make_hourly(eia_coal) From 041c69dcae003606c8a20cb200276193aaada83f Mon Sep 17 00:00:00 2001 From: Kamran Date: Tue, 1 Sep 2026 17:36:11 -0700 Subject: [PATCH 13/13] fix: clamp EIA electricity price query to >=2001 for wy2000 exports pricing Third EIA data floor: retail-sales electricity prices begin 2001, so the exports 'wholesale' path in add_extra_components crashed for weather year 2000 on an empty payload. No index relabel needed -- _build_cost_timeseries already maps the data onto investment periods. Co-Authored-By: Claude Fable 5 --- workflow/scripts/external_regions.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/workflow/scripts/external_regions.py b/workflow/scripts/external_regions.py index b4afa8aeb..6038258db 100644 --- a/workflow/scripts/external_regions.py +++ b/workflow/scripts/external_regions.py @@ -148,7 +148,16 @@ def calc_import_export_costs(n: pypsa.Network, carrier: str) -> float: def load_import_export_costs(eia_api: str, year: int) -> pd.DataFrame: """Loads fuel costs from EIA.""" - return FuelCosts(fuel="electricity", year=year, api=eia_api).get_data() + # EIA retail-sales electricity prices begin 2001; earlier years return an + # empty payload that crashes format_data. No date shifting is needed: the + # downstream _build_cost_timeseries relabels the index onto the network's + # investment periods anyway. + data_year = max(year, 2001) + if data_year != year: + logger.warning( + f"No EIA electricity prices before 2001; using {data_year} prices for year {year}", + ) + return FuelCosts(fuel="electricity", year=data_year, api=eia_api).get_data() def format_import_export_costs(n: pypsa.Network, fuel_costs: pd.DataFrame) -> pd.DataFrame: