A greedy procedure never looks back, which makes it fast and usually wrong. Dijkstra’s procedure for shortest paths is the rare one where never looking back is exactly right.
New Julia packages used
HTTP is Julia’s package for the HTTP protocol that web pages and web services use. As a client it sends a request, such as getting a web address, and returns the server’s response; it can also run a web server, and it supports HTTP/2, server-sent events and WebSockets.
JSON3 reads and writes JSON, the text format in which most web services answer. Reading a JSON document gives an object whose fields are reached with a dot or brackets, like a Julia struct or dictionary, and Julia values can be written back out as JSON.
New Logjam functions used
addconnectors: Add connector links from demand points to a road network.
cropnetwork: Crop a road network to a bounding box defined by coordinates.
prune_reindex: Prune network to common vertices and reindex node IDs sequentially.
segcost: Calculate the cost of each segment in a location sequence.
shortestpaths: Compute shortest path distances and parent pointers from the first n nodes.
thin: Remove degree-2 nodes from a road network, merging their incident links.
tracepath: Reconstruct a shortest-path node sequence from a Dijkstra parent vector.
x2ln: Convert graph edges to NaN-separated line coordinates for plotting.
For many network applications in logistics, the network is really just a way of talking about the problem: underneath, the problem is a math program, and nothing about its being a network is used beyond that way of thinking about it. The assignment and transportation problems of Lecture 4.2, among them assigning a ride-hailing service’s cars to passengers wanting a ride, are of this kind. The shortest path problem is different: it is the one problem in this course that really is a network problem.
One of its main applications is to take an actual road network and find the shortest path through it. Lecture 2.2, Sec. 7, found great circle distances to be fairly accurate as long as the distances are not within about 50 miles. Inside a metro area they are probably not the best, especially where a major river runs through it with only two or three bridges, which become a bottleneck: a point one mile away on the other side of the river may be three miles from the nearest bridge. Over such short distances the disparity can be large, and the same holds wherever waterways, islands or canals shape the roads. A road network brings those features into the problem, which is one good reason to use it. The other is that what is usually wanted is not the shortest distance but the shortest travel time. The route Google Maps recommends first is almost always the one with the shortest travel time, not necessarily the shortest distance. Distance is used as the weight on each arc until Sec. 5, which compares using travel time as the weight with using distance.
Great circle distances are not accurate over shorter (< 50 mi) distances, or in areas with waterways that restrict road travel, and there the actual distance using the road network is more accurate. And a traveler will usually want to find the shortest travel time path, not the shortest distance path.
The terms network and graph can be used interchangeably. The terms graph/vertex/edge are typically used when the focus is on properties of the mathematical object, and the terms network/node/arc are used when the focus is on the properties of a real-world system. The terms graph and digraph (short for directed graph) are used in mathematics to describe networks in which the arcs are either all undirected or directed, respectively. In many network applications, both directed and undirected arcs are used (e.g., one- and two-way roads), and the arcs have associated weights that represent, for example, the distance or travel time of the road being represented. In road network data, such as the Oak Ridge and FAF5 networks of Sec. 4, an arc is called a link: a road segment between two nodes, one-way or two-way, so that a two-way link becomes two directed arcs in a digraph. Table 1 sets the three sets of terms side by side.
Table 1: Network, graph and digraph terms.
Network
Graph
Digraph
Node
Vertex
Vertex
Arc/link (undirected/directed)
Edge (undirected)
Edge (directed)
In Julia, both graphs and digraphs can be created. If a network with a mix of directed and undirected arcs is needed, then a digraph can be used, with two directed arcs in opposite directions used to represent an undirected arc.
There are different network representations that make it easier to operate on the network for specific applications, and it is usually possible to translate from one representation to another. Given a network with n nodes and m arcs:
Arc list:m-element array of the source nodes, m-element array of destination nodes, and an m-element array of the weights of each arc in the network.
Adjacency matrix:n \times n matrix of arc weights, with weights of 0 indicating no arc; each undirected arc between nodes i and j is represented both as element (i, j) and as element (j, i), so the adjacency matrix of a graph is symmetric.
Show the code that draws this figure
# Code block 1: the three-node digraph, drawn from its arc listXY3 = [0.00.0; 2.00.0; 1.01.5] # node positionsr3 =0.22# node radius, data unitsc3 =vec(sum(XY3, dims =1)) ./3# the centroid, to place labelsfig =Figure(size = (300, 250))ax =Axis(fig[1, 1]; aspect =DataAspect())hidedecorations!(ax); hidespines!(ax)for (i, j, w) inzip([1, 3, 3], [2, 1, 2], [1, 3, 2]) p, q = XY3[i, :], XY3[j, :] u = (q .- p) ./hypot((q .- p)...) t, h = p .+ r3 .* u, q .- (r3 +0.015) .* u # outline to outlinearrows2d!(ax, [t[1]], [t[2]], [h[1] - t[1]], [h[2] - t[2]]; color = colorant"#252525", shaftwidth =1.5, tipwidth =9, tiplength =9) m = (p .+ q) ./2 m = m .+0.2.* (m .- c3) ./hypot((m .- c3)...) # outside the arctext!(ax, m[1], m[2]; text =string(w), fontsize =14, align = (:center, :center), color = colorant"#252525")endfor k in1:3poly!(ax, Circle(Point2f(XY3[k, 1], XY3[k, 2]), r3); color =:white, strokecolor = colorant"#252525", strokewidth =1.8)text!(ax, XY3[k, 1], XY3[k, 2]; text =string(k), fontsize =15, align = (:center, :center), color = colorant"#252525")endxlims!(ax, -0.4, 2.4); ylims!(ax, -0.45, 1.85)fig
Figure 1: A digraph of three nodes and three arcs, each arc labeled with its weight.
For example, the digraph of Fig. 1 has the arc list and the adjacency matrix W of Eq. 1:
The arc list has one row for each arc, and W has one row and one column for each node.
Complete bipartite digraph
A complete bipartite directed graph (or digraph), as in Fig. 2, can represent suppliers shipping to multiple distribution centers (DCs) by a single mode of transport.
Show the code that draws this figure
# Code block 2: a network of two groups, suppliers 1-3 and DCs 4-5supc = colorant"#d21f26"# suppliersdcc = colorant"#1e8a3c"# DCsink = colorant"#252525"r =0.24# node radius, data unitsXY = [0.02.0; 0.01.0; 0.00.0; 2.41.6; 2.40.4] # node positions# One arc from node i to node j, outline to outline: kind :one (an arrow),# :two (a two-way arc, no arrowhead) or :absent (dashed); off shifts it# sideways, for a pair of parallel arcsfunctionarc!(ax, i, j, lab; kind =:one, off =0.0, lpos =0.5, lside =1) p, q = XY[i, :], XY[j, :] u = (q .- p) ./hypot((q .- p)...) nrm = [-u[2], u[1]] # unit normal s =asin(clamp(off / r, -1, 1)) # where a shifted arc meets t = p .+ r *cos(s) .* u .+ off .* nrm # the circle h = q .- (r *cos(s) +0.015) .* u .+ off .* nrmif kind ==:absentlines!(ax, [t[1], h[1]], [t[2], h[2]]; color = ink, linewidth =1, linestyle =:dash)elseif kind ==:twolines!(ax, [t[1], h[1] +0.015u[1]], [t[2], h[2] +0.015u[2]]; color = ink, linewidth =1.5)elsearrows2d!(ax, [t[1]], [t[2]], [h[1] - t[1]], [h[2] - t[2]]; color = ink, shaftwidth =1.5, tipwidth =9, tiplength =9)endif !isempty(lab) m = t .+ lpos .* (h .- t)text!(ax, m[1], m[2]; text = lab, color = ink, fontsize =14, align = (:center, :center), offset = (12lside * nrm[1], 12lside * nrm[2]))endreturnnothingendfunctiontwogroup(arcs; w =340, h =300) fig =Figure(size = (w, h)) ax =Axis(fig[1, 1]; aspect =DataAspect())hidedecorations!(ax); hidespines!(ax)for (i, j, lab, kw) in arcsarc!(ax, i, j, lab; kw...)endfor k in1:size(XY, 1) # the nodes, over the arcs x, y = XY[k, :] c = k <=3 ? supc : dccpoly!(ax, Circle(Point2f(x, y), r); color =:white, strokecolor = c, strokewidth =1.8)text!(ax, x, y; text =string(k), color = c, fontsize =15, align = (:center, :center))endxlims!(ax, -0.4, 2.8); ylims!(ax, -0.4, 2.4)return figendtwogroup([(1, 4, "6", (lpos =0.25,)), (1, 5, "10", (lpos =0.2, lside =-1)), (2, 4, "0", (lpos =0.3,)), (2, 5, "13", (lpos =0.3, lside =-1)), (3, 4, "9", (lpos =0.2,)), (3, 5, "16", (lpos =0.3, lside =-1))])
Figure 2: A complete bipartite digraph: three suppliers (red) each ship to both of two DCs (green), one arc per pair, each labeled with its weight.
The weighted adjacency matrix of Fig. 2 has a row and a column for each of its five nodes:
The function SimpleWeightedDiGraph is used to create a directed network using an adjacency matrix, where each element (i, j) in the matrix represents an arc directed from a starting node (i) to an ending node (j):
# Code block 4: a digraph from the adjacency matrixusingGraphs, SimpleWeightedGraphsg =SimpleWeightedDiGraph(W) # Create digraph from adjacency matrixprintln.(edges(g));
Edge 1 => 4 with weight 6.0
Edge 1 => 5 with weight 10.0
Edge 2 => 5 with weight 13.0
Edge 3 => 4 with weight 9.0
Edge 3 => 5 with weight 16.0
The problem is that arc 2 => 4, of 0 weight, is not included.
The network of Fig. 2 is a common one in logistics. If nodes 1 to 3 were DCs and nodes 4 and 5 customers, it would probably not be a distribution network, since a distribution network would usually allocate each customer to one DC, and here every DC serves both customers. It is more likely two manufacturing plants that both draw on the same three major suppliers of the components that go into their product. Just looking at the network view of a problem gives some insight into what is possible, and the weights could represent, for example, the cost of supplying along each arc.
The zero on arc 2 → 4 is the one thing to note. It does not mean that the arc is unusable. If the weight is a distance, a zero means that supplier 2 is located at plant 4, in the same building or next door, so in many logistics applications the zero is relevant and the arc exists. A lot of software, though, has assumed that a zero means the absence of an arc, because it represents a graph as a sparse matrix, and in a sparse matrix a zero is the absence of an entry. Math solvers such as CPLEX and Gurobi have done the same: given an arc whose cost is zero, a solver could eliminate that variable from the model as an arc not wanted, and in a big model the cause can take a day or more to find. The Julia package SimpleWeightedGraphs, which stores a graph as a sparse matrix, puts the problem front and center in its documentation. A distance matrix with zeros works fine for a calculation that never builds a graph; the problem arises only once the matrix becomes a weighted adjacency matrix, stored as a sparse array.
In SimpleWeightedGraphs, an arc of zero weight is lost when the graph is built from a full matrix, as in Code block 4, when the arc is added with add_edge!, and when arithmetic is done on the graph’s sparse weight matrix; the package’s documentation suggests a very small weight in its place.1 A graph built in one go from its arc list keeps an arc of zero weight as it is, with no small value standing in for it, and every network in the rest of this lecture is built that way:
# Code block 5: a digraph from its arc listi = [1, 1, 2, 2, 3, 3] # source node of each arcj = [4, 5, 4, 5, 4, 5] # destination node of each arcw = [6., 10., 0., 13., 9., 16.] # weight of each arcg =SimpleWeightedDiGraph(i, j, w) # Create digraph using arc listsprintln.(edges(g));
Edge 1 => 4 with weight 6.0
Edge 1 => 5 with weight 10.0
Edge 2 => 4 with weight 0.0
Edge 2 => 5 with weight 13.0
Edge 3 => 4 with weight 9.0
Edge 3 => 5 with weight 16.0
# Code block 6: the adjacency matrix, a view of the digraph@showhas_edge(g, 2, 4) # Arc 2 => 4 is in the digraphMatrix(adjacency_matrix(g))
1
Since the adjacency matrix is stored as a sparse matrix, convert it to a regular (dense) matrix for display purposes.
The digraph has arc 2 => 4, yet its adjacency matrix shows a 0 in that position, the same as where there is no arc at all: the arc list carries what the matrix cannot.
Bipartite graph
A bipartite graph has one- or two-way connections between nodes in two groups, as in Fig. 3, and its arcs can be listed in an arc list matrix, one row per arc.
Figure 3: A bipartite graph with one- and two-way connections: arc 1–4 runs both ways at weight 6, 3 → 5 and 5 → 3 have different weights, and there is no arc between 3 and 4 (dashed).
Fig. 3 is a variation on Fig. 2. It is still bipartite, with two different sets of nodes, but it is not complete: the dashed line means there is no arc from 3 to 4. Some of its arcs also run in both directions, so it combines directed and undirected arcs. Its weighted adjacency matrix shows a zero for the missing arc from 3 to 4, and here that zero really should be a zero, meaning no arc.
An alternative to the weighted adjacency matrix is the arc list matrix, which simply lists the source and destination of each arc along with its weight; the same information can be kept as three separate vectors of sources, destinations and weights. Its first column is the source of each arc, its second the destination and its third the cost on the arc. The arcs from 3 to 5 and from 5 to 3 are two different arcs with different values, so each has its own row. The undirected arc between 1 and 4 appears twice, once from 4 to 1 and once from 1 to 4, both with value 6. That duplication is the big representational difference between undirected and directed: Julia, like most graph software, builds a network as one or the other, and an undirected arc in a directed graph is represented by duplicating it in both directions with the same weight. It matters for road networks, which have both two-way and one-way roads.
# Code block 8: the bipartite graph, a digraph from its arc listi = [4, 5, 1, 2, 1, 2, 3] # arc 1-4 listed both waysj = [1, 3, 4, 4, 5, 5, 5]w = [6., 3., 6., 0., 10., 13., 16.]g =SimpleWeightedDiGraph(i, j, w)println.(edges(g));
Edge 1 => 4 with weight 6.0
Edge 1 => 5 with weight 10.0
Edge 2 => 4 with weight 0.0
Edge 2 => 5 with weight 13.0
Edge 3 => 5 with weight 16.0
Edge 4 => 1 with weight 6.0
Edge 5 => 3 with weight 3.0
A network whose arcs are all two-way can instead be built as a graph. Taking the arcs of Code block 5 as an example, the only change is that each arc in the network is undirected from the starting node (i) to the ending node (j). The function SimpleWeightedGraph is used to create an undirected network, as opposed to using SimpleWeightedDiGraph to create a directed network:
# Code block 9: the arcs of Code block 5, undirectedg =SimpleWeightedGraph([1, 1, 2, 2, 3, 3], [4, 5, 4, 5, 4, 5], [6., 10., 0., 13., 9., 16.]) # Create graphMatrix(adjacency_matrix(g)) # Both directions represented for each arc
A multigraph has multiple connections between the same pair of nodes, such as multiple modes of transport, as in Fig. 4. A simple graph does not have multiple connections.
Figure 4: A multigraph: two two-way arcs join 1 and 4, at weights 6 and 18, and two arcs run from 3 to 5, at weights 16 and 3, as two modes of transport might.
The package used for these networks is called SimpleWeightedGraphs, and the term simple comes up quite a bit. It is there because graph theory has another concept, the multigraph, in which two nodes can be joined by more than one arc. In Fig. 4, nodes 1 and 4 are joined by two undirected arcs with different weights, and two directed arcs run from 3 to 5, again with different weights. An adjacency matrix cannot hold two values for one pair of nodes: it holds one, so a choice has to be made of which to put there. A multigraph therefore cannot be represented by a weighted adjacency matrix, and the only graphs that can be are simple graphs, which do not have multiple connections.
By default the parallel arcs are added together: 6 + 18 = 24 between 1 and 4, and 16 + 3 = 19 from 3 to 5, weights that belong to no arc in the network. When parallel arcs are reduced to one, the shorter is the one to keep, and that takes a single keyword:
# Code block 12: keep the shorter of each pair of parallel arcsMatrix(weights(SimpleWeightedDiGraph(i, j, w; combine = min)))
Sec. 4 returns to this rule, since a road network often has two links between the same pair of nodes.
The parallel arcs are summed by default because a sparse matrix adds up whatever is entered more than once for the same element, and for most logistics applications a sum does not represent what actually happens. There are contexts, though, where a multigraph is what is wanted. Of the two arcs between nodes 1 and 4, the 18 might be a costlier mode of transport and the 6 a cheaper one with limited capacity, so that only so much flow can use the cheaper mode and beyond that the costlier one is needed. Representing that takes lower and upper bounds on the arcs, which the minimum cost network flow of Lecture 4.2 provides, and it can only be done with an arc list, never with an adjacency matrix. That is one of the reasons for having several representations of a network: some cannot represent everything, and in other cases one is just more convenient than another, as the adjacency matrix often is even though it cannot represent everything.
2. Dijkstra’s procedure
Dijkstra’s procedure is used to find the shortest path in a network. It is a greedy procedure that can quickly find the (global) optimal solution.
A small network can be solved almost exhaustively by looking at it, but on larger networks solving by hand does not go well, and a systematic procedure is needed. Dijkstra’s shortest path procedure is greedy in the sense that it always makes whatever the next best choice is. A greedy procedure usually gives a reasonable solution but not an optimal one. Dijkstra’s procedure has the nice feature of being guaranteed optimal, so wherever it applies it should be used. It was developed in 1959, which is surprisingly recent for so basic a problem. Other shortest path procedures existed before it, but they were computationally more difficult, and this one is about the easiest to implement.
It is already implemented, so the reason to see how it works is not to implement it but to see another kind of greedy procedure, after the ADD and DROP procedures of Lecture 2.4 and before the ones in the Routing lectures, Lectures 5.1 to 5.3. Exposure to their logic gives a feel for what is computationally feasible in a logistics problem. When there is no proper way to solve a problem, as a linear program or otherwise, some kind of greedy procedure can quite often be devised. Its answer is not likely to be optimal, but it is better than having no answer at all, and many companies would be happy just to have a feasible solution. With the trucks, drivers and customers they have, some cannot find even a feasible one, never mind an optimal one.
The Dijkstra procedure requires that arcs have all nonnegative lengths. It is a greedy label setting algorithm since a step to the final solution is made as each node is labeled.
The shortest path problem is stated as a model in Model 1, and its formulation is the procedure itself, written as pseudocode.
minimize: total length of a path from an origin node to a destination node
solve for:
(a) path, a sequence of arcs from the origin to the destination, each starting at the node where the previous one ended.
subject to:
(a) network: the path uses only arcs of the network, each in its direction.
assumptions:
(a) every arc length is nonnegative;
(b) length of a path is the sum of the lengths of its arcs.
Model 1: Shortest path
Model 1 formulation: Shortest path
algorithm dijkstra;{ input: n × n matrix W of arc lengths, ∞ where no arc; origin s }{ output: distance d(i) and predecessor pred(i) of each node i }begin S := ∅; { no node labeled yet } d(i) := ∞ for each i ∈ {1, …, n}; { no node reached yet } d(s) := 0, pred(s) := 0; { s reached, from no node } while |S| < n do { until every node is labeled } begin i := arg min{ d(j) : j ∈ S̄ }; { nearest node not yet in S } S := S ∪ {i}; { label i: d(i) is now final } for each j with Wᵢⱼ < ∞ do { each node an arc from i reaches } if d(j) > d(i) + Wᵢⱼ then { shorter by way of i? } d(j) := d(i) + Wᵢⱼ, pred(j) := i; { then update j } end; return d, pred;end;
Neither the pseudocode nor its Julia implementation below needs to be changed or adapted, and the implementation can be treated as a black box. The one thing to take from them is that the whole procedure is not hundreds of lines of code. It is easy to implement in almost any programming language, even as a macro in a spreadsheet, which has no shortest path procedure of its own.
The critical step is line 10 of the pseudocode, which finds, among all of the unlabeled nodes, the one at the minimum distance. The \arg\min says only to find the minimum; it does not say what data structure to use to find it, but the choice makes a difference to the complexity. With n nodes, looking through every unlabeled node every time makes the procedure O(n^2). Routing from Raleigh to Los Angeles that way would probably take even a fast phone a minute or so, whereas keeping the unlabeled nodes in a priority queue, or in the more specialized Fibonacci heap, takes a couple of seconds, which is a big difference in performance. The procedure is so well known that an AI assistant, asked to write it and then to make it more efficient, will add a priority queue at once. Sec. 3 compares the time complexity of the shortest path procedures.
The O(\cdot) notation, read “order of,” states how the work a procedure does grows with the size of its problem, here the number of nodes n and of arcs m, leaving out constant factors and smaller terms. A procedure that is O(n^2) does work that, for large n, grows no faster than a constant times n^2, so doubling the number of nodes roughly quadruples the work. With the constants left out, the notation compares procedures rather than computers or programs: on a road network, where m is a small multiple of n, an O(m \log n) procedure eventually does less work than an O(n^2) one, whatever the constants. A procedure whose work grows as a power of the size, such as n^2 or n^4, is called polynomial, and one whose work can grow as 2^n is called exponential, which in time outgrows any power.2
Model 1 implementation: Shortest path
# Model: shortest path, by Dijkstra's procedurefunctiondijkstra(W, s) n =size(W, 1) # Number of nodes d =fill(Inf, n) # Distance from s pred =zeros(Int, n) # Predecessor of each node S =falses(n) # Labeled nodes: none yet d[s] =0.0# Distance from s to s = 0whilecount(S) < n# Find the unlabeled node at the minimum distance (SLOW step!) i =0for j in1:nif !S[j] && (i ==0|| d[j] < d[i]) i = jendend S[i] =true# Label node i# Update each node that an arc from i reaches more cheaplyfor j in1:nif d[i] + W[i, j] < d[j] d[j] = d[i] + W[i, j] pred[j] = iendendendreturn (d = d, pred = pred)end
1
The distance to each node from s, infinite until the node is reached. Infinity sometimes causes trouble in code, and the sum of all the arc lengths plus one then serves as infinity just as well.
2
i == 0 means that no unlabeled node has been found yet in this pass, so the first one met becomes the candidate, and each later one replaces it only if its distance is smaller.
dijkstra (generic function with 1 method)
The implementation is pretty easy to follow, and should a modification to the procedure ever be needed, it is not difficult to get into the details.
How the minimum is found in lines 10 to 15 is what separates a real implementation of Dijkstra’s procedure from a pretend one. Looking through every unlabeled node makes the procedure O(n^2): with 10 nodes the work is on the order of 100, and with 1,000 nodes on the order of a million, which becomes time consuming. A route planner that takes three or four minutes to route a ten-mile trip is very likely doing exactly that. Keeping the unlabeled nodes in a priority queue, such as a Fibonacci heap, makes finding the minimum much faster, which is what lets a route from Raleigh to California come back at once. A shortest path on a network of a thousand or more nodes is a good benchmark of whether an implementation is efficient.
Timed on one computer, on the whole FAF5 road network of Sec. 4, the heap-based procedure of Graphs.jl finds the shortest paths from one node to all of its nearly 345,000 nodes in about 0.05 s, and MATLAB’s built-in procedure in about 0.03 s. On a 2,476-node piece of the same network, the linear search of the implementation of Model 1 takes about 0.03 s, and the same loop in MATLAB about 0.04 s. Written instead with whole-array operations, as MATLAB code usually is, the procedure takes about 0.05 s in MATLAB but about 0.13 s in Julia, which builds new arrays at every step: Julia is fast when a loop is written out element by element, as it is here. The heap-based procedure handles the whole network in about twice the time the linear search takes on a piece about 140 times smaller.
Example 1: Six-node network
Determine the shortest path from node 1 to node 6 in the network of Fig. 5, by hand, then in Julia, and then with its arcs directed.
Show the code that draws this figure
# Code block 13: the six-node networkredc = colorant"#d21f26"# labels and labeled nodesP6 = [0.00.0; 1.00.58; 1.0-0.58; 2.00.58; 2.0-0.58; 3.00.0]A6 = [124; 132; 231; 245; 348; 3510; 452; 466; 563]r6 =0.15# node radius, data unitsside = [1, -1, -1, 1, 1, 1, -1, 1, -1] # side of each arc for its lengthunder = [true, false, true, false, true, true] # label below the node?# The network on ax: node k labeled labels[k], the nodes in done circled# in red, the arcs of path drawn in red, arrows if directedfunctionnet6!(ax; labels =fill("", size(P6, 1)), done =Int[], path =Int[], directed =false, fs =13)onpath(i, j) =any(k -> (path[k], path[k+1]) == (i, j) || (!directed && (path[k], path[k+1]) == (j, i)),1:length(path)-1)for (k, (i, j, w)) inenumerate(eachrow(A6)) p, q = P6[i, :], P6[j, :] u = (q .- p) ./hypot((q .- p)...) t, h = p .+ r6 .* u, q .- r6 .* u # outline to outline c, lw =onpath(i, j) ? (redc, 3.0) : (ink, 1.2)if directedarrows2d!(ax, [t[1]], [t[2]], [h[1] - t[1]], [h[2] - t[2]]; color = c, shaftwidth = lw, tipwidth =9, tiplength =9)elselines!(ax, [t[1], h[1]], [t[2], h[2]]; color = c, linewidth = lw)end m = (p .+ q) ./2.+0.11* side[k] .* [-u[2], u[1]]text!(ax, m[1], m[2]; text =string(w), color = ink, fontsize = fs, align = (:center, :center))endfor k in1:size(P6, 1) x, y = P6[k, :]poly!(ax, Circle(Point2f(x, y), r6); color =:white, strokecolor = k in done ? redc : ink, strokewidth = k in done ? 2.5:1.2)text!(ax, x, y; text =string(k), color = ink, fontsize = fs, align = (:center, :center))isempty(labels[k]) &&continue dy = under[k] ? -(r6 +0.1) : r6 +0.1text!(ax, x, y + dy; text = labels[k], color = redc, fontsize = fs, font =:bold, align = (:center, :center))endxlims!(ax, -0.35, 3.35); ylims!(ax, -0.95, 0.95)returnnothingendfig =Figure(size = (420, 240))ax =Axis(fig[1, 1]; aspect =DataAspect())hidedecorations!(ax); hidespines!(ax)net6!(ax)text!(ax, [-0.25, 3.25], [0.0, 0.0]; text = [L"s", L"t"], color = redc, fontsize =17, align = (:center, :center))fig
Figure 5: The six-node network, each arc labeled with its length and traveled in either direction; s marks the origin and t the destination.
Example 1(a): By hand
Determine the shortest path from node 1 to node 6, and its length, by labeling the nodes one at a time.
Dijkstra’s procedure is not a heuristic: it gives the optimal path. Before applying it, though, it is always a good idea to calibrate, providing a baseline. Taking the nodes as cities and each arc as a road that can be traveled in both directions, the question is the shortest total distance from city 1 to city 6. The network is about as big as can be solved by looking at it, and keeping a running total along the likely paths turns up 14 and then 13. Nothing better than 13 turns up, and 13 turns out to be the minimum. An algorithm does not usually work the way the eye does, though.
Dijkstra’s procedure first marks the starting node, node 1, and the destination or target node, node 6, as s and t in Fig. 5. It keeps a running total of the current cost from node 1 to every other node, and initially every one is set to infinity, even node 1 itself. A node circled in red has been labeled: the cost of reaching it is known, and it never needs to be looked at again. Each label has two arguments, the distance first and then the predecessor, the node it was reached from, which is what allows the path to it to be reconstructed. Node 1 is labeled 0,0, a distance of 0 and a predecessor of 0, since it is reached from no other node.
The procedure then loops, as Fig. 6 shows step by step. From the labeled nodes, at first only node 1, it looks at the nodes reachable by a direct arc, nodes 2 and 3, and determines for each the minimum cost of reaching it. Node 3 can be reached from node 1 at a cost of 2, so its infinity is replaced by 2,1, and node 2 at a cost of 4, so its infinity becomes 4,1. The next step is the computationally most intensive one, and it is what makes Dijkstra’s procedure slow without the right data structure: among all of the unlabeled nodes, at costs of 2, 4 and infinity, pick the minimum. Looking through every node every time is the bottleneck, and priority queues, Fibonacci heaps and other data structures do it much more efficiently. Node 3, at 2, is labeled, and its cost can never get lower.
The procedure repeats from node 3. Node 2 can now be reached at a cost of 3, the 2 to reach node 3 plus the arc of 1, which beats the 4 coming from node 1, so its 4,1 is replaced by 3,3. Nodes 4 and 5 are reached for the first time, at 10,3 and 12,3. The cost is always cumulative from node 1, and the predecessor says where the path came from. That is the key idea of Dijkstra’s procedure: since the nodes are labeled one at a time and never revisited, each arc only needs to be looked at once. Node 2, at 3, is now the lowest unlabeled node and is labeled next, which changes node 4 to 8,2, three plus five. Node 4, at 8, is labeled next, changing node 5 to 10,4 and reaching node 6 for the first time, at 14,4. The target has been reached, at a cost of 14, but the procedure is not finished: node 5, at 10, is labeled next, and node 6 then becomes 13,5, which is where the 13 comes from.
The predecessors then retrace the path, starting at the end. Node 6 was reached from node 5, node 5 from node 4, node 4 from node 2, node 2 from node 3, and node 3 from node 1, which is the way back to the starting point.
Committing to the 13 found by eye before the procedure runs is the Prior check of Lecture 1.2, and the procedure confirms it.
Show the code that draws this figure
# Code block 14: the loop of dijkstra, keeping the labels after each stepfunctiondijkstra_steps(W, s) n =size(W, 1) d =fill(Inf, n); pred =zeros(Int, n); S =falses(n); d[s] =0.0 steps = []whilecount(S) < n i =0for j in1:nif !S[j] && (i ==0|| d[j] < d[i]) i = jendend S[i] =truefor j in1:nif d[i] + W[i, j] < d[j] d[j] = d[i] + W[i, j] pred[j] = iendendpush!(steps, (i = i, d =copy(d), pred =copy(pred)))endreturn stepsendn6 =size(P6, 1) # number of nodesW6 =fill(Inf, n6, n6) # arc lengths, both directionsfor (i, j, w) ineachrow(A6) W6[i, j] = W6[j, i] = wendsteps =dijkstra_steps(W6, 1)lbl(d, p) =isinf(d) ? "∞":string(round(Int, d), ",", p)p6 = [6] # the path, traced back from 6while p6[1] !=1pushfirst!(p6, steps[end].pred[p6[1]])endfig =Figure(size = (687, 620))for (k, st) inenumerate(steps) axk =Axis(fig[(k -1) ÷2+1, (k -1) %2+1]; aspect =DataAspect(), title ="Step $k: label node $(st.i)", titlesize =13, titlefont =:regular)hidedecorations!(axk); hidespines!(axk)net6!(axk; labels = [lbl(st.d[v], st.pred[v]) for v in1:n6], done = [s.i for s in steps[1:k]], path = k ==length(steps) ? p6 :Int[], fs =12)endcolgap!(fig.layout, 6); rowgap!(fig.layout, 4)fig
Figure 6: Dijkstra’s procedure on the six-node network, one frame per step. Each step labels the unlabeled node at the smallest distance, circling it in red, and updates the nodes it reaches; a label reads distance, predecessor, and ∞ marks a node not yet reached. The last frame draws the shortest path in red.
Table 2 sets out the same labels as a table, one row per step, which is the form a hand solution takes.
Show the code that builds the table
# Code block 15: the labels after each step, node by nodeprintln("| Step | Labeled | ", # <1>join(["Node $v" for v in1:n6], " | "), " |")println("|:-:|:-:|", repeat(":-:|", n6))for (k, st) inenumerate(steps) c = [lbl(st.d[v], st.pred[v]) for v in1:n6] c[st.i] ="**"* c[st.i] *"**"# <2>println("| ", k, " | ", st.i, " | ", join(c, " | "), " |")endprintln("\n: Labels (distance, predecessor) after each step of ", # <3>"Dijkstra's procedure from node 1, the node labeled at that ","step in bold. {#tbl-dijkstra-steps .fit}\n")
Table 2: Labels (distance, predecessor) after each step of Dijkstra’s procedure from node 1, the node labeled at that step in bold.
Step
Labeled
Node 1
Node 2
Node 3
Node 4
Node 5
Node 6
1
1
0,0
4,1
2,1
∞
∞
∞
2
3
0,0
3,3
2,1
10,3
12,3
∞
3
2
0,0
3,3
2,1
8,2
12,3
∞
4
4
0,0
3,3
2,1
8,2
10,4
14,4
5
5
0,0
3,3
2,1
8,2
10,4
13,5
6
6
0,0
3,3
2,1
8,2
10,4
13,5
The block prints a Markdown table, one println for each row, and the page renders what it prints as a table, since the cell’s output is marked asis.
The label of the node labeled at this step is wrapped in **, Markdown’s bold.
A line that opens with a colon after a table is the table’s caption in Markdown, and the {#tbl-...} at its end gives the table its number and the anchor the text refers to.
The last row of Table 2 gives the predecessor of each node, from node 1 to node 6:
Table 2 displays all of the intermediate labeling steps, as the loop of the implementation computes them: the node labeled at each step, the current distance estimate of every node, and the predecessors. Initially every distance is infinity except at the starting node, and every predecessor is 0. The path is recovered from the last set of predecessors by starting at the end and going backwards, as above, until the predecessor 0 of node 1 is reached, which is the signal to stop. The predecessors matter in addition to the path because from them the shortest path from node 1 to any of the other nodes can be determined: five different paths, all from a single vector, an economy over listing five separate paths.
Only node 1 is labeled at the start. The destination is just marked, and the procedure keeps working until it labels the destination; as soon as it does, it can stop, since a label is never revised, even with many nodes still unlabeled. On the whole U.S. road network, a route from Raleigh to Charlotte can stop as soon as Charlotte is labeled, with no need to consider all the nodes out west in California, and at that point the labels give the shortest distance from the starting point to each intermediate node labeled. Had there been a direct arc of length 5 from node 1 to node 6, for example, node 6 would have been labeled right after nodes 1, 3 and 2, since its 5 would be less than the distances then held by nodes 4 and 5, and the procedure would have stopped with those two left unlabeled. An arc of length zero causes no trouble either, as where two facilities are at the same location: the procedure would definitely use it.
Example 1(b): In Julia
Determine the same path with the dijkstra_shortest_paths function of Graphs.jl and with the implementation of Model 1, and determine whether both agree with Ex. 1(a).
The network can be created by defining the starting and ending nodes for each arc along with each arc’s weight:
# Code block 16: the six-node network, from its arc listi = [1, 1, 2, 2, 3, 3, 4, 4, 5] # Index of starting node of each arcj = [2, 3, 3, 4, 4, 5, 5, 6, 6] # Index of ending node of each arc# Weight (e.g., cost, distance, time) of each arcw = [4, 2, 1, 5, 8, 10, 2, 6, 3]g =SimpleWeightedGraph(i, j, w) # Create undirected network
{6, 9} undirected simple Int64 graph with Int64 weights
The shortest paths are found first using Dijkstra’s procedure from Graphs.jl (very efficient compared to the implementation of Model 1, called after it):
ds =dijkstra_shortest_paths(g, 1); # Code block 17: from node 1 to all
ds.dists # Code block 18: distance from node 1 to all nodes
6-element Vector{Int64}:
0
3
2
8
10
13
ds.parents # Code block 19: node predecessors (also called parents)
6-element Vector{Int64}:
0
3
1
2
4
5
p =enumerate_paths(ds, 6) # Code block 20: path from node 1 to node 6
6-element Vector{Int64}:
1
3
2
4
5
6
ds.dists[6] # Code block 21: distance from node 1 to 6
13
The implementation of Model 1 operates on the adjacency matrix representation of the network. Its only inputs are the adjacency matrix and the starting node.
W =Array(float(adjacency_matrix(g))) # Code block 22: adjacency matrix
# Code block 23: non-arcs to InfW = [has_edge(g, i, j) ? W[i, j] :Inf for i in1:nv(g), j in1:nv(g)]
1
A 0 in the adjacency matrix means no arc, which the implementation would read as an arc that costs nothing, so each pair with no arc gets Inf and is never selected. has_edge keeps a real arc of length 0 as it is, such as the arc 2 => 4 of Code blocks 4 to 6.
ds.dists, ds.parents # Code block 25: using dijkstra_shortest_paths(g, 1)
([0, 3, 2, 8, 10, 13], [0, 3, 1, 2, 4, 5])
The hand trace of Ex. 1(a), the implementation of Model 1 and Graphs.jl give the same distances and predecessors, and so the same path of length 13, which is the Triangulate check of Lecture 1.2.
Example 1(c): Directed arcs
Determine the shortest path from node 1 to node 6 when each arc can be traveled only from its starting node to its ending node.
The only change is that each arc in the network is directed from the starting node (i) to the ending node (j). The function SimpleWeightedDiGraph is used to create a directed network, as opposed to using SimpleWeightedGraph to create an undirected network:
{6, 9} directed simple Int64 graph with Int64 weights
# Code block 27: the shortest path from node 1 to node 6, directedds =dijkstra_shortest_paths(g, 1)p =enumerate_paths(ds, 6)println("Min cost path: ", p, ", Cost: ", ds.dists[6])
Min cost path: [1, 2, 4, 5, 6], Cost: 14
The undirected path of Ex. 1(a) and 1(b) was [1, 3, 2, 4, 5, 6] at a cost of 13. In the directed graph, there is no edge from 3 to 2, resulting in the different path and higher cost, both shown in Fig. 7.
Show the code that draws this figure
# Code block 28: the directed network and its shortest pathfig =Figure(size = (420, 240))ax =Axis(fig[1, 1]; aspect =DataAspect())hidedecorations!(ax); hidespines!(ax)net6!(ax; directed =true, path = p)fig
Figure 7: The six-node network with each arc directed from its lower-numbered node to its higher-numbered one. The shortest path from node 1 to node 6, in red, no longer uses the arc between nodes 2 and 3, which runs only from 2 to 3.
3. Other shortest-path procedures
The shortest path problem can also be viewed as a minimum cost network flow problem and solved as a linear program, the form Lecture 4.2 takes up. That is worth being aware of as an additional modeling capability: a shortest path problem with additional considerations, such as bounds on its arcs or nodes, can include them when it is modeled as a minimum cost network flow, which Dijkstra’s procedure cannot do. A standard shortest path problem, though, is solved directly by Dijkstra’s procedure, as long as no arc length is negative, assumption (a) of Model 1. With negative arc lengths, other shortest path procedures are needed. They typically do not arise in logistics engineering applications outside of project scheduling, which this course does not cover. They are described briefly here, to show where Dijkstra’s procedure sits among them.
Dijkstra’s procedure is the one needed for most road settings, and it works well. The one situation where it does not work is where the arc lengths are a mix of positive and negative. Dijkstra’s procedure is a label setting procedure; the other techniques are called label correcting, and the most famous of them is the Floyd-Warshall algorithm, which is O(n^3) where Dijkstra’s procedure is O(n^2), and which works when some of the arcs are negative. It is fair to ask why such a procedure deserves attention when Dijkstra’s has always sufficed, and in what case an arc length would ever be negative. One rule of thumb is that if something exists and is talked about, there is a good chance it solves some problem that Dijkstra’s procedure does not, and here the problem is project scheduling. A task that cannot start until a minimum time has passed, and has to be finished by some maximum time, is modeled as a network with a negative arc, and that network needs a label correcting procedure: for any significant project scheduling, Floyd-Warshall is the way to go.
Both a label setting procedure like Dijkstra’s and a label correcting procedure like Floyd-Warshall keep a distance label and a predecessor at every node. What is different is that a node’s label becomes permanent when the node is labeled in a label setting procedure, while in a label correcting procedure every label is temporary and provides only an upper bound. Permanently labeling the node is safe only because, when no arc is negative, a node labeled later can never reach it earlier or more cheaply. The label correcting procedure only stops when no arc can improve any label. Since a label can change many times, the label correcting procedure can be slower than Dijkstra’s procedure. Also, if there is a cycle of negative total length in a graph, a label correcting procedure will never terminate, and so any procedure must have a method of detecting the negative cycle.
Negative arcs are used in project scheduling to represent maximum lags between activities. A company scheduling the maintenance of oil refineries is one example. Assume that each refinery has to be shut down every three years or so for maintenance for about three months, with everything torn down and hundreds of workers on site, and saving one day of the time the maintenance requires is worth millions of dollars. Project scheduling represents the lags allowable between activities as arcs. An activity has to occur no sooner than one day after some event, say, but no later than three days after it. That is a minimum and a maximum lag, and the minimum is a positive arc length while the maximum is a negative one.
Each activity’s earliest start time is the length of the longest path leading to it, so every arc can only express a “no earlier than” constraint on the activity it points to. A maximum lag is a “no later than” constraint, so it has to be restated from the other side: if the activity can start at most three days after the event, then the event can occur at the earliest three days before the activity. That restated constraint becomes an arc pointing back from the activity to the event, with a length of minus three.
Lags like these are relative time windows: once a task is complete, the next has to be done within one hour, say, perhaps because of some chemical process, and when the first will be done is not known in advance. A relative time window has to be solved with a shortest path procedure with negative arcs, where Dijkstra’s procedure does not work and Floyd-Warshall is used. The time windows of vehicle routing, the subject of Lecture 5.3, are absolute instead, such as a delivery window that opens at 1 p.m., and an absolute time window can be handled by a simple scanning procedure, with no shortest path procedure needed. Table 3 sets out the time complexity of these procedures for n nodes and m arcs, beside that of a linear program and of the Hungarian procedure for the transportation problem of Lecture 4.2.
Table 3: Time complexity of shortest path procedures.
Procedure
Problem
Time Complexity
Simplex
LP
O(2^n)
Ellipsoid
LP
O(n^4)
Hungarian
Transportation
O(n^3)
Floyd-Warshall
Shortest Path with Negative Arcs (all-pairs)
O(n^3)
Bellman-Ford
Shortest Path with Negative Arcs (single-src)
O(m \cdot n)
Dijkstra (linear min)
Shortest Path with Nonnegative Arcs
O(n^2)
Dijkstra (Fibonacci heap)
Shortest Path with Nonnegative Arcs
O(m \log n)
Number of nodes
n
Number of arcs
m
The following compares the time complexity of determining the road distance from Raleigh, NC, to Portland, OR, using the U.S. Interstate Highway network:
# Code block 29: the complexities, for the Interstate networkn =862# Number of nodes in Interstate road networkm =1223# Number of arcs@show n^4# LP for shortest path using ellipsoid@show n^2# Dijkstra using linear search for minimum@showm*log(n) # Dijkstra using Fibonacci heap for minimum@show (n^4)/(m*log(n)) # Ratio LP to Fibonacci@show (n^2)/(m*log(n)); # Ratio linear search to Fibonacci
n ^ 4 = 552114385936
n ^ 2 = 743044
m * log(n) = 8266.569196021695
n ^ 4 / (m * log(n)) = 6.6788818050625674e7
n ^ 2 / (m * log(n)) = 89.88541466000086
The arithmetic makes the table concrete. For the six nodes of Ex. 1, Dijkstra’s procedure as the implementation of Model 1 runs it, looking through all of the unlabeled nodes for the minimum, does work proportional to n^2 = 36. Selecting the minimum unlabeled node is the one step a Fibonacci heap speeds up, reducing it from n to \log n, and n \log n is only 10.75. A linear program, at best n^4, would be 1,296. Either Dijkstra procedure is a lot more efficient than a linear program, which for shortest paths is worth using only to add other features, since such problems explode at realistic sizes. A more realistic network has about 70,000 nodes, roughly the road nodes in the United States that trucks can travel on, leaving out the residential streets. With a Fibonacci heap the figure is under a million, which is fast; the naive n^2 procedure is just under 5 billion, which is quite time-consuming; and a linear program is far larger than anyone would consider solving. Any realistic routing procedure needs the Fibonacci heap’s capability.
On the Interstate network of Code block 29, the linear search does about 90 times the work of the heap.
NoteA* algorithm
The A* algorithm adds to Dijkstra’s procedure a heuristic lower bound (LB) estimate of each node’s remaining distance to the destination.3 It is used in AI search for all types of applications (tic-tac-toe, chess).4 In path planning applications, great circle distance from each node to destination could be used as an LB estimate of the remaining distance.
Finding the shortest path from Raleigh to Dallas with Dijkstra’s procedure starts at Raleigh, where roads leave in all directions, 360 degrees around it, and keeps labeling road and intersection nodes until it labels Dallas and can stop. The labeled nodes, green in Fig. 8, cover the eastern half of the United States: to get to Dallas the procedure labels half the country, which is not terribly efficient. The A* algorithm adds to each node’s distance from Dijkstra’s procedure an estimate of a lower bound on its remaining distance to the destination, here the great circle distance from the node directly to Dallas, which can be computed very quickly. For A* the estimate has to be a lower bound: it can never overestimate, or give a middle estimate, and always has to underestimate the actual value. The great circle distance is perfect for this, since it is the distance as the crow flies and a road distance is always a little greater. Every time the procedure updates a node’s cost it adds the node’s great circle distance, as Eq. 4 below shows, and that works like a cutting plane. Of two nodes, if the Dijkstra distance plus the lower bound of one already exceeds the other, the first can be eliminated without losing anything, since the bound never overestimates. So everything outside the ellipse in Fig. 8 never has to be looked at. The same idea applies to games: a heuristic that puts a value on a position, such as two bishops and a knight, can be used as long as it never overestimates.
For each node i it reaches, A* adds the two distances:
d_{A^{\star}}(\text{Raleigh}, \text{Dallas}) = d_{\text{dijk}}(\text{Raleigh}, i) + d_{GC}(i, \text{Dallas}), \quad \text{for each node } i
\tag{4}
where
d_{A^{\star}}(\text{Raleigh}, \text{Dallas})
=A* estimate of the distance from Raleigh to Dallas by way of node i (mi)
d_{\text{dijk}}(\text{Raleigh}, i)
=distance from Raleigh to node i found by Dijkstra’s procedure (mi)
d_{GC}(i, \text{Dallas})
=great circle distance from node i to Dallas, a lower bound on the road distance (mi).
Show the code that draws this figure
# Code block 30: Dijkstra against A*, Raleigh to DallasusingGeoMakieN, L =prune_reindex(faf5nodes(), faf5links())nN =nrow(N)# Both directions of each link, one-way links honored, and of two parallel# links the shorter kept (Sec. 4 explains both)fw, bw = L.DIR .!=-1, L.DIR .!=1gr =SimpleWeightedDiGraph([L.SRC[fw]; L.DST[bw]], [L.DST[fw]; L.SRC[bw]], [L.DIST[fw]; L.DIST[bw]]; combine = min)ral =loc2lonlat("Raleigh"; state =:NC)dal =loc2lonlat("Dallas"; state =:TX)LL = [N.LON N.LAT]closest(c) =argmin(dgc.(eachrow(LL), Ref([c.LON, c.LAT])))s, t =closest(ral), closest(dal)dsr =dijkstra_shortest_paths(gr, s)dRD = dsr.dists[t] # Raleigh to Dallas, mih =dgc.(eachrow(LL), Ref(LL[t, :])) # great circle to Dallas, midijk =findall(<(dRD), dsr.dists) # labeled by Dijkstraastar =findall(v -> dsr.dists[v] + h[v] < dRD, 1:nN) # by A*fig, ax =makemap(N.LON[dijk], N.LAT[dijk])resize!(fig, 687, 430)scatter!(ax, N.LON[dijk], N.LAT[dijk]; color = (dcc, 0.6), markersize =1.5)scatter!(ax, N.LON[astar], N.LAT[astar]; color = (supc, 0.8), markersize =2)scatter!(ax, [ral.LON, dal.LON], [ral.LAT, dal.LAT]; color = ink, strokecolor =:white, strokewidth =1, markersize =10)fig
Figure 8: Raleigh to Dallas on the national road network. Green: the nodes Dijkstra’s procedure labels before it reaches Dallas. Red: the far smaller set A* labels, with the great circle distance to Dallas as its lower bound. Raleigh and Dallas are the dark dots, in the east and the west.
On the road network of Sec. 4, the shortest route from Raleigh to Dallas is 1,158.2 mi, and the great circle distance is 1,053.9 mi, below it as a lower bound must be: a Bounds check of Lecture 1.2. Dijkstra’s procedure labels 215,602 nodes before it reaches Dallas, and A* only 36,975, about 17% as many.
Commercial route planners are faster still. How Google Maps finds its routes has not been published, but the documented planners, among them Bing Maps and the open-source OSRM used in Ex. 3, process the road network once in advance, so that each route then takes milliseconds or less, and on a detailed road network a great circle lower bound speeds up A* little, if at all.5
4. Road networks
Road data typically cannot just be downloaded and used as is: some processing is needed to make it into a more useful form. The first modification is thinning, which, like the others, is something like an algorithm, a series of conditions applied to the data.
Each of the modifications, thinning, pruning and reindexing, and adding connectors, is a function in Logjam. The thinning below works on the Interstate links of a national road network, whose two CSV files are linked under Companion script at the top of the page. That network is derived from FAF5, the national road network that Logjam’s faf5nodes and faf5links return and Ex. 2 below uses. FAF5 keeps each side of a divided highway as its own one-way link and each interchange as a set of ramps, so its Interstates in the lower 48 states are reduced to a single two-way road between each pair of neighboring interchanges, with the nodes along the road kept.6
# Code block 31: the Interstate network, from its two CSV filesusingCSVdfN =DataFrame(CSV.File("data/interstate_nodes.csv"))dfL =DataFrame(CSV.File("data/interstate_links.csv"))nrow(dfN), nrow(dfL), maximum(dfN.IDX)
(34746, 35165, 1949186)
Each node keeps its index in the whole road network, which is why the largest index, 1,949,186, is far above the number of nodes, 34,746.
Thin
To thin a network, remove all degree-2 nodes from the network, and add the cost of both arcs incident to each degree-2 node to make one arc of the two. If that results in multiple arcs between a pair of nodes, keep the one of minimum cost.
# Code block 32: a schematic of I-40 around Raleigh, thinnedXY0 = [1.33-6.70; 1.27-7.89; 1.67-7.28; 2.97-7.78; 4.05-7.27;4.45-6.68; 4.83-7.91; 5.53-6.62; 5.62-7.66; 6.41-7.47;6.59-6.79; 7.07-7.86] # where each node is drawnA0 = [135; 235; 3420; 4515; 5610; 6815; 8915; 5715;7910; 91010; 10115; 10125] # arc list: i, j, costN0 =DataFrame(IDX =1:size(XY0, 1), X = XY0[:, 1], Y = XY0[:, 2])L0 =DataFrame(SRC = A0[:, 1], DST = A0[:, 2], DIST = A0[:, 3])tN0, tL0 =thin(N0, L0)prt(tL0)
# Code block 33: the schematic, before and after thinningbluec = colorant"#1f6fb5"# the nodesrd, rk =0.07, 0.15# node dot and ring radiiside0 = [0, 0, -1, -1, 1, 1, 1, -1, -1, -1, 0, 0] # side of each label# The two ends of the segment from the outline of the circle of radius# rp at p to the outline of the circle of radius rq at qfunctionrim(p, q, rp, rq) u = (q .- p) ./hypot((q .- p)...)return p .+ rp .* u, q .- rq .* uendfig =Figure(size = (620, 230))ax =Axis(fig[1, 1]; aspect =DataAspect())hidedecorations!(ax); hidespines!(ax)for (k, (i, j, w)) inenumerate(eachrow(A0)) # the network p, q = XY0[i, :], XY0[j, :] t, h =rim(p, q, rd, rd)lines!(ax, [t[1], h[1]], [t[2], h[2]]; color = ink, linewidth =1.2) side0[k] ==0&&continue# the spurs go unlabeled u = (q .- p) ./hypot((q .- p)...) m = (p .+ q) ./2.+0.17* side0[k] .* [-u[2], u[1]]text!(ax, m[1], m[2]; text =string(w), color = ink, fontsize =13, align = (:center, :center))endfor (i, j) inzip(tL0.SRC, tL0.DST) # what thinning leaves t, h =rim(XY0[i, :], XY0[j, :], rk, rk)lines!(ax, [t[1], h[1]], [t[2], h[2]]; color = (supc, 0.85), linewidth =2.5)endkept =Set([tL0.SRC; tL0.DST])for k in1:nrow(N0) c =Point2f(XY0[k, 1], XY0[k, 2])poly!(ax, Circle(c, rd); color = bluec) k in kept &&poly!(ax, Circle(c, rk); color =:transparent, strokecolor = supc, strokewidth =2.5)endtext!(ax, 3.9, -7.5; text =string(maximum(tL0.DIST)), color = supc, fontsize =15, font =:bold, align = (:center, :center))fig
Figure 9: A network before and after thinning. Every degree-2 node between the two junctions goes, and the five arcs that remain, in red, join the six nodes circled in red, one arc of cost 70 replacing all of the arcs between the junctions.
The schematic of Fig. 9 is the actual road segment around Raleigh: I-40 runs between its junctions with I-85 on the left and I-95 on the right, and the Beltline around Raleigh forms the loop between them. Thinning first eliminates all degree-2 nodes from the network, the nodes with just one arc in and one arc out. The two junction nodes are of degree three, and all of the others between them are of degree two. Removing a degree-2 node replaces its two arcs with a single arc whose cost is their sum, so arcs of 10 and 15 become a single arc of 25. Where that leaves two arcs between the same pair of nodes, such as the two sides of the Beltline, only the minimum cost one is kept: a traveler would always take the shorter, so the longer can be forgotten. Applying the procedure repeatedly, until nothing else needs to be thinned, leaves the junction of I-40 and I-85 joined to the junction of I-40 and I-95 by a single arc with a total distance of 70. What is interesting is that the entire Beltline around Raleigh is eliminated. Whether that is good from a practical point of view is a fair question, since Raleigh is a whole city: to use it, there has to be a way of identifying where Raleigh is, of embedding the location of a city in the network, which is the third procedure, adding connectors. Thinning ignores all of that. Fig. 10 is the more detailed version, with all of the actual nodes around Raleigh. It includes I-540, which runs across the north of Raleigh between two different junctions, so thinning keeps it as one arc between them.
Show the code that draws this map
# Code block 34: the Interstates around Raleigh, before and after thinningRRal = [-79.261835.2497; -78.196436.1762] # the rectangle: lon, lat# The nodes inside the rectangle, and the links between themrN, rL =cropnetwork(dfN, dfL, RRal[:, 1], RRal[:, 2]; xexpand =0, yexpand =0)trN, trL =thin(rN, rL)fig, ax =makemap(rN.LON, rN.LAT; doRoadbkgd =false, xexpand =0.05, yexpand =0.05)resize!(fig, 400, 420)gR =SimpleWeightedGraph(rL.SRC, rL.DST, rL.DIST)lines!(ax, x2ln(gR, rN.LON), x2ln(gR, rN.LAT); color = supc, linewidth =1)scatter!(ax, rN.LON, rN.LAT; color =:transparent, strokecolor = supc, strokewidth =1, markersize =5)for (i, j) inzip(trL.SRC, trL.DST)lines!(ax, rN.LON[[i, j]], rN.LAT[[i, j]]; color = (dcc, 0.8), linewidth =4)endscatter!(ax, rN.LON[trN.IDX], rN.LAT[trN.IDX]; color = ink, markersize =6)# The data carries no road names, so the routes are labeled by handlabs = [("I-85", -78.755, 36.07, (:left, :center)), ("I-40", -79.02, 35.905, (:right, :center)), ("I-540", -78.68, 35.94, (:center, :bottom)), ("I-440", -78.69, 35.735, (:center, :top)), ("I-40", -78.53, 35.58, (:left, :center)), ("I-95", -78.32, 35.46, (:center, :top)), ("Raleigh", -78.70, 35.665, (:center, :top)), ("Durham", -78.89, 36.052, (:center, :bottom))]for (t, x, y, a) in labstext!(ax, x, y; text = t, fontsize =12, color = ink, align = a)endfig
Figure 10: Thinned Interstates around Raleigh. The Interstate links inside a rectangle around the city, each node a red circle, and in green the links thinning leaves, between the nodes in black. I-440 is the Beltline, and I-540 runs across the north of the city.
Around Raleigh, thinning takes the 228 Interstate links inside the rectangle to 12. Thinning the whole Interstate network leaves each remaining node with its original index, so the nodes are then numbered from 1 by prune_reindex, the function of the next modification:
# Code block 35: thin the Interstate network, then reindex its nodestN, tL =thin(dfN, dfL)@shownrow(tN), nrow(tL), maximum(tN.IDX)tN, tL =prune_reindex(tN, tL)@shownrow(tN), nrow(tL), maximum(tN.IDX);
# Code block 36: the Interstate network, before and after thinningfunctionimap(N, L, c) g =SimpleWeightedGraph(L.SRC, L.DST, L.DIST) fig, ax =makemap(N.LON, N.LAT; doRoadbkgd =false, xexpand =0.02, yexpand =0.02)resize!(fig, 340, 240) ax.title =string(commas(nrow(N)), " nodes, ",commas(nrow(L)), " links") ax.titlesize, ax.titlefont =13, :regularlines!(ax, x2ln(g, N.LON), x2ln(g, N.LAT); color = c, linewidth =0.6)return figendimap(prune_reindex(dfN, dfL)..., supc) # indices from 1, for x2lnimap(tN, tL, ink) # Code block 37: the thinned network
(a) Not thinned
(b) Thinned
Figure 11: The Interstate network before and after thinning.
Thinning takes the Interstate network from 35,165 links to 1,024. Fig. 11 (a) is the network as it is, and Fig. 11 (b) the network thinned.
Taking just the links flagged as Interstate, with their nodes, gives the network on the left of Fig. 11. It is worth asking why it wiggles so much, rather than running in straight lines like the one on the right. The wiggles are the nodes along each Interstate where one FAF5 link meets the next, and they follow the course of the road. The problem is that they add about 34,000 arcs, a fairly big data set, and every one is extra work for a shortest path procedure such as Dijkstra’s working along the road. They are redundant when all that is wanted is to get from the beginning of a stretch to its end, which a straight line does, and thinning gives the network on the right.
Julia is fast enough that thinning is now optional. In MATLAB, finding shortest paths was slow enough that thinning really helped. Now the actual road network can be used as it is, with all of its arcs, and there are advantages in not having to thin it at all.
Prune and reindex
To prune and reindex a network is to extract the portion of its graph with only those nodes and/or arcs that satisfy some condition, numbering the nodes that remain from 1.
Pruning and reindexing are needed even when thinning is not. Given a road network for the entire United States, a problem about North Carolina or another smaller region wants just a portion of it. In Fig. 12, suppose the criterion is arcs with weights less than 35. Eliminating all of the arcs of 35 or more leaves the set of nodes shown under “Arcs < 35”: nodes 5 and 7 lose all of their arcs, and it is easy to throw them away, while the others stay. The one problem is that node 6 of the original network is still numbered 6, after nodes 1 to 4. That needs to be corrected, because a node labeled 6 makes the representation assume six nodes: the adjacency matrix would be 6 × 6, and the graph would be incomplete, with no node 5. So node 6 is reindexed as node 5, starred in the figure, giving a five-node graph. The same is done for a region of the road network, the nodes inside the dashed rectangle, which have to be relabeled 1, 2, 3 and 4. Pruning takes a subset of the arcs and nodes, and reindexing renumbers the nodes so that their indices run from 1 up to however many nodes the graph has, which is necessary for all of the graph procedures to work correctly.
# Code block 38: seven cities near Greensboro, pruned two waysgso =loc2lonlat("Greensboro"; state =:NC)C7 =filter(r -> r.POP >30_000&&dgc([gso.LON, gso.LAT], [r.LON, r.LAT]) <50, usplace())sort!(C7, [:ST, :NAME]) # nodes 1 to 7# Arc list: neighboring cities, each cost about 1.2 x great circle (mi)A7 = [1231; 1325; 1436; 1742; 2464; 25100; 2760;3414; 3629; 3749; 4543; 4621; 5639; 6769]N7 =DataFrame(IDX =1:nrow(C7), LON = C7.LON, LAT = C7.LAT)L7 =DataFrame(SRC = A7[:, 1], DST = A7[:, 2], DIST = A7[:, 3])println(join(string.(1:4, " ", C7.NAME[1:4]), ", "))println(join(string.(5:7, " ", C7.NAME[5:7]), ", "))aN, aL =prune_reindex(N7, filter(r -> r.DIST <35, L7)) # arcs < 35R7 = [-80.608535.5897; -79.697136.2184] # a rectangle: lon, latcN, cL =cropnetwork(N7, L7, R7[:, 1], R7[:, 2]; xexpand =0, yexpand =0) # nodes inside itprt(aL)prt(cL)
1
Sorting by state and name fixes the order of the cities, so the node numbers in the arc list A7 refer to the right cities.
2
cropnetwork widens the rectangle by default, by an expansion factor of 0.1 on each axis; a factor of 0 keeps exactly the rectangle.
# Code block 39: the seven-city network and its two subgraphs, drawn# Lon and lat drawn to scale at 36 degrees northxy(lon, lat) = [lon *cosd(36), lat]# The network on ax: nodes as numbered circles, each arc from outline to# outline and labeled with its cost, on side s (1 or -1) and fraction f# of the way along it, as given for an arc (i, j) in labfunctioncitynet!(ax, N, L; r =0.034, fs =12, labels =true, lab =Dict())for (i, j, w) inzip(L.SRC, L.DST, L.DIST) p, q =xy(N.LON[i], N.LAT[i]), xy(N.LON[j], N.LAT[j]) u = (q .- p) ./hypot((q .- p)...) t, h = p .+ r .* u, q .- r .* ulines!(ax, [t[1], h[1]], [t[2], h[2]]; color = ink, linewidth =1) labels ||continue s, f =get(lab, (i, j), (1, 0.5)) m = p .+ f .* (q .- p) .+0.045s .* [-u[2], u[1]]text!(ax, m[1], m[2]; text =string(w), color = ink, fontsize = fs, align = (:center, :center))endfor k in1:nrow(N) x, y =xy(N.LON[k], N.LAT[k])poly!(ax, Circle(Point2f(x, y), r); color =:white, strokecolor = bluec, strokewidth =1.5)text!(ax, x, y; text =string(k), color = bluec, fontsize = fs -1, align = (:center, :center))endend# A subgraph moved dx and dy (in drawn units) to its place in the figureshift(N, dx, dy) =DataFrame(IDX = N.IDX, LON = N.LON .+ dx /cosd(36), LAT = N.LAT .+ dy)fig =Figure(size = (687, 580))ax =Axis(fig[1, 1]; aspect =DataAspect())hidedecorations!(ax); hidespines!(ax)citynet!(ax, N7, L7; lab =Dict((1, 3) => (-1, 0.5), (3, 4) => (-1, 0.62), (3, 6) => (-1, 0.5), (2, 5) => (-1, 0.5), (4, 5) => (-1, 0.5), (1, 4) => (1, 0.3)))(x0, y0), (x1, y1) =xy(R7[1, 1], R7[1, 2]), xy(R7[2, 1], R7[2, 2])lines!(ax, [x0, x1, x1, x0, x0], [y0, y0, y1, y1, y0]; color = supc, linestyle =:dash, linewidth =1.8)text!(ax, xy(-80.50, 36.62)...; text ="(a) Network", fontsize =14, align = (:left, :center))# Each subgraph is derived from (a), so an arrow runs from (a) to eachdyb, dxc =-0.85, 1.55# where the two subgraphs are putarrow!(p, q) =arrows2d!(ax, [p[1]], [p[2]], [q[1] - p[1]], [q[2] - p[2]]; color = (ink, 0.7), shaftwidth =1.5, tipwidth =11, tiplength =11)arrow!(xy(-79.55, 35.76), xy(-79.55, 35.47)) # (a) to (b), downarrow!(xy(-78.96, 36.15), xy(-78.49, 36.15)) # (a) to (c), righta2 =shift(aN, 0.0, dyb) # (b) arcs < 35, belowx5, y5 =xy(a2.LON[5], a2.LAT[5])scatter!(ax, [x5], [y5]; marker =:star5, markersize =36, color = (supc, 0.85)) # the renumbered nodecitynet!(ax, a2, aL; lab =Dict((1, 3) => (-1, 0.5), (3, 4) => (-1, 0.62), (3, 5) => (-1, 0.5)))text!(ax, xy(-79.55, 35.39)...; text ="(b) Arcs < 35", fontsize =14, align = (:center, :center))c3 =shift(cN, dxc, 0.2) # (c) nodes in the rectanglecitynet!(ax, c3, cL; lab =Dict((1, 2) => (1, 0.5)))text!(ax, xy(-80.16+ dxc /cosd(36), 36.45)...; text ="(c) Nodes in rectangle", fontsize =14, align = (:center, :center))fig
Figure 12: Pruning and reindexing. (a) The seven cities joined by arcs costed at about 1.2 times their great circle distance, with a dashed rectangle. (b) The subgraph of the arcs of (a) under 35, its renumbered node starred. (c) The subgraph of the nodes of (a) inside the rectangle. Each subgraph’s nodes are numbered from 1.
Pruning by an attribute such as a state can go wrong at the edge of a region. Suppose the ten most populous cities in Texas need a ten by ten matrix of shortest path distances, as Ex. 2 below builds for the ten North Carolina cities of over 100,000 people. Texas has a big chunk missing in its northwest corner, beside the Panhandle, and the most obvious route between some of its cities goes outside the state, probably through Oklahoma. Taking only the roads in Texas eliminates that route, and the path has to go a longer way around to stay in Texas, where cutting through Oklahoma would have been fine. The way to account for that is to restrict the nodes instead of the links: take the longitude and latitude of all ten cities and put a bounding box around them, not just fitting them but with some extra expansion, and use it to restrict the nodes of the road network. That keeps not just the Texas roads but every road that fits in the rectangle, a larger road network, and the result goes through Oklahoma, which is more realistic. Logjam’s cropnetwork does exactly that, as it did for the rectangle in Fig. 12.
Add connectors
To add connectors, given new nodes, add arcs that connect the new nodes to the existing nodes in a graph and to each other. The distance of a connector arc is the great circle distance times a circuity factor for connectors, 1.3 by default. Each new node is connected to the three existing nodes at the corners of the triangle that contains it, in a triangulation of the existing nodes.
The last modification starts from points rather than roads. A problem is quite often about points such as the street address of a factory or a distribution center, the address of a customer, or the center of population of a city or a ZIP code: demand points with a longitude and latitude, and there is no reason to expect any of them to fall on a node of the road network. These new nodes are typically what the problem is actually about, and they have to be connected, by creating additional links from each new node to the existing road network. These are roads that do not actually exist, virtual roads, added so that the question of the shortest distance from a customer to a plant can be asked meaningfully. A path then runs from one new node along a connector to the road network, along the road network, and along another connector to its destination; between two new nodes close enough together, it need not get on the road network at all. Since the longitude and latitude of each new node are known, the great circle distance from it to the road nodes can be calculated, but that is a straight line, not reality. The circuity factor of 1.2 from Lecture 2.2, Sec. 7, could serve as a default, but connectors lie within a small region, where a higher value, such as 1.3 or even 1.4, tends to work better. The road network is less dense in less populated areas, so a point may not be near a major road, and travel to it may not be very direct, perhaps along a lot of back roads; the higher factor accounts for that. It is just a matter of what works reasonably well in practice.
# Code block 40: six new nodes, connected to the seven-city networkx1 = [-79.990, -80.308, -79.629, -79.899, -79.554, -79.143] # lony1 = [35.921, 35.736, 35.894, 36.206, 36.232, 36.321] # latn1 =length(x1)N8, L8 =addconnectors(N7, L7, x1, y1)isc =isequal.(L8.SOURCE, "CONNECTOR") # the connectorsnnx =count(isc .& (L8.DST .> n1)) # new node to existingnnn =count(isc .& (L8.DST .<= n1)) # new node to new nodeprt(sort(L8[isc .& (L8.SRC .==1), 1:3], :DST)) # those of new node 1
1
A road link has no SOURCE value, and == with a missing value gives missing rather than false; isequal gives false, so the result can be used to select rows.
# Code block 41: the connectors, drawn with Code block 39's citynet!fig =Figure(size = (687, 560))ax =Axis(fig[1, 1]; aspect =DataAspect())hidedecorations!(ax); hidespines!(ax)P8 = [N8.LON .*cosd(36) N8.LAT] # every node, drawn as xy() doesrn =0.02# a new node's radiushs =Any[]for (i, j) inzip(L8.SRC[isc], L8.DST[isc]) # the connectors, first p, q = P8[i, :], P8[j, :] newnew = j <= n1 u = (q .- p) ./hypot((q .- p)...) t = p .+ rn .* u # outline to outline h = q .- (newnew ? rn :0.034) .* upush!(hs, lines!(ax, [t[1], h[1]], [t[2], h[2]]; color = newnew ? bluec : supc, linestyle = newnew ? :dashdot ::dash, linewidth =1.8))endcitynet!(ax, N7, L7; labels =false, fs =14) # the cities, over themfor k in1:n1 # the new nodes, filledpoly!(ax, Circle(Point2f(P8[k, 1], P8[k, 2]), rn); color = bluec, strokecolor =:white, strokewidth =1)endk1 =findfirst(>(n1), L8.DST[isc]) # a connector of each kindk2 =findfirst(<=(n1), L8.DST[isc])sc =MarkerElement(color = bluec, marker =:circle, markersize =14)axislegend(ax, [hs[k1], hs[k2], sc], ["New to existing node", "New to new node", "New node"]; position =:lt, labelsize =14)fig
Figure 13: Connectors added to the seven-city network of Fig. 12. Each of the six new nodes is joined to the cities at the corners of the triangle that contains it, two for the node outside every triangle (red, dashed), and to its neighboring new nodes (blue, dash-dot).
Fig. 13 shows the result. addconnectors numbers the new nodes first, 1 to 6, and shifts the seven cities to 7 to 13. Of its connectors, 17 join a new node to a city and 11 join two new nodes.
NoteDelaunay triangulation
A Delaunay triangulation finds the triangulation of points that maximizes the minimum angle of any triangle, and it captures proximity relationships.
The triangulation that adding connectors uses is the Delaunay triangulation. A triangulation of a set of points divides the region they span into shapes that are all triangles, and any set of points can be triangulated. What makes a triangulation a Delaunay triangulation is that it is the single best one, the triangulation that maximizes the minimum angle of any of its triangles. It would be easy to triangulate with long, skinny triangles, two long sides and one very short one, but that is not a good triangulation. It is remarkable that the Delaunay triangulation can be calculated very efficiently. It also makes a good default road network. If the points are the population centers of a set of cities, a Delaunay triangulation is pretty close to the roads that connect them, unless a major waterway or something similar is in the way: two big cities near each other probably have a road connecting them, and that is true of every arc. So, given a set of points without an actual road network, calculate the Delaunay triangulation; it is usually a pretty good approximation. Traveling west in North Carolina from Raleigh to Asheville, for example, the direct way by distance is due west, but I-40 goes well north, up to Greensboro and Winston-Salem, because the road network follows the big cities, and the actual road distance pretty much follows that sequence of cities rather than the great circle. A Delaunay triangulation, with its arcs connecting Greensboro and Winston-Salem, captures that, which makes it better than the great circle distance. So there are three levels: the great circle distance; a Delaunay triangulation, a little better; and the actual road network, better yet.
Fig. 14 shows the Delaunay triangulation of the seven cities of Fig. 12, whose arcs are the ones Code block 38 lists, beside another triangulation of the same cities, in which the arcs from High Point (4) to Burlington (1) and to Greensboro (3) are replaced by arcs from Chapel Hill (2) to Greensboro (3) and to Winston-Salem (6). The two replacement arcs make long, skinny triangles.
Show the code that draws this figure
# Code block 42: the Delaunay triangulation and another triangulation# Each triangulation as its triangles, three cities to a rowT7 = [124; 127; 134; 137; 245; 346; 367; 456] # DelaunayTf = [123; 127; 137; 236; 246; 245; 367; 456] # 2 flippedP7 = [xy(N7.LON[k], N7.LAT[k]) for k in1:nrow(N7)] # drawn positions# The smallest angle of any triangle of T, in degreesfunctionminangle(T) a =Inffor t ineachrow(T), m in1:3 p, q, r = P7[t[m]], P7[t[mod1(m +1, 3)]], P7[t[mod1(m +2, 3)]] u, v = q .- p, r .- p a =min(a, acosd(sum(u .* v) / (hypot(u...) *hypot(v...))))endreturn aend# The arcs of a triangulation, from its trianglesarcs(T) =unique([minmax(t[a], t[b]) for t ineachrow(T) for (a, b) in ((1, 2), (2, 3), (1, 3))])fig =Figure(size = (687, 330))for (k, (T, nm)) inenumerate([(T7, "Delaunay"), (Tf, "Another")]) red = k ==2 ? [(2, 3), (2, 6)] :Tuple{Int, Int}[] # replacements A =setdiff(arcs(T), red) a =round(minangle(T), digits =1) ax =Axis(fig[1, k]; aspect =DataAspect(), titlesize =14, titlefont =:regular, title =string(nm, ": smallest angle ", a, "°"))hidedecorations!(ax); hidespines!(ax)citynet!(ax, N7, DataFrame(SRC =first.(A), DST =last.(A), DIST =0); labels =false, fs =13)for (i, j) in red # outline to outline, in red u = (P7[j] .- P7[i]) ./hypot((P7[j] .- P7[i])...) t, h = P7[i] .+0.034.* u, P7[j] .-0.034.* ulines!(ax, [t[1], h[1]], [t[2], h[2]]; color = supc, linewidth =2.5)endendfig
Figure 14: Two triangulations of the seven cities of Fig. 12. Left, the Delaunay triangulation; right, the same cities with the arcs from High Point (4) to Burlington (1) and to Greensboro (3) replaced by the two arcs in red, from Chapel Hill (2) to Greensboro (3) and to Winston-Salem (6). Each panel gives the smallest angle of any of its triangles.
Inside the procedure that adds connectors, the Delaunay triangulation is how the decisions are made. A Delaunay triangulation can be queried: for a particular point, it tells which triangle the point is inside, which is exactly what is needed to connect the point to the network. The Delaunay triangulation is critical to adding connectors.
National road network
Finding a national road network to compute on is something of a long story. The Interstate data of earlier offerings of this course came from the Oak Ridge National Highway Network, last updated in January 2008.
Oak Ridge National Laboratory built that network primarily for vehicle routing and scheduling, military movements included, with funding from the U.S. Army Forces Command, the National Guard Bureau, the Federal Highway Administration, the Military Traffic Management Command, the Bureau of Transportation Statistics and the Department of Energy.7 The national road network of this course is now FAF5, the highway network of version 5 of the Freight Analysis Framework (FAF), which the Bureau of Transportation Statistics and the Federal Highway Administration produce and Oak Ridge develops. It descends from the National Highway Planning Network that Oak Ridge built for the Federal Highway Administration in the early 1990s, and was built new for the FAF 5.0 release in 2022, taking in every road of the earlier FAF networks and of the National Highway System. It is routable as published, and Logjam returns it as faf5nodes and faf5links.8
The functions faf5nodes and faf5links return node and link data, respectively, for the U.S. FAF5 road network. Each node is defined by an index (IDX) and its longitude (LON) and latitude (LAT), and each link is defined by its starting (source) node index (SRC), ending (destination) node index (DST), its distance (DIST) in miles, and other link attributes like state and county Federal Information Processing Standards (FIPS) codes (STFIP and COFIP), its direction (DIR), and its speed and travel time each way. The network is directed, because some of its links are one-way.
# Code block 43: the FAF5 road networkN, L =prune_reindex(faf5nodes(), faf5links()) # nodes numbered from 1@shownrow(N), nrow(L)# A link's direction: 1 one-way SRC to DST, -1 one-way DST to SRC, 0 both@showcount(==(1), L.DIR)@showcount(==(-1), L.DIR)@showcount(==(0), L.DIR)nodepair =minmax.(L.SRC, L.DST) # the two nodes each link joinsnshare =length(nodepair) -length(unique(nodepair)) # a pair repeated
A divided highway is two one-way links, one for each direction of travel, and 1,427 links join a pair of nodes that another link already joins.
Getting the roads of Wake County works just like the earlier work with data frames. Besides each link’s source and destination, the i and j of its arc, and its distance in miles, the weight when distance is minimized, the national road data carries the state and county FIPS codes as attributes, an urban area code, and the speed of travel on the link, which is really useful: with the distance it gives a travel time, which can be used as the weight instead. There is also the name of the road and other information, so there are many attributes that can be used. Pruning and reindexing with the state FIPS code of North Carolina, 37, and the county FIPS code of Wake County, 183, leaves all of the roads of Wake County, and a graph can be created from them. Drawn, as in Fig. 15, they are the major roads of the county, the Interstates and major highways, without any of the residential streets. Two points then stand for two facilities, between which the shortest distance is wanted. Adding connectors joins each of the two blue points of Fig. 15 to three nodes of the road network, in red, and the shortest path from the first new node to the second can then be determined: it starts at node 1, the first new point, runs through road nodes, and ends at node 2, the other new point. Its length can be compared with the great circle distance and with Google Maps, which in 2024 gave about 12.9 miles.
Use STFIP == 37 (NC) and COFIP == 37183 (Wake County) to select FAF5 links, and prune_reindex to create node and link data frames that just contain Wake County roads.
# Code block 44: the roads of Wake County, pruned and reindexed# FAF5's county code includes the state's: 37 and 183 make 37183NW, LW =prune_reindex(N, filter(r -> r.STFIP ==37&& r.COFIP ==37183, L))nrow(NW), nrow(LW)
(833, 1172)
# Code block 45: two points connected, and the shortest path between themXY2 = [-78.7586452972295935.74329140697632;-78.7035264437191335.866083946889844]NW2, LW2 =addconnectors(NW, LW, XY2[:, 1], XY2[:, 2]; add_nf_nf =false)fw, bw = LW2.DIR .!=-1, LW2.DIR .!=1# arcs each waygW =SimpleWeightedDiGraph([LW2.SRC[fw]; LW2.DST[bw]], [LW2.DST[fw]; LW2.SRC[bw]], [LW2.DIST[fw]; LW2.DIST[bw]]; combine = min)D, P =shortestpaths(gW, 2) # between the two new nodespW =tracepath(P[1], 1, 2) # the path from node 1 to node 2D
1
A connector between two new nodes is only an estimate, 1.3 times their great circle distance, so when the road distance between them is what is wanted, the connectors are added with add_nf_nf = false, which joins the new nodes to the network alone.
2
A graph holds one arc for each pair of nodes and direction, so the graph of a road network is built by two rules: each link gives an arc in each direction it can be traveled (fw marks the links that can be traveled from SRC to DST, and bw those that can be traveled back), and of two arcs joining the same nodes in the same direction, the shorter is kept, the rule for parallel arcs in Sec. 1, which combine = min applies. Built from the arc list, a link of zero length stays an arc.
3
shortestpaths(g, k) runs Dijkstra’s procedure from each of the first k nodes and returns the distances between them and the predecessors; here the first 2 nodes are the two new points, since addconnectors numbers the new nodes first.
2×2 Matrix{Float64}:
0.0 12.6502
12.3379 0.0
Show the code that draws this map
# Code block 46: Wake County's roads, the two points and their pathgw =SimpleWeightedGraph(LW.SRC, LW.DST, LW.DIST; combine = min)cn = LW2.DST[isequal.(LW2.SOURCE, "CONNECTOR")] # the road nodes joined# The map zoomed to the two points, their connectors and their pathzm =unique([pW; cn])fig, ax =makemap(NW2.LON[zm], NW2.LAT[zm]; doRoadbkgd =false, xexpand =0.9, yexpand =0.12)resize!(fig, 500, 520)lines!(ax, x2ln(gw, NW.LON), x2ln(gw, NW.LAT); color = (ink, 0.6), linewidth =0.8)# The connectors, dashed: virtual roads that stand for streets not in FAF5for r ineachrow(LW2[isequal.(LW2.SOURCE, "CONNECTOR"), :])lines!(ax, [XY2[r.SRC, 1], NW2.LON[r.DST]], [XY2[r.SRC, 2], NW2.LAT[r.DST]]; color = bluec, linestyle =:dash, linewidth =1.5)end# The path: solid along roads, dashed along its two connectorslines!(ax, NW2.LON[pW[2:end-1]], NW2.LAT[pW[2:end-1]]; color = (supc, 0.8), linewidth =3)for s in (pW[1:2], pW[end-1:end])lines!(ax, NW2.LON[s], NW2.LAT[s]; color = (supc, 0.8), linewidth =3, linestyle =:dash)endscatter!(ax, NW2.LON[cn], NW2.LAT[cn]; color = supc, markersize =9)scatter!(ax, XY2[:, 1], XY2[:, 2]; color = bluec, markersize =14, strokecolor =:white, strokewidth =1)text!(ax, XY2[:, 1], XY2[:, 2]; text = ["1", "2"], fontsize =15, color = bluec, font =:bold, align = [(:left, :center), (:left, :bottom)], offset = [(9, 0), (7, 5)])fig
Figure 15: The roads of Wake County around two new points (blue, numbered as addconnectors numbers them), their connectors (dashed) to the road nodes they join (red dots), and the shortest path between them (red), dashed where it runs along a connector.
By road, the second point is 12.65 mi from the first, and the first 12.34 mi from the second, since one-way links make the two directions differ. The great circle distance is 9.03 mi, and 1.2 times it, 10.83 mi, falls short of the road distance over so short a trip, as Lecture 2.2, Sec. 7, found for distances under about 50 miles. Comparing the road distance with Google Maps’ 12.9 miles is the Landmark check of Lecture 1.2.
Example 2: Road distances between 100K NC cities
Determine the road distances between the ten North Carolina cities of over 100,000 people, those of Lecture 2.2, Sec. 6, and compare the distance from Cary to Charlotte with its estimate from the great circle distance.
# Code block 47: the NC cities over 100K, from the Logjam gazetteerdf =filter(r -> (r.STFIP ==st2fips(:NC)) && (r.POP >100_000),usplace())select!(df, :NAME, :LAT, :LON, :POP)prt(df)
NAME LAT LON POP
───────────────────────────────────────
Cary 35.78 -78.82 174,721
Charlotte 35.21 -80.83 874,579
Concord 35.39 -80.64 105,240
Durham 35.98 -78.90 283,506
Fayetteville 35.08 -78.97 208,501
Greensboro 36.10 -79.83 299,035
High Point 35.99 -79.99 114,059
Raleigh 35.83 -78.64 467,665
Wilmington 34.21 -77.89 115,451
Winston-Salem 36.10 -80.26 249,545
# Code block 48: their road distances, a 10 x 10 matrixn =nrow(df)Nc, Lc =cropnetwork(N, L, df.LON, df.LAT) # a box around the citiesNc, Lc =addconnectors(Nc, Lc, df.LON, df.LAT; add_nf_nf =false)fw, bw = Lc.DIR .!=-1, Lc.DIR .!=1# as in Code block 45g =SimpleWeightedDiGraph([Lc.SRC[fw]; Lc.DST[bw]], [Lc.DST[fw]; Lc.SRC[bw]], [Lc.DIST[fw]; Lc.DIST[bw]]; combine = min)D, P =shortestpaths(g, n) # from each city to eachround.(D, digits =1)
Code block 49 gives the distance from Cary (1) to Charlotte (2), by road and by 1.2 times the great circle distance, the circuity factor of Lecture 2.2, Sec. 7:
# Code block 49: Cary to Charlotte, by road and by great circle x 1.2@show D[1, 2]@showdgc([df.LON[1], df.LAT[1]], [df.LON[2], df.LAT[2]]) *1.2# Is every road distance at least its great circle distance?G = [dgc([df.LON[i], df.LAT[i]], [df.LON[j], df.LAT[j]]) for i in1:n, j in1:n]@showall(D .>= G);
No road distance falls below its great circle distance, which is a lower bound on it: the Bounds check of Lecture 1.2. The matrix is not quite symmetric, since one-way links make the distance from Charlotte back to Cary 131.1 mi.
Cary to Charlotte is 131.7 mi by road, against 143.9 mi estimated as 1.2 times the great circle distance
Show the code that draws this map
# Code block 50: the cities, their road network, and Cary to CharlottegL =SimpleWeightedGraph(Lc.SRC, Lc.DST, Lc.DIST; combine = min)p12 =tracepath(P[1], 1, 2) # Cary to Charlottefig, ax =makemap(Nc.LON, Nc.LAT; doRoadbkgd =false, xexpand =0.06, yexpand =0.03)resize!(fig, 687, 400)lines!(ax, x2ln(gL, Nc.LON), x2ln(gL, Nc.LAT); color = (ink, 0.3), linewidth =0.5)lines!(ax, Nc.LON[p12], Nc.LAT[p12]; color = (supc, 0.85), linewidth =3)scatter!(ax, df.LON, df.LAT; color = bluec, markersize =10)# aligntext's label placement, moved for the labels the route would crossal =aligntext(df.LON, df.LAT; offsetamt =6)A, O =copy(al[1].second), copy(al[2].second)fix =Dict("Cary"=> ((:center, :top), (0, -7)),"Concord"=> ((:left, :top), (6, -6)),"Charlotte"=> ((:left, :top), (6, -6)),"Wilmington"=> ((:right, :top), (-6, -6)))for (k, nm) inenumerate(df.NAME)haskey(fix, nm) && ((A[k], O[k]) = fix[nm])endtext!(ax, df.LON, df.LAT; text = df.NAME, fontsize =12, align = A, offset = O)fig
Figure 16: The road network in a box around North Carolina’s ten cities of over 100,000 people, with the shortest route by distance from Cary to Charlotte in red.
The ten cities get the same analysis as Wake County’s two points, on the road network around them, and Dijkstra’s procedure then gives a ten by ten matrix: only the distances between pairs of the ten cities are wanted, not those between the individual road nodes. From Cary, node 1, to Charlotte, node 2, the road distance, connectors included, is below 1.2 times the great circle distance, and in 2024 Google Maps gave 162 miles, higher still and about 30 miles more than the road distance found then. That is a significant difference, and worth asking why. Minimizing distance, the procedure takes the route in red in Fig. 16, which is the shortest route by distance. The routes Google Maps recommended, on I-40 or on US 64 and US 421, were both a lot longer, because Google Maps is looking at travel time: it finds the route with the shortest travel time and does not even consider the shortest distance, whereas the only weight given here was distance.
The procedure did what it is supposed to do: it minimized distance, using whichever roads in the network give the shortest route. To duplicate Google Maps, the distance would be translated into time, which would more likely give the I-40 route, and that is the subject of Sec. 5.
5. Travel time
From Cary to Charlotte, Google Maps recommended a much longer route than the shortest one of Ex. 2 because it was minimizing time: it could travel a lot faster on I-40 and I-85 than on US 64 through all of the smaller towns. Minimizing time instead of distance changes only the weight on each arc, and Logjam’s functions are built to allow it. prune_reindex and addconnectors require no fixed column names: they work with the first three columns of the node and link data frames by position, so that different data sources can be used without renaming every column. For the nodes, the first column is the index, followed by the longitude and the latitude; for the links, the first is the source node, the second the destination, and the third the weight, labeled here as distance. To minimize time, the third column would hold a time instead, labeled however is convenient, and pruning and the rest would still work; generically it could be called the weight.
The next step converts the road distances to travel times, since most often the path wanted is not the shortest in distance but the shortest in travel time. The conversion rests on an estimate of the average speed of travel along each link. In the 2018 version of this analysis the speeds were estimates of likely reasonable values, which could be adjusted, and which a function would take as parameters, set by how optimistic or pessimistic one is about actual travel times. Rural Interstates were taken at 75 mph, five above a posted limit of 70. Urban roads that are not Interstates were taken at 25 mph: their posted speed is probably 35 to 45 mph, and the 25 accounts for stopping at red lights. An alternative is to include an expected delay explicitly at every node passed through, for its red light, with a speed between nodes of perhaps 45 mph; taking an average speed that includes a reasonable amount of stopping time is the simpler choice. Each link’s time is its distance divided by its speed, and that time replaces the distance as the link’s weight. A connector, from a facility such as a distribution center, a manufacturing facility or a customer location, was given an average speed of about 20 mph. Connectors joining two facilities directly were ignored: as long as the road network is reasonably dense, travel from one facility to another can be assumed to use it. Only on a sparse network such as the Interstates alone, with areas of the country that have no Interstate, would they be relevant.
The other two speeds were 65 mph on urban Interstates and 50 mph on rural roads that are not Interstates.
An estimate of travel time over a long trip can be remarkably accurate, to within a few minutes even with heavy traffic and construction along the way, unless there is a large delay such as an accident. The reason is the law of large numbers. A long trip is made up of many small road segments, each an opportunity for its estimate to be a little over or a little under, and when the estimates of all of the segments are added up, the overestimates and underestimates tend to balance out. That holds only if the errors on the segments are independent. Segments far apart from each other can reliably be assumed to be, unless an accident were so large that its delay extended over the entire trip.
Connectors pose one problem for speed. Their distance is estimated, as the great circle distance times 1.3, the circuity taken for connectors, but there is no speed for them at all. From a home to a major road of the highway system, the trip is not going to average 60 mph down a residential street: the actual travel speed is perhaps 35 mph at best, and averaging in the time spent waiting at stoplights gives something like 25 to 30 mph. The same holds on highways that are not Interstates, such as US 64 heading west out of Raleigh, with all of its red lights. The network model could be extended to model the nodes themselves, with a cost at each node such as the average length of its signal; that is possible for an application that needs it, though it is not needed here.
FAF5 needs no table of speeds, because each of its links carries a free-flow travel time in each direction, AB_TIME from SRC to DST and BA_TIME back, in minutes. A few links carry no time, and for those the time comes from the link’s posted speed, SPEED; a connector, which has neither, is taken at 20 mph, as in 2018. Every link then has a time each way, and the travel times between the ten cities of Ex. 2 follow on the same network:
# Code block 51: travel times between the cities of Ex. 2, in minutescon =isequal.(Lc.SOURCE, "CONNECTOR")# A link's minutes: FAF5's free-flow time t, or, where a link has none,# 60 d / v at its posted speed v; a connector at 20 mphmins(t, v, d) = t >0 ? t :60d / vtAB = [c ? 60d /20:mins(t, v, d) for (c, t, v, d) inzip(con, Lc.AB_TIME, Lc.SPEED, Lc.DIST)]tBA = [c ? 60d /20:mins(t, v, d) for (c, t, v, d) inzip(con, Lc.BA_TIME, Lc.SPEED, Lc.DIST)]fw, bw = Lc.DIR .!=-1, Lc.DIR .!=1# as in Code block 45gT =SimpleWeightedDiGraph([Lc.SRC[fw]; Lc.DST[bw]], [Lc.DST[fw]; Lc.SRC[bw]], [tAB[fw]; tBA[bw]]; combine = min)T, PT =shortestpaths(gT, n) # minutesround.(T)
1
Only the times from city to city are wanted, so Dijkstra’s procedure is run from each of the n cities, the first n nodes, and only the first n entries of each result are kept, the n by n matrix T.
Dijkstra’s procedure then runs on the travel times from each city: in 2018, with 16 cities, a 16 by 16 matrix of the shortest travel times, computed in 0.3 s. The whole road network is used to get those times and then thrown away, keeping just the small matrix of times, as Code block 51 keeps the 10 by 10 matrix T. Dijkstra’s procedure also returns the predecessors, which are much larger, since they cover every node of the network, and from them the paths can be extracted. The shortest travel time route is not necessarily the route that is shortest in distance. So the distance of each fastest route cannot be found by running Dijkstra’s procedure on distance; instead, it is added up along the route, as Code block 52 does for every pair of cities.
# Code block 52: each route measured the other way, along its pathWd, Wt =weights(g), weights(gT) # miles and minutes, arc by arc# segcost gives the weight of each link along a pathDT = [i == j ? 0.0:sum(segcost(tracepath(PT[i], i, j), Wd)) for i in1:n, j in1:n] # miles on each fastest routeTD = [i == j ? 0.0:sum(segcost(tracepath(P[i], i, j), Wt)) for i in1:n, j in1:n] # minutes on each shortest routeir =findfirst(==("Raleigh"), df.NAME)prt(DataFrame("From"=> df.NAME[[1, ir]], "To"=> df.NAME[2],"Shortest, mi"=>round.(D[[1, ir], 2], digits =1),"its min"=>round.(TD[[1, ir], 2]),"Fastest, min"=>round.(T[[1, ir], 2]),"its mi"=>round.(DT[[1, ir], 2], digits =1)))
1
The path of the shortest travel time route from city i to city j is taken from the predecessors PT by tracepath, and the distances of its links are added along it, for every pair of cities, to give the distance of each shortest travel time route.
From To Shortest, mi its min Fastest, min its mi
─────────────────────────────────────────────────────────────────
Cary Charlotte 131.70 176 172 134.90
Raleigh Charlotte 143.50 184 178 169.70
For 85 of the 90 ordered pairs of cities, the fastest route is longer than the shortest. From Raleigh to Charlotte, the fastest route is 26.2 mi longer than the shortest and 7 min faster. From Cary to Charlotte, on FAF5’s free-flow times, the fastest route is only 3.1 mi longer than the shortest.
Example 3: High Point to Goldsboro
Determine the shortest route and the fastest route by road from High Point to Goldsboro, and the distance and travel time of each, taking FAF5’s free-flow times on its links and 20 mph on the connectors.
# Code block 53: High Point to Goldsboro, shortest and fastesthp =loc2lonlat("High Point"; state =:NC)gb =loc2lonlat("Goldsboro"; state =:NC)X, Y = [hp.LON, gb.LON], [hp.LAT, gb.LAT]N3, L3 =cropnetwork(N, L, X, Y; xexpand =0.2, yexpand =0.5)N3, L3 =addconnectors(N3, L3, X, Y; add_nf_nf =false)c3 =isequal.(L3.SOURCE, "CONNECTOR")t3AB = [c ? 60d /20:mins(t, v, d) # minutes, as in Code block 51 for (c, t, v, d) inzip(c3, L3.AB_TIME, L3.SPEED, L3.DIST)]t3BA = [c ? 60d /20:mins(t, v, d) for (c, t, v, d) inzip(c3, L3.BA_TIME, L3.SPEED, L3.DIST)]fw, bw = L3.DIR .!=-1, L3.DIR .!=1I3, J3 = [L3.SRC[fw]; L3.DST[bw]], [L3.DST[fw]; L3.SRC[bw]]g3 =SimpleWeightedDiGraph(I3, J3, [L3.DIST[fw]; L3.DIST[bw]]; combine = min) # milesg3T =SimpleWeightedDiGraph(I3, J3, [t3AB[fw]; t3BA[bw]]; combine = min) # minutespD =tracepath(dijkstra_shortest_paths(g3, 1).parents, 1, 2) # shortestpF =tracepath(dijkstra_shortest_paths(g3T, 1).parents, 1, 2) # fastestmiD =sum.(segcost.([pD, pF], Ref(weights(g3)))) # miles of eachminT =sum.(segcost.([pD, pF], Ref(weights(g3T)))) # minutes of each# hm, defined with the page's setup, writes minutes as hours and minutesprt(DataFrame("Route"=> ["Shortest", "Fastest"],"mi"=>round.(miD, digits =1), "Time"=>hm.(minT)))
Route mi Time
──────────────────────────────
Shortest 137.60 2 h 55 min
Fastest 151.00 2 h 44 min
Shortest route 137.6 mi in 2 h 55 min; fastest route 151.0 mi in 2 h 44 min
Show the code that draws this map
# Code block 54: the two routes, High Point to Goldsborog3u =SimpleWeightedGraph(L3.SRC, L3.DST, L3.DIST; combine = min)fig, ax =makemap(N3.LON, N3.LAT; doRoadbkgd =false, xexpand =0.02, yexpand =0.04)resize!(fig, 687, 360)lines!(ax, x2ln(g3u, N3.LON), x2ln(g3u, N3.LAT); color = (ink, 0.3), linewidth =0.5)lD =lines!(ax, N3.LON[pD], N3.LAT[pD]; color = (supc, 0.85), linewidth =3)lF =lines!(ax, N3.LON[pF], N3.LAT[pF]; color = (dcc, 0.85), linewidth =3)scatter!(ax, X, Y; color = bluec, markersize =11)text!(ax, X, Y; text = ["High Point", "Goldsboro"], fontsize =12, align = [(:right, :bottom), (:right, :top)], offset = [(-6, 6), (-8, -4)])axislegend(ax, [lD, lF], ["Shortest", "Fastest"]; position =:rt, labelsize =12)fig
Figure 17: High Point to Goldsboro on the road network: the shortest route (red) and the fastest route (green).
The new nodes are numbered first, so the cities come first in the arc list and every road node has a higher number, which is what lets the path between two cities be read directly. The 2018 version of this example, on the Oak Ridge roads, drew the fastest route from High Point to Goldsboro on the map with the distance and the estimated time in its title. Instead of going through downtown Raleigh, which the short-distance path would do, the fastest route goes around the Beltline, then continues along I-40 to US 70 and over to Goldsboro. It was estimated at 143.5 miles and 2 hours 28 minutes. Google Maps then gave 2 hours 33 minutes for its shortest travel time route, within five minutes of the estimate, at 144 miles, almost exactly the same distance. The estimated times were probably close to what Google Maps uses, even though the speeds were never designed to match it: an assumed average speed must be fairly close to Google’s own. A match is not always this close, but this is the underlying mechanics of what Google Maps does when it solves a problem like this one.
# Code block 55: a landmark, OSRM's routes on OpenStreetMap datausingHTTP, JSON3a =loc2lonlat("High Point"; state =:NC) # its LON and LATb =loc2lonlat("Goldsboro"; state =:NC)url ="https://router.project-osrm.org/route/v1/driving/"*"$(a.LON),$(a.LAT);$(b.LON),$(b.LAT)"*"?alternatives=true&overview=false"r = JSON3.read(HTTP.get(url).body) # the router's answerosrm =DataFrame(DIST = [rt.distance for rt in r.routes] ./1609.344, TIME = [rt.duration for rt in r.routes] ./60)
1
OSRM returns each route’s distance in meters and its duration in seconds, so dividing by 1,609.344 m to the mile and 60 s to the minute gives miles and minutes.
Table 4: OSRM’s routes from High Point to Goldsboro, fastest first, as the request of Code block 55 returned them on 2026-10-09.
Lecture 2.7, Sec. 3, anticipated an application of set covering, Model 4 in Lecture 2.7, to a road network: the fewest intersections such that every intersection is within a stated reach of a selected one, the selected intersections being truck terminals. Locating transshipment terminals adds one condition, that each terminal lie within a day’s drive of another so that a driver can carry freight between them, and with that condition in the formulation from the start the problem is Model 2, a special case of set covering.
minimize: number of terminals
solve for:
(a) intersections to select as terminals.
subject to:
(a) Coverage: every intersection is within half a day’s drive of a terminal; (b) Linehaul: every terminal is within a day’s drive of another terminal.
assumptions:
(a) travel time between two intersections is the time of the shortest path between them;
(b) any intersection can be a terminal.
=\begin{cases} 1, & \text{if intersection } i \text{ is a terminal} \\ 0, & \text{otherwise} \end{cases}
t_{ij}
=travel time from intersection i to intersection j (h)
t_{\max}
=longest time a driver may drive in a day (h)
N
=intersections
=\{1, \ldots, n\}.
Model 2 states two constraints, and both are rows. Row (a), coverage, is the row of Eq. 9 in Lecture 2.7 with every cost one, intersection j being covered by a terminal at i when t_{ij} \le t_{\max}/2, half a day’s drive. Row (b), linehaul, is the row Model 2 adds: an intersection can be selected as a terminal only if another terminal is selected within t_{\max} of it. Its sum is over the set \{k \in N \setminus \{i\} : t_{ik} \le t_{\max}\}, read “the k in N other than i such that t_{ik} \le t_{\max}”: the colon reads “such that,” as it does under the sum of row (a), and N \setminus \{i\} is the set difference of N and \{i\}, the members of N that are not in \{i\}, so the sum counts the terminals within t_{\max} of i other than i itself, and the implementation writes the set difference as setdiff(N, [i]). Leaving i out matters: since t_{ii} = 0, the sum would otherwise include x_i, and the row would hold whatever was selected.
Model 2 implementation: Transshipment terminals
# Model: transshipment terminals, set covering with a linehaul conditionusingJuMP, HiGHSfunctionterminals(T, tmax) N =1:size(T, 1) model =Model(HiGHS.Optimizer)@variable(model, x[i ∈ N], Bin)@objective(model, Min, sum(x[i] for i ∈ N))@constraint(model, coverage[j ∈ N], # (a)sum(x[i] for i ∈ N if T[i, j] ≤ tmax /2) ≥1)@constraint(model, linehaul[i ∈ N], # (b) x[i] ≤sum(x[k] for k ∈setdiff(N, [i])if T[i, k] ≤ tmax))set_silent(model)optimize!(model)returnfindall(==(1.0), snapvals(value.(x)))end
terminals (generic function with 1 method)
Example 4: Transshipment terminals
A national less-than-truckload (LTL) carrier would like to know what would be the minimum number of transshipment terminals needed to cover the continental U.S. and their locations. All of their existing local pickup and delivery terminals are located near the intersection of two or more Interstate highways. Since hours-of-service (HOS) requirements limit a driver to eleven hours of driving per day (not including loading/unloading time), they would like each of their local terminals to be within five and a half hours travel time of a transshipment terminal (so that one driver can complete one roundtrip in a day) and, for linehaul operations between the transshipment terminals, they would like for all of the transshipment terminals to be within eleven hours of each other.
This is Model 2 with t_{\max} = 11 h, on the Interstate network thinned in Sec. 4, with travel on it taken at an average of 60 mph.
The solution takes three steps: find the intersections, compute the travel time between every pair of them, and solve Model 2, after which both of its conditions are checked on the terminals it selects.
An intersection is a node of the thinned Interstate network of Code block 35 that is not the end of an Interstate, a node of degree 1. Removing an end can leave a new end behind, so the ends are removed until none is left, and the network is then thinned again:
# Code block 56: the Interstate intersectionsiN, iL = tN, tL # the thinned network of Code block 35whiletrue# remove the ends until none is left deg = [count(==(k), [iL.SRC; iL.DST]) for k in1:nrow(iN)]all(>(1), deg) &&break# no node of degree 1 is leftglobal iN, iL =prune_reindex(iN, filter(r -> deg[r.SRC] >1&& deg[r.DST] >1, iL))endiN, iL =prune_reindex(thin(iN, iL)...) # remove new degree-2 nodesnrow(iN)
1
A while loop has a scope of its own, so a variable assigned inside it is local to the loop unless it is declared otherwise. global makes this line update the iN and iL defined above the loop, so that each pass starts from the network the last pass left. Without it, in a script run as a file, such as the companion .jl, Julia would warn that the assignment is ambiguous, treat iN and iL as new variables local to the loop, and stop with an error at the loop’s first line, which reads iN before the loop has given it a value. Typed at the REPL, the loop would update them without the keyword, which does no harm there.
593
Dividing each link’s miles by 60 mph gives its hours, and the shortest paths then give the hours between every pair of intersections, the t_{ij} of Eq. 5. The implementation of Model 2 then selects the terminals for an eleven-hour day, which puts every intersection within five and a half hours of a terminal, so that one driver can drive there and back:
# Code block 57: travel times, then the terminals for an 11-hour daygI =SimpleWeightedGraph(iL.SRC, iL.DST, iL.DIST ./60; combine = min)T, _ =shortestpaths(gI, weights(gI), nv(gI)) # hours, all pairsidx =terminals(T, 11)length(idx)
1
The last argument is how many of the first nodes Dijkstra’s procedure is run from, and nv(gI), the number of nodes, runs it from every node, so T holds the hours between every pair of intersections.
17
The two constraints are hard limits, so the terminals selected can be given the Bounds check of Lecture 1.2: how far any intersection lies from its nearest terminal, against five and a half hours for constraint (a), and how far any terminal lies from its nearest other terminal, against eleven for constraint (b):
# Code block 58: checking constraints (a) and (b) on the terminalstoterm = [minimum(T[i, j] for i in idx) for j in1:nv(gI)]near = [minimum(T[a, b] for b in idx if b != a) for a in idx]round(maximum(toterm), digits =2), round(maximum(near), digits =2)
(5.47, 10.77)
Every intersection is within 5.47 h of a terminal, and every terminal within 10.77 h of another.
17 transshipment terminals put every one of the 593 Interstate intersections within five and a half hours, and each terminal within eleven hours of another
Show the code that draws this map
# Code block 59: the network and its terminalsfig, ax =makemap(iN.LON, iN.LAT; doRoadbkgd =false, xexpand =0.02, yexpand =0.02)resize!(fig, 687, 430)lines!(ax, x2ln(gI, iN.LON), x2ln(gI, iN.LAT); color = ink, linewidth =0.6)scatter!(ax, iN.LON[idx], iN.LAT[idx]; color =:transparent, strokecolor = supc, strokewidth =2, markersize =14)fig
Figure 18: The Interstate network after its ends are dropped, and the transshipment terminals (red) that put every intersection within five and a half hours of one and each terminal within eleven hours of another.
R. K. Ahuja, T. L. Magnanti and J. B. Orlin, Network Flows: Theory, Algorithms, and Applications, Prentice Hall, 1993, Sec. 3.2, “Complexity Analysis.”↩︎
P. E. Hart, N. J. Nilsson and B. Raphael, “A formal basis for the heuristic determination of minimum cost paths,” IEEE Transactions on Systems Science and Cybernetics, 4(2), 1968, pp. 100–107, where A* was introduced.↩︎
S. J. Russell and P. Norvig, Artificial Intelligence: A Modern Approach, 4th ed., Pearson, 2021, Sec. 3.5.2, “A* search.”↩︎
H. Bast, D. Delling, A. Goldberg, M. Müller-Hannemann, T. Pajor, P. Sanders, D. Wagner and R. F. Werneck, “Route planning in transportation networks,” in Algorithm Engineering, Springer, 2016, pp. 19–80, Secs. 2.2 and 6.↩︎
Bureau of Transportation Statistics, National Transportation Atlas Database, FAF5 Network Links, data dictionary doi:10.21949/1529027.↩︎
