This file shows how to solve the shortest path problem for the SEERVADA PARK prototype example from Hillier and Lieberman, Section 9.1.
from pyomo.environ import *
from pyomo.opt import *
opt = solvers.SolverFactory("glpk")
We define the node set and the distance values on each arc of the network:
N = ['O', 'A', 'B', 'C', 'D', 'E', 'T']
D = {('O','A'):2,
('O','B'):5,
('O','C'):4,
('A','D'):7,
('A','B'):2,
('B','D'):4,
('B','E'):3,
('B','C'):1,
('C','E'):4,
('D','E'):1,
('D','T'):5,
('E','T'):7}
We will need the list of all arcs later. We could hand-code it, but it's easier to extract it from D
as the list of key values.
A = list(D.keys())
The SEERVADA PARK road network as described in the example is an undirected network. However, in the linear programming formulation of the shortest path problem, we can only use directed arcs. Thus, for each arc in D
which we listed above, we need to add an arc of the same distance with opposite orientation. To save us some writing, we add the backlinks with two lines of Python code:
for (i,j) in A:
D[(j,i)]=D[(i,j)]
print(D) # Just to check what our new D looks like
{('O', 'A'): 2, ('O', 'B'): 5, ('O', 'C'): 4, ('A', 'D'): 7, ('A', 'B'): 2, ('B', 'D'): 4, ('B', 'E'): 3, ('B', 'C'): 1, ('C', 'E'): 4, ('D', 'E'): 1, ('D', 'T'): 5, ('E', 'T'): 7, ('A', 'O'): 2, ('B', 'O'): 5, ('C', 'O'): 4, ('D', 'A'): 7, ('B', 'A'): 2, ('D', 'B'): 4, ('E', 'B'): 3, ('C', 'B'): 1, ('E', 'C'): 4, ('E', 'D'): 1, ('T', 'D'): 5, ('T', 'E'): 7}
After augmenting the distance directory, we also need to update the list of arcs.
A = list(D.keys())
Now we can set up the model and the decision variable x
, the flow on the network:
model = ConcreteModel()
model.x = Var(A, within=NonNegativeReals)
For the shortest path problem, we assume an inflow of 1 at the source node O
and an outflow of 1 at destination node T
. At all other nodes, inflow equals outflow.
def flow_rule(model, n):
InFlow = sum(model.x[i,j] for (i,j) in A if j==n)
OutFlow = sum(model.x[i,j] for (i,j) in A if i==n)
if n == 'O':
return OutFlow - InFlow == 1 # b at origin = 1
elif n == 'T':
return OutFlow - InFlow == -1 # b at target = -1
else:
return OutFlow - InFlow == 0 # b at all other nodes = 0
model.flow = Constraint(N, rule=flow_rule)
The objective is minimizing distance traveled.
model.distance = Objective(expr = sum(model.x[a]*D[a] for a in A), sense=minimize)
Now we solve and look at the results:
results = opt.solve(model)
model.x.get_values()
{('O', 'A'): 1.0, ('O', 'B'): 0.0, ('O', 'C'): 0.0, ('A', 'D'): 0.0, ('A', 'B'): 1.0, ('B', 'D'): 1.0, ('B', 'E'): 0.0, ('B', 'C'): 0.0, ('C', 'E'): 0.0, ('D', 'E'): 0.0, ('D', 'T'): 1.0, ('E', 'T'): 0.0, ('A', 'O'): 0.0, ('B', 'O'): 0.0, ('C', 'O'): 0.0, ('D', 'A'): 0.0, ('B', 'A'): 0.0, ('D', 'B'): 0.0, ('E', 'B'): 0.0, ('C', 'B'): 0.0, ('E', 'C'): 0.0, ('E', 'D'): 0.0, ('T', 'D'): 0.0, ('T', 'E'): 0.0}
model.distance.expr()
13.0
So the shortest path is O-A-B-D-T
with a distance of 13. As you will recognize immediately in this simple example, there is an alternative shortest path O-A-B-E-D-T
; Pyomo does not report such alternative solutions.