fully included halloran et al 2023, cleanup heat, bev args

rules/build_electricity.smk

@@ -277,7 +277,6 @@ rule build_temperature_profiles_halloran:
     input:
         regions=RESOURCES + "regions_onshore_elec_s{simpl}_eso.geojson",
         population="data/population_layout/GB_residential_population_2011_1_km.tif",
-        # cutout="cutouts/"+ CDIR + config["heating"]["cutout"] + ".nc",
         cutout="cutouts/" + CDIR + config["atlite"]["default_cutout"] + ".nc",
     output:
         temp_air=RESOURCES + "temp_air_elec_s{simpl}_eso_halloran.nc",
@@ -293,7 +292,6 @@ rule build_heating_profiles_halloran:
     input:
         regions=RESOURCES + "regions_onshore_elec_s{simpl}_eso.geojson",
         population="data/population_layout/GB_residential_population_2011_1_km.tif",
-        # cutout = "cutouts/"+ CDIR + config["heating"]["cutout"] + ".nc",
         cutout="cutouts/" + CDIR + config["atlite"]["default_cutout"] + ".nc",
     output:
         profile_air_source_heating=RESOURCES + "load_air_source_heating_elec_s{simpl}_eso.nc",
@@ -750,7 +748,7 @@ rule prepare_network:
         capacity_constraints=RESOURCES + "fes_capacity_constraints_{fes}_{year}.csv",
         profile_boiler_heating="data/heat_load_profile_BDEW.csv",
         profile_air_source_heating=RESOURCES + "load_air_source_heating_elec_s{simpl}_eso.nc",
-        profile_ground_source_heating=RESOURCES + "load_ground_source_heating_elec_s{simpl}_eso.nc"
+        profile_ground_source_heating=RESOURCES + "load_ground_source_heating_elec_s{simpl}_eso.nc",
         clustered_pop_layout=RESOURCES + "pop_layout_elec_s{simpl}_eso.csv",
         heat_demand=RESOURCES + "heat_demand_total_elec_s{simpl}_eso.nc",
         cop_air_total=RESOURCES + "cop_air_total_elec_s{simpl}_eso.nc",

scripts/build_heating_profiles_halloran.py

@@ -254,12 +254,16 @@ def heat_demand_watson(cutout, source, **params):
     cutout_rio = cutout.data
     cutout_rio = cutout_rio.rio.write_crs('EPSG:4326')
     # transform to same CRS and resolution as cutout
-    population_match = population.rio.reproject_match(cutout_rio,
-                                                      resampling = rasterio.enums.Resampling.sum)
+    population_match = population.rio.reproject_match(
+        cutout_rio,
+        resampling = rasterio.enums.Resampling.sum
+        )
     population_match = population_match.squeeze().drop('band')
     population_match = population_match.where(population_match>0.)
+
     # population_match = population_match.fillna(0.)
     households = population_match/2.4 # England and Wales average household size
+
     # change large negative values to NaN-- may need to change to 0
     households = households.where(households>0.)
     households = households.fillna(0.)

scripts/prepare_network.py

@@ -585,26 +585,17 @@ def add_gas_ccs(n, costs):
 
 
 # adapted from `add_heat` method in `scripts/prepare_sector_network.py`
-def add_heat_pump_load(
-    n,
-    heat_demand_file,
-    ashp_cop_file,
-    energy_totals_file,
-    profile_boiler_intraday,
-    profile_air_source,
-    profile_ground_source,
-    scenario,
-    year,
-    flexopts,
-    flex_config,
-    ):
+def add_heat_pump_load(n, flexopts, costs):
 
-    year = int(year)
+    year = int(snakemake.wildcards.year)
+    scenario = snakemake.wildcards.fes
 
-    intraday_profiles = pd.read_csv(profile_boiler_intraday, index_col=0)
+    flex_config = snakemake.params.flex_config
+
+    intraday_profiles = pd.read_csv(snakemake.input["profile_boiler_heating"], index_col=0)
 
     daily_space_heat_demand = (
-        xr.open_dataarray(heat_demand_file)
+        xr.open_dataarray(snakemake.input["heat_demand"])
         .to_pandas()
         .reindex(index=n.snapshots, method="ffill")
     )
@@ -616,7 +607,7 @@ def add_heat_pump_load(
 
     daily_space_heat_demand = daily_space_heat_demand[nodes]
 
-    pop_weighted_energy_totals = pd.read_csv(energy_totals_file, index_col=0)
+    pop_weighted_energy_totals = pd.read_csv(snakemake.input["energy_totals"], index_col=0)
 
     sectors = ["residential"]
     uses = ["water", "space"]
@@ -648,33 +639,56 @@ def add_heat_pump_load(
         for region in daily_space_heat_demand.columns
     })
 
-    cop_air = xr.open_dataarray(ashp_cop_file).to_dataframe().iloc[:,0].unstack()
-    cop_air = cop_air[nodes]
+    # Assuming that profiles for boilers are the same as for resistive heaters
+    profile_resistive = heat_demand
 
-    cop = cop.rename(
-        columns={old: [col for col in heat_demand.columns if col in old][0]
-        for old in cop.columns}
-    )
+    get_profile = lambda file: xr.open_dataarray(file).to_pandas().transpose()[nodes]
+
+    profile_ashp = get_profile(snakemake.input["profile_air_source_heating"])
+    profile_gshp = get_profile(snakemake.input["profile_ground_source_heating"])
 
-    # get electricity demand by dividing heat demand through cop
-    heat_demand = heat_demand.divide(cop)
+    get_cop = lambda f: xr.open_dataarray(f).to_dataframe().iloc[:,0].unstack()[nodes]
 
-    # scale according to scenario
-    # get number of elec load through residential heating
+    cop_ashp = get_cop(snakemake.input["cop_air_total"])
+    cop_gshp = get_cop(snakemake.input["cop_ground_total"])
+    cop_resistive = costs.at["decentral resistive heater", "efficiency"]
 
-    hp_load_base, hp_load_future = (
+    heat_tech_shares = get_heating_shares(scenario, year)
+    heat_tech_shares /= heat_tech_shares.sum()
+
+    elec_heat_today, elec_heat_future = (
         get_electric_heat_demand(scenario, year, n.snapshots[0].year)
     )
 
-    heat_tech_shares = get_heating_shares(scenario, year)
+    def build_total_profile(profiles, shares, cops, total):
+
+        assert np.allclose(sum(shares), 1.), "Heating shares do not sum to 1."
+
+        process = lambda profile, share, cop: (profile / profile.sum().sum() * share).divide(cop)
+        profiles = list(map(process, profiles, shares, cops))
 
-    future_heat_demand = heat_demand / heat_demand.sum().sum() * hp_load_future
-    future_heat_demand.columns = heat_demand_spatial.nodes
+        scale = lambda profile: profile / sum([p.sum().sum() for p in profiles]) * total
+        profiles = list(map(scale, profiles))
 
-    # estimate demand subsumed in the general electricity demand, and remove it
-    base_heat_demand = heat_demand / heat_demand.sum().sum() * hp_load_base
-    base_heat_demand.columns = nodes
-    n.loads_t.p_set[nodes] -= base_heat_demand.reindex(nodes, axis=1)
+        return (p := pd.concat(profiles, axis=1)).T.groupby(p.columns).sum().T
+
+    assert heat_tech_shares.index.tolist() == ["Direct electric", "GSHP", "ASHP"], "Wrong order in heat shares."
+
+    profile_future = build_total_profile(
+        [profile_resistive, profile_gshp, profile_ashp],
+        heat_tech_shares.tolist(),
+        [cop_resistive, cop_gshp, cop_ashp],
+        elec_heat_future,
+    )
+
+    profile_today = build_total_profile(
+        [profile_resistive, profile_gshp, profile_ashp],
+        heat_tech_shares.tolist(),
+        [cop_resistive, cop_gshp, cop_ashp],
+        elec_heat_today,
+    )
+
+    n.loads_t.p_set[nodes] -= profile_today
 
     n.madd(
         "Bus",
@@ -689,15 +703,15 @@ def add_heat_pump_load(
         heat_demand_spatial.nodes,
         bus=heat_demand_spatial.nodes,
         carrier="heat demand",
-        p_set=future_heat_demand,
+        p_set=profile_future.set_axis(heat_demand_spatial.nodes, axis=1),
     )
 
     n.madd(
         "Link",
         heat_pumps_spatial.nodes,
         bus0=nodes,
         bus1=heat_demand_spatial.nodes,
-        carrier="heat pump",
+        carrier="electric heating",
         p_nom_extendable=True,
     )
 
@@ -749,7 +763,7 @@ def add_heat_pump_load(
         store_use_window = charging_window
 
         daily_p_max = (
-            future_heat_demand 
+            profile_future
             .groupby(pd.Grouper(freq="h"))
             .max()
             .reindex(s, method="ffill")
@@ -762,7 +776,7 @@ def add_heat_pump_load(
         p_min_pu = - (daily_p_max / p_nom).mul(discharging_window, axis=0)
 
         daily_e_max = (
-            future_heat_demand
+            profile_future
             .rolling(shift_size)
             .sum()
             .shift(-shift_size)
@@ -812,7 +826,7 @@ def add_heat_pump_load(
         max_hours = flex_config["water_tank_max_hours"]
 
         p_nom = (
-            future_heat_demand
+            profile_future
             .rolling(max_hours)
             .mean()
             .max()
@@ -831,7 +845,7 @@ def add_heat_pump_load(
         )
 
 
-def add_bev(n, transport_config, flex_config, flexopts):
+def add_bev(n, flexopts):
     """
     Adds BEV load and respective stores units;
     adapted from `add_land_transport` method in `scripts/prepare_sector_network.py`
@@ -840,6 +854,9 @@ def add_bev(n, transport_config, flex_config, flexopts):
     year = int(snakemake.wildcards.year)
     scenario = snakemake.wildcards.fes
 
+    transport_config = snakemake.params.sector_config
+    flex_config = snakemake.params.flex_config
+
     nodes = spatial.nodes
 
     if "go" in flexopts:
@@ -1972,26 +1989,10 @@ def prepare_validation(n, config):
     add_biogas(n, other_costs)
 
     logger.info("Adding heat pump load.")
-    add_heat_pump_load(
-        n,
-        snakemake.input["heat_demand"],
-        snakemake.input["cop_air_total"],
-        snakemake.input["energy_totals"],
-        snakemake.input["profile_boiler_heating"],
-        snakemake.input["profile_air_source_heating"],
-        snakemake.input["profile_ground_source_heating"],
-        snakemake.wildcards.fes,
-        snakemake.wildcards.year,
-        flexopts,
-        snakemake.params.flex_config,
-    )
+    add_heat_pump_load(n, flexopts, other_costs)
 
     logger.info("Adding BEV load.")
-    add_bev(n,
-        snakemake.params.sector_config,
-        snakemake.params.flex_config,
-        flexopts,
-    )
+    add_bev(n, flexopts)
 
     logger.info("Adding battery storage.")
     add_batteries(n, elec_costs, opts=opts)