Include today's district heating share for myopic optimisation

Snakefile

@@ -159,6 +159,7 @@ rule build_energy_totals:
         co2="data/eea/UNFCCC_v23.csv",
         swiss="data/switzerland-sfoe/switzerland-new_format.csv",
         idees="data/jrc-idees-2015",
+        district_heat_share='data/district_heat_share.csv',
         eurostat=input_eurostat
     output:
         energy_name='resources/energy_totals.csv',

config.default.yaml

@@ -141,8 +141,16 @@ existing_capacities:
 
 
 sector:
-  central: true
-  central_fraction: 0.6
+  district_heating:
+    potential: 0.6  # maximum fraction of urban demand which can be supplied by district heating
+     # increase of today's district heating demand to potential maximum district heating share
+     # progress = 0 means today's district heating share, progress = 1 means maximum fraction of urban demand is supplied by district heating
+    progress:
+      2020: 0.0
+      2030: 0.3
+      2040: 0.6
+      2050: 1.0
+    district_heating_loss: 0.15
   bev_dsm_restriction_value: 0.75  #Set to 0 for no restriction on BEV DSM
   bev_dsm_restriction_time: 7  #Time at which SOC of BEV has to be dsm_restriction_value
   transport_heating_deadband_upper: 20.
@@ -151,7 +159,6 @@ sector:
   ICE_upper_degree_factor: 1.6
   EV_lower_degree_factor: 0.98
   EV_upper_degree_factor: 0.63
-  district_heating_loss: 0.15
   bev_dsm: true #turns on EV battery
   bev_availability: 0.5  #How many cars do smart charging
   bev_energy: 0.05  #average battery size in MWh

data/district_heat_share.csv

@@ -0,0 +1,34 @@
+country,share to satisfy heat demand (residential) in percent,capacity[MWth]
+AT,14,11200
+BG,16,6162
+BA,8,
+HR,6.3,2221
+CZ,40,
+DK,65,
+FI,38,23390
+FR,5,
+DE,13.8,
+HU,7.92875588637399,8549
+IS,90,8079000
+IE,0.8,
+IT,3,8727
+LV,73,2254
+LT,56,
+MK,23.7745607009008,636
+NO,4,3400
+PL,42,54912
+PT,0.070754716981132,34
+RS,25,5821
+SI,8.86,1739
+ES,0.251589260787732,1273
+SE,50.4,
+UK,2,
+BY,70,
+EE,52,5406
+KO,3,207
+RO,23,9962
+SK,54,15000
+NL,4,9800
+CH,4,2792
+AL,0,
+ME,0,

doc/data.csv

@@ -25,3 +25,6 @@ Comparative level investment,comparative_level_investment.csv,Eurostat,https://e
 Electricity taxes,electricity_taxes_eu.csv,Eurostat,https://appsso.eurostat.ec.europa.eu/nui/show.do?dataset=nrg_pc_204&lang=en
 Building topologies and corresponding standard values,tabula-calculator-calcsetbuilding.csv,unknown,https://episcope.eu/fileadmin/tabula/public/calc/tabula-calculator.xlsx
 Retrofitting thermal envelope costs for Germany,retro_cost_germany.csv,unkown,https://www.iwu.de/forschung/handlungslogiken/kosten-energierelevanter-bau-und-anlagenteile-bei-modernisierung/
+District heating most countries,jrc-idees-2015/,CC BY 4.0,https://ec.europa.eu/jrc/en/potencia/jrc-idees,,
+District heating missing countries,district_heat_share.csv,unkown,https://www.euroheat.org/knowledge-hub/country-profiles,,
+

scripts/build_energy_totals.py

@@ -212,6 +212,12 @@ def idees_per_country(ct, year):
     assert df.index[47] == "Electricity"
     ct_totals["electricity residential"] = df[47]
 
+    assert df.index[46] == "Derived heat"
+    ct_totals["Derived heat residential"] = df[46]
+
+    assert df.index[50] == 'Thermal uses'
+    ct_totals["Thermal uses residential"] = df[50]
+
     # services
 
     df = pd.read_excel(fn_services, "SER_hh_fec", index_col=0)[year]
@@ -239,6 +245,12 @@ def idees_per_country(ct, year):
     assert df.index[50] == "Electricity"
     ct_totals["electricity services"] = df[50]
 
+    assert df.index[49] == "Derived heat"
+    ct_totals["Derived heat services"] = df[49]
+
+    assert df.index[53] == 'Thermal uses'
+    ct_totals["Thermal uses services"] = df[53]
+
     # transport
 
     df = pd.read_excel(fn_transport, "TrRoad_ene", index_col=0)[year]
@@ -342,6 +354,7 @@ def build_idees(countries, year):
     with mp.Pool(processes=nprocesses) as pool:
         totals_list = list(tqdm(pool.imap(func, countries), **tqdm_kwargs))
 
+
     totals = pd.concat(totals_list, axis=1)
 
     # convert ktoe to TWh
@@ -351,6 +364,13 @@ def build_idees(countries, year):
     # convert TWh/100km to kWh/km
     totals.loc["passenger car efficiency"] *= 10
 
+    # district heating share
+    district_heat = totals.loc[["Derived heat residential",
+                                "Derived heat services"]].sum()
+    total_heat = totals.loc[["Thermal uses residential",
+                             "Thermal uses services"]].sum()
+    totals.loc["district heat share"] = district_heat.div(total_heat)
+
     return totals.T
 
 
@@ -493,7 +513,7 @@ def build_energy_totals(countries, eurostat, swiss, idees):
 
     for purpose in ["passenger", "freight"]:
         attrs = [f"total domestic aviation {purpose}", f"total international aviation {purpose}"]
-        df.loc[missing, f"total aviation {purpose}"] = df.loc[missing, attrs].sum(axis=1)   
+        df.loc[missing, f"total aviation {purpose}"] = df.loc[missing, attrs].sum(axis=1)
 
     if "BA" in df.index:
         # fill missing data for BA (services and road energy data)
@@ -502,6 +522,14 @@ def build_energy_totals(countries, eurostat, swiss, idees):
         ratio = df.at["BA", "total residential"] / df.at["RS", "total residential"]
         df.loc['BA', missing] = ratio * df.loc["RS", missing]
 
+    # Missing district heating share
+    dh_share = pd.read_csv(snakemake.input.district_heat_share,
+                           index_col=0, usecols=[0, 1])
+    # make conservative assumption and take minimum from both data sets
+    df["district heat share"] = (pd.concat([df["district heat share"],
+                                            dh_share.reindex(index=df.index)/100],
+                                           axis=1).min(axis=1))
+
     return df
 
 

scripts/prepare_sector_network.py

@@ -489,8 +489,7 @@ def add_dac(n, costs):
     efficiency3 = -(costs.at['direct air capture', 'heat-input'] - costs.at['direct air capture', 'compression-heat-output'])
 
     n.madd("Link",
-        locations,
-        suffix=" DAC",
+        heat_buses.str.replace(" heat", " DAC"),
         bus0="co2 atmosphere",
         bus1=spatial.co2.df.loc[locations, "nodes"].values,
         bus2=locations.values,
@@ -636,6 +635,9 @@ def prepare_data(n):
 
     nodal_energy_totals = energy_totals.loc[pop_layout.ct].fillna(0.)
     nodal_energy_totals.index = pop_layout.index
+    # district heat share not weighted by population
+    district_heat_share = round(nodal_energy_totals["district heat share"],
+                                ndigits=2)
     nodal_energy_totals = nodal_energy_totals.multiply(pop_layout.fraction, axis=0)
 
     # copy forward the daily average heat demand into each hour, so it can be multipled by the intraday profile
@@ -758,7 +760,7 @@ def prepare_data(n):
     )
 
 
-    return nodal_energy_totals, heat_demand, ashp_cop, gshp_cop, solar_thermal, transport, avail_profile, dsm_profile, nodal_transport_data
+    return nodal_energy_totals, heat_demand, ashp_cop, gshp_cop, solar_thermal, transport, avail_profile, dsm_profile, nodal_transport_data, district_heat_share
 
 
 # TODO checkout PyPSA-Eur script
@@ -1336,11 +1338,10 @@ def add_heat(n, costs):
 
     sectors = ["residential", "services"]
 
-    nodes = create_nodes_for_heat_sector()
 
-    #NB: must add costs of central heating afterwards (EUR 400 / kWpeak, 50a, 1% FOM from Fraunhofer ISE)
+    nodes, dist_fraction, urban_fraction = create_nodes_for_heat_sector()
 
-    urban_fraction = options['central_fraction'] * pop_layout["urban"] / pop_layout[["urban", "rural"]].sum(axis=1)
+    #NB: must add costs of central heating afterwards (EUR 400 / kWpeak, 50a, 1% FOM from Fraunhofer ISE)
 
     # exogenously reduce space heat demand
     if options["reduce_space_heat_exogenously"]:
@@ -1372,15 +1373,22 @@ def add_heat(n, costs):
         ## Add heat load
 
         for sector in sectors:
+            # heat demand weighting
             if "rural" in name:
                 factor = 1 - urban_fraction[nodes[name]]
-            elif "urban" in name:
-                factor = urban_fraction[nodes[name]]
+            elif "urban central" in name:
+                factor = dist_fraction[nodes[name]]
+            elif "urban decentral" in name:
+                factor = urban_fraction[nodes[name]] - \
+                    dist_fraction[nodes[name]]
+            else:
+                raise NotImplementedError(f" {name} not in " f"heat systems: {heat_systems}")
+
             if sector in name:
                 heat_load = heat_demand[[sector + " water",sector + " space"]].groupby(level=1,axis=1).sum()[nodes[name]].multiply(factor)
 
         if name == "urban central":
-            heat_load = heat_demand.groupby(level=1,axis=1).sum()[nodes[name]].multiply(urban_fraction[nodes[name]] * (1 + options['district_heating_loss']))
+            heat_load = heat_demand.groupby(level=1,axis=1).sum()[nodes[name]].multiply(factor * (1 + options['district_heating']['district_heating_loss']))
 
         n.madd("Load",
             nodes[name],
@@ -1661,23 +1669,42 @@ def create_nodes_for_heat_sector():
     # urban are areas with high heating density
     # urban can be split into district heating (central) and individual heating (decentral)
 
+    ct_urban = pop_layout.urban.groupby(pop_layout.ct).sum()
+    # distribution of urban population within a country
+    pop_layout["urban_ct_fraction"] = pop_layout.urban / pop_layout.ct.map(ct_urban.get)
+
     sectors = ["residential", "services"]
 
     nodes = {}
+    urban_fraction = pop_layout.urban / pop_layout[["rural", "urban"]].sum(axis=1)
+
     for sector in sectors:
         nodes[sector + " rural"] = pop_layout.index
-
-        if options["central"]:
-            # TODO: this looks hardcoded, move to config
-            urban_decentral_ct = pd.Index(["ES", "GR", "PT", "IT", "BG"])
-            nodes[sector + " urban decentral"] = pop_layout.index[pop_layout.ct.isin(urban_decentral_ct)]
-        else:
-            nodes[sector + " urban decentral"] = pop_layout.index
-
-    # for central nodes, residential and services are aggregated
-    nodes["urban central"] = pop_layout.index.symmetric_difference(nodes["residential urban decentral"])
-
-    return nodes
+        nodes[sector + " urban decentral"] = pop_layout.index
+
+    # maximum potential of urban demand covered by district heating
+    central_fraction = options['district_heating']["potential"]
+    # district heating share at each node
+    dist_fraction_node = district_heat_share * pop_layout["urban_ct_fraction"] / pop_layout["fraction"]
+    nodes["urban central"] = dist_fraction_node.index
+    # if district heating share larger than urban fraction -> set urban
+    # fraction to district heating share
+    urban_fraction = pd.concat([urban_fraction, dist_fraction_node],
+                               axis=1).max(axis=1)
+    # difference of max potential and today's share of district heating
+    diff = (urban_fraction * central_fraction) - dist_fraction_node
+    progress = get(options["district_heating"]["potential"], investment_year)
+    dist_fraction_node += diff * progress
+    print("************************************")
+    print(
+        "the current DH share compared to the maximum possible is increased \
+           \n by a progress factor of ",
+        progress,
+        "resulting DH share: ",
+        dist_fraction_node)
+    print("**********************************")
+
+    return nodes, dist_fraction_node, urban_fraction
 
 
 def add_biomass(n, costs):
@@ -1993,7 +2020,7 @@ def add_industry(n, costs):
 
     if options["oil_boilers"]:
 
-        nodes_heat = create_nodes_for_heat_sector()
+        nodes_heat = create_nodes_for_heat_sector()[0]
 
         for name in ["residential rural", "services rural", "residential urban decentral", "services urban decentral"]:
 
@@ -2191,17 +2218,18 @@ def limit_individual_line_extension(n, maxext):
     hvdc = n.links.index[n.links.carrier == 'DC']
     n.links.loc[hvdc, 'p_nom_max'] = n.links.loc[hvdc, 'p_nom'] + maxext
 
-
+#%%
 if __name__ == "__main__":
     if 'snakemake' not in globals():
         from helper import mock_snakemake
         snakemake = mock_snakemake(
             'prepare_sector_network',
             simpl='',
-            clusters="45",
+            opts="",
+            clusters="37",
             lv=1.0,
-            opts='',
             sector_opts='Co2L0-168H-T-H-B-I-solar3-dist1',
+            planning_horizons="2020",
             planning_horizons="2030",
         )
 
@@ -2254,10 +2282,10 @@ def limit_individual_line_extension(n, maxext):
         if o == "biomasstransport":
             options["biomass_transport"] = True
 
-    nodal_energy_totals, heat_demand, ashp_cop, gshp_cop, solar_thermal, transport, avail_profile, dsm_profile, nodal_transport_data = prepare_data(n)
+    nodal_energy_totals, heat_demand, ashp_cop, gshp_cop, solar_thermal, transport, avail_profile, dsm_profile, nodal_transport_data, district_heat_share = prepare_data(n)
 
     if "nodistrict" in opts:
-        options["central"] = False
+        options["district_heating"]["progress"] = 0.0
 
     if "T" in opts:
         add_land_transport(n, costs)