Improve exogenous industry demand

CHANGELOG.md

@@ -1,4 +1,6 @@
 # Changelog
+- Added an option to source industry energy demand from UBA MWMS (Projektionsbericht 2025) for the years 2025-2035
+- renamed some scripts
 - Added an option to source mobility demand from UBA MWMS (Projektionsbericht 2025) for the years 2025-2035
 - Renamed functions and script for exogenous mobility demand
 - Improved the transport demand data, added an option to source 2020 and 2025 data from AGEB instead of Aladin

Snakefile

@@ -555,6 +555,10 @@ rule modify_prenetwork:
         bev_charge_rate=config_provider("sector", "bev_charge_rate"),
         bev_energy=config_provider("sector", "bev_energy"),
         bev_dsm_availability=config_provider("sector", "bev_dsm_availability"),
+        uba_for_industry=config_provider("iiasa_database", "uba_for_industry"),
+        scale_industry_non_energy=config_provider(
+            "iiasa_database", "scale_industry_non_energy"
+        ),
     input:
         costs_modifications="ariadne-data/costs_{planning_horizons}-modifications.csv",
         network=resources(
@@ -575,13 +579,20 @@ rule modify_prenetwork:
         industrial_demand=resources(
             "industrial_energy_demand_base_s_{clusters}_{planning_horizons}.csv"
         ),
+        industrial_production_per_country_tomorrow=resources(
+            "industrial_production_per_country_tomorrow_{planning_horizons}-modified.csv"
+        ),
+        industry_sector_ratios=resources(
+            "industry_sector_ratios_{planning_horizons}.csv"
+        ),
         pop_weighted_energy_totals=resources(
             "pop_weighted_energy_totals_s_{clusters}.csv"
         ),
         shipping_demand=resources("shipping_demand_s_{clusters}.csv"),
         regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"),
         regions_offshore=resources("regions_offshore_base_s_{clusters}.geojson"),
         offshore_connection_points="ariadne-data/offshore_connection_points.csv",
+        new_industrial_energy_demand="ariadne-data/UBA_Projektionsbericht2025_Abbildung31_MWMS.csv",
     output:
         network=resources(
             "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}_final.nc"
@@ -595,7 +606,7 @@ rule modify_prenetwork:
         "scripts/pypsa-de/modify_prenetwork.py"
 
 
-ruleorder: modify_industry_demand > build_industrial_production_per_country_tomorrow
+ruleorder: modify_industry_production > build_industrial_production_per_country_tomorrow
 
 
 rule modify_existing_heating:
@@ -656,7 +667,7 @@ rule build_existing_chp_de:
         "scripts/pypsa-de/build_existing_chp_de.py"
 
 
-rule modify_industry_demand:
+rule modify_industry_production:
     params:
         reference_scenario=config_provider("iiasa_database", "reference_scenario"),
     input:
@@ -671,9 +682,9 @@ rule modify_industry_demand:
     resources:
         mem_mb=1000,
     log:
-        logs("modify_industry_demand_{planning_horizons}.log"),
+        logs("modify_industry_production_{planning_horizons}.log"),
     script:
-        "scripts/pypsa-de/modify_industry_demand.py"
+        "scripts/pypsa-de/modify_industry_production.py"
 
 
 rule build_wasserstoff_kernnetz:

ariadne-data/UBA_Projektionsbericht2025_Abbildung31_MWMS.csv

@@ -0,0 +1,6 @@
+carrier,2025,2030,2035
+fossil,324,258,191
+industry electricity,211,234,249
+solid biomass for industry,31,35,31
+H2 for industry,0,6,42
+low-temperature heat for industry,48,59,63

config/config.de.yaml

@@ -4,7 +4,8 @@
 
 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#run
 run:
-  prefix: 20250807_merge_july
+  prefix: 20250716_improve_industry_demand
+
   name:
   # - ExPol
   - KN2045_Mix
@@ -45,6 +46,8 @@ iiasa_database:
   region: Deutschland
   ageb_for_mobility: true # In 2020 use AGEB data for final energy demand and KBA for vehicles
   uba_for_mobility: false # For 2025–2035 use MWMS scenario from UBA Projektionsbericht 2025
+  uba_for_industry: false # For 2025–2035 use MWMS scenario from UBA Projektionsbericht 2025
+  scale_industry_non_energy: false # Scale non-energy industry demand directly proportional to energy demand
 
 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#foresight
 foresight: myopic

scripts/pypsa-de/export_ariadne_variables.py

@@ -1818,9 +1818,17 @@ def get_secondary_energy(n, region, _industry_demand):
         axis=0,
     ).sum()
     mwh_coal_per_mwh_coke = 1.366
+    coke_fraction = (
+        industry_demand.get("coke")
+        * mwh_coal_per_mwh_coke
+        / (
+            industry_demand.get("coke") * mwh_coal_per_mwh_coke
+            + industry_demand.get("coal")
+        )
+    )
     # Coke is added as a coal demand, so we need to convert back to units of coke for secondary energy
     var["Secondary Energy|Solids|Coal"] = var["Secondary Energy|Solids"] = (
-        industry_demand.get("coke", 0) / mwh_coal_per_mwh_coke
+        sum_load(n, "coal for industry", region) * coke_fraction / mwh_coal_per_mwh_coke
     )
 
     biomass_usage = (
@@ -1990,14 +1998,17 @@ def get_final_energy(
     # !: Pypsa-eur does not strictly distinguish between energy and
     # non-energy use
 
-    var["Final Energy|Industry|Electricity"] = industry_demand.get("electricity")
-    # or use: sum_load(n, "industry electricity", region)
+    var["Final Energy|Industry|Electricity"] = sum_load(
+        n, "industry electricity", region
+    )
     # electricity is not used for non-energy purposes
     var["Final Energy|Industry excl Non-Energy Use|Electricity"] = var[
         "Final Energy|Industry|Electricity"
     ]
 
-    var["Final Energy|Industry|Heat"] = industry_demand.get("low-temperature heat")
+    var["Final Energy|Industry|Heat"] = sum_load(
+        n, "low-temperature heat for industry", region
+    )
     # heat is not used for non-energy purposes
     var["Final Energy|Industry excl Non-Energy Use|Heat"] = var[
         "Final Energy|Industry|Heat"
@@ -2009,7 +2020,7 @@ def get_final_energy(
     # var["Final Energy|Industry|Geothermal"] = \
     # Not implemented
 
-    var["Final Energy|Industry|Gases"] = industry_demand.get("methane")
+    var["Final Energy|Industry|Gases"] = sum_load(n, "gas for industry", region)
 
     for gas_type in gas_fractions.index:
         var[f"Final Energy|Industry|Gases|{gas_type}"] = (
@@ -2031,7 +2042,7 @@ def get_final_energy(
     # var["Final Energy|Industry|Power2Heat"] = \
     # Q: misleading description
 
-    var["Final Energy|Industry|Hydrogen"] = industry_demand.get("hydrogen")
+    var["Final Energy|Industry|Hydrogen"] = sum_load(n, "H2 for industry", region)
     # subtract non-energy used hydrogen from total hydrogen demand
     var["Final Energy|Industry excl Non-Energy Use|Hydrogen"] = (
         var["Final Energy|Industry|Hydrogen"]
@@ -2075,16 +2086,29 @@ def get_final_energy(
 
     # var["Final Energy|Industry|Other"] = \
 
-    var["Final Energy|Industry|Solids|Biomass"] = industry_demand.get("solid biomass")
+    var["Final Energy|Industry|Solids|Biomass"] = sum_load(
+        n, "solid biomass for industry", region
+    )
     var["Final Energy|Industry excl Non-Energy Use|Solids|Biomass"] = var[
         "Final Energy|Industry|Solids|Biomass"
     ]
 
     mwh_coal_per_mwh_coke = 1.366
-    # Coke is added as a coal demand, so we need to convert back to units of coke for final energy
+    coke_fraction = (
+        industry_demand.get("coke")
+        * mwh_coal_per_mwh_coke
+        / (
+            industry_demand.get("coke") * mwh_coal_per_mwh_coke
+            + industry_demand.get("coal")
+        )
+    )
+    # Contains coke demand, which is a coal product
+    # Here coke is considered a secondary energy source
     var["Final Energy|Industry|Solids|Coal"] = (
-        industry_demand.get("coal")
-        + industry_demand.get("coke") / mwh_coal_per_mwh_coke
+        sum_load(n, "coal for industry", region) * (1 - coke_fraction)
+        + sum_load(n, "coal for industry", region)
+        * coke_fraction
+        / mwh_coal_per_mwh_coke
     )
     var["Final Energy|Industry excl Non-Energy Use|Solids|Coal"] = var[
         "Final Energy|Industry|Solids|Coal"
@@ -2572,10 +2596,10 @@ def get_final_energy(
     return var * MWh2TWh
 
 
-def get_emissions(n, region, _energy_totals, industry_demand):
+def get_emissions(n, region, _energy_totals, _industry_demand):
     energy_totals = _energy_totals.loc[region[0:2]]
 
-    industry_DE = industry_demand.filter(
+    industry_demand = _industry_demand.filter(
         like=region,
         axis=0,
     ).sum()
@@ -2881,8 +2905,22 @@ def get_emissions(n, region, _energy_totals, industry_demand):
     )  # considered 0 anyways
 
     mwh_coal_per_mwh_coke = 1.366  # from eurostat energy balance
-    # 0.3361 t/MWh, 1e-6 to convert to Mt
-    coking_emissions = industry_DE.coke * (mwh_coal_per_mwh_coke - 1) * 0.3361 * t2Mt
+    coke_fraction = (
+        industry_demand.get("coke")
+        * mwh_coal_per_mwh_coke
+        / (
+            industry_demand.get("coke") * mwh_coal_per_mwh_coke
+            + industry_demand.get("coal")
+        )
+    )
+    # 0.3361 t_CO2/MWh
+    coking_emissions = (
+        sum_load(n, "coal for industry", region)
+        * coke_fraction
+        * (mwh_coal_per_mwh_coke - 1)
+        * 0.3361
+        * t2Mt
+    )
     var["Emissions|Gross Fossil CO2|Energy|Demand|Industry"] = (
         co2_emissions.reindex(
             [

scripts/pypsa-de/modify_industry_demand.py → scripts/pypsa-de/modify_industry_production.py

@@ -23,7 +23,7 @@
 if __name__ == "__main__":
     if "snakemake" not in globals():
         snakemake = mock_snakemake(
-            "modify_industry_demand",
+            "modify_industry_production",
             simpl="",
             clusters=22,
             opts="",

scripts/pypsa-de/modify_prenetwork.py

@@ -1255,6 +1255,120 @@ def scale_capacity(n, scaling):
                 ]
 
 
+def modify_industry_demand(
+    n,
+    year,
+    industry_energy_demand_file,
+    industry_production_file,
+    sector_ratios_file,
+    scale_non_energy=False,
+):
+    logger.info("Modifying industry demand in Germany.")
+
+    industry_production = pd.read_csv(
+        industry_production_file,
+        index_col="kton/a",
+    ).rename_axis("country")
+
+    sector_ratios = pd.read_csv(
+        sector_ratios_file,
+        header=[0, 1],
+        index_col=0,
+    ).rename_axis("carrier")
+
+    new_demand = pd.read_csv(
+        industry_energy_demand_file,
+        index_col=0,
+    )[str(year)].mul(1e6)
+
+    subcategories = ["HVC", "Methanol", "Chlorine", "Ammonia"]
+    carrier = ["hydrogen", "methane", "naphtha"]
+
+    ip = industry_production.loc["DE", subcategories]  # kt/a
+    sr = sector_ratios["DE"].loc[carrier, subcategories]  # MWh/tMaterial
+    _non_energy = sr.multiply(ip).sum(axis=1) * 1e3
+
+    non_energy = pd.Series(
+        {
+            "industry electricity": 0.0,
+            "low-temperature heat for industry": 0.0,
+            "solid biomass for industry": 0.0,
+            "H2 for industry": _non_energy["hydrogen"],
+            "coal for industry": 0.0,
+            "gas for industry": _non_energy["methane"],
+            "naphtha for industry": _non_energy["naphtha"],
+        }
+    )
+
+    _industry_loads = [
+        "solid biomass for industry",
+        "gas for industry",
+        "H2 for industry",
+        "industry methanol",
+        "naphtha for industry",
+        "low-temperature heat for industry",
+        "industry electricity",
+        "coal for industry",
+    ]
+    industry_loads = n.loads.query(
+        f"carrier in {_industry_loads} and bus.str.startswith('DE')"
+    )
+
+    if scale_non_energy:
+        new_demand_without_non_energy = new_demand.sum()
+        pypsa_industry_without_non_energy = (
+            industry_loads.p_set.sum() * 8760 - non_energy.sum()
+        )
+        non_energy_scaling_factor = (
+            new_demand_without_non_energy / pypsa_industry_without_non_energy
+        )
+        logger.info(
+            f"Scaling non-energy use by {non_energy_scaling_factor:.2f} to match UBA data."
+        )
+        non_energy_corrected = non_energy * non_energy_scaling_factor
+    else:
+        non_energy_corrected = non_energy
+
+    for carrier in [
+        "industry electricity",
+        "H2 for industry",
+        "solid biomass for industry",
+        "low-temperature heat for industry",
+    ]:
+        loads_i = n.loads.query(
+            f"carrier == '{carrier}' and bus.str.startswith('DE')"
+        ).index
+        logger.info(
+            f"Total load of {carrier} in DE before scaling: {n.loads.loc[loads_i, 'p_set'].sum() * 8760:.2f} MWh/a"
+        )
+        total_load = industry_loads.p_set.loc[loads_i].sum() * 8760
+        scaling_factor = (
+            new_demand[carrier] + non_energy_corrected[carrier]
+        ) / total_load
+        n.loads.loc[loads_i, "p_set"] *= scaling_factor
+        logger.info(
+            f"Total load of {carrier} in DE after scaling: {n.loads.loc[loads_i, 'p_set'].sum() * 8760:.2f} MWh/a"
+        )
+
+    # Fossil fuels are aggregated in UBA MWMS but have to be scaled separately
+    fossil_loads = industry_loads.query("carrier.str.contains('gas|coal|naphtha')")
+    fossil_totals = (
+        fossil_loads[["p_set", "carrier"]].groupby("carrier").p_set.sum() * 8760
+    )
+    fossil_energy = fossil_totals - non_energy[fossil_totals.index]
+    fossil_energy_corrected = fossil_energy * new_demand["fossil"] / fossil_energy.sum()
+    fossil_totals_corrected = (
+        fossil_energy_corrected + non_energy_corrected[fossil_totals.index]
+    )
+    for carrier in fossil_totals.index:
+        loads_i = fossil_loads.query(
+            f"carrier == '{carrier}' and bus.str.startswith('DE')"
+        ).index
+        n.loads.loc[loads_i, "p_set"] *= (
+            fossil_totals_corrected[carrier] / fossil_totals[carrier]
+        )
+
+
 if __name__ == "__main__":
     if "snakemake" not in globals():
         snakemake = mock_snakemake(
@@ -1337,4 +1451,18 @@ def scale_capacity(n, scaling):
 
     sanitize_custom_columns(n)
 
+    if snakemake.params.uba_for_industry and current_year >= 2025:
+        if current_year >= 2040:
+            logger.error(
+                "The UBA for industry data is only available for 2025, 2030 and 2035. Please check your config."
+            )
+        modify_industry_demand(
+            n,
+            current_year,
+            snakemake.input.new_industrial_energy_demand,
+            snakemake.input.industrial_production_per_country_tomorrow,
+            snakemake.input.industry_sector_ratios,
+            scale_non_energy=snakemake.params.scale_industry_non_energy,
+        )
+
     n.export_to_netcdf(snakemake.output.network)