diff --git a/Snakefile b/Snakefile index 1e53dc76..981c340d 100644 --- a/Snakefile +++ b/Snakefile @@ -13,6 +13,7 @@ configfile: "config.yaml" wildcard_constraints: weather_year="[0-9]{4}|", + capacity_year="[0-9]{4}|", lv="[a-z0-9\.]+", simpl="[a-zA-Z0-9]*", clusters="[0-9]+m?", @@ -680,3 +681,39 @@ if config["foresight"] == "myopic": resources: mem_mb=config['solving']['mem'] benchmark: RDIR + "/benchmarks/solve_network/elec{weather_year}_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}" script: "scripts/solve_network.py" + + +rule solve_operations_network: + input: + pre=RDIR + "/prenetworks/elec{weather_year}_s_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}.nc", + post=RDIR + "/postnetworks/elec{capacity_year}_s_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}.nc", + output: "elec{capacity_year}_s_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}_{weather_year}.nc" + shadow: "shallow" + log: + solver=RDIR + "/logs/elec{capacity_year}_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}_{weather_year}_solver.log", + python=RDIR + "/logs/elec{capacity_year}_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}_{weather_year}_python.log", + threads: 4 + resources: mem_mb=10000 + benchmark: RDIR + "/benchmarks/solve_operations_network/elec{capacity_year}_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}_{weather_year}" + script: "scripts/solve_operations_network.py" + + +def solved_previous_year(wildcards): + previous_year = int(wildcards.weather_year) - 1 + return RDIR + "/postnetworks/elec" + previous_year + "_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}.nc" + + +rule solve_operations_network_myopic: + input: + pre=RDIR + "/prenetworks/elec{weather_year}_s_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}.nc", + post=RDIR + "/postnetworks/elec{capacity_year}_s_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}.nc", + previous=solved_previous_year + output: "elec{capacity_year}_s_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}_{weather_year}_myopic.nc" + shadow: "shallow" + log: + solver=RDIR + "/logs/elec{capacity_year}_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}_{weather_year}_solver.log", + python=RDIR + "/logs/elec{capacity_year}_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}_{weather_year}_python.log", + threads: 4 + resources: mem_mb=10000 + benchmark: RDIR + "/benchmarks/solve_operations_network_myopic/elec{capacity_year}_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}_{weather_year}" + script: "scripts/solve_operations_network_myopic.py" diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 19ef5f52..4dd4d969 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -71,7 +71,7 @@ def prepare_network(n, solve_opts=None): df.where(df>solve_opts['clip_p_max_pu'], other=0., inplace=True) if solve_opts.get('load_shedding'): - n.add("Carrier", "Load") + n.add("Carrier", "load") n.madd("Generator", n.buses.index, " load", bus=n.buses.index, carrier='load', @@ -246,7 +246,7 @@ def extra_functionality(n, snapshots): add_co2_sequestration_limit(n, snapshots) -def solve_network(n, config, opts='', **kwargs): +def solve_network(n, config, opts='', snapshots=None, **kwargs): solver_options = config['solving']['solver'].copy() solver_name = solver_options.pop('name') cf_solving = config['solving']['options'] @@ -259,12 +259,15 @@ def solve_network(n, config, opts='', **kwargs): n.config = config n.opts = opts + if snapshots is None: + snapshots = n.snapshots + if cf_solving.get('skip_iterations', False): - network_lopf(n, solver_name=solver_name, solver_options=solver_options, + network_lopf(n, snapshots, solver_name=solver_name, solver_options=solver_options, extra_functionality=extra_functionality, keep_shadowprices=keep_shadowprices, **kwargs) else: - ilopf(n, solver_name=solver_name, solver_options=solver_options, + ilopf(n, snapshots, solver_name=solver_name, solver_options=solver_options, track_iterations=track_iterations, min_iterations=min_iterations, max_iterations=max_iterations, @@ -295,7 +298,7 @@ if __name__ == "__main__": if tmpdir is not None: from pathlib import Path Path(tmpdir).mkdir(parents=True, exist_ok=True) - opts = snakemake.wildcards.opts.split('-') + opts = snakemake.wildcards.sector_opts.split('-') solve_opts = snakemake.config['solving']['options'] fn = getattr(snakemake.log, 'memory', None) diff --git a/scripts/solve_operations_network.py b/scripts/solve_operations_network.py new file mode 100644 index 00000000..2021c7c9 --- /dev/null +++ b/scripts/solve_operations_network.py @@ -0,0 +1,126 @@ +"""Solve operations network.""" + + +import pypsa +import numpy as np + +from solve_network import solve_network, prepare_network +from helper import override_component_attrs + +import logging +logger = logging.getLogger(__name__) +pypsa.pf.logger.setLevel(logging.WARNING) + + +def set_parameters_from_optimized(n, n_optim): + lines_typed_i = n.lines.index[n.lines.type != ''] + n.lines.loc[lines_typed_i, 'num_parallel'] = \ + n_optim.lines['num_parallel'].reindex(lines_typed_i, fill_value=0.) + n.lines.loc[lines_typed_i, 's_nom'] = ( + np.sqrt(3) * n.lines['type'].map(n.line_types.i_nom) * + n.lines.bus0.map(n.buses.v_nom) * n.lines.num_parallel) + + lines_untyped_i = n.lines.index[n.lines.type == ''] + for attr in ('s_nom', 'r', 'x'): + n.lines.loc[lines_untyped_i, attr] = \ + n_optim.lines[attr].reindex(lines_untyped_i, fill_value=0.) + n.lines['s_nom_extendable'] = False + + links_dc_i = n.links.index[n.links.p_nom_extendable] + n.links.loc[links_dc_i, 'p_nom'] = \ + n_optim.links['p_nom_opt'].reindex(links_dc_i, fill_value=0.) + n.links.loc[links_dc_i, 'p_nom_extendable'] = False + + gen_extend_i = n.generators.index[n.generators.p_nom_extendable] + n.generators.loc[gen_extend_i, 'p_nom'] = \ + n_optim.generators['p_nom_opt'].reindex(gen_extend_i, fill_value=0.) + n.generators.loc[gen_extend_i, 'p_nom_extendable'] = False + + stor_units_extend_i = n.storage_units.index[n.storage_units.p_nom_extendable] + n.storage_units.loc[stor_units_extend_i, 'p_nom'] = \ + n_optim.storage_units['p_nom_opt'].reindex(stor_units_extend_i, fill_value=0.) + n.storage_units.loc[stor_units_extend_i, 'p_nom_extendable'] = False + + stor_extend_i = n.stores.index[n.stores.e_nom_extendable] + n.stores.loc[stor_extend_i, 'e_nom'] = \ + n_optim.stores['e_nom_opt'].reindex(stor_extend_i, fill_value=0.) + n.stores.loc[stor_extend_i, 'e_nom_extendable'] = False + + return n + + +def remove_unused_components(n, threshold=50): + logger.info("Remove assets that are barely used to speed things up.") + + for c in n.iterate_components({"Store", "Link", "Generator"}): + attr = "e_nom" if c.name == "Store" else "p_nom" + to_remove = c.df.loc[c.df[attr] < threshold].index + logger.info(f"Removing barely used {c.name}s:\n{to_remove}") + n.mremove(c.name, to_remove) + + return n + + +def add_load_shedding(n, voll=1e4): + logger.info("Add load shedding to all buses.") + + if "load" in n.generators.carrier.unique(): + to_remove = n.generators.query("carrier == 'load'").index + logger.info(f"Removing pre-existing load shedding:\n{to_remove}") + n.mremove("Generator", to_remove) + + n.madd("Generator", n.buses.index, + suffix=" load", + bus=n.buses.index, + carrier='load', + marginal_cost=voll, + p_nom=1e6 + ) + + return n + + +if __name__ == "__main__": + if 'snakemake' not in globals(): + from helper import mock_snakemake + snakemake = mock_snakemake( + 'solve_operations_network', + capacity_year=1952, + simpl='', + opts='', + clusters=37, + lv=2.0, + sector_opts='Co2L0-25H-T-H-B-I-A', + planning_horizons=2030, + weather_year=2013 + ) + + logging.basicConfig(filename=snakemake.log.python, + level=snakemake.config['logging_level']) + + tmpdir = snakemake.config['solving'].get('tmpdir') + if tmpdir is not None: + from pathlib import Path + Path(tmpdir).mkdir(parents=True, exist_ok=True) + + overrides = override_component_attrs(snakemake.input.overrides) + + n = pypsa.Network(snakemake.input.pre, override_component_attrs=overrides) + n_post = pypsa.Network(snakemake.input.post, override_component_attrs=overrides) + n = set_parameters_from_optimized(n, n_post) + del n_post + + n = remove_unused_components(n) + n = add_load_shedding(n) + + opts = snakemake.wildcards.sector_opts.split('-') + solve_opts = snakemake.config['solving']['options'] + solve_opts['skip_iterations'] = True + + n = prepare_network(n, solve_opts) + + n = solve_network(n, config=snakemake.config, opts=opts, + solver_dir=tmpdir, + solver_logfile=snakemake.log.solver) + + n.export_to_netcdf(snakemake.output[0]) diff --git a/scripts/solve_operations_network_myopic.py b/scripts/solve_operations_network_myopic.py new file mode 100644 index 00000000..a0e30b0a --- /dev/null +++ b/scripts/solve_operations_network_myopic.py @@ -0,0 +1,143 @@ +"""Solve myopic operations network.""" + + +import pypsa +import pandas as pd + +from solve_network import solve_network, prepare_network +from solve_operations_network import set_parameters_from_optimized, remove_unused_components, add_load_shedding +from helper import override_component_attrs + +import logging +logger = logging.getLogger(__name__) +pypsa.pf.logger.setLevel(logging.WARNING) + + +def prepare_myopic(n, config, store_soc, storage_unit_soc): + + n.stores.e_cyclic = False + n.storage_units.cyclic_state_of_charge = False + + biomass_stores = n.stores.carrier.str.isin(["solid biomass", "biogas"]) + biomass_potential = n.stores.loc[biomass_stores, "e_initial"] + + n.stores.e_initial = store_soc + n.storage_units.state_of_charge_initial = storage_unit_soc + + n.remove("GlobalConstraint", "CO2Limit") + n.stores.at["co2 atmosphere", "marginal_cost"] = -config["co2_price"] + + # handle co2 sequestration + assert sum(n.stores.carriers == "co2 stored") == 1, "Myopic operation not implemented for spatially resolved CO2 sequestration." + n.stores.at["co2 stored", 'e_nom'] = config['co2_sequestration_limit'] * 1e6 # t/a + + # reset co2 emissions + n.stores.loc[n.stores.carrier == 'co2 stored', "e_initial"] = 0. + n.stores.at["co2 atmosphere", "e_initial"] = 0. + + # replenish fossil gas and oil with 1000 TWh each + fossil_stores = n.stores.carrier.str.isin(["gas", "oil"]) + n.stores.loc[fossil_stores, 'e_initial'] = 1e9 + n.stores.loc[fossil_stores, 'e_nom'] = 10e9 + + # replenish annual solid biomass and biogas potentials + n.stores.loc[biomass_stores, "e_initial"] = biomass_potential + + # set storage bidding prices + # TODO bidding prices for Fischer-Tropsch, Methanation for H2? + bidding_prices = config["bidding_prices"] + for c in n.iterate_components({"Store", "Link", "StorageUnit"}): + c.df.marginal_cost.update(c.df.carrier.map(bidding_prices).dropna()) + + # deduct industry solid biomass + assert sum(n.stores.carriers == "solid biomass") == 1, "Myopic operation not implemented for spatially resolved solid biomass." + n.stores.at["EU solid biomass", "e_initial"] -= n.loads.at["solid biomass for industry", "p_set"] * 8760 + n.remove("Load", "solid biomass for industry") + + return n + + +def solve_network_myopic(n, config, opts='', **kwargs): + + rolling_horizon = config["operations"]["rolling_horizon"] + + freq = int(pd.infer_freq(n.snapshots)[:-1]) + window = rolling_horizon["window"] * 24 // freq + overlap = rolling_horizon["overlap"] * 24 // freq + kept = window - overlap + length = len(n.snapshots) + + assert kept > 0, f"Overlap ({overlap} days) must be smaller than windows ({window} days)." + + for i in range(length // kept): + + snapshots = n.snapshots[i * kept:(i + 1) * kept + overlap] + logger.info(f"Optimising operations from {snapshots[0]} to {snapshots[-1]}") + + n = solve_network(n, config, opts=opts, snapshots=snapshots, **kwargs) + + last_kept = n.snapshots[(i + 1) * kept - 1] + logger.info(f"Setting initial SOCs from {last_kept} for next iteration.\n") + + n.stores.e_initial = n.stores_t.e.loc[last_kept] + n.storage_units.state_of_charge_initial = n.storage_units_t.state_of_charge.loc[last_kept] + + return n + + +if __name__ == "__main__": + if 'snakemake' not in globals(): + from helper import mock_snakemake + snakemake = mock_snakemake( + 'solve_operations_network_myopic', + capacity_year=1952, + simpl='', + opts='', + clusters=37, + lv=2.0, + sector_opts='Co2L0-25H-T-H-B-I-A', + planning_horizons=2030, + weather_year=2013 + ) + + logging.basicConfig(filename=snakemake.log.python, + level=snakemake.config['logging_level']) + + tmpdir = snakemake.config['solving'].get('tmpdir') + if tmpdir is not None: + from pathlib import Path + Path(tmpdir).mkdir(parents=True, exist_ok=True) + + config = snakemake.config["operations"] + overrides = override_component_attrs(snakemake.input.overrides) + + n = pypsa.Network(snakemake.input.pre, override_component_attrs=overrides) + + n_post = pypsa.Network(snakemake.input.post, override_component_attrs=overrides) + n = set_parameters_from_optimized(n, n_post) + del n_post + + n_previous = pypsa.Network(snakemake.input.previous, override_component_attrs=overrides) + store_soc = n_previous.stores_t.e.iloc[-1] + storage_unit_soc = n_previous.storage_units_t.state_of_charge.iloc[-1] + del n_previous + + n = remove_unused_components(n) + n = add_load_shedding(n) + n = prepare_myopic(n, config, store_soc, storage_unit_soc) + + opts = snakemake.wildcards.sector_opts.split('-') + solve_opts = snakemake.config['solving']['options'] + solve_opts['skip_iterations'] = True + + n = prepare_network(n, solve_opts) + + n = solve_network_myopic( + n, + config=snakemake.config, + opts=opts, + solver_dir=tmpdir, + solver_logfile=snakemake.log.solver + ) + + n.export_to_netcdf(snakemake.output[0])