Include today's district heating share for myopic optimisation

.gitignore

@@ -28,7 +28,7 @@ gurobi.log
 /data/.nfs*
 /data/Industrial_Database.csv
 /data/retro/tabula-calculator-calcsetbuilding.csv
-/data
+/data/nuts*
 *.org
 
 *.nc

Snakefile

@@ -323,6 +323,7 @@ rule prepare_sector_network:
         transport_name='resources/transport_data.csv',
         traffic_data_KFZ = "data/emobility/KFZ__count",
         traffic_data_Pkw = "data/emobility/Pkw__count",
+        dh_share = "data/district_heat_share.csv",
         biomass_potentials='resources/biomass_potentials.csv',
         heat_profile="data/heat_load_profile_BDEW.csv",
         costs=CDIR + "costs_{planning_horizons}.csv",

config.default.yaml

@@ -143,6 +143,12 @@ existing_capacities:
 sector:
   central: true
   central_fraction: 0.6
+  district_heating_increase: true
+  dh_strength:
+    2020: 0   # starting at today's share
+    2030: 0.2
+    2040: 0.5
+    2050: 1    # maximum possible share defined in central fraction
   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.
@@ -335,7 +341,7 @@ solving:
 
 plotting:
   map:
-    boundaries: [-11, 30, 34, 71] 
+    boundaries: [-11, 30, 34, 71]
     color_geomap:
       ocean: white
       land: whitesmoke

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,

scripts/prepare_sector_network.py

@@ -75,7 +75,7 @@ def co2_emissions_year(countries, opts, year):
     co2_emissions = co2_totals.loc[countries, sectors].sum().sum()
 
     # convert MtCO2 to GtCO2
-    co2_emissions *= 0.001  
+    co2_emissions *= 0.001
 
     return co2_emissions
 
@@ -102,17 +102,17 @@ def build_carbon_budget(o, fn):
 
     #emissions at the beginning of the path (last year available 2018)
     e_0 = co2_emissions_year(countries, opts, year=2018)
-    
+
     #emissions in 2019 and 2020 assumed equal to 2018 and substracted
     carbon_budget -= 2 * e_0
-    
+
     planning_horizons = snakemake.config['scenario']['planning_horizons']
     t_0 = planning_horizons[0]
 
     if "be" in o:
 
         # final year in the path
-        t_f = t_0 + (2 * carbon_budget / e_0).round(0) 
+        t_f = t_0 + (2 * carbon_budget / e_0).round(0)
 
         def beta_decay(t):
             cdf_term = (t - t_0) / (t_f - t_0)
@@ -328,7 +328,7 @@ def add_co2_tracking(n, options):
 
 def add_dac(n, costs):
 
-    heat_carriers = ["urban central heat", "services urban decentral heat"]
+    heat_carriers = ["urban central heat"]
     heat_buses = n.buses.index[n.buses.carrier.isin(heat_carriers)]
     locations = n.buses.location[heat_buses]
 
@@ -818,7 +818,7 @@ def insert_gas_distribution_costs(n, costs):
     # TODO options?
 
     f_costs = options['gas_distribution_grid_cost_factor']
-    
+
     print("Inserting gas distribution grid with investment cost factor of", f_costs)
 
     capital_cost = costs.loc['electricity distribution grid']["fixed"] * f_costs
@@ -827,7 +827,7 @@ def insert_gas_distribution_costs(n, costs):
     gas_b = n.links.index[n.links.carrier.str.contains("gas boiler") &
                           (~n.links.carrier.str.contains("urban central"))]
     n.links.loc[gas_b, "capital_cost"] += capital_cost
-    
+
     # micro CHPs
     mchp = n.links.index[n.links.carrier.str.contains("micro gas")]
     n.links.loc[mchp,  "capital_cost"] += capital_cost
@@ -1075,7 +1075,7 @@ def add_land_transport(n, costs):
             suffix=" EV battery",
             carrier="Li ion"
         )
-        
+
         p_set = electric_share * (transport[nodes] + cycling_shift(transport[nodes], 1) + cycling_shift(transport[nodes], 2)) / 3
 
         n.madd("Load",
@@ -1086,8 +1086,8 @@ def add_land_transport(n, costs):
             p_set=p_set
         )
 
-        
-        p_nom = nodal_transport_data["number cars"] * options.get("bev_charge_rate", 0.011) * electric_share  
+
+        p_nom = nodal_transport_data["number cars"] * options.get("bev_charge_rate", 0.011) * electric_share
 
         n.madd("Link",
             nodes,
@@ -1119,7 +1119,7 @@ def add_land_transport(n, costs):
 
     if electric_share > 0 and options["bev_dsm"]:
 
-        e_nom = nodal_transport_data["number cars"] * options.get("bev_energy", 0.05) * options["bev_availability"] * electric_share 
+        e_nom = nodal_transport_data["number cars"] * options.get("bev_energy", 0.05) * options["bev_availability"] * electric_share
 
         n.madd("Store",
             nodes,
@@ -1179,11 +1179,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"]:
@@ -1199,7 +1198,7 @@ def add_heat(n, costs):
         "services urban decentral",
         "urban central"
     ]
-    
+
     for name in heat_systems:
 
         name_type = "central" if name == "urban central" else "decentral"
@@ -1215,15 +1214,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:
+                factor = None
+
             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_loss']))
 
         n.madd("Load",
             nodes[name],
@@ -1281,16 +1287,16 @@ def add_heat(n, costs):
                 p_nom_extendable=True
             )
 
-            
+
             if isinstance(options["tes_tau"], dict):
                 tes_time_constant_days = options["tes_tau"][name_type]
             else:
                 logger.warning("Deprecated: a future version will require you to specify 'tes_tau' ",
                                "for 'decentral' and 'central' separately.")
                 tes_time_constant_days = options["tes_tau"] if name_type == "decentral" else 180.
-            
+
             # conversion from EUR/m^3 to EUR/MWh for 40 K diff and 1.17 kWh/m^3/K
-            capital_cost = costs.at[name_type + ' water tank storage', 'fixed'] / 0.00117 / 40 
+            capital_cost = costs.at[name_type + ' water tank storage', 'fixed'] / 0.00117 / 40
 
             n.madd("Store",
                 nodes[name] + f" {name} water tanks",
@@ -1503,24 +1509,55 @@ def create_nodes_for_heat_sector():
     # rural are areas with low heating density and individual heating
     # 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()
+    pop_layout["urban_ct_fraction"] = pop_layout["urban"] / \
+            pop_layout["ct"].map(ct_urban.get)
+    # todays district heating share per country
+    dist_heat_share_ct = pd.read_csv(snakemake.input.dh_share, index_col=0,
+                                  usecols=[0,1]).dropna()/100
+    dist_heat_share = pop_layout.ct.map(dist_heat_share_ct["district heating share"])
+
+
     sectors = ["residential", "services"]
-    
+
     nodes = {}
+    urban_fraction = pop_layout["urban"] / \
+                     (pop_layout[["urban", "rural"]].sum(axis=1))
+
     for sector in sectors:
         nodes[sector + " rural"] = pop_layout.index
+        nodes[sector + " urban decentral"] = pop_layout.index
+
+    if options["central"] and not options['district_heating_increase']:
+        central_fraction = options['central_fraction']
+        dist_fraction = central_fraction * urban_fraction
+        nodes["urban central"] = dist_fraction.index
+
+    if options['district_heating_increase']:  # take current district heating share
+        dist_fraction = dist_heat_share * \
+            pop_layout["urban_ct_fraction"] / pop_layout["fraction"]
+        nodes["urban central"] = dist_fraction.index
+        # if district heating share larger than urban fraction -> set urban
+        # fraction to district heating share
+        urban_fraction = pd.concat([urban_fraction, dist_fraction],
+                                   axis=1).max(axis=1)
+        diff = urban_fraction - dist_fraction
+        dist_fraction += diff * get(options["dh_strength"], investment_year)
+        print("************************************")
+        print(
+            "the current DH share compared to the maximum possible is increased \
+               \n by a factor of ",
+            get(options["dh_strength"], investment_year),
+            "resulting DH share: ",
+            dist_fraction)
+        print("**********************************")
 
-        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
+    else:
+        dist_fraction = urban_fraction * 0
+        nodes["urban central"] = dist_fraction.index
+
+    return nodes, dist_fraction, urban_fraction
 
 
 def add_biomass(n, costs):
@@ -1754,7 +1791,7 @@ def add_industry(n, costs):
     if shipping_hydrogen_share < 1:
 
         shipping_oil_share = 1 - shipping_hydrogen_share
-        
+
         p_set = shipping_oil_share * nodal_energy_totals.loc[nodes, all_navigation].sum(axis=1) * 1e6 / 8760.
 
         n.madd("Load",
@@ -1764,7 +1801,7 @@ def add_industry(n, costs):
             carrier="shipping oil",
             p_set=p_set
         )
-        
+
         co2 = shipping_oil_share * nodal_energy_totals.loc[nodes, all_navigation].sum().sum() * 1e6 / 8760 * costs.at["oil", "CO2 intensity"]
 
         n.add("Load",
@@ -1783,7 +1820,7 @@ def add_industry(n, costs):
         )
 
     if "EU oil Store" not in n.stores.index:
-        
+
         #could correct to e.g. 0.001 EUR/kWh * annuity and O&M
         n.add("Store",
             "EU oil Store",
@@ -1805,7 +1842,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"]:
 
@@ -1950,7 +1987,7 @@ def add_waste_heat(n):
 
 
 def decentral(n):
-    """Removes the electricity transmission system."""    
+    """Removes the electricity transmission system."""
     n.lines.drop(n.lines.index, inplace=True)
     n.links.drop(n.links.index[n.links.carrier.isin(["DC", "B2B"])], inplace=True)
 
@@ -2002,17 +2039,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=48,
+            opts="",
+            clusters="37",
             lv=1.0,
             sector_opts='Co2L0-168H-T-H-B-I-solar3-dist1',
-            planning_horizons=2020,
+            planning_horizons="2020",
         )
 
     logging.basicConfig(level=snakemake.config['logging_level'])
@@ -2027,7 +2065,7 @@ def limit_individual_line_extension(n, maxext):
     n = pypsa.Network(snakemake.input.network, override_component_attrs=overrides)
 
     pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0)
-    Nyears = n.snapshot_weightings.generators.sum() / 8760
+    Nyears = n.snapshot_weightings.sum() / 8760
 
     costs = prepare_costs(snakemake.input.costs,
                           snakemake.config['costs']['USD2013_to_EUR2013'],
@@ -2038,7 +2076,7 @@ def limit_individual_line_extension(n, maxext):
     patch_electricity_network(n)
 
     if snakemake.config["foresight"] == 'myopic':
-        
+
         add_lifetime_wind_solar(n, costs)
 
         conventional = snakemake.config['existing_capacities']['conventional_carriers']