harmonive solve_network across workflow

config.default.yaml

@@ -629,15 +629,7 @@ solving:
     track_iterations: false
     min_iterations: 4
     max_iterations: 6
-    keep_shadowprices:
+    seed: 123
-    - Bus
-    - Line
-    - Link
-    - Transformer
-    - GlobalConstraint
-    - Generator
-    - Store
-    - StorageUnit
 
   solver:
     name: gurobi

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",
+        network=RESULTS + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc",
-        optimized=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

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"

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"

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):

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"):
+    if "lv_limit" in n.global_constraints.index:
-        metrics.at["line_volume_limit", label] = n.line_volume_limit
+        metrics.at["line_volume_limit", label] = n.global_constraints.at[
-        metrics.at["line_volume_shadow", label] = n.line_volume_limit_dual
+            "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"]

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
+Solves sector-coupled network.
-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 <https://pypsa.readthedocs.io/en/latest/optimal_power_flow.html#linear-optimal-power-flow>`_.
-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 <https://arxiv.org/abs/1907.10548>`_), *16th International Conference on the European Energy Market*, 2019. `arXiv:1907.10548 <https://arxiv.org/abs/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 os
 import re
-from pathlib import Path
 
 import numpy as np
 import pandas as pd
 import pypsa
-from _helpers import configure_logging
+import xarray as xr
-from linopy import merge
+from _helpers import (
-from pypsa.descriptors import get_switchable_as_dense as get_as_dense
+    configure_logging,
-from pypsa.optimization.abstract import optimize_transmission_expansion_iteratively
+    override_component_attrs,
-from pypsa.optimization.optimize import optimize
+    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 solve_opts.get("load_shedding"):
-    if 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_minmax = pd.read_csv(
-    agg_p_nom_limits = os.path.join(
+        config["electricity"]["agg_p_nom_limits"], index_col=[0, 1]
-        pypsa_eur_path, config["electricity"].get("agg_p_nom_limits")
     )
-    try:
+    logger.info("Adding generation capacity constraints per carrier and country")
-        agg_p_nom_minmax = pd.read_csv(agg_p_nom_limits, index_col=list(range(2)))
+    p_nom = n.model["Generator-p_nom"]
-    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 per carrier generation capacity constraints for " "individual countries"
-    )
-    capacity_variable = n.model["Generator-p_nom"]
 
-    lhs = []
+    gens = n.generators.query("p_nom_extendable").rename_axis(index="Generator-ext")
-    ext_carriers = n.generators.query("p_nom_extendable").carrier.unique()
+    grouper = [gens.bus.map(n.buses.country), gens.carrier]
-    for c in ext_carriers:
+    grouper = xr.DataArray(pd.MultiIndex.from_arrays(grouper), dims=["Generator-ext"])
-        ext_carrier = n.generators.query("p_nom_extendable and carrier == @c")
+    lhs = p_nom.groupby(grouper).sum().rename(bus="country")
-        country_grouper = (
-            ext_carrier.bus.map(n.buses.country)
-            .rename_axis("Generator-ext")
-            .rename("country")
-        )
-        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)
 
+    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 >= min_matrix, name="agg_p_nom_min", mask=min_matrix.notnull()
+            lhs.sel(group=index) >= minimum.loc[index], name="agg_p_nom_min"
         )
+
+    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 <= max_matrix, name="agg_p_nom_max", mask=max_matrix.notnull()
+            lhs.sel(group=index) <= maximum.loc[index], name="agg_p_nom_max"
         )
 
 
 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)
+        ggrouper = n.generators.bus.map(n.buses.country).to_xarray()
-        lgrouper = n.loads.bus.map(n.buses.country)
+        lgrouper = n.loads.bus.map(n.buses.country).to_xarray()
-        sgrouper = n.storage_units.bus.map(n.buses.country)
+        sgrouper = n.storage_units.bus.map(n.buses.country).to_xarray()
     else:
-        ggrouper = n.generators.bus
+        ggrouper = n.generators.bus.to_xarray()
-        lgrouper = n.loads.bus
+        lgrouper = n.loads.bus.to_xarray()
-        sgrouper = n.storage_units.bus
+        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))
+        (p * (n.snapshot_weightings.generators * scaling))
-        .groupby_sum(ggrouper)
+        .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))
+            (spillage * (-n.snapshot_weightings.stores * scaling))
-            .groupby_sum(sgrouper)
+            .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")
+    ext_carrier_i = xr.DataArray(ext_i.carrier.rename_axis("Generator-ext"))
-    lhs = capacity_variable.groupby_sum(ext_carrier_i)
+    lhs = p_nom.groupby(ext_carrier_i).sum()
-    rhs = mincaps[lhs.coords["carrier"].values].rename_axis("carrier")
+    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"]
+    p_nom = n.model["Generator-p_nom"].loc[ext_gens_i]
-    ext_cap_var = capacity_variable.sel({"Generator-ext": ext_gens_i})
+    lhs = p_nom.sum()
-    lhs = ext_cap_var.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 = lhs + (p_nom_vres * (-EPSILON_VRES * capacity_factor)).sum()
-            lhs,
-            (renewable_capacity_variables * (-EPSILON_VRES * capacity_factor)).sum(
-                ["Generator"]
-            ),
-        )
 
     # 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(
+    lhs = n.model.constraints["Generator-fix-p-upper"].lhs
-        dispatch * 1,
+    lhs = lhs + reserve.loc[:, lhs.coords["Generator-fix"]].drop("Generator")
-        reserve * 1,
+    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:
+    lhs = n.model.constraints["Generator-ext-p-upper"].lhs
-        capacity_variable = n.model["Generator-p_nom"]
+    lhs = lhs + reserve.loc[:, lhs.coords["Generator-ext"]].drop("Generator")
-        lhs = merge(
+    rhs = n.model.constraints["Generator-ext-p-upper"].rhs
-            lhs,
+    n.model.add_constraints(lhs >= rhs, name="Generator-ext-p-upper-reserve")
-            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)
 
 
 def add_battery_constraints(n):
     """
-    Add constraints to ensure that the ratio between the charger and
+    Add constraint ensuring that charger = discharger, i.e.
-    discharger.
-
     1 * charger_size - efficiency * discharger_size = 0
     """
-    nodes = n.buses.index[n.buses.carrier == "battery"]
+    if not n.links.p_nom_extendable.any():
-    if nodes.empty:
         return
-    vars_link = n.model["Link-p_nom"]
+
-    eff = n.links.loc[nodes + " discharger", "efficiency"]
+    discharger_bool = n.links.index.str.contains("battery discharger")
-    lhs = merge(
+    charger_bool = n.links.index.str.contains("battery charger")
-        vars_link.sel({"Link-ext": nodes + " charger"}) * 1,
+
-        vars_link.sel({"Link-ext": nodes + " discharger"}) * -eff,
+    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()
+    set_of_options = config["solving"]["solver"]["options"]
-    solver_name = solver_options.pop("name")
+    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(
+        status, condition = n.optimize(
-            n,
             solver_name=solver_name,
-            solver_options=solver_options,
             extra_functionality=extra_functionality,
+            **solver_options,
             **kwargs,
         )
     else:
-        optimize_transmission_expansion_iteratively(
+        status, condition = n.optimize.optimize_transmission_expansion_iteratively(
-            n,
             solver_name=solver_name,
-            solver_options=solver_options,
             track_iterations=track_iterations,
             min_iterations=min_iterations,
             max_iterations=max_iterations,
             extra_functionality=extra_functionality,
+            **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",
+            ll="v1.5",
-            opts="Co2L-BAU-24H",
+            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")
+    opts = (snakemake.wildcards.opts + "-" + snakemake.wildcards.sector_opts).split("-")
-    if tmpdir is not None:
+    opts = [o for o in opts if o != ""]
-        Path(tmpdir).mkdir(parents=True, exist_ok=True)
-    opts = snakemake.wildcards.opts.split("-")
     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])
+        if "overrides" in snakemake.input.keys():
-        n = prepare_network(n, solve_opts)
+            overrides = override_component_attrs(snakemake.input.overrides)
-        n = solve_network(
+            n = pypsa.Network(
-            n,
+                snakemake.input.network, override_component_attrs=overrides
-            snakemake.config,
-            opts,
-            solver_dir=tmpdir,
-            solver_logfile=snakemake.log.solver,
             )
+        else:
+            n = pypsa.Network(snakemake.input.network)
+
+        n = prepare_network(n, solve_opts, config=snakemake.config)
+
+        n = solve_network(
+            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])
 

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",
+            ll="v1.5",
-            opts="Co2L-BAU-24H",
+            sector_opts="",
+            planning_horizons="",
         )
+
     configure_logging(snakemake)
+    update_config_with_sector_opts(snakemake.config, snakemake.wildcards.sector_opts)
 
-    tmpdir = snakemake.config["solving"].get("tmpdir")
+    opts = (snakemake.wildcards.opts + "-" + snakemake.wildcards.sector_opts).split("-")
-    if tmpdir is not None:
+    opts = [o for o in opts if o != ""]
-        Path(tmpdir).mkdir(parents=True, exist_ok=True)
+    solve_opts = snakemake.config["solving"]["options"]
 
-    n = pypsa.Network(snakemake.input.unprepared)
+    np.random.seed(solve_opts.get("seed", 123))
-    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
 
     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:
-        n = solve_network(
+            overrides = override_component_attrs(snakemake.input.overrides)
-            n,
+            n = pypsa.Network(
-            snakemake.config,
+                snakemake.input.network, override_component_attrs=overrides
-            opts,
-            solver_dir=tmpdir,
-            solver_logfile=snakemake.log.solver,
             )
+        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, 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])
 

scripts/solve_sector_network.py

@@ -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))