Merge pull request #903 from PyPSA/cluster-network-replace-pyomo

Cluster network replace pyomo
This commit is contained in:
Fabian Hofmann 2024-01-31 13:31:59 +01:00 committed by GitHub
commit 0bb70b42ef
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
4 changed files with 27 additions and 40 deletions

View File

@ -53,16 +53,6 @@ jobs:
run: | run: |
echo -ne "url: ${CDSAPI_URL}\nkey: ${CDSAPI_TOKEN}\n" > ~/.cdsapirc echo -ne "url: ${CDSAPI_URL}\nkey: ${CDSAPI_TOKEN}\n" > ~/.cdsapirc
- name: Add solver to environment
run: |
echo -e "- glpk\n- ipopt<3.13.3" >> envs/environment.yaml
if: ${{ matrix.os }} == 'windows-latest'
- name: Add solver to environment
run: |
echo -e "- glpk\n- ipopt" >> envs/environment.yaml
if: ${{ matrix.os }} != 'windows-latest'
- name: Setup micromamba - name: Setup micromamba
uses: mamba-org/setup-micromamba@v1 uses: mamba-org/setup-micromamba@v1
with: with:

View File

@ -60,6 +60,11 @@ Upcoming Release
* The rule ``plot_network`` has been split into separate rules for plotting * The rule ``plot_network`` has been split into separate rules for plotting
electricity, hydrogen and gas networks. electricity, hydrogen and gas networks.
* To determine the optimal topology to meet the number of clusters, the workflow used pyomo in combination with ``ipopt`` or ``gurobi``. This dependency has been replaced by using ``linopy`` in combination with ``scipopt`` or ``gurobi``. The environment file has been updated accordingly.
* The ``highs`` solver was added to the default environment file.
PyPSA-Eur 0.9.0 (5th January 2024) PyPSA-Eur 0.9.0 (5th January 2024)
================================== ==================================

View File

@ -35,8 +35,9 @@ dependencies:
- netcdf4 - netcdf4
- networkx - networkx
- scipy - scipy
- glpk
- shapely>=2.0 - shapely>=2.0
- pyomo - pyscipopt
- matplotlib - matplotlib
- proj - proj
- fiona - fiona
@ -47,7 +48,6 @@ dependencies:
- tabula-py - tabula-py
- pyxlsb - pyxlsb
- graphviz - graphviz
- ipopt
# Keep in conda environment when calling ipython # Keep in conda environment when calling ipython
- ipython - ipython
@ -60,3 +60,4 @@ dependencies:
- pip: - pip:
- tsam>=2.3.1 - tsam>=2.3.1
- highspy

View File

@ -122,14 +122,15 @@ Exemplary unsolved network clustered to 37 nodes:
""" """
import logging import logging
import os
import warnings import warnings
from functools import reduce from functools import reduce
import geopandas as gpd import geopandas as gpd
import linopy
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
import pandas as pd import pandas as pd
import pyomo.environ as po
import pypsa import pypsa
import seaborn as sns import seaborn as sns
from _helpers import configure_logging, update_p_nom_max from _helpers import configure_logging, update_p_nom_max
@ -214,7 +215,7 @@ def get_feature_for_hac(n, buses_i=None, feature=None):
return feature_data return feature_data
def distribute_clusters(n, n_clusters, focus_weights=None, solver_name="cbc"): def distribute_clusters(n, n_clusters, focus_weights=None, solver_name="scip"):
""" """
Determine the number of clusters per country. Determine the number of clusters per country.
""" """
@ -254,31 +255,22 @@ def distribute_clusters(n, n_clusters, focus_weights=None, solver_name="cbc"):
L.sum(), 1.0, rtol=1e-3 L.sum(), 1.0, rtol=1e-3
), f"Country weights L must sum up to 1.0 when distributing clusters. Is {L.sum()}." ), f"Country weights L must sum up to 1.0 when distributing clusters. Is {L.sum()}."
m = po.ConcreteModel() m = linopy.Model()
clusters = m.add_variables(
def n_bounds(model, *n_id): lower=1, upper=N, coords=[L.index], name="n", integer=True
return (1, N[n_id])
m.n = po.Var(list(L.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] - L.loc[i] * n_clusters) ** 2 for i in L.index),
sense=po.minimize,
) )
m.add_constraints(clusters.sum() == n_clusters, name="tot")
opt = po.SolverFactory(solver_name) # leave out constant in objective (L * n_clusters) ** 2
if solver_name == "appsi_highs" or not opt.has_capability("quadratic_objective"): m.objective = (clusters * clusters - 2 * clusters * L * n_clusters).sum()
logger.warning( if solver_name == "gurobi":
f"The configured solver `{solver_name}` does not support quadratic objectives. Falling back to `ipopt`." 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") solver_name = "scip"
m.solve(solver_name=solver_name)
results = opt.solve(m) return m.solution["n"].to_series().astype(int)
assert (
results["Solver"][0]["Status"] == "ok"
), f"Solver returned non-optimally: {results}"
return pd.Series(m.n.get_values(), index=L.index).round().astype(int)
def busmap_for_n_clusters( def busmap_for_n_clusters(
@ -372,7 +364,7 @@ def busmap_for_n_clusters(
return ( return (
n.buses.groupby(["country", "sub_network"], group_keys=False) n.buses.groupby(["country", "sub_network"], group_keys=False)
.apply(busmap_for_country) .apply(busmap_for_country, include_groups=False)
.squeeze() .squeeze()
.rename("busmap") .rename("busmap")
) )
@ -385,7 +377,7 @@ def clustering_for_n_clusters(
aggregate_carriers=None, aggregate_carriers=None,
line_length_factor=1.25, line_length_factor=1.25,
aggregation_strategies=dict(), aggregation_strategies=dict(),
solver_name="cbc", solver_name="scip",
algorithm="hac", algorithm="hac",
feature=None, feature=None,
extended_link_costs=0, extended_link_costs=0,
@ -462,7 +454,6 @@ if __name__ == "__main__":
params = snakemake.params params = snakemake.params
solver_name = snakemake.config["solving"]["solver"]["name"] solver_name = snakemake.config["solving"]["solver"]["name"]
solver_name = "appsi_highs" if solver_name == "highs" else solver_name
n = pypsa.Network(snakemake.input.network) n = pypsa.Network(snakemake.input.network)