Linopy transition

Snakefile

@@ -560,6 +560,7 @@ rule add_electricity:
 
 rule simplify_network:
     params:
+        aggregation_strategies=config["cluster_options"]["aggregation_strategies"],
         renewable=config["renewable"],
         geo_crs=config["crs"]["geo_crs"],
         cluster_options=config["cluster_options"],

config.default.yaml

@@ -65,7 +65,7 @@ cluster_options:
     remove_stubs_across_borders: true
     p_threshold_drop_isolated: 20 # [MW] isolated buses are being discarded if bus mean power is below the specified threshold
     p_threshold_merge_isolated: 300 # [MW] isolated buses are being merged into a single isolated bus if a bus mean power is below the specified threshold
-    s_threshold_fetch_isolated: 0.05 # [-] a share of the national load for merging an isolated network into a backbone network
+    s_threshold_fetch_isolated: false # [-] a share of the national load for merging an isolated network into a backbone network
   cluster_network:
     algorithm: kmeans
     feature: solar+onwind-time
@@ -396,8 +396,6 @@ monte_carlo:
       type: beta
       args: [0.5, 2]
 
-
-
 solving:
   options:
     formulation: kirchhoff
@@ -411,15 +409,63 @@ solving:
     #nhours: 10
   solver:
     name: gurobi
-    threads: 4
-    method: 2 # barrier (=ipm)
-    crossover: 0
-    BarConvTol: 1.e-5
-    FeasibilityTol: 1.e-6
-    AggFill: 0
-    PreDual: 0
-    GURO_PAR_BARDENSETHRESH: 200
-
+    options: gurobi-default
+
+  solver_options:
+    highs-default:
+      # refer to https://ergo-code.github.io/HiGHS/options/definitions.html#solver
+      threads: 4
+      solver: "ipm"
+      run_crossover: "off"
+      small_matrix_value: 1e-6
+      large_matrix_value: 1e9
+      primal_feasibility_tolerance: 1e-5
+      dual_feasibility_tolerance: 1e-5
+      ipm_optimality_tolerance: 1e-4
+      parallel: "on"
+      random_seed: 123
+    gurobi-default:
+      threads: 4
+      method: 2 # barrier
+      crossover: 0
+      BarConvTol: 1.e-6
+      Seed: 123
+      AggFill: 0
+      PreDual: 0
+      GURO_PAR_BARDENSETHRESH: 200
+    gurobi-numeric-focus:
+      name: gurobi
+      NumericFocus: 3       # Favour numeric stability over speed
+      method: 2             # barrier
+      crossover: 0          # do not use crossover
+      BarHomogeneous: 1     # Use homogeneous barrier if standard does not converge
+      BarConvTol: 1.e-5
+      FeasibilityTol: 1.e-4
+      OptimalityTol: 1.e-4
+      ObjScale: -0.5
+      threads: 8
+      Seed: 123
+    gurobi-fallback:        # Use gurobi defaults
+      name: gurobi
+      crossover: 0
+      method: 2             # barrier
+      BarHomogeneous: 1     # Use homogeneous barrier if standard does not converge
+      BarConvTol: 1.e-5
+      FeasibilityTol: 1.e-5
+      OptimalityTol: 1.e-5
+      Seed: 123
+      threads: 8
+    cplex-default:
+      threads: 4
+      lpmethod: 4 # barrier
+      solutiontype: 2 # non basic solution, ie no crossover
+      barrier.convergetol: 1.e-4 #5
+      feasopt.tolerance: 1.e-5 #6
+    cbc-default: {} # Used in CI
+    glpk-default: {} # Used in CI
+
+  mem_mb: 30000 #memory in MB; 20 GB enough for 50+B+I+H2; 100 GB for 181+B+I+H2
+  runtime: 6h #runtime in humanfriendly style https://humanfriendly.readthedocs.io/en/latest/
 
 plotting:
   map:

config.tutorial.yaml

@@ -78,7 +78,7 @@ cluster_options:
     remove_stubs_across_borders: true
     p_threshold_drop_isolated: 20 # [MW] isolated buses are being discarded if bus mean power is below the specified threshold
     p_threshold_merge_isolated: 300 # [MW] isolated buses are being merged into a single isolated bus if a bus mean power is below the specified threshold
-    s_threshold_fetch_isolated: 0.05 # [-] a share of the national load for merging an isolated network into a backbone network
+    s_threshold_fetch_isolated: false # [-] a share of the national load for merging an isolated network into a backbone network
   cluster_network:
     algorithm: kmeans
     feature: solar+onwind-time
@@ -405,6 +405,63 @@ solving:
     #nhours: 10
   solver:
     name: glpk
+    options: glpk-default
+
+  solver_options:
+    highs-default:
+      # refer to https://ergo-code.github.io/HiGHS/options/definitions.html#solver
+      threads: 4
+      solver: "ipm"
+      run_crossover: "off"
+      small_matrix_value: 1e-6
+      large_matrix_value: 1e9
+      primal_feasibility_tolerance: 1e-5
+      dual_feasibility_tolerance: 1e-5
+      ipm_optimality_tolerance: 1e-4
+      parallel: "on"
+      random_seed: 123
+    gurobi-default:
+      threads: 4
+      method: 2 # barrier
+      crossover: 0
+      BarConvTol: 1.e-6
+      Seed: 123
+      AggFill: 0
+      PreDual: 0
+      GURO_PAR_BARDENSETHRESH: 200
+    gurobi-numeric-focus:
+      name: gurobi
+      NumericFocus: 3       # Favour numeric stability over speed
+      method: 2             # barrier
+      crossover: 0          # do not use crossover
+      BarHomogeneous: 1     # Use homogeneous barrier if standard does not converge
+      BarConvTol: 1.e-5
+      FeasibilityTol: 1.e-4
+      OptimalityTol: 1.e-4
+      ObjScale: -0.5
+      threads: 8
+      Seed: 123
+    gurobi-fallback:        # Use gurobi defaults
+      name: gurobi
+      crossover: 0
+      method: 2             # barrier
+      BarHomogeneous: 1     # Use homogeneous barrier if standard does not converge
+      BarConvTol: 1.e-5
+      FeasibilityTol: 1.e-5
+      OptimalityTol: 1.e-5
+      Seed: 123
+      threads: 8
+    cplex-default:
+      threads: 4
+      lpmethod: 4 # barrier
+      solutiontype: 2 # non basic solution, ie no crossover
+      barrier.convergetol: 1.e-4 #5
+      feasopt.tolerance: 1.e-5 #6
+    cbc-default: {} # Used in CI
+    glpk-default: {} # Used in CI
+
+  mem_mb: 30000 #memory in MB; 20 GB enough for 50+B+I+H2; 100 GB for 181+B+I+H2
+  runtime: 6h #runtime in humanfriendly style https://humanfriendly.readthedocs.io/en/latest/
 
 
 plotting:

doc/release_notes.rst

@@ -24,6 +24,8 @@ E.g. if a new rule becomes available describe how to use it `snakemake -j1 run_t
 
 * Add an option to use csv format for custom demand imports. `PR #995 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/995>`__
 
+* Re-implement network optimisation using linopy methods. `PR #796 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/796>`__
+
 **Minor Changes and bug-fixing**
 
 * Minor bug-fixing to run the cluster wildcard min `PR #1019 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/1019>`__

envs/environment.yaml

@@ -12,7 +12,7 @@ dependencies:
 - pip
 - mamba   # esp for windows build
 
-- pypsa>=0.24, <0.25
+- pypsa
 # - atlite>=0.2.4  # until https://github.com/PyPSA/atlite/issues/244 is not merged
 - dask
 - powerplantmatching>=0.5.7
@@ -27,6 +27,7 @@ dependencies:
 - memory_profiler
 - ruamel.yaml<=0.17.26
 - pytables
+- pyscipopt # added to compy with the quadratic objective requirement of the clustering script
 - lxml
 - numpy
 - pandas
@@ -77,7 +78,7 @@ dependencies:
 
 # Default solver for tests (required for CI)
 - glpk
-- ipopt<3.13.3
+- ipopt #<3.13.3
 - gurobi
 
 - pip:

scripts/_helpers.py

@@ -453,31 +453,6 @@ def dlProgress(count, blockSize, totalSize, roundto=roundto):
     urllib.request.urlretrieve(url, file, reporthook=dlProgress, data=data)
 
 
-def get_aggregation_strategies(aggregation_strategies):
-    """
-    Default aggregation strategies that cannot be defined in .yaml format must
-    be specified within the function, otherwise (when defaults are passed in
-    the function's definition) they get lost when custom values are specified
-    in the config.
-    """
-    import numpy as np
-
-    # to handle the new version of PyPSA.
-    try:
-        from pypsa.clustering.spatial import _make_consense
-    except Exception:
-        # TODO: remove after new release and update minimum pypsa version
-        from pypsa.clustering.spatial import _make_consense
-
-    bus_strategies = dict(country=_make_consense("Bus", "country"))
-    bus_strategies.update(aggregation_strategies.get("buses", {}))
-
-    generator_strategies = {"build_year": lambda x: 0, "lifetime": lambda x: np.inf}
-    generator_strategies.update(aggregation_strategies.get("generators", {}))
-
-    return bus_strategies, generator_strategies
-
-
 def mock_snakemake(rulename, **wildcards):
     """
     This function is expected to be executed from the "scripts"-directory of "

scripts/base_network.py

@@ -513,6 +513,18 @@ def base_network(
 
     if hvdc_as_lines_config:
         lines = pd.concat([lines_ac, lines_dc])
+        lines.drop(
+            labels=[
+                "bus0_lon",
+                "bus0_lat",
+                "bus1_lon",
+                "bus1_lat",
+                "bus_0_coors",
+                "bus_1_coors",
+            ],
+            axis=1,
+            inplace=True,
+        )
         n.import_components_from_dataframe(lines, "Line")
     else:
         lines_dc = _set_electrical_parameters_links(links_config, lines_dc)
@@ -522,6 +534,18 @@ def base_network(
             axis=1,
             result_type="reduce",
         )
+        lines_ac.drop(
+            labels=[
+                "bus0_lon",
+                "bus0_lat",
+                "bus1_lon",
+                "bus1_lat",
+                "bus_0_coors",
+                "bus_1_coors",
+            ],
+            axis=1,
+            inplace=True,
+        )
         n.import_components_from_dataframe(lines_ac, "Line")
         # The columns which names starts with "bus" are mixed up with the third-bus specification
         # when executing additional_linkports()

scripts/cluster_network.py

@@ -125,15 +125,14 @@
 from functools import reduce
 
 import geopandas as gpd
+import linopy
 import numpy as np
 import pandas as pd
-import pyomo.environ as po
 import pypsa
 from _helpers import (
     REGION_COLS,
     configure_logging,
     create_logger,
-    get_aggregation_strategies,
     sets_path_to_root,
     update_p_nom_max,
 )
@@ -336,47 +335,22 @@ def distribute_clusters(
         distribution_factor.sum(), 1.0, rtol=1e-3
     ), f"Country weights L must sum up to 1.0 when distributing clusters. Is {distribution_factor.sum()}."
 
-    m = po.ConcreteModel()
-
-    def n_bounds(model, *n_id):
-        """
-        Create a function that makes a bound pair for pyomo.
-
-        Use n_bounds(model, n_id) if N is Single-Index
-        Use n_bounds(model, *n_id) if N is Multi-Index
-        Example: https://pyomo.readthedocs.io/en/stable/pyomo_modeling_components/Variables.html
-
-        Returns
-        -------
-        bounds = A function (or Python object) that gives a (lower,upper) bound pair i.e.(1,10) for the variable
-        """
-        return (1, N[n_id])
-
-    m.n = po.Var(list(distribution_factor.index), bounds=n_bounds, domain=po.Integers)
-    m.tot = po.Constraint(expr=(po.summation(m.n) == n_clusters))
-    m.objective = po.Objective(
-        expr=sum(
-            (m.n[i] - distribution_factor.loc[i] * n_clusters) ** 2
-            for i in distribution_factor.index
-        ),
-        sense=po.minimize,
+    m = linopy.Model()
+    clusters = m.add_variables(
+        lower=1, upper=N, coords=[L.index], name="n", integer=True
     )
-
-    opt = po.SolverFactory(solver_name)
-    if not opt.has_capability("quadratic_objective"):
-        logger.warning(
-            f"The configured solver `{solver_name}` does not support quadratic objectives. Falling back to `ipopt`."
+    m.add_constraints(clusters.sum() == n_clusters, name="tot")
+    # leave out constant in objective (L * n_clusters) ** 2
+    m.objective = (clusters * clusters - 2 * clusters * L * n_clusters).sum()
+    if solver_name == "gurobi":
+        logging.getLogger("gurobipy").propagate = False
+    elif solver_name != "scip":
+        logger.info(
+            f"The configured solver `{solver_name}` does not support quadratic objectives. Falling back to `scip`."
         )
-        opt = po.SolverFactory("ipopt")
-
-    results = opt.solve(m)
-    assert (
-        results["Solver"][0]["Status"] == "ok"
-    ), f"Solver returned non-optimally: {results}"
-
-    return (
-        pd.Series(m.n.get_values(), index=distribution_factor.index).round().astype(int)
-    )
+        solver_name = "scip"
+    m.solve(solver_name=solver_name)
+    return m.solution["n"].to_series().astype(int)
 
 
 def busmap_for_gadm_clusters(inputs, n, gadm_level, geo_crs, country_list):
@@ -577,9 +551,10 @@ def clustering_for_n_clusters(
     extended_link_costs=0,
     focus_weights=None,
 ):
-    bus_strategies, generator_strategies = get_aggregation_strategies(
-        aggregation_strategies
-    )
+    line_strategies = aggregation_strategies.get("lines", dict())
+    line_strategies.update({"geometry": "first", "bounds": "first"})
+    generator_strategies = aggregation_strategies.get("generators", dict())
+    one_port_strategies = aggregation_strategies.get("one_ports", dict())
 
     if not isinstance(custom_busmap, pd.Series):
         if alternative_clustering:
@@ -605,12 +580,20 @@ def clustering_for_n_clusters(
     clustering = get_clustering_from_busmap(
         n,
         busmap,
-        bus_strategies=bus_strategies,
         aggregate_generators_weighted=True,
         aggregate_generators_carriers=aggregate_carriers,
         aggregate_one_ports=["Load", "StorageUnit"],
         line_length_factor=line_length_factor,
+        line_strategies=line_strategies,
         generator_strategies=generator_strategies,
+        bus_strategies={
+            "lat": "mean",
+            "lon": "mean",
+            "tag_substation": "first",
+            "tag_area": "first",
+            "country": "first",
+        },
+        one_port_strategies=one_port_strategies,
         scale_link_capital_costs=False,
     )