# -*- coding: utf-8 -*- <<<<<<< HEAD # SPDX-FileCopyrightText: : 2017-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`. """ import logging import re from pathlib import Path import numpy as np import pandas as pd import pypsa from _helpers import configure_logging 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 vresutils.benchmark import memory_logger logger = logging.getLogger(__name__) def prepare_network(n, solve_opts): if "clip_p_max_pu" in solve_opts: for df in (n.generators_t.p_max_pu, n.storage_units_t.inflow): df.where(df > solve_opts["clip_p_max_pu"], other=0.0, inplace=True) load_shedding = solve_opts.get("load_shedding") if load_shedding: n.add("Carrier", "load", color="#dd2e23", nice_name="Load shedding") buses_i = n.buses.query("carrier == 'AC'").index if not np.isscalar(load_shedding): 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, carrier="load", sign=1e-3, # Adjust sign to measure p and p_nom in kW instead of MW marginal_cost=load_shedding, ======= """ Solve network. """ import logging import numpy as np import pypsa from helper 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 not "seq" 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 >>>>>>> pypsa-eur-sec/master p_nom=1e9, # kW ) if solve_opts.get("noisy_costs"): <<<<<<< HEAD for t in n.iterate_components(n.one_port_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: ======= 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) >>>>>>> pypsa-eur-sec/master t.df["marginal_cost"] += 1e-2 + 2e-3 * ( np.random.random(len(t.df)) - 0.5 ) for t in n.iterate_components(["Line", "Link"]): <<<<<<< HEAD ======= np.random.seed(123) >>>>>>> pypsa-eur-sec/master 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 <<<<<<< HEAD return n def add_CCL_constraints(n, config): agg_p_nom_limits = 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']." ) logger.info( "Adding per carrier generation capacity constraints for " "individual countries" ) 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, } ) .dropna(subset=["p_nom"]) .groupby(["country", "carrier"]) .p_nom.apply(join_exprs) ) 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): 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) else: ggrouper = n.generators.bus lgrouper = n.loads.bus sgrouper = n.storage_units.bus load = ( n.snapshot_weightings.generators @ n.loads_t.p_set.groupby(lgrouper, axis=1).sum() ) inflow = ( n.snapshot_weightings.stores @ n.storage_units_t.inflow.groupby(sgrouper, axis=1).sum() ) inflow = inflow.reindex(load.index).fillna(0.0) rhs = scaling * (level * load - inflow) lhs_gen = ( linexpr( (n.snapshot_weightings.generators * scaling, get_var(n, "Generator", "p").T) ) .T.groupby(ggrouper, axis=1) .apply(join_exprs) ) if not n.storage_units_t.inflow.empty: lhs_spill = ( linexpr( ( -n.snapshot_weightings.stores * scaling, get_var(n, "StorageUnit", "spill").T, ) ) .T.groupby(sgrouper, axis=1) .apply(join_exprs) ) lhs_spill = lhs_spill.reindex(lhs_gen.index).fillna("") lhs = lhs_gen + lhs_spill else: lhs = lhs_gen define_constraints(n, lhs, ">=", rhs, "equity", "min") def add_BAU_constraints(n, config): mincaps = pd.Series(config["electricity"]["BAU_mincapacities"]) lhs = ( linexpr((1, get_var(n, "Generator", "p_nom"))) .groupby(n.generators.carrier) .apply(join_exprs) ) define_constraints(n, lhs, ">=", mincaps[lhs.index], "Carrier", "bau_mincaps") def add_SAFE_constraints(n, config): peakdemand = ( 1.0 + config["electricity"]["SAFE_reservemargin"] ) * n.loads_t.p_set.sum(axis=1).max() conv_techs = config["plotting"]["conv_techs"] 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") def add_operational_reserve_margin_constraint(n, config): reserve_config = config["electricity"]["operational_reserve"] EPSILON_LOAD = reserve_config["epsilon_load"] EPSILON_VRES = reserve_config["epsilon_vres"] CONTINGENCY = reserve_config["contingency"] # Reserve Variables reserve = get_var(n, "Generator", "r") lhs = linexpr((1, reserve)).sum(1) # 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) # Total demand at t demand = n.loads_t.p_set.sum(1) # VRES potential of non extendable generators capacity_factor = n.generators_t.p_max_pu[vres_i.difference(ext_i)] renewable_capacity = n.generators.p_nom[vres_i.difference(ext_i)] potential = (capacity_factor * renewable_capacity).sum(1) # Right-hand-side rhs = EPSILON_LOAD * demand + EPSILON_VRES * potential + CONTINGENCY define_constraints(n, lhs, ">=", rhs, "Reserve margin") def update_capacity_constraint(n): 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") 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)) 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="" ) rhs = (p_max_pu[fix_i] * capacity_fixed).reindex(columns=gen_i, fill_value=0) define_constraints(n, lhs, "<=", rhs, "Generators", "updated_capacity_constraint") 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. """ define_variables(n, 0, np.inf, "Generator", "r", axes=[sns, n.generators.index]) add_operational_reserve_margin_constraint(n, config) update_capacity_constraint(n) def add_battery_constraints(n): nodes = n.buses.index[n.buses.carrier == "battery"] if nodes.empty or ("Link", "p_nom") not in n.variables.index: 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, ), ) define_constraints(n, lhs, "=", 0, "Link", "charger_ratio") 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) ======= 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) >>>>>>> pypsa-eur-sec/master 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: <<<<<<< HEAD network_lopf( n, solver_name=solver_name, solver_options=solver_options, **kwargs ) else: ilopf( n, solver_name=solver_name, solver_options=solver_options, track_iterations=track_iterations, min_iterations=min_iterations, max_iterations=max_iterations, **kwargs ) return n if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake snakemake = mock_snakemake( "solve_network", simpl="", clusters="5", ll="v1.5", opts="" ) configure_logging(snakemake) tmpdir = snakemake.config["solving"].get("tmpdir") if tmpdir is not None: Path(tmpdir).mkdir(parents=True, exist_ok=True) opts = snakemake.wildcards.opts.split("-") ======= 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 helper import mock_snakemake snakemake = mock_snakemake( "solve_network_myopic", simpl="", opts="", clusters="45", lv=1.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("-") >>>>>>> pypsa-eur-sec/master solve_opts = snakemake.config["solving"]["options"] fn = getattr(snakemake.log, "memory", None) with memory_logger(filename=fn, interval=30.0) as mem: <<<<<<< HEAD n = pypsa.Network(snakemake.input[0]) n = prepare_network(n, solve_opts) n = solve_network( n, snakemake.config, opts, extra_functionality=extra_functionality, solver_dir=tmpdir, solver_logfile=snakemake.log.solver, ) ======= 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"] >>>>>>> pypsa-eur-sec/master 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))