From 08c4643b690dce71c4eb11bbac79a8e075f8da4f Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 18 Jul 2023 18:15:40 +0300 Subject: [PATCH 01/57] Initial update comment --- scripts/solve_network.py | 238 +++++++++++++++++++++------------------ 1 file changed, 130 insertions(+), 108 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index e8f54a00f..95dc3f521 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -86,16 +86,10 @@ import pandas as pd import pypsa from _helpers import configure_logging +from linopy import merge from pypsa.descriptors import get_switchable_as_dense as get_as_dense -from pypsa.linopf import ( - define_constraints, - define_variables, - get_var, - ilopf, - join_exprs, - linexpr, - network_lopf, -) +from pypsa.optimization.abstract import optimize_transmission_expansion_iteratively +from pypsa.optimization.optimize import optimize from vresutils.benchmark import memory_logger logger = logging.getLogger(__name__) @@ -166,29 +160,32 @@ def add_CCL_constraints(n, config): ) gen_country = n.generators.bus.map(n.buses.country) - # cc means country and carrier - p_nom_per_cc = ( - pd.DataFrame( - { - "p_nom": linexpr((1, get_var(n, "Generator", "p_nom"))), - "country": gen_country, - "carrier": n.generators.carrier, - } + capacity_variable = n.model["Generator-p_nom"] + + lhs = [] + ext_carriers = n.generators.query("p_nom_extendable").carrier.unique() + for c in ext_carriers: + ext_carrier = n.generators.query("p_nom_extendable and carrier == @c") + country_grouper = ( + ext_carrier.bus.map(n.buses.country) + .rename_axis("Generator-ext") + .rename("country") ) - .dropna(subset=["p_nom"]) - .groupby(["country", "carrier"]) - .p_nom.apply(join_exprs) + ext_carrier_per_country = capacity_variable.loc[ + country_grouper.index + ].groupby_sum(country_grouper) + lhs.append(ext_carrier_per_country) + lhs = merge(lhs, dim=pd.Index(ext_carriers, name="carrier")) + + min_matrix = agg_p_nom_minmax["min"].to_xarray().unstack().reindex_like(lhs) + max_matrix = agg_p_nom_minmax["max"].to_xarray().unstack().reindex_like(lhs) + + n.model.add_constraints( + lhs >= min_matrix, name="agg_p_nom_min", mask=min_matrix.notnull() + ) + n.model.add_constraints( + lhs <= max_matrix, name="agg_p_nom_max", mask=max_matrix.notnull() ) - minimum = agg_p_nom_minmax["min"].dropna() - if not minimum.empty: - minconstraint = define_constraints( - n, p_nom_per_cc[minimum.index], ">=", minimum, "agg_p_nom", "min" - ) - maximum = agg_p_nom_minmax["max"].dropna() - if not maximum.empty: - maxconstraint = define_constraints( - n, p_nom_per_cc[maximum.index], "<=", maximum, "agg_p_nom", "max" - ) def add_EQ_constraints(n, o, scaling=1e-1): @@ -212,76 +209,89 @@ def add_EQ_constraints(n, o, scaling=1e-1): ) inflow = inflow.reindex(load.index).fillna(0.0) rhs = scaling * (level * load - inflow) + dispatch_variable = n.model["Generator-p"].T lhs_gen = ( - linexpr( - (n.snapshot_weightings.generators * scaling, get_var(n, "Generator", "p").T) - ) - .T.groupby(ggrouper, axis=1) - .apply(join_exprs) + (dispatch_variable * (n.snapshot_weightings.generators * scaling)) + .groupby_sum(ggrouper) + .sum("snapshot") ) - lhs_spill = ( - linexpr( - ( - -n.snapshot_weightings.stores * scaling, - get_var(n, "StorageUnit", "spill").T, - ) + if not n.storage_units_t.inflow.empty: + spillage_variable = n.model["StorageUnit-spill"] + lhs_spill = ( + (spillage_variable * (-n.snapshot_weightings.stores * scaling)) + .groupby_sum(sgrouper) + .sum("snapshot") ) - .T.groupby(sgrouper, axis=1) - .apply(join_exprs) - ) - lhs_spill = lhs_spill.reindex(lhs_gen.index).fillna("") - lhs = lhs_gen + lhs_spill - define_constraints(n, lhs, ">=", rhs, "equity", "min") + lhs = merge(lhs_gen, lhs_spill) + else: + lhs = lhs_gen + n.model.add_constraints(lhs >= rhs, name="equity_min") def add_BAU_constraints(n, config): mincaps = pd.Series(config["electricity"]["BAU_mincapacities"]) - lhs = ( - linexpr((1, get_var(n, "Generator", "p_nom"))) - .groupby(n.generators.carrier) - .apply(join_exprs) - ) - define_constraints(n, lhs, ">=", mincaps[lhs.index], "Carrier", "bau_mincaps") + # lhs = ( + # linexpr((1, get_var(n, "Generator", "p_nom"))) + # .groupby(n.generators.carrier) + # .apply(join_exprs) + # ) + # define_constraints(n, lhs, ">=", mincaps[lhs.index], "Carrier", "bau_mincaps") + capacity_variable = n.model["Generator-p_nom"] + ext_i = n.generators.query("p_nom_extendable") + ext_carrier_i = ext_i.carrier.rename_axis("Generator-ext") + lhs = capacity_variable.groupby_sum(ext_carrier_i) + rhs = mincaps[lhs.coords["carrier"].values].rename_axis("carrier") + n.model.add_constraints(lhs >= rhs, name="bau_mincaps") def add_SAFE_constraints(n, config): - peakdemand = ( - 1.0 + config["electricity"]["SAFE_reservemargin"] - ) * n.loads_t.p_set.sum(axis=1).max() + peakdemand = n.loads_t.p_set.sum(axis=1).max() + margin = 1.0 + config["electricity"]["SAFE_reservemargin"] + reserve_margin = peakdemand * margin conv_techs = config["plotting"]["conv_techs"] + ext_gens_i = n.generators.query("carrier in @conv_techs & p_nom_extendable").index + capacity_variable = n.model["Generator-p_nom"] + ext_cap_var = capacity_variable.sel({"Generator-ext": ext_gens_i}) + lhs = ext_cap_var.sum() exist_conv_caps = n.generators.query( "~p_nom_extendable & carrier in @conv_techs" ).p_nom.sum() - ext_gens_i = n.generators.query("carrier in @conv_techs & p_nom_extendable").index - lhs = linexpr((1, get_var(n, "Generator", "p_nom")[ext_gens_i])).sum() - rhs = peakdemand - exist_conv_caps - define_constraints(n, lhs, ">=", rhs, "Safe", "mintotalcap") + rhs = reserve_margin - exist_conv_caps + n.model.add_constraints(lhs >= rhs, name="safe_mintotalcap") -def add_operational_reserve_margin_constraint(n, config): +def add_operational_reserve_margin_constraint(n, sns, config): reserve_config = config["electricity"]["operational_reserve"] EPSILON_LOAD = reserve_config["epsilon_load"] EPSILON_VRES = reserve_config["epsilon_vres"] CONTINGENCY = reserve_config["contingency"] # Reserve Variables - reserve = get_var(n, "Generator", "r") - lhs = linexpr((1, reserve)).sum(1) + n.model.add_variables( + 0, np.inf, coords=[sns, n.generators.index], name="Generator-r" + ) + reserve = n.model["Generator-r"] + lhs = reserve.sum("Generator") # Share of extendable renewable capacities ext_i = n.generators.query("p_nom_extendable").index vres_i = n.generators_t.p_max_pu.columns if not ext_i.empty and not vres_i.empty: capacity_factor = n.generators_t.p_max_pu[vres_i.intersection(ext_i)] - renewable_capacity_variables = get_var(n, "Generator", "p_nom")[ - vres_i.intersection(ext_i) - ] - lhs += linexpr( - (-EPSILON_VRES * capacity_factor, renewable_capacity_variables) - ).sum(1) + renewable_capacity_variables = ( + n.model["Generator-p_nom"] + .sel({"Generator-ext": vres_i.intersection(ext_i)}) + .rename({"Generator-ext": "Generator"}) + ) + lhs = merge( + lhs, + (renewable_capacity_variables * (-EPSILON_VRES * capacity_factor)).sum( + ["Generator"] + ), + ) - # Total demand at t - demand = n.loads_t.p.sum(1) + # Total demand per t + demand = n.loads_t.p_set.sum(1) # VRES potential of non extendable generators capacity_factor = n.generators_t.p_max_pu[vres_i.difference(ext_i)] @@ -291,7 +301,8 @@ def add_operational_reserve_margin_constraint(n, config): # Right-hand-side rhs = EPSILON_LOAD * demand + EPSILON_VRES * potential + CONTINGENCY - define_constraints(n, lhs, ">=", rhs, "Reserve margin") + # define_constraints(n, lhs, ">=", rhs, "Reserve margin") + n.model.add_constraints(lhs >= rhs, name="reserve_margin") def update_capacity_constraint(n): @@ -299,24 +310,27 @@ def update_capacity_constraint(n): ext_i = n.generators.query("p_nom_extendable").index fix_i = n.generators.query("not p_nom_extendable").index - dispatch = get_var(n, "Generator", "p") - reserve = get_var(n, "Generator", "r") - - capacity_fixed = n.generators.p_nom[fix_i] - + dispatch = n.model["Generator-p"] + reserve = n.model["Generator-r"] p_max_pu = get_as_dense(n, "Generator", "p_max_pu") + capacity_fixed = n.generators.p_nom[fix_i] - lhs = linexpr((1, dispatch), (1, reserve)) + lhs = merge( + dispatch * 1, + reserve * 1, + ) if not ext_i.empty: - capacity_variable = get_var(n, "Generator", "p_nom") - lhs += linexpr((-p_max_pu[ext_i], capacity_variable)).reindex( - columns=gen_i, fill_value="" + capacity_variable = n.model["Generator-p_nom"] + lhs = merge( + lhs, + capacity_variable.rename({"Generator-ext": "Generator"}) * -p_max_pu[ext_i], ) - rhs = (p_max_pu[fix_i] * capacity_fixed).reindex(columns=gen_i, fill_value=0) - - define_constraints(n, lhs, "<=", rhs, "Generators", "updated_capacity_constraint") + rhs = (p_max_pu[fix_i] * capacity_fixed).reindex(columns=gen_i) + n.model.add_constraints( + lhs <= rhs, name="gen_updated_capacity_constraint", mask=rhs.notnull() + ) def add_operational_reserve_margin(n, sns, config): @@ -325,26 +339,25 @@ def add_operational_reserve_margin(n, sns, config): https://genxproject.github.io/GenX/dev/core/#Reserves. """ - define_variables(n, 0, np.inf, "Generator", "r", axes=[sns, n.generators.index]) - - add_operational_reserve_margin_constraint(n, config) - + # define_variables(n, 0, np.inf, "Generator", "r", axes=[sns, n.generators.index]) + # add_operational_reserve_margin_constraint(n, config) + add_operational_reserve_margin_constraint(n, sns, config) update_capacity_constraint(n) def add_battery_constraints(n): nodes = n.buses.index[n.buses.carrier == "battery"] - if nodes.empty or ("Link", "p_nom") not in n.variables.index: + # if nodes.empty or ("Link", "p_nom") not in n.variables.index: + if nodes.empty: return - link_p_nom = get_var(n, "Link", "p_nom") - lhs = linexpr( - (1, link_p_nom[nodes + " charger"]), - ( - -n.links.loc[nodes + " discharger", "efficiency"].values, - link_p_nom[nodes + " discharger"].values, - ), + vars_link = n.model["Link-p_nom"] + eff = n.links.loc[nodes + " discharger", "efficiency"] + lhs = merge( + vars_link.sel({"Link-ext": nodes + " charger"}) * 1, + # for some reasons, eff is one element longer as compared with vars_link + vars_link.sel({"Link-ext": nodes + " discharger"}) * -eff[0], ) - define_constraints(n, lhs, "=", 0, "Link", "charger_ratio") + n.model.add_constraints(lhs == 0, name="link_charger_ratio") def extra_functionality(n, snapshots): @@ -384,16 +397,21 @@ def solve_network(n, config, opts="", **kwargs): n.config = config n.opts = opts - if cf_solving.get("skip_iterations", False): - network_lopf( + skip_iterations = cf_solving.get("skip_iterations", False) + if not n.lines.s_nom_extendable.any(): + skip_iterations = True + logger.info("No expandable lines found. Skipping iterative solving.") + + if skip_iterations: + optimize( n, solver_name=solver_name, solver_options=solver_options, extra_functionality=extra_functionality, - **kwargs + **kwargs, ) else: - ilopf( + optimize_transmission_expansion_iteratively( n, solver_name=solver_name, solver_options=solver_options, @@ -401,7 +419,7 @@ def solve_network(n, config, opts="", **kwargs): min_iterations=min_iterations, max_iterations=max_iterations, extra_functionality=extra_functionality, - **kwargs + **kwargs, ) return n @@ -429,15 +447,19 @@ def solve_network(n, config, opts="", **kwargs): fn = getattr(snakemake.log, "memory", None) with memory_logger(filename=fn, interval=30.0) as mem: n = pypsa.Network(snakemake.input[0]) - if snakemake.config["augmented_line_connection"].get("add_to_snakefile"): - n.lines.loc[ - n.lines.index.str.contains("new"), "s_nom_min" - ] = snakemake.config["augmented_line_connection"].get("min_expansion") + # TODO Double-check handling the augmented case + # if snakemake.config["augmented_line_connection"].get("add_to_snakefile"): + # n.lines.loc[ + # n.lines.index.str.contains("new"), "s_nom_min" + # ] = snakemake.config["augmented_line_connection"].get("min_expansion") n = prepare_network(n, solve_opts) n = solve_network( n, - config=snakemake.config, - opts=opts, + # TODO The initial version looks clearer... + # config=snakemake.config, + # opts=opts, + snakemake.config, + opts, solver_dir=tmpdir, solver_logfile=snakemake.log.solver, ) From 28e11680ed7ef71c28a36a0520508a9817598237 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 22 Jul 2023 19:29:02 +0300 Subject: [PATCH 02/57] Align with PyPSA-Eur approach --- scripts/solve_network.py | 43 +++++++++++++++------------------------- 1 file changed, 16 insertions(+), 27 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 95dc3f521..75d5ffab0 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -95,7 +95,7 @@ logger = logging.getLogger(__name__) -def prepare_network(n, solve_opts): +def prepare_network(n, solve_opts, config=None): if "clip_p_max_pu" in solve_opts: for df in (n.generators_t.p_max_pu, n.storage_units_t.inflow): df.where(df > solve_opts["clip_p_max_pu"], other=0.0, inplace=True) @@ -438,32 +438,21 @@ def solve_network(n, config, opts="", **kwargs): ) configure_logging(snakemake) - tmpdir = snakemake.config["solving"].get("tmpdir") - if tmpdir is not None: - Path(tmpdir).mkdir(parents=True, exist_ok=True) + # tmpdir = snakemake.config["solving"].get("tmpdir") + # if tmpdir is not None: + # Path(tmpdir).mkdir(parents=True, exist_ok=True) opts = snakemake.wildcards.opts.split("-") solve_opts = snakemake.config["solving"]["options"] - fn = getattr(snakemake.log, "memory", None) - with memory_logger(filename=fn, interval=30.0) as mem: - n = pypsa.Network(snakemake.input[0]) - # TODO Double-check handling the augmented case - # if snakemake.config["augmented_line_connection"].get("add_to_snakefile"): - # n.lines.loc[ - # n.lines.index.str.contains("new"), "s_nom_min" - # ] = snakemake.config["augmented_line_connection"].get("min_expansion") - n = prepare_network(n, solve_opts) - n = solve_network( - n, - # TODO The initial version looks clearer... - # config=snakemake.config, - # opts=opts, - snakemake.config, - opts, - solver_dir=tmpdir, - solver_logfile=snakemake.log.solver, - ) - n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) - n.export_to_netcdf(snakemake.output[0]) - - logger.info("Maximum memory usage: {}".format(mem.mem_usage)) + n = pypsa.Network(snakemake.input[0]) + # TODO Double-check handling the augmented case + # if snakemake.config["augmented_line_connection"].get("add_to_snakefile"): + # n.lines.loc[ + # n.lines.index.str.contains("new"), "s_nom_min" + # ] = snakemake.config["augmented_line_connection"].get("min_expansion") + n = prepare_network(n, solve_opts, config=solve_opts) + n = solve_network( + n, config=snakemake.config, opts=opts, solver_logfile=snakemake.log.solver + ) + n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) + n.export_to_netcdf(snakemake.output[0]) From 4bb53276aa08c2583723b91f1b819653d760c4ba Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 17 Nov 2023 09:35:39 +0000 Subject: [PATCH 03/57] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- scripts/solve_network.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 3e075c258..59e6186e5 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -229,9 +229,8 @@ def add_EQ_constraints(n, o, scaling=1e-1): def add_BAU_constraints(n, config): - ext_c = n.generators.query("p_nom_extendable").carrier.unique() - + mincaps = pd.Series( config["electricity"].get("BAU_mincapacities", {key: 0 for key in ext_c}) ) @@ -584,7 +583,7 @@ def solve_network(n, config, opts="", **kwargs): n = pypsa.Network(snakemake.input[0]) # TODO Double-check handling the augmented case - #if snakemake.params.augmented_line_connection.get("add_to_snakefile"): + # if snakemake.params.augmented_line_connection.get("add_to_snakefile"): # n.lines.loc[ # n.lines.index.str.contains("new"), "s_nom_min" # ] = snakemake.params.augmented_line_connection.get("min_expansion") From 7ca9ee55169bb57d86e3b228bfe61e4f9e3b0752 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 15 Feb 2024 18:46:12 +0300 Subject: [PATCH 04/57] Fix merge --- scripts/solve_network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 420e4b25a..601e978d3 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -85,7 +85,7 @@ import numpy as np import pandas as pd import pypsa -from _helpers import configure_logging +from _helpers import configure_logging, create_logger from linopy import merge from pypsa.descriptors import get_switchable_as_dense as get_as_dense from pypsa.optimization.abstract import optimize_transmission_expansion_iteratively From 4c2d3461791cd77b2c6f164b670696f6af50f839 Mon Sep 17 00:00:00 2001 From: ekatef Date: Mon, 26 Feb 2024 01:34:16 +0300 Subject: [PATCH 05/57] Replace aggregategenerators with aggregateoneport --- scripts/simplify_network.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py index a150cd757..6bcd81e2a 100644 --- a/scripts/simplify_network.py +++ b/scripts/simplify_network.py @@ -102,7 +102,6 @@ from add_electricity import load_costs from cluster_network import cluster_regions, clustering_for_n_clusters from pypsa.clustering.spatial import ( - aggregategenerators, aggregateoneport, busmap_by_stubs, get_clustering_from_busmap, @@ -278,8 +277,12 @@ def replace_components(n, c, df, pnl): _, generator_strategies = get_aggregation_strategies(aggregation_strategies) carriers = set(n.generators.carrier) - set(exclude_carriers) - generators, generators_pnl = aggregategenerators( - n, busmap, carriers=carriers, custom_strategies=generator_strategies + generators, generators_pnl = aggregateoneport( + n, + busmap, + "Generator", + carriers=carriers, + custom_strategies=generator_strategies, ) replace_components(n, "Generator", generators, generators_pnl) From 97b24e6dc0e328c68fa1b748349fd4be7e864cfd Mon Sep 17 00:00:00 2001 From: ekatef Date: Mon, 26 Feb 2024 01:39:55 +0300 Subject: [PATCH 06/57] Add aggregation strategies as a parameter --- Snakefile | 1 + scripts/simplify_network.py | 14 +++++++++----- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/Snakefile b/Snakefile index 459834375..035e9eddd 100644 --- a/Snakefile +++ b/Snakefile @@ -555,6 +555,7 @@ rule add_electricity: rule simplify_network: params: + aggregation_strategies=config["cluster_options"]["aggregation_strategies"], renewable=config["renewable"], geo_crs=config["crs"]["geo_crs"], cluster_options=config["cluster_options"], diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py index 6bcd81e2a..4bfda9e8a 100644 --- a/scripts/simplify_network.py +++ b/scripts/simplify_network.py @@ -874,7 +874,7 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()): snakemake.params.config_links, snakemake.output, exclude_carriers, - aggregation_strategies, + snakemake.params.aggregation_strategies, ) busmaps = [trafo_map, simplify_links_map] @@ -891,12 +891,14 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()): hvdc_as_lines, lines_length_factor, snakemake.output, - aggregation_strategies=aggregation_strategies, + aggregation_strategies=snakemake.params.aggregation_strategies, ) busmaps.append(stub_map) if cluster_config.get("to_substations", False): - n, substation_map = aggregate_to_substations(n, aggregation_strategies) + n, substation_map = aggregate_to_substations( + n, snakemake.params.aggregation_strategies + ) busmaps.append(substation_map) # treatment of outliers (nodes without a profile for considered carrier): @@ -928,7 +930,9 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()): logger.info( f"clustering preparation (hac): aggregating {len(buses_i)} buses of type {carrier}." ) - n, busmap_hac = aggregate_to_substations(n, aggregation_strategies, buses_i) + n, busmap_hac = aggregate_to_substations( + n, params.aggregation_strategies, buses_i + ) busmaps.append(busmap_hac) if snakemake.wildcards.simpl: @@ -958,7 +962,7 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()): solver_name, cluster_config.get("algorithm", "hac"), cluster_config.get("feature", None), - aggregation_strategies, + aggregation_strategies=snakemake.params.aggregation_strategies, ) busmaps.append(cluster_map) From 004a8506cf43670307ea7331527707106801f006 Mon Sep 17 00:00:00 2001 From: ekatef Date: Mon, 26 Feb 2024 01:41:12 +0300 Subject: [PATCH 07/57] Re-define aggregation strategies --- scripts/simplify_network.py | 37 ++++++++++++++++++++----------------- 1 file changed, 20 insertions(+), 17 deletions(-) diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py index 4bfda9e8a..ea9efbd0c 100644 --- a/scripts/simplify_network.py +++ b/scripts/simplify_network.py @@ -274,7 +274,7 @@ def replace_components(n, c, df, pnl): _adjust_capital_costs_using_connection_costs(n, connection_costs_to_bus, output) - _, generator_strategies = get_aggregation_strategies(aggregation_strategies) + generator_strategies = aggregation_strategies["generators"] carriers = set(n.generators.carrier) - set(exclude_carriers) generators, generators_pnl = aggregateoneport( @@ -590,9 +590,9 @@ def aggregate_to_substations(n, aggregation_strategies=dict(), buses_i=None): if not dist.empty: busmap.loc[buses_i] = dist.idxmin(1) - bus_strategies, generator_strategies = get_aggregation_strategies( - aggregation_strategies - ) + line_strategies = aggregation_strategies.get("lines", dict()) + generator_strategies = aggregation_strategies.get("generators", dict()) + one_port_strategies = aggregation_strategies.get("one_ports", dict()) clustering = get_clustering_from_busmap( n, @@ -602,7 +602,9 @@ def aggregate_to_substations(n, aggregation_strategies=dict(), buses_i=None): aggregate_generators_carriers=None, aggregate_one_ports=["Load", "StorageUnit"], line_length_factor=1.0, + line_strategies=line_strategies, generator_strategies=generator_strategies, + one_port_strategies=one_port_strategies, scale_link_capital_costs=False, ) return clustering.network, busmap @@ -804,19 +806,20 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()): if (busmap.index == busmap).all(): return n, n.buses.index.to_series() - bus_strategies, generator_strategies = get_aggregation_strategies( - aggregation_strategies - ) + line_strategies = aggregation_strategies.get("lines", dict()) + generator_strategies = aggregation_strategies.get("generators", dict()) + one_port_strategies = aggregation_strategies.get("one_ports", dict()) clustering = get_clustering_from_busmap( n, busmap, - bus_strategies=bus_strategies, aggregate_generators_weighted=True, aggregate_generators_carriers=None, aggregate_one_ports=["Load", "StorageUnit"], line_length_factor=1.0, + line_strategies=line_strategies, generator_strategies=generator_strategies, + one_port_strategies=one_port_strategies, scale_link_capital_costs=False, ) @@ -846,14 +849,14 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()): "exclude_carriers", [] ) hvdc_as_lines = snakemake.params.electricity["hvdc_as_lines"] - aggregation_strategies = snakemake.params.cluster_options.get( - "aggregation_strategies", {} - ) - # translate str entries of aggregation_strategies to pd.Series functions: - aggregation_strategies = { - p: {k: getattr(pd.Series, v) for k, v in aggregation_strategies[p].items()} - for p in aggregation_strategies.keys() - } + # aggregation_strategies = snakemake.params.cluster_options.get( + # "aggregation_strategies", {} + # ) + ## translate str entries of aggregation_strategies to pd.Series functions: + # aggregation_strategies = { + # p: {k: getattr(pd.Series, v) for k, v in aggregation_strategies[p].items()} + # for p in aggregation_strategies.keys() + # } n, trafo_map = simplify_network_to_base_voltage(n, linetype, base_voltage) Nyears = n.snapshot_weightings.objective.sum() / 8760 @@ -989,7 +992,7 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()): n, merged_nodes_map = merge_isolated_nodes( n, threshold=p_threshold_merge_isolated, - aggregation_strategies=aggregation_strategies, + aggregation_strategies=snakemake.params.aggregation_strategies, ) busmaps.append(merged_nodes_map) From f19c69cc511f572c93420418552499736a55aa8a Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 27 Feb 2024 01:08:50 +0300 Subject: [PATCH 08/57] Update the environment --- envs/environment.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/envs/environment.yaml b/envs/environment.yaml index b65376cb7..06767535b 100644 --- a/envs/environment.yaml +++ b/envs/environment.yaml @@ -12,7 +12,7 @@ dependencies: - pip - mamba # esp for windows build -- pypsa>=0.24, <0.25 +- pypsa>=0.24, <0.26 # - atlite>=0.2.4 # until https://github.com/PyPSA/atlite/issues/244 is not merged - dask - powerplantmatching>=0.5.7 @@ -77,7 +77,7 @@ dependencies: # Default solver for tests (required for CI) - glpk -- ipopt<3.13.3 +- ipopt #<3.13.3 - gurobi - pip: From 3d72eee87685cec333810eaf2c320647767cbe9b Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 27 Feb 2024 01:09:20 +0300 Subject: [PATCH 09/57] A temporal fix for merge_isolated --- scripts/simplify_network.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py index ea9efbd0c..2d5a38fa3 100644 --- a/scripts/simplify_network.py +++ b/scripts/simplify_network.py @@ -988,13 +988,13 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()): p_threshold_merge_isolated = cluster_config.get("p_threshold_merge_isolated", False) n = drop_isolated_nodes(n, threshold=p_threshold_drop_isolated) - if p_threshold_merge_isolated: - n, merged_nodes_map = merge_isolated_nodes( - n, - threshold=p_threshold_merge_isolated, - aggregation_strategies=snakemake.params.aggregation_strategies, - ) - busmaps.append(merged_nodes_map) + # if p_threshold_merge_isolated: + # n, merged_nodes_map = merge_isolated_nodes( + # n, + # threshold=p_threshold_merge_isolated, + # aggregation_strategies=snakemake.params.aggregation_strategies, + # ) + # busmaps.append(merged_nodes_map) n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) n.export_to_netcdf(snakemake.output.network) From f94362b678f5bcce62d52c0a134ee455798c249f Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 27 Feb 2024 01:11:45 +0300 Subject: [PATCH 10/57] Replace pyomo methods with linopy ones when clustering --- scripts/cluster_network.py | 54 ++++++++++---------------------------- 1 file changed, 14 insertions(+), 40 deletions(-) diff --git a/scripts/cluster_network.py b/scripts/cluster_network.py index 6b559deba..63afcea95 100644 --- a/scripts/cluster_network.py +++ b/scripts/cluster_network.py @@ -128,7 +128,6 @@ import matplotlib.pyplot as plt import numpy as np import pandas as pd -import pyomo.environ as po import pypsa import seaborn as sns import shapely @@ -339,47 +338,22 @@ def distribute_clusters( distribution_factor.sum(), 1.0, rtol=1e-3 ), f"Country weights L must sum up to 1.0 when distributing clusters. Is {distribution_factor.sum()}." - m = po.ConcreteModel() - - def n_bounds(model, *n_id): - """ - Create a function that makes a bound pair for pyomo. - - Use n_bounds(model, n_id) if N is Single-Index - Use n_bounds(model, *n_id) if N is Multi-Index - Example: https://pyomo.readthedocs.io/en/stable/pyomo_modeling_components/Variables.html - - Returns - ------- - bounds = A function (or Python object) that gives a (lower,upper) bound pair i.e.(1,10) for the variable - """ - return (1, N[n_id]) - - m.n = po.Var(list(distribution_factor.index), bounds=n_bounds, domain=po.Integers) - m.tot = po.Constraint(expr=(po.summation(m.n) == n_clusters)) - m.objective = po.Objective( - expr=sum( - (m.n[i] - distribution_factor.loc[i] * n_clusters) ** 2 - for i in distribution_factor.index - ), - sense=po.minimize, + m = linopy.Model() + clusters = m.add_variables( + lower=1, upper=N, coords=[L.index], name="n", integer=True ) - - opt = po.SolverFactory(solver_name) - if not opt.has_capability("quadratic_objective"): - logger.warning( - f"The configured solver `{solver_name}` does not support quadratic objectives. Falling back to `ipopt`." + m.add_constraints(clusters.sum() == n_clusters, name="tot") + # leave out constant in objective (L * n_clusters) ** 2 + m.objective = (clusters * clusters - 2 * clusters * L * n_clusters).sum() + if solver_name == "gurobi": + logging.getLogger("gurobipy").propagate = False + elif solver_name != "scip": + logger.info( + f"The configured solver `{solver_name}` does not support quadratic objectives. Falling back to `scip`." ) - opt = po.SolverFactory("ipopt") - - results = opt.solve(m) - assert ( - results["Solver"][0]["Status"] == "ok" - ), f"Solver returned non-optimally: {results}" - - return ( - pd.Series(m.n.get_values(), index=distribution_factor.index).round().astype(int) - ) + solver_name = "scip" + m.solve(solver_name=solver_name) + return m.solution["n"].to_series().astype(int) def busmap_for_gadm_clusters(inputs, n, gadm_level, geo_crs, country_list): From 6b302f834ac1b263c61ab0baff1e1a86b0c10c19 Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 27 Feb 2024 01:14:15 +0300 Subject: [PATCH 11/57] Update aggregation strategies --- scripts/cluster_network.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/scripts/cluster_network.py b/scripts/cluster_network.py index 63afcea95..2edf92a0d 100644 --- a/scripts/cluster_network.py +++ b/scripts/cluster_network.py @@ -554,9 +554,9 @@ def clustering_for_n_clusters( extended_link_costs=0, focus_weights=None, ): - bus_strategies, generator_strategies = get_aggregation_strategies( - aggregation_strategies - ) + line_strategies = aggregation_strategies.get("lines", dict()) + generator_strategies = aggregation_strategies.get("generators", dict()) + one_port_strategies = aggregation_strategies.get("one_ports", dict()) if not isinstance(custom_busmap, pd.Series): if alternative_clustering: @@ -582,12 +582,13 @@ def clustering_for_n_clusters( clustering = get_clustering_from_busmap( n, busmap, - bus_strategies=bus_strategies, aggregate_generators_weighted=True, aggregate_generators_carriers=aggregate_carriers, aggregate_one_ports=["Load", "StorageUnit"], line_length_factor=line_length_factor, + line_strategies=line_strategies, generator_strategies=generator_strategies, + one_port_strategies=one_port_strategies, scale_link_capital_costs=False, ) From 28877d36d1bec4d6a642d2b7f657f2185869445d Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 29 Feb 2024 14:54:06 +0300 Subject: [PATCH 12/57] A temporal fix for the solver to keep the environment same --- scripts/cluster_network.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/scripts/cluster_network.py b/scripts/cluster_network.py index 2edf92a0d..f94eb8537 100644 --- a/scripts/cluster_network.py +++ b/scripts/cluster_network.py @@ -125,6 +125,7 @@ from functools import reduce import geopandas as gpd +import linopy import matplotlib.pyplot as plt import numpy as np import pandas as pd @@ -351,7 +352,8 @@ def distribute_clusters( logger.info( f"The configured solver `{solver_name}` does not support quadratic objectives. Falling back to `scip`." ) - solver_name = "scip" + # solver_name = "scip" + solver_name = "ipopt" m.solve(solver_name=solver_name) return m.solution["n"].to_series().astype(int) From 060b973d829e8e10fba8dfd9208910ed5c623179 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 29 Feb 2024 14:55:03 +0300 Subject: [PATCH 13/57] Update solving setup --- scripts/solve_network.py | 58 +++++++++++++++++++++++----------------- 1 file changed, 33 insertions(+), 25 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 57f68c5d7..208bfe605 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -518,42 +518,47 @@ def extra_functionality(n, snapshots): add_battery_constraints(n) -def solve_network(n, config, opts="", **kwargs): - solver_options = config["solving"]["solver"].copy() - solver_name = solver_options.pop("name") - cf_solving = config["solving"]["options"] - track_iterations = cf_solving.get("track_iterations", False) - min_iterations = cf_solving.get("min_iterations", 4) - max_iterations = cf_solving.get("max_iterations", 6) +def solve_network(n, config, solving, **kwargs): + set_of_options = solving["solver"]["options"] + cf_solving = solving["options"] + + kwargs["solver_options"] = ( + solving["solver_options"][set_of_options] if set_of_options else {} + ) + kwargs["solver_name"] = solving["solver"]["name"] # add to network for extra_functionality n.config = config - n.opts = opts skip_iterations = cf_solving.get("skip_iterations", False) if not n.lines.s_nom_extendable.any(): skip_iterations = True logger.info("No expandable lines found. Skipping iterative solving.") + # status, condition = run_glpk(**kwargs) + + status, condition = n.optimize.optimize_transmission_expansion_iteratively(**kwargs) + if skip_iterations: - optimize( - n, - solver_name=solver_name, - solver_options=solver_options, - extra_functionality=extra_functionality, - **kwargs, - ) + status, condition = n.optimize(**kwargs) else: - optimize_transmission_expansion_iteratively( - n, - solver_name=solver_name, - solver_options=solver_options, - track_iterations=track_iterations, - min_iterations=min_iterations, - max_iterations=max_iterations, - extra_functionality=extra_functionality, - **kwargs, + kwargs["track_iterations"] = (cf_solving.get("track_iterations", False),) + kwargs["min_iterations"] = (cf_solving.get("min_iterations", 4),) + kwargs["max_iterations"] = (cf_solving.get("max_iterations", 6),) + status, condition = n.optimize.optimize_transmission_expansion_iteratively( + **kwargs + ) + + if status != "ok": # and not rolling_horizon: + logger.warning( + f"Solving status '{status}' with termination condition '{condition}'" ) + if "infeasible" in condition: + labels = n.model.compute_infeasibilities() + logger.info(f"Labels:\n{labels}") + n.model.print_infeasibilities() + raise RuntimeError("Solving status 'infeasible'") + return n @@ -585,7 +590,10 @@ def solve_network(n, config, opts="", **kwargs): # ] = snakemake.params.augmented_line_connection.get("min_expansion") n = prepare_network(n, solve_opts, config=solve_opts) n = solve_network( - n, config=snakemake.config, opts=opts, solver_logfile=snakemake.log.solver + n, + config=snakemake.config, + solving=snakemake.params.solving, + solver_logfile=snakemake.log.solver, ) n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) n.export_to_netcdf(snakemake.output[0]) From e953d5ab284854e4fd94384981d6dc5f0c166f57 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 29 Feb 2024 14:56:32 +0300 Subject: [PATCH 14/57] Update solver settings in the config --- config.default.yaml | 57 +++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 53 insertions(+), 4 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index 254caccd1..48f04a8cb 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -386,15 +386,64 @@ solving: #nhours: 10 solver: name: gurobi + options: gurobi-default + +solver_options: + highs-default: + # refer to https://ergo-code.github.io/HiGHS/options/definitions.html#solver + threads: 4 + solver: "ipm" + run_crossover: "off" + small_matrix_value: 1e-6 + large_matrix_value: 1e9 + primal_feasibility_tolerance: 1e-5 + dual_feasibility_tolerance: 1e-5 + ipm_optimality_tolerance: 1e-4 + parallel: "on" + random_seed: 123 + gurobi-default: threads: 4 - method: 2 # barrier (=ipm) + method: 2 # barrier crossover: 0 - BarConvTol: 1.e-5 - FeasibilityTol: 1.e-6 + BarConvTol: 1.e-6 + Seed: 123 AggFill: 0 PreDual: 0 GURO_PAR_BARDENSETHRESH: 200 - + gurobi-numeric-focus: + name: gurobi + NumericFocus: 3 # Favour numeric stability over speed + method: 2 # barrier + crossover: 0 # do not use crossover + BarHomogeneous: 1 # Use homogeneous barrier if standard does not converge + BarConvTol: 1.e-5 + FeasibilityTol: 1.e-4 + OptimalityTol: 1.e-4 + ObjScale: -0.5 + threads: 8 + Seed: 123 + gurobi-fallback: # Use gurobi defaults + name: gurobi + crossover: 0 + method: 2 # barrier + BarHomogeneous: 1 # Use homogeneous barrier if standard does not converge + BarConvTol: 1.e-5 + FeasibilityTol: 1.e-5 + OptimalityTol: 1.e-5 + Seed: 123 + threads: 8 + cplex-default: + threads: 4 + lpmethod: 4 # barrier + solutiontype: 2 # non basic solution, ie no crossover + barrier.convergetol: 1.e-5 + feasopt.tolerance: 1.e-6 + copt-default: + Threads: 8 + LpMethod: 2 + Crossover: 0 + cbc-default: {} # Used in CI + glpk-default: {} # Used in CI plotting: map: From 5fd80606a7375bb44c415c4bb619764043d92cd1 Mon Sep 17 00:00:00 2001 From: ekatef Date: Wed, 6 Mar 2024 15:17:09 +0300 Subject: [PATCH 15/57] Update environment --- envs/environment.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/envs/environment.yaml b/envs/environment.yaml index 06767535b..e1a1ff955 100644 --- a/envs/environment.yaml +++ b/envs/environment.yaml @@ -27,6 +27,7 @@ dependencies: - memory_profiler - ruamel.yaml<=0.17.26 - pytables +- pyscipopt # added to compy with the quadratic objective requirement of the clustering script - lxml - numpy - pandas From f86da66b2d3209af1e617f87bfc727fd45fd3a1a Mon Sep 17 00:00:00 2001 From: ekatef Date: Wed, 6 Mar 2024 15:17:34 +0300 Subject: [PATCH 16/57] Update a solver for clustering --- scripts/cluster_network.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/scripts/cluster_network.py b/scripts/cluster_network.py index f94eb8537..78f19a060 100644 --- a/scripts/cluster_network.py +++ b/scripts/cluster_network.py @@ -352,8 +352,7 @@ def distribute_clusters( logger.info( f"The configured solver `{solver_name}` does not support quadratic objectives. Falling back to `scip`." ) - # solver_name = "scip" - solver_name = "ipopt" + solver_name = "scip" m.solve(solver_name=solver_name) return m.solution["n"].to_series().astype(int) From a2a48f72678ab0979b72b878f609b4266a517c97 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 4 Apr 2024 16:39:47 +0300 Subject: [PATCH 17/57] Define bus aggregation strategies --- scripts/cluster_network.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/scripts/cluster_network.py b/scripts/cluster_network.py index 78f19a060..1a3b360da 100644 --- a/scripts/cluster_network.py +++ b/scripts/cluster_network.py @@ -589,6 +589,13 @@ def clustering_for_n_clusters( line_length_factor=line_length_factor, line_strategies=line_strategies, generator_strategies=generator_strategies, + bus_strategies={ + "lat": "mean", + "lon": "mean", + "tag_substation": "first", + "tag_area": "first", + "country": "first", + }, one_port_strategies=one_port_strategies, scale_link_capital_costs=False, ) From 2be7fe1f14f96b233a6a62750e4c3bb379e9252e Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 4 Apr 2024 17:29:27 +0300 Subject: [PATCH 18/57] Update environment --- envs/environment.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/envs/environment.yaml b/envs/environment.yaml index 308fc7c4e..907d7be5e 100644 --- a/envs/environment.yaml +++ b/envs/environment.yaml @@ -12,7 +12,7 @@ dependencies: - pip - mamba # esp for windows build -- pypsa>=0.24, <0.26 +- pypsa>=0.25, <0.26 # - atlite>=0.2.4 # until https://github.com/PyPSA/atlite/issues/244 is not merged - dask - powerplantmatching>=0.5.7 From 377455e9ec2be6bc3407d7ba3a4be588e2be2bb0 Mon Sep 17 00:00:00 2001 From: ekatef Date: Wed, 17 Apr 2024 00:00:45 +0300 Subject: [PATCH 19/57] Get back merge of isolated nodes --- scripts/simplify_network.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py index c9a120f78..49061a823 100644 --- a/scripts/simplify_network.py +++ b/scripts/simplify_network.py @@ -986,15 +986,14 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()): 0.0, cluster_config.get("p_threshold_drop_isolated", 0.0) ) p_threshold_merge_isolated = cluster_config.get("p_threshold_merge_isolated", False) - n = drop_isolated_nodes(n, threshold=p_threshold_drop_isolated) - # if p_threshold_merge_isolated: - # n, merged_nodes_map = merge_isolated_nodes( - # n, - # threshold=p_threshold_merge_isolated, - # aggregation_strategies=snakemake.params.aggregation_strategies, - # ) - # busmaps.append(merged_nodes_map) + if p_threshold_merge_isolated: + n, merged_nodes_map = merge_isolated_nodes( + n, + threshold=p_threshold_merge_isolated, + aggregation_strategies=snakemake.params.aggregation_strategies, + ) + busmaps.append(merged_nodes_map) n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) n.export_to_netcdf(snakemake.output.network) From a01f188e1f0763e2549c0454f98ca73ecd4a296b Mon Sep 17 00:00:00 2001 From: ekatef Date: Wed, 17 Apr 2024 00:04:30 +0300 Subject: [PATCH 20/57] Add keys to aggregation strategies --- scripts/cluster_network.py | 1 + scripts/simplify_network.py | 17 ++++++++++++++++- 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/scripts/cluster_network.py b/scripts/cluster_network.py index 1a3b360da..405e0a0a5 100644 --- a/scripts/cluster_network.py +++ b/scripts/cluster_network.py @@ -556,6 +556,7 @@ def clustering_for_n_clusters( focus_weights=None, ): line_strategies = aggregation_strategies.get("lines", dict()) + line_strategies.update({"geometry": "first", "bounds": "first"}) generator_strategies = aggregation_strategies.get("generators", dict()) one_port_strategies = aggregation_strategies.get("one_ports", dict()) diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py index 49061a823..8b36490fb 100644 --- a/scripts/simplify_network.py +++ b/scripts/simplify_network.py @@ -591,18 +591,25 @@ def aggregate_to_substations(n, aggregation_strategies=dict(), buses_i=None): busmap.loc[buses_i] = dist.idxmin(1) line_strategies = aggregation_strategies.get("lines", dict()) + line_strategies.update({"geometry": "first", "bounds": "first"}) generator_strategies = aggregation_strategies.get("generators", dict()) one_port_strategies = aggregation_strategies.get("one_ports", dict()) clustering = get_clustering_from_busmap( n, busmap, - bus_strategies=bus_strategies, aggregate_generators_weighted=True, aggregate_generators_carriers=None, aggregate_one_ports=["Load", "StorageUnit"], line_length_factor=1.0, line_strategies=line_strategies, + bus_strategies={ + "lat": "mean", + "lon": "mean", + "tag_substation": "first", + "tag_area": "first", + "country": "first", + }, generator_strategies=generator_strategies, one_port_strategies=one_port_strategies, scale_link_capital_costs=False, @@ -807,6 +814,7 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()): return n, n.buses.index.to_series() line_strategies = aggregation_strategies.get("lines", dict()) + line_strategies.update({"geometry": "first", "bounds": "first"}) generator_strategies = aggregation_strategies.get("generators", dict()) one_port_strategies = aggregation_strategies.get("one_ports", dict()) @@ -818,6 +826,13 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()): aggregate_one_ports=["Load", "StorageUnit"], line_length_factor=1.0, line_strategies=line_strategies, + bus_strategies={ + "lat": "mean", + "lon": "mean", + "tag_substation": "first", + "tag_area": "first", + "country": "first", + }, generator_strategies=generator_strategies, one_port_strategies=one_port_strategies, scale_link_capital_costs=False, From 644c2ea65521ad10a866dc03bc3c079360ec01ca Mon Sep 17 00:00:00 2001 From: ekatef Date: Wed, 17 Apr 2024 00:21:58 +0300 Subject: [PATCH 21/57] Update the tutorial config --- config.tutorial.yaml | 57 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) diff --git a/config.tutorial.yaml b/config.tutorial.yaml index 243d96ce9..242c4e245 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -384,6 +384,63 @@ solving: #nhours: 10 solver: name: glpk + options: glpk-default + + solver_options: + highs-default: + # refer to https://ergo-code.github.io/HiGHS/options/definitions.html#solver + threads: 4 + solver: "ipm" + run_crossover: "off" + small_matrix_value: 1e-6 + large_matrix_value: 1e9 + primal_feasibility_tolerance: 1e-5 + dual_feasibility_tolerance: 1e-5 + ipm_optimality_tolerance: 1e-4 + parallel: "on" + random_seed: 123 + gurobi-default: + threads: 4 + method: 2 # barrier + crossover: 0 + BarConvTol: 1.e-6 + Seed: 123 + AggFill: 0 + PreDual: 0 + GURO_PAR_BARDENSETHRESH: 200 + gurobi-numeric-focus: + name: gurobi + NumericFocus: 3 # Favour numeric stability over speed + method: 2 # barrier + crossover: 0 # do not use crossover + BarHomogeneous: 1 # Use homogeneous barrier if standard does not converge + BarConvTol: 1.e-5 + FeasibilityTol: 1.e-4 + OptimalityTol: 1.e-4 + ObjScale: -0.5 + threads: 8 + Seed: 123 + gurobi-fallback: # Use gurobi defaults + name: gurobi + crossover: 0 + method: 2 # barrier + BarHomogeneous: 1 # Use homogeneous barrier if standard does not converge + BarConvTol: 1.e-5 + FeasibilityTol: 1.e-5 + OptimalityTol: 1.e-5 + Seed: 123 + threads: 8 + cplex-default: + threads: 4 + lpmethod: 4 # barrier + solutiontype: 2 # non basic solution, ie no crossover + barrier.convergetol: 1.e-4 #5 + feasopt.tolerance: 1.e-5 #6 + cbc-default: {} # Used in CI + glpk-default: {} # Used in CI + + mem: 30000 #memory in MB; 20 GB enough for 50+B+I+H2; 100 GB for 181+B+I+H2 + walltime: "12:00:00" plotting: From d210ce10b74d84429d5a20ef0f13a1d1b09c7deb Mon Sep 17 00:00:00 2001 From: ekatef Date: Wed, 17 Apr 2024 01:42:48 +0300 Subject: [PATCH 22/57] Lift version constraints for pypsa --- envs/environment.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/envs/environment.yaml b/envs/environment.yaml index 907d7be5e..50ca9369d 100644 --- a/envs/environment.yaml +++ b/envs/environment.yaml @@ -12,7 +12,7 @@ dependencies: - pip - mamba # esp for windows build -- pypsa>=0.25, <0.26 +- pypsa # - atlite>=0.2.4 # until https://github.com/PyPSA/atlite/issues/244 is not merged - dask - powerplantmatching>=0.5.7 From c4016e257f22b73017c29e197394db81d5572099 Mon Sep 17 00:00:00 2001 From: Ekaterina Date: Thu, 9 May 2024 20:17:23 +0300 Subject: [PATCH 23/57] Add Davide's suggestion for calling the solver log Co-authored-by: Davide Fioriti <67809479+davide-f@users.noreply.github.com> --- scripts/solve_network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 208bfe605..60af9a36d 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -593,7 +593,7 @@ def solve_network(n, config, solving, **kwargs): n, config=snakemake.config, solving=snakemake.params.solving, - solver_logfile=snakemake.log.solver, + log_fn=snakemake.log.solver, ) n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) n.export_to_netcdf(snakemake.output[0]) From 5991ca723e3cb3cf7671aa87cf079f6ad5e77bdd Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 13 May 2024 07:05:50 +0000 Subject: [PATCH 24/57] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- scripts/simplify_network.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py index 79804dc74..5bfb5a11a 100644 --- a/scripts/simplify_network.py +++ b/scripts/simplify_network.py @@ -1134,9 +1134,9 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()): ) p_threshold_merge_isolated = cluster_config.get("p_threshold_merge_isolated", False) s_threshold_fetch_isolated = cluster_config.get("s_threshold_fetch_isolated", False) - + n = drop_isolated_nodes(n, threshold=p_threshold_drop_isolated) - + if p_threshold_merge_isolated: n, merged_nodes_map = merge_isolated_nodes( n, From 32674dbcd0fb01a3548c951c4dcc6f83239cc0f2 Mon Sep 17 00:00:00 2001 From: ekatef Date: Mon, 13 May 2024 10:06:33 +0300 Subject: [PATCH 25/57] Remove unprecified columns from lines component --- scripts/base_network.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/scripts/base_network.py b/scripts/base_network.py index fe91ada84..208535e20 100644 --- a/scripts/base_network.py +++ b/scripts/base_network.py @@ -517,6 +517,18 @@ def base_network( if hvdc_as_lines_config: lines = pd.concat([lines_ac, lines_dc]) + lines.drop( + labels=[ + "bus0_lon", + "bus0_lat", + "bus1_lon", + "bus1_lat", + "bus_0_coors", + "bus_1_coors", + ], + axis=1, + inplace=True, + ) n.import_components_from_dataframe(lines, "Line") else: lines_dc = _set_electrical_parameters_links(links_config, lines_dc) @@ -526,6 +538,18 @@ def base_network( axis=1, result_type="reduce", ) + lines_ac.drop( + labels=[ + "bus0_lon", + "bus0_lat", + "bus1_lon", + "bus1_lat", + "bus_0_coors", + "bus_1_coors", + ], + axis=1, + inplace=True, + ) n.import_components_from_dataframe(lines_ac, "Line") # The columns which names starts with "bus" are mixed up with the third-bus specification # when executing additional_linkports() From 7026712c98ba3492f260c9386ee6b7282424262d Mon Sep 17 00:00:00 2001 From: Davide Fioriti <67809479+davide-f@users.noreply.github.com> Date: Thu, 16 May 2024 10:10:28 +0200 Subject: [PATCH 26/57] Fix typo aggregation_strategies --- scripts/simplify_network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py index 5bfb5a11a..c9567d427 100644 --- a/scripts/simplify_network.py +++ b/scripts/simplify_network.py @@ -1149,7 +1149,7 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()): n, fetched_nodes_map = merge_into_network( n, threshold=s_threshold_fetch_isolated, - aggregation_strategies=aggregation_strategies, + aggregation_strategies=snakemake.params.aggregation_strategies, ) busmaps.append(fetched_nodes_map) From c444005d32ce0d09ff2bd3746b45ea106f4091fe Mon Sep 17 00:00:00 2001 From: Davide Fioriti Date: Thu, 16 May 2024 17:41:55 +0200 Subject: [PATCH 27/57] Remove get_aggregation_strategies --- scripts/_helpers.py | 25 ------------------------- scripts/cluster_network.py | 1 - scripts/simplify_network.py | 20 ++++++++++++++------ 3 files changed, 14 insertions(+), 32 deletions(-) diff --git a/scripts/_helpers.py b/scripts/_helpers.py index 4fdf000b9..7c736acd0 100644 --- a/scripts/_helpers.py +++ b/scripts/_helpers.py @@ -431,31 +431,6 @@ def dlProgress(count, blockSize, totalSize, roundto=roundto): urllib.request.urlretrieve(url, file, reporthook=dlProgress, data=data) -def get_aggregation_strategies(aggregation_strategies): - """ - Default aggregation strategies that cannot be defined in .yaml format must - be specified within the function, otherwise (when defaults are passed in - the function's definition) they get lost when custom values are specified - in the config. - """ - import numpy as np - - # to handle the new version of PyPSA. - try: - from pypsa.clustering.spatial import _make_consense - except Exception: - # TODO: remove after new release and update minimum pypsa version - from pypsa.clustering.spatial import _make_consense - - bus_strategies = dict(country=_make_consense("Bus", "country")) - bus_strategies.update(aggregation_strategies.get("buses", {})) - - generator_strategies = {"build_year": lambda x: 0, "lifetime": lambda x: np.inf} - generator_strategies.update(aggregation_strategies.get("generators", {})) - - return bus_strategies, generator_strategies - - def mock_snakemake(rulename, **wildcards): """ This function is expected to be executed from the "scripts"-directory of " diff --git a/scripts/cluster_network.py b/scripts/cluster_network.py index c04d626ba..09e4ffd49 100644 --- a/scripts/cluster_network.py +++ b/scripts/cluster_network.py @@ -133,7 +133,6 @@ REGION_COLS, configure_logging, create_logger, - get_aggregation_strategies, sets_path_to_root, update_p_nom_max, ) diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py index c9567d427..469a5d5cb 100644 --- a/scripts/simplify_network.py +++ b/scripts/simplify_network.py @@ -96,7 +96,6 @@ from _helpers import ( configure_logging, create_logger, - get_aggregation_strategies, update_p_nom_max, ) from add_electricity import load_costs @@ -859,19 +858,28 @@ def merge_into_network(n, threshold, aggregation_strategies=dict()): if (busmap.index == busmap).all(): return n, n.buses.index.to_series() - bus_strategies, generator_strategies = get_aggregation_strategies( - aggregation_strategies - ) - + line_strategies = aggregation_strategies.get("lines", dict()) + line_strategies.update({"geometry": "first", "bounds": "first"}) + generator_strategies = aggregation_strategies.get("generators", dict()) + one_port_strategies = aggregation_strategies.get("one_ports", dict()) + clustering = get_clustering_from_busmap( n, busmap, - bus_strategies=bus_strategies, aggregate_generators_weighted=True, aggregate_generators_carriers=None, aggregate_one_ports=["Load", "StorageUnit"], line_length_factor=1.0, + line_strategies=line_strategies, + bus_strategies={ + "lat": "mean", + "lon": "mean", + "tag_substation": "first", + "tag_area": "first", + "country": "first", + }, generator_strategies=generator_strategies, + one_port_strategies=one_port_strategies, scale_link_capital_costs=False, ) From a7db9bb0bf67b919efc731d7dca6db4f1e8b0f6e Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 16 May 2024 15:45:13 +0000 Subject: [PATCH 28/57] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- scripts/simplify_network.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py index 469a5d5cb..413b48556 100644 --- a/scripts/simplify_network.py +++ b/scripts/simplify_network.py @@ -93,11 +93,7 @@ import pandas as pd import pypsa import scipy as sp -from _helpers import ( - configure_logging, - create_logger, - update_p_nom_max, -) +from _helpers import configure_logging, create_logger, update_p_nom_max from add_electricity import load_costs from cluster_network import cluster_regions, clustering_for_n_clusters from pypsa.clustering.spatial import ( @@ -862,7 +858,7 @@ def merge_into_network(n, threshold, aggregation_strategies=dict()): line_strategies.update({"geometry": "first", "bounds": "first"}) generator_strategies = aggregation_strategies.get("generators", dict()) one_port_strategies = aggregation_strategies.get("one_ports", dict()) - + clustering = get_clustering_from_busmap( n, busmap, From 2074e18e8c11b2fc354abe69e5dc82a68d6a263e Mon Sep 17 00:00:00 2001 From: ekatef Date: Wed, 5 Jun 2024 22:07:33 +0300 Subject: [PATCH 29/57] Update the default config --- config.default.yaml | 111 +++++++++++++++++++++----------------------- 1 file changed, 54 insertions(+), 57 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index ca8ebbb1c..7fdd09afc 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -376,8 +376,6 @@ monte_carlo: type: beta args: [0.5, 2] - - solving: options: formulation: kirchhoff @@ -393,62 +391,61 @@ solving: name: gurobi options: gurobi-default -solver_options: - highs-default: + solver_options: + highs-default: # refer to https://ergo-code.github.io/HiGHS/options/definitions.html#solver - threads: 4 - solver: "ipm" - run_crossover: "off" - small_matrix_value: 1e-6 - large_matrix_value: 1e9 - primal_feasibility_tolerance: 1e-5 - dual_feasibility_tolerance: 1e-5 - ipm_optimality_tolerance: 1e-4 - parallel: "on" - random_seed: 123 - gurobi-default: - threads: 4 - method: 2 # barrier - crossover: 0 - BarConvTol: 1.e-6 - Seed: 123 - AggFill: 0 - PreDual: 0 - GURO_PAR_BARDENSETHRESH: 200 - gurobi-numeric-focus: - name: gurobi - NumericFocus: 3 # Favour numeric stability over speed - method: 2 # barrier - crossover: 0 # do not use crossover - BarHomogeneous: 1 # Use homogeneous barrier if standard does not converge - BarConvTol: 1.e-5 - FeasibilityTol: 1.e-4 - OptimalityTol: 1.e-4 - ObjScale: -0.5 - threads: 8 - Seed: 123 - gurobi-fallback: # Use gurobi defaults - name: gurobi - crossover: 0 - method: 2 # barrier - BarHomogeneous: 1 # Use homogeneous barrier if standard does not converge - BarConvTol: 1.e-5 - FeasibilityTol: 1.e-5 - OptimalityTol: 1.e-5 - Seed: 123 - threads: 8 - cplex-default: - threads: 4 - lpmethod: 4 # barrier - solutiontype: 2 # non basic solution, ie no crossover - barrier.convergetol: 1.e-5 - feasopt.tolerance: 1.e-6 - copt-default: - Threads: 8 - LpMethod: 2 - Crossover: 0 - cbc-default: {} # Used in CI - glpk-default: {} # Used in CI + threads: 4 + solver: "ipm" + run_crossover: "off" + small_matrix_value: 1e-6 + large_matrix_value: 1e9 + primal_feasibility_tolerance: 1e-5 + dual_feasibility_tolerance: 1e-5 + ipm_optimality_tolerance: 1e-4 + parallel: "on" + random_seed: 123 + gurobi-default: + threads: 4 + method: 2 # barrier + crossover: 0 + BarConvTol: 1.e-6 + Seed: 123 + AggFill: 0 + PreDual: 0 + GURO_PAR_BARDENSETHRESH: 200 + gurobi-numeric-focus: + name: gurobi + NumericFocus: 3 # Favour numeric stability over speed + method: 2 # barrier + crossover: 0 # do not use crossover + BarHomogeneous: 1 # Use homogeneous barrier if standard does not converge + BarConvTol: 1.e-5 + FeasibilityTol: 1.e-4 + OptimalityTol: 1.e-4 + ObjScale: -0.5 + threads: 8 + Seed: 123 + gurobi-fallback: # Use gurobi defaults + name: gurobi + crossover: 0 + method: 2 # barrier + BarHomogeneous: 1 # Use homogeneous barrier if standard does not converge + BarConvTol: 1.e-5 + FeasibilityTol: 1.e-5 + OptimalityTol: 1.e-5 + Seed: 123 + threads: 8 + cplex-default: + threads: 4 + lpmethod: 4 # barrier + solutiontype: 2 # non basic solution, ie no crossover + barrier.convergetol: 1.e-4 #5 + feasopt.tolerance: 1.e-5 #6 + cbc-default: {} # Used in CI + glpk-default: {} # Used in CI + + mem: 30000 #memory in MB; 20 GB enough for 50+B+I+H2; 100 GB for 181+B+I+H2 + walltime: "12:00:00" plotting: map: From f932c6a256d1f84d40035d3fea77db4c58f3d99f Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 6 Jun 2024 16:32:27 +0300 Subject: [PATCH 30/57] Remove debug parts --- scripts/solve_network.py | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 9d0d49ec8..d61e29f7d 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -526,17 +526,13 @@ def solve_network(n, config, solving, **kwargs): ) kwargs["solver_name"] = solving["solver"]["name"] - # add to network for extra_functionality - n.config = config - skip_iterations = cf_solving.get("skip_iterations", False) if not n.lines.s_nom_extendable.any(): skip_iterations = True logger.info("No expandable lines found. Skipping iterative solving.") - # status, condition = run_glpk(**kwargs) - - status, condition = n.optimize.optimize_transmission_expansion_iteratively(**kwargs) + # add to network for extra_functionality + n.config = config if skip_iterations: status, condition = n.optimize(**kwargs) @@ -575,9 +571,6 @@ def solve_network(n, config, solving, **kwargs): ) configure_logging(snakemake) - # tmpdir = snakemake.params.solving.get("tmpdir") - # if tmpdir is not None: - # Path(tmpdir).mkdir(parents=True, exist_ok=True) opts = snakemake.wildcards.opts.split("-") solve_opts = snakemake.params.solving["options"] From be6905db94cf1dc7709f330061743aa61258026f Mon Sep 17 00:00:00 2001 From: ekatef Date: Fri, 14 Jun 2024 18:52:09 +0300 Subject: [PATCH 31/57] Remove outdated comments --- scripts/solve_network.py | 9 --------- 1 file changed, 9 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index d61e29f7d..d9451bc16 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -229,12 +229,6 @@ def add_EQ_constraints(n, o, scaling=1e-1): def add_BAU_constraints(n, config): mincaps = pd.Series(config["electricity"]["BAU_mincapacities"]) - # lhs = ( - # linexpr((1, get_var(n, "Generator", "p_nom"))) - # .groupby(n.generators.carrier) - # .apply(join_exprs) - # ) - # define_constraints(n, lhs, ">=", mincaps[lhs.index], "Carrier", "bau_mincaps") capacity_variable = n.model["Generator-p_nom"] ext_i = n.generators.query("p_nom_extendable") ext_carrier_i = ext_i.carrier.rename_axis("Generator-ext") @@ -310,7 +304,6 @@ def add_operational_reserve_margin_constraint(n, sns, config): # Right-hand-side rhs = EPSILON_LOAD * demand + EPSILON_VRES * potential + CONTINGENCY - # define_constraints(n, lhs, ">=", rhs, "Reserve margin") n.model.add_constraints(lhs >= rhs, name="reserve_margin") @@ -348,8 +341,6 @@ def add_operational_reserve_margin(n, sns, config): https://genxproject.github.io/GenX/dev/core/#Reserves. """ - # define_variables(n, 0, np.inf, "Generator", "r", axes=[sns, n.generators.index]) - # add_operational_reserve_margin_constraint(n, config) add_operational_reserve_margin_constraint(n, sns, config) update_capacity_constraint(n) From fef6951801d93dc499731a938266b196dee471fb Mon Sep 17 00:00:00 2001 From: ekatef Date: Fri, 14 Jun 2024 18:56:33 +0300 Subject: [PATCH 32/57] Update implementation of the constraints --- scripts/solve_network.py | 46 ++++++++++++++++------------------------ 1 file changed, 18 insertions(+), 28 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index d9451bc16..94274c17e 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -229,35 +229,27 @@ def add_EQ_constraints(n, o, scaling=1e-1): def add_BAU_constraints(n, config): mincaps = pd.Series(config["electricity"]["BAU_mincapacities"]) - capacity_variable = n.model["Generator-p_nom"] + p_nom = n.model["Generator-p_nom"] ext_i = n.generators.query("p_nom_extendable") - ext_carrier_i = ext_i.carrier.rename_axis("Generator-ext") - lhs = capacity_variable.groupby_sum(ext_carrier_i) - rhs = mincaps[lhs.coords["carrier"].values].rename_axis("carrier") + ext_carrier_i = xr.DataArray(ext_i.carrier.rename_axis("Generator-ext")) + lhs = p_nom.groupby(ext_carrier_i).sum() + rhs = mincaps[lhs.indexes["carrier"]].rename_axis("carrier") n.model.add_constraints(lhs >= rhs, name="bau_mincaps") - maxcaps = pd.Series( - config["electricity"].get("BAU_maxcapacities", {key: np.inf for key in ext_c}) - ) - lhs = ( - linexpr((1, get_var(n, "Generator", "p_nom"))) - .groupby(n.generators.carrier) - .apply(join_exprs) - ) - define_constraints(n, lhs, "<=", maxcaps[lhs.index], "Carrier", "bau_maxcaps") - def add_SAFE_constraints(n, config): peakdemand = n.loads_t.p_set.sum(axis=1).max() margin = 1.0 + config["electricity"]["SAFE_reservemargin"] reserve_margin = peakdemand * margin - conv_techs = config["plotting"]["conv_techs"] - ext_gens_i = n.generators.query("carrier in @conv_techs & p_nom_extendable").index + conventional_carriers = config["electricity"]["conventional_carriers"] + ext_gens_i = n.generators.query( + "carrier in @conventional_carriers & p_nom_extendable" + ).index capacity_variable = n.model["Generator-p_nom"] - ext_cap_var = capacity_variable.sel({"Generator-ext": ext_gens_i}) - lhs = ext_cap_var.sum() + p_nom = n.model["Generator-p_nom"].loc[ext_gens_i] + lhs = p_nom.sum() exist_conv_caps = n.generators.query( - "~p_nom_extendable & carrier in @conv_techs" + "~p_nom_extendable & carrier in @conventional_carriers" ).p_nom.sum() rhs = reserve_margin - exist_conv_caps n.model.add_constraints(lhs >= rhs, name="safe_mintotalcap") @@ -283,7 +275,7 @@ def add_operational_reserve_margin_constraint(n, sns, config): capacity_factor = n.generators_t.p_max_pu[vres_i.intersection(ext_i)] renewable_capacity_variables = ( n.model["Generator-p_nom"] - .sel({"Generator-ext": vres_i.intersection(ext_i)}) + .loc[vres_i.intersection(ext_i)] .rename({"Generator-ext": "Generator"}) ) lhs = merge( @@ -294,12 +286,12 @@ def add_operational_reserve_margin_constraint(n, sns, config): ) # Total demand per t - demand = n.loads_t.p_set.sum(1) + 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(1) + potential = (capacity_factor * renewable_capacity).sum(axis=1) # Right-hand-side rhs = EPSILON_LOAD * demand + EPSILON_VRES * potential + CONTINGENCY @@ -314,9 +306,10 @@ def update_capacity_constraint(n): dispatch = n.model["Generator-p"] reserve = n.model["Generator-r"] - p_max_pu = get_as_dense(n, "Generator", "p_max_pu") capacity_fixed = n.generators.p_nom[fix_i] + p_max_pu = get_as_dense(n, "Generator", "p_max_pu") + lhs = merge( dispatch * 1, reserve * 1, @@ -324,12 +317,9 @@ def update_capacity_constraint(n): if not ext_i.empty: capacity_variable = n.model["Generator-p_nom"] - lhs = merge( - lhs, - capacity_variable.rename({"Generator-ext": "Generator"}) * -p_max_pu[ext_i], - ) + lhs = dispatch + reserve - capacity_variable * xr.DataArray(p_max_pu[ext_i]) - rhs = (p_max_pu[fix_i] * capacity_fixed).reindex(columns=gen_i) + rhs = (p_max_pu[fix_i] * capacity_fixed).reindex(columns=gen_i, fill_value=0) n.model.add_constraints( lhs <= rhs, name="gen_updated_capacity_constraint", mask=rhs.notnull() ) From a6ee0a5bf12336710208888c552bb83ecad2c47e Mon Sep 17 00:00:00 2001 From: ekatef Date: Fri, 14 Jun 2024 19:03:06 +0300 Subject: [PATCH 33/57] Maintain consistency in naming --- scripts/solve_network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 94274c17e..3d700a5f9 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -466,7 +466,7 @@ def add_RES_constraints(n, res_share): # todo: for links (lhs_charge and lhs_discharge), account for snapshot weightings lhs = lhs_gen + lhs_dispatch + lhs_store + lhs_charge + lhs_discharge - define_constraints(n, lhs, "=", rhs, "RES share") + define_constraints(n, lhs, "=", rhs, "res_share") def extra_functionality(n, snapshots): From 577464cd2544bf8cd3163cdd56d0821c58843995 Mon Sep 17 00:00:00 2001 From: ekatef Date: Fri, 14 Jun 2024 22:34:31 +0300 Subject: [PATCH 34/57] Update doc-strings --- scripts/solve_network.py | 112 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 107 insertions(+), 5 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 3d700a5f9..dadaceac1 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -52,15 +52,15 @@ linear optimal power flow (plus investment planning) is provided in the `documentation of PyPSA `_. -The optimization is based on the ``pyomo=False`` setting in the :func:`network.lopf` and :func:`pypsa.linopf.ilopf` function. -Additionally, some extra constraints specified in :mod:`prepare_network` are added. +The optimization is based on the :func:`network.optimize` function. +Additionally, some extra constraints specified in :mod:`prepare_network` and :mod:`solve_network` are added. Solving the network in multiple iterations is motivated through the dependence of transmission line capacities and impedances on values of corresponding flows. As lines are expanded their electrical parameters change, which renders the optimisation bilinear even if the power flow equations are linearized. To retain the computational advantage of continuous linear programming, a sequential linear programming technique is used, where in between iterations the line impedances are updated. -Details (and errors made through this heuristic) are discussed in the paper +Details (and errors introduced through this heuristic) are discussed in the paper - Fabian Neumann and Tom Brown. `Heuristics for Transmission Expansion Planning in Low-Carbon Energy System Models `_), *16th International Conference on the European Energy Market*, 2019. `arXiv:1907.10548 `_. @@ -144,6 +144,25 @@ def prepare_network(n, solve_opts, config=None): def add_CCL_constraints(n, config): + """ + Add CCL (country & carrier limit) constraint to the network. + + Add minimum and maximum levels of generator nominal capacity per carrier + for individual countries. Opts and path for agg_p_nom_minmax.csv must be defined + in config.yaml. Default file is available at data/agg_p_nom_minmax.csv. + + Parameters + ---------- + n : pypsa.Network + config : dict + + Example + ------- + scenario: + opts: [Co2L-CCL-24H] + electricity: + agg_p_nom_limits: data/agg_p_nom_minmax.csv + """ agg_p_nom_limits = config["electricity"].get("agg_p_nom_limits") try: @@ -188,6 +207,28 @@ def add_CCL_constraints(n, config): def add_EQ_constraints(n, o, scaling=1e-1): + """ + Add equity constraints to the network. + + Currently this is only implemented for the electricity sector only. + + Opts must be specified in the config.yaml. + + Parameters + ---------- + n : pypsa.Network + o : str + + Example + ------- + scenario: + opts: [Co2L-EQ0.7-24h] + + Require each country or node to on average produce a minimal share + of its total electricity consumption itself. Example: EQ0.7c demands each country + to produce on average at least 70% of its consumption; EQ0.7 demands + each node to produce on average at least 70% of its consumption. + """ float_regex = "[0-9]*\.?[0-9]+" level = float(re.findall(float_regex, o)[0]) if o[-1] == "c": @@ -228,6 +269,30 @@ def add_EQ_constraints(n, o, scaling=1e-1): def add_BAU_constraints(n, config): + """ + Add a per-carrier minimal overall capacity. + + BAU_mincapacities and opts must be adjusted in the config.yaml. + + Parameters + ---------- + n : pypsa.Network + 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 + """ mincaps = pd.Series(config["electricity"]["BAU_mincapacities"]) p_nom = n.model["Generator-p_nom"] ext_i = n.generators.query("p_nom_extendable") @@ -238,6 +303,25 @@ def add_BAU_constraints(n, config): def add_SAFE_constraints(n, config): + """ + Add a capacity reserve margin of a certain fraction above the peak demand. + Renewable generators and storage do not contribute. Ignores network. + + Parameters + ---------- + n : pypsa.Network + config : dict + + Example + ------- + config.yaml requires to specify opts: + + scenario: + opts: [Co2L-SAFE-24h] + electricity: + SAFE_reservemargin: 0.1 + Which sets a reserve margin of 10% above the peak demand. + """ peakdemand = n.loads_t.p_set.sum(axis=1).max() margin = 1.0 + config["electricity"]["SAFE_reservemargin"] reserve_margin = peakdemand * margin @@ -329,13 +413,31 @@ def add_operational_reserve_margin(n, sns, config): """ Build reserve margin constraints based on the formulation given in https://genxproject.github.io/GenX/dev/core/#Reserves. - """ + Parameters + ---------- + n : pypsa.Network + sns: pd.DatetimeIndex + config : dict + + Example: + -------- + config.yaml requires to specify operational_reserve: + operational_reserve: # like https://genxproject.github.io/GenX/dev/core/#Reserves + activate: true + epsilon_load: 0.02 # percentage of load at each snapshot + epsilon_vres: 0.02 # percentage of VRES at each snapshot + contingency: 400000 # MW + """ add_operational_reserve_margin_constraint(n, sns, config) update_capacity_constraint(n) def add_battery_constraints(n): + """ + Add constraint ensuring that charger = discharger, i.e. + 1 * charger_size - efficiency * discharger_size = 0 + """ nodes = n.buses.index[n.buses.carrier == "battery"] # if nodes.empty or ("Link", "p_nom") not in n.variables.index: if nodes.empty: @@ -472,7 +574,7 @@ def add_RES_constraints(n, res_share): def extra_functionality(n, snapshots): """ Collects supplementary constraints which will be passed to - ``pypsa.linopf.network_lopf``. + ``pypsa.optimization.optimize``. If you want to enforce additional custom constraints, this is a good location to add them. The arguments ``opts`` and ``snakemake.config`` are expected to be attached to the network. From 417bf62ca6361a23d6270c34a64da9e4f7afa176 Mon Sep 17 00:00:00 2001 From: ekatef Date: Fri, 14 Jun 2024 23:15:24 +0300 Subject: [PATCH 35/57] Update equity constraint --- scripts/solve_network.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index dadaceac1..c45506901 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -249,20 +249,24 @@ def add_EQ_constraints(n, o, scaling=1e-1): ) inflow = inflow.reindex(load.index).fillna(0.0) rhs = scaling * (level * load - inflow) - dispatch_variable = n.model["Generator-p"].T + dispatch_variable = n.model["Generator-p"] lhs_gen = ( (dispatch_variable * (n.snapshot_weightings.generators * scaling)) - .groupby_sum(ggrouper) + .groupby(ggrouper.to_xarray()) + .sum() .sum("snapshot") ) + # TODO: double check that this is really needed, why do have to subtract the spillage if not n.storage_units_t.inflow.empty: spillage_variable = n.model["StorageUnit-spill"] lhs_spill = ( (spillage_variable * (-n.snapshot_weightings.stores * scaling)) .groupby_sum(sgrouper) + .groupby(sgrouper.to_xarray()) + .sum() .sum("snapshot") ) - lhs = merge(lhs_gen, lhs_spill) + lhs = lhs_gen + lhs_spill else: lhs = lhs_gen n.model.add_constraints(lhs >= rhs, name="equity_min") From 6b060d064f1066b73988514394e3b24dfd4d4a30 Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 25 Jun 2024 00:02:00 +0300 Subject: [PATCH 36/57] Fix typos --- scripts/solve_network.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index c45506901..505472563 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -463,10 +463,10 @@ def add_RES_constraints(n, res_share): cgrouper = n.links.bus0.map(n.buses.country) logger.warning( - "The add_RES_constraints functionality is still work in progress. " - "Unexpected results might be incurred, particularly if " + "The add_RES_constraints() is still work in progress. " + "Unexpected results might be incurre, particularly if " "temporal clustering is applied or if an unexpected change of technologies " - "is subject to the obtimisation." + "is subject to the optimisation." ) load = ( From 26849ab2127e7228ed2f9c01e6a21d90cbee8640 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 30 Jun 2024 00:15:25 +0300 Subject: [PATCH 37/57] Remove hard-coding for renewable technologies --- scripts/solve_network.py | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 505472563..a62eb5e74 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -456,7 +456,7 @@ def add_battery_constraints(n): n.model.add_constraints(lhs == 0, name="link_charger_ratio") -def add_RES_constraints(n, res_share): +def add_RES_constraints(n, res_share, config): lgrouper = n.loads.bus.map(n.buses.country) ggrouper = n.generators.bus.map(n.buses.country) sgrouper = n.storage_units.bus.map(n.buses.country) @@ -476,15 +476,7 @@ def add_RES_constraints(n, res_share): rhs = res_share * load - res_techs = [ - "solar", - "onwind", - "offwind-dc", - "offwind-ac", - "battery", - "hydro", - "ror", - ] + res_techs = config.renewable_carriers charger = ["H2 electrolysis", "battery charger"] discharger = ["H2 fuel cell", "battery discharger"] @@ -597,7 +589,7 @@ def extra_functionality(n, snapshots): for o in opts: if "RES" in o: res_share = float(re.findall("[0-9]*\.?[0-9]+$", o)[0]) - add_RES_constraints(n, res_share) + add_RES_constraints(n, res_share, config) for o in opts: if "EQ" in o: add_EQ_constraints(n, o) From 60faf6316946435c05e745eb8aa84f2f7c4614ea Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 30 Jun 2024 00:29:13 +0300 Subject: [PATCH 38/57] Revise naming --- scripts/solve_network.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index a62eb5e74..af37318e7 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -476,12 +476,12 @@ def add_RES_constraints(n, res_share, config): rhs = res_share * load - res_techs = config.renewable_carriers + renew_techs = config.renewable_carriers charger = ["H2 electrolysis", "battery charger"] discharger = ["H2 fuel cell", "battery discharger"] - gens_i = n.generators.query("carrier in @res_techs").index - stores_i = n.storage_units.query("carrier in @res_techs").index + gens_i = n.generators.query("carrier in @renew_techs").index + stores_i = n.storage_units.query("carrier in @renew_techs").index charger_i = n.links.query("carrier in @charger").index discharger_i = n.links.query("carrier in @discharger").index From f54d1f686080212619a5b9a491ad731842ba0489 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 30 Jun 2024 00:29:44 +0300 Subject: [PATCH 39/57] Move sign under the common code chunks --- scripts/solve_network.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index af37318e7..7aa5f9c75 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -514,7 +514,7 @@ def add_RES_constraints(n, res_share, config): ( linexpr( ( - -n.snapshot_weightings.stores, + n.snapshot_weightings.stores, get_var(n, "StorageUnit", "p_store")[stores_i].T, ) ) @@ -532,7 +532,7 @@ def add_RES_constraints(n, res_share, config): ( linexpr( ( - -n.snapshot_weightings.stores, + n.snapshot_weightings.stores, get_var(n, "Link", "p")[charger_i].T, ) ) @@ -560,9 +560,8 @@ def add_RES_constraints(n, res_share, config): .fillna("") ) - # signs of resp. terms are coded in the linexpr. # todo: for links (lhs_charge and lhs_discharge), account for snapshot weightings - lhs = lhs_gen + lhs_dispatch + lhs_store + lhs_charge + lhs_discharge + lhs = lhs_gen + lhs_dispatch - lhs_store - lhs_charge + lhs_discharge define_constraints(n, lhs, "=", rhs, "res_share") From cbea30eb6a2c000a93493d85ad9c3adaa30da35a Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 30 Jun 2024 00:39:17 +0300 Subject: [PATCH 40/57] Add a missed efficiency --- scripts/solve_network.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 7aa5f9c75..d9ea28d92 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -532,7 +532,9 @@ def add_RES_constraints(n, res_share, config): ( linexpr( ( - n.snapshot_weightings.stores, + n.snapshot_weightings.stores.apply( + lambda r: r * n.links.loc[charger_i].efficiency + ), get_var(n, "Link", "p")[charger_i].T, ) ) From 6a7e770b9a76449bd9acf07d310f7fb15c5b3d73 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 30 Jun 2024 00:54:06 +0300 Subject: [PATCH 41/57] Select expressions --- scripts/solve_network.py | 83 ++++++++++++++++------------------------ 1 file changed, 33 insertions(+), 50 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index d9ea28d92..f74cd747c 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -494,70 +494,53 @@ def add_RES_constraints(n, res_share, config): .apply(join_exprs) ) - # StorageUnits - lhs_dispatch = ( + store_disp_expr = linexpr( ( - linexpr( - ( - n.snapshot_weightings.stores, - get_var(n, "StorageUnit", "p_dispatch")[stores_i].T, - ) - ) - .T.groupby(sgrouper, axis=1) - .apply(join_exprs) + n.snapshot_weightings.stores, + get_var(n, "StorageUnit", "p_dispatch")[stores_i].T, + ) + ) + store_expr = linexpr( + ( + n.snapshot_weightings.stores, + get_var(n, "StorageUnit", "p_store")[stores_i].T, ) + ) + charge_expr = linexpr( + ( + n.snapshot_weightings.stores.apply( + lambda r: r * n.links.loc[charger_i].efficiency + ), + get_var(n, "Link", "p")[charger_i].T, + ) + ) + discharge_expr = linexpr( + ( + n.snapshot_weightings.stores.apply( + lambda r: r * n.links.loc[discharger_i].efficiency + ), + get_var(n, "Link", "p")[discharger_i].T, + ) + ) + + # StorageUnits + lhs_dispatch = ( + (store_disp_expr.T.groupby(sgrouper, axis=1).apply(join_exprs)) .reindex(lhs_gen.index) .fillna("") ) - lhs_store = ( - ( - linexpr( - ( - n.snapshot_weightings.stores, - get_var(n, "StorageUnit", "p_store")[stores_i].T, - ) - ) - .T.groupby(sgrouper, axis=1) - .apply(join_exprs) - ) + (store_expr.T.groupby(sgrouper, axis=1).apply(join_exprs)) .reindex(lhs_gen.index) .fillna("") ) - - # Stores (or their resp. Link components) - # Note that the variables "p0" and "p1" currently do not exist. - # Thus, p0 and p1 must be derived from "p" (which exists), taking into account the link efficiency. lhs_charge = ( - ( - linexpr( - ( - n.snapshot_weightings.stores.apply( - lambda r: r * n.links.loc[charger_i].efficiency - ), - get_var(n, "Link", "p")[charger_i].T, - ) - ) - .T.groupby(cgrouper, axis=1) - .apply(join_exprs) - ) + (charge_expr.T.groupby(cgrouper, axis=1).apply(join_exprs)) .reindex(lhs_gen.index) .fillna("") ) - lhs_discharge = ( - ( - linexpr( - ( - n.snapshot_weightings.stores.apply( - lambda r: r * n.links.loc[discharger_i].efficiency - ), - get_var(n, "Link", "p")[discharger_i], - ) - ) - .groupby(cgrouper, axis=1) - .apply(join_exprs) - ) + (discharge_expr.T.groupby(cgrouper, axis=1).apply(join_exprs)) .reindex(lhs_gen.index) .fillna("") ) From e7bd8bb9f223a51a573c705a019fa64b7571b26c Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 30 Jun 2024 01:05:13 +0300 Subject: [PATCH 42/57] Put repetitions into a function --- scripts/solve_network.py | 34 ++++++++++++---------------------- 1 file changed, 12 insertions(+), 22 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index f74cd747c..7b0ab3941 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -523,29 +523,19 @@ def add_RES_constraints(n, res_share, config): ) ) - # StorageUnits - lhs_dispatch = ( - (store_disp_expr.T.groupby(sgrouper, axis=1).apply(join_exprs)) - .reindex(lhs_gen.index) - .fillna("") - ) - lhs_store = ( - (store_expr.T.groupby(sgrouper, axis=1).apply(join_exprs)) - .reindex(lhs_gen.index) - .fillna("") - ) - lhs_charge = ( - (charge_expr.T.groupby(cgrouper, axis=1).apply(join_exprs)) - .reindex(lhs_gen.index) - .fillna("") - ) - lhs_discharge = ( - (discharge_expr.T.groupby(cgrouper, axis=1).apply(join_exprs)) - .reindex(lhs_gen.index) - .fillna("") - ) + def form_lhs_definition(lin_expression, ngrouper, target_index=lhs_gen.index): + lhs_form = ( + (lin_expression.T.groupby(ngrouper, axis=1).apply(join_exprs)) + .reindex(target_index) + .fillna("") + ) + return lhs_form + + lhs_dispatch = form_lhs_definition(store_disp_expr, sgrouper) + lhs_store = form_lhs_definition(store_expr, sgrouper) + lhs_charge = form_lhs_definition(charge_expr, cgrouper) + lhs_discharge = form_lhs_definition(discharge_expr, cgrouper) - # todo: for links (lhs_charge and lhs_discharge), account for snapshot weightings lhs = lhs_gen + lhs_dispatch - lhs_store - lhs_charge + lhs_discharge define_constraints(n, lhs, "=", rhs, "res_share") From 0e50477bf8ec6d3a871e180d66635e871801ff3e Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 30 Jun 2024 01:05:46 +0300 Subject: [PATCH 43/57] Add a supplementary variable --- scripts/solve_network.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 7b0ab3941..e104a8e97 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -481,10 +481,13 @@ def add_RES_constraints(n, res_share, config): discharger = ["H2 fuel cell", "battery discharger"] gens_i = n.generators.query("carrier in @renew_techs").index + stores_i = n.storage_units.query("carrier in @renew_techs").index charger_i = n.links.query("carrier in @charger").index discharger_i = n.links.query("carrier in @discharger").index + stores_t_weights = n.snapshot_weightings.stores + # Generators lhs_gen = ( linexpr( @@ -496,29 +499,25 @@ def add_RES_constraints(n, res_share, config): store_disp_expr = linexpr( ( - n.snapshot_weightings.stores, + stores_t_weights, get_var(n, "StorageUnit", "p_dispatch")[stores_i].T, ) ) store_expr = linexpr( ( - n.snapshot_weightings.stores, + stores_t_weights, get_var(n, "StorageUnit", "p_store")[stores_i].T, ) ) charge_expr = linexpr( ( - n.snapshot_weightings.stores.apply( - lambda r: r * n.links.loc[charger_i].efficiency - ), + stores_t_weights.apply(lambda r: r * n.links.loc[charger_i].efficiency), get_var(n, "Link", "p")[charger_i].T, ) ) discharge_expr = linexpr( ( - n.snapshot_weightings.stores.apply( - lambda r: r * n.links.loc[discharger_i].efficiency - ), + stores_t_weights.apply(lambda r: r * n.links.loc[discharger_i].efficiency), get_var(n, "Link", "p")[discharger_i].T, ) ) From d014acc4836e31417d8b5e0b4a986d5bd738b8e2 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 30 Jun 2024 01:09:40 +0300 Subject: [PATCH 44/57] Replace get_var in add_res_constraints --- scripts/solve_network.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index e104a8e97..a162435c2 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -490,9 +490,7 @@ def add_RES_constraints(n, res_share, config): # Generators lhs_gen = ( - linexpr( - (n.snapshot_weightings.generators, get_var(n, "Generator", "p")[gens_i].T) - ) + linexpr((n.snapshot_weightings.generators, n.model["Generator-p"][gens_i])) .T.groupby(ggrouper, axis=1) .apply(join_exprs) ) @@ -500,25 +498,25 @@ def add_RES_constraints(n, res_share, config): store_disp_expr = linexpr( ( stores_t_weights, - get_var(n, "StorageUnit", "p_dispatch")[stores_i].T, + n.model["StorageUnit-p_dispatch"][stores_i], ) ) store_expr = linexpr( ( stores_t_weights, - get_var(n, "StorageUnit", "p_store")[stores_i].T, + n.model["StorageUnit-p_store"][stores_i], ) ) charge_expr = linexpr( ( stores_t_weights.apply(lambda r: r * n.links.loc[charger_i].efficiency), - get_var(n, "Link", "p")[charger_i].T, + n.model["Link-p"][charger_i], ) ) discharge_expr = linexpr( ( stores_t_weights.apply(lambda r: r * n.links.loc[discharger_i].efficiency), - get_var(n, "Link", "p")[discharger_i].T, + n.model["Link-p"][discharger_i], ) ) From 56ce52e2b43318f09430b6e201b1e128aecab6ed Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 30 Jun 2024 01:12:53 +0300 Subject: [PATCH 45/57] Replace a constraint definition for add_res_constraints --- scripts/solve_network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index a162435c2..39385aecd 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -535,7 +535,7 @@ def form_lhs_definition(lin_expression, ngrouper, target_index=lhs_gen.index): lhs = lhs_gen + lhs_dispatch - lhs_store - lhs_charge + lhs_discharge - define_constraints(n, lhs, "=", rhs, "res_share") + n.model.add_constraints(lhs == rhs, name="res_share") def extra_functionality(n, snapshots): From c1155c3d6b27738cc4b20baac510fcece4c3d0fa Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 30 Jun 2024 01:39:39 +0300 Subject: [PATCH 46/57] Fix typo --- scripts/solve_network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 39385aecd..bcefaf199 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -476,7 +476,7 @@ def add_RES_constraints(n, res_share, config): rhs = res_share * load - renew_techs = config.renewable_carriers + renew_techs = config["electricity"]["renewable_carriers"] charger = ["H2 electrolysis", "battery charger"] discharger = ["H2 fuel cell", "battery discharger"] From a3cbd04b9955d9eb6bd480361ac291053876a6e0 Mon Sep 17 00:00:00 2001 From: ekatef Date: Mon, 1 Jul 2024 11:31:51 +0300 Subject: [PATCH 47/57] Add initialization of the model --- scripts/solve_network.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index bcefaf199..ea4ce1038 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -625,6 +625,8 @@ def solve_network(n, config, solving, **kwargs): solve_opts = snakemake.params.solving["options"] n = pypsa.Network(snakemake.input[0]) + # needed to get `n.model` property + n.optimize.create_model() # TODO Double-check handling the augmented case # if snakemake.params.augmented_line_connection.get("add_to_snakefile"): # n.lines.loc[ From 671dbb53a241b4910f6d386f06c7f5df810da1b6 Mon Sep 17 00:00:00 2001 From: ekatef Date: Mon, 1 Jul 2024 11:32:08 +0300 Subject: [PATCH 48/57] Add a TODO --- scripts/solve_network.py | 1 + 1 file changed, 1 insertion(+) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index ea4ce1038..864a083f5 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -458,6 +458,7 @@ def add_battery_constraints(n): def add_RES_constraints(n, res_share, config): lgrouper = n.loads.bus.map(n.buses.country) + # TODO drop load ggrouper = n.generators.bus.map(n.buses.country) sgrouper = n.storage_units.bus.map(n.buses.country) cgrouper = n.links.bus0.map(n.buses.country) From 78ca7d08c0e6fc07a1bca9c965c0e816edbc7cdc Mon Sep 17 00:00:00 2001 From: ekatef Date: Mon, 1 Jul 2024 15:55:41 +0300 Subject: [PATCH 49/57] Update linear expressions --- scripts/solve_network.py | 34 +++++++++------------------------- 1 file changed, 9 insertions(+), 25 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 864a083f5..584e4cb90 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -491,34 +491,18 @@ def add_RES_constraints(n, res_share, config): # Generators lhs_gen = ( - linexpr((n.snapshot_weightings.generators, n.model["Generator-p"][gens_i])) - .T.groupby(ggrouper, axis=1) - .apply(join_exprs) - ) + n.model["Generator-p"].loc[:, gens_i] * n.snapshot_weightings.generators + ) # .groupby(ggrouper).apply(join_exprs) - store_disp_expr = linexpr( - ( - stores_t_weights, - n.model["StorageUnit-p_dispatch"][stores_i], - ) - ) - store_expr = linexpr( - ( - stores_t_weights, - n.model["StorageUnit-p_store"][stores_i], - ) + store_disp_expr = ( + n.model["StorageUnit-p_dispatch"].loc[:, stores_i] * stores_t_weights ) - charge_expr = linexpr( - ( - stores_t_weights.apply(lambda r: r * n.links.loc[charger_i].efficiency), - n.model["Link-p"][charger_i], - ) + store_expr = n.model["StorageUnit-p_store"].loc[:, stores_i] * stores_t_weights + charge_expr = n.model["Link-p"].loc[:, charger_i] * stores_t_weights.apply( + lambda r: r * n.links.loc[charger_i].efficiency ) - discharge_expr = linexpr( - ( - stores_t_weights.apply(lambda r: r * n.links.loc[discharger_i].efficiency), - n.model["Link-p"][discharger_i], - ) + discharge_expr = n.model["Link-p"].loc[:, discharger_i] * stores_t_weights.apply( + lambda r: r * n.links.loc[discharger_i].efficiency ) def form_lhs_definition(lin_expression, ngrouper, target_index=lhs_gen.index): From 564fd05db29e46450b6b7bffa90807fe5064bcc3 Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 2 Jul 2024 23:29:44 +0300 Subject: [PATCH 50/57] Update constraints --- scripts/solve_network.py | 38 ++++++++++++++++++++++++-------------- 1 file changed, 24 insertions(+), 14 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 584e4cb90..a7fdb247b 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -491,8 +491,10 @@ def add_RES_constraints(n, res_share, config): # Generators lhs_gen = ( - n.model["Generator-p"].loc[:, gens_i] * n.snapshot_weightings.generators - ) # .groupby(ggrouper).apply(join_exprs) + (n.model["Generator-p"].loc[:, gens_i] * n.snapshot_weightings.generators) + # .groupby(ggrouper.to_xarray()) + .sum() + ) store_disp_expr = ( n.model["StorageUnit-p_dispatch"].loc[:, stores_i] * stores_t_weights @@ -505,18 +507,26 @@ def add_RES_constraints(n, res_share, config): lambda r: r * n.links.loc[discharger_i].efficiency ) - def form_lhs_definition(lin_expression, ngrouper, target_index=lhs_gen.index): - lhs_form = ( - (lin_expression.T.groupby(ngrouper, axis=1).apply(join_exprs)) - .reindex(target_index) - .fillna("") - ) - return lhs_form - - lhs_dispatch = form_lhs_definition(store_disp_expr, sgrouper) - lhs_store = form_lhs_definition(store_expr, sgrouper) - lhs_charge = form_lhs_definition(charge_expr, cgrouper) - lhs_discharge = form_lhs_definition(discharge_expr, cgrouper) + lhs_dispatch = ( + store_disp_expr + # .groupby(sgrouper) + .sum() + ) + lhs_store = ( + store_expr + # .groupby(sgrouper) + .sum() + ) + lhs_charge = ( + charge_expr + # .groupby(cgrouper) + .sum() + ) + lhs_discharge = ( + discharge_expr + # .groupby(cgrouper) + .sum() + ) lhs = lhs_gen + lhs_dispatch - lhs_store - lhs_charge + lhs_discharge From 0ea21255771bd487661ee119cf76070f47ca4366 Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 2 Jul 2024 23:32:28 +0300 Subject: [PATCH 51/57] Add TODO --- scripts/solve_network.py | 1 + 1 file changed, 1 insertion(+) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index a7fdb247b..773c45283 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -490,6 +490,7 @@ def add_RES_constraints(n, res_share, config): stores_t_weights = n.snapshot_weightings.stores # Generators + # TODO restore grouping by countries un-commenting calls of groupby() lhs_gen = ( (n.model["Generator-p"].loc[:, gens_i] * n.snapshot_weightings.generators) # .groupby(ggrouper.to_xarray()) From 194ff0ef2c2ee26502d05b80ccebba60df9952a5 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 21 Jul 2024 17:04:14 +0200 Subject: [PATCH 52/57] Add a release note --- doc/release_notes.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/release_notes.rst b/doc/release_notes.rst index acf736f5a..aa7850775 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -24,6 +24,8 @@ E.g. if a new rule becomes available describe how to use it `snakemake -j1 run_t * Add an option to use csv format for custom demand imports. `PR #995 `__ +* Re-implement network optimisation using linopy methods. `PR #796 `__ + **Minor Changes and bug-fixing** * Minor bug-fixing to run the cluster wildcard min `PR #1019 `__ From e90799175951fe589188110022b99f7d277e1009 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 25 Jul 2024 16:32:10 +0200 Subject: [PATCH 53/57] Revise keywords following humanfriendly conventions --- config.default.yaml | 4 ++-- config.tutorial.yaml | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index c404d4484..7b4a953f6 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -464,8 +464,8 @@ solving: cbc-default: {} # Used in CI glpk-default: {} # Used in CI - mem: 30000 #memory in MB; 20 GB enough for 50+B+I+H2; 100 GB for 181+B+I+H2 - walltime: "12:00:00" + mem_mb: 30000 #memory in MB; 20 GB enough for 50+B+I+H2; 100 GB for 181+B+I+H2 + runtime: 6h #runtime in humanfriendly style https://humanfriendly.readthedocs.io/en/latest/ plotting: map: diff --git a/config.tutorial.yaml b/config.tutorial.yaml index 63035ae65..52d9e71b2 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -460,8 +460,8 @@ solving: cbc-default: {} # Used in CI glpk-default: {} # Used in CI - mem: 30000 #memory in MB; 20 GB enough for 50+B+I+H2; 100 GB for 181+B+I+H2 - walltime: "12:00:00" + mem_mb: 30000 #memory in MB; 20 GB enough for 50+B+I+H2; 100 GB for 181+B+I+H2 + runtime: 6h #runtime in humanfriendly style https://humanfriendly.readthedocs.io/en/latest/ plotting: From 6ec28f19b36c1fe2056077a380f4abcfff519888 Mon Sep 17 00:00:00 2001 From: ekatef Date: Mon, 12 Aug 2024 19:49:21 +0200 Subject: [PATCH 54/57] Account for Davide's comment --- scripts/solve_network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 773c45283..2138b6623 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -94,7 +94,7 @@ logger = create_logger(__name__) -def prepare_network(n, solve_opts, config=None): +def prepare_network(n, solve_opts, config): if "clip_p_max_pu" in solve_opts: for df in (n.generators_t.p_max_pu, n.storage_units_t.inflow): df.where(df > solve_opts["clip_p_max_pu"], other=0.0, inplace=True) From f0c95cb95ea14d826e4660bc9b58e2ea4aad325e Mon Sep 17 00:00:00 2001 From: ekatef Date: Mon, 12 Aug 2024 19:54:04 +0200 Subject: [PATCH 55/57] Fix a typo --- scripts/solve_network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 2138b6623..c70e978e6 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -465,7 +465,7 @@ def add_RES_constraints(n, res_share, config): logger.warning( "The add_RES_constraints() is still work in progress. " - "Unexpected results might be incurre, particularly if " + "Unexpected results might be incurred, particularly if " "temporal clustering is applied or if an unexpected change of technologies " "is subject to the optimisation." ) From 4b5d8daa402daff77c4f65444512abdaf5389e3c Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 13 Aug 2024 10:26:51 +0200 Subject: [PATCH 56/57] Add a comment on the approach used for the reserve margin constraint --- scripts/solve_network.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index c70e978e6..1e5ba5a75 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -344,6 +344,14 @@ def add_SAFE_constraints(n, config): def add_operational_reserve_margin_constraint(n, sns, config): + """ + Build reserve margin constraints based on the formulation + as suggested in GenX + https://energy.mit.edu/wp-content/uploads/2017/10/Enhanced-Decision-Support-for-a-Changing-Electricity-Landscape.pdf + It implies that the reserve margin also accounts for optimal + dispatch of distributed energy resources (DERs) and demand response + which is a novel feature of GenX. + """ reserve_config = config["electricity"]["operational_reserve"] EPSILON_LOAD = reserve_config["epsilon_load"] EPSILON_VRES = reserve_config["epsilon_vres"] From 1affa86e67f9f7dd0343a03e266b5c4d07a14599 Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 13 Aug 2024 17:19:44 +0200 Subject: [PATCH 57/57] Switch-off fetching of isolated networks --- config.default.yaml | 2 +- config.tutorial.yaml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index 7b4a953f6..fb1bbffbc 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -65,7 +65,7 @@ cluster_options: remove_stubs_across_borders: true p_threshold_drop_isolated: 20 # [MW] isolated buses are being discarded if bus mean power is below the specified threshold p_threshold_merge_isolated: 300 # [MW] isolated buses are being merged into a single isolated bus if a bus mean power is below the specified threshold - s_threshold_fetch_isolated: 0.05 # [-] a share of the national load for merging an isolated network into a backbone network + s_threshold_fetch_isolated: false # [-] a share of the national load for merging an isolated network into a backbone network cluster_network: algorithm: kmeans feature: solar+onwind-time diff --git a/config.tutorial.yaml b/config.tutorial.yaml index 52d9e71b2..f0960dd72 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -78,7 +78,7 @@ cluster_options: remove_stubs_across_borders: true p_threshold_drop_isolated: 20 # [MW] isolated buses are being discarded if bus mean power is below the specified threshold p_threshold_merge_isolated: 300 # [MW] isolated buses are being merged into a single isolated bus if a bus mean power is below the specified threshold - s_threshold_fetch_isolated: 0.05 # [-] a share of the national load for merging an isolated network into a backbone network + s_threshold_fetch_isolated: false # [-] a share of the national load for merging an isolated network into a backbone network cluster_network: algorithm: kmeans feature: solar+onwind-time