08: Tutorial: Using a random road graph model to understand road networks robustness to link failures
The purpose of this tutorial is to address some of the queries that I received about the model and data found in the following research paper:
Sohouenou, P. Y. R., Christidis, P., Christodoulou, A., Neves, L. A. C., & Lo Presti, D. (2020). Using a random road graph model to understand road networks robustness to link failures. International Journal of Critical Infrastructure Protection DOI: https://doi.org/10.1016/j.ijcip.2020.100353
Considering that the code was separated in several scripts and partly run on a High-Performance Computer (with its specific requirements), there is no point in sharing the raw code. Instead, I wrote this tutorial, which shows how to generate a GREREC network, acquire road network data, and perform a robustness analysis of these networks. The aim is to help anyone reading this to adapt the principles/code presented here to their needs and prefered programming language.
The tutorial is organised in three sections:
- The GREREC model (= random road network generator)
- The road network data
- The robustness analysis
As an indication, Section 1 and 3 were performed using an R 3.4.4 code, while Section 2 was performed using a Python 3.7 code.
Disclaimer: I’m neither a computer scientist nor a software engineer. Therefore, the code shown below is just a way to get the job done. I'm sure that it could be more organised and optimised.
NB: Please if you use this tutorial for your research don't forget to cite the article (above) ^^.
GREREC model (= random road network generator)¶
As stated in the paper, the GREREC model was developed to randomly generate a variety of abstract networks presenting the topological and operational characteristics of real-road networks. For details about the model characteristics and relevance to road networks, I direct you to the paper.
GREREC stands for Grid network with Random Edges and Regional Edge Costs. The procedure used to generate a graph with the GREREC model requires five steps:
The five steps/functions behind the generator¶
# Required libraries
library(igraph) # for network/graph analysis
library(dplyr) # for fast data manipulation
Step 1: Generate a graph with a rectangular grid topology¶
The dimensions m and n of the rectangle (i.e. the number of nodes per row and columns respectively) are the only parameters necessary to define the grid. The graph generated has nm nodes and the vertex on the i-th column and j-th row is denoted as $v_{i,j}$.
The code used to create this graph is shown in the cell below. I used a function that generates an adjacency matrix, which corresponds to the graph. As you will see in the code below, for the sake of simplicity, the nodes are given a label that is a number in [1,nm] depending on their position rather than the label $v_{i,j}$.
# Function to compute a Grid network of size n*m
# The nodes of the grid are called 1 -> n*m, from left to right and bottom to top
# The functon creates an adjacency matrix that can be transformed into a graph using an igraph function (see Section 1.2 Example)
Grid_adjacency <- function(m, n = NULL)
{
if (missing(n))
{
stop("n missing")
}
else
{
A = matrix(0, m * n, m * n)
for (i in 1:(m * n))
{
up = i - m
down = i + m
left = i - 1
right = i + 1
if (up > 0)
A[i, up] = 1
if (down <= (m * n))
A[i, down] = 1
if (left %% m != 0)
A[i, left] = 1
if (i %% m != 0)
A[i, right] = 1
}
}
A
}
Step 2: Check and remove the existing edges by the order “left to right, bottom to top” with probability (1-p)¶
The code used for this step is shown in the cell below. I wrote a function — Obstacle(Adj_mat, p) — that randomly removes some edges from the adjacency matrix of a grid graph, created with the previous function.
# Function to keep the horizontal and vertical edges with a probability p in [0,1]
# it modifies the adjacency matrix of the grid network (rk: the matrix is no longer symetric)
Obstacle <- function(Adj_mat,p) {
x <- nrow(Adj_mat)
if (p==1) { return(Adj_mat) # do nothing i.e. all the edges are kept
} else for(i in 1:x){
for(j in 1:x){
if (Adj_mat[i,j] == 1 && runif(1,0,1) < (1-p)) {Adj_mat[i,j]=0}
}
}
return(Adj_mat)
}
Step 3: For each vertex $v_{i,j}$ where both i and j are odd numbers, generate the four diagonal edges departing from $v_{i,j}$ with probability q.¶
The code used for this step is shown in the cell below. I wrote a function — Shortcuts_1(Adj_mat,m,n,q) — that randomly build shortcuts (=diagonals) in the matrix of the grid graph, created with the previous function. To develop this function I draw several examples on a paper. I encourage you to do the same to understand it.
# Function to build shortcuts with a probability q in [0,1]
Shortcuts_1 <- function(Adj_mat, m, n, q) # shortcuts departing from the nodes v(i,j) where both i and j are odd numbers
{
if (q==0) { return(Adj_mat) }
else
{
for(k in seq(1, n-1, by=2)) # select row 1, 3, 5 etc.(counting from 0)
{
for( i in seq( k*m + 1, (k+1)*m, by=1) ) # go trough the nodes in these rows
{
if ( (i-m) %% 2 == 0 ) # and select the even nodes that are in column 1, 3, etc.( => nodes for diag departure)
{
if (i %% m == 0) # if i on the right rim (only 2 diagonals possible)
{
diag_left_down = i - m - 1
diag_left_up = i + m - 1
if (diag_left_down <= (m * n) && runif(1,0,1) < q)
Adj_mat[i, diag_left_down] = 1
if (diag_left_up <= (m * n) && runif(1,0,1) < q )
Adj_mat[i, diag_left_up] = 1
} else {
diag_left_down = i - m - 1
diag_left_up = i + m - 1
diag_right_up = i + m + 1
diag_right_down = i - m + 1
if ( diag_left_down <= (m * n) && runif(1,0,1) < q)
Adj_mat[i, diag_left_down] = 1
if (diag_left_up <= (m * n) && runif(1,0,1) < q )
Adj_mat[i, diag_left_up] = 1
if ( diag_right_up <= (m * n) && runif(1,0,1) < q )
Adj_mat[i, diag_right_up] = 1
if ( diag_right_down <= (m * n) && runif(1,0,1) < q )
Adj_mat[i, diag_right_down] = 1
}
}
}
}
return(Adj_mat)
}
}
Step 4: For each vertex $v_{i,j}$ where i is an even number (regardless of j), generate the four diagonal edges departing from $v_{i,j}$ with probability q providing that the diagonal is not intersecting an existing one.¶
The function, Shortcuts_2(Adj_mat,m,n,q), used for this step is shown in the cell below. This function should be applied after Shortcuts_1(Adj_mat,m,n,q).
# Rk: Shortcuts_2 depends on the result of Shortcuts_1 -> builds a shortcut providing that the crossing diagonal does not exist
Shortcuts_2 <- function(Adj_mat, m, n, q) # shortcuts departing from the nodes v(i,j) where i is even and j is odd
{
if (q==0) { return(Adj_mat) }
else
{
for(k in seq(1, n-1, by=2)) # select row 1, 3, 5 etc.(counting from 0)
{
for( i in seq( k*m + 1, (k+1)*m, by=1) ) # go trough the nodes in these rows
{
if ( (i-m) %% 2 == 1 ) # and select the even nodes that are in column 0, 2, etc.( => nodes for diag departure)
{
if ( (i-1) %% m == 0) # if we are on the left rim (only 2 diagonals possible)
{
diag_right_up = i + m + 1
diag_right_down = i - m + 1
if ( diag_right_up <= (m * n) && Adj_mat[i+1, diag_right_up - 1 ] == 0 && runif(1,0,1) < q )
Adj_mat[i, diag_right_up] = 1
if ( diag_right_down <= (m * n) && Adj_mat[i+1, diag_right_down - 1] == 0 && runif(1,0,1) < q )
Adj_mat[i, diag_right_down] = 1
} else if ( i %% m == 0) #if we are on the right rim
{
diag_left_down = i - m - 1
diag_left_up = i + m - 1
if ( diag_left_down <= (m * n) && Adj_mat[ i-1, diag_left_down + 1] == 0 && runif(1,0,1) < q)
Adj_mat[i, diag_left_down] = 1
if (diag_left_up <= (m * n) && Adj_mat[i-1, diag_left_up + 1] == 0 && runif(1,0,1) < q )
Adj_mat[i, diag_left_up] = 1
} else {
diag_left_down = i - m - 1
diag_left_up = i + m - 1
diag_right_up = i + m + 1
diag_right_down = i - m + 1
if ( diag_left_down <= (m * n) && Adj_mat[ i-1, diag_left_down + 1] == 0 && runif(1,0,1) < q)
Adj_mat[i, diag_left_down] = 1
if (diag_left_up <= (m * n) && Adj_mat[i-1, diag_left_up + 1] == 0 && runif(1,0,1) < q )
Adj_mat[i, diag_left_up] = 1
if ( diag_right_up <= (m * n) && Adj_mat[i+1, diag_right_up - 1] == 0 && runif(1,0,1) < q )
Adj_mat[i, diag_right_up] = 1
if ( diag_right_down <= (m * n) && Adj_mat[i+1, diag_right_down - 1] == 0 && runif(1,0,1) < q )
Adj_mat[i, diag_right_down] = 1
}
}
}
}
return(Adj_mat)
}
}
Step 5: Randomly assign a travel cost to the edges by the order ”left to right, bottom to top” with the following rule: all the links departing from the same node have the same cost of travel.¶
The code used for this step is shown in the cell below. The function (n,max_value, min_value) generates an array of link cost values, that can be used to assign costs to the graph created with the previous functions.
# Function to generate link costs as normally distributed random numbers in the discrete interval [min_value, max_value]
Normal_link_costs <- function(n, max_value, min_value) # n is the number of values generated
{
x <- rnorm(n)
max_x <- max(x)
min_x <- min(x)
x <- (x - min_x) / (max_x - min_x)
x <- x * (max_value - min_value)
x <- x + min_value
x <- round(x)
return(x)
}
Example¶
The functions presented above can be combined to generate a GREREC graph. To this end, one need the choose some parameters to input into the model. I propose to use $n=m=4$, $p=0.9$ and $q=0.3$ in this example.
Graph generation¶
m <- 4
n <- 4
p <- 0.9
q <- 0.3
### Network generation ###
Matrix_grid <- Grid_adjacency(m, n)
Matrix_with_obstacles <- Obstacle(Matrix_grid, p)
Matrix_road_network <- Shortcuts_1(Matrix_with_obstacles, m, n, q)
Matrix_road_network <- Shortcuts_2(Matrix_road_network, m, n, q)
Matrix_road_network <- pmax(Matrix_road_network, t(Matrix_road_network)) # forces the matrix to be symetric
Matrix_road_network[ lower.tri(Matrix_road_network) ] <- 0 # puts zeros into the lower part of the matrix for an undirected analysis
## Link cost assignment
Grid_dim <- max(n,m)
N_nodes <- n*m
Link_costs <- Normal_link_costs (N_nodes, 1, Grid_dim) # normal distribution stretched between min (=1) and max (=Grid_dim)
Matrix_road_network <- Matrix_road_network * Link_costs
# Transform the adjacency matrix into a graph using an igraph function
G_road_network <- graph_from_adjacency_matrix(Matrix_road_network, mode = c("upper"), weighted = T)
# Check the connectedness of the graph generated. Remember that depending on the parameters the graph generated can be unconnected (see paper for the explanations)
if (is.connected(G_road_network) == F) {return("not connected")} # We only analysed networks that are connected in the paper !!
Can we visualise the graph?¶
Visualising the graph is rather tricky. I used the plot function of igraph (see the cell below). The problem is that the function does not always realize that the graph is planar and can be plotted in 2D. It usually works well for small networks ($n,m <5$) using layout = layout_with_fr(G_road_network, dim =2), OR layout = layout.grid, edge.width=E(G_road_network)$weight).
I obtained this result by testing several layouts. Based on my online research, I understood that the complexity of the problem exponentially increased with the graph size, hence no algorithm can currently automatically find that a large graph can be plotted in 2D — without coordinates of course.
plot(G_road_network, main = 'A GREREC network', vertex.label.color = "black", vertex.label.dist = 1.5,
vertex.size = 4, vertex.shape = 'circle', edge.curved = 0, frame = 'T', margin = c(0.05,0.05,0.05,0.05),
layout = layout_with_fr(G_road_network, dim=2), edge.width = E(G_road_network)$weight )
# #layout = layout.grid, edge.width=E(G_road_network)$weight) #other layout option
Ta da!
Real road network data¶
In the paper, we also analyse graphs extracted from road network datasets.
To acquire the samples, we used the Python package OSMnx (https://osmnx.readthedocs.io/en/stable/) to download drivable street network data within the chosen boundaries from OpenStreetMap(https://www.openstreetmap.org/) and automatically process them into length-weighted nonplanar graphs. I won't expand on this here as there is already a great tutorial on this (here: https://geoffboeing.com/2016/11/osmnx-python-street-networks/) written by the author of the package himself!
Robustness analysis¶
The purpose of the paper was to analyse the robustness of the GREREC and real networks to find patterns. In the paper, I considered 161 and 30 GREREC and real networks, respectively, using Quasi-monte Carlo simulations and high-performance computing to handle the complexity of the simulations. As there is little chance that you have access to the same computer configuration and try to tackle the same problem, there is no point in including this in the tutorial. This section of the tutorial will instead provide the code and explanations required to evaluate the robustness of the example (above) to one scenario only. This will be done in four steps.
The robustness indicator ($RO$) used in the paper measures the demand-weighted average increase in travel time (TT) along the Origin-Destination pairs of the network.
$RO = \sum \limits _{w} k_{w}*\Big({ 1 + \frac{TT_{d}^w - TT_{0}^w}{TT_{0}^w} \Big) ^{-1}}$
where $w$ and $k_w$ are an OD pair and the associated weighting factor, respectively. $k_w$ is the ratio between the travel demand on $w$ and the total travel demand. $TT_0^{w}$ and $TT_d^{w}$ are the undisrupted and disrupted travel times, respectively.
In the paper, the OD pairs were considered to be of equal importance (i.e. $k_w = 1/N_{OD}$, $N_{OD}$ being the number of OD pairs) and $TT^w$ is the travel cost on the least-path cost for $w$. The funtion below can be used to compute the main term of the sum, which measures the relative increase in travel time along an OD pair.
OD_TT_increase <- function(x,y) {return ( 1/( 1 + (y-x)/x) ) } # function to compute the increase in TT along an OD pair.
Step 1: Randomly define some origin/destination nodes¶
In the graphs representing road networks, some of the nodes do not serve as origin-destination points. To take this specificity into account, in the paper we assumed that the bottom-left ($v_{0,0}$) and top-right ($v_{mn,mn}$) nodes of the original grid (see the image above) were OD points and randomly selected additional OD points in the network with the probability $r$. When $r=0$ only these two nodes were considered as OD points, while when $r=1$ all the nodes in the network served as OD points.
The function below can be used to randomly select some nodes in a graph.
# Function to select certain nodes as Origin-Destination points with probability r in [0,1]
# If r = 0 only two nodes of the rim are selected (node 1 and node m*n), if r = 1 all OD pairs are selected
# the function returns a list of Origin nodes
ODpairs_selection <- function(r) {
Origins <- c(1,m*n)
for(i in 2:(m*n - 1) ){
if( runif(1,0,1) < r ) { Origins <- c(Origins, i) }
}
Origins <- sort(Origins)
return(Origins)
}
r = 0.8 # lets use r = 0.8
Origins <- ODpairs_selection(r)
Destinations <- Origins # origins are also destinations
Step 2: Evaluate the travel time in the undisrupted network.¶
The travel times bewteen the origins and destinations can be computed using the function distances() of igraph and stored in a dataframe, as shown below.
shortestpaths <- distances(G_road_network, v = Origins, to = Destinations, algorithm = "dijkstra") # Igraph function to calculate the length of the shortest paths
dimnames(shortestpaths) <- list(O = Origins, D = Destinations)
OD_data <- as.data.frame(as.table(shortestpaths)) # the data is stored in a data frame
OD_data <- OD_data[!(OD_data$O == OD_data$D),] # removes the shortest path from O to O
names(OD_data)[3] <- "TT_0"
Step 3: Select a link to delete¶
Two methods can be used to delete a link in a graph:
- Delete a link from the original graph using the function delete_edges() of igraph
- Write a 0 in one of the non-zero entries of the adjacency matrix of the original graph and build a new graph.
I tested both and found that option 2 was faster. Therefore, I used option 2.
In this example, we will delete the link between node 1 and 2 (see the network plotted above). You can print the matrix using the print() function to see its non-zero entries, which are the links of the graph.
Matrix_road_network_d <- Matrix_road_network
Matrix_road_network_d[1,2] <- 0
G_road_network_d <- graph_from_adjacency_matrix(Matrix_road_network_d, mode = c("upper"), weighted = T )
Step 4: Compute the robustness indicator¶
To this end, the travel time along the OD pairs of the disrupted network should be computed and used to compute the network robustness indicator($RO$).
shortestpaths_d <- distances(G_road_network_d, v = Origins, to = Destinations, algorithm = "dijkstra")
dimnames(shortestpaths_d) <- list(O = Origins, D = Destinations)
OD_data_d <- as.data.frame(as.table(shortestpaths_d))
OD_data_d <- OD_data_d[!(OD_data_d$O == OD_data_d$D),] # deletes the shortest path that starts and arrives to the same node
OD_data_d <- bind_cols(OD_data, OD_data_d[3])
OD_data_d[5] <- mapply(OD_TT_increase, OD_data_d[3], OD_data_d[4]) # calculates relative increase in TT along each OD pairs
Network_robustness <- colMeans(OD_data_d[5]) # calculates the network robustness (equal weights for all OD pairs)
print(Network_robustness)
V5 0.9974359
The result should be Network_robustness = 0.974.