Merge pull request #282 from PyPSA/ci-review

Transition to `linopy` implementation

config.default.yaml

@@ -400,6 +400,10 @@ solving:
 
   solver:
     name: gurobi
+    options: gurobi-default
+
+  solver_options:
+    gurobi-default:
       threads: 4
       method: 2 # barrier
       crossover: 0
@@ -408,14 +412,38 @@ solving:
       AggFill: 0
       PreDual: 0
       GURO_PAR_BARDENSETHRESH: 200
-    #FeasibilityTol: 1.e-6
+      seed: 10              # Consistent seed for all plattforms
+    gurobi-numeric-focus:
+      name: gurobi
+      NumericFocus: 3       # Favour numeric stability over speed
+      method: 2             # barrier
+      crossover: 0          # do not use crossover
+      BarHomogeneous: 1     # Use homogeneous barrier if standard does not converge
+      BarConvTol: 1.e-5
+      FeasibilityTol: 1.e-4
+      OptimalityTol: 1.e-4
+      ObjScale: -0.5
+      threads: 8
+      Seed: 123
+    gurobi-fallback:        # Use gurobi defaults
+      name: gurobi
+      crossover: 0
+      method: 2             # barrier
+      BarHomogeneous: 1     # Use homogeneous barrier if standard does not converge
+      BarConvTol: 1.e-5
+      FeasibilityTol: 1.e-5
+      OptimalityTol: 1.e-5
+      Seed: 123
+      threads: 8
+    cplex-default:
+      threads: 4
+      lpmethod: 4 # barrier
+      solutiontype: 2 # non basic solution, ie no crossover
+      barrier_convergetol: 1.e-5
+      feasopt_tolerance: 1.e-6
+
+    cbc-default: {} # Used in CI 
 
-    #name: cplex
-    #threads: 4
-    #lpmethod: 4 # barrier
-    #solutiontype: 2 # non basic solution, ie no crossover
-    #barrier_convergetol: 1.e-5
-    #feasopt_tolerance: 1.e-6
   mem: 30000 #memory in MB; 20 GB enough for 50+B+I+H2; 100 GB for 181+B+I+H2
 
 

scripts/solve_network.py

@@ -1,16 +1,9 @@
 """Solve network."""
 
 import pypsa
-
 import numpy as np
-import pandas as pd
-
-from pypsa.linopt import get_var, linexpr, define_constraints
-
-from pypsa.linopf import network_lopf, ilopf
 
 from vresutils.benchmark import memory_logger
-
 from helper import override_component_attrs, update_config_with_sector_opts
 
 import logging
@@ -64,22 +57,36 @@ def _add_land_use_constraint_m(n):
     n.generators.p_nom_max.clip(lower=0, inplace=True)
 
 
-def prepare_network(n, solve_opts=None):
+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., 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
-               # intersect between macroeconomic and surveybased
-               # willingness to pay
-               # http://journal.frontiersin.org/article/10.3389/fenrg.2015.00055/full
                p_nom=1e9 # kW
         )
 
@@ -103,177 +110,138 @@ def prepare_network(n, solve_opts=None):
     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")
 
-    chargers_b = n.links.carrier.str.contains("battery charger")
+    dischargers_ext= n.links[discharger_bool].query("p_nom_extendable").index
-    chargers = n.links.index[chargers_b & n.links.p_nom_extendable]
+    chargers_ext= n.links[charger_bool].query("p_nom_extendable").index
-    dischargers = chargers.str.replace("charger", "discharger")
 
-    if chargers.empty or ('Link', 'p_nom') not in n.variables.index:
+    eff = n.links.efficiency[dischargers_ext].values
-        return
+    lhs = n.model["Link-p_nom"].loc[chargers_ext] - n.model["Link-p_nom"].loc[dischargers_ext] * eff
 
-    link_p_nom = get_var(n, "Link", "p_nom")
+    n.model.add_constraints(lhs == 0, name="Link-charger_ratio")
-
-    lhs = linexpr((1,link_p_nom[chargers]),
-                  (-n.links.loc[dischargers, "efficiency"].values,
-                   link_p_nom[dischargers].values))
-
-    define_constraints(n, lhs, "=", 0, 'Link', 'charger_ratio')
 
 
 def add_chp_constraints(n):
 
-    electric_bool = (n.links.index.str.contains("urban central")
+    electric = (n.links.index.str.contains("urban central")
                 & n.links.index.str.contains("CHP")
                 & n.links.index.str.contains("electric"))
-    heat_bool = (n.links.index.str.contains("urban central")
+    heat = (n.links.index.str.contains("urban central")
             & n.links.index.str.contains("CHP")
             & n.links.index.str.contains("heat"))
 
-    electric = n.links.index[electric_bool]
+    electric_ext = n.links[electric].query("p_nom_extendable").index
-    heat = n.links.index[heat_bool]
+    heat_ext = n.links[heat].query("p_nom_extendable").index
 
-    electric_ext = n.links.index[electric_bool & n.links.p_nom_extendable]
+    electric_fix = n.links[electric].query("~p_nom_extendable").index
-    heat_ext = n.links.index[heat_bool & n.links.p_nom_extendable]
+    heat_fix = n.links[heat].query("~p_nom_extendable").index
 
-    electric_fix = n.links.index[electric_bool & ~n.links.p_nom_extendable]
+    p = n.model["Link-p"] # dimension: [time, link]
-    heat_fix = n.links.index[heat_bool & ~n.links.p_nom_extendable]
-
-    link_p = get_var(n, "Link", "p")
 
+    # 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"]
 
-        link_p_nom = get_var(n, "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')
 
-        #ratio of output heat to electricity set by p_nom_ratio
+        rename = {"Link-ext": "Link"}
-        lhs = linexpr((n.links.loc[electric_ext, "efficiency"]
+        lhs = p.loc[:, electric_ext] + p.loc[:, heat_ext] - p_nom.rename(rename).loc[electric_ext]
-                       *n.links.loc[electric_ext, "p_nom_ratio"],
+        n.model.add_constraints(lhs <= 0, name='chplink-top_iso_fuel_line_ext')
-                       link_p_nom[electric_ext]),
-                      (-n.links.loc[heat_ext, "efficiency"].values,
-                       link_p_nom[heat_ext].values))
 
-        define_constraints(n, lhs, "=", 0, 'chplink', 'fix_p_nom_ratio')
-
-        #top_iso_fuel_line for extendable
-        lhs = linexpr((1,link_p[heat_ext]),
-                      (1,link_p[electric_ext].values),
-                      (-1,link_p_nom[electric_ext].values))
-
-        define_constraints(n, lhs, "<=", 0, 'chplink', 'top_iso_fuel_line_ext')
-
-    if not electric_fix.empty:
 
     # top_iso_fuel_line for fixed
-        lhs = linexpr((1,link_p[heat_fix]),
+    if not electric_fix.empty:
-                      (1,link_p[electric_fix].values))
+        lhs = p.loc[:, electric_fix] + p.loc[:, heat_fix]
-
+        rhs = n.links.p_nom[electric_fix]
-        rhs = n.links.loc[electric_fix, "p_nom"].values
+        n.model.add_constraints(lhs <= rhs, name='chplink-top_iso_fuel_line_fix')
-
-        define_constraints(n, lhs, "<=", rhs, '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')
 
-        #backpressure
-        lhs = linexpr((n.links.loc[electric, "c_b"].values
-                       *n.links.loc[heat, "efficiency"],
-                       link_p[heat]),
-                      (-n.links.loc[electric, "efficiency"].values,
-                       link_p[electric].values))
-
-        define_constraints(n, lhs, "<=", 0, 'chplink', 'backpressure')
-
-def basename(x):
-     return x.split("-2")[0]
 
 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
+    if h2_retrofitted_i.empty or gas_pipes_i.empty:
-
-    link_p_nom = get_var(n, "Link", "p_nom")
-
-    CH4_per_H2 = 1 / n.config["sector"]["H2_retrofit_capacity_per_CH4"]
-    fr = "H2 pipeline retrofitted"
-    to = "gas pipeline"
-
-    pipe_capacity = n.links.loc[gas_pipes_i, 'p_nom'].rename(basename)
-
-    lhs = linexpr(
-        (CH4_per_H2, link_p_nom.loc[h2_retrofitted_i].rename(index=lambda x: x.replace(fr, to))),
-        (1, link_p_nom.loc[gas_pipes_i])
-    )
-
-    lhs.rename(basename, inplace=True)
-    define_constraints(n, lhs, "=", pipe_capacity, 'Link', 'pipe_retrofit')
-
-
-def add_co2_sequestration_limit(n, sns):
-
-    co2_stores = n.stores.loc[n.stores.carrier=='co2 stored'].index
-
-    if co2_stores.empty or ('Store', 'e') not in n.variables.index:
         return
 
-    vars_final_co2_stored = get_var(n, 'Store', 'e').loc[sns[-1], co2_stores]
+    p_nom = n.model["Link-p_nom"]
 
-    lhs = linexpr((1, vars_final_co2_stored)).sum()
+    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")
 
-    limit = n.config["sector"].get("co2_sequestration_potential", 200) * 1e6
+    n.model.add_constraints(lhs == rhs, name='Link-pipe_retrofit')
-    for o in opts:
-        if not "seq" in o: continue
-        limit = float(o[o.find("seq")+3:]) * 1e6
-        break
-
-    name = 'co2_sequestration_limit'
-    sense = "<="
-
-    n.add("GlobalConstraint", name, sense=sense, constant=limit,
-          type=np.nan, carrier_attribute=np.nan)
-
-    define_constraints(n, lhs, sense, limit, 'GlobalConstraint',
-                       'mu', axes=pd.Index([name]), spec=name)
 
 
 def extra_functionality(n, snapshots):
     add_battery_constraints(n)
     add_pipe_retrofit_constraint(n)
-    add_co2_sequestration_limit(n, snapshots)
 
 
-def solve_network(n, config, opts='', **kwargs):
+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 {}
-    cf_solving = config['solving']['options']
+    solver_name = config['solving']['solver']['name']
-    track_iterations = cf_solving.get('track_iterations', False)
+    cf_solving = config["solving"]["options"]
-    min_iterations = cf_solving.get('min_iterations', 4)
+    track_iterations = cf_solving.get("track_iterations", False)
-    max_iterations = cf_solving.get('max_iterations', 6)
+    min_iterations = cf_solving.get("min_iterations", 4)
-    keep_shadowprices = cf_solving.get('keep_shadowprices', True)
+    max_iterations = cf_solving.get("max_iterations", 6)
 
     # add to network for extra_functionality
     n.config = config
     n.opts = opts
 
-    if cf_solving.get('skip_iterations', False):
+    skip_iterations = cf_solving.get("skip_iterations", False)
-        network_lopf(n, solver_name=solver_name, solver_options=solver_options,
+    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,
-                     keep_shadowprices=keep_shadowprices, **kwargs)
+            **solver_options,
+            **kwargs,
+        )
     else:
-        ilopf(n, solver_name=solver_name, solver_options=solver_options,
+        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,
-              keep_shadowprices=keep_shadowprices,
+            **solver_options,
-              **kwargs)
+            **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
@@ -281,10 +249,10 @@ if __name__ == "__main__":
             'solve_network',
             simpl='',
             opts="",
-            clusters="37",
+            clusters="5",
             lv=1.0,
-            sector_opts='168H-T-H-B-I-A-solar+p3-dist1',
+            sector_opts='Co2L0-3H-T-H-B-I-A-solar+p3-dist1',
-            planning_horizons="2030",
+            planning_horizons="2050",
         )
 
     logging.basicConfig(filename=snakemake.log.python,
@@ -305,11 +273,9 @@ if __name__ == "__main__":
         overrides = override_component_attrs(snakemake.input.overrides)
         n = pypsa.Network(snakemake.input.network, override_component_attrs=overrides)
 
-        n = prepare_network(n, solve_opts)
+        n = prepare_network(n, solve_opts, config=snakemake.config)
 
-        n = solve_network(n, config=snakemake.config, opts=opts,
+        n = solve_network(n, config=snakemake.config, opts=opts, log_fn=snakemake.log.solver)
-                          solver_dir=tmpdir,
-                          solver_logfile=snakemake.log.solver)
 
         if "lv_limit" in n.global_constraints.index:
             n.line_volume_limit = n.global_constraints.at["lv_limit", "constant"]