validation: add structure for price and crossborder comparison

config/test/config.validation.yaml

@@ -0,0 +1,12 @@
+# SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: CC0-1.0
+
+run:
+  name: "validation-test" # use this to keep track of runs with different settings
+
+scenario:
+  clusters:  # number of nodes in Europe, any integer between 37 (1 node per country-zone) and several hundred
+  - 37
+  opts: # only relevant for PyPSA-Eur
+  - 'Ept-12h'

rules/collect.smk

@@ -84,6 +84,7 @@ rule validate_elec_networks:
         ),
         expand(
             RESULTS
-            + "figures/.validation_plots_elec_s{simpl}_{clusters}_ec_l{ll}_{opts}",
+            + "figures/.validation_{kind}_plots_elec_s{simpl}_{clusters}_ec_l{ll}_{opts}",
-            **config["scenario"]
+            **config["scenario"],
+            kind=["production", "prices", "cross_border"]
         ),

rules/common.smk

@@ -16,7 +16,7 @@ def memory(w):
             factor *= int(m.group(1)) / 8760
             break
     if w.clusters.endswith("m") or w.clusters.endswith("c"):
-        return int(factor * (35000 + 180 * int(w.clusters[:-1])))
+        return int(factor * (35000 + 600 * int(w.clusters[:-1])))
     elif w.clusters == "all":
         return int(factor * (18000 + 180 * 4000))
     else:

rules/validate.smk

@@ -2,6 +2,14 @@
 #
 # SPDX-License-Identifier: MIT
 
+PRODUCTION_PLOTS = [
+    "production_bar",
+    "production_deviation_bar",
+    "seasonal_operation_area",
+]
+CROSS_BORDER_PLOTS = []
+PRICES_PLOTS = ["price_bar"]
+
 
 rule build_electricity_production:
     """
@@ -21,10 +29,45 @@ rule build_electricity_production:
         "../scripts/build_electricity_production.py"
 
 
-PLOTS = ["production_bar", "production_deviation_bar", "seasonal_operation_area"]
+rule build_cross_border_flows:
+    """
+    This rule builds the cross-border flows from ENTSO-E data.
+    The data is used for validation of the optimization results.
+    """
+    params:
+        snapshots=config["snapshots"],
+        countries=config["countries"],
+    input:
+        network=RESOURCES + "networks/base.nc",
+    output:
+        RESOURCES + "historical_cross_border_flows.csv",
+    log:
+        LOGS + "build_cross_border_flows.log",
+    resources:
+        mem_mb=5000,
+    script:
+        "../scripts/build_cross_border_flows.py"
 
 
-rule plot_electricity_production:
+rule build_electricity_prices:
+    """
+    This rule builds the electricity prices from ENTSO-E data.
+    The data is used for validation of the optimization results.
+    """
+    params:
+        snapshots=config["snapshots"],
+        countries=config["countries"],
+    output:
+        RESOURCES + "historical_electricity_prices.csv",
+    log:
+        LOGS + "build_electricity_prices.log",
+    resources:
+        mem_mb=5000,
+    script:
+        "../scripts/build_electricity_prices.py"
+
+
+rule plot_validation_electricity_production:
     input:
         network=RESULTS + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc",
         electricity_production=RESOURCES + "historical_electricity_production.csv",
@@ -32,9 +75,41 @@ rule plot_electricity_production:
         **{
             plot: RESULTS
             + f"figures/validation_{plot}_elec_s{{simpl}}_{{clusters}}_ec_l{{ll}}_{{opts}}.pdf"
-            for plot in PLOTS
+            for plot in PRODUCTION_PLOTS
         },
         plots_touch=RESULTS
-        + "figures/.validation_plots_elec_s{simpl}_{clusters}_ec_l{ll}_{opts}",
+        + "figures/.validation_production_plots_elec_s{simpl}_{clusters}_ec_l{ll}_{opts}",
     script:
-        "../scripts/plot_electricity_production.py"
+        "../scripts/plot_validation_electricity_production.py"
+
+
+rule plot_validation_cross_border_flows:
+    input:
+        network=RESULTS + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc",
+        cross_border_flows=RESOURCES + "historical_cross_border_flows.csv",
+    output:
+        **{
+            plot: RESULTS
+            + f"figures/validation_{plot}_elec_s{{simpl}}_{{clusters}}_ec_l{{ll}}_{{opts}}.pdf"
+            for plot in CROSS_BORDER_PLOTS
+        },
+        plots_touch=RESULTS
+        + "figures/.validation_cross_border_plots_elec_s{simpl}_{clusters}_ec_l{ll}_{opts}",
+    script:
+        "../scripts/plot_validation_cross_border_flows.py"
+
+
+rule plot_validation_electricity_prices:
+    input:
+        network=RESULTS + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc",
+        electricity_prices=RESOURCES + "historical_electricity_prices.csv",
+    output:
+        **{
+            plot: RESULTS
+            + f"figures/validation_{plot}_elec_s{{simpl}}_{{clusters}}_ec_l{{ll}}_{{opts}}.pdf"
+            for plot in PRICES_PLOTS
+        },
+        plots_touch=RESULTS
+        + "figures/.validation_prices_plots_elec_s{simpl}_{clusters}_ec_l{ll}_{opts}",
+    script:
+        "../scripts/plot_validation_electricity_prices.py"

scripts/base_network.py

@@ -337,7 +337,7 @@ def _load_lines_from_eg(buses, eg_lines):
     )
 
     lines["length"] /= 1e3
-
+    lines["carrier"] = "AC"
     lines = _remove_dangling_branches(lines, buses)
 
     return lines

scripts/build_cross_border_flows.py

@@ -0,0 +1,65 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+
+import logging
+
+import pandas as pd
+import pypsa
+from _helpers import configure_logging
+from entsoe import EntsoePandasClient
+from entsoe.exceptions import InvalidBusinessParameterError, NoMatchingDataError
+from requests import HTTPError
+
+logger = logging.getLogger(__name__)
+
+if __name__ == "__main__":
+    if "snakemake" not in globals():
+        from _helpers import mock_snakemake
+
+        snakemake = mock_snakemake("build_cross_border_flows")
+    configure_logging(snakemake)
+
+    api_key = snakemake.config["private"]["keys"]["entsoe_api"]
+    client = EntsoePandasClient(api_key=api_key)
+
+    n = pypsa.Network(snakemake.input.network)
+    start = pd.Timestamp(snakemake.params.snapshots["start"], tz="Europe/Brussels")
+    end = pd.Timestamp(snakemake.params.snapshots["end"], tz="Europe/Brussels")
+
+    branches = n.branches().query("carrier in ['AC', 'DC']")
+    c = n.buses.country
+    branch_countries = pd.concat([branches.bus0.map(c), branches.bus1.map(c)], axis=1)
+    branch_countries = branch_countries.query("bus0 != bus1")
+    branch_countries = branch_countries.apply(sorted, axis=1, result_type="broadcast")
+    country_pairs = branch_countries.drop_duplicates().reset_index(drop=True)
+
+    flows = []
+    unavailable_borders = []
+    for from_country, to_country in country_pairs.values:
+        try:
+            flow_directed = client.query_crossborder_flows(
+                from_country, to_country, start=start, end=end
+            )
+            flow_reverse = client.query_crossborder_flows(
+                to_country, from_country, start=start, end=end
+            )
+            flow = (flow_directed - flow_reverse).rename(
+                f"{from_country} - {to_country}"
+            )
+            flow = flow.tz_localize(None).resample("1h").mean()
+            flow = flow.loc[start.tz_localize(None) : end.tz_localize(None)]
+            flows.append(flow)
+        except (HTTPError, NoMatchingDataError, InvalidBusinessParameterError):
+            unavailable_borders.append(f"{from_country}-{to_country}")
+
+    if unavailable_borders:
+        logger.warning(
+            "Historical electricity cross-border flows for countries"
+            f" {', '.join(unavailable_borders)} not available."
+        )
+
+    flows = pd.concat(flows, axis=1)
+    flows.to_csv(snakemake.output[0])

scripts/build_electricity_prices.py

@@ -0,0 +1,52 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+
+import logging
+
+import pandas as pd
+from _helpers import configure_logging
+from entsoe import EntsoePandasClient
+from entsoe.exceptions import NoMatchingDataError
+
+logger = logging.getLogger(__name__)
+
+if __name__ == "__main__":
+    if "snakemake" not in globals():
+        from _helpers import mock_snakemake
+
+        snakemake = mock_snakemake("build_cross_border_flows")
+    configure_logging(snakemake)
+
+    api_key = snakemake.config["private"]["keys"]["entsoe_api"]
+    client = EntsoePandasClient(api_key=api_key)
+
+    start = pd.Timestamp(snakemake.params.snapshots["start"], tz="Europe/Brussels")
+    end = pd.Timestamp(snakemake.params.snapshots["end"], tz="Europe/Brussels")
+
+    countries = snakemake.params.countries
+
+    prices = []
+    unavailable_countries = []
+
+    for country in countries:
+        country_code = country
+
+        try:
+            gen = client.query_day_ahead_prices(country, start=start, end=end)
+            gen = gen.tz_localize(None).resample("1h").mean()
+            gen = gen.loc[start.tz_localize(None) : end.tz_localize(None)]
+            prices.append(gen)
+        except NoMatchingDataError:
+            unavailable_countries.append(country)
+
+    if unavailable_countries:
+        logger.warning(
+            f"Historical electricity prices for countries {', '.join(unavailable_countries)} not available."
+        )
+
+    keys = [c for c in countries if c not in unavailable_countries]
+    prices = pd.concat(prices, keys=keys, axis=1)
+    prices.to_csv(snakemake.output[0])

scripts/plot_validation_electricity_prices.py

@@ -0,0 +1,57 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+
+import matplotlib.pyplot as plt
+import pandas as pd
+import pypsa
+import seaborn as sns
+from _helpers import configure_logging
+from pypsa.statistics import get_bus_and_carrier
+
+sns.set_theme("paper", style="whitegrid")
+
+if __name__ == "__main__":
+    if "snakemake" not in globals():
+        from _helpers import mock_snakemake
+
+        snakemake = mock_snakemake(
+            "plot_electricity_prices",
+            simpl="",
+            opts="Ept-12h",
+            clusters="37",
+            ll="v1.0",
+        )
+    configure_logging(snakemake)
+
+    n = pypsa.Network(snakemake.input.network)
+    n.loads.carrier = "load"
+
+    historic = pd.read_csv(
+        snakemake.input.electricity_prices,
+        index_col=0,
+        header=[0, 1],
+        parse_dates=True,
+    )
+
+    if len(historic.index) > len(n.snapshots):
+        historic = historic.resample(n.snapshots.inferred_freq).mean().loc[n.snapshots]
+
+    optimized = n.buses_t.marginal_price.groupby(n.buses.country, axis=1).mean()
+
+    data = pd.concat([historic, optimized], keys=["Historic", "Optimized"], axis=1)
+    data.columns.names = ["Kind", "Country"]
+
+    # %% total production per carrier
+    fig, ax = plt.subplots(figsize=(6, 6))
+
+    df = data.mean().unstack().T
+    df.plot.barh(ax=ax, xlabel="Electricity Price [€/MWh]", ylabel="")
+    ax.grid(axis="y")
+    fig.savefig(snakemake.output.price_bar, bbox_inches="tight")
+
+    # touch file
+    with open(snakemake.output.plots_touch, "a"):
+        pass

scripts/plot_validation_electricity_production.py

@@ -45,8 +45,6 @@ if __name__ == "__main__":
         header=[0, 1],
         parse_dates=True,
     )
-    historic = historic.drop("Other renewable", axis=1, level=1)
-    historic = historic.drop("Marine", axis=1, level=1)
 
     colors = n.carriers.set_index("nice_name").color.where(
         lambda s: s != "", "lightgrey"