From 4cb21f05ec4c316e00935a1a4d9cb6b905dc795b Mon Sep 17 00:00:00 2001 From: Max Parzen Date: Thu, 9 Mar 2023 13:24:25 +0000 Subject: [PATCH 01/10] add Linopy to PyPSA --- scripts/solve_network.py | 348 ++++++++++++++++++++++++++------------- 1 file changed, 237 insertions(+), 111 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 51891275..9757597f 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -78,6 +78,7 @@ Details (and errors made through this heuristic) are discussed in the paper """ import logging +import os import re from pathlib import Path @@ -85,16 +86,10 @@ import numpy as np 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__) @@ -148,47 +143,90 @@ def prepare_network(n, solve_opts): def add_CCL_constraints(n, config): - agg_p_nom_limits = config["electricity"].get("agg_p_nom_limits") + """ + 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 + """ + pypsa_eur_path = os.path.dirname(os.getcwd()) + agg_p_nom_limits = os.path.join( + pypsa_eur_path, config["electricity"].get("agg_p_nom_limits") + ) try: agg_p_nom_minmax = pd.read_csv(agg_p_nom_limits, index_col=list(range(2))) except IOError: logger.exception( "Need to specify the path to a .csv file containing " - "aggregate capacity limits per country in " - "config['electricity']['agg_p_nom_limit']." + "aggregate capacity limits per country. " + "Path specified in config['electricity']['agg_p_nom_limit']. " + f"Currently read path is 'pypsa-eur/{agg_p_nom_limits}'." ) logger.info( "Adding per carrier generation capacity constraints for " "individual countries" ) + capacity_variable = n.model["Generator-p_nom"] - 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, - } + 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): + """ + Add equality constraints to the network. + + 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 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": @@ -209,78 +247,142 @@ 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") ) if not n.storage_units_t.inflow.empty: + spillage_variable = n.model["StorageUnit-spill"] lhs_spill = ( - linexpr( - ( - -n.snapshot_weightings.stores * scaling, - get_var(n, "StorageUnit", "spill").T, - ) - ) - .T.groupby(sgrouper, axis=1) - .apply(join_exprs) + (spillage_variable * (-n.snapshot_weightings.stores * scaling)) + .groupby_sum(sgrouper) + .sum("snapshot") ) - lhs_spill = lhs_spill.reindex(lhs_gen.index).fillna("") - lhs = lhs_gen + lhs_spill + lhs = merge(lhs_gen, lhs_spill) else: lhs = lhs_gen - define_constraints(n, lhs, ">=", rhs, "equity", "min") + n.model.add_constraints(lhs >= rhs, name="equity_min") 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"]) - 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() + """ + 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 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): + """ + Define minimum operational reserve margin for a given snapshot. + + Parameters + ---------- + n : pypsa.Network + 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 + """ 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 + # Total demand per t demand = n.loads_t.p_set.sum(1) # VRES potential of non extendable generators @@ -291,59 +393,76 @@ 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") + n.model.add_constraints(lhs >= rhs, name="reserve_margin") def update_capacity_constraint(n): + """ + Update the capacity constraint to include the new capacity variables. + + Parameters + ---------- + n : pypsa.Network + """ gen_i = n.generators.index 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") - + 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 = 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): """ 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 """ - 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): + """ + Add constraints to ensure that the ratio between the charger and + discharger. + + 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: 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, + vars_link.sel({"Link-ext": nodes + " discharger"}) * -eff, ) - define_constraints(n, lhs, "=", 0, "Link", "charger_ratio") + n.model.add_constraints(lhs == 0, name="link_charger_ratio") def extra_functionality(n, snapshots): @@ -373,11 +492,8 @@ def extra_functionality(n, snapshots): def solve_network(n, config, opts="", **kwargs): - set_of_options = config["solving"]["solver"]["options"] - solver_options = ( - config["solving"]["solver_options"][set_of_options] if set_of_options else {} - ) - solver_name = config["solving"]["solver"]["name"] + 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) @@ -393,19 +509,25 @@ def solve_network(n, config, opts="", **kwargs): logger.info("No expandable lines found. Skipping iterative solving.") if skip_iterations: - network_lopf( - n, solver_name=solver_name, solver_options=solver_options, **kwargs + optimize( + n, + solver_name=solver_name, + solver_options=solver_options, + extra_functionality=extra_functionality, + **kwargs, ) else: - ilopf( + 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, - **kwargs + extra_functionality=extra_functionality, + **kwargs, ) + return n @@ -413,8 +535,13 @@ if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake + os.chdir(os.path.dirname(os.path.abspath(__file__))) snakemake = mock_snakemake( - "solve_network", simpl="", clusters="5", ll="v1.5", opts="" + "solve_network", + simpl="", + clusters="5", + ll="copt", + opts="Co2L-BAU-24H", ) configure_logging(snakemake) @@ -432,11 +559,10 @@ if __name__ == "__main__": n, snakemake.config, opts, - extra_functionality=extra_functionality, 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)) + logger.info("Maximum memory usage: {}".format(mem.mem_usage)) \ No newline at end of file From a9b09e4ae4fa03725575c6fcf3ab2ad3a81c992d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 9 Mar 2023 13:27:06 +0000 Subject: [PATCH 02/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- 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 9757597f..c7e16de2 100755 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -565,4 +565,4 @@ if __name__ == "__main__": 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)) \ No newline at end of file + logger.info("Maximum memory usage: {}".format(mem.mem_usage)) From 2baf8f5ba41ad9b885d9ef554dfaa33ff048a299 Mon Sep 17 00:00:00 2001 From: Fabian Date: Thu, 9 Mar 2023 22:51:56 +0100 Subject: [PATCH 03/10] harmonive solve_network across workflow --- config.default.yaml | 10 +- rules/solve_electricity.smk | 9 +- rules/solve_myopic.smk | 2 +- rules/solve_overnight.smk | 2 +- scripts/_helpers.py | 7 +- scripts/make_summary.py | 10 +- scripts/solve_network.py | 589 +++++++++++++++------------- scripts/solve_operations_network.py | 94 ++--- scripts/solve_sector_network.py | 365 ----------------- 9 files changed, 371 insertions(+), 717 deletions(-) mode change 100755 => 100644 scripts/solve_network.py delete mode 100644 scripts/solve_sector_network.py diff --git a/config.default.yaml b/config.default.yaml index 23e066b3..a1a26965 100755 --- a/config.default.yaml +++ b/config.default.yaml @@ -629,15 +629,7 @@ solving: track_iterations: false min_iterations: 4 max_iterations: 6 - keep_shadowprices: - - Bus - - Line - - Link - - Transformer - - GlobalConstraint - - Generator - - Store - - StorageUnit + seed: 123 solver: name: gurobi diff --git a/rules/solve_electricity.smk b/rules/solve_electricity.smk index bda5a0a4..8ddeca92 100644 --- a/rules/solve_electricity.smk +++ b/rules/solve_electricity.smk @@ -5,9 +5,9 @@ rule solve_network: input: - RESOURCES + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc", + network=RESOURCES + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc", output: - RESULTS + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc", + network=RESULTS + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc", log: solver=normpath( LOGS + "solve_network/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}_solver.log" @@ -31,10 +31,9 @@ rule solve_network: rule solve_operations_network: input: - unprepared=RESOURCES + "networks/elec_s{simpl}_{clusters}_ec.nc", - optimized=RESULTS + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc", + network=RESULTS + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc", output: - RESULTS + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}_op.nc", + network=RESULTS + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}_op.nc", log: solver=normpath( LOGS diff --git a/rules/solve_myopic.smk b/rules/solve_myopic.smk index b0085d7d..041bee84 100644 --- a/rules/solve_myopic.smk +++ b/rules/solve_myopic.smk @@ -103,4 +103,4 @@ rule solve_sector_network_myopic: conda: "../envs/environment.yaml" script: - "../scripts/solve_sector_network.py" + "../scripts/solve_network.py" diff --git a/rules/solve_overnight.smk b/rules/solve_overnight.smk index 0dc6ad17..c39662ec 100644 --- a/rules/solve_overnight.smk +++ b/rules/solve_overnight.smk @@ -35,4 +35,4 @@ rule solve_sector_network: conda: "../envs/environment.yaml" script: - "../scripts/solve_sector_network.py" + "../scripts/solve_network.py" diff --git a/scripts/_helpers.py b/scripts/_helpers.py index f7f1c557..2e74c75d 100644 --- a/scripts/_helpers.py +++ b/scripts/_helpers.py @@ -339,11 +339,12 @@ def mock_snakemake(rulename, configfiles=[], **wildcards): kwargs = ( dict(rerun_triggers=[]) if parse(sm.__version__) > Version("7.7.0") else {} ) - workflow = sm.Workflow(snakefile, **kwargs) - workflow.include(snakefile) - if isinstance(configfiles, str): configfiles = [configfiles] + + workflow = sm.Workflow(snakefile, overwrite_configfiles=configfiles, **kwargs) + workflow.include(snakefile) + if configfiles: for f in configfiles: if not os.path.exists(f): diff --git a/scripts/make_summary.py b/scripts/make_summary.py index 7aff83fb..e3281de2 100644 --- a/scripts/make_summary.py +++ b/scripts/make_summary.py @@ -442,9 +442,13 @@ def calculate_metrics(n, label, metrics): ["line_volume_AC", "line_volume_DC"], label ].sum() - if hasattr(n, "line_volume_limit"): - metrics.at["line_volume_limit", label] = n.line_volume_limit - metrics.at["line_volume_shadow", label] = n.line_volume_limit_dual + if "lv_limit" in n.global_constraints.index: + metrics.at["line_volume_limit", label] = n.global_constraints.at[ + "lv_limit", "constant" + ] + metrics.at["line_volume_shadow", label] = n.global_constraints.at[ + "lv_limit", "mu" + ] if "CO2Limit" in n.global_constraints.index: metrics.at["co2_shadow", label] = n.global_constraints.at["CO2Limit", "mu"] diff --git a/scripts/solve_network.py b/scripts/solve_network.py old mode 100755 new mode 100644 index c7e16de2..3ad364be --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -1,127 +1,158 @@ # -*- coding: utf-8 -*- -# SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors +# SPDX-FileCopyrightText: : 2020-2023 The PyPSA-Eur Authors # # SPDX-License-Identifier: MIT - """ -Solves linear optimal power flow for a network iteratively while updating -reactances. - -Relevant Settings ------------------ - -.. code:: yaml - - solving: - tmpdir: - options: - formulation: - clip_p_max_pu: - load_shedding: - noisy_costs: - nhours: - min_iterations: - max_iterations: - skip_iterations: - track_iterations: - solver: - name: - -.. seealso:: - Documentation of the configuration file ``config.yaml`` at - :ref:`electricity_cf`, :ref:`solving_cf`, :ref:`plotting_cf` - -Inputs ------- - -- ``networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc``: confer :ref:`prepare` - -Outputs -------- - -- ``results/networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc``: Solved PyPSA network including optimisation results - - .. image:: img/results.png - :scale: 40 % - -Description ------------ - -Total annual system costs are minimised with PyPSA. The full formulation of the -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. - -Solving the network in multiple iterations is motivated through the dependence of transmission line capacities and impedances. -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 - -- 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 `_. - -.. warning:: - Capital costs of existing network components are not included in the objective function, - since for the optimisation problem they are just a constant term (no influence on optimal result). - - Therefore, these capital costs are not included in ``network.objective``! - - If you want to calculate the full total annual system costs add these to the objective value. - -.. tip:: - The rule :mod:`solve_all_networks` runs - for all ``scenario`` s in the configuration file - the rule :mod:`solve_network`. +Solves sector-coupled network. """ import logging -import os import re -from pathlib import Path import numpy as np 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.optimization.abstract import optimize_transmission_expansion_iteratively -from pypsa.optimization.optimize import optimize +import xarray as xr +from _helpers import ( + configure_logging, + override_component_attrs, + update_config_with_sector_opts, +) from vresutils.benchmark import memory_logger logger = logging.getLogger(__name__) +pypsa.pf.logger.setLevel(logging.WARNING) -def prepare_network(n, solve_opts): +def add_land_use_constraint(n, config): + if "m" in snakemake.wildcards.clusters: + _add_land_use_constraint_m(n, config) + else: + _add_land_use_constraint(n, config) + + +def _add_land_use_constraint(n): + # warning: this will miss existing offwind which is not classed AC-DC and has carrier 'offwind' + + for carrier in ["solar", "onwind", "offwind-ac", "offwind-dc"]: + ext_i = (n.generators.carrier == carrier) & ~n.generators.p_nom_extendable + existing = ( + n.generators.loc[ext_i, "p_nom"] + .groupby(n.generators.bus.map(n.buses.location)) + .sum() + ) + existing.index += " " + carrier + "-" + snakemake.wildcards.planning_horizons + n.generators.loc[existing.index, "p_nom_max"] -= existing + + # check if existing capacities are larger than technical potential + existing_large = n.generators[ + n.generators["p_nom_min"] > n.generators["p_nom_max"] + ].index + if len(existing_large): + logger.warning( + f"Existing capacities larger than technical potential for {existing_large},\ + adjust technical potential to existing capacities" + ) + n.generators.loc[existing_large, "p_nom_max"] = n.generators.loc[ + existing_large, "p_nom_min" + ] + + n.generators.p_nom_max.clip(lower=0, inplace=True) + + +def _add_land_use_constraint_m(n, config): + # if generators clustering is lower than network clustering, land_use accounting is at generators clusters + + planning_horizons = config["scenario"]["planning_horizons"] + grouping_years = config["existing_capacities"]["grouping_years"] + current_horizon = snakemake.wildcards.planning_horizons + + for carrier in ["solar", "onwind", "offwind-ac", "offwind-dc"]: + existing = n.generators.loc[n.generators.carrier == carrier, "p_nom"] + ind = list( + set( + [ + i.split(sep=" ")[0] + " " + i.split(sep=" ")[1] + for i in existing.index + ] + ) + ) + + previous_years = [ + str(y) + for y in planning_horizons + grouping_years + if y < int(snakemake.wildcards.planning_horizons) + ] + + for p_year in previous_years: + ind2 = [ + i for i in ind if i + " " + carrier + "-" + p_year in existing.index + ] + sel_current = [i + " " + carrier + "-" + current_horizon for i in ind2] + sel_p_year = [i + " " + carrier + "-" + p_year for i in ind2] + n.generators.loc[sel_current, "p_nom_max"] -= existing.loc[ + sel_p_year + ].rename(lambda x: x[:-4] + current_horizon) + + n.generators.p_nom_max.clip(lower=0, inplace=True) + + +def add_co2_sequestration_limit(n, limit=200): + """ + Add a global constraint on the amount of Mt CO2 that can be sequestered. + """ + n.carriers.loc["co2 stored", "co2_absorptions"] = -1 + n.carriers.co2_absorptions = n.carriers.co2_absorptions.fillna(0) + + limit = limit * 1e6 + for o in opts: + if "seq" not in o: + continue + limit = float(o[o.find("seq") + 3 :]) * 1e6 + break + + n.add( + "GlobalConstraint", + "co2_sequestration_limit", + sense="<=", + constant=limit, + type="primary_energy", + carrier_attribute="co2_absorptions", + ) + + +def prepare_network(n, solve_opts=None, config=None): if "clip_p_max_pu" in solve_opts: - for df in (n.generators_t.p_max_pu, n.storage_units_t.inflow): + for df in ( + n.generators_t.p_max_pu, + n.generators_t.p_min_pu, # TODO: check if this can be removed + n.storage_units_t.inflow, + ): df.where(df > solve_opts["clip_p_max_pu"], other=0.0, inplace=True) - load_shedding = solve_opts.get("load_shedding") - if load_shedding: + if solve_opts.get("load_shedding"): + # intersect between macroeconomic and surveybased willingness to pay + # http://journal.frontiersin.org/article/10.3389/fenrg.2015.00055/full + # TODO: retrieve color and nice name from config n.add("Carrier", "load", color="#dd2e23", nice_name="Load shedding") buses_i = n.buses.query("carrier == 'AC'").index if not np.isscalar(load_shedding): + # TODO: do not scale via sign attribute (use Eur/MWh instead of Eur/kWh) load_shedding = 1e2 # Eur/kWh - # intersect between macroeconomic and surveybased - # willingness to pay - # http://journal.frontiersin.org/article/10.3389/fenrg.2015.00055/full) + n.madd( "Generator", buses_i, " load", - bus=buses_i, + bus=n.buses.index, carrier="load", sign=1e-3, # Adjust sign to measure p and p_nom in kW instead of MW - marginal_cost=load_shedding, + marginal_cost=load_shedding, # Eur/kWh p_nom=1e9, # kW ) if solve_opts.get("noisy_costs"): - for t in n.iterate_components(n.one_port_components): + for t in n.iterate_components(): # if 'capital_cost' in t.df: # t.df['capital_cost'] += 1e1 + 2.*(np.random.random(len(t.df)) - 0.5) if "marginal_cost" in t.df: @@ -139,6 +170,13 @@ def prepare_network(n, solve_opts): n.set_snapshots(n.snapshots[:nhours]) n.snapshot_weightings[:] = 8760.0 / nhours + if config["foresight"] == "myopic": + add_land_use_constraint(n, config) + + if n.stores.carrier.eq("co2 stored").any(): + limit = config["sector"].get("co2_sequestration_potential", 200) + add_co2_sequestration_limit(n, limit=limit) + return n @@ -162,53 +200,37 @@ def add_CCL_constraints(n, config): electricity: agg_p_nom_limits: data/agg_p_nom_minmax.csv """ - pypsa_eur_path = os.path.dirname(os.getcwd()) - agg_p_nom_limits = os.path.join( - pypsa_eur_path, config["electricity"].get("agg_p_nom_limits") + agg_p_nom_minmax = pd.read_csv( + config["electricity"]["agg_p_nom_limits"], index_col=[0, 1] ) - try: - agg_p_nom_minmax = pd.read_csv(agg_p_nom_limits, index_col=list(range(2))) - except IOError: - logger.exception( - "Need to specify the path to a .csv file containing " - "aggregate capacity limits per country. " - "Path specified in config['electricity']['agg_p_nom_limit']. " - f"Currently read path is 'pypsa-eur/{agg_p_nom_limits}'." + logger.info("Adding generation capacity constraints per carrier and country") + p_nom = n.model["Generator-p_nom"] + + gens = n.generators.query("p_nom_extendable").rename_axis(index="Generator-ext") + grouper = [gens.bus.map(n.buses.country), gens.carrier] + grouper = xr.DataArray(pd.MultiIndex.from_arrays(grouper), dims=["Generator-ext"]) + lhs = p_nom.groupby(grouper).sum().rename(bus="country") + + minimum = xr.DataArray(agg_p_nom_minmax["min"].dropna()).rename(dim_0="group") + index = minimum.indexes["group"].intersection(lhs.indexes["group"]) + if not index.empty: + n.model.add_constraints( + lhs.sel(group=index) >= minimum.loc[index], name="agg_p_nom_min" ) - logger.info( - "Adding per carrier generation capacity constraints for " "individual countries" - ) - 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") + maximum = xr.DataArray(agg_p_nom_minmax["max"].dropna()).rename(dim_0="group") + index = maximum.indexes["group"].intersection(lhs.indexes["group"]) + if not index.empty: + n.model.add_constraints( + lhs.sel(group=index) <= maximum.loc[index], name="agg_p_nom_max" ) - 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() - ) def add_EQ_constraints(n, o, scaling=1e-1): """ - Add equality constraints to the network. + Add equity constraints to the network. + + Currently this is only implemented for the electricity sector only. Opts must be specified in the config.yaml. @@ -223,20 +245,21 @@ def add_EQ_constraints(n, o, scaling=1e-1): opts: [Co2L-EQ0.7-24H] Require each country or node to on average produce a minimal share - of its total consumption itself. Example: EQ0.7c demands each country + 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. """ + # TODO: Generalize to cover myopic and other sectors? float_regex = "[0-9]*\.?[0-9]+" level = float(re.findall(float_regex, o)[0]) if o[-1] == "c": - ggrouper = n.generators.bus.map(n.buses.country) - lgrouper = n.loads.bus.map(n.buses.country) - sgrouper = n.storage_units.bus.map(n.buses.country) + ggrouper = n.generators.bus.map(n.buses.country).to_xarray() + lgrouper = n.loads.bus.map(n.buses.country).to_xarray() + sgrouper = n.storage_units.bus.map(n.buses.country).to_xarray() else: - ggrouper = n.generators.bus - lgrouper = n.loads.bus - sgrouper = n.storage_units.bus + ggrouper = n.generators.bus.to_xarray() + lgrouper = n.loads.bus.to_xarray() + sgrouper = n.storage_units.bus.to_xarray() load = ( n.snapshot_weightings.generators @ n.loads_t.p_set.groupby(lgrouper, axis=1).sum() @@ -247,20 +270,23 @@ 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 + p = n.model["Generator-p"] lhs_gen = ( - (dispatch_variable * (n.snapshot_weightings.generators * scaling)) - .groupby_sum(ggrouper) + (p * (n.snapshot_weightings.generators * scaling)) + .groupby(ggrouper) + .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"] + spillage = n.model["StorageUnit-spill"] lhs_spill = ( - (spillage_variable * (-n.snapshot_weightings.stores * scaling)) - .groupby_sum(sgrouper) + (spillage * (-n.snapshot_weightings.stores * scaling)) + .groupby(sgrouper) + .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") @@ -292,14 +318,16 @@ def add_BAU_constraints(n, config): OCGT bus 1 + OCGT bus 2 + ... > 100000 """ 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() + index = mincaps.index.intersection(lhs.indexes["carrier"]) + rhs = mincaps[index].rename_axis("carrier") n.model.add_constraints(lhs >= rhs, name="bau_mincaps") +# TODO: think about removing or make per country def add_SAFE_constraints(n, config): """ Add a capacity reserve margin of a certain fraction above the peak demand. @@ -323,11 +351,11 @@ 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 + # TODO: do not take this from the plotting config! 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() + 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.sum() @@ -335,13 +363,15 @@ def add_SAFE_constraints(n, config): n.model.add_constraints(lhs >= rhs, name="safe_mintotalcap") -def add_operational_reserve_margin_constraint(n, sns, config): +def add_operational_reserve_margin(n, sns, config): """ - Define minimum operational reserve margin for a given snapshot. + 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: @@ -370,130 +400,148 @@ def add_operational_reserve_margin_constraint(n, sns, config): 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 = ( + p_nom_vres = ( 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( - lhs, - (renewable_capacity_variables * (-EPSILON_VRES * capacity_factor)).sum( - ["Generator"] - ), - ) + lhs = lhs + (p_nom_vres * (-EPSILON_VRES * capacity_factor)).sum() # Total demand per t - demand = n.loads_t.p_set.sum(1) + demand = n.loads_t.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 n.model.add_constraints(lhs >= rhs, name="reserve_margin") - -def update_capacity_constraint(n): - """ - Update the capacity constraint to include the new capacity variables. - - Parameters - ---------- - n : pypsa.Network - """ - gen_i = n.generators.index - ext_i = n.generators.query("p_nom_extendable").index - fix_i = n.generators.query("not p_nom_extendable").index - - 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 = merge( - dispatch * 1, - reserve * 1, - ) + lhs = n.model.constraints["Generator-fix-p-upper"].lhs + lhs = lhs + reserve.loc[:, lhs.coords["Generator-fix"]].drop("Generator") + rhs = n.model.constraints["Generator-fix-p-upper"].rhs + n.model.add_constraints(lhs <= rhs, name="Generator-fix-p-upper-reserve") - 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], - ) - - 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): - """ - 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 - """ - add_operational_reserve_margin_constraint(n, sns, config) - update_capacity_constraint(n) + lhs = n.model.constraints["Generator-ext-p-upper"].lhs + lhs = lhs + reserve.loc[:, lhs.coords["Generator-ext"]].drop("Generator") + rhs = n.model.constraints["Generator-ext-p-upper"].rhs + n.model.add_constraints(lhs >= rhs, name="Generator-ext-p-upper-reserve") def add_battery_constraints(n): """ - Add constraints to ensure that the ratio between the charger and - discharger. - + 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: + if not n.links.p_nom_extendable.any(): return - vars_link = n.model["Link-p_nom"] - eff = n.links.loc[nodes + " discharger", "efficiency"] - lhs = merge( - vars_link.sel({"Link-ext": nodes + " charger"}) * 1, - vars_link.sel({"Link-ext": nodes + " discharger"}) * -eff, + + discharger_bool = n.links.index.str.contains("battery discharger") + charger_bool = n.links.index.str.contains("battery charger") + + dischargers_ext = n.links[discharger_bool].query("p_nom_extendable").index + chargers_ext = n.links[charger_bool].query("p_nom_extendable").index + + eff = n.links.efficiency[dischargers_ext].values + lhs = ( + n.model["Link-p_nom"].loc[chargers_ext] + - n.model["Link-p_nom"].loc[dischargers_ext] * eff ) - n.model.add_constraints(lhs == 0, name="link_charger_ratio") + + n.model.add_constraints(lhs == 0, name="Link-charger_ratio") + + +def add_chp_constraints(n): + electric = ( + n.links.index.str.contains("urban central") + & n.links.index.str.contains("CHP") + & n.links.index.str.contains("electric") + ) + heat = ( + n.links.index.str.contains("urban central") + & n.links.index.str.contains("CHP") + & n.links.index.str.contains("heat") + ) + + electric_ext = n.links[electric].query("p_nom_extendable").index + heat_ext = n.links[heat].query("p_nom_extendable").index + + electric_fix = n.links[electric].query("~p_nom_extendable").index + heat_fix = n.links[heat].query("~p_nom_extendable").index + + p = n.model["Link-p"] # dimension: [time, link] + + # output ratio between heat and electricity and top_iso_fuel_line for extendable + if not electric_ext.empty: + p_nom = n.model["Link-p_nom"] + + lhs = ( + p_nom.loc[electric_ext] + * (n.links.p_nom_ratio * n.links.efficiency)[electric_ext].values + - p_nom.loc[heat_ext] * n.links.efficiency[heat_ext].values + ) + n.model.add_constraints(lhs == 0, name="chplink-fix_p_nom_ratio") + + rename = {"Link-ext": "Link"} + lhs = ( + p.loc[:, electric_ext] + + p.loc[:, heat_ext] + - p_nom.rename(rename).loc[electric_ext] + ) + n.model.add_constraints(lhs <= 0, name="chplink-top_iso_fuel_line_ext") + + # top_iso_fuel_line for fixed + if not electric_fix.empty: + lhs = p.loc[:, electric_fix] + p.loc[:, heat_fix] + rhs = n.links.p_nom[electric_fix] + n.model.add_constraints(lhs <= rhs, name="chplink-top_iso_fuel_line_fix") + + # back-pressure + if not electric.empty: + lhs = ( + p.loc[:, heat] * (n.links.efficiency[heat] * n.links.c_b[electric].values) + - p.loc[:, electric] * n.links.efficiency[electric] + ) + n.model.add_constraints(lhs <= rhs, name="chplink-backpressure") + + +def add_pipe_retrofit_constraint(n): + """ + Add constraint for retrofitting existing CH4 pipelines to H2 pipelines. + """ + gas_pipes_i = n.links.query("carrier == 'gas pipeline' and p_nom_extendable").index + h2_retrofitted_i = n.links.query( + "carrier == 'H2 pipeline retrofitted' and p_nom_extendable" + ).index + + if h2_retrofitted_i.empty or gas_pipes_i.empty: + return + + p_nom = n.model["Link-p_nom"] + + CH4_per_H2 = 1 / n.config["sector"]["H2_retrofit_capacity_per_CH4"] + lhs = p_nom.loc[gas_pipes_i] + CH4_per_H2 * p_nom.loc[h2_retrofitted_i] + rhs = n.links.p_nom[gas_pipes_i].rename_axis("Link-ext") + + n.model.add_constraints(lhs == rhs, name="Link-pipe_retrofit") def extra_functionality(n, snapshots): - """ - Collects supplementary constraints which will be passed to - ``pypsa.linopf.network_lopf``. - - 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. - """ - opts = n.opts - config = n.config - if "BAU" in opts and n.generators.p_nom_extendable.any(): - add_BAU_constraints(n, config) - if "SAFE" in opts and n.generators.p_nom_extendable.any(): - add_SAFE_constraints(n, config) - if "CCL" in opts and n.generators.p_nom_extendable.any(): - add_CCL_constraints(n, config) - reserve = config["electricity"].get("operational_reserve", {}) - if reserve.get("activate"): - add_operational_reserve_margin(n, snapshots, config) - for o in opts: - if "EQ" in o: - add_EQ_constraints(n, o) add_battery_constraints(n) + add_pipe_retrofit_constraint(n) def solve_network(n, config, opts="", **kwargs): - solver_options = config["solving"]["solver"].copy() - solver_name = solver_options.pop("name") + set_of_options = config["solving"]["solver"]["options"] + solver_options = ( + config["solving"]["solver_options"][set_of_options] if set_of_options else {} + ) + solver_name = config["solving"]["solver"]["name"] cf_solving = config["solving"]["options"] track_iterations = cf_solving.get("track_iterations", False) min_iterations = cf_solving.get("min_iterations", 4) @@ -509,25 +557,30 @@ def solve_network(n, config, opts="", **kwargs): logger.info("No expandable lines found. Skipping iterative solving.") if skip_iterations: - optimize( - n, + status, condition = n.optimize( solver_name=solver_name, - solver_options=solver_options, extra_functionality=extra_functionality, + **solver_options, **kwargs, ) else: - optimize_transmission_expansion_iteratively( - n, + status, condition = n.optimize.optimize_transmission_expansion_iteratively( solver_name=solver_name, - solver_options=solver_options, track_iterations=track_iterations, min_iterations=min_iterations, max_iterations=max_iterations, extra_functionality=extra_functionality, + **solver_options, **kwargs, ) + if status != "ok": + logger.warning( + f"Solving status '{status}' with termination condition '{condition}'" + ) + if "infeasible" in condition: + raise RuntimeError("Solving status 'infeasible'") + return n @@ -535,33 +588,41 @@ if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake - os.chdir(os.path.dirname(os.path.abspath(__file__))) snakemake = mock_snakemake( - "solve_network", + "solve_sector_network", + configfiles="test/config.overnight.yaml", simpl="", + opts="", clusters="5", - ll="copt", - opts="Co2L-BAU-24H", + ll="v1.5", + sector_opts="CO2L0-24H-T-H-B-I-A-solar+p3-dist1", + planning_horizons="2030", ) configure_logging(snakemake) + update_config_with_sector_opts(snakemake.config, snakemake.wildcards.sector_opts) - tmpdir = snakemake.config["solving"].get("tmpdir") - if tmpdir is not None: - Path(tmpdir).mkdir(parents=True, exist_ok=True) - opts = snakemake.wildcards.opts.split("-") + opts = (snakemake.wildcards.opts + "-" + snakemake.wildcards.sector_opts).split("-") + opts = [o for o in opts if o != ""] solve_opts = snakemake.config["solving"]["options"] + np.random.seed(solve_opts.get("seed", 123)) + fn = getattr(snakemake.log, "memory", None) with memory_logger(filename=fn, interval=30.0) as mem: - n = pypsa.Network(snakemake.input[0]) - n = prepare_network(n, solve_opts) + if "overrides" in snakemake.input.keys(): + overrides = override_component_attrs(snakemake.input.overrides) + n = pypsa.Network( + snakemake.input.network, override_component_attrs=overrides + ) + else: + n = pypsa.Network(snakemake.input.network) + + n = prepare_network(n, solve_opts, config=snakemake.config) + n = solve_network( - n, - snakemake.config, - opts, - solver_dir=tmpdir, - solver_logfile=snakemake.log.solver, + n, config=snakemake.config, opts=opts, log_fn=snakemake.log.solver ) + n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) n.export_to_netcdf(snakemake.output[0]) diff --git a/scripts/solve_operations_network.py b/scripts/solve_operations_network.py index e58a6716..5b16b4a6 100644 --- a/scripts/solve_operations_network.py +++ b/scripts/solve_operations_network.py @@ -46,98 +46,60 @@ Description """ import logging -from pathlib import Path import numpy as np import pypsa -from _helpers import configure_logging +from _helpers import ( + configure_logging, + override_component_attrs, + update_config_with_sector_opts, +) from solve_network import prepare_network, solve_network from vresutils.benchmark import memory_logger logger = logging.getLogger(__name__) -def set_parameters_from_optimized(n, n_optim): - lines_typed_i = n.lines.index[n.lines.type != ""] - n.lines.loc[lines_typed_i, "num_parallel"] = n_optim.lines["num_parallel"].reindex( - lines_typed_i, fill_value=0.0 - ) - n.lines.loc[lines_typed_i, "s_nom"] = ( - np.sqrt(3) - * n.lines["type"].map(n.line_types.i_nom) - * n.lines.bus0.map(n.buses.v_nom) - * n.lines.num_parallel - ) - - lines_untyped_i = n.lines.index[n.lines.type == ""] - for attr in ("s_nom", "r", "x"): - n.lines.loc[lines_untyped_i, attr] = n_optim.lines[attr].reindex( - lines_untyped_i, fill_value=0.0 - ) - n.lines["s_nom_extendable"] = False - - links_dc_i = n.links.index[n.links.p_nom_extendable] - n.links.loc[links_dc_i, "p_nom"] = n_optim.links["p_nom_opt"].reindex( - links_dc_i, fill_value=0.0 - ) - n.links.loc[links_dc_i, "p_nom_extendable"] = False - - gen_extend_i = n.generators.index[n.generators.p_nom_extendable] - n.generators.loc[gen_extend_i, "p_nom"] = n_optim.generators["p_nom_opt"].reindex( - gen_extend_i, fill_value=0.0 - ) - n.generators.loc[gen_extend_i, "p_nom_extendable"] = False - - stor_units_extend_i = n.storage_units.index[n.storage_units.p_nom_extendable] - n.storage_units.loc[stor_units_extend_i, "p_nom"] = n_optim.storage_units[ - "p_nom_opt" - ].reindex(stor_units_extend_i, fill_value=0.0) - n.storage_units.loc[stor_units_extend_i, "p_nom_extendable"] = False - - stor_extend_i = n.stores.index[n.stores.e_nom_extendable] - n.stores.loc[stor_extend_i, "e_nom"] = n_optim.stores["e_nom_opt"].reindex( - stor_extend_i, fill_value=0.0 - ) - n.stores.loc[stor_extend_i, "e_nom_extendable"] = False - - return n - - if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake snakemake = mock_snakemake( "solve_operations_network", + configfiles="test/config.test1.yaml", simpl="", + opts="", clusters="5", - ll="copt", - opts="Co2L-BAU-24H", + ll="v1.5", + sector_opts="", + planning_horizons="", ) + configure_logging(snakemake) + update_config_with_sector_opts(snakemake.config, snakemake.wildcards.sector_opts) - tmpdir = snakemake.config["solving"].get("tmpdir") - if tmpdir is not None: - Path(tmpdir).mkdir(parents=True, exist_ok=True) + opts = (snakemake.wildcards.opts + "-" + snakemake.wildcards.sector_opts).split("-") + opts = [o for o in opts if o != ""] + solve_opts = snakemake.config["solving"]["options"] - n = pypsa.Network(snakemake.input.unprepared) - n_optim = pypsa.Network(snakemake.input.optimized) - n = set_parameters_from_optimized(n, n_optim) - del n_optim - - opts = snakemake.wildcards.opts.split("-") - snakemake.config["solving"]["options"]["skip_iterations"] = False + np.random.seed(solve_opts.get("seed", 123)) fn = getattr(snakemake.log, "memory", None) with memory_logger(filename=fn, interval=30.0) as mem: - n = prepare_network(n, snakemake.config["solving"]["options"]) + if "overrides" in snakemake.input: + overrides = override_component_attrs(snakemake.input.overrides) + n = pypsa.Network( + snakemake.input.network, override_component_attrs=overrides + ) + else: + n = pypsa.Network(snakemake.input.network) + + n.optimize.fix_optimal_capacities() + n = prepare_network(n, solve_opts, config=snakemake.config) n = solve_network( - n, - snakemake.config, - opts, - solver_dir=tmpdir, - solver_logfile=snakemake.log.solver, + n, config=snakemake.config, opts=opts, log_fn=snakemake.log.solver ) + n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) n.export_to_netcdf(snakemake.output[0]) diff --git a/scripts/solve_sector_network.py b/scripts/solve_sector_network.py deleted file mode 100644 index 8bc8b4f6..00000000 --- a/scripts/solve_sector_network.py +++ /dev/null @@ -1,365 +0,0 @@ -# -*- coding: utf-8 -*- -# SPDX-FileCopyrightText: : 2020-2023 The PyPSA-Eur Authors -# -# SPDX-License-Identifier: MIT -""" -Solves sector-coupled network. -""" - -import logging - -import numpy as np -import pypsa -from _helpers import override_component_attrs, update_config_with_sector_opts -from vresutils.benchmark import memory_logger - -logger = logging.getLogger(__name__) -pypsa.pf.logger.setLevel(logging.WARNING) - - -def add_land_use_constraint(n): - if "m" in snakemake.wildcards.clusters: - _add_land_use_constraint_m(n) - else: - _add_land_use_constraint(n) - - -def _add_land_use_constraint(n): - # warning: this will miss existing offwind which is not classed AC-DC and has carrier 'offwind' - - for carrier in ["solar", "onwind", "offwind-ac", "offwind-dc"]: - ext_i = (n.generators.carrier == carrier) & ~n.generators.p_nom_extendable - existing = ( - n.generators.loc[ext_i, "p_nom"] - .groupby(n.generators.bus.map(n.buses.location)) - .sum() - ) - existing.index += " " + carrier + "-" + snakemake.wildcards.planning_horizons - n.generators.loc[existing.index, "p_nom_max"] -= existing - - # check if existing capacities are larger than technical potential - existing_large = n.generators[ - n.generators["p_nom_min"] > n.generators["p_nom_max"] - ].index - if len(existing_large): - logger.warning( - f"Existing capacities larger than technical potential for {existing_large},\ - adjust technical potential to existing capacities" - ) - n.generators.loc[existing_large, "p_nom_max"] = n.generators.loc[ - existing_large, "p_nom_min" - ] - - n.generators.p_nom_max.clip(lower=0, inplace=True) - - -def _add_land_use_constraint_m(n): - # if generators clustering is lower than network clustering, land_use accounting is at generators clusters - - planning_horizons = snakemake.config["scenario"]["planning_horizons"] - grouping_years = snakemake.config["existing_capacities"]["grouping_years"] - current_horizon = snakemake.wildcards.planning_horizons - - for carrier in ["solar", "onwind", "offwind-ac", "offwind-dc"]: - existing = n.generators.loc[n.generators.carrier == carrier, "p_nom"] - ind = list( - set( - [ - i.split(sep=" ")[0] + " " + i.split(sep=" ")[1] - for i in existing.index - ] - ) - ) - - previous_years = [ - str(y) - for y in planning_horizons + grouping_years - if y < int(snakemake.wildcards.planning_horizons) - ] - - for p_year in previous_years: - ind2 = [ - i for i in ind if i + " " + carrier + "-" + p_year in existing.index - ] - sel_current = [i + " " + carrier + "-" + current_horizon for i in ind2] - sel_p_year = [i + " " + carrier + "-" + p_year for i in ind2] - n.generators.loc[sel_current, "p_nom_max"] -= existing.loc[ - sel_p_year - ].rename(lambda x: x[:-4] + current_horizon) - - n.generators.p_nom_max.clip(lower=0, inplace=True) - - -def add_co2_sequestration_limit(n, limit=200): - """ - Add a global constraint on the amount of Mt CO2 that can be sequestered. - """ - n.carriers.loc["co2 stored", "co2_absorptions"] = -1 - n.carriers.co2_absorptions = n.carriers.co2_absorptions.fillna(0) - - limit = limit * 1e6 - for o in opts: - if "seq" not in o: - continue - limit = float(o[o.find("seq") + 3 :]) * 1e6 - break - - n.add( - "GlobalConstraint", - "co2_sequestration_limit", - sense="<=", - constant=limit, - type="primary_energy", - carrier_attribute="co2_absorptions", - ) - - -def prepare_network(n, solve_opts=None, config=None): - if "clip_p_max_pu" in solve_opts: - for df in ( - n.generators_t.p_max_pu, - n.generators_t.p_min_pu, - n.storage_units_t.inflow, - ): - df.where(df > solve_opts["clip_p_max_pu"], other=0.0, inplace=True) - - if solve_opts.get("load_shedding"): - # intersect between macroeconomic and surveybased willingness to pay - # http://journal.frontiersin.org/article/10.3389/fenrg.2015.00055/full - n.add("Carrier", "Load") - n.madd( - "Generator", - n.buses.index, - " load", - bus=n.buses.index, - carrier="load", - sign=1e-3, # Adjust sign to measure p and p_nom in kW instead of MW - marginal_cost=1e2, # Eur/kWh - p_nom=1e9, # kW - ) - - if solve_opts.get("noisy_costs"): - for t in n.iterate_components(): - # if 'capital_cost' in t.df: - # t.df['capital_cost'] += 1e1 + 2.*(np.random.random(len(t.df)) - 0.5) - if "marginal_cost" in t.df: - np.random.seed(174) - t.df["marginal_cost"] += 1e-2 + 2e-3 * ( - np.random.random(len(t.df)) - 0.5 - ) - - for t in n.iterate_components(["Line", "Link"]): - np.random.seed(123) - t.df["capital_cost"] += ( - 1e-1 + 2e-2 * (np.random.random(len(t.df)) - 0.5) - ) * t.df["length"] - - if solve_opts.get("nhours"): - nhours = solve_opts["nhours"] - n.set_snapshots(n.snapshots[:nhours]) - n.snapshot_weightings[:] = 8760.0 / nhours - - if snakemake.config["foresight"] == "myopic": - add_land_use_constraint(n) - - if n.stores.carrier.eq("co2 stored").any(): - limit = config["sector"].get("co2_sequestration_potential", 200) - add_co2_sequestration_limit(n, limit=limit) - - return n - - -def add_battery_constraints(n): - """ - Add constraint ensuring that charger = discharger: - 1 * charger_size - efficiency * discharger_size = 0 - """ - discharger_bool = n.links.index.str.contains("battery discharger") - charger_bool = n.links.index.str.contains("battery charger") - - dischargers_ext = n.links[discharger_bool].query("p_nom_extendable").index - chargers_ext = n.links[charger_bool].query("p_nom_extendable").index - - eff = n.links.efficiency[dischargers_ext].values - lhs = ( - n.model["Link-p_nom"].loc[chargers_ext] - - n.model["Link-p_nom"].loc[dischargers_ext] * eff - ) - - n.model.add_constraints(lhs == 0, name="Link-charger_ratio") - - -def add_chp_constraints(n): - electric = ( - n.links.index.str.contains("urban central") - & n.links.index.str.contains("CHP") - & n.links.index.str.contains("electric") - ) - heat = ( - n.links.index.str.contains("urban central") - & n.links.index.str.contains("CHP") - & n.links.index.str.contains("heat") - ) - - electric_ext = n.links[electric].query("p_nom_extendable").index - heat_ext = n.links[heat].query("p_nom_extendable").index - - electric_fix = n.links[electric].query("~p_nom_extendable").index - heat_fix = n.links[heat].query("~p_nom_extendable").index - - p = n.model["Link-p"] # dimension: [time, link] - - # output ratio between heat and electricity and top_iso_fuel_line for extendable - if not electric_ext.empty: - p_nom = n.model["Link-p_nom"] - - lhs = ( - p_nom.loc[electric_ext] - * (n.links.p_nom_ratio * n.links.efficiency)[electric_ext].values - - p_nom.loc[heat_ext] * n.links.efficiency[heat_ext].values - ) - n.model.add_constraints(lhs == 0, name="chplink-fix_p_nom_ratio") - - rename = {"Link-ext": "Link"} - lhs = ( - p.loc[:, electric_ext] - + p.loc[:, heat_ext] - - p_nom.rename(rename).loc[electric_ext] - ) - n.model.add_constraints(lhs <= 0, name="chplink-top_iso_fuel_line_ext") - - # top_iso_fuel_line for fixed - if not electric_fix.empty: - lhs = p.loc[:, electric_fix] + p.loc[:, heat_fix] - rhs = n.links.p_nom[electric_fix] - n.model.add_constraints(lhs <= rhs, name="chplink-top_iso_fuel_line_fix") - - # back-pressure - if not electric.empty: - lhs = ( - p.loc[:, heat] * (n.links.efficiency[heat] * n.links.c_b[electric].values) - - p.loc[:, electric] * n.links.efficiency[electric] - ) - n.model.add_constraints(lhs <= rhs, name="chplink-backpressure") - - -def add_pipe_retrofit_constraint(n): - """ - Add constraint for retrofitting existing CH4 pipelines to H2 pipelines. - """ - gas_pipes_i = n.links.query("carrier == 'gas pipeline' and p_nom_extendable").index - h2_retrofitted_i = n.links.query( - "carrier == 'H2 pipeline retrofitted' and p_nom_extendable" - ).index - - if h2_retrofitted_i.empty or gas_pipes_i.empty: - return - - p_nom = n.model["Link-p_nom"] - - CH4_per_H2 = 1 / n.config["sector"]["H2_retrofit_capacity_per_CH4"] - lhs = p_nom.loc[gas_pipes_i] + CH4_per_H2 * p_nom.loc[h2_retrofitted_i] - rhs = n.links.p_nom[gas_pipes_i].rename_axis("Link-ext") - - n.model.add_constraints(lhs == rhs, name="Link-pipe_retrofit") - - -def extra_functionality(n, snapshots): - add_battery_constraints(n) - add_pipe_retrofit_constraint(n) - - -def solve_network(n, config, opts="", **kwargs): - set_of_options = config["solving"]["solver"]["options"] - solver_options = ( - config["solving"]["solver_options"][set_of_options] if set_of_options else {} - ) - solver_name = config["solving"]["solver"]["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) - - # 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.") - - if skip_iterations: - status, condition = n.optimize( - solver_name=solver_name, - extra_functionality=extra_functionality, - **solver_options, - **kwargs, - ) - else: - status, condition = n.optimize.optimize_transmission_expansion_iteratively( - solver_name=solver_name, - track_iterations=track_iterations, - min_iterations=min_iterations, - max_iterations=max_iterations, - extra_functionality=extra_functionality, - **solver_options, - **kwargs, - ) - - if status != "ok": - logger.warning( - f"Solving status '{status}' with termination condition '{condition}'" - ) - - return n - - -# %% -if __name__ == "__main__": - if "snakemake" not in globals(): - from _helpers import mock_snakemake - - snakemake = mock_snakemake( - "solve_network_myopic", - simpl="", - opts="", - clusters="45", - ll="v1.0", - sector_opts="8760H-T-H-B-I-A-solar+p3-dist1", - planning_horizons="2020", - ) - - logging.basicConfig( - filename=snakemake.log.python, level=snakemake.config["logging"]["level"] - ) - - update_config_with_sector_opts(snakemake.config, snakemake.wildcards.sector_opts) - - tmpdir = snakemake.config["solving"].get("tmpdir") - if tmpdir is not None: - from pathlib import Path - - Path(tmpdir).mkdir(parents=True, exist_ok=True) - opts = snakemake.wildcards.sector_opts.split("-") - solve_opts = snakemake.config["solving"]["options"] - - fn = getattr(snakemake.log, "memory", None) - with memory_logger(filename=fn, interval=30.0) as mem: - overrides = override_component_attrs(snakemake.input.overrides) - n = pypsa.Network(snakemake.input.network, override_component_attrs=overrides) - - n = prepare_network(n, solve_opts, config=snakemake.config) - - n = solve_network( - n, config=snakemake.config, opts=opts, log_fn=snakemake.log.solver - ) - - if "lv_limit" in n.global_constraints.index: - n.line_volume_limit = n.global_constraints.at["lv_limit", "constant"] - n.line_volume_limit_dual = n.global_constraints.at["lv_limit", "mu"] - - 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)) From 4f0ddf2e95ab5a4b8b52a4564ff72c388a493669 Mon Sep 17 00:00:00 2001 From: Fabian Date: Fri, 10 Mar 2023 14:12:22 +0100 Subject: [PATCH 04/10] fix Nyears scaling --- scripts/build_transport_demand.py | 4 +- scripts/make_summary.py | 2 +- scripts/plot_network.py | 16 +++--- scripts/prepare_sector_network.py | 81 ++++++++++++++++--------------- scripts/solve_network.py | 2 +- test/config.myopic.yaml | 12 +++++ test/config.overnight.yaml | 14 +++++- test/config.test1.yaml | 13 +++++ 8 files changed, 94 insertions(+), 50 deletions(-) diff --git a/scripts/build_transport_demand.py b/scripts/build_transport_demand.py index f6a4b30b..a684036d 100644 --- a/scripts/build_transport_demand.py +++ b/scripts/build_transport_demand.py @@ -83,7 +83,7 @@ def build_transport_demand(traffic_fn, airtemp_fn, nodes, nodal_transport_data): ) transport = ( - (transport_shape.multiply(energy_totals_transport) * 1e6 * Nyears) + (transport_shape.multiply(energy_totals_transport) * 1e6 * nyears) .divide(efficiency_gain * ice_correction) .multiply(1 + dd_EV) ) @@ -181,7 +181,7 @@ if __name__ == "__main__": snapshots = pd.date_range(freq="h", **snakemake.config["snapshots"], tz="UTC") - Nyears = 1 + nyears = len(snapshots) / 8760 nodal_transport_data = build_nodal_transport_data( snakemake.input.transport_data, pop_layout diff --git a/scripts/make_summary.py b/scripts/make_summary.py index e3281de2..5d1ee3dc 100644 --- a/scripts/make_summary.py +++ b/scripts/make_summary.py @@ -699,7 +699,7 @@ if __name__ == "__main__": for planning_horizon in snakemake.config["scenario"]["planning_horizons"] } - Nyears = 1 + Nyears = len(pd.date_range(freq="h", **snakemake.config["snapshots"])) / 8760 costs_db = prepare_costs( snakemake.input.costs, diff --git a/scripts/plot_network.py b/scripts/plot_network.py index 2110718b..76d00c82 100644 --- a/scripts/plot_network.py +++ b/scripts/plot_network.py @@ -23,7 +23,7 @@ from make_summary import assign_carriers from plot_summary import preferred_order, rename_techs from pypsa.plot import add_legend_circles, add_legend_lines, add_legend_patches -plt.style.use(["ggplot", "matplotlibrc"]) +# plt.style.use(["ggplot", "matplotlibrc"]) def rename_techs_tyndp(tech): @@ -394,8 +394,7 @@ def plot_h2_map(network, regions): link_widths=link_widths_retro, branch_components=["Link"], ax=ax, - color_geomap=False, - boundaries=map_opts["boundaries"], + **map_opts, ) regions.plot( @@ -922,11 +921,11 @@ if __name__ == "__main__": snakemake = mock_snakemake( "plot_network", simpl="", - clusters="181", - ll="vopt", opts="", - sector_opts="Co2L0-730H-T-H-B-I-A-solar+p3-linemaxext10", - planning_horizons="2050", + clusters="5", + ll="v1.5", + sector_opts="CO2L0-1H-T-H-B-I-A-solar+p3-dist1", + planning_horizons="2030", ) logging.basicConfig(level=snakemake.config["logging"]["level"]) @@ -938,6 +937,9 @@ if __name__ == "__main__": map_opts = snakemake.config["plotting"]["map"] + if map_opts["boundaries"] is None: + map_opts["boundaries"] = regions.total_bounds[[0, 2, 1, 3]] + [-1, 1, -1, 1] + plot_map( n, components=["generators", "links", "stores", "storage_units"], diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 4f0640cd..470867f2 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -676,7 +676,7 @@ def add_dac(n, costs): ) -def add_co2limit(n, Nyears=1.0, limit=0.0): +def add_co2limit(n, nyears=1.0, limit=0.0): logger.info(f"Adding CO2 budget limit as per unit of 1990 levels of {limit}") countries = snakemake.config["countries"] @@ -688,7 +688,7 @@ def add_co2limit(n, Nyears=1.0, limit=0.0): co2_limit = co2_totals.loc[countries, sectors].sum().sum() - co2_limit *= limit * Nyears + co2_limit *= limit * nyears n.add( "GlobalConstraint", @@ -732,7 +732,7 @@ def cycling_shift(df, steps=1): return df -def prepare_costs(cost_file, config, Nyears): +def prepare_costs(cost_file, config, nyears): # set all asset costs and other parameters costs = pd.read_csv(cost_file, index_col=[0, 1]).sort_index() @@ -750,7 +750,7 @@ def prepare_costs(cost_file, config, Nyears): return annuity(v["lifetime"], v["discount rate"]) + v["FOM"] / 100 costs["fixed"] = [ - annuity_factor(v) * v["investment"] * Nyears for i, v in costs.iterrows() + annuity_factor(v) * v["investment"] * nyears for i, v in costs.iterrows() ] return costs @@ -1557,8 +1557,7 @@ def add_land_transport(n, costs): co2 = ( ice_share / ice_efficiency - * transport[nodes].sum().sum() - / 8760 + * transport[nodes].sum(axis=1) * costs.at["oil", "CO2 intensity"] ) @@ -2333,11 +2332,13 @@ def add_industry(n, costs): logger.info("Add industrial demand") nodes = pop_layout.index + nsnapshots = n.snapshot_weightings.generators.sum() + nyears = nsnapshots / 8760 # 1e6 to convert TWh to MWh industrial_demand = ( pd.read_csv(snakemake.input.industrial_demand, index_col=0) * 1e6 - ) + ) * nyears n.madd( "Bus", @@ -2352,10 +2353,10 @@ def add_industry(n, costs): industrial_demand.loc[spatial.biomass.locations, "solid biomass"].rename( index=lambda x: x + " solid biomass for industry" ) - / 8760 + / nsnapshots ) else: - p_set = industrial_demand["solid biomass"].sum() / 8760 + p_set = industrial_demand["solid biomass"].sum() / nsnapshots n.madd( "Load", @@ -2402,7 +2403,7 @@ def add_industry(n, costs): unit="MWh_LHV", ) - gas_demand = industrial_demand.loc[nodes, "methane"] / 8760.0 + gas_demand = industrial_demand.loc[nodes, "methane"] / nsnapshots if options["gas_network"]: spatial_gas_demand = gas_demand.rename(index=lambda x: x + " gas for industry") @@ -2454,7 +2455,7 @@ def add_industry(n, costs): suffix=" H2 for industry", bus=nodes + " H2", carrier="H2 for industry", - p_set=industrial_demand.loc[nodes, "hydrogen"] / 8760, + p_set=industrial_demand.loc[nodes, "hydrogen"] / nsnapshots, ) shipping_hydrogen_share = get(options["shipping_hydrogen_share"], investment_year) @@ -2470,11 +2471,11 @@ def add_industry(n, costs): domestic_navigation = pop_weighted_energy_totals.loc[ nodes, "total domestic navigation" ].squeeze() - international_navigation = pd.read_csv( - snakemake.input.shipping_demand, index_col=0 - ).squeeze() + international_navigation = ( + pd.read_csv(snakemake.input.shipping_demand, index_col=0).squeeze() * nyears + ) all_navigation = domestic_navigation + international_navigation - p_set = all_navigation * 1e6 / 8760 + p_set = all_navigation * 1e6 / nsnapshots if shipping_hydrogen_share: oil_efficiency = options.get( @@ -2681,7 +2682,7 @@ def add_industry(n, costs): ) demand_factor = options.get("HVC_demand_factor", 1) - p_set = demand_factor * industrial_demand.loc[nodes, "naphtha"].sum() / 8760 + p_set = demand_factor * industrial_demand.loc[nodes, "naphtha"].sum() / nsnapshots if demand_factor != 1: logger.warning(f"Changing HVC demand by {demand_factor*100-100:+.2f}%.") @@ -2699,7 +2700,7 @@ def add_industry(n, costs): demand_factor * pop_weighted_energy_totals.loc[nodes, all_aviation].sum(axis=1).sum() * 1e6 - / 8760 + / nsnapshots ) if demand_factor != 1: logger.warning(f"Changing aviation demand by {demand_factor*100-100:+.2f}%.") @@ -2718,7 +2719,8 @@ def add_industry(n, costs): co2_release = ["naphtha for industry", "kerosene for aviation"] co2 = ( n.loads.loc[co2_release, "p_set"].sum() * costs.at["oil", "CO2 intensity"] - - industrial_demand.loc[nodes, "process emission from feedstock"].sum() / 8760 + - industrial_demand.loc[nodes, "process emission from feedstock"].sum() + / nsnapshots ) n.add( @@ -2741,12 +2743,13 @@ def add_industry(n, costs): for node in nodes ], carrier="low-temperature heat for industry", - p_set=industrial_demand.loc[nodes, "low-temperature heat"] / 8760, + p_set=industrial_demand.loc[nodes, "low-temperature heat"] / nsnapshots, ) # remove today's industrial electricity demand by scaling down total electricity demand for ct in n.buses.country.dropna().unique(): # TODO map onto n.bus.country + loads_i = n.loads.index[ (n.loads.index.str[:2] == ct) & (n.loads.carrier == "electricity") ] @@ -2765,7 +2768,7 @@ def add_industry(n, costs): suffix=" industry electricity", bus=nodes, carrier="industry electricity", - p_set=industrial_demand.loc[nodes, "electricity"] / 8760, + p_set=industrial_demand.loc[nodes, "electricity"] / nsnapshots, ) n.madd( @@ -2782,10 +2785,10 @@ def add_industry(n, costs): -industrial_demand.loc[nodes, sel] .sum(axis=1) .rename(index=lambda x: x + " process emissions") - / 8760 + / nsnapshots ) else: - p_set = -industrial_demand.loc[nodes, sel].sum(axis=1).sum() / 8760 + p_set = -industrial_demand.loc[nodes, sel].sum(axis=1).sum() / nsnapshots # this should be process emissions fossil+feedstock # then need load on atmosphere for feedstock emissions that are currently going to atmosphere via Link Fischer-Tropsch demand @@ -2829,10 +2832,10 @@ def add_industry(n, costs): industrial_demand.loc[spatial.ammonia.locations, "ammonia"].rename( index=lambda x: x + " NH3" ) - / 8760 + / nsnapshots ) else: - p_set = industrial_demand["ammonia"].sum() / 8760 + p_set = industrial_demand["ammonia"].sum() / nsnapshots n.madd( "Load", @@ -2884,6 +2887,7 @@ def add_agriculture(n, costs): logger.info("Add agriculture, forestry and fishing sector.") nodes = pop_layout.index + nsnapshots = n.snapshot_weightings.generators.sum() # electricity @@ -2895,7 +2899,7 @@ def add_agriculture(n, costs): carrier="agriculture electricity", p_set=pop_weighted_energy_totals.loc[nodes, "total agriculture electricity"] * 1e6 - / 8760, + / nsnapshots, ) # heat @@ -2908,7 +2912,7 @@ def add_agriculture(n, costs): carrier="agriculture heat", p_set=pop_weighted_energy_totals.loc[nodes, "total agriculture heat"] * 1e6 - / 8760, + / nsnapshots, ) # machinery @@ -2944,7 +2948,7 @@ def add_agriculture(n, costs): / efficiency_gain * machinery_nodal_energy * 1e6 - / 8760, + / nsnapshots, ) if oil_share > 0: @@ -2953,14 +2957,14 @@ def add_agriculture(n, costs): ["agriculture machinery oil"], bus=spatial.oil.nodes, carrier="agriculture machinery oil", - p_set=oil_share * machinery_nodal_energy.sum() * 1e6 / 8760, + p_set=oil_share * machinery_nodal_energy.sum() * 1e6 / nsnapshots, ) co2 = ( oil_share * machinery_nodal_energy.sum() * 1e6 - / 8760 + / nsnapshots * costs.at["oil", "CO2 intensity"] ) @@ -3229,19 +3233,19 @@ def set_temporal_aggregation(n, opts, solver_name): return n -# %% if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake snakemake = mock_snakemake( "prepare_sector_network", + configfiles="test/config.overnight.yaml", simpl="", opts="", - clusters="37", + clusters="5", ll="v1.5", - sector_opts="cb40ex0-365H-T-H-B-I-A-solar+p3-dist1", - planning_horizons="2020", + sector_opts="CO2L0-24H-T-H-B-I-A-solar+p3-dist1", + planning_horizons="2030", ) logging.basicConfig(level=snakemake.config["logging"]["level"]) @@ -3258,16 +3262,17 @@ if __name__ == "__main__": n = pypsa.Network(snakemake.input.network, override_component_attrs=overrides) pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0) - Nyears = n.snapshot_weightings.generators.sum() / 8760 + nsnapshots = n.snapshot_weightings.objective.sum() + nyears = nsnapshots / 8760 costs = prepare_costs( snakemake.input.costs, snakemake.config["costs"], - Nyears, + nyears, ) - pop_weighted_energy_totals = pd.read_csv( - snakemake.input.pop_weighted_energy_totals, index_col=0 + pop_weighted_energy_totals = ( + pd.read_csv(snakemake.input.pop_weighted_energy_totals, index_col=0) * nyears ) patch_electricity_network(n) @@ -3369,7 +3374,7 @@ if __name__ == "__main__": limit = float(limit.replace("p", ".").replace("m", "-")) break logger.info(f"Add CO2 limit from {limit_type}") - add_co2limit(n, Nyears, limit) + add_co2limit(n, nyears, limit) for o in opts: if not o[:10] == "linemaxext": diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 3ad364be..99458f7e 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -31,7 +31,7 @@ def add_land_use_constraint(n, config): _add_land_use_constraint(n, config) -def _add_land_use_constraint(n): +def _add_land_use_constraint(n, config): # warning: this will miss existing offwind which is not classed AC-DC and has carrier 'offwind' for carrier in ["solar", "onwind", "offwind-ac", "offwind-dc"]: diff --git a/test/config.myopic.yaml b/test/config.myopic.yaml index 71c05f3c..9136d20d 100644 --- a/test/config.myopic.yaml +++ b/test/config.myopic.yaml @@ -58,3 +58,15 @@ solving: name: glpk options: glpk-default mem: 4000 + +plotting: + map: + boundaries: + eu_node_location: + x: -5.5 + y: 46. + costs_max: 1000 + costs_threshold: 0.0000001 + energy_max: + energy_min: + energy_threshold: 0.000001 diff --git a/test/config.overnight.yaml b/test/config.overnight.yaml index ec9bd502..318dc52a 100644 --- a/test/config.overnight.yaml +++ b/test/config.overnight.yaml @@ -16,7 +16,7 @@ scenario: clusters: - 5 sector_opts: - - CO2L0-24H-T-H-B-I-A-solar+p3-dist1 + - CO2L0-1H-T-H-B-I-A-solar+p3-dist1 planning_horizons: - 2030 @@ -59,3 +59,15 @@ solving: name: glpk options: glpk-default mem: 4000 + +plotting: + map: + boundaries: + eu_node_location: + x: -5.5 + y: 46. + costs_max: 1000 + costs_threshold: 0.0000001 + energy_max: + energy_min: + energy_threshold: 0.000001 diff --git a/test/config.test1.yaml b/test/config.test1.yaml index d95858de..6798e38c 100755 --- a/test/config.test1.yaml +++ b/test/config.test1.yaml @@ -65,3 +65,16 @@ solving: solver: name: glpk options: "glpk-default" + + +plotting: + map: + boundaries: + eu_node_location: + x: -5.5 + y: 46. + costs_max: 1000 + costs_threshold: 0.0000001 + energy_max: + energy_min: + energy_threshold: 0.000001 From f8d178c7147a22d8b10ae3e43aefa3d58703998a Mon Sep 17 00:00:00 2001 From: Fabian Date: Fri, 10 Mar 2023 14:26:31 +0100 Subject: [PATCH 05/10] prepare_sector: revert taking mean config.overnight: revert hourly resolution --- scripts/prepare_sector_network.py | 6 ++++-- test/config.overnight.yaml | 2 +- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 470867f2..f90f608f 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -1409,6 +1409,7 @@ def add_land_transport(n, costs): # TODO options? logger.info("Add land transport") + nsnapshots = n.snapshot_weightings.generators.sum() transport = pd.read_csv( snakemake.input.transport_demand, index_col=0, parse_dates=True @@ -1557,7 +1558,8 @@ def add_land_transport(n, costs): co2 = ( ice_share / ice_efficiency - * transport[nodes].sum(axis=1) + * transport[nodes].sum().sum() + / nsnapshots * costs.at["oil", "CO2 intensity"] ) @@ -3262,7 +3264,7 @@ if __name__ == "__main__": n = pypsa.Network(snakemake.input.network, override_component_attrs=overrides) pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0) - nsnapshots = n.snapshot_weightings.objective.sum() + nsnapshots = n.snapshot_weightings.generators.sum() nyears = nsnapshots / 8760 costs = prepare_costs( diff --git a/test/config.overnight.yaml b/test/config.overnight.yaml index 318dc52a..06cc5fd6 100644 --- a/test/config.overnight.yaml +++ b/test/config.overnight.yaml @@ -16,7 +16,7 @@ scenario: clusters: - 5 sector_opts: - - CO2L0-1H-T-H-B-I-A-solar+p3-dist1 + - CO2L0-24H-T-H-B-I-A-solar+p3-dist1 planning_horizons: - 2030 From 99a9ea8beb8088332d9098f778a466b98179e623 Mon Sep 17 00:00:00 2001 From: Fabian Date: Fri, 10 Mar 2023 15:58:53 +0100 Subject: [PATCH 06/10] adress review comments --- scripts/plot_network.py | 2 +- scripts/prepare_sector_network.py | 53 +++++++++++++++---------------- scripts/solve_network.py | 22 +++++++++++++ 3 files changed, 49 insertions(+), 28 deletions(-) diff --git a/scripts/plot_network.py b/scripts/plot_network.py index 76d00c82..f5063ce8 100644 --- a/scripts/plot_network.py +++ b/scripts/plot_network.py @@ -23,7 +23,7 @@ from make_summary import assign_carriers from plot_summary import preferred_order, rename_techs from pypsa.plot import add_legend_circles, add_legend_lines, add_legend_patches -# plt.style.use(["ggplot", "matplotlibrc"]) +plt.style.use(["ggplot", "matplotlibrc"]) def rename_techs_tyndp(tech): diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index f90f608f..da6eab72 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -1409,7 +1409,7 @@ def add_land_transport(n, costs): # TODO options? logger.info("Add land transport") - nsnapshots = n.snapshot_weightings.generators.sum() + nhours = n.snapshot_weightings.generators.sum() transport = pd.read_csv( snakemake.input.transport_demand, index_col=0, parse_dates=True @@ -1559,7 +1559,7 @@ def add_land_transport(n, costs): ice_share / ice_efficiency * transport[nodes].sum().sum() - / nsnapshots + / nhours * costs.at["oil", "CO2 intensity"] ) @@ -2334,8 +2334,8 @@ def add_industry(n, costs): logger.info("Add industrial demand") nodes = pop_layout.index - nsnapshots = n.snapshot_weightings.generators.sum() - nyears = nsnapshots / 8760 + nhours = n.snapshot_weightings.generators.sum() + nyears = nhours / 8760 # 1e6 to convert TWh to MWh industrial_demand = ( @@ -2355,10 +2355,10 @@ def add_industry(n, costs): industrial_demand.loc[spatial.biomass.locations, "solid biomass"].rename( index=lambda x: x + " solid biomass for industry" ) - / nsnapshots + / nhours ) else: - p_set = industrial_demand["solid biomass"].sum() / nsnapshots + p_set = industrial_demand["solid biomass"].sum() / nhours n.madd( "Load", @@ -2405,7 +2405,7 @@ def add_industry(n, costs): unit="MWh_LHV", ) - gas_demand = industrial_demand.loc[nodes, "methane"] / nsnapshots + gas_demand = industrial_demand.loc[nodes, "methane"] / nhours if options["gas_network"]: spatial_gas_demand = gas_demand.rename(index=lambda x: x + " gas for industry") @@ -2457,7 +2457,7 @@ def add_industry(n, costs): suffix=" H2 for industry", bus=nodes + " H2", carrier="H2 for industry", - p_set=industrial_demand.loc[nodes, "hydrogen"] / nsnapshots, + p_set=industrial_demand.loc[nodes, "hydrogen"] / nhours, ) shipping_hydrogen_share = get(options["shipping_hydrogen_share"], investment_year) @@ -2477,7 +2477,7 @@ def add_industry(n, costs): pd.read_csv(snakemake.input.shipping_demand, index_col=0).squeeze() * nyears ) all_navigation = domestic_navigation + international_navigation - p_set = all_navigation * 1e6 / nsnapshots + p_set = all_navigation * 1e6 / nhours if shipping_hydrogen_share: oil_efficiency = options.get( @@ -2684,7 +2684,7 @@ def add_industry(n, costs): ) demand_factor = options.get("HVC_demand_factor", 1) - p_set = demand_factor * industrial_demand.loc[nodes, "naphtha"].sum() / nsnapshots + p_set = demand_factor * industrial_demand.loc[nodes, "naphtha"].sum() / nhours if demand_factor != 1: logger.warning(f"Changing HVC demand by {demand_factor*100-100:+.2f}%.") @@ -2702,7 +2702,7 @@ def add_industry(n, costs): demand_factor * pop_weighted_energy_totals.loc[nodes, all_aviation].sum(axis=1).sum() * 1e6 - / nsnapshots + / nhours ) if demand_factor != 1: logger.warning(f"Changing aviation demand by {demand_factor*100-100:+.2f}%.") @@ -2721,8 +2721,7 @@ def add_industry(n, costs): co2_release = ["naphtha for industry", "kerosene for aviation"] co2 = ( n.loads.loc[co2_release, "p_set"].sum() * costs.at["oil", "CO2 intensity"] - - industrial_demand.loc[nodes, "process emission from feedstock"].sum() - / nsnapshots + - industrial_demand.loc[nodes, "process emission from feedstock"].sum() / nhours ) n.add( @@ -2745,7 +2744,7 @@ def add_industry(n, costs): for node in nodes ], carrier="low-temperature heat for industry", - p_set=industrial_demand.loc[nodes, "low-temperature heat"] / nsnapshots, + p_set=industrial_demand.loc[nodes, "low-temperature heat"] / nhours, ) # remove today's industrial electricity demand by scaling down total electricity demand @@ -2770,7 +2769,7 @@ def add_industry(n, costs): suffix=" industry electricity", bus=nodes, carrier="industry electricity", - p_set=industrial_demand.loc[nodes, "electricity"] / nsnapshots, + p_set=industrial_demand.loc[nodes, "electricity"] / nhours, ) n.madd( @@ -2787,10 +2786,10 @@ def add_industry(n, costs): -industrial_demand.loc[nodes, sel] .sum(axis=1) .rename(index=lambda x: x + " process emissions") - / nsnapshots + / nhours ) else: - p_set = -industrial_demand.loc[nodes, sel].sum(axis=1).sum() / nsnapshots + p_set = -industrial_demand.loc[nodes, sel].sum(axis=1).sum() / nhours # this should be process emissions fossil+feedstock # then need load on atmosphere for feedstock emissions that are currently going to atmosphere via Link Fischer-Tropsch demand @@ -2834,10 +2833,10 @@ def add_industry(n, costs): industrial_demand.loc[spatial.ammonia.locations, "ammonia"].rename( index=lambda x: x + " NH3" ) - / nsnapshots + / nhours ) else: - p_set = industrial_demand["ammonia"].sum() / nsnapshots + p_set = industrial_demand["ammonia"].sum() / nhours n.madd( "Load", @@ -2889,7 +2888,7 @@ def add_agriculture(n, costs): logger.info("Add agriculture, forestry and fishing sector.") nodes = pop_layout.index - nsnapshots = n.snapshot_weightings.generators.sum() + nhours = n.snapshot_weightings.generators.sum() # electricity @@ -2901,7 +2900,7 @@ def add_agriculture(n, costs): carrier="agriculture electricity", p_set=pop_weighted_energy_totals.loc[nodes, "total agriculture electricity"] * 1e6 - / nsnapshots, + / nhours, ) # heat @@ -2914,7 +2913,7 @@ def add_agriculture(n, costs): carrier="agriculture heat", p_set=pop_weighted_energy_totals.loc[nodes, "total agriculture heat"] * 1e6 - / nsnapshots, + / nhours, ) # machinery @@ -2950,7 +2949,7 @@ def add_agriculture(n, costs): / efficiency_gain * machinery_nodal_energy * 1e6 - / nsnapshots, + / nhours, ) if oil_share > 0: @@ -2959,14 +2958,14 @@ def add_agriculture(n, costs): ["agriculture machinery oil"], bus=spatial.oil.nodes, carrier="agriculture machinery oil", - p_set=oil_share * machinery_nodal_energy.sum() * 1e6 / nsnapshots, + p_set=oil_share * machinery_nodal_energy.sum() * 1e6 / nhours, ) co2 = ( oil_share * machinery_nodal_energy.sum() * 1e6 - / nsnapshots + / nhours * costs.at["oil", "CO2 intensity"] ) @@ -3264,8 +3263,8 @@ if __name__ == "__main__": n = pypsa.Network(snakemake.input.network, override_component_attrs=overrides) pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0) - nsnapshots = n.snapshot_weightings.generators.sum() - nyears = nsnapshots / 8760 + nhours = n.snapshot_weightings.generators.sum() + nyears = nhours / 8760 costs = prepare_costs( snakemake.input.costs, diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 99458f7e..1fa27b32 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -532,6 +532,28 @@ def add_pipe_retrofit_constraint(n): def extra_functionality(n, snapshots): + """ + Collects supplementary constraints which will be passed to + ``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. + """ + opts = n.opts + config = n.config + if "BAU" in opts and n.generators.p_nom_extendable.any(): + add_BAU_constraints(n, config) + if "SAFE" in opts and n.generators.p_nom_extendable.any(): + add_SAFE_constraints(n, config) + if "CCL" in opts and n.generators.p_nom_extendable.any(): + add_CCL_constraints(n, config) + reserve = config["electricity"].get("operational_reserve", {}) + if reserve.get("activate"): + add_operational_reserve_margin(n, snapshots, config) + for o in opts: + if "EQ" in o: + add_EQ_constraints(n, o) add_battery_constraints(n) add_pipe_retrofit_constraint(n) From 30b8bb3c941e83e29ef7c6c0e2f825a68c966b41 Mon Sep 17 00:00:00 2001 From: Fabian Hofmann Date: Fri, 10 Mar 2023 16:00:52 +0100 Subject: [PATCH 07/10] Update scripts/solve_network.py Co-authored-by: Fabian Neumann --- 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 1fa27b32..a36e9edf 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -# SPDX-FileCopyrightText: : 2020-2023 The PyPSA-Eur Authors +# SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors # # SPDX-License-Identifier: MIT """ From 5fe4032e473715a46738a3175ee8a30a82053846 Mon Sep 17 00:00:00 2001 From: Fabian Date: Fri, 10 Mar 2023 16:07:00 +0100 Subject: [PATCH 08/10] solve_network: reinsert docstring --- scripts/solve_network.py | 67 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 65 insertions(+), 2 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index a36e9edf..e22b1a6d 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -3,9 +3,72 @@ # # SPDX-License-Identifier: MIT """ -Solves sector-coupled network. -""" +Solves optimal operation and capacity for a network with the option to +iteratively optimize while updating line reactances. +This script is used for optimizing the electrical network as well as the +sector coupled network. + +Relevant Settings +----------------- +.. code:: yaml + solving: + options: + formulation: + clip_p_max_pu: + load_shedding: + noisy_costs: + nhours: + min_iterations: + max_iterations: + skip_iterations: + track_iterations: + solver: + name: + options: + +.. seealso:: + Documentation of the configuration file ``config.yaml`` at + :ref:`electricity_cf`, :ref:`solving_cf`, :ref:`plotting_cf` + +Inputs +------ +- ``networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc``: confer :ref:`prepare` + +Outputs +------- +- ``results/networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc``: Solved PyPSA network including optimisation results + .. image:: ../img/results.png + :scale: 40 % + +Description +----------- +Total annual system costs are minimised with PyPSA. The full formulation of the +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. +Solving the network in multiple iterations is motivated through the dependence of transmission line capacities and impedances. +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 + +- 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 `_. + +.. warning:: + Capital costs of existing network components are not included in the objective function, + since for the optimisation problem they are just a constant term (no influence on optimal result). + Therefore, these capital costs are not included in ``network.objective``! + If you want to calculate the full total annual system costs add these to the objective value. + +.. tip:: + The rule :mod:`solve_all_networks` runs + for all ``scenario`` s in the configuration file + the rule :mod:`solve_network`. +""" import logging import re From e9bf639f2d20ea20c28960695f2c30895adef7d0 Mon Sep 17 00:00:00 2001 From: Fabian Date: Fri, 10 Mar 2023 16:10:01 +0100 Subject: [PATCH 09/10] solve_network: add line breaks in docstring --- scripts/solve_network.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index e22b1a6d..f670ab5d 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -2,6 +2,7 @@ # SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors # # SPDX-License-Identifier: MIT + """ Solves optimal operation and capacity for a network with the option to iteratively optimize while updating line reactances. @@ -11,7 +12,9 @@ sector coupled network. Relevant Settings ----------------- + .. code:: yaml + solving: options: formulation: @@ -33,22 +36,28 @@ Relevant Settings Inputs ------ + - ``networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc``: confer :ref:`prepare` Outputs ------- + - ``results/networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc``: Solved PyPSA network including optimisation results + .. image:: ../img/results.png :scale: 40 % Description ----------- + Total annual system costs are minimised with PyPSA. The full formulation of the 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. + Solving the network in multiple iterations is motivated through the dependence of transmission line capacities and impedances. As lines are expanded their electrical parameters change, which renders the optimisation bilinear even if the power flow equations are linearized. @@ -59,12 +68,15 @@ Details (and errors made 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 `_. .. warning:: + Capital costs of existing network components are not included in the objective function, since for the optimisation problem they are just a constant term (no influence on optimal result). + Therefore, these capital costs are not included in ``network.objective``! If you want to calculate the full total annual system costs add these to the objective value. .. tip:: + The rule :mod:`solve_all_networks` runs for all ``scenario`` s in the configuration file the rule :mod:`solve_network`. From ebc5a993388c31f170a7a56002d0e2f478d48506 Mon Sep 17 00:00:00 2001 From: Fabian Date: Fri, 10 Mar 2023 16:12:02 +0100 Subject: [PATCH 10/10] follow up: fix img path --- 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 f670ab5d..ec27d04a 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -44,7 +44,7 @@ Outputs - ``results/networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc``: Solved PyPSA network including optimisation results - .. image:: ../img/results.png + .. image:: img/results.png :scale: 40 % Description