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 so-called 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]])
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.10/x64/lib/python3.11/site-packages/networkx/classes/graph.py:1134, in Graph.remove_edge(self, u, v)
1133 try:
-> 1134 del self._adj[u][v]
1135 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.10/x64/lib/python3.11/site-packages/networkx/classes/graph.py:1138, in Graph.remove_edge(self, u, v)
1136 del self._adj[v][u]
1137 except KeyError as err:
-> 1138 raise NetworkXError(f"The edge {u}-{v} is not in the graph") from err
1139 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")
Note
For a full list of arguments of the function see https://networkx.org/documentation/stable/reference/generated/networkx.drawing.nx_pylab.draw_networkx.html
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
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 so-called 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.
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.
Show code cell content
G = nx.Graph()
Show 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)),
]
Show code cell content
G.add_edges_from(edges)
Show 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 ofnetworkx
and the keyword argumentedge_color
. Plot the network on a figure with a size of 10x10 inches. Change thewidth
of the edges to the value 4.
Show code cell content
G = subgraphs[3] # iterate here
Show code cell content
len(G.nodes)
Show code cell content
len(G.edges)
Show code cell content
nx.is_planar(G)
Show code cell content
pd.Series({k: v for k, v in G.degree}).mean()
Show code cell content
len(nx.cycle_basis(G))
Show code cell content
cycle_length = pd.Series([len(c) for c in nx.cycle_basis(G)])
Show code cell content
cycle_length.plot.hist(bins=100)
Show code cell content
# alternative
pd.Series([len(c) for c in nx.cycle_basis(G)]).value_counts().sort_index().plot()
Show code cell content
cycle_length.describe()
Show code cell content
nx.adjacency_matrix(G).todense()
Show code cell content
K = nx.incidence_matrix(G).todense()
K
Show 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)
Show code cell content
H = B.dot(K.T)
Show code cell content
L = K.dot(H)
Show code cell content
L_inv = np.linalg.pinv(L)
Show code cell content
PTDF = H.dot(L_inv)
Show code cell content
imbalance = np.zeros(len(G.nodes))
imbalance[0] = 1
imbalance[-1] = -1
Show code cell content
flows = PTDF.dot(imbalance)
Show code cell content
abs_flows = np.abs(flows)
Show 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)