Note

You can download this example as a Jupyter notebook or start it in interactive mode.

Network Clustering

In this example, we show how pypsa can deal with spatial clustering of networks.

[1]:
import pypsa
import re
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import pandas as pd
from pypsa.networkclustering import get_clustering_from_busmap, busmap_by_kmeans
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [1], line 1
----> 1 import pypsa
      2 import re
      3 import numpy as np

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/__init__.py:10
      1 # -*- coding: utf-8 -*-
      4 """
      5 Python for Power Systems Analysis (PyPSA)
      6
      7 Grid calculation library.
      8 """
---> 10 from pypsa import (
     11     components,
     12     contingency,
     13     descriptors,
     14     examples,
     15     geo,
     16     io,
     17     linopf,
     18     linopt,
     19     networkclustering,
     20     opf,
     21     opt,
     22     optimization,
     23     pf,
     24     plot,
     25 )
     26 from pypsa.components import Network, SubNetwork
     28 __version__ = "0.21.2"

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/components.py:50
     37 from pypsa.io import (
     38     export_to_csv_folder,
     39     export_to_hdf5,
   (...)
     47     import_series_from_dataframe,
     48 )
     49 from pypsa.opf import network_lopf, network_opf
---> 50 from pypsa.optimization.optimize import OptimizationAccessor
     51 from pypsa.pf import (
     52     calculate_B_H,
     53     calculate_dependent_values,
   (...)
     62     sub_network_pf,
     63 )
     64 from pypsa.plot import iplot, plot

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/optimization/__init__.py:7
      1 #!/usr/bin/env python3
      2 # -*- coding: utf-8 -*-
      3 """
      4 Build optimisation problems from PyPSA networks with Linopy.
      5 """
----> 7 from pypsa.optimization import abstract, constraints, optimize, variables
      8 from pypsa.optimization.optimize import create_model

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/optimization/constraints.py:9
      6 import logging
      8 import pandas as pd
----> 9 from linopy.expressions import LinearExpression, merge
     10 from numpy import arange, cumsum, inf, nan, roll
     11 from scipy import sparse

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/__init__.py:9
      1 #!/usr/bin/env python3
      2 # -*- coding: utf-8 -*-
      3 """
      4 Created on Wed Mar 10 11:03:06 2021.
      5
      6 @author: fabulous
      7 """
----> 9 from linopy import model, remote
     10 from linopy.expressions import merge
     11 from linopy.io import read_netcdf

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/model.py:22
     20 from linopy import solvers
     21 from linopy.common import best_int, replace_by_map
---> 22 from linopy.constraints import (
     23     AnonymousConstraint,
     24     AnonymousScalarConstraint,
     25     Constraints,
     26 )
     27 from linopy.eval import Expr
     28 from linopy.expressions import LinearExpression, ScalarLinearExpression

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/constraints.py:21
     18 from scipy.sparse import coo_matrix
     19 from xarray import DataArray, Dataset
---> 21 from linopy import expressions, variables
     22 from linopy.common import (
     23     _merge_inplace,
     24     has_assigned_model,
   (...)
     27     replace_by_map,
     28 )
     31 class Constraint(DataArray):

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/expressions.py:23
     20 from xarray.core.dataarray import DataArrayCoordinates
     21 from xarray.core.groupby import _maybe_reorder, peek_at
---> 23 from linopy import constraints, variables
     24 from linopy.common import as_dataarray
     27 def exprwrap(method, *default_args, **new_default_kwargs):

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/variables.py:398
    393     roll = varwrap(DataArray.roll)
    395     rolling = varwrap(DataArray.rolling)
--> 398 @dataclass(repr=False)
    399 class Variables:
    400     """
    401     A variables container used for storing multiple variable arrays.
    402     """
    404     labels: Dataset = Dataset()

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/dataclasses.py:1211, in dataclass.<locals>.wrap(cls)
   1210 def wrap(cls):
-> 1211     return _process_class(cls, init, repr, eq, order, unsafe_hash,
   1212                           frozen, match_args, kw_only, slots,
   1213                           weakref_slot)

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/dataclasses.py:959, in _process_class(cls, init, repr, eq, order, unsafe_hash, frozen, match_args, kw_only, slots, weakref_slot)
    956         kw_only = True
    957     else:
    958         # Otherwise it's a field of some type.
--> 959         cls_fields.append(_get_field(cls, name, type, kw_only))
    961 for f in cls_fields:
    962     fields[f.name] = f

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/dataclasses.py:816, in _get_field(cls, a_name, a_type, default_kw_only)
    812 # For real fields, disallow mutable defaults.  Use unhashable as a proxy
    813 # indicator for mutability.  Read the __hash__ attribute from the class,
    814 # not the instance.
    815 if f._field_type is _FIELD and f.default.__class__.__hash__ is None:
--> 816     raise ValueError(f'mutable default {type(f.default)} for field '
    817                      f'{f.name} is not allowed: use default_factory')
    819 return f

ValueError: mutable default <class 'xarray.core.dataset.Dataset'> for field labels is not allowed: use default_factory
[2]:
n = pypsa.examples.scigrid_de()
n.lines["type"] = np.nan  # delete the 'type' specifications to make this example easier
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [2], line 1
----> 1 n = pypsa.examples.scigrid_de()
      2 n.lines["type"] = np.nan  # delete the 'type' specifications to make this example easier

NameError: name 'pypsa' is not defined

The important information that pypsa needs for spatial clustering is in the busmap. It contains the mapping of which buses should be grouped together, similar to the groupby groups as we know it from pandas.

You can either calculate a busmap from the provided clustering algorithms or you can create/use your own busmap.

Cluster by custom busmap

Let’s start with creating our own. In the following, we group all buses together which belong to the same operator. Buses which do not have a specific operator just stay on its own.

[3]:
groups = n.buses.operator.apply(lambda x: re.split(" |,|;", x)[0])
busmap = groups.where(groups != "", n.buses.index)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [3], line 1
----> 1 groups = n.buses.operator.apply(lambda x: re.split(" |,|;", x)[0])
      2 busmap = groups.where(groups != "", n.buses.index)

NameError: name 'n' is not defined

Now we cluster the network based on the busmap.

[4]:
C = get_clustering_from_busmap(n, busmap)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [4], line 1
----> 1 C = get_clustering_from_busmap(n, busmap)

NameError: name 'get_clustering_from_busmap' is not defined

C is a Clustering object which contains all important information. Among others, the new network is now stored in that Clustering object.

[5]:
nc = C.network
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [5], line 1
----> 1 nc = C.network

NameError: name 'C' is not defined

We have a look at the original and the clustered topology

[6]:
fig, (ax, ax1) = plt.subplots(
    1, 2, subplot_kw={"projection": ccrs.EqualEarth()}, figsize=(12, 12)
)
plot_kwrgs = dict(bus_sizes=1e-3, line_widths=0.5)
n.plot(ax=ax, title="original", **plot_kwrgs)
nc.plot(ax=ax1, title="clustered by operator", **plot_kwrgs)
fig.tight_layout()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [6], line 1
----> 1 fig, (ax, ax1) = plt.subplots(
      2     1, 2, subplot_kw={"projection": ccrs.EqualEarth()}, figsize=(12, 12)
      3 )
      4 plot_kwrgs = dict(bus_sizes=1e-3, line_widths=0.5)
      5 n.plot(ax=ax, title="original", **plot_kwrgs)

NameError: name 'plt' is not defined

Looks a bit messy as over 120 buses do not have assigned operators.

Clustering by busmap created from K-means

Let’s now make a clustering based on the kmeans algorithm. Therefore we calculate the busmap from a non-weighted kmeans clustering.

[7]:
weighting = pd.Series(1, n.buses.index)
busmap2 = busmap_by_kmeans(n, bus_weightings=weighting, n_clusters=50)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [7], line 1
----> 1 weighting = pd.Series(1, n.buses.index)
      2 busmap2 = busmap_by_kmeans(n, bus_weightings=weighting, n_clusters=50)

NameError: name 'pd' is not defined

We use this new kmeans-based busmap to create a new clustered method.

[8]:
C2 = get_clustering_from_busmap(n, busmap2)
nc2 = C2.network
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [8], line 1
----> 1 C2 = get_clustering_from_busmap(n, busmap2)
      2 nc2 = C2.network

NameError: name 'get_clustering_from_busmap' is not defined

Again, let’s plot the networks to compare:

[9]:
fig, (ax, ax1) = plt.subplots(
    1, 2, subplot_kw={"projection": ccrs.EqualEarth()}, figsize=(12, 12)
)
plot_kwrgs = dict(bus_sizes=1e-3, line_widths=0.5)
n.plot(ax=ax, title="original", **plot_kwrgs)
nc2.plot(ax=ax1, title="clustered by kmeans", **plot_kwrgs)
fig.tight_layout()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [9], line 1
----> 1 fig, (ax, ax1) = plt.subplots(
      2     1, 2, subplot_kw={"projection": ccrs.EqualEarth()}, figsize=(12, 12)
      3 )
      4 plot_kwrgs = dict(bus_sizes=1e-3, line_widths=0.5)
      5 n.plot(ax=ax, title="original", **plot_kwrgs)

NameError: name 'plt' is not defined

There are other clustering algorithms in the pipeline of pypsa as the hierarchical clustering which performs better than the kmeans. Also the get_clustering_from_busmap function supports various arguments on how components in the network should be aggregated.