Merge branch 'main' into refining-scenarios

Author: lindnemi

.gitignore

@@ -31,6 +31,7 @@ doc/_build
 
 config/config.yaml
 config/scenarios.yaml
+config/scenarios.automated.yaml
 
 config.yaml
 config/config.yaml
@@ -82,3 +83,8 @@ zenodo.org
 globalenergymonitor.org/wp-content/uploads/2023/07/Europe-Gas-Tracker-2023-03-v3.xlsx
 fnb-gas.de
 raw.githubusercontent.com
+cutouts
+
+
+# custom local files
+local

Snakefile

@@ -581,6 +581,9 @@ rule plot_ariadne_variables:
         transmission_investment_csv=RESULTS + "ariadne/transmission_investment.csv",
         trassenlaenge_csv=RESULTS + "ariadne/trassenlaenge.csv",
         Kernnetz_Investment_plot=RESULTS + "ariadne/Kernnetz_Investment_plot.png",
+        elec_trade=RESULTS + "ariadne/elec-trade-DE.png",
+        h2_trade=RESULTS + "ariadne/h2-trade-DE.png",
+        trade_balance=RESULTS + "ariadne/trade-balance-DE.png",
     log:
         RESULTS + "logs/plot_ariadne_variables.log",
     script:

config/config.default.yaml

@@ -909,7 +909,7 @@ solving:
     skip_iterations: true
     rolling_horizon: false
     seed: 123
-    custom_extra_functionality: "../data/custom_extra_functionality.py"
+    custom_extra_functionality: "../scripts/pypsa-de/additional_functionality.py"
     # io_api: "direct"  # Increases performance but only supported for the highs and gurobi solvers
     # options that go into the optimize function
     track_iterations: false

rules/solve_myopic.smk

@@ -114,9 +114,6 @@ rule solve_sector_network_myopic:
             "sector", "co2_sequestration_potential", default=200
         ),
         custom_extra_functionality=input_custom_extra_functionality,
-        # custom_extra_functionality=os.path.join(
-        #     os.path.dirname(workflow.snakefile), "scripts/pypsa-de/additional_functionality.py"
-        # ),
         energy_year=config_provider("energy", "energy_totals_year"),
     input:
         network=RESULTS

scripts/pypsa-de/additional_functionality.py

@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
-
 import logging
+import sys
 
 import pandas as pd
 from xarray import DataArray
@@ -51,7 +51,7 @@ def add_capacity_limits(n, investment_year, limits_capacity, sense="maximum"):
 
                 lhs = nom.sum()
 
-                cname = f"capacity_{sense}-{ct}-{c.name}-{carrier.replace(' ','-')}"
+                cname = f"capacity_{sense}-{ct}-{c.name}-{carrier.replace(' ', '-')}"
 
                 if cname in n.global_constraints.index:
                     logger.warning(
@@ -163,7 +163,7 @@ def h2_import_limits(n, investment_year, limits_volume_max):
     for ct in limits_volume_max["h2_import"]:
         limit = limits_volume_max["h2_import"][ct][investment_year] * 1e6
 
-        logger.info(f"limiting H2 imports in {ct} to {limit/1e6} TWh/a")
+        logger.info(f"limiting H2 imports in {ct} to {limit / 1e6} TWh/a")
         pipeline_carrier = [
             "H2 pipeline",
             "H2 pipeline (Kernnetz)",
@@ -234,7 +234,7 @@ def h2_production_limits(n, investment_year, limits_volume_min, limits_volume_ma
         limit_upper = limits_volume_max["electrolysis"][ct][investment_year] * 1e6
 
         logger.info(
-            f"limiting H2 electrolysis in DE between {limit_lower/1e6} and {limit_upper/1e6} TWh/a"
+            f"limiting H2 electrolysis in DE between {limit_lower / 1e6} and {limit_upper / 1e6} TWh/a"
         )
 
         production = n.links[
@@ -283,7 +283,7 @@ def electricity_import_limits(n, investment_year, limits_volume_max):
     for ct in limits_volume_max["electricity_import"]:
         limit = limits_volume_max["electricity_import"][ct][investment_year] * 1e6
 
-        logger.info(f"limiting electricity imports in {ct} to {limit/1e6} TWh/a")
+        logger.info(f"limiting electricity imports in {ct} to {limit / 1e6} TWh/a")
 
         incoming_line = n.lines.index[
             (n.lines.carrier == "AC")
@@ -386,7 +386,7 @@ def add_co2limit_country(n, limit_countries, snakemake, debug=False):
             ]
 
             logger.info(
-                f"For {ct} adding following link carriers to port {port} CO2 constraint: {n.links.loc[links,'carrier'].unique()}"
+                f"For {ct} adding following link carriers to port {port} CO2 constraint: {n.links.loc[links, 'carrier'].unique()}"
             )
 
             if port == "0":
@@ -614,7 +614,7 @@ def add_h2_derivate_limit(n, investment_year, limits_volume_max):
     for ct in limits_volume_max["h2_derivate_import"]:
         limit = limits_volume_max["h2_derivate_import"][ct][investment_year] * 1e6
 
-        logger.info(f"limiting H2 derivate imports in {ct} to {limit/1e6} TWh/a")
+        logger.info(f"limiting H2 derivate imports in {ct} to {limit / 1e6} TWh/a")
 
         incoming = n.links.loc[
             [

scripts/pypsa-de/export_ariadne_variables.py

@@ -2331,11 +2331,49 @@ def get_final_energy(
     # straightforward for the other categories
     # !!! TODO this assert is temporarily disabled because of https://github.com/PyPSA/pypsa-eur/issues/985
 
+    central_heat_withdrawal = (
+        n.statistics.withdrawal(
+            bus_carrier=["urban central heat"],
+            **kwargs,
+        )
+        .filter(
+            like=region,
+        )
+        .groupby("carrier")
+        .sum()
+        .multiply(MWh2PJ)
+    )
+
+    var["Final Energy|Carbon Dioxide Removal|Heat"] = decentral_heat_withdrawal.get(
+        "DAC", 0
+    ) + central_heat_withdrawal.get("DAC", 0)
+
+    electricity = (
+        n.statistics.withdrawal(
+            bus_carrier="AC",
+            **kwargs,
+        )
+        .filter(
+            like=region,
+        )
+        .groupby("carrier")
+        .sum()
+        .multiply(MWh2PJ)
+    )
+
+    var["Final Energy|Carbon Dioxide Removal|Electricity"] = electricity.get("DAC", 0)
+
+    var["Final Energy|Carbon Dioxide Removal"] = (
+        var["Final Energy|Carbon Dioxide Removal|Electricity"]
+        + var["Final Energy|Carbon Dioxide Removal|Heat"]
+    )
+
     var["Final Energy|Electricity"] = (
         var["Final Energy|Agriculture|Electricity"]
         + var["Final Energy|Residential and Commercial|Electricity"]
         + var["Final Energy|Transportation|Electricity"]
         + var["Final Energy|Industry excl Non-Energy Use|Electricity"]
+        + var["Final Energy|Carbon Dioxide Removal|Electricity"]
     )
 
     var["Final Energy|Solids"] = (
@@ -2385,6 +2423,7 @@ def get_final_energy(
         var["Final Energy|Agriculture|Heat"]
         + var["Final Energy|Residential and Commercial|Heat"]
         + var["Final Energy|Industry excl Non-Energy Use|Heat"]
+        + var["Final Energy|Carbon Dioxide Removal|Heat"]
     )
     # var["Final Energy|Solar"] = \
     var["Final Energy|Hydrogen"] = (
@@ -2410,30 +2449,6 @@ def get_final_energy(
 
     var["Final Energy|Waste"] = waste_withdrawal.get("HVC to air", 0)
 
-    var["Final Energy|Carbon Dioxide Removal|Heat"] = decentral_heat_withdrawal.get(
-        "DAC", 0
-    )
-
-    electricity = (
-        n.statistics.withdrawal(
-            bus_carrier="AC",
-            **kwargs,
-        )
-        .filter(
-            like=region,
-        )
-        .groupby("carrier")
-        .sum()
-        .multiply(MWh2PJ)
-    )
-
-    var["Final Energy|Carbon Dioxide Removal|Electricity"] = electricity.get("DAC", 0)
-
-    var["Final Energy|Carbon Dioxide Removal"] = (
-        var["Final Energy|Carbon Dioxide Removal|Electricity"]
-        + var["Final Energy|Carbon Dioxide Removal|Heat"]
-    )
-
     var["Final Energy incl Non-Energy Use incl Bunkers"] = (
         var["Final Energy|Industry"]
         + var["Final Energy|Residential and Commercial"]
@@ -2897,6 +2912,7 @@ def get_emissions(n, region, _energy_totals, industry_demand):
         + var["Emissions|Gross Fossil CO2|Energy|Supply|Hydrogen"]
         + var["Emissions|Gross Fossil CO2|Energy|Supply|Liquids"]
         + var["Emissions|Gross Fossil CO2|Energy|Supply|Solids"]
+        + var["Emissions|Gross Fossil CO2|Energy|Supply|Gases"]
     )
 
     var["Emissions|Gross Fossil CO2|Energy"] = (
@@ -2927,12 +2943,6 @@ def get_emissions(n, region, _energy_totals, industry_demand):
         var["Emissions|CO2|Energy"] + var["Emissions|CO2|Industrial Processes"]
     )
 
-    var["Emissions|Gross Fossil CO2|Energy|Supply"] = (
-        var["Emissions|Gross Fossil CO2|Energy|Supply|Electricity"]
-        + var["Emissions|Gross Fossil CO2|Energy|Supply|Heat"]
-        + var["Emissions|Gross Fossil CO2|Energy|Supply|Hydrogen"]
-    )
-
     emission_difference = var["Emissions|CO2"] - (
         var["Emissions|CO2|Energy and Industrial Processes"]
         + var["Emissions|CO2|Energy|Demand|Bunkers"]
@@ -4174,72 +4184,92 @@ def get_trade(n, region):
     var = pd.Series()
 
     def get_export_import_links(n, region, carriers):
-        exporting = n.links.index[
+        # note: links can also used bidirectional if efficiency=1 (e.g. "H2 pipeline retrofitted")
+        outgoing = n.links.index[
             (n.links.carrier.isin(carriers))
             & (n.links.bus0.str[:2] == region)
             & (n.links.bus1.str[:2] != region)
         ]
 
-        importing = n.links.index[
+        incoming = n.links.index[
             (n.links.carrier.isin(carriers))
             & (n.links.bus0.str[:2] != region)
             & (n.links.bus1.str[:2] == region)
         ]
 
         exporting_p = (
-            n.links_t.p0.loc[:, exporting]
+            # if p0 > 0 (=clip(lower=0)) system is withdrawing from bus0 (DE) and feeding into bus1 (non-DE) -> export
+            n.links_t.p0.loc[:, outgoing]
             .clip(lower=0)
             .multiply(n.snapshot_weightings.generators, axis=0)
             .values.sum()
-            - n.links_t.p0.loc[:, importing]
-            .clip(upper=0)
+            +
+            # if p1 > 0 system is withdrawing from bus1 (DE) and feeding into bus0 (non-DE) -> export
+            n.links_t.p1.loc[:, incoming]
+            .clip(lower=0)
             .multiply(n.snapshot_weightings.generators, axis=0)
             .values.sum()
         )
+
         importing_p = (
-            n.links_t.p0.loc[:, importing]
-            .clip(lower=0)
+            # if p1 < 0 (=clip(upper=0)) system is feeding into bus1 (DE) and withdrawing from bus0 (non-DE) -> import (with negative sign here)
+            n.links_t.p1.loc[:, incoming]
+            .clip(upper=0)
             .multiply(n.snapshot_weightings.generators, axis=0)
             .values.sum()
-            - n.links_t.p0.loc[:, exporting]
+            * -1
+            +
+            # if p0 < 0 (=clip(upper=0)) system is feeding into bus0 (DE) and withdrawing from bus1 (non-DE) -> import (with negative sign here)
+            n.links_t.p0.loc[:, outgoing]
             .clip(upper=0)
             .multiply(n.snapshot_weightings.generators, axis=0)
             .values.sum()
+            * -1
         )
 
         return exporting_p, importing_p
 
     # Trade|Secondary Energy|Electricity|Volume
-    exporting_ac = n.lines.index[
+    outgoing_ac = n.lines.index[
         (n.lines.carrier == "AC")
         & (n.lines.bus0.str[:2] == region)
         & (n.lines.bus1.str[:2] != region)
     ]
 
-    importing_ac = n.lines.index[
+    incoming_ac = n.lines.index[
         (n.lines.carrier == "AC")
         & (n.lines.bus0.str[:2] != region)
         & (n.lines.bus1.str[:2] == region)
     ]
+
     exporting_p_ac = (
-        n.lines_t.p0.loc[:, exporting_ac]
+        # if p0 > 0 (=clip(lower=0)) system is withdrawing from bus0 (DE) and feeding into bus1 (non-DE) -> export
+        n.lines_t.p0.loc[:, outgoing_ac]
         .clip(lower=0)
         .multiply(n.snapshot_weightings.generators, axis=0)
         .values.sum()
-        - n.lines_t.p0.loc[:, importing_ac]
-        .clip(upper=0)
+        +
+        # if p1 > 0 system is withdrawing from bus1 (DE) and feeding into bus0 (non-DE) -> export
+        n.lines_t.p1.loc[:, incoming_ac]
+        .clip(lower=0)
         .multiply(n.snapshot_weightings.generators, axis=0)
         .values.sum()
     )
+
     importing_p_ac = (
-        n.lines_t.p0.loc[:, importing_ac]
-        .clip(lower=0)
+        # if p1 < 0 (=clip(upper=0)) system is feeding into bus1 (DE) and withdrawing from bus0 (non-DE) -> import (with negative sign here)
+        n.lines_t.p1.loc[:, incoming_ac]
+        .clip(upper=0)
         .multiply(n.snapshot_weightings.generators, axis=0)
         .values.sum()
-        - n.lines_t.p0.loc[:, exporting_ac]
+        * -1
+        +
+        # if p0 < 0 (=clip(upper=0)) system is feeding into bus0 (DE) and withdrawing from bus1 (non-DE) -> import (with negative sign here)
+        n.lines_t.p0.loc[:, outgoing_ac]
         .clip(upper=0)
         .multiply(n.snapshot_weightings.generators, axis=0)
         .values.sum()
+        * -1
     )
 
     exports_dc, imports_dc = get_export_import_links(n, region, ["DC"])
@@ -5141,6 +5171,23 @@ def get_data(
     var["Investment|Energy Supply|Hydrogen|Fossil"] = var[
         "Investment|Energy Supply|Hydrogen|Gas"
     ]
+    # For internal use only and translated directly to TWh
+    var["Demand|Electricity"] = (
+        var.reindex(
+            [
+                "Secondary Energy|Electricity|Storage Losses",
+                "Secondary Energy Input|Electricity|Heat",
+                "Secondary Energy Input|Electricity|Hydrogen",
+                "Secondary Energy Input|Electricity|Liquids",
+                "Final Energy|Industry|Electricity",
+                "Final Energy|Agriculture|Electricity",
+                "Final Energy|Residential and Commercial|Electricity",
+                "Final Energy|Transportation|Electricity",
+                "Final Energy|Carbon Dioxide Removal|Electricity",
+            ]
+        ).sum()
+        / 3.6
+    )
 
     data = []
     for v in var.index: