diff --git a/config/config.default.yaml b/config/config.default.yaml index d9f692065..b95e3e3b0 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -382,6 +382,8 @@ sector: 2050: 1.0 district_heating_loss: 0.15 cluster_heat_buses: true + endogenous_transport: true + electrobiofuels: true bev_dsm_restriction_value: 0.75 bev_dsm_restriction_time: 7 transport_heating_deadband_upper: 20. @@ -773,21 +775,22 @@ solving: threads: 4 method: 2 # barrier crossover: 0 - BarConvTol: 1.e-6 + BarConvTol: 1.e-5 + FeasibilityTol: 1.e-4 + OptimalityTol: 1.e-4 Seed: 123 AggFill: 0 PreDual: 0 GURO_PAR_BARDENSETHRESH: 200 gurobi-numeric-focus: - name: gurobi - NumericFocus: 3 # Favour numeric stability over speed + # NumericFocus: 3 # Favour numeric stability over speed method: 2 # barrier crossover: 0 # do not use crossover BarHomogeneous: 1 # Use homogeneous barrier if standard does not converge BarConvTol: 1.e-5 FeasibilityTol: 1.e-4 OptimalityTol: 1.e-4 - ObjScale: -0.5 + # ObjScale: -0.5 threads: 8 Seed: 123 gurobi-fallback: # Use gurobi defaults @@ -913,6 +916,7 @@ plotting: gas pipeline: '#ebbca0' gas pipeline new: '#a87c62' # oil + electrobiofuels: '#a1b5ae' oil: '#c9c9c9' imported oil: '#a3a3a3' oil boiler: '#adadad' @@ -980,6 +984,7 @@ plotting: BEV charger: '#baf238' V2G: '#e5ffa8' land transport EV: '#baf238' + land transport demand: '#38baf2' Li ion: '#baf238' # hot water storage water tanks: '#e69487' diff --git a/config/test/config.perfect.yaml b/config/test/config.perfect.yaml index 5d77c9c51..118916aff 100644 --- a/config/test/config.perfect.yaml +++ b/config/test/config.perfect.yaml @@ -18,7 +18,7 @@ scenario: clusters: - 5 sector_opts: - - 8760h-T-H-B-I-A-dist1 + - 24h-T-H-B-I-A-dist1 planning_horizons: - 2030 - 2040 diff --git a/doc/release_notes.rst b/doc/release_notes.rst index ad1bbf4fa..07c64df4b 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -10,6 +10,12 @@ Release Notes Upcoming Release ================ + +* fix bug in land transport for temperature correction of ICEs and fuel cell cars + +* restructure land transport, demand is now attached to one single node, the +different car types (ICE, EV, fuel cell) are modelled as links + * Upgrade to Snakemake v8.5+. This version is the new minimum version required. To upgrade an existing environment, run ``conda install -c bioconda snakemake-minimal">=8.5"`` and ``pip install snakemake-storage-plugin-http`` diff --git a/rules/build_electricity.smk b/rules/build_electricity.smk index 24f328eb8..f3451eed0 100644 --- a/rules/build_electricity.smk +++ b/rules/build_electricity.smk @@ -310,7 +310,7 @@ rule build_renewable_profiles: benchmarks("build_renewable_profiles_{technology}") threads: config["atlite"].get("nprocesses", 4) resources: - mem_mb=config["atlite"].get("nprocesses", 4) * 5000, + mem_mb=config["atlite"].get("nprocesses", 4) * 10000, wildcard_constraints: technology="(?!hydro).*", # Any technology other than hydro conda: @@ -481,7 +481,7 @@ rule simplify_network: benchmarks("simplify_network/elec_s{simpl}") threads: 1 resources: - mem_mb=12000, + mem_mb=30000, conda: "../envs/environment.yaml" script: diff --git a/rules/build_sector.smk b/rules/build_sector.smk index 8d41e893b..ad8b3b15c 100644 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -323,6 +323,31 @@ rule build_biomass_potentials: "../scripts/build_biomass_potentials.py" +rule build_existing_car_ages: + input: + ACEA_report=storage( + "https://www.acea.auto/files/ACEA-Report-Vehicles-on-European-roads-.pdf", + keep_local=True, + ), + output: + car_ages=resources( + "car_ages.csv" + ), + truck_ages=resources( + "truck_ages.csv" + ), + threads: 1 + resources: + mem_mb=1000, + log: + logs("build_existing_car_ages.log"), + benchmark: + benchmarks("build_existing_car_ages") + conda: + "../envs/environment.yaml" + script: + "../scripts/build_existing_car_ages.py" + rule build_biomass_transport_costs: input: transport_cost_data=storage( @@ -735,6 +760,7 @@ rule build_transport_demand: transport_data=resources("transport_data.csv"), traffic_data_KFZ="data/bundle-sector/emobility/KFZ__count", traffic_data_Pkw="data/bundle-sector/emobility/Pkw__count", + traffic_data_Lkw="data/bundle-sector/emobility/Lkw__count", temp_air_total=resources("temp_air_total_elec_s{simpl}_{clusters}.nc"), output: transport_demand=resources("transport_demand_s{simpl}_{clusters}.csv"), diff --git a/rules/solve_myopic.smk b/rules/solve_myopic.smk index 57b8a9d3f..599bf1f1f 100644 --- a/rules/solve_myopic.smk +++ b/rules/solve_myopic.smk @@ -10,6 +10,8 @@ rule add_existing_baseyear: existing_capacities=config_provider("existing_capacities"), costs=config_provider("costs"), input: + car_ages=resources("car_ages.csv"), + truck_ages=resources("truck_ages.csv"), network=RESULTS + "prenetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc", powerplants=resources("powerplants.csv"), diff --git a/rules/solve_perfect.smk b/rules/solve_perfect.smk index a565d9788..f570cfdce 100644 --- a/rules/solve_perfect.smk +++ b/rules/solve_perfect.smk @@ -8,6 +8,8 @@ rule add_existing_baseyear: existing_capacities=config_provider("existing_capacities"), costs=config_provider("costs"), input: + car_ages=resources("car_ages.csv"), + truck_ages=resources("truck_ages.csv"), network=RESULTS + "prenetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc", powerplants=resources("powerplants.csv"), @@ -113,7 +115,8 @@ rule solve_sector_network_perfect: + "postnetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_brownfield_all_years.nc", threads: solver_threads resources: - mem_mb=config_provider("solving", "mem"), + mem_mb=config_provider("solving", "mem_mb"), + runtime=config_provider("solving", "runtime", default="24h"), shadow: "shallow" log: diff --git a/scripts/_helpers.py b/scripts/_helpers.py index a1504c3c7..b244beb12 100644 --- a/scripts/_helpers.py +++ b/scripts/_helpers.py @@ -570,9 +570,18 @@ def update_config_from_wildcards(config, w, inplace=True): if "CCL" in opts: config["solving"]["constraints"]["CCL"] = True - + + co2_budget = { + "1p5": 34.2 , # 25.7 # Budget in Gt CO2 for 1.5 for Europe, global 420 Gt, assuming per capita share + "1p6": 43.259666 , # 35 # Budget in Gt CO2 for 1.6 for Europe, global 580 Gt + "1p7": 51.4, # 45 # Budget in Gt CO2 for 1.7 for Europe, global 800 Gt + "2p0": 69.778 , # 73.9 # Budget in Gt CO2 for 2 for Europe, global 1170 Gt + } eq_value = get_opt(opts, r"^EQ+\d*\.?\d+(c|)") for o in opts: + m = re.match(r"^\d+p\d$", o, re.IGNORECASE) + if m is not None: + config["co2_budget"]= co2_budget[m.group(0)] if eq_value is not None: config["solving"]["constraints"]["EQ"] = eq_value elif "EQ" in o: diff --git a/scripts/add_existing_baseyear.py b/scripts/add_existing_baseyear.py index 0bbe19f0e..852f52079 100644 --- a/scripts/add_existing_baseyear.py +++ b/scripts/add_existing_baseyear.py @@ -21,7 +21,12 @@ update_config_from_wildcards, ) from add_electricity import sanitize_carriers -from prepare_sector_network import cluster_heat_buses, define_spatial, prepare_costs +from prepare_sector_network import ( + cluster_heat_buses, + define_spatial, + get, + prepare_costs, +) logger = logging.getLogger(__name__) cc = coco.CountryConverter() @@ -55,6 +60,67 @@ def add_build_year_to_new_assets(n, baseyear): c.pnl[attr] = c.pnl[attr].rename(columns=rename) +def add_existing_land_transport(baseyear, options): + # today ICE capacity assuming all internal combustion + share = get(options["land_transport_ice_share"], baseyear) + ice_i = n.links[n.links.carrier == "land transport oil"].index + p_nom = n.links.loc[ice_i, "p_nom"] / share + efficiency = n.links_t.efficiency[ice_i] + p_max_pu = n.links_t.p_max_pu[ice_i] + + car_ages = pd.read_csv(snakemake.input.car_ages, index_col=[0]).iloc[:,:-2] + car_ages.columns = car_ages.columns.astype(int) + # group data in 5 years interval + interval = 5 + # mapping forward (mapping backward would be year//5*5) + group_mapping = {year: year//interval*interval+4 for year in car_ages.columns} + + grouped = car_ages.T.groupby(group_mapping).sum().T + + pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0) + + grouped = (grouped.reindex(pop_layout.ct) + .fillna(grouped.mean()).set_index(pop_layout.index)) + + for build_year in grouped.columns: + df = n.links.loc[ice_i] + df = df[df.lifetime + build_year > baseyear] + if df.empty: + continue + share = grouped[build_year] + df["build_year"] = build_year + df["p_nom"] = p_nom.mul(share.values) + df["p_nom_extendable"] = False + df.rename( + index=lambda x: x.replace(f"-{baseyear}", f"-{build_year}"), inplace=True + ) + profile = p_max_pu.rename( + columns=lambda x: x.replace(f"-{baseyear}", f"-{build_year}")) + eff = efficiency.rename( + columns=lambda x: x.replace(f"-{baseyear}", f"-{build_year}")) + + n.madd( + "Link", + df.index, + bus0=df.bus0, + bus1=df.bus1, + bus2=df.bus2, + carrier=df.carrier, + efficiency=eff, + capital_cost=df.capital_cost, + marginal_cost=df.marginal_cost, + efficiency2=df.efficiency2, + p_nom_extendable=False, + p_nom=df.p_nom, + p_min_pu=profile, + p_max_pu=profile, + build_year=df.build_year, + lifetime=df.lifetime, + ) + + n.links.loc[ice_i, "p_nom"] = 0 + + def add_existing_renewables(df_agg): """ Append existing renewables to the df_agg pd.DataFrame with the conventional @@ -553,13 +619,13 @@ def add_heating_capacities_installed_before_baseyear( snakemake = mock_snakemake( "add_existing_baseyear", - # configfiles="config/test/config.myopic.yaml", + # configfiles="config/config.myopic.yaml", simpl="", clusters="37", ll="v1.0", opts="", - sector_opts="8760-T-H-B-I-A-dist1", - planning_horizons=2020, + sector_opts="49sn-T-H-B-I-A-dist1", + planning_horizons=2025, ) configure_logging(snakemake) @@ -619,6 +685,9 @@ def add_heating_capacities_installed_before_baseyear( if options.get("cluster_heat_buses", False): cluster_heat_buses(n) + if options["endogenous_transport"]: + add_existing_land_transport(baseyear, options) + n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) sanitize_carriers(n, snakemake.config) diff --git a/scripts/build_energy_totals.py b/scripts/build_energy_totals.py index 1ffc4ae2c..e365c0a4e 100644 --- a/scripts/build_energy_totals.py +++ b/scripts/build_energy_totals.py @@ -240,9 +240,10 @@ def idees_per_country(ct, year, base_dir): ct_totals["total agriculture"] = df[row] # transport - + df = pd.read_excel(fn_transport, "TrRoad_ene", index_col=0)[year] - + + # energy consumptiob by fuel (ktoe) ct_totals["total road"] = df["by fuel (EUROSTAT DATA)"] ct_totals["electricity road"] = df["Electricity"] @@ -269,9 +270,13 @@ def idees_per_country(ct, year, base_dir): row = "Heavy duty vehicles (Diesel oil incl. biofuels)" ct_totals["total heavy duty road freight"] = df[row] - + + # vehicle efficiency (kgoe/100km) assert df.index[61] == "Passenger cars" ct_totals["passenger car efficiency"] = df.iloc[61] + + assert df.index[81] == "Heavy duty vehicles" + ct_totals["heavy duty efficiency"] = df.iloc[81] df = pd.read_excel(fn_transport, "TrRail_ene", index_col=0)[year] @@ -330,9 +335,25 @@ def idees_per_country(ct, year, base_dir): ct_totals["total domestic navigation"] = df["by fuel (EUROSTAT DATA)"] df = pd.read_excel(fn_transport, "TrRoad_act", index_col=0)[year] - + + # total number of light duty vehicles assert df.index[85] == "Passenger cars" - ct_totals["passenger cars"] = df.iloc[85] + + ct_totals["Number Passenger cars"] = df.iloc[85] + + assert df.index[84] == "Powered 2-wheelers" + ct_totals['Number Powered 2-wheelers'] = df.iloc[84] + + assert df.index[99] == "Light duty vehicles" + ct_totals['Number Light duty vehicles'] = df.iloc[99] + + # total number of heavy duty vehicles + + assert df.index[92] == "Motor coaches, buses and trolley buses" + ct_totals["Number Motor coaches, buses and trolley buses"] = df.iloc[92] + + assert df.index[105] == 'Heavy duty vehicles' + ct_totals['Number Heavy duty vehicles'] = df.iloc[105] return pd.Series(ct_totals, name=ct) @@ -356,11 +377,12 @@ def build_idees(countries, year): totals = pd.concat(totals_list, axis=1) # convert ktoe to TWh - exclude = totals.index.str.fullmatch("passenger cars") + exclude = totals.index.str.contains("Number") totals.loc[~exclude] *= 11.63 / 1e3 # convert TWh/100km to kWh/km totals.loc["passenger car efficiency"] *= 10 + totals.loc["heavy duty efficiency"] *= 10 return totals.T @@ -368,7 +390,13 @@ def build_idees(countries, year): def build_energy_totals(countries, eurostat, swiss, idees): eurostat_fuels = {"electricity": "Electricity", "total": "Total all products"} - to_drop = ["passenger cars", "passenger car efficiency"] + to_drop = ["Number Passenger cars", + "Number Powered 2-wheelers", + "Number Light duty vehicles", + "Number Motor coaches, buses and trolley buses", + "Number Heavy duty vehicles", + "passenger car efficiency", + "heavy duty efficiency"] df = idees.reindex(countries).drop(to_drop, axis=1) eurostat_countries = eurostat.index.levels[0] @@ -671,34 +699,39 @@ def build_transport_data(countries, population, idees): transport_data = pd.DataFrame(index=countries) # collect number of cars - - transport_data["number cars"] = idees["passenger cars"] + car_cols = ["Number Passenger cars", + "Number Powered 2-wheelers", + "Number Light duty vehicles", + "Number Motor coaches, buses and trolley buses", + "Number Heavy duty vehicles"] + + transport_data[car_cols] = idees[car_cols] # CH from http://ec.europa.eu/eurostat/statistics-explained/index.php/Passenger_cars_in_the_EU#Luxembourg_has_the_highest_number_of_passenger_cars_per_inhabitant if "CH" in countries: - transport_data.at["CH", "number cars"] = 4.136e6 - - missing = transport_data.index[transport_data["number cars"].isna()] - if not missing.empty: - logger.info( - f"Missing data on cars from:\n{list(missing)}\nFilling gaps with averaged data." - ) - - cars_pp = transport_data["number cars"] / population - transport_data.loc[missing, "number cars"] = cars_pp.mean() * population + transport_data.at["CH", "Number Passenger cars"] = 4.136e6 + + for col in car_cols: + missing = transport_data.index[transport_data[col].isna()] + if not missing.empty: + logger.info( + f"Missing data on {col} from:\n{list(missing)}\nFilling gaps with averaged data." + ) + cars_pp = transport_data[col] / population + transport_data.loc[missing, col] = cars_pp.mean() * population # collect average fuel efficiency in kWh/km - - transport_data["average fuel efficiency"] = idees["passenger car efficiency"] - - missing = transport_data.index[transport_data["average fuel efficiency"].isna()] - if not missing.empty: - logger.info( - f"Missing data on fuel efficiency from:\n{list(missing)}\nFilling gapswith averaged data." - ) - - fill_values = transport_data["average fuel efficiency"].mean() - transport_data.loc[missing, "average fuel efficiency"] = fill_values + for vehicle_type in ["passenger car", "heavy duty"]: + transport_data[f"average fuel efficiency {vehicle_type}"] = idees[f"{vehicle_type} efficiency"] + + missing = transport_data.index[transport_data[f"average fuel efficiency {vehicle_type}"].isna()] + if not missing.empty: + logger.info( + f"Missing data on fuel efficiency from:\n{list(missing)}\nFilling gapswith averaged data." + ) + + fill_values = transport_data[f"average fuel efficiency {vehicle_type}"].mean() + transport_data.loc[missing, f"average fuel efficiency {vehicle_type}"] = fill_values return transport_data @@ -827,7 +860,7 @@ def rescale_idees_from_eurostat( return energy - +#%% if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake diff --git a/scripts/build_existing_car_ages.py b/scripts/build_existing_car_ages.py new file mode 100644 index 000000000..fc360e5ac --- /dev/null +++ b/scripts/build_existing_car_ages.py @@ -0,0 +1,101 @@ +# -*- coding: utf-8 -*- +# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors +# +# SPDX-License-Identifier: MIT +""" +Read in existing car and truck ages from ACEA report (2024) +""" + + +import tabula +import numpy as np +import country_converter as coco +import pandas as pd +import logging +logger = logging.getLogger(__name__) +from _helpers import configure_logging, set_scenario_config + +def distirbute_cars(row): + if row.name=="EU": return row + i = 0 + for j, number in enumerate(row): + if np.isnan(number): + i += 1 + else: + distributed = number / (i+1) + row.iloc[j-i:j+1] = distributed + i = 0 + return row + +def distribute_older_ages(ds, number_of_years=8): + "Distribute columne <10 years based on averaged data" + share_last_year = ds.fillna(0).loc['>10 years'] + ds.drop('>10 years', inplace=True) + mean = ds.fillna(0).mean() + # convert index to int + ds.index = ds.index.astype(int) + years_i = np.arange(2012, 2012-number_of_years, -1) + av = share_last_year/mean + distributed = pd.Series(data=0, index=years_i, dtype="float64") + for i, year in enumerate(years_i): + if i+1 < av: + distributed[year] = mean + elif (i == int(av)): + distributed[year] = (av - i) * mean + else: + distributed[year] = 0 + + if share_last_year - distributed.iloc[:-1].sum()> mean: + distributed.iloc[-1] = share_last_year - distributed.iloc[:-1].sum() + + return pd.concat([ds, distributed]) + +def read_ages_from_pdf(fn, page, name="CARS BY AGE"): + tables = tabula.read_pdf(fn, pages=page) + table = tables[0] + table.set_index(name, inplace=True) + + table.columns = table.iloc[1,:] + table = table.iloc[2:,:] + + table = table.applymap(lambda x: np.nan if isinstance(x, str) and '–' in x else x) + table = table.apply(lambda x: (x.str.replace(",","")).astype(float)) + + # convert into to iso-code + cc = coco.CountryConverter() + table.index = cc.convert(table.index, to="iso2") + table.drop('not found', inplace=True) + + table = table.apply(lambda x: distirbute_cars(x), axis=1) + + # get shares instead of total numbers + table.iloc[:, :-2] = table.iloc[:, :-2].div(table["Total"], axis=0) + table.rename(columns={'(in years)': "Average age"}, inplace=True) + df = table.iloc[:,:-2] + # swedish data is missing, filling with average data + df = df.fillna(df.mean()) + final_table = df.apply(distribute_older_ages, axis=1) + + return pd.concat([final_table, table.iloc[:, -2:]], axis=1) + + +#%% +if __name__ == "__main__": + if "snakemake" not in globals(): + from _helpers import mock_snakemake + snakemake = mock_snakemake( + "build_existing_car_ages", + ) + + configure_logging(snakemake) + set_scenario_config(snakemake) + fn = snakemake.input.ACEA_report + # car ages + car_ages = read_ages_from_pdf(fn, page=11, name="CARS BY AGE") + + # trucks ages + truck_ages = read_ages_from_pdf(fn, page=13, name="TRUCKS BY AGE") + + + car_ages.to_csv(snakemake.output.car_ages) + truck_ages.to_csv(snakemake.output.truck_ages) diff --git a/scripts/build_transport_demand.py b/scripts/build_transport_demand.py index de561e3f8..ace11cbf6 100644 --- a/scripts/build_transport_demand.py +++ b/scripts/build_transport_demand.py @@ -19,41 +19,59 @@ def build_nodal_transport_data(fn, pop_layout): + # get numbers of car and fuel efficiency per country transport_data = pd.read_csv(fn, index_col=0) + # TODO check what is going wrong in LU + # transport_data.loc['LU', 'average fuel efficiency heavy duty'] = transport_data.loc['SI', 'average fuel efficiency heavy duty'] + + # break number of cars down to nodal level based on population density nodal_transport_data = transport_data.loc[pop_layout.ct].fillna(0.0) + nodal_transport_data.index = pop_layout.index - nodal_transport_data["number cars"] = ( - pop_layout["fraction"] * nodal_transport_data["number cars"] + car_cols = ["Number Passenger cars", + "Number Powered 2-wheelers", + "Number Light duty vehicles", + "Number Motor coaches, buses and trolley buses", + "Number Heavy duty vehicles"] + + nodal_transport_data[car_cols] = ( + nodal_transport_data[car_cols].mul(pop_layout["fraction"], axis=0) ) - nodal_transport_data.loc[ - nodal_transport_data["average fuel efficiency"] == 0.0, - "average fuel efficiency", - ] = transport_data["average fuel efficiency"].mean() - + # fill missing fuel efficiency [kWh/km] with average data + eff_cols = ['average fuel efficiency passenger car', + 'average fuel efficiency heavy duty'] + for col in eff_cols: + nodal_transport_data.loc[ + nodal_transport_data[col] == 0.0, + col, + ] = transport_data[col].mean() + return nodal_transport_data -def build_transport_demand(traffic_fn, airtemp_fn, nodes, nodal_transport_data): - ## Get overall demand curve for all vehicles +def get_shape(traffic_fn): traffic = pd.read_csv(traffic_fn, skiprows=2, usecols=["count"]).squeeze("columns") + # create annual profile take account time zone + summer time transport_shape = generate_periodic_profiles( dt_index=snapshots, nodes=nodes, weekly_profile=traffic.values, ) transport_shape = transport_shape / transport_shape.sum() - - # electric motors are more efficient, so alter transport demand - - plug_to_wheels_eta = options["bev_plug_to_wheel_efficiency"] - battery_to_wheels_eta = plug_to_wheels_eta * options["bev_charge_efficiency"] - - efficiency_gain = ( - nodal_transport_data["average fuel efficiency"] / battery_to_wheels_eta - ) + + return transport_shape + +def build_transport_demand(traffic_fn_Pkw, traffic_fn_Lkw, + airtemp_fn, nodes, nodal_transport_data): + """ + Returns transport demand per bus in unit kinetic energy. + """ + transport_shape_light = get_shape(traffic_fn_Pkw) + + transport_shape_heavy = get_shape(traffic_fn_Lkw) # get heating demand for correction to demand time series temperature = xr.open_dataarray(airtemp_fn).to_pandas() @@ -67,30 +85,41 @@ def build_transport_demand(traffic_fn, airtemp_fn, nodes, nodal_transport_data): options["ICE_upper_degree_factor"], ) - dd_EV = transport_degree_factor( - temperature, - options["transport_heating_deadband_lower"], - options["transport_heating_deadband_upper"], - options["EV_lower_degree_factor"], - options["EV_upper_degree_factor"], - ) - # divide out the heating/cooling demand from ICE totals - # and multiply back in the heating/cooling demand for EVs - ice_correction = (transport_shape * (1 + dd_ICE)).sum() / transport_shape.sum() - - energy_totals_transport = ( - pop_weighted_energy_totals["total road"] - + pop_weighted_energy_totals["total rail"] - - pop_weighted_energy_totals["electricity rail"] - ) - - return ( - (transport_shape.multiply(energy_totals_transport) * 1e6 * nyears) - .divide(efficiency_gain * ice_correction) - .multiply(1 + dd_EV) - ) - + ice_correction_light = (transport_shape_light * (1 + dd_ICE)).sum() / transport_shape_light.sum() + ice_correction_heavy = (transport_shape_light * (1 + dd_ICE)).sum() / transport_shape_heavy.sum() + + light_duty_cols = ['total two-wheel', 'total passenger cars', + 'total light duty road freight'] + + heavy_duty_cols = ['total other road passenger', # motor coaches, buses, trolley buses + 'total heavy duty road freight'] + + light_duty = pop_weighted_energy_totals[light_duty_cols].sum(axis=1) + + heavy_duty = pop_weighted_energy_totals[heavy_duty_cols].sum(axis=1) + + rail = pop_weighted_energy_totals["total rail"] - pop_weighted_energy_totals["electricity rail"] + + def get_demand(transport_shape, energy_totals_transport, nyears, + fuel_efficiency, ice_correction, name): + + demand = (transport_shape.multiply(energy_totals_transport) * 1e6 * nyears).divide( + fuel_efficiency * ice_correction + ) + + return pd.concat([demand], keys=[name], axis=1) + + demand_light = get_demand(transport_shape_light, light_duty, nyears, + # convert 1 kWh/km = 0.1 MWh/ 100 km + 0.1*nodal_transport_data["average fuel efficiency passenger car"], + ice_correction_light, name="light") + demand_heavy = get_demand(transport_shape_heavy, (heavy_duty + rail), nyears, + 0.1*nodal_transport_data["average fuel efficiency heavy duty"], + ice_correction_heavy, name="heavy") + + return pd.concat([demand_light, demand_heavy], axis=1) + def transport_degree_factor( temperature, @@ -125,11 +154,14 @@ def bev_availability_profile(fn, snapshots, nodes, options): """ Derive plugged-in availability for passenger electric vehicles. """ + # car count in typical week traffic = pd.read_csv(fn, skiprows=2, usecols=["count"]).squeeze("columns") - + # maximum share plugged-in availability for passenger electric vehicles avail_max = options["bev_avail_max"] + # average share plugged-in availability for passenger electric vehicles avail_mean = options["bev_avail_mean"] + # linear scaling, highest when traffic is lowest, decreases if traffic increases avail = avail_max - (avail_max - avail_mean) * (traffic - traffic.min()) / ( traffic.mean() - traffic.min() ) @@ -150,6 +182,8 @@ def bev_availability_profile(fn, snapshots, nodes, options): def bev_dsm_profile(snapshots, nodes, options): dsm_week = np.zeros((24 * 7,)) + # assuming that at a certain time ("bev_dsm_restriction_time") EVs have to + # be charged to a minimum value (defined in bev_dsm_restriction_value) dsm_week[(np.arange(0, 7, 1) * 24 + options["bev_dsm_restriction_time"])] = options[ "bev_dsm_restriction_value" ] @@ -161,6 +195,7 @@ def bev_dsm_profile(snapshots, nodes, options): ) +# %% if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake @@ -168,7 +203,7 @@ def bev_dsm_profile(snapshots, nodes, options): snakemake = mock_snakemake( "build_transport_demand", simpl="", - clusters=48, + clusters=37, ) configure_logging(snakemake) set_scenario_config(snakemake) @@ -192,7 +227,8 @@ def bev_dsm_profile(snapshots, nodes, options): ) transport_demand = build_transport_demand( - snakemake.input.traffic_data_KFZ, + snakemake.input.traffic_data_Pkw, + snakemake.input.traffic_data_Lkw, snakemake.input.temp_air_total, nodes, nodal_transport_data, diff --git a/scripts/make_summary.py b/scripts/make_summary.py index 18642afc5..e083bc672 100644 --- a/scripts/make_summary.py +++ b/scripts/make_summary.py @@ -294,6 +294,7 @@ def calculate_energy(n, label, energy): ) # remove values where bus is missing (bug in nomopyomo) no_bus = c.df.index[c.df["bus" + port] == ""] + if totals.empty: continue totals.loc[no_bus] = float( n.component_attrs[c.name].loc["p" + port, "default"] ) @@ -651,9 +652,13 @@ def make_summaries(networks_dict): df = {output: pd.DataFrame(columns=columns, dtype=float) for output in outputs} for label, filename in networks_dict.items(): logger.info(f"Make summary for scenario {label}, using {filename}") - - n = pypsa.Network(filename) - + try: + n = pypsa.Network(filename) + except FileNotFoundError: + logger.info(f"{label} not yet solved.") + continue + + if not hasattr(n, "objective"): continue assign_carriers(n) assign_locations(n) @@ -667,7 +672,7 @@ def to_csv(df): for key in df: df[key].to_csv(snakemake.output[key]) - +#%% if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake diff --git a/scripts/prepare_perfect_foresight.py b/scripts/prepare_perfect_foresight.py index f7e8495ec..812f68d07 100644 --- a/scripts/prepare_perfect_foresight.py +++ b/scripts/prepare_perfect_foresight.py @@ -18,6 +18,7 @@ update_config_from_wildcards, ) from add_existing_baseyear import add_build_year_to_new_assets +from prepare_sector_network import adjust_transport_temporal_agg from pypsa.descriptors import expand_series from pypsa.io import import_components_from_dataframe from six import iterkeys @@ -165,7 +166,8 @@ def concat_networks(years): year = years[i] network = pypsa.Network(network_path) adjust_electricity_grid(network, year, years) - add_build_year_to_new_assets(network, year) + if not i == 0: + add_build_year_to_new_assets(network, year) # static ---------------------------------- for component in network.iterate_components( @@ -334,7 +336,7 @@ def set_carbon_constraints(n): n.add( "GlobalConstraint", "carbon_neutral", - type="co2_limit", + type="co2_atmosphere", carrier_attribute="co2_emissions", sense="<=", constant=0, @@ -497,8 +499,8 @@ def apply_time_segmentation_perfect( simpl="", opts="", clusters="37", - ll="v1.5", - sector_opts="1p7-4380H-T-H-B-I-A-dist1", + ll="v1.0", + sector_opts="1p5-730H-T-H-B-I-A-dist1", ) configure_logging(snakemake) set_scenario_config(snakemake) @@ -520,6 +522,7 @@ def apply_time_segmentation_perfect( segments = snakemake.params.time_resolution if isinstance(segments, (int, float)): n = apply_time_segmentation_perfect(n, segments, solver_name=solver_name) + adjust_transport_temporal_agg(n) # adjust global constraints lv limit if the same for all years n = adjust_lvlimit(n) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 2e8bf6fdb..bf5892d90 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -24,6 +24,7 @@ ) from add_electricity import calculate_annuity, sanitize_carriers, sanitize_locations from build_energy_totals import build_co2_totals, build_eea_co2, build_eurostat_co2 +from build_transport_demand import transport_degree_factor from networkx.algorithms import complement from networkx.algorithms.connectivity.edge_augmentation import k_edge_augmentation from prepare_network import maybe_adjust_costs_and_potentials @@ -1496,90 +1497,122 @@ def add_storage_and_grids(n, costs): ) -def add_land_transport(n, costs): - # TODO options? +def check_land_transport_shares(shares): + # Sums up the shares, ignoring None values + total_share = sum(filter(None, shares)) + if total_share != 1: + logger.warning( + f"Total land transport shares sum up to {total_share:.2%}," + "corresponding to increased or decreased demand assumptions." + ) - logger.info("Add land transport") - transport = pd.read_csv( - snakemake.input.transport_demand, index_col=0, parse_dates=True - ) - number_cars = pd.read_csv(snakemake.input.transport_data, index_col=0)[ - "number cars" - ] - avail_profile = pd.read_csv( - snakemake.input.avail_profile, index_col=0, parse_dates=True - ) - dsm_profile = pd.read_csv( - snakemake.input.dsm_profile, index_col=0, parse_dates=True +def get_temp_efficency( + car_efficiency, + temperature, + deadband_lw, + deadband_up, + degree_factor_lw, + degree_factor_up, +): + """ + Correct temperature depending on heating and cooling for respective car + type. + """ + # temperature correction for EVs + dd = transport_degree_factor( + temperature, + deadband_lw, + deadband_up, + degree_factor_lw, + degree_factor_up, ) - fuel_cell_share = get(options["land_transport_fuel_cell_share"], investment_year) - electric_share = get(options["land_transport_electric_share"], investment_year) - ice_share = get(options["land_transport_ice_share"], investment_year) + temp_eff = 1 / (1 + dd) - total_share = fuel_cell_share + electric_share + ice_share - if total_share != 1: - logger.warning( - f"Total land transport shares sum up to {total_share:.2%}, corresponding to increased or decreased demand assumptions." - ) + return car_efficiency * temp_eff - logger.info(f"FCEV share: {fuel_cell_share*100}%") - logger.info(f"EV share: {electric_share*100}%") - logger.info(f"ICEV share: {ice_share*100}%") - nodes = pop_layout.index +def add_EVs( + n, + nodes, + avail_profile, + dsm_profile, + p_set, + electric_share, + number_cars, + temperature, +): - if electric_share > 0: - n.add("Carrier", "Li ion") + n.add("Carrier", "Li ion") - n.madd( - "Bus", - nodes, - suffix=" EV battery", - location=nodes, - carrier="Li ion", - unit="MWh_el", - ) + n.madd( + "Bus", + nodes, + suffix=" EV battery", + location=nodes, + carrier="Li ion", + unit="MWh_el", + ) + + # https://ev-database.org/de/cheatsheet/energy-consumption-electric-car + # average energy consumption 188 Wh/km = 0.188 kWh/km = 18.8 kWh/100 km = 0.0188 MWh/ 100 km + # 1/0.188 -> 1 MWh = 53.19 * 100 km + costs.at["Battery electric (passenger cars)", "efficiency"] = 53.19 + car_efficiency = costs.at["Battery electric (passenger cars)", "efficiency"] - p_set = ( - electric_share - * ( - transport[nodes] - + cycling_shift(transport[nodes], 1) - + cycling_shift(transport[nodes], 2) - ) - / 3 - ) + # temperature corrected efficiency + efficiency = get_temp_efficency( + car_efficiency, + temperature, + options["transport_heating_deadband_lower"], + options["transport_heating_deadband_upper"], + options["EV_lower_degree_factor"], + options["EV_upper_degree_factor"], + ) + suffix = " land transport EV" - n.madd( - "Load", - nodes, - suffix=" land transport EV", - bus=nodes + " EV battery", - carrier="land transport EV", - p_set=p_set, - ) + p_shifted = (p_set + cycling_shift(p_set, 1) + cycling_shift(p_set, 2)) / 3 - p_nom = number_cars * options.get("bev_charge_rate", 0.011) * electric_share + cyclic_eff = p_set.div(p_shifted) - n.madd( - "Link", - nodes, - suffix=" BEV charger", - bus0=nodes, - bus1=nodes + " EV battery", - p_nom=p_nom, - carrier="BEV charger", - p_max_pu=avail_profile[nodes], - efficiency=options.get("bev_charge_efficiency", 0.9), - # These were set non-zero to find LU infeasibility when availability = 0.25 - # p_nom_extendable=True, - # p_nom_min=p_nom, - # capital_cost=1e6, #i.e. so high it only gets built where necessary - ) + efficiency *= cyclic_eff + + p_nom = electric_share * p_set.div(efficiency).max() + + profile = p_set.div(efficiency) / p_set.div(efficiency).max() + + n.madd( + "Link", + nodes, + suffix=suffix, + bus0=nodes + " EV battery", + bus1=nodes + " land transport", + carrier="land transport EV", + efficiency=efficiency, + p_min_pu=profile, + p_max_pu=profile, + p_nom=p_nom, + p_nom_extendable=False, + lifetime=1, + ) - if electric_share > 0 and options["v2g"]: + p_nom = number_cars * options.get("bev_charge_rate", 0.011) * electric_share + + n.madd( + "Link", + nodes, + suffix=" BEV charger", + bus0=nodes, + bus1=nodes + " EV battery", + p_nom=p_nom, + carrier="BEV charger", + p_max_pu=avail_profile[nodes], + lifetime=1, + efficiency=options.get("bev_charge_efficiency", 0.9), + ) + + if options["v2g"]: n.madd( "Link", nodes, @@ -1589,10 +1622,11 @@ def add_land_transport(n, costs): p_nom=p_nom, carrier="V2G", p_max_pu=avail_profile[nodes], + lifetime=1, efficiency=options.get("bev_charge_efficiency", 0.9), ) - if electric_share > 0 and options["bev_dsm"]: + if options["bev_dsm"]: e_nom = ( number_cars * options.get("bev_energy", 0.05) @@ -1610,63 +1644,270 @@ def add_land_transport(n, costs): e_nom=e_nom, e_max_pu=1, e_min_pu=dsm_profile[nodes], - ) + lifetime=1, + ) + +def add_electrobiofuels(n, nodes): + + print('Adding electrobiofuels') + efuel_scale_factor = costs.at['BtL', 'C stored'] + n.madd("Link", + nodes + " electrobiofuels", + bus0=spatial.biomass.nodes, + bus1=spatial.oil.nodes, + bus2=spatial.h2.nodes, + bus3="co2 atmosphere", + carrier="electrobiofuels", + lifetime=costs.at['electrobiofuels', 'lifetime'], + efficiency=costs.at['electrobiofuels', 'efficiency-biomass'], + efficiency2=-costs.at['electrobiofuels', 'efficiency-hydrogen'], + efficiency3=-costs.at['solid biomass', 'CO2 intensity'] + costs.at['BtL', 'CO2 stored'] * (1 - costs.at['Fischer-Tropsch', 'capture rate']), + p_nom_extendable=True, + capital_cost=costs.at['BtL', 'fixed'] * costs.at['electrobiofuels', 'efficiency-biomass'] \ + + efuel_scale_factor * costs.at['Fischer-Tropsch', 'fixed'] * costs.at['electrobiofuels', 'efficiency-hydrogen'], + marginal_cost=costs.at['BtL', 'VOM'] * costs.at['electrobiofuels', 'efficiency-biomass'] \ + + efuel_scale_factor * costs.at['Fischer-Tropsch', 'VOM'] * costs.at['electrobiofuels', 'efficiency-hydrogen'] + ) + + +def add_fuel_cell_cars(n, nodes, p_set, fuel_cell_share, temperature): + + # https://h2-mobility.de/wp-content/uploads/2021/02/H2M_Flottenpapier_English_20180822.pdf + # assume in average 1kg_H2 per 100 km -> 1kg_H2 = 33 kWh_H2 (LHV) + # 1 MWh_H2 = 30.003 100km + + car_efficiency = 30.003 # options["transport_fuel_cell_efficiency"] + + # temperature corrected efficiency + efficiency = get_temp_efficency( + car_efficiency, + temperature, + options["transport_heating_deadband_lower"], + options["transport_heating_deadband_upper"], + options["ICE_lower_degree_factor"], + options["ICE_upper_degree_factor"], + ) + + suffix = " land transport fuel cell" + + p_nom = fuel_cell_share * p_set.div(efficiency).max() + + profile = p_set.div(efficiency) / p_set.div(efficiency).max() - if fuel_cell_share > 0: - n.madd( - "Load", - nodes, - suffix=" land transport fuel cell", - bus=nodes + " H2", - carrier="land transport fuel cell", - p_set=fuel_cell_share - / options["transport_fuel_cell_efficiency"] - * transport[nodes], - ) + n.madd( + "Link", + nodes, + suffix=suffix, + bus0=spatial.h2.nodes, + bus1=nodes + " land transport", + carrier="land transport fuel cell", + efficiency=efficiency, + p_nom_extendable=False, + p_nom=p_nom, + p_min_pu=profile, + p_max_pu=profile, + lifetime=1, + ) - if ice_share > 0: - add_carrier_buses(n, "oil") - ice_efficiency = options["transport_internal_combustion_efficiency"] +def add_ice_cars(n, nodes, p_set, ice_share, temperature): - p_set_land_transport_oil = ( - ice_share - / ice_efficiency - * transport[nodes].rename(columns=lambda x: x + " land transport oil") - ) + add_carrier_buses(n, "oil") + + # average consumption 7 liter per 100 km + # 0.008889 MWh_petrol = 1 liter + # 0.062223 MWh_petrol / 100 km + # 1 MWh_petrol = 16.0712 + + car_efficiency = 16.0712 # options["transport_internal_combustion_efficiency"] - if not options["regional_oil_demand"]: - p_set_land_transport_oil = p_set_land_transport_oil.sum(axis=1).to_frame( - name="EU land transport oil" - ) + # temperature corrected efficiency + efficiency = get_temp_efficency( + car_efficiency, + temperature, + options["transport_heating_deadband_lower"], + options["transport_heating_deadband_upper"], + options["ICE_lower_degree_factor"], + options["ICE_upper_degree_factor"], + ) + suffix = " land transport ICE" - n.madd( - "Bus", - spatial.oil.land_transport, - location=spatial.oil.demand_locations, - carrier="land transport oil", - unit="land transport", - ) + p_nom = ice_share * p_set.div(efficiency).max() - n.madd( - "Load", - spatial.oil.land_transport, - bus=spatial.oil.land_transport, - carrier="land transport oil", - p_set=p_set_land_transport_oil, - ) + profile = p_set.div(efficiency) / p_set.div(efficiency).max() - n.madd( - "Link", - spatial.oil.land_transport, - bus0=spatial.oil.nodes, - bus1=spatial.oil.land_transport, - bus2="co2 atmosphere", - carrier="land transport oil", - efficiency2=costs.at["oil", "CO2 intensity"], - p_nom_extendable=True, + n.madd( + "Link", + nodes, + suffix=suffix, + bus0=spatial.oil.nodes, + bus1=nodes + " land transport", + bus2=["co2 atmosphere"], + carrier="land transport oil", + efficiency=efficiency, + efficiency2=costs.at["oil", "CO2 intensity"], + p_nom_extendable=False, + p_nom=p_nom, + p_min_pu=profile, + p_max_pu=profile, + lifetime=1, + ) + + +def adjust_endogenous_transport(n): + + logger.info("Assume endogenous land transport") + carrier = [ + "land transport EV", + "land transport oil", + "land transport fuel cell", + "BEV charger", + "V2G", + ] + + links_i = n.links[n.links.carrier.isin(carrier)].index + n.links.loc[links_i, "p_nom_extendable"] = True + n.links.loc[links_i, "lifetime"] = 15 + + store_carrier = ["battery storage", "Li ion"] + store_i = n.stores[n.stores.carrier.isin(store_carrier)].index + n.stores.loc[store_i, "e_nom_extendable"] = True + + # costs todo + # assume here for all of Europe + # average driving distance 15 000 km /year and car = 150 100km /year + # EV -------------------------- + cost_EV = costs.loc["Battery electric (passenger cars)", "fixed"] / 150 + # FCE ---------------------------- + cost_FCE = costs.loc["Hydrogen fuel cell (passenger cars)", "fixed"] / 150 + # ICE --------------------------------------------------------- + cost_ICE = costs.at["Liquid fuels ICE (passenger cars)", "fixed"] / 150 + # cost in unit input depending on car type + costs_car_type = { + "land transport EV": cost_EV, + "land transport fuel cell": cost_FCE, + "land transport oil": cost_ICE, + } + + # add dummy generator only needed for solving with glpk with higher solver tolerance + n.add("Carrier", "dummy transport", color="#dd2e23", nice_name="Dummy transport") + buses_i = n.buses[n.buses.carrier == "land transport demand"].index + n.madd( + "Generator", + buses_i, + " load", + bus=buses_i, + carrier="load", + marginal_cost=1e9, + p_nom=1e5, + ) + + # n.madd( + # "Generator", + # buses_i, + # " load negative", + # bus=buses_i, + # carrier="load", + # marginal_cost=1e9, + # p_nom=1e5, + # p_max_pu=0, + # p_min_pu=-1, + # sign=-1, + # ) + + for car_type, cost in costs_car_type.items(): + car_i = n.links[n.links.carrier == car_type].index + n.links.loc[car_i, "capital_cost"] = cost + + +def add_land_transport(n, costs): + # TODO options? + + logger.info("Add land transport") + + # read in transport demand in units 100 km + transport = pd.read_csv( + snakemake.input.transport_demand, index_col=[0], header=[0,1], + parse_dates=True + )["light"] + car_cols = ['Number Passenger cars', 'Number Powered 2-wheelers', + 'Number Light duty vehicles', + 'Number Motor coaches, buses and trolley buses', + 'Number Heavy duty vehicles'] + number_cars = pd.read_csv(snakemake.input.transport_data, index_col=0, + )[ + car_cols + ].sum(axis=1) + avail_profile = pd.read_csv( + snakemake.input.avail_profile, index_col=0, parse_dates=True + ) + dsm_profile = pd.read_csv( + snakemake.input.dsm_profile, index_col=0, parse_dates=True + ) + + # exogenous share of passenger car type + engine_types = ["fuel_cell", "electric", "ice"] + shares = pd.Series() + for engine in engine_types: + shares[engine] = get(options[f"land_transport_{engine}_share"], investment_year) + logger.info(f"{engine} share: {shares[engine]*100}%") + + check_land_transport_shares(shares) + + endogenous = options["endogenous_transport"] + nodes = spatial.nodes + + # Add load for transport demand + n.add("Carrier", "land transport demand") + + n.madd( + "Bus", + nodes, + location=nodes, + suffix=" land transport", + carrier="land transport demand", + unit="100 km", + ) + + p_set = transport[nodes] + + # add demand + n.madd( + "Load", + nodes, + suffix=" land transport", + bus=nodes + " land transport", + carrier="land transport demand", + p_set=p_set, + ) + + # temperature for correction factor for heating/cooling + temperature = xr.open_dataarray(snakemake.input.temp_air_total).to_pandas() + + if shares["electric"] > 0 or endogenous: + add_EVs( + n, + nodes, + avail_profile, + dsm_profile, + p_set, + shares["electric"], + number_cars, + temperature, ) + if shares["fuel_cell"] > 0 or endogenous: + add_fuel_cell_cars(n, nodes, p_set, shares["fuel_cell"], temperature) + + if shares["ice"] > 0 or endogenous: + add_ice_cars(n, nodes, p_set, shares["ice"], temperature) + + if endogenous: + adjust_endogenous_transport(n) + + if options['electrobiofuels']: + add_electrobiofuels(n, nodes) + def build_heat_demand(n): heat_demand_shape = ( @@ -2915,7 +3156,7 @@ def add_industry(n, costs): p_set = ( demand_factor - * pop_weighted_energy_totals.loc[nodes, all_aviation].sum(axis=1) + * pop_weighted_energy_totals.loc[nodes, all_aviation].replace(np.inf, 0).sum(axis=1) * 1e6 / nhours ).rename(lambda x: x + " kerosene for aviation") @@ -3560,19 +3801,47 @@ def lossy_bidirectional_links(n, carrier, efficiencies={}): ) +def adjust_transport_temporal_agg(n): + + engine_types = { + "fuel_cell": "land transport fuel cell", + "electric": "land transport EV", + "ice": "land transport oil", + } + + p_set = n.loads_t.p_set.loc[:, n.loads.carrier == "land transport demand"] + + for engine, carrier in engine_types.items(): + + links_i = n.links[n.links.carrier == carrier].index + if links_i.empty: + continue + + share = get(options[f"land_transport_{engine}_share"], investment_year) + efficiency = n.links_t.efficiency.loc[:, links_i] + p_set.columns = efficiency.columns + p_nom = share * p_set.div(efficiency).max() + profile = p_set.div(efficiency) / p_set.div(efficiency).max() + + n.links.loc[links_i, "p_nom"] = p_nom + n.links_t.p_max_pu[links_i] = profile + n.links_t.p_min_pu[links_i] = profile + + +# %% if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake snakemake = mock_snakemake( "prepare_sector_network", - configfiles="test/config.overnight.yaml", + # configfiles="test/config.overnight.yaml", simpl="", opts="", clusters="37", ll="v1.0", - sector_opts="CO2L0-24H-T-H-B-I-A-dist1", - planning_horizons="2030", + sector_opts="25sn-T-H-B-I-A-dist1", + planning_horizons="2050", ) configure_logging(snakemake) @@ -3658,6 +3927,8 @@ def lossy_bidirectional_links(n, carrier, efficiencies={}): resolution = snakemake.params.time_resolution n = set_temporal_aggregation(n, resolution, solver_name) + adjust_transport_temporal_agg(n) + co2_budget = snakemake.params.co2_budget if isinstance(co2_budget, str) and co2_budget.startswith("cb"): fn = "results/" + snakemake.params.RDIR + "/csvs/carbon_budget_distribution.csv" diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 7e53e6063..6629571a4 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -353,7 +353,7 @@ def prepare_network( buses_i = n.buses.index if not np.isscalar(load_shedding): # TODO: do not scale via sign attribute (use Eur/MWh instead of Eur/kWh) - load_shedding = 1e2 # Eur/kWh + load_shedding = 1e4 # Eur/kWh n.madd( "Generator", @@ -820,6 +820,56 @@ def add_co2_atmosphere_constraint(n, snapshots): n.model.add_constraints(lhs <= rhs, name=f"GlobalConstraint-{name}") +def add_endogenous_transport_constraints(n, snapshots): + """ + Add constraints to relate number of EVs to EV charger, V2G and DSM. + """ + + # get index TODO only extendable + link_ext = n.links[n.links.p_nom_extendable] + ev_i = link_ext[link_ext.carrier == "land transport EV"].index + bev_i = link_ext[link_ext.carrier == "BEV charger"].index + v2g_i = link_ext[link_ext.carrier == "V2G"].index + bev_dsm_i = n.stores[ + (n.stores.carrier == "battery storage") & n.stores.e_nom_extendable + ].index + + if ev_i.empty: + return + # factor + f = ( + n.links.loc[ev_i, "p_nom"] + .rename(n.links.bus0) + .div(n.links.loc[bev_i, "p_nom"].rename(n.links.bus1)) + ) + + # variables + link_p_nom = n.model.variables.Link_p_nom + + # constraint for BEV charger + lhs = link_p_nom.loc[ev_i] - (link_p_nom.loc[bev_i] * f.values) + n.model.add_constraints(lhs == 0, name="p_nom-EV-BEV") + + if not v2g_i.empty: + # constraint for V2G + lhs = link_p_nom.loc[ev_i] - (link_p_nom.loc[v2g_i] * f.values) + n.model.add_constraints(lhs == 0, name="p_nom-EV-V2G") + + if not bev_dsm_i.empty: + # factor + f = ( + n.links.loc[ev_i, "p_nom"] + .rename(n.links.bus0) + .div(n.stores.loc[bev_dsm_i, "e_nom"].rename(n.links.bus1)) + ) + + store_e_nom = n.model.variables.Store_e_nom + + # constraint for DSM + lhs = link_p_nom.loc[ev_i] - (store_e_nom.loc[bev_dsm_i] * f.values) + n.model.add_constraints(lhs == 0, name="e_nom-EV-DSM") + + def extra_functionality(n, snapshots): """ Collects supplementary constraints which will be passed to @@ -855,6 +905,8 @@ def extra_functionality(n, snapshots): else: add_co2_atmosphere_constraint(n, snapshots) + if n.config["sector"]["endogenous_transport"]: + add_endogenous_transport_constraints(n, snapshots) if snakemake.params.custom_extra_functionality: source_path = snakemake.params.custom_extra_functionality assert os.path.exists(source_path), f"{source_path} does not exist" @@ -922,19 +974,20 @@ def solve_network(n, config, solving, **kwargs): return n +# %% if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake snakemake = mock_snakemake( "solve_sector_network", - configfiles="../config/test/config.perfect.yaml", + # configfiles="../config/test/config.perfect.yaml", simpl="", opts="", clusters="37", ll="v1.0", - sector_opts="CO2L0-1H-T-H-B-I-A-dist1", - planning_horizons="2030", + sector_opts="Co2L0-25sn-T-H-B-I-A-dist1", + planning_horizons="2050", ) configure_logging(snakemake) set_scenario_config(snakemake)