Introduction to networkx#

NetworkX: Network Analysis in Python

NetworkX is a Python package for the creation, manipulation, and study of the structure, dynamics, and functions of complex networks.

Main features:

  • Data structures for graphs, digraphs, and multigraphs

  • Many standard graph algorithms

  • Network structure and analysis measures

Note

Documentation for this package is available at https://networkx.github.io/.

Note

If you have not yet set up Python on your computer, you can execute this tutorial in your browser via Google Colab. Click on the rocket in the top right corner and launch “Colab”. If that doesn’t work download the .ipynb file and import it in Google Colab.

Then install the following packages by executing the following command in a Jupyter cell at the top of the notebook.

!pip install numpy networkx pandas matplotlib

Let’s perform some network analysis on this simple graph:

Network Analysis with numpy.#

import numpy as np

Say we want to calculate the Laplacian \(L\) of this graph based on its incidence matrix \(K\), which is an \(N\times N\) matrix defined as \(L=KK^\top\) for an undirected graph. The Laplacian matrix of a graph is a representation that captures the connectivity and structure of the graph by quantifying the difference between the degree of each vertex and the adjacency relationships with its neighbouring vertices. We first need to write down the incidence matrix \(K\) as a np.array. Let’s also use the convention that edges are oriented such that they are directed at the node with the higher label value (i.e. from node 1 to node 2, not vice versa).

K = np.array(
    [
        [1, -1, 0, 0],
        [1, 0, -1, 0],
        [0, 1, -1, 0],
        [0, 0, 1, -1],
    ]
).T
K
array([[ 1,  1,  0,  0],
       [-1,  0,  1,  0],
       [ 0, -1, -1,  1],
       [ 0,  0,  0, -1]])

To calculate the Laplacian matrix, we form a dot product of the incidence matrix with its transpose:

L = K.dot(K.T)
L
array([[ 2, -1, -1,  0],
       [-1,  2, -1,  0],
       [-1, -1,  3, -1],
       [ 0,  0, -1,  1]])

This is all fine for small graphs, but inconvient for larger graphs. Let’s take the help some Python packages have to offer…

Making our life easier with networkx#

First, let’s import the library. It is commonly imported under the alias nx.

import networkx as nx

This is how we can create an empty graph with no nodes and no edges.

G = nx.Graph()

Nodes#

We can add one node at a time,

G.add_node(1)

with attributes

G.add_node(2, country="DE")

or add nodes from a list

G.add_nodes_from([3, 4])

We can also add nodes along with node attributes if your container yields 2-tuples of the form (node, node_attribute_dict):

G.add_nodes_from(
    [
        (5, {"color": "red"}),
        (6, {"color": "green"}),
    ]
)

Edges#

G can also be grown by adding one edge at a time,

G.add_edge(1, 2)

even with attributes

G.add_edge(3, 4, weight=2)

or by adding a list of edges,

G.add_edges_from([(1, 3), (2, 5)])

or as a 3-tuple with 2 nodes followed by an edge attribute dictionary

G.add_edges_from([(2, 3, {"weight": 3})])

Examining elements of a graph#

We can examine the nodes and edges.

G.nodes
NodeView((1, 2, 3, 4, 5, 6))
G.number_of_nodes()
6
G.edges
EdgeView([(1, 2), (1, 3), (2, 5), (2, 3), (3, 4)])
G.number_of_edges()
5

Accessing graph elements#

Access an edge:

G.edges[2, 3]
{'weight': 3}

Access an attribute of an edge:

G.edges[2, 3]["weight"]
3

Find all neighbours of node 1:

G[1]
AtlasView({2: {}, 3: {}})

Removing elements#

One can remove nodes and edges from the graph in a similar fashion to adding. Use methods G.remove_node(), G.remove_nodes_from(), G.remove_edge() and G.remove_edges_from(), e.g.

G.remove_node(5)
G.remove_edge(2, 5)
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.9/x64/lib/python3.11/site-packages/networkx/classes/graph.py:1127, in Graph.remove_edge(self, u, v)
   1126 try:
-> 1127     del self._adj[u][v]
   1128     if u != v:  # self-loop needs only one entry removed

KeyError: 5

The above exception was the direct cause of the following exception:

NetworkXError                             Traceback (most recent call last)
Cell In[23], line 1
----> 1 G.remove_edge(2, 5)

File /opt/hostedtoolcache/Python/3.11.9/x64/lib/python3.11/site-packages/networkx/classes/graph.py:1131, in Graph.remove_edge(self, u, v)
   1129         del self._adj[v][u]
   1130 except KeyError as err:
-> 1131     raise NetworkXError(f"The edge {u}-{v} is not in the graph") from err
   1132 nx._clear_cache(self)

NetworkXError: The edge 2-5 is not in the graph

NB: Removing a node will also remove all adjacent edges!

You can remove all nodes and edges with

# G.clear()

Visualising graphs#

A basic drawing function for networks is also available in networkx

nx.draw(G)

with options for labeling graph elements

nx.draw(G, with_labels=True)

and integration to matplotlib

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(5, 5))
nx.draw(G, with_labels=True, ax=ax, node_color="green", font_weight="bold")
plt.savefig("tmp.png")

Analysing graphs#

The networkx library comes with many functions to analyse graphs. Here are a few examples we’re going to need for linearised power flow calculations in electricity transmission networks:

Are all nodes in the network connected with each other?

nx.is_connected(G)

What are the components that are connected / isolated?

list(nx.connected_components(G))

Is the network planar? I.e. can the graph be drawn such that edges don’t cross?

nx.is_planar(G)

What is the frequency of degrees in the network?

nx.degree_histogram(G)
import pandas as pd

pd.Series(nx.degree_histogram(G))

What is the adjacency matrix? (Careful, networkx will yield a weighted adjacency matrix by default!)

A = nx.adjacency_matrix(G, weight=None).todense()
A

What is the incidence matrix? (Careful, networkx will yield a incidence matrix without orientation by default!)

nx.incidence_matrix(G, oriented=True).todense()

What is the Laplacian matrix? (Careful, networkx will yield a weighted adjacency matrix by default!)

L = nx.laplacian_matrix(G, weight=None).todense()
L

What is the degree matrix?

G.degree()

Here, we need to do some work to produce a diagonal matrix from the output of G.degree():

degrees = [val for node, val in G.degree()]
D = np.diag(degrees)
D

Let’s do a few checks using different formulas for the Laplacian matrix:

\[L = KK^\top\]
K.dot(K.T)
assert (K.dot(K.T) == L).all()
\[L = D-A\]
D - A
assert (D - A == L).all()

Find a cycle basis (i.e. a collection of independent cycles through which all other cycles can be represented through linear combination):

nx.cycle_basis(G)

This function returns a list of sequences. Each sequence indicates a series of nodes to traverse for the respective cycle.

Linearised power flow in the European transmission network#

In this example, we are going to load a slightly simplified dataset of the European high-voltage transmission network. In this dataset, HVDC links and back-to-back converters have been left out, and the data only shows AC transmission lines.

url = "https://tubcloud.tu-berlin.de/s/898dEKqG6XEDqqS/download/nodes.csv"
nodes = pd.read_csv(url, index_col=0)
nodes.head(5)
url = "https://tubcloud.tu-berlin.de/s/FmFrJkiWpg2QcQA/download/edges.csv"
edges = pd.read_csv(url, index_col=0)
edges.head(5)
edges

networkx provides a utility function nx.from_pandas_edgelist() to build a network from edges listed in a pandas.DataFrame:

N = nx.from_pandas_edgelist(edges, "bus0", "bus1", edge_attr=["x_pu", "s_nom"])

We can get some basic info about the graph:

print(N)
degrees = [val for node, val in N.degree()]
np.mean(degrees)

The nodes DataFrame provides us with the coordinates of the graph’s nodes. To make the nx.draw() function use these coordinates, we need to bring them into a particular format.

{nodeA: (x, y),
 nodeB: (x, y),
 nodeC: (x, y)}
pos = nodes.apply(tuple, axis=1).to_dict()

Let’s just look at the first 5 elements of the dictionary to check:

{k: pos[k] for k in list(pos.keys())[:5]}

Now, we can draw the European transmission network:

fig, ax = plt.subplots(figsize=(10, 10))
nx.draw(N, pos=pos, node_size=0)

You can already see that not all parts of the Network are connected with each other. Ireland, Great Britain, Scandinavia, the Baltics and some Islands in the Mediterranean are not connected to the continental grid. At least not via AC transmission lines. They are through HVDC links and back-to-back converters. These subgraphs denote the different synchronous zones of the European transmission network:

len(list(nx.connected_components(N)))

Let’s build subgraphs for the synchronous zones by iterating over the connected components:

subgraphs = []
for c in nx.connected_components(N):
    subgraphs.append(N.subgraph(c).copy())

We can now color-code them in the network plot:

fig, ax = plt.subplots(figsize=(10, 10))
colors = ["red", "blue", "green", "orange", "teal", "cyan", "black"]
for i, sub in enumerate(subgraphs):
    sub_pos = {k: v for k, v in pos.items() if k in sub.nodes}
    nx.draw(sub, pos=sub_pos, node_size=0, edge_color=colors[i])

Let’s checkout the synchronous zone of Ireland. We want to compute the Power Transfer Distribution Factor (PTDF) matrix for this area.

The PTDF matrix measures the sensitivity of power flows in each transmission line relative to incremental changes in nodal power injections or withdrawals throughout the electricity network. It is a matrix representation showing how changes in power injections or withdrawals at various nodes in a power grid affect the flow in each of its transmission lines. It is an alternative formulation for linearised power flow to the cycle-based approache from the lecture.

\[f_\ell = \frac{1}{x_\ell}\sum_{i,j} K_{i\ell} (L^{-1})_{ij} p_j\]
\[f_\ell = \sum_j \text{PTDF}_{\ell j} p_j\]
IE = subgraphs[4]

Incidence matrix:

K = nx.incidence_matrix(IE, oriented=True).todense()
K
plt.imshow(K, cmap="RdBu")

Unweighted Laplacian:

L = nx.laplacian_matrix(IE, weight=None).todense()
L
plt.imshow(L, cmap="viridis")

If we want to calculate the weighted Laplacian \(L = KBK^\top\), we also need the diagonal matrix \(B\) for the susceptances (i.e inverse of reactances; \(1/x_\ell\)):

x_pu = nx.get_edge_attributes(IE, "x_pu").values()
b = [1 / x for x in x_pu]
B = np.diag(b)
plt.imshow(B, cmap="Blues", vmax=30000)

Now, we can calculate the weighted Laplacian:

H = B.dot(K.T)
L = K.dot(H)
plt.imshow(L, cmap="RdBu", vmax=2e4, vmin=-2e4)

The weighted Laplacian of the network is not invertible, but we can use the Moore Penrose pseudo-inverse

L_inv = np.linalg.pinv(L)

or apply a trick by disregarding the first column and row of the weighted Laplacian for inversion:

n_nodes = L.shape[0]

L_inv = np.linalg.inv(L[1:, 1:])
L_inv = np.hstack((np.zeros((n_nodes - 1, 1)), L_inv))
L_inv = np.vstack((np.zeros(n_nodes), L_inv))
plt.imshow(L_inv, cmap="viridis")

Now, we have all the ingredients to calculate the PTDF matrix

PTDF = H.dot(L_inv)
plt.imshow(PTDF, cmap="RdBu", vmin=-1, vmax=1)

This matrix, we can now use to calculate the (linearised) power flow in transmission lines based on nodal power imbalances. Here’s an example:

imbalance = np.linspace(-10, 10, len(IE.nodes))
imbalance
imbalance.sum()
flows = PTDF.dot(imbalance)
flows

Exercise 1#

Build the following graph in networkx using the functions to add graphs. The weight of each edge should correspond to the sum of the node labels it connects.

Hide code cell content
G = nx.Graph()
Hide code cell content
edges = [
    (0, 1, dict(weight=1)),
    (0, 2, dict(weight=2)),
    (0, 3, dict(weight=3)),
    (0, 4, dict(weight=4)),
    (1, 2, dict(weight=3)),
    (1, 4, dict(weight=5)),
    (2, 3, dict(weight=5)),
    (3, 4, dict(weight=7)),
]
Hide code cell content
G.add_edges_from(edges)
Hide code cell content
nx.draw(G, with_labels=True)

Exercise 2#

For each subgraph (each represents a synchronous zone), determine

  • the number of transmission lines (aka edges) and buses (aka nodes)

  • whether the network is planar

  • the average number of transmission lines connecting to a bus

  • the number of cycles forming the cycle basis

  • the histogram of the length of the cycles (i.e. number of edges per cycle) in the cycle basis

  • the average length of the cycles in the cycle basis?

  • the adjacency matrix

  • the incidence matrix

  • the weighted Laplacian (weighted by line susceptance = 1 / line reactance)

  • the PTDF matrix

  • the network flows if the first bus generates 1 unit of power and the last bus consumes 1 unit of power

  • Show the absolute flows in a plot using the .draw function of networkx and the keyword argument edge_color. Plot the network on a figure with a size of 10x10 inches. Change the width of the edges to the value 4.

Hide code cell content
G = subgraphs[3]  # iterate here
Hide code cell content
len(G.nodes)
Hide code cell content
len(G.edges)
Hide code cell content
nx.is_planar(G)
Hide code cell content
pd.Series({k: v for k, v in G.degree}).mean()
Hide code cell content
len(nx.cycle_basis(G))
Hide code cell content
cycle_length = pd.Series([len(c) for c in nx.cycle_basis(G)])
Hide code cell content
cycle_length.plot.hist(bins=100)
Hide code cell content
# alternative
pd.Series([len(c) for c in nx.cycle_basis(G)]).value_counts().sort_index().plot()
Hide code cell content
cycle_length.describe()
Hide code cell content
nx.adjacency_matrix(G).todense()
Hide code cell content
K = nx.incidence_matrix(G).todense()
K
Hide code cell content
x_pu = nx.get_edge_attributes(G, "x_pu").values()
b = [1 / x for x in x_pu]
B = np.diag(b)
Hide code cell content
H = B.dot(K.T)
Hide code cell content
L = K.dot(H)
Hide code cell content
L_inv = np.linalg.pinv(L)
Hide code cell content
PTDF = H.dot(L_inv)
Hide code cell content
imbalance = np.zeros(len(G.nodes))
imbalance[0] = 1
imbalance[-1] = -1
Hide code cell content
flows = PTDF.dot(imbalance)
Hide code cell content
abs_flows = np.abs(flows)
Hide code cell content
sub_pos = {k: v for k, v in pos.items() if k in G.nodes}
fig, ax = plt.subplots(figsize=(10, 10))
nx.draw(G, pos=sub_pos, node_size=0, edge_color=abs_flows, width=4)