From 9717cb0707514fe5cba8b81a6384a0678183ed6c Mon Sep 17 00:00:00 2001 From: Fabian Hofmann Date: Mon, 10 Feb 2025 10:02:06 +0100 Subject: [PATCH] Consistent function scope (#1537) * start removing snakemake from functions * move snakemake objects out of prepare_sector_functions * continue with add_existing_baseyear * continue with prepare_perfect_foresight * apply it to add_electricity * continue with solve network * add pylint check to reusable functions * iron out follow up issues * adjust pylint check --- .github/workflows/test.yaml | 5 + envs/environment.yaml | 1 + scripts/add_brownfield.py | 46 +- scripts/add_electricity.py | 212 ++++++- scripts/add_existing_baseyear.py | 216 ++++--- scripts/prepare_network.py | 12 +- scripts/prepare_perfect_foresight.py | 236 +++++--- scripts/prepare_sector_network.py | 808 +++++++++++++++++++++++---- scripts/solve_network.py | 248 +++++--- 9 files changed, 1420 insertions(+), 364 deletions(-) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 6c06ee118..610c23b96 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -85,6 +85,11 @@ jobs: conda env update -n pypsa-eur -f ${{ env.env_file }} echo "Run conda list" && conda list + - name: Run pylint check on scripts + # check for undefined variables to reuse functions across scripts + run: | + pylint --disable=all --enable=E0601 --output-format=parseable scripts/add_* scripts/prepare_* scripts/solve_* + - name: Run snakemake test workflows run: | make test diff --git a/envs/environment.yaml b/envs/environment.yaml index e46a59639..4de83565d 100644 --- a/envs/environment.yaml +++ b/envs/environment.yaml @@ -59,6 +59,7 @@ dependencies: - ipython - pre-commit - ruff +- pylint - pip: - gurobipy diff --git a/scripts/add_brownfield.py b/scripts/add_brownfield.py index 0e91f7571..50a9503f3 100644 --- a/scripts/add_brownfield.py +++ b/scripts/add_brownfield.py @@ -23,7 +23,32 @@ idx = pd.IndexSlice -def add_brownfield(n, n_p, year): +def add_brownfield( + n, + n_p, + year, + h2_retrofit=False, + h2_retrofit_capacity_per_ch4=None, + capacity_threshold=None, +): + """ + Add brownfield capacity from previous network. + + Parameters + ---------- + n : pypsa.Network + Network to add brownfield to + n_p : pypsa.Network + Previous network to get brownfield from + year : int + Planning year + h2_retrofit : bool + Whether to allow hydrogen pipeline retrofitting + h2_retrofit_capacity_per_ch4 : float + Ratio of hydrogen to methane capacity for pipeline retrofitting + capacity_threshold : float + Threshold for removing assets with low capacity + """ logger.info(f"Preparing brownfield for the year {year}") # electric transmission grid set optimised capacities of previous as minimum @@ -49,11 +74,9 @@ def add_brownfield(n, n_p, year): & c.df.index.str.contains("heat") ] - threshold = snakemake.params.threshold_capacity - if not chp_heat.empty: threshold_chp_heat = ( - threshold + capacity_threshold * c.df.efficiency[chp_heat.str.replace("heat", "electric")].values * c.df.p_nom_ratio[chp_heat.str.replace("heat", "electric")].values / c.df.efficiency[chp_heat].values @@ -67,7 +90,7 @@ def add_brownfield(n, n_p, year): c.name, c.df.index[ (c.df[f"{attr}_nom_extendable"] & ~c.df.index.isin(chp_heat)) - & (c.df[f"{attr}_nom_opt"] < threshold) + & (c.df[f"{attr}_nom_opt"] < capacity_threshold) ], ) @@ -85,7 +108,7 @@ def add_brownfield(n, n_p, year): n.import_series_from_dataframe(c.pnl[tattr], c.name, tattr) # deal with gas network - if snakemake.params.H2_retrofit: + if h2_retrofit: # subtract the already retrofitted from the maximum capacity h2_retrofitted_fixed_i = n.links[ (n.links.carrier == "H2 pipeline retrofitted") @@ -118,7 +141,7 @@ def add_brownfield(n, n_p, year): pipe_capacity = n.links.loc[gas_pipes_i, "p_nom"] fr = "H2 pipeline retrofitted" to = "gas pipeline" - CH4_per_H2 = 1 / snakemake.params.H2_retrofit_capacity_per_CH4 + CH4_per_H2 = 1 / h2_retrofit_capacity_per_ch4 already_retrofitted.index = already_retrofitted.index.str.replace(fr, to) remaining_capacity = ( pipe_capacity @@ -301,7 +324,14 @@ def update_heat_pump_efficiency(n: pypsa.Network, n_p: pypsa.Network, year: int) update_heat_pump_efficiency(n, n_p, year) - add_brownfield(n, n_p, year) + add_brownfield( + n, + n_p, + year, + h2_retrofit=snakemake.params.H2_retrofit, + h2_retrofit_capacity_per_ch4=snakemake.params.H2_retrofit_capacity_per_CH4, + capacity_threshold=snakemake.params.threshold_capacity, + ) disable_grid_expansion_if_limit_hit(n) diff --git a/scripts/add_electricity.py b/scripts/add_electricity.py index 35f797c22..9f27fc9dd 100755 --- a/scripts/add_electricity.py +++ b/scripts/add_electricity.py @@ -70,15 +70,46 @@ logger = logging.getLogger(__name__) -def normed(s): +def normed(s: pd.Series) -> pd.Series: + """ + Normalize a pandas Series by dividing each element by the sum of all elements. + + Parameters + ---------- + s : pd.Series + Input series to normalize + + Returns + ------- + pd.Series + Normalized series where all elements sum to 1 + """ return s / s.sum() -def calculate_annuity(n, r): +def calculate_annuity(n: float, r: float | pd.Series) -> float | pd.Series: """ - Calculate the annuity factor for an asset with lifetime n years and. + Calculate the annuity factor for an asset with lifetime n years and discount rate r. + + The annuity factor is used to calculate the annual payment required to pay off a loan + over n years at interest rate r. For example, annuity(20, 0.05) * 20 = 1.6. - discount rate of r, e.g. annuity(20, 0.05) * 20 = 1.6 + Parameters + ---------- + n : float + Lifetime of the asset in years + r : float | pd.Series + Discount rate (interest rate). Can be a single float or a pandas Series of rates. + + Returns + ------- + float | pd.Series + Annuity factor. Returns a float if r is float, or pd.Series if r is pd.Series. + + Examples + -------- + >>> calculate_annuity(20, 0.05) + 0.08024258718774728 """ if isinstance(r, pd.Series): return pd.Series(1 / n, index=r.index).where( @@ -343,6 +374,20 @@ def attach_load( busmap_fn: str, scaling: float = 1.0, ) -> None: + """ + Attach load data to the network. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network to attach the load data to. + load_fn : str + Path to the load data file. + busmap_fn : str + Path to the busmap file. + scaling : float, optional + Scaling factor for the load data, by default 1.0. + """ load = ( xr.open_dataarray(load_fn).to_dataframe().squeeze(axis=1).unstack(level="time") ) @@ -363,6 +408,20 @@ def set_transmission_costs( line_length_factor: float = 1.0, link_length_factor: float = 1.0, ) -> None: + """ + Set the transmission costs for lines and links in the network. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network to set the transmission costs for. + costs : pd.DataFrame + DataFrame containing the cost data. + line_length_factor : float, optional + Factor to scale the line length, by default 1.0. + link_length_factor : float, optional + Factor to scale the link length, by default 1.0. + """ n.lines["capital_cost"] = ( n.lines["length"] * line_length_factor @@ -396,12 +455,32 @@ def set_transmission_costs( def attach_wind_and_solar( n: pypsa.Network, costs: pd.DataFrame, - input_profiles: str, + profile_filenames: dict, carriers: list | set, extendable_carriers: list | set, line_length_factor: float = 1.0, landfall_lengths: dict = None, ) -> None: + """ + Attach wind and solar generators to the network. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network to attach the generators to. + costs : pd.DataFrame + DataFrame containing the cost data. + profile_filenames : dict + Dictionary containing the paths to the wind and solar profiles. + carriers : list | set + List of renewable energy carriers to attach. + extendable_carriers : list | set + List of extendable renewable energy carriers. + line_length_factor : float, optional + Factor to scale the line length, by default 1.0. + landfall_lengths : dict, optional + Dictionary containing the landfall lengths for offshore wind, by default None. + """ add_missing_carriers(n, carriers) if landfall_lengths is None: @@ -413,7 +492,7 @@ def attach_wind_and_solar( landfall_length = landfall_lengths.get(car, 0.0) - with xr.open_dataset(getattr(input_profiles, "profile_" + car)) as ds: + with xr.open_dataset(profile_filenames["profile_" + car]) as ds: if ds.indexes["bus"].empty: continue @@ -460,16 +539,40 @@ def attach_wind_and_solar( def attach_conventional_generators( - n, - costs, - ppl, - conventional_carriers, - extendable_carriers, - conventional_params, - conventional_inputs, - unit_commitment=None, - fuel_price=None, + n: pypsa.Network, + costs: pd.DataFrame, + ppl: pd.DataFrame, + conventional_carriers: list, + extendable_carriers: dict, + conventional_params: dict, + conventional_inputs: dict, + unit_commitment: pd.DataFrame = None, + fuel_price: pd.DataFrame = None, ): + """ + Attach conventional generators to the network. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network to attach the generators to. + costs : pd.DataFrame + DataFrame containing the cost data. + ppl : pd.DataFrame + DataFrame containing the power plant data. + conventional_carriers : list + List of conventional energy carriers. + extendable_carriers : dict + Dictionary of extendable energy carriers. + conventional_params : dict + Dictionary of conventional generator parameters. + conventional_inputs : dict + Dictionary of conventional generator inputs. + unit_commitment : pd.DataFrame, optional + DataFrame containing unit commitment data, by default None. + fuel_price : pd.DataFrame, optional + DataFrame containing fuel price data, by default None. + """ carriers = list(set(conventional_carriers) | set(extendable_carriers["Generator"])) ppl = ppl.query("carrier in @carriers") @@ -532,7 +635,7 @@ def attach_conventional_generators( # Values affecting generators of technology k country-specific # First map generator buses to countries; then map countries to p_max_pu values = pd.read_csv( - snakemake.input[f"conventional_{carrier}_{attr}"], index_col=0 + conventional_inputs[f"conventional_{carrier}_{attr}"], index_col=0 ).iloc[:, 0] bus_values = n.buses.country.map(values) n.generators.update( @@ -543,7 +646,35 @@ def attach_conventional_generators( n.generators.loc[idx, attr] = values -def attach_hydro(n, costs, ppl, profile_hydro, hydro_capacities, carriers, **params): +def attach_hydro( + n: pypsa.Network, + costs: pd.DataFrame, + ppl: pd.DataFrame, + profile_hydro: str, + hydro_capacities: str, + carriers: list, + **params, +): + """ + Attach hydro generators and storage units to the network. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network to attach the hydro units to. + costs : pd.DataFrame + DataFrame containing the cost data. + ppl : pd.DataFrame + DataFrame containing the power plant data. + profile_hydro : str + Path to the hydro profile data. + hydro_capacities : str + Path to the hydro capacities data. + carriers : list + List of hydro energy carriers. + **params : + Additional parameters for hydro units. + """ add_missing_carriers(n, carriers) add_co2_emissions(n, costs, carriers) @@ -716,8 +847,12 @@ def attach_OPSD_renewables(n: pypsa.Network, tech_map: dict[str, list[str]]) -> def estimate_renewable_capacities( - n: pypsa.Network, year: int, tech_map: dict, expansion_limit: bool, countries: list -) -> None: + n: pypsa.Network, + year: int, + tech_map: dict, + expansion_limit: bool, + countries: list, +): """ Estimate a different between renewable capacities in the network and reported country totals from IRENASTAT dataset. Distribute the difference @@ -778,7 +913,26 @@ def estimate_renewable_capacities( ) -def attach_storageunits(n, costs, extendable_carriers, max_hours): +def attach_storageunits( + n: pypsa.Network, + costs: pd.DataFrame, + extendable_carriers: dict, + max_hours: dict, +): + """ + Attach storage units to the network. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network to attach the storage units to. + costs : pd.DataFrame + DataFrame containing the cost data. + extendable_carriers : dict + Dictionary of extendable energy carriers. + max_hours : dict + Dictionary of maximum hours for storage units. + """ carriers = extendable_carriers["StorageUnit"] n.add("Carrier", carriers) @@ -809,7 +963,23 @@ def attach_storageunits(n, costs, extendable_carriers, max_hours): ) -def attach_stores(n, costs, extendable_carriers): +def attach_stores( + n: pypsa.Network, + costs: pd.DataFrame, + extendable_carriers: dict, +): + """ + Attach stores to the network. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network to attach the stores to. + costs : pd.DataFrame + DataFrame containing the cost data. + extendable_carriers : dict + Dictionary of extendable energy carriers. + """ carriers = extendable_carriers["Store"] n.add("Carrier", carriers) diff --git a/scripts/add_existing_baseyear.py b/scripts/add_existing_baseyear.py index 91e01939e..fa6b78fb6 100644 --- a/scripts/add_existing_baseyear.py +++ b/scripts/add_existing_baseyear.py @@ -31,13 +31,16 @@ spatial = SimpleNamespace() -def add_build_year_to_new_assets(n, baseyear): +def add_build_year_to_new_assets(n: pypsa.Network, baseyear: int) -> None: """ + Add build year to new assets in the network. + Parameters ---------- n : pypsa.Network + Network to modify baseyear : int - year in which optimized assets are built + Year in which optimized assets are built """ # Give assets with lifetimes and no build year the build year baseyear for c in n.iterate_components(["Link", "Generator", "Store"]): @@ -57,14 +60,33 @@ def add_build_year_to_new_assets(n, baseyear): c.pnl[attr] = c.pnl[attr].rename(columns=rename) -def add_existing_renewables(df_agg, costs): +def add_existing_renewables( + n: pypsa.Network, + costs: pd.DataFrame, + df_agg: pd.DataFrame, + countries: list[str], +) -> None: """ - Append existing renewables to the df_agg pd.DataFrame with the conventional - power plants. + Add existing renewable capacities to conventional power plant data. + + Parameters + ---------- + df_agg : pd.DataFrame + DataFrame containing conventional power plant data + costs : pd.DataFrame + Technology cost data with 'lifetime' column indexed by technology + n : pypsa.Network + Network containing topology and generator data + countries : list + List of country codes to consider + + Returns + ------- + None + Modifies df_agg in-place """ tech_map = {"solar": "PV", "onwind": "Onshore", "offwind-ac": "Offshore"} - countries = snakemake.config["countries"] # noqa: F841 irena = pm.data.IRENASTAT().powerplant.convert_country_to_alpha2() irena = irena.query("Country in @countries") irena = irena.groupby(["Technology", "Country", "Year"]).Capacity.sum() @@ -121,23 +143,41 @@ def add_existing_renewables(df_agg, costs): df_agg.at[name, "bus"] = node -def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, baseyear): +def add_power_capacities_installed_before_baseyear( + n: pypsa.Network, + costs: pd.DataFrame, + grouping_years: list[int], + baseyear: int, + powerplants_file: str, + countries: list[str], + capacity_threshold: float, + lifetime_values: dict[str, float], +) -> None: """ + Add power generation capacities installed before base year. + Parameters ---------- n : pypsa.Network - grouping_years : - intervals to group existing capacities - costs : - to read lifetime to estimate YearDecomissioning + Network to modify + costs : pd.DataFrame + Technology costs + grouping_years : list + Intervals to group existing capacities baseyear : int + Base year for analysis + powerplants_file : str + Path to powerplants CSV file + countries : list + List of countries to consider + capacity_threshold : float + Minimum capacity threshold + lifetime_values : dict + Default values for missing data """ - logger.debug( - f"Adding power capacities installed before {baseyear} from" - " powerplants_s_{clusters}.csv" - ) + logger.debug(f"Adding power capacities installed before {baseyear}") - df_agg = pd.read_csv(snakemake.input.powerplants, index_col=0) + df_agg = pd.read_csv(powerplants_file, index_col=0) rename_fuel = { "Hard Coal": "coal", @@ -177,15 +217,16 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas mean = df_agg.loc[biomass_i, "DateIn"].mean() df_agg.loc[biomass_i, "DateIn"] = df_agg.loc[biomass_i, "DateIn"].fillna(int(mean)) # Fill missing DateOut - dateout = ( - df_agg.loc[biomass_i, "DateIn"] - + snakemake.params.costs["fill_values"]["lifetime"] - ) + dateout = df_agg.loc[biomass_i, "DateIn"] + lifetime_values["lifetime"] df_agg.loc[biomass_i, "DateOut"] = df_agg.loc[biomass_i, "DateOut"].fillna(dateout) # include renewables in df_agg - add_existing_renewables(df_agg, costs) - + add_existing_renewables( + df_agg=df_agg, + costs=costs, + n=n, + countries=countries, + ) # drop assets which are already phased out / decommissioned phased_out = df_agg[df_agg["DateOut"] < baseyear].index df_agg.drop(phased_out, inplace=True) @@ -237,9 +278,7 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas # capacity is the capacity in MW at each node for this capacity = df.loc[grouping_year, generator] capacity = capacity[~capacity.isna()] - capacity = capacity[ - capacity > snakemake.params.existing_capacities["threshold_capacity"] - ] + capacity = capacity[capacity > capacity_threshold] suffix = "-ac" if generator == "offwind" else "" name_suffix = f" {generator}{suffix}-{grouping_year}" name_suffix_by = f" {generator}{suffix}-{baseyear}" @@ -373,7 +412,13 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas ] -def get_efficiency(heat_system, carrier, nodes, heating_efficiencies, costs): +def get_efficiency( + heat_system: HeatSystem, + carrier: str, + nodes: pd.Index, + efficiencies: dict[str, float], + costs: pd.DataFrame, +) -> pd.Series | float: """ Computes the heating system efficiency based on the sector and carrier type. @@ -385,7 +430,7 @@ def get_efficiency(heat_system, carrier, nodes, heating_efficiencies, costs): The type of fuel or energy carrier (e.g., 'gas', 'oil'). nodes : pandas.Series A pandas Series containing node information used to match the heating efficiency data. - heating_efficiencies : dict + efficiencies : dict A dictionary containing efficiency values for different carriers and sectors. costs : pandas.DataFrame A DataFrame containing boiler cost and efficiency data for different heating systems. @@ -407,10 +452,10 @@ def get_efficiency(heat_system, carrier, nodes, heating_efficiencies, costs): efficiency = costs.at[boiler_costs_name, "efficiency"] elif heat_system.sector.value == "residential": key = f"{carrier} residential space efficiency" - efficiency = nodes.str[:2].map(heating_efficiencies[key]) + efficiency = nodes.str[:2].map(efficiencies[key]) elif heat_system.sector.value == "services": key = f"{carrier} services space efficiency" - efficiency = nodes.str[:2].map(heating_efficiencies[key]) + efficiency = nodes.str[:2].map(efficiencies[key]) else: logger.warning(f"{heat_system} not defined.") @@ -419,41 +464,68 @@ def get_efficiency(heat_system, carrier, nodes, heating_efficiencies, costs): def add_heating_capacities_installed_before_baseyear( n: pypsa.Network, - baseyear: int, - grouping_years: list, - cop: dict, - time_dep_hp_cop: bool, costs: pd.DataFrame, + baseyear: int, + grouping_years: list[int], + existing_capacities: pd.DataFrame, + heat_pump_cop: xr.DataArray, + heat_pump_source_types: dict[str, list[str]], + efficiency_file: str, + use_time_dependent_cop: bool, default_lifetime: int, - existing_heating: pd.DataFrame, -): + energy_totals_year: int, + capacity_threshold: float, + use_electricity_distribution_grid: bool, +) -> None: """ + Add heating capacities installed before base year. + Parameters ---------- n : pypsa.Network - baseyear : last year covered in the existing capacities database - grouping_years : intervals to group existing capacities - linear decommissioning of heating capacities from 2020 to 2045 is - currently assumed heating capacities split between residential and - services proportional to heating load in both 50% capacities - in rural buses 50% in urban buses - cop: xr.DataArray - DataArray with time-dependent coefficients of performance (COPs) heat pumps. Coordinates are heat sources (see config), heat system types (see :file:`scripts/enums/HeatSystemType.py`), nodes and snapshots. - time_dep_hp_cop: bool - If True, time-dependent (dynamic) COPs are used for heat pumps + Network to modify + costs : pd.DataFrame + Technology costs + baseyear : int + Base year for analysis + grouping_years : list + Intervals to group capacities + heat_pump_cop : xr.DataArray + Heat pump coefficients of performance + use_time_dependent_cop : bool + Use time-dependent COPs + heating_default_lifetime : int + Default lifetime for heating systems + existing_capacities : pd.DataFrame + Existing heating capacity distribution + heat_pump_source_types : dict + Heat pump sources by system type + efficiency_file : str + Path to heating efficiencies file + energy_totals_year : int + Year for energy totals + capacity_threshold : float + Minimum capacity threshold + use_electricity_distribution_grid : bool + Whether to use electricity distribution grid """ logger.debug(f"Adding heating capacities installed before {baseyear}") - for heat_system in existing_heating.columns.get_level_values(0).unique(): + # Load heating efficiencies + heating_efficiencies = pd.read_csv(efficiency_file, index_col=[1, 0]).loc[ + energy_totals_year + ] + + for heat_system in existing_capacities.columns.get_level_values(0).unique(): heat_system = HeatSystem(heat_system) nodes = pd.Index( n.buses.location[n.buses.index.str.contains(f"{heat_system} heat")] ) - if (not heat_system == HeatSystem.URBAN_CENTRAL) and options[ - "electricity_distribution_grid" - ]: + if ( + not heat_system == HeatSystem.URBAN_CENTRAL + ) and use_electricity_distribution_grid: nodes_elec = nodes + " low voltage" else: nodes_elec = nodes @@ -485,20 +557,18 @@ def add_heating_capacities_installed_before_baseyear( for ratio, grouping_year in zip(ratios, valid_grouping_years): # Add heat pumps - for heat_source in snakemake.params.heat_pump_sources[ - heat_system.system_type.value - ]: + for heat_source in heat_pump_source_types[heat_system.system_type.value]: costs_name = heat_system.heat_pump_costs_name(heat_source) efficiency = ( - cop.sel( + heat_pump_cop.sel( heat_system=heat_system.system_type.value, heat_source=heat_source, name=nodes, ) .to_pandas() .reindex(index=n.snapshots) - if time_dep_hp_cop + if use_time_dependent_cop else costs.at[costs_name, "efficiency"] ) @@ -512,7 +582,7 @@ def add_heating_capacities_installed_before_baseyear( efficiency=efficiency, capital_cost=costs.at[costs_name, "efficiency"] * costs.at[costs_name, "fixed"], - p_nom=existing_heating.loc[ + p_nom=existing_capacities.loc[ nodes, (heat_system.value, f"{heat_source} heat pump") ] * ratio @@ -537,7 +607,9 @@ def add_heating_capacities_installed_before_baseyear( * costs.at[heat_system.resistive_heater_costs_name, "fixed"] ), p_nom=( - existing_heating.loc[nodes, (heat_system.value, "resistive heater")] + existing_capacities.loc[ + nodes, (heat_system.value, "resistive heater") + ] * ratio / costs.at[heat_system.resistive_heater_costs_name, "efficiency"] ), @@ -564,7 +636,7 @@ def add_heating_capacities_installed_before_baseyear( * costs.at[heat_system.gas_boiler_costs_name, "fixed"] ), p_nom=( - existing_heating.loc[nodes, (heat_system.value, "gas boiler")] + existing_capacities.loc[nodes, (heat_system.value, "gas boiler")] * ratio / costs.at[heat_system.gas_boiler_costs_name, "efficiency"] ), @@ -589,7 +661,7 @@ def add_heating_capacities_installed_before_baseyear( capital_cost=costs.at[heat_system.oil_boiler_costs_name, "efficiency"] * costs.at[heat_system.oil_boiler_costs_name, "fixed"], p_nom=( - existing_heating.loc[nodes, (heat_system.value, "oil boiler")] + existing_capacities.loc[nodes, (heat_system.value, "oil boiler")] * ratio / costs.at[heat_system.oil_boiler_costs_name, "efficiency"] ), @@ -610,13 +682,13 @@ def add_heating_capacities_installed_before_baseyear( ) # delete links with capacities below threshold - threshold = snakemake.params.existing_capacities["threshold_capacity"] n.remove( "Link", [ index for index in n.links.index.to_list() - if str(grouping_year) in index and n.links.p_nom[index] < threshold + if str(grouping_year) in index + and n.links.p_nom[index] < capacity_threshold ], ) @@ -661,7 +733,14 @@ def add_heating_capacities_installed_before_baseyear( grouping_years_power = snakemake.params.existing_capacities["grouping_years_power"] grouping_years_heat = snakemake.params.existing_capacities["grouping_years_heat"] add_power_capacities_installed_before_baseyear( - n, grouping_years_power, costs, baseyear + n=n, + costs=costs, + grouping_years=grouping_years_power, + baseyear=baseyear, + powerplants_file=snakemake.input.powerplants, + countries=snakemake.config["countries"], + capacity_threshold=snakemake.params.existing_capacities["threshold_capacity"], + lifetime_values=snakemake.params.costs["fill_values"], ) if options["heating"]: @@ -672,19 +751,26 @@ def add_heating_capacities_installed_before_baseyear( add_heating_capacities_installed_before_baseyear( n=n, + costs=costs, baseyear=baseyear, grouping_years=grouping_years_heat, - cop=xr.open_dataarray(snakemake.input.cop_profiles), - time_dep_hp_cop=options["time_dep_hp_cop"], - costs=costs, + heat_pump_cop=xr.open_dataarray(snakemake.input.cop_profiles), + use_time_dependent_cop=options["time_dep_hp_cop"], default_lifetime=snakemake.params.existing_capacities[ "default_heating_lifetime" ], - existing_heating=pd.read_csv( + existing_capacities=pd.read_csv( snakemake.input.existing_heating_distribution, header=[0, 1], index_col=0, ), + heat_pump_source_types=snakemake.params.heat_pump_sources, + efficiency_file=snakemake.input.heating_efficiencies, + energy_totals_year=snakemake.params["energy_totals_year"], + capacity_threshold=snakemake.params.existing_capacities[ + "threshold_capacity" + ], + use_electricity_distribution_grid=options["electricity_distribution_grid"], ) # Set defaults for missing missing values diff --git a/scripts/prepare_network.py b/scripts/prepare_network.py index 175558bee..4ac98a4b7 100755 --- a/scripts/prepare_network.py +++ b/scripts/prepare_network.py @@ -119,8 +119,8 @@ def add_emission_prices(n, emission_prices={"co2": 0.0}, exclude_co2=False): n.storage_units["marginal_cost"] += su_ep -def add_dynamic_emission_prices(n): - co2_price = pd.read_csv(snakemake.input.co2_price, index_col=0, parse_dates=True) +def add_dynamic_emission_prices(n, fn): + co2_price = pd.read_csv(fn, index_col=0, parse_dates=True) co2_price = co2_price[~co2_price.index.duplicated()] co2_price = co2_price.reindex(n.snapshots).ffill().bfill() @@ -182,13 +182,13 @@ def set_transmission_limit(n, kind, factor, costs, Nyears=1): return n -def average_every_nhours(n, offset): +def average_every_nhours(n, offset, drop_leap_day=False): logger.info(f"Resampling the network to {offset}") m = n.copy(with_time=False) snapshot_weightings = n.snapshot_weightings.resample(offset).sum() sns = snapshot_weightings.index - if snakemake.params.drop_leap_day: + if drop_leap_day: sns = sns[~((sns.month == 2) & (sns.day == 29))] m.set_snapshots(snapshot_weightings.index) m.snapshot_weightings = snapshot_weightings @@ -315,7 +315,7 @@ def set_line_nom_max( time_resolution = snakemake.params.time_resolution is_string = isinstance(time_resolution, str) if is_string and time_resolution.lower().endswith("h"): - n = average_every_nhours(n, time_resolution) + n = average_every_nhours(n, time_resolution, snakemake.params.drop_leap_day) # segments with package tsam if is_string and time_resolution.lower().endswith("seg"): @@ -336,7 +336,7 @@ def set_line_nom_max( logger.info( "Setting time dependent emission prices according spot market price" ) - add_dynamic_emission_prices(n) + add_dynamic_emission_prices(n, snakemake.input.co2_price) elif emission_prices["enable"]: add_emission_prices( n, dict(co2=snakemake.params.costs["emission_prices"]["co2"]) diff --git a/scripts/prepare_perfect_foresight.py b/scripts/prepare_perfect_foresight.py index 91c09827e..90d5edc47 100644 --- a/scripts/prepare_perfect_foresight.py +++ b/scripts/prepare_perfect_foresight.py @@ -25,37 +25,63 @@ # helper functions --------------------------------------------------- -def get_missing(df, n, c): +def get_missing(df: pd.DataFrame, n: pypsa.Network, c: str) -> pd.DataFrame: """ - Get in network n missing assets of df for component c. + Get missing assets in network n compared to df for component c. - Input: - df: pandas DataFrame, static values of pypsa components - n : pypsa Network to which new assets should be added - c : string, pypsa component.list_name (e.g. "generators") + Parameters + ---------- + df : pd.DataFrame + Static values of pypsa components + n : pypsa.Network + Network to which new assets should be added + c : str + pypsa component.list_name (e.g. "generators") - Return: - pd.DataFrame with static values of missing assets + Returns + ------- + pd.DataFrame + Static values of missing assets """ df_final = getattr(n, c) missing_i = df.index.difference(df_final.index) return df.loc[missing_i] -def get_social_discount(t, r=0.01): +def get_social_discount(t: int, r: float = 0.01) -> float: """ - Calculate for a given time t and social discount rate r [per unit] the - social discount. + Calculate social discount for given time and rate. + + Parameters + ---------- + t : int + Time period in years + r : float, default 0.01 + Social discount rate per unit + + Returns + ------- + float + Social discount factor """ return 1 / (1 + r) ** t -def get_investment_weighting(time_weighting, r=0.01): +def get_investment_weighting(time_weighting: pd.Series, r: float = 0.01) -> pd.Series: """ - Define cost weighting. + Define cost weighting for investment periods. - Returns cost weightings depending on the the time_weighting - (pd.Series) and the social discountrate r + Parameters + ---------- + time_weighting : pd.Series + Time weightings for each period + r : float, default 0.01 + Social discount rate per unit + + Returns + ------- + pd.Series + Cost weightings for each investment period """ end = time_weighting.cumsum() start = time_weighting.cumsum().shift().fillna(0) @@ -67,7 +93,7 @@ def get_investment_weighting(time_weighting, r=0.01): ) -def add_year_to_constraints(n, baseyear): +def add_year_to_constraints(n: pypsa.Network, baseyear: int) -> None: """ Add investment period to global constraints and rename index. @@ -83,7 +109,7 @@ def add_year_to_constraints(n, baseyear): c.df.rename(index=lambda x: x + "-" + str(baseyear), inplace=True) -def hvdc_transport_model(n): +def hvdc_transport_model(n: pypsa.Network) -> None: """ Convert AC lines to DC links for multi-decade optimisation with line expansion. @@ -120,18 +146,23 @@ def hvdc_transport_model(n): ) -def adjust_electricity_grid(n, year, years): +def adjust_electricity_grid(n: pypsa.Network, year: int, years: list[int]) -> None: """ - Add carrier to lines. Replace AC lines with DC links in case of line - expansion. Add lifetime to DC links in case of line expansion. + Adjust electricity grid for multi-decade optimization. Parameters ---------- - n : pypsa.Network + n : pypsa.Network + Network to adjust year : int - year in which optimized assets are built - years: list - investment periods + Year in which optimized assets are built + years : list[int] + List of investment periods + + Returns + ------- + None + Modifies network in place """ n.lines["carrier"] = "AC" links_i = n.links[n.links.carrier == "DC"].index @@ -145,22 +176,24 @@ def adjust_electricity_grid(n, year, years): # -------------------------------------------------------------------- -def concat_networks(years): +def concat_networks(years: list[int], network_paths: list[str]) -> pypsa.Network: """ - Concat given pypsa networks and adds build_year. + Concat given pypsa networks and add build years. - Return: - n : pypsa.Network for the whole planning horizon - """ + Parameters + ---------- + years : list[int] + List of years representing investment periods + network_paths : list[str] + List of paths to network files for each investment period - # input paths of sector coupling networks - network_paths = [snakemake.input.brownfield_network] + [ - snakemake.input[f"network_{year}"] for year in years[1:] - ] - # final concatenated network + Returns + ------- + pypsa.Network + Network for the whole planning horizon + """ n = pypsa.Network() - # iterate over single year networks and concat to perfect foresight network for i, network_path in enumerate(network_paths): year = years[i] network = pypsa.Network(network_path) @@ -241,10 +274,22 @@ def concat_networks(years): return n -def adjust_stores(n): +def adjust_stores(n: pypsa.Network) -> None: """ - Make sure that stores still behave cyclic over one year and not whole - modelling horizon. + Adjust store behavior for multi-year optimization. + + Ensures stores behave cyclically over one year and not whole modeling horizon. + Sets appropriate flags for different types of stores. + + Parameters + ---------- + n : pypsa.Network + Network to adjust + + Returns + ------- + pypsa.Network + Network with adjusted store settings """ # cyclic constraint cyclic_i = n.stores[n.stores.e_cyclic].index @@ -265,13 +310,28 @@ def adjust_stores(n): co2_i = n.stores[n.stores.carrier.isin(e_initial_store)].index n.stores.loc[co2_i, "e_initial_per_period"] = True - return n - -def set_phase_out(n, carrier, ct, phase_out_year): +def set_phase_out( + n: pypsa.Network, carrier: list[str], ct: str, phase_out_year: int +) -> None: """ - Set planned phase outs for given carrier,country (ct) and planned year of - phase out (phase_out_year). + Set planned phase outs for given assets. + + Parameters + ---------- + n : pypsa.Network + Network to adjust + carrier : list[str] + List of carrier names to phase out + ct : str + Two-letter country code + phase_out_year : int + Year when phase out should be complete + + Returns + ------- + None + Modifies network in place """ df = n.links[(n.links.carrier.isin(carrier)) & (n.links.bus1.str[:2] == ct)] # assets which are going to be phased out before end of their lifetime @@ -281,7 +341,7 @@ def set_phase_out(n, carrier, ct, phase_out_year): n.links.loc[assets_i, "lifetime"] = (phase_out_year - build_year).astype(float) -def set_all_phase_outs(n): +def set_all_phase_outs(n: pypsa.Network) -> None: # TODO move this to a csv or to the config planned = [ (["nuclear"], "DE", 2022), @@ -308,13 +368,28 @@ def set_all_phase_outs(n): n.remove("Link", remove_i) -def set_carbon_constraints(n): +def set_carbon_constraints( + n: pypsa.Network, co2_budget: float | None, sector_opts: str +) -> None: """ Add global constraints for carbon emissions. + + Parameters + ---------- + n : pypsa.Network + Network to add constraints to + co2_budget : float or None + CO2 budget in Gt CO2; if float, converted to t CO2 + sector_opts : str + String of sector options separated by "-" + + Returns + ------- + pypsa.Network + Network with carbon constraints added """ - budget = snakemake.config["co2_budget"] - if budget and isinstance(budget, float): - budget *= 1e9 # convert to t CO2 + if co2_budget and isinstance(co2_budget, float): + budget = co2_budget * 1e9 # convert to t CO2 logger.info(f"add carbon budget of {budget}") n.add( @@ -342,7 +417,7 @@ def set_carbon_constraints(n): ) # set minimum CO2 emission constraint to avoid too fast reduction - if "co2min" in snakemake.wildcards.sector_opts.split("-"): + if "co2min" in sector_opts.split("-"): emissions_1990 = 4.53693 emissions_2019 = 3.344096 target_2030 = 0.45 * emissions_1990 @@ -361,10 +436,8 @@ def set_carbon_constraints(n): constant=co2min * 1e9 * time_weightings, ) - return n - -def adjust_lvlimit(n): +def adjust_lvlimit(n: pypsa.Network) -> None: """ Convert global constraints for single investment period to one uniform if all attributes stay the same. @@ -379,20 +452,16 @@ def adjust_lvlimit(n): n.remove(c, remove_i) n.add(c, glc.index, **glc) - return n - -def adjust_CO2_glc(n): +def adjust_CO2_glc(n: pypsa.Network) -> None: c = "GlobalConstraint" glc_name = "CO2Limit" glc_type = "primary_energy" mask = (n.df(c).index.str.contains(glc_name)) & (n.df(c).type == glc_type) n.df(c).loc[mask, "type"] = "co2_limit" - return n - -def add_H2_boilers(n): +def add_H2_boilers(n: pypsa.Network) -> None: """ Gas boilers can be retrofitted to run with H2. @@ -423,18 +492,24 @@ def add_H2_boilers(n): def apply_time_segmentation_perfect( - n, segments, solver_name="cbc", overwrite_time_dependent=True -): + n: pypsa.Network, segments: int, solver_name: str = "cbc" +) -> None: """ - Aggregating time series to segments with different lengths. - - Input: - n: pypsa Network - segments: (int) number of segments in which the typical period should be - subdivided - solver_name: (str) name of solver - overwrite_time_dependent: (bool) overwrite time dependent data of pypsa network - with typical time series created by tsam + Aggregate time series to segments with different lengths. + + Parameters + ---------- + n : pypsa.Network + Network to segment + segments : int + Number of segments for typical period subdivision + solver_name : str, default "cbc" + Name of solver to use for segmentation + + Returns + ------- + pypsa.Network + Network with segmented time series """ try: import tsam.timeseriesaggregation as tsam @@ -485,10 +560,8 @@ def apply_time_segmentation_perfect( n.set_snapshots(sn_weightings.index) n.snapshot_weightings = n.snapshot_weightings.mul(sn_weightings, axis=0) - return n - -def update_heat_pump_efficiency(n: pypsa.Network, years: list[int]): +def update_heat_pump_efficiency(n: pypsa.Network, years: list[int]) -> None: """ Update the efficiency of heat pumps from previous years to current year (e.g. 2030 heat pumps receive 2040 heat pump COPs in 2030). @@ -545,20 +618,23 @@ def update_heat_pump_efficiency(n: pypsa.Network, years: list[int]): ) # concat prepared networks of planning horizon to single network ------------ - n = concat_networks(years) + network_paths = [snakemake.input.brownfield_network] + [ + snakemake.input[f"network_{year}"] for year in years[1:] + ] + n = concat_networks(years, network_paths) # temporal aggregate solver_name = snakemake.config["solving"]["solver"]["name"] segments = snakemake.params.time_resolution if isinstance(segments, (int, float)): - n = apply_time_segmentation_perfect(n, segments, solver_name=solver_name) + apply_time_segmentation_perfect(n, segments, solver_name=solver_name) # adjust global constraints lv limit if the same for all years - n = adjust_lvlimit(n) + adjust_lvlimit(n) # adjust global constraints CO2 limit - n = adjust_CO2_glc(n) + adjust_CO2_glc(n) # adjust stores to multi period investment - n = adjust_stores(n) + adjust_stores(n) # set phase outs set_all_phase_outs(n) @@ -567,7 +643,11 @@ def update_heat_pump_efficiency(n: pypsa.Network, years: list[int]): add_H2_boilers(n) # set carbon constraints - n = set_carbon_constraints(n) + set_carbon_constraints( + n, + co2_budget=snakemake.config["co2_budget"], + sector_opts=snakemake.wildcards.sector_opts, + ) # update meta n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 552bdc0e3..84fb562e7 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -265,7 +265,16 @@ def co2_emissions_year( # TODO: move to own rule with sector-opts wildcard? -def build_carbon_budget(o, input_eurostat, fn, emissions_scope, input_co2, options): +def build_carbon_budget( + o, + input_eurostat, + fn, + emissions_scope, + input_co2, + options, + countries, + planning_horizons, +): """ Distribute carbon budget following beta or exponential transition path. """ @@ -279,8 +288,6 @@ def build_carbon_budget(o, input_eurostat, fn, emissions_scope, input_co2, optio carbon_budget = float(o[o.find("cb") + 2 : o.find("ex")]) r = float(o[o.find("ex") + 2 :]) - countries = snakemake.params.countries - e_1990 = co2_emissions_year( countries, input_eurostat, @@ -300,7 +307,6 @@ def build_carbon_budget(o, input_eurostat, fn, emissions_scope, input_co2, optio year=2018, ) - planning_horizons = snakemake.params.planning_horizons if not isinstance(planning_horizons, list): planning_horizons = [planning_horizons] t_0 = planning_horizons[0] @@ -407,12 +413,27 @@ def make_index(c): def update_wind_solar_costs( n: pypsa.Network, costs: pd.DataFrame, - line_length_factor: int | float = 1, + profiles: dict[str, str], landfall_lengths: dict = None, + line_length_factor: int | float = 1, ) -> None: """ Update costs for wind and solar generators added with pypsa-eur to those cost in the planning year. + + Parameters + ---------- + n : pypsa.Network + Network to update generator costs + costs : pd.DataFrame + Cost assumptions DataFrame + line_length_factor : int | float, optional + Factor to multiply line lengths by, by default 1 + landfall_lengths : dict, optional + Dictionary of landfall lengths per technology, by default None + profiles : dict[str, str] + Dictionary mapping technology names to profile file paths + e.g. {'offwind-dc': 'path/to/profile.nc'} """ if landfall_lengths is None: @@ -429,13 +450,14 @@ def update_wind_solar_costs( ] # for offshore wind, need to calculated connection costs - for connection in ["dc", "ac", "float"]: - tech = "offwind-" + connection + for key, fn in profiles.items(): + tech = key[len("profile_") :] landfall_length = landfall_lengths.get(tech, 0.0) + if tech not in n.generators.carrier.values: continue - profile = snakemake.input["profile_offwind-" + connection] - with xr.open_dataset(profile) as ds: + + with xr.open_dataset(fn) as ds: # if-statement for compatibility with old profiles if "year" in ds.indexes: ds = ds.sel(year=ds.year.min(), drop=True) @@ -555,14 +577,22 @@ def add_carrier_buses(n, carrier, nodes=None): # TODO: PyPSA-Eur merge issue -def remove_elec_base_techs(n): +def remove_elec_base_techs(n: pypsa.Network, carriers_to_keep: dict) -> None: """ Remove conventional generators (e.g. OCGT) and storage units (e.g. batteries and H2) from base electricity-only network, since they're added here differently using links. + + Parameters + ---------- + n : pypsa.Network + Network to remove components from + carriers_to_keep : dict + Dictionary specifying which carriers to keep for each component type + e.g. {'Generator': ['hydro'], 'StorageUnit': ['PHS']} """ - for c in n.iterate_components(snakemake.params.pypsa_eur): - to_keep = snakemake.params.pypsa_eur[c.name] + for c in n.iterate_components(carriers_to_keep): + to_keep = carriers_to_keep[c.name] to_remove = pd.Index(c.df.carrier.unique()).symmetric_difference(to_keep) if to_remove.empty: continue @@ -582,10 +612,12 @@ def remove_non_electric_buses(n): n.buses = n.buses[n.buses.carrier.isin(["AC", "DC"])] -def patch_electricity_network(n, costs, landfall_lengths): - remove_elec_base_techs(n) +def patch_electricity_network(n, costs, carriers_to_keep, profiles, landfall_lengths): + remove_elec_base_techs(n, carriers_to_keep) remove_non_electric_buses(n) - update_wind_solar_costs(n, costs, landfall_lengths=landfall_lengths) + update_wind_solar_costs( + n, costs, landfall_lengths=landfall_lengths, profiles=profiles + ) n.loads["carrier"] = "electricity" n.buses["location"] = n.buses.index n.buses["unit"] = "MWh_el" @@ -605,7 +637,44 @@ def add_eu_bus(n, x=-5.5, y=46): n.add("Carrier", "none") -def add_co2_tracking(n, costs, options): +def add_co2_tracking(n, costs, options, sequestration_potential_file=None): + """ + Add CO2 tracking components to the network including atmospheric CO2, + CO2 storage, and sequestration infrastructure. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network container object + costs : pd.DataFrame + Cost assumptions for different technologies, must include + 'CO2 storage tank' with 'fixed' cost column + options : dict + Configuration options containing at least: + - regional_co2_sequestration_potential: dict with keys + - enable: bool + - max_size: float + - years_of_storage: float + - co2_sequestration_cost: float + - co2_sequestration_lifetime: float + - co2_vent: bool + sequestration_potential_file : str, optional + Path to CSV file containing regional CO2 sequestration potentials. + Required if options['regional_co2_sequestration_potential']['enable'] is True. + + Returns + ------- + None + Modifies the network object in-place by adding CO2-related components. + + Notes + ----- + Adds several components to track CO2: + - Atmospheric CO2 store + - CO2 storage tanks + - CO2 sequestration infrastructure + - Optional CO2 venting facilities + """ # minus sign because opposite to how fossil fuels used: # CH4 burning puts CH4 down, atmosphere up n.add("Carrier", "co2", co2_emissions=-1.0) @@ -666,13 +735,16 @@ def add_co2_tracking(n, costs, options): ) if options["regional_co2_sequestration_potential"]["enable"]: + if sequestration_potential_file is None: + raise ValueError( + "sequestration_potential_file must be provided when " + "regional_co2_sequestration_potential is enabled" + ) upper_limit = ( options["regional_co2_sequestration_potential"]["max_size"] * 1e3 ) # Mt annualiser = options["regional_co2_sequestration_potential"]["years_of_storage"] - e_nom_max = pd.read_csv( - snakemake.input.sequestration_potential, index_col=0 - ).squeeze() + e_nom_max = pd.read_csv(sequestration_potential_file, index_col=0).squeeze() e_nom_max = ( e_nom_max.reindex(spatial.co2.locations) .fillna(0.0) @@ -710,7 +782,36 @@ def add_co2_tracking(n, costs, options): ) -def add_co2_network(n, costs): +def add_co2_network(n, costs, co2_network_cost_factor=1.0): + """ + Add CO2 transport network to the PyPSA network. + + Creates a CO2 pipeline network with both onshore and submarine pipeline segments, + considering different costs for each type. The network allows bidirectional flow + and is extendable. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network container object + costs : pd.DataFrame + Cost assumptions for different technologies. Must contain entries for + 'CO2 pipeline' and 'CO2 submarine pipeline' with 'fixed' and 'lifetime' + columns + co2_network_cost_factor : float, optional + Factor to scale the capital costs of the CO2 network, default 1.0 + + Returns + ------- + None + Modifies the network object in-place by adding CO2 pipeline links + + Notes + ----- + The function creates bidirectional CO2 pipeline links between nodes, with costs + depending on the underwater fraction of the pipeline. The network topology is + created using the create_network_topology helper function. + """ logger.info("Adding CO2 network.") co2_links = create_network_topology(n, "CO2 pipeline ") @@ -728,8 +829,7 @@ def add_co2_network(n, costs): * co2_links.length ) capital_cost = cost_onshore + cost_submarine - cost_factor = snakemake.config["sector"]["co2_network_cost_factor"] - capital_cost *= cost_factor + capital_cost *= co2_network_cost_factor n.add( "Link", @@ -1040,15 +1140,45 @@ def add_dac(n, costs): ) -def add_co2limit(n, options, nyears=1.0, limit=0.0): - logger.info(f"Adding CO2 budget limit as per unit of 1990 levels of {limit}") +def add_co2limit(n, options, co2_totals_file, countries, nyears, limit): + """ + Add a global CO2 emissions constraint to the network. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network container object + options : dict + Dictionary of options determining which sectors to consider for emissions + co2_totals_file : str + Path to CSV file containing historical CO2 emissions data in Mt + (megatonnes) per country and sector + countries : list + List of country codes to consider for the CO2 limit + nyears : float, optional + Number of years for the CO2 budget, by default 1.0 + limit : float, optional + CO2 limit as a fraction of 1990 levels, by default 0.0 - countries = snakemake.params.countries + Returns + ------- + None + The function modifies the network object in-place by adding a global + CO2 constraint. + + Notes + ----- + The function reads historical CO2 emissions data, calculates a total limit + based on the specified countries and sectors, and adds a global constraint + to the network. The limit is calculated as a fraction of historical emissions + multiplied by the number of years. + """ + logger.info(f"Adding CO2 budget limit as per unit of 1990 levels of {limit}") sectors = determine_emission_sectors(options) # convert Mt to tCO2 - co2_totals = 1e6 * pd.read_csv(snakemake.input.co2_totals_name, index_col=0) + co2_totals = 1e6 * pd.read_csv(co2_totals_file, index_col=0) co2_limit = co2_totals.loc[countries, sectors].sum().sum() @@ -1349,7 +1479,74 @@ def add_electricity_grid_connection(n, costs): ] -def add_storage_and_grids(n, costs): +def add_storage_and_grids( + n, + costs, + pop_layout, + h2_cavern_file, + cavern_types, + clustered_gas_network_file, + gas_input_nodes_file, + spatial, + options, +): + """ + Add storage and grid infrastructure to the network including hydrogen, gas, and battery systems. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network container object + costs : pd.DataFrame + Technology cost assumptions + pop_layout : pd.DataFrame + Population layout with index of locations/nodes + h2_cavern_file : str + Path to CSV file containing hydrogen cavern storage potentials + cavern_types : list + List of underground storage types to consider + clustered_gas_network_file : str, optional + Path to CSV file containing gas network data + gas_input_nodes_file : str, optional + Path to CSV file containing gas input nodes data + spatial : object, optional + Object containing spatial information about nodes and their locations + options : dict, optional + Dictionary of configuration options. Defaults to empty dict if not provided. + Key options include: + - hydrogen_fuel_cell : bool + - hydrogen_turbine : bool + - hydrogen_underground_storage : bool + - gas_network : bool + - H2_retrofit : bool + - H2_network : bool + - methanation : bool + - coal_cc : bool + - SMR_cc : bool + - SMR : bool + - marginal_cost_storage : float + - min_part_load_methanation : float + - cc_fraction : float + logger : logging.Logger, optional + Logger for output messages. If None, no logging is performed. + + Returns + ------- + None + The function modifies the network object in-place by adding various + storage and grid components. + + Notes + ----- + This function adds multiple types of storage and grid infrastructure: + - Hydrogen infrastructure (electrolysis, fuel cells, storage) + - Gas network infrastructure + - Battery storage systems + - Carbon capture and conversion facilities (if enabled in options) + """ + # Set defaults + options = options or {} + logger.info("Add hydrogen storage") nodes = pop_layout.index @@ -1406,8 +1603,7 @@ def add_storage_and_grids(n, costs): lifetime=costs.at["OCGT", "lifetime"], ) - cavern_types = snakemake.params.sector["hydrogen_underground_storage_locations"] - h2_caverns = pd.read_csv(snakemake.input.h2_cavern, index_col=0) + h2_caverns = pd.read_csv(h2_cavern_file, index_col=0) if ( not h2_caverns.empty @@ -1457,8 +1653,7 @@ def add_storage_and_grids(n, costs): ) if options["gas_network"] or options["H2_retrofit"]: - fn = snakemake.input.clustered_gas_network - gas_pipes = pd.read_csv(fn, index_col=0) + gas_pipes = pd.read_csv(clustered_gas_network_file, index_col=0) if options["gas_network"]: logger.info( @@ -1499,8 +1694,7 @@ def add_storage_and_grids(n, costs): # remove fossil generators where there is neither # production, LNG terminal, nor entry-point beyond system scope - fn = snakemake.input.gas_input_nodes_simplified - gas_input_nodes = pd.read_csv(fn, index_col=0) + gas_input_nodes = pd.read_csv(gas_input_nodes_file, index_col=0) unique = gas_input_nodes.index.unique() gas_i = n.generators.carrier == "gas" @@ -1939,36 +2133,86 @@ def add_ice_cars(n, p_set, ice_share, temperature): ) -def add_land_transport(n, costs): - logger.info("Add land transport") +def add_land_transport( + n, + costs, + transport_demand_file, + transport_data_file, + avail_profile_file, + dsm_profile_file, + temp_air_total_file, + options, + investment_year, + nodes, + logger=None, +): + """ + Add land transport demand and infrastructure to the network. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network container object + costs : pd.DataFrame + Cost assumptions for different technologies + transport_demand_file : str + Path to CSV file containing transport demand in driven km [100 km] + transport_data_file : str + Path to CSV file containing number of cars per region + avail_profile_file : str + Path to CSV file containing availability profiles + dsm_profile_file : str + Path to CSV file containing demand-side management profiles + temp_air_total_file : str + Path to netCDF file containing air temperature data + options : dict + Dictionary containing configuration options, specifically: + - land_transport_fuel_cell_share + - land_transport_electric_share + - land_transport_ice_share + investment_year : int + Year for which to get the transport shares + nodes : list-like + List of spatial nodes to consider + logger : logging.Logger, optional + Logger instance for output. If None, logging is skipped. + + Returns + ------- + None + Modifies the network object in-place by adding transport-related + components and their properties. + + Notes + ----- + The function adds different types of land transport (electric vehicles, + fuel cell vehicles, and internal combustion engines) to the network + based on specified shares and profiles. + """ + if logger: + logger.info("Add land transport") # read in transport demand in units driven km [100 km] - 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 - ) + transport = pd.read_csv(transport_demand_file, index_col=0, parse_dates=True) + number_cars = pd.read_csv(transport_data_file, index_col=0)["number cars"] + avail_profile = pd.read_csv(avail_profile_file, index_col=0, parse_dates=True) + dsm_profile = pd.read_csv(dsm_profile_file, 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}%") + share_key = f"land_transport_{engine}_share" + shares[engine] = get(options[share_key], investment_year) + if logger: + logger.info(f"{engine} share: {shares[engine] * 100}%") check_land_transport_shares(shares) - p_set = transport[spatial.nodes] + p_set = transport[nodes] # temperature for correction factor for heating/cooling - temperature = xr.open_dataarray(snakemake.input.temp_air_total).to_pandas() + temperature = xr.open_dataarray(temp_air_total_file).to_pandas() if shares["electric"] > 0: add_EVs( @@ -1988,11 +2232,39 @@ def add_land_transport(n, costs): add_ice_cars(n, p_set, shares["ice"], temperature) -def build_heat_demand(n): +def build_heat_demand( + n, hourly_heat_demand_file, pop_weighted_energy_totals, heating_efficiencies +): + """ + Build heat demand time series and adjust electricity load to account for electric heating. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network container object + hourly_heat_demand_file : str + Path to netCDF file containing hourly heat demand data + pop_weighted_energy_totals : pd.DataFrame + Population-weighted energy totals containing columns for total and + electricity consumption for different sectors and uses + heating_efficiencies : dict + Dictionary mapping sector and use combinations to their heating efficiencies + + Returns + ------- + pd.DataFrame + Heat demand time series with hierarchical columns for different sectors + and uses (residential/services, water/space) + + Notes + ----- + The function: + - Constructs heat demand profiles for different sectors and uses + - Adjusts the electricity load profiles by subtracting electric heating + - Modifies the network object in-place by updating n.loads_t.p_set + """ heat_demand_shape = ( - xr.open_dataset(snakemake.input.hourly_heat_demand_total) - .to_dataframe() - .unstack(level=1) + xr.open_dataset(hourly_heat_demand_file).to_dataframe().unstack(level=1) ) sectors = [sector.value for sector in HeatSector] @@ -2031,29 +2303,95 @@ def build_heat_demand(n): def add_heat( n: pypsa.Network, costs: pd.DataFrame, - cop: xr.DataArray, - direct_heat_source_utilisation_profile: xr.DataArray, + cop_profiles_file: str, + direct_heat_source_utilisation_profile_file: str, + hourly_heat_demand_total_file: str, + district_heat_share_file: str, + solar_thermal_total_file: str, + retro_cost_file: str, + floor_area_file: str, + heat_source_profile_files: dict[str, str], + params: dict, + pop_weighted_energy_totals: pd.DataFrame, + heating_efficiencies: pd.DataFrame, + pop_layout: pd.DataFrame, + spatial: object, + options: dict, + investment_year: int, ): """ - Add heat sector to the network. + Add heat sector to the network including heat demand, heat pumps, storage, and conversion technologies. Parameters ---------- - n (pypsa.Network): The PyPSA network object. - costs (pd.DataFrame): DataFrame containing cost information. - cop (xr.DataArray): DataArray containing coefficient of performance (COP) values. + n : pypsa.Network + The PyPSA network object + costs : pd.DataFrame + DataFrame containing cost information for different technologies + cop_profiles_file : str + Path to NetCDF file containing coefficient of performance (COP) values for heat pumps + direct_heat_source_utilisation_profile_file : str + Path to NetCDF file containing direct heat source utilisation profiles + hourly_heat_demand_total_file : str + Path to CSV file containing hourly heat demand data + district_heat_share_file : str + Path to CSV file containing district heating share information + solar_thermal_total_file : str + Path to NetCDF file containing solar thermal generation data + retro_cost_file : str + Path to CSV file containing retrofitting costs + floor_area_file : str + Path to CSV file containing floor area data + heat_source_profile_files : dict[str, str] + Dictionary mapping heat source names to their data file paths + params : dict + Dictionary containing parameters including: + - heat_pump_sources + - heat_utilisation_potentials + - direct_utilisation_heat_sources + pop_weighted_energy_totals : pd.DataFrame + Population-weighted energy totals by region + heating_efficiencies : pd.DataFrame + Heating system efficiencies + pop_layout : pd.DataFrame + Population layout data with columns for fraction and country + spatial : object + Object containing spatial data with attributes for different carriers (gas, co2, etc.) + options : dict + Dictionary containing configuration options for heat sector components + investment_year : int + Year for which to get the heat sector components and costs Returns ------- - None + None + Modifies the network object in-place by adding heat sector components + + Notes + ----- + The function adds various heat sector components to the network including: + - Heat demand for different sectors (residential, services) + - Heat pumps with different heat sources + - Thermal energy storage if enabled + - Gas boilers if enabled + - Solar thermal if enabled + - Combined heat and power (CHP) plants if enabled + - Building retrofitting options if enabled """ logger.info("Add heat sector") sectors = [sector.value for sector in HeatSector] - heat_demand = build_heat_demand(n) + heat_demand = build_heat_demand( + n, + hourly_heat_demand_total_file, + pop_weighted_energy_totals, + heating_efficiencies, + ) - district_heat_info = pd.read_csv(snakemake.input.district_heat_share, index_col=0) + cop = xr.open_dataarray(cop_profiles_file) + direct_heat_profile = xr.open_dataarray(direct_heat_source_utilisation_profile_file) + district_heat_info = pd.read_csv(district_heat_share_file, index_col=0) dist_fraction = district_heat_info["district fraction of node"] urban_fraction = district_heat_info["urban fraction"] @@ -2068,7 +2406,7 @@ def add_heat( if options["solar_thermal"]: solar_thermal = ( - xr.open_dataarray(snakemake.input.solar_thermal_total) + xr.open_dataarray(solar_thermal_total_file) .to_pandas() .reindex(index=n.snapshots) ) @@ -2147,9 +2485,7 @@ def add_heat( ) ## Add heat pumps - for heat_source in snakemake.params.heat_pump_sources[ - heat_system.system_type.value - ]: + for heat_source in params.heat_pump_sources[heat_system.system_type.value]: costs_name_heat_pump = heat_system.heat_pump_costs_name(heat_source) cop_heat_pump = ( cop.sel( @@ -2163,10 +2499,10 @@ def add_heat( else costs.at[costs_name_heat_pump, "efficiency"] ) - if heat_source in snakemake.params.heat_utilisation_potentials: + if heat_source in params.heat_utilisation_potentials: # get potential p_max_source = pd.read_csv( - snakemake.input[heat_source], + heat_source_profile_files[heat_source], index_col=0, ).squeeze()[nodes] @@ -2180,7 +2516,7 @@ def add_heat( carrier=heat_carrier, ) - if heat_source in snakemake.params.direct_utilisation_heat_sources: + if heat_source in params.direct_utilisation_heat_sources: capital_cost = ( costs.at[ heat_system.heat_source_costs_name(heat_source), "fixed" @@ -2223,10 +2559,10 @@ def add_heat( lifetime=costs.at[costs_name_heat_pump, "lifetime"], ) - if heat_source in snakemake.params.direct_utilisation_heat_sources: + if heat_source in params.direct_utilisation_heat_sources: # 1 if source temperature exceeds forward temperature, 0 otherwise: efficiency_direct_utilisation = ( - direct_heat_source_utilisation_profile.sel( + direct_heat_profile.sel( heat_source=heat_source, name=nodes, ) @@ -2459,13 +2795,13 @@ def add_heat( # demand 'dE' [per unit of original heat demand] for each country and # different retrofitting strengths [additional insulation thickness in m] retro_data = pd.read_csv( - snakemake.input.retro_cost, + retro_cost_file, index_col=[0, 1], skipinitialspace=True, header=[0, 1], ) # heated floor area [10^6 * m^2] per country - floor_area = pd.read_csv(snakemake.input.floor_area, index_col=[0, 1]) + floor_area_file = pd.read_csv(floor_area_file, index_col=[0, 1]) n.add("Carrier", "retrofitting") @@ -2504,7 +2840,7 @@ def add_heat( # get floor aread at node and region (urban/rural) in m^2 floor_area_node = ( - pop_layout.loc[node].fraction * floor_area.loc[ct, "value"] * 10**6 + pop_layout.loc[node].fraction * floor_area_file.loc[ct, "value"] * 10**6 ).loc[sec] * f # total heat demand at node [MWh] demand = n.loads_t.p_set[name] @@ -2591,10 +2927,68 @@ def add_methanol(n, costs): add_methanol_reforming_cc(n, costs) -def add_biomass(n, costs): +def add_biomass( + n, + costs, + options, + spatial, + cf_industry, + pop_layout, + biomass_potentials_file, + biomass_transport_costs_file=None, +): + """ + Add biomass-related components to the PyPSA network. + + This function adds various biomass-related components including biogas, + solid biomass, municipal solid waste, biomass transport, and different + biomass conversion technologies (CHP, boilers, BtL, BioSNG, etc.). + + Parameters + ---------- + n : pypsa.Network + The PyPSA network container object + costs : pd.DataFrame + DataFrame containing technology cost assumptions + options : dict + Dictionary of configuration options including keys like: + - gas_network : bool + - biomass_transport : bool + - biomass_spatial : bool + - municipal_solid_waste : bool + - biomass_to_liquid : bool + - etc. + spatial : object + Object containing spatial information about different carriers (gas, biomass, etc.) + cf_industry : dict + Dictionary containing industrial sector configuration + pop_layout : pd.DataFrame + DataFrame containing population layout information + biomass_potentials_file : str + Path to CSV file containing biomass potentials data + biomass_transport_costs_file : str, optional + Path to CSV file containing biomass transport costs data. + Required if biomass_transport or biomass_spatial options are True. + + Returns + ------- + None + The function modifies the network object in-place by adding + biomass-related components. + + Notes + ----- + The function adds various types of biomass-related components depending + on the options provided, including: + - Biogas and solid biomass generators + - Municipal solid waste if enabled + - Biomass transport infrastructure + - Biomass conversion technologies (CHP, boilers, BtL, BioSNG) + - Carbon capture options for different processes + """ logger.info("Add biomass") - biomass_potentials = pd.read_csv(snakemake.input.biomass_potentials, index_col=0) + biomass_potentials = pd.read_csv(biomass_potentials_file, index_col=0) # need to aggregate potentials if gas not nodally resolved if options["gas_network"]: @@ -2865,9 +3259,7 @@ def add_biomass(n, costs): if options["biomass_transport"]: # add biomass transport - transport_costs = pd.read_csv( - snakemake.input.biomass_transport_costs, index_col=0 - ) + transport_costs = pd.read_csv(biomass_transport_costs_file, index_col=0) transport_costs = transport_costs.squeeze() biomass_transport = create_network_topology( n, "biomass transport ", bidirectional=False @@ -2909,9 +3301,7 @@ def add_biomass(n, costs): elif options["biomass_spatial"]: # add artificial biomass generators at nodes which include transport costs - transport_costs = pd.read_csv( - snakemake.input.biomass_transport_costs, index_col=0 - ) + transport_costs = pd.read_csv(biomass_transport_costs_file, index_col=0) transport_costs = transport_costs.squeeze() bus_transport_costs = spatial.biomass.nodes.to_series().apply( lambda x: transport_costs[x[:2]] @@ -3231,7 +3621,69 @@ def add_biomass(n, costs): ) -def add_industry(n, costs): +def add_industry( + n, + costs, + industrial_demand_file, + shipping_demand_file, + pop_layout, + pop_weighted_energy_totals, + options, + spatial, + cf_industry, + investment_year, +): + """ + Add industry, shipping, aviation, and their corresponding carrier buses to the network. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network container object + costs : pd.DataFrame + Costs data including carbon capture, fuel costs, etc. + industrial_demand_file : str + Path to CSV file containing industrial demand data + shipping_demand_file : str + Path to CSV file containing shipping demand data + pop_layout : pd.DataFrame + Population layout data with index of nodes + pop_weighted_energy_totals : pd.DataFrame + Population-weighted energy totals including aviation and navigation data + options : dict + Dictionary of configuration options including: + - biomass_spatial + - biomass_transport + - gas_network + - methanol configuration + - regional_oil_demand + - shipping shares (hydrogen, methanol, oil) + - and others + spatial : object + Object containing spatial configuration for different carriers + (biomass, gas, oil, methanol, etc.) + cf_industry : dict + Industry-specific configuration parameters + investment_year : int + Year for which investment costs should be considered + HeatSystem : Enum + Enumeration defining different heat system types + + Returns + ------- + None + The function modifies the network object in-place by adding + industry-related components. + + Notes + ----- + This function adds multiple components to the network including: + - Industrial demand for various carriers + - Shipping and aviation infrastructure + - Carbon capture facilities + - Heat systems + - Process emission handling + """ logger.info("Add industrial demand") # add oil buses for shipping, aviation and naptha for industry add_carrier_buses(n, "oil") @@ -3243,9 +3695,7 @@ def add_industry(n, costs): nyears = nhours / 8760 # 1e6 to convert TWh to MWh - industrial_demand = ( - pd.read_csv(snakemake.input.industrial_demand, index_col=0) * 1e6 - ) * nyears + industrial_demand = pd.read_csv(industrial_demand_file, index_col=0) * 1e6 * nyears n.add( "Bus", @@ -3442,8 +3892,7 @@ def add_industry(n, costs): nodes, ["total domestic navigation"] ].squeeze() international_navigation = ( - pd.read_csv(snakemake.input.shipping_demand, index_col=0).squeeze(axis=1) - * nyears + pd.read_csv(shipping_demand_file, index_col=0).squeeze(axis=1) * nyears ) all_navigation = domestic_navigation + international_navigation p_set = all_navigation * 1e6 / nhours @@ -4393,13 +4842,62 @@ def lossy_bidirectional_links(n, carrier, efficiencies={}): ) -def add_enhanced_geothermal(n, egs_potentials, egs_overlap, costs): +def add_enhanced_geothermal( + n, + costs, + costs_config, + egs_potentials, + egs_overlap, + egs_config, + egs_capacity_factors=None, +): """ - Adds EGS potential to model. + Add Enhanced Geothermal System (EGS) potential to the network model. - Built in scripts/build_egs_potentials.py - """ + Parameters + ---------- + n : pypsa.Network + The PyPSA network container object. + egs_potentials : str + Path to CSV file containing EGS potential data. + egs_overlap : str + Path to CSV file defining overlap between gridded geothermal potential + estimation and bus regions. + costs : pd.DataFrame + Technology cost assumptions including fields for lifetime, FOM, investment, + and efficiency parameters. + egs_config : dict + Configuration for enhanced geothermal systems with keys: + - var_cf : bool + Whether to use time-varying capacity factors + - flexible : bool + Whether to add flexible operation using geothermal reservoir + - max_hours : float + Maximum hours of storage if flexible + - max_boost : float + Maximum power boost factor if flexible + costs_config : dict + General cost configuration containing: + - fill_values : dict + With key 'discount rate' for financial calculations + egs_capacity_factors : str, optional + Path to CSV file with time-varying capacity factors. + Required if egs_config['var_cf'] is True. + Returns + ------- + None + Modifies the network object in-place by adding EGS components. + + Notes + ----- + Implements EGS with 2020 CAPEX from Aghahosseini et al 2021. + The function adds various components to the network: + - Geothermal heat generators and buses + - Organic Rankine Cycle links for electricity generation + - District heating links if urban central heat exists + - Optional storage units for flexible operation + """ if len(spatial.geothermal_heat.nodes) > 1: logger.warning( "'add_enhanced_geothermal' not implemented for multiple geothermal nodes." @@ -4417,9 +4915,6 @@ def add_enhanced_geothermal(n, egs_potentials, egs_overlap, costs): ) logger.info("[EGS] district heat distribution part -> 'geothermal district heat'") - egs_config = snakemake.params["sector"]["enhanced_geothermal"] - costs_config = snakemake.config["costs"] - # matrix defining the overlap between gridded geothermal potential estimation, and bus regions overlap = pd.read_csv(egs_overlap, index_col=0) overlap.columns = overlap.columns.astype(int) @@ -4478,9 +4973,7 @@ def add_enhanced_geothermal(n, egs_potentials, egs_overlap, costs): ) if egs_config["var_cf"]: - efficiency = pd.read_csv( - snakemake.input.egs_capacity_factors, parse_dates=True, index_col=0 - ) + efficiency = pd.read_csv(egs_capacity_factors, parse_dates=True, index_col=0) logger.info("Adding Enhanced Geothermal with time-varying capacity factors.") else: efficiency = 1.0 @@ -4639,12 +5132,18 @@ def add_enhanced_geothermal(n, egs_potentials, egs_overlap, costs): ) pop_weighted_energy_totals.update(pop_weighted_heat_totals) + carriers_to_keep = snakemake.params.pypsa_eur + profiles = { + key: snakemake.input[key] + for key in snakemake.input.keys() + if key.startswith("profile") + } landfall_lengths = { tech: settings["landfall_length"] for tech, settings in snakemake.params.renewable.items() if "landfall_length" in settings.keys() } - patch_electricity_network(n, costs, landfall_lengths) + patch_electricity_network(n, costs, carriers_to_keep, profiles, landfall_lengths) fn = snakemake.input.heating_efficiencies year = int(snakemake.params["energy_totals_year"]) @@ -4661,27 +5160,78 @@ def add_enhanced_geothermal(n, egs_potentials, egs_overlap, costs): add_eu_bus(n) - add_co2_tracking(n, costs, options) + add_co2_tracking( + n, + costs, + options, + sequestration_potential_file=snakemake.input.sequestration_potential, + ) add_generation(n, costs) - add_storage_and_grids(n, costs) + add_storage_and_grids( + n=n, + costs=costs, + pop_layout=pop_layout, + h2_cavern_file=snakemake.input.h2_cavern, + cavern_types=snakemake.params.sector["hydrogen_underground_storage_locations"], + clustered_gas_network_file=snakemake.input.clustered_gas_network, + gas_input_nodes_file=snakemake.input.gas_input_nodes_simplified, + spatial=spatial, + options=options, + ) if options["transport"]: - add_land_transport(n, costs) + add_land_transport( + n=n, + costs=costs, + transport_demand_file=snakemake.input.transport_demand, + transport_data_file=snakemake.input.transport_data, + avail_profile_file=snakemake.input.avail_profile, + dsm_profile_file=snakemake.input.dsm_profile, + temp_air_total_file=snakemake.input.temp_air_total, + options=options, + investment_year=investment_year, + nodes=spatial.nodes, + logger=logger, + ) if options["heating"]: add_heat( n=n, costs=costs, - cop=xr.open_dataarray(snakemake.input.cop_profiles), - direct_heat_source_utilisation_profile=xr.open_dataarray( - snakemake.input.direct_heat_source_utilisation_profiles - ), + cop_profiles_file=snakemake.input.cop_profiles, + direct_heat_source_utilisation_profile_file=snakemake.input.direct_heat_source_utilisation_profiles, + hourly_heat_demand_total_file=snakemake.input.hourly_heat_demand_total, + district_heat_share_file=snakemake.input.district_heat_share, + solar_thermal_total_file=snakemake.input.solar_thermal_total, + retro_cost_file=snakemake.input.retro_cost, + floor_area_file=snakemake.input.floor_area, + heat_source_profile_files={ + source: snakemake.input[source] + for source in snakemake.params.heat_utilisation_potentials + if source in snakemake.input.keys() + }, + params=snakemake.params, + pop_weighted_energy_totals=pop_weighted_energy_totals, + heating_efficiencies=heating_efficiencies, + pop_layout=pop_layout, + spatial=spatial, + options=options, + investment_year=investment_year, ) if options["biomass"]: - add_biomass(n, costs) + add_biomass( + n=n, + costs=costs, + options=options, + spatial=spatial, + cf_industry=cf_industry, + pop_layout=pop_layout, + biomass_potentials_file=snakemake.input.biomass_potentials, + biomass_transport_costs_file=snakemake.input.biomass_transport_costs, + ) if options["ammonia"]: add_ammonia(n, costs) @@ -4690,7 +5240,18 @@ def add_enhanced_geothermal(n, egs_potentials, egs_overlap, costs): add_methanol(n, costs) if options["industry"]: - add_industry(n, costs) + add_industry( + n=n, + costs=costs, + industrial_demand_file=snakemake.input.industrial_demand, + shipping_demand_file=snakemake.input.shipping_demand, + pop_layout=pop_layout, + pop_weighted_energy_totals=pop_weighted_energy_totals, + options=options, + spatial=spatial, + cf_industry=cf_industry, + investment_year=investment_year, + ) if options["heating"]: add_waste_heat(n) @@ -4708,7 +5269,13 @@ def add_enhanced_geothermal(n, egs_potentials, egs_overlap, costs): remove_h2_network(n) if options["co2_network"]: - add_co2_network(n, costs) + add_co2_network( + n, + costs, + co2_network_cost_factor=snakemake.config["sector"][ + "co2_network_cost_factor" + ], + ) if options["allam_cycle_gas"]: add_allam_gas(n, costs) @@ -4730,12 +5297,21 @@ def add_enhanced_geothermal(n, egs_potentials, egs_overlap, costs): emissions_scope, input_co2, options, + snakemake.params.countries, + snakemake.params.planning_horizons, ) co2_cap = pd.read_csv(fn, index_col=0).squeeze() limit = co2_cap.loc[investment_year] else: limit = get(co2_budget, investment_year) - add_co2limit(n, options, nyears, limit) + add_co2limit( + n, + options, + snakemake.input.co2_totals_name, + snakemake.params.countries, + nyears, + limit, + ) maxext = snakemake.params["lines"]["max_extension"] if maxext is not None: @@ -4747,7 +5323,13 @@ def add_enhanced_geothermal(n, egs_potentials, egs_overlap, costs): if options["enhanced_geothermal"].get("enable", False): logger.info("Adding Enhanced Geothermal Systems (EGS).") add_enhanced_geothermal( - n, snakemake.input["egs_potentials"], snakemake.input["egs_overlap"], costs + n, + costs=costs, + costs_config=snakemake.config["costs"], + egs_potentials=snakemake.input["egs_potentials"], + egs_overlap=snakemake.input["egs_overlap"], + egs_config=snakemake.params["sector"]["enhanced_geothermal"], + egs_capacity_factors="path/to/capacity_factors.csv", ) if options["gas_distribution_grid"]: diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 92c4d9e9d..b3fe65006 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -31,6 +31,7 @@ import os import re import sys +from typing import Any import numpy as np import pandas as pd @@ -55,9 +56,19 @@ class ObjectiveValueError(Exception): pass -def add_land_use_constraint_perfect(n): +def add_land_use_constraint_perfect(n: pypsa.Network) -> None: """ Add global constraints for tech capacity limit. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network instance + + Returns + ------- + pypsa.Network + Network with added land use constraints """ logger.info("Add land-use constraint for perfect foresight") @@ -115,10 +126,23 @@ def check_p_min_p_max(p_nom_max): bus = df_carrier.bus n.buses.loc[bus, name] = df_carrier.p_nom_max.values - return n +def add_land_use_constraint(n: pypsa.Network, planning_horizon: str) -> None: + """ + Add land use constraints for renewable energy potential. -def add_land_use_constraint(n): + Parameters + ---------- + n : pypsa.Network + The PyPSA network instance + planning_horizon : str + The planning horizon year as string + + Returns + ------- + pypsa.Network + Modified PyPSA network with constraints added + """ # warning: this will miss existing offwind which is not classed AC-DC and has carrier 'offwind' for carrier in [ @@ -136,7 +160,7 @@ def add_land_use_constraint(n): .groupby(n.generators.bus.map(n.buses.location)) .sum() ) - existing.index += " " + carrier + "-" + snakemake.wildcards.planning_horizons + existing.index += f" {carrier}-{planning_horizon}" n.generators.loc[existing.index, "p_nom_max"] -= existing # check if existing capacities are larger than technical potential @@ -155,7 +179,7 @@ def add_land_use_constraint(n): n.generators["p_nom_max"] = n.generators["p_nom_max"].clip(lower=0) -def add_solar_potential_constraints(n, config): +def add_solar_potential_constraints(n: pypsa.Network, config: dict) -> None: """ Add constraint to make sure the sum capacity of all solar technologies (fixed, tracking, ets. ) is below the region potential. @@ -213,7 +237,11 @@ def add_solar_potential_constraints(n, config): n.model.add_constraints(lhs <= rhs, name="solar_potential") -def add_co2_sequestration_limit(n, limit_dict): +def add_co2_sequestration_limit( + n: pypsa.Network, + limit_dict: dict[str, float], + planning_horizons: str, +) -> None: """ Add a global constraint on the amount of Mt CO2 that can be sequestered. """ @@ -228,7 +256,7 @@ def add_co2_sequestration_limit(n, limit_dict): ) names = limit.index else: - limit = get(limit_dict, int(snakemake.wildcards.planning_horizons)) + limit = get(limit_dict, int(planning_horizons)) periods = [np.nan] names = pd.Index(["co2_sequestration_limit"]) @@ -243,7 +271,7 @@ def add_co2_sequestration_limit(n, limit_dict): ) -def add_carbon_constraint(n, snapshots): +def add_carbon_constraint(n: pypsa.Network, snapshots: pd.DatetimeIndex) -> None: glcs = n.global_constraints.query('type == "co2_atmosphere"') if glcs.empty: return @@ -269,7 +297,7 @@ def add_carbon_constraint(n, snapshots): n.model.add_constraints(lhs <= rhs, name=f"GlobalConstraint-{name}") -def add_carbon_budget_constraint(n, snapshots): +def add_carbon_budget_constraint(n: pypsa.Network, snapshots: pd.DatetimeIndex) -> None: glcs = n.global_constraints.query('type == "Co2Budget"') if glcs.empty: return @@ -296,12 +324,11 @@ def add_carbon_budget_constraint(n, snapshots): n.model.add_constraints(lhs <= rhs, name=f"GlobalConstraint-{name}") -def add_max_growth(n): +def add_max_growth(n: pypsa.Network, opts: dict) -> None: """ Add maximum growth rates for different carriers. """ - opts = snakemake.params["sector"]["limit_max_growth"] # take maximum yearly difference between investment periods since historic growth is per year factor = n.investment_period_weightings.years.max() * opts["factor"] for carrier in opts["max_growth"].keys(): @@ -318,10 +345,10 @@ def add_max_growth(n): ) n.carriers.loc[carrier, "max_relative_growth"] = max_r_per_period - return n - -def add_retrofit_gas_boiler_constraint(n, snapshots): +def add_retrofit_gas_boiler_constraint( + n: pypsa.Network, snapshots: pd.DatetimeIndex +) -> None: """ Allow retrofitting of existing gas boilers to H2 boilers. """ @@ -366,13 +393,34 @@ def add_retrofit_gas_boiler_constraint(n, snapshots): def prepare_network( - n, - solve_opts=None, - config=None, - foresight=None, - planning_horizons=None, - co2_sequestration_potential=None, -): + n: pypsa.Network, + solve_opts: dict, + foresight: str, + planning_horizons: str, + co2_sequestration_potential: dict[str, float], + limit_max_growth: dict[str, Any] | None = None, +) -> None: + """ + Prepare network with various constraints and modifications. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network instance + solve_opts : Dict + Dictionary of solving options containing clip_p_max_pu, load_shedding etc. + foresight : str + Planning foresight type ('myopic' or 'perfect') + planning_horizons : str + List of planning horizon years + co2_sequestration_potential : Dict[str, float] + CO2 sequestration potential constraints by year + + Returns + ------- + pypsa.Network + Modified PyPSA network with added constraints + """ if "clip_p_max_pu" in solve_opts: for df in ( n.generators_t.p_max_pu, @@ -440,21 +488,21 @@ def prepare_network( n.snapshot_weightings[:] = 8760.0 / nhours if foresight == "myopic": - add_land_use_constraint(n) + add_land_use_constraint(n, planning_horizons) if foresight == "perfect": - n = add_land_use_constraint_perfect(n) - if snakemake.params["sector"]["limit_max_growth"]["enable"]: - n = add_max_growth(n) + add_land_use_constraint_perfect(n) + if limit_max_growth is not None and limit_max_growth["enable"]: + add_max_growth(n, limit_max_growth) if n.stores.carrier.eq("co2 sequestered").any(): limit_dict = co2_sequestration_potential - add_co2_sequestration_limit(n, limit_dict=limit_dict) - - return n + add_co2_sequestration_limit( + n, limit_dict=limit_dict, planning_horizons=planning_horizons + ) -def add_CCL_constraints(n, config): +def add_CCL_constraints(n: pypsa.Network, config: dict, planning_horizons: str) -> None: """ Add CCL (country & carrier limit) constraint to the network. @@ -466,6 +514,7 @@ def add_CCL_constraints(n, config): ---------- n : pypsa.Network config : dict + planning_horizons : str Example ------- @@ -476,7 +525,7 @@ def add_CCL_constraints(n, config): """ agg_p_nom_minmax = pd.read_csv( config["solving"]["agg_p_nom_limits"]["file"], index_col=[0, 1], header=[0, 1] - )[snakemake.wildcards.planning_horizons] + )[planning_horizons] logger.info("Adding generation capacity constraints per carrier and country") p_nom = n.model["Generator-p_nom"] @@ -496,8 +545,7 @@ def add_CCL_constraints(n, config): index="Generator-cst" ) gens_cst = gens_cst[ - (gens_cst["build_year"] + gens_cst["lifetime"]) - >= int(snakemake.wildcards.planning_horizons) + (gens_cst["build_year"] + gens_cst["lifetime"]) >= int(planning_horizons) ] if config["solving"]["agg_p_nom_limits"]["agg_offwind"]: gens_cst = gens_cst.replace(rename_offwind) @@ -608,30 +656,16 @@ def add_EQ_constraints(n, o, scaling=1e-1): n.model.add_constraints(lhs >= rhs, name="equity_min") -def add_BAU_constraints(n, config): +def add_BAU_constraints(n: pypsa.Network, config: dict) -> None: """ - Add a per-carrier minimal overall capacity. - - BAU_mincapacities and opts must be adjusted in the config.yaml. + Add business-as-usual (BAU) constraints for minimum capacities. Parameters ---------- n : pypsa.Network + PyPSA network instance config : dict - - Example - ------- - scenario: - opts: [Co2L-BAU-24h] - electricity: - BAU_mincapacities: - solar: 0 - onwind: 0 - OCGT: 100000 - offwind-ac: 0 - offwind-dc: 0 - Which sets minimum expansion across all nodes e.g. in Europe to 100GW. - OCGT bus 1 + OCGT bus 2 + ... > 100000 + Configuration dictionary containing BAU minimum capacities """ mincaps = pd.Series(config["electricity"]["BAU_mincapacities"]) p_nom = n.model["Generator-p_nom"] @@ -725,18 +759,18 @@ def add_operational_reserve_margin(n, sns, config): p_nom_vres * (-EPSILON_VRES * xr.DataArray(capacity_factor)) ).sum("Generator") - # Total demand per t - demand = get_as_dense(n, "Load", "p_set").sum(axis=1) + # Total demand per t + demand = get_as_dense(n, "Load", "p_set").sum(axis=1) - # VRES potential of non extendable generators - capacity_factor = n.generators_t.p_max_pu[vres_i.difference(ext_i)] - renewable_capacity = n.generators.p_nom[vres_i.difference(ext_i)] - potential = (capacity_factor * renewable_capacity).sum(axis=1) + # VRES potential of non extendable generators + capacity_factor = n.generators_t.p_max_pu[vres_i.difference(ext_i)] + renewable_capacity = n.generators.p_nom[vres_i.difference(ext_i)] + potential = (capacity_factor * renewable_capacity).sum(axis=1) - # Right-hand-side - rhs = EPSILON_LOAD * demand + EPSILON_VRES * potential + CONTINGENCY + # Right-hand-side + rhs = EPSILON_LOAD * demand + EPSILON_VRES * potential + CONTINGENCY - n.model.add_constraints(lhs >= rhs, name="reserve_margin") + n.model.add_constraints(lhs >= rhs, name="reserve_margin") # additional constraint that capacity is not exceeded gen_i = n.generators.index @@ -918,7 +952,20 @@ def add_co2_atmosphere_constraint(n, snapshots): n.model.add_constraints(lhs <= rhs, name=f"GlobalConstraint-{name}") -def extra_functionality(n, snapshots): +def extra_functionality( + n: pypsa.Network, + snapshots: pd.DatetimeIndex, +) -> None: + """ + Add custom constraints and functionality. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network instance with config and params attributes + snapshots : pd.DatetimeIndex + Simulation timesteps + """ """ Collects supplementary constraints which will be passed to ``pypsa.optimization.optimize``. @@ -934,7 +981,7 @@ def extra_functionality(n, snapshots): if constraints["SAFE"] and n.generators.p_nom_extendable.any(): add_SAFE_constraints(n, config) if constraints["CCL"] and n.generators.p_nom_extendable.any(): - add_CCL_constraints(n, config) + add_CCL_constraints(n, config, n.params.planning_horizons) reserve = config["electricity"].get("operational_reserve", {}) if reserve.get("activate"): @@ -970,10 +1017,25 @@ def extra_functionality(n, snapshots): module_name = os.path.splitext(os.path.basename(source_path))[0] module = importlib.import_module(module_name) custom_extra_functionality = getattr(module, module_name) - custom_extra_functionality(n, snapshots, snakemake) + custom_extra_functionality(n, snapshots, snakemake) # pylint: disable=E0601 -def check_objective_value(n, solving): +def check_objective_value(n: pypsa.Network, solving: dict) -> None: + """ + Check if objective value matches expected value within tolerance. + + Parameters + ---------- + n : pypsa.Network + Network with solved objective + solving : Dict + Dictionary containing objective checking parameters + + Raises + ------ + ObjectiveValueError + If objective value differs from expected value beyond tolerance + """ check_objective = solving["check_objective"] if check_objective["enable"]: atol = check_objective["atol"] @@ -986,7 +1048,48 @@ def check_objective_value(n, solving): ) -def solve_network(n, config, params, solving, **kwargs): +def solve_network( + n: pypsa.Network, + config: dict, + params: dict, + solving: dict, + rule_name: str | None = None, + **kwargs, +) -> None: + """ + Solve network optimization problem. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network instance + config : Dict + Configuration dictionary containing solver settings + params : Dict + Dictionary of solving parameters + solving : Dict + Dictionary of solving options and configuration + rule_name : str, optional + Name of the snakemake rule being executed + **kwargs + Additional keyword arguments passed to the solver + + Returns + ------- + n : pypsa.Network + Solved network instance + status : str + Solution status + condition : str + Termination condition + + Raises + ------ + RuntimeError + If solving status is infeasible + ObjectiveValueError + If objective value differs from expected value + """ set_of_options = solving["solver"]["options"] cf_solving = solving["options"] @@ -1016,7 +1119,7 @@ def solve_network(n, config, params, solving, **kwargs): n.config = config n.params = params - if rolling_horizon and snakemake.rule == "solve_operations_network": + if rolling_horizon and rule_name == "solve_operations_network": kwargs["horizon"] = cf_solving.get("horizon", 365) kwargs["overlap"] = cf_solving.get("overlap", 0) n.optimize.optimize_with_rolling_horizon(**kwargs) @@ -1047,8 +1150,6 @@ def solve_network(n, config, params, solving, **kwargs): n.model.print_infeasibilities() raise RuntimeError("Solving status 'infeasible'") - return n - # %% if __name__ == "__main__": @@ -1073,14 +1174,15 @@ def solve_network(n, config, params, solving, **kwargs): np.random.seed(solve_opts.get("seed", 123)) n = pypsa.Network(snakemake.input.network) + planning_horizons = snakemake.wildcards.get("planning_horizons", None) - n = prepare_network( + prepare_network( n, - solve_opts, - config=snakemake.config, + solve_opts=snakemake.params.solving["options"], foresight=snakemake.params.foresight, - planning_horizons=snakemake.params.planning_horizons, + planning_horizons=planning_horizons, co2_sequestration_potential=snakemake.params["co2_sequestration_potential"], + limit_max_growth=snakemake.params.get("sector", {}).get("limit_max_growth"), ) logging_frequency = snakemake.config.get("solving", {}).get( @@ -1089,12 +1191,12 @@ def solve_network(n, config, params, solving, **kwargs): with memory_logger( filename=getattr(snakemake.log, "memory", None), interval=logging_frequency ) as mem: - n = solve_network( + solve_network( n, config=snakemake.config, params=snakemake.params, solving=snakemake.params.solving, - log_fn=snakemake.log.solver, + rule_name=snakemake.rule, ) logger.info(f"Maximum memory usage: {mem.mem_usage}")