Computers & Operations Research 40 (2013) 378–384

Contents lists available at SciVerse ScienceDirect

Computers & Operations Research

journal homepage: www.elsevier.com/locate/caor

On an exact method for the constrained shortest path problem

Leonardo Lozano, Andrés L. Medaglia n
Centro para la Optimización y Probabilidad Aplicada (COPA), Departamento de Ingenierı́a Industrial, Universidad de los Andes, Cr 1E No. 19A-10, ML711, Bogotá, Colombia

a r t i c l e i n f o a b s t r a c t

Available online 16 July 2012 The constrained shortest path (CSP) is a well known NP-Hard problem. Besides from its straightforward
Keywords: application as a network problem, the CSP is also used as a building block under column-generation
Constrained shortest path solution methods for crew scheduling and crew rostering problems. We propose an exact solution
Shortest path problem method for the CSP capable of handling large-scale networks in a reasonable amount of time. We
Column generation compared our approach with three different state-of-the-art algorithms for the CSP and found optimal
solutions on networks with up to 40,000 nodes and 800,000 arcs. We extended the algorithm to
effectively solve the auxiliary problems of a multi-activity shift scheduling problem and a bus rapid
transit route design problem tackled with column generation. We obtained significant speedups against
alternative column generation schemes that solve the auxiliary problem with state-of-the-art
commercial (linear) optimizers. We also present a first parallel version of our algorithm that shows
promising results.
& 2012 Elsevier Ltd. All rights reserved.

1. Introduction or label-correcting algorithms. Perhaps the most distinctive

feature of these algorithms is that they can be extremely fast
Consider a directed graph G ¼ ðN ,AÞ, where N ¼ fv1 , . . . , for reasonably sized networks; however, they might fail to scale
vi , . . . ,vn g represents the set of nodes and well for very large networks due to the ‘‘curse of dimensionality’’
 
A ¼ ði,jÞ9vi A N , vj A N ,ia j the set of arcs. For each arc ði,jÞ A A of DP, more specifically, the number of labels that need to be
the nonnegative weights, cij and tij, represent the cost and stored can be huge and, consequently, impossible to manage. One
resource consumption incurred by traversing the arc. The con- of the first algorithms in this category is the one proposed by
strained shortest path (CSP) problem consists of finding a path P Joksch [22], lately extended by Dumitrescu and Boland [13] by
from a start node vs A N to an end node ve A N that minimizes the including preprocessing techniques. More recently, Zhu and
total cost, subject to not exceeding a maximum resource con- Wilhelm [34,35] developed a three-stage approach specialized
sumption T. The CSP is known to be NP-Hard even for the case of for column generation subproblems that transforms the CSP into
one resource [16]. a shortest path problem and solves it using a labeling method.
Aside from its straightforward application (i.e., shortest path The second major solution approach for the CSP is based on the
subject to a time limit), the CSP and its variants naturally arise as use of Lagrangian relaxation to solve the integer programming
auxiliary problems in column generation schemes for air cargo formulation of the problem. The efficiency of this approach relies
planning and routing [11], flight planning [18], crew pairing on the effectiveness of the underlying unconstrained shortest
[23,25], tail assignment problem in aircraft scheduling [20], path algorithms. Handler and Zang [21] solved the CSP by using a
day-to-day crew operation [32], and crew rostering problems k-th shortest path algorithm; that is, they identify k paths, sort
[15], among others. them by length, and evaluate them successively until they find
Some variants of the CSP include non-additive objective the first path that satisfies the (side) constraint. The same authors
functions [27], probability functions defined over the network also proposed a Lagrangian relaxation approach that can be seen
[2], forbidden paths [33], and replenishment arcs [30]. as a generalization of the k-shortest path algorithm for the CSP,
Solution strategies for the CSP can be classified into one of two that is, a method that searches in a direction that (linearly)
main categories: Dynamic Programming (DP) and Lagrangian combines the information of the original objective and the (side)
relaxation. Methods based on DP are also known as label-setting constraint. Santos et al. [29] extended this idea by refining the
search direction based on the relative tightness of the (side)
constraint. More recently, Carlyle et al. [5] used Lagrangian
n
relaxation combined with enumeration of near-shortest paths.
Corresponding author. Tel.: þ57 1 3394949x2880; fax: þ57 1 3324321.
E-mail address: [email protected] (A.L. Medaglia). This research aims to expand the body of knowledge of the
URL: http://wwwprof.uniandes.edu.co/~amedagli (A.L. Medaglia). CSP. Specifically, it proposes an exact algorithm that consistently

0305-0548/$ - see front matter & 2012 Elsevier Ltd. All rights reserved.
http://dx.doi.org/10.1016/j.cor.2012.07.008
 L. Lozano, A.L. Medaglia / Computers & Operations Research 40 (2013) 378–384 379

outperforms three state-of-the-art algorithms, namely, those by Algorithm 2. Pulse function.

Santos et al. [29], Zhu and Wilhelm [35], and Dumitrescu and
Input: vk, c, t, P
Boland [13], over a vast set of instances from the literature. In
Output: void
addition and most importantly, the proposed approach provides a
1: if checkDominanceðvk ,c,tÞ ¼ false then fsee Section 3:1g
flexible solution framework that can be easily parallelized and
2: if checkFeasibilityðvk ,tÞ ¼ true then fsee Section 3:2g
extended to include multiple side constraints and negative costs.
3: if checkBoundsðvk ,cÞ ¼ false then fsee Section 3:3g
This paper is organized as follows. Section 2 presents a broad
overview of the proposed pulse algorithm, focusing on its intui-
4: P 0 ’P [ fvk g
tion. Section 3 provides detailed description of the algorithm’s 5: for vi A N ðvk Þ do
components. Section 4 presents our computational experience 6: c0 ’c þ cki
with the proposed algorithm against the state-of-the-art algo- 7: t 0 ’t þ t ki
rithms for the CSP. Section 5 illustrates our experience while 8: pulseðvi ,c0 ,t 0 ,P 0 Þ
embedding the algorithm within column generation schemes. 9: end for
Finally, Section 6 concludes the paper and outlines future 10: end if
extensions. 11: end if
12: end if

2. The pulse algorithm: intuition and overview Note that the way the pulse algorithm explores the graph
follows a depth-first search (DFS) exploration scheme truncated
The idea behind the pulse algorithm is very simple, almost by the pruning strategies. An illustrative animation of the algo-
naive, yet very powerful. The algorithm is based on the idea of rithm is available at http://hdl.handle.net/1992/1144.
propagating pulses through a network from a start node vs A N to
an end node ve A N . As a pulse traverses the network from node to
node, it builds a partial path P including the nodes already 3. Pruning strategies
visited, the cumulative objective function cðPÞ and the cumulative
resource consumption tðPÞ. Each pulse that reaches the final node The pulse algorithm heavily depends on the strength of its
ve contains all the information for a feasible path P from vs to ve. pruning strategies. This section outlines the different strategies
If nothing prevents the pulses from propagating, the algorithm implemented for the CSP. It is noteworthy that given that all
completely enumerates all possible paths from vs to ve, ensuring weights defined over the arcs are nonnegative, the optimal
that the optimal path P n is always found. At the core of the path must be an elementary path that contains no cycles.
algorithm lies the ability to (effectively and aggressively) prune Whenever the function pulse is invoked on a node vk, the
pulses (i.e., prevent their propagation), without jeopardizing the algorithm verifies if visiting node vk creates a cycle in the partial
optimal path. This idea is shared with other algorithms like path P associated with the pulse. This verification is performed in
branch and bound, where an implicit enumeration is performed constant time using a binary vector that keeps a record of the
with relative efficiency. Similarly, the strength of the pulse nodes already visited. Paths containing cycles are not considered
algorithm depends on the pruning strategies; in our case for the in the search.
CSP, pruning by dominance, bounds, and infeasibility.
Algorithm 1 presents the pseudocode of the pulse algorithm.
Lines 1 through 3 initialize the partial path P, the initial objective 3.1. Pruning by dominance
function c0, and the initial resource consumption t0. Line 4 runs a
one-to-all shortest path algorithm to find the minimum resource Dominance relationships can be defined over partial paths as
consumption for every node (see Section 3.2). Line 5 invokes the done in DP-based algorithms that extend Bellman-Ford’s cf. [24].
recursive function pulse starting the propagation from node vs. Let P 1 and P 2 be two partial paths to a given node vk. If
Line 6 returns the optimal path P n found in the recursion. cðP 1 Þ 4 cðP 2 Þ and tðP 1 Þ 4tðP 2 Þ then P 1 is said to be strongly
dominated by P 2 . If cðP 1 Þ 4 cðP 2 Þ and tðP 1 Þ ZtðP 2 Þ or
Algorithm 1. Pulse algorithm. cðP 1 Þ Z cðP 2 Þ and tðP 1 Þ 4 tðP 2 Þ then P 1 is said to be (weakly)
Input: G, vs, ve, T dominated by P 2 . If cðP 1 Þ ¼ cðP 2 Þ, tðP 1 Þ ¼ tðP 2 Þ, and P 1 a P 2 then
P 1 and P 2 are alternative paths to node vk.
Output: Optimal path P n
Given these relations, we use labels to check for strong
1: P’fg
dominance, (weak) dominance, or the existence of alternative
2: c0 ’0
paths. For each node vk A N we define a limited memory
3: t 0 ’0  
Lðvk Þ ¼ ðckl ,t kl Þ9l ¼ 1, . . . ,Q where ckl and tkl are the cost and
4: initializationðG,TÞ fsee Section 3:3g
resource consumption labels, respectively; and Q denotes the
5: pulseðvs ,c0 ,t 0 ,PÞ
total number of labels (i.e., memory size). At first glance, it seems
6: return P n reasonable to store all the values for cðPÞ and tðPÞ corresponding
to the paths that have passed through a given node vk, so that it is
Algorithm 2 shows the body of the recursive function pulse, possible to prune any incoming pulse in case it is dominated.
 
where N ðvk Þ ¼ vi A N 9ðk,iÞ A A is the set of head nodes of the However, this exhaustive approach would require huge amounts
outgoing arcs from node vk. Each time the function pulse is of memory to store the values for all partial paths in large
invoked on the final node ve, the information for the best path networks, ultimately, forcing the method to fail in terms of
known is stored globally and the pulse is no longer propagated. scalability. In contrast, by choosing a limited memory of Q labels
Lines 1 through 3 use the different pruning strategies in order to (per node), it is possible that some dominated paths may
prune the incoming pulse. If the pulse is not pruned, line 4 adds inefficiently pass through node vk, but still the optimal solution
the current node to the partial path and lines 5 through 9 propa- is never pruned (just dominated paths).
gate the pulse by recursively invoking the function pulse on all In contrast to a classical DP approach, it is noteworthy that
nodes vi A N ðvk Þ. these labels are not constantly expanding, nor the optimality
380 L. Lozano, A.L. Medaglia / Computers & Operations Research 40 (2013) 378–384

depends on storing them all. We evaluated the impact of different known objective function value. Given that all costs are nonnega-
storage rules and values of Q on the computational performance tive, if a partial path P exceeds the primal bound, that is cðPÞ Zc,
of the pulse algorithm. We tried the following storage rules: (1) to this path can be safely pruned because another path with a lower
randomly replace a pair of labels each time a new nondominated (or equal) total cost has already been found. Furthermore, we
path is found (random rule); (2) to always keep the labels with the derived a bound on the minimum cost cðvk Þ from any node vk to
lowest cost and resource consumption (elitist rule); and (3) to use the final node ve, following a similar procedure to the one used for
a distance metric to rank the labels and keep always a diversified tðvk Þ. Again, it is necessary to reverse the network and use a one-
(heterogeneous) label set (diversified rule). For the latter rule, we to-all shortest path algorithm to minimize the cost from the first
have borrowed the crowding distance from NSGA-II [9], where it node v0s ¼ ve in the reversed network to any other node vk; this is
is used as a parameter-less approach to preserve a diversified equivalent to find the minimum cost cðvk Þ from any node vk to ve
population in the approximate Pareto frontier. We defer until in the original network. As this shortest path algorithm does not
Section 4 the discussion regarding the computational efficiency of take into account the resource consumption, the values cðvk Þ are a
these rules. In summary, these labels prune dominated partial dual bound on the real minimum cost subject to the resource
paths, yet they consume a reasonable (and limited) amount of constraint. Given a partial path P to node vk, if cðPÞ þ cðvk Þ Zc,
memory. then path P can be safely pruned because the algorithm has
already found a complete feasible path that is better (or equal)
3.2. Pruning by infeasibility than any complete path that includes the partial path P. The
procedure to find all the dual bounds cðvk Þ is done as part of the
The idea of this strategy is to prune pulses at an early stage, as initialization of the algorithm. The use of the dual bound strength-
soon as it is known that the traveling pulse will not be able to ens the algorithm at the expense of a single run of a one-to-all
reach the final node meeting the (side) constraint. To do so, it is shortest path algorithm. It is worth mentioning, that similar
necessary to calculate the minimum resource consumption from bounds have been used by An-based shortest path algorithms [26].
any node vk to the final node ve. To obtain this information
we reverse the original network (i.e., change the direction of
all arcs) creating 0 0 4. Computational experiments for the CSP
  a new directed graph G ¼ ðN ,A Þ where
A0 ¼ ðj,iÞ9ði,jÞ A A and the weights are set to c0ji ¼ cij and t 0ji ¼ t ij
for ðj,iÞ A A0 . The new start node is the original end node v0s ¼ ve . In this section we compare the pulse algorithm against the
Next, we run a one-to-all shortest path algorithm to find the algorithms by Santos et al. [29], Zhu and Wilhelm [35], and
minimum resource consumption from the initial node to all other Dumitrescu and Boland [13] on a wide range of benchmark
nodes in graph G0 . This is equivalent to find the minimum problems from the literature. The proposed algorithm was imple-
resource consumption tðvk Þ from any node vk to the final node mented in Java, compiled using Eclipse SDK version 3.4.2, and the
ve in the original network. Then, we set the maximum resource experiments executed on a computer with an Intel Mobile Core
consumption tðvk Þ for a partial path up to any given node vk to 2 Duo @ 2.4 GHz CPU P8600 with 512 MB of RAM allocated to the
tðvk Þ ¼ Ttðvk Þ. That is, if for a partial path P to node vk it is memory heap size of the Java Virtual Machine on Windows XP
known that tðPÞ 4 tðvk Þ, then partial path P is infeasible because Professional. The shortest path problems were solved using our
tðPÞ þ tðvk Þ 4 T. The procedure to find all the tðvk Þ values is done Dijkstra’s double buckets implementation adapted from Cher-
as part of the initialization of the algorithm. kassky et al. [6].
This pruning strategy shares the same spirit of similar algo- Regarding the pruning by dominance strategy, a mixture
rithms that have been used to reduce the size of the original between the random and elitist storage rules (with three labels)
network, but it does it in a much simpler way. In contrast to the was proven to be both fast and effective. At each node we keep a
preprocessing phases of the algorithms by Dumitrescu and Boland first nondominated label that is overwritten every time that a
[13] and Zhu and Wilhelm [35], where they use iterative proce- partial path exhibits a lower cost than the one stored, a second
dures that calculate forward and backward shortest paths for all nondominated label that is overwritten every time that a partial
nodes and modify the original graph by deleting nodes and arcs, path exhibits a lower resource consumption than the one stored,
we only use a backward one-to-all shortest path algorithm for and a third nondominated label that is replaced randomly.
obtaining the lower bounds on the resource consumption. Gro-
micho et al. [19] also propose a similar idea to restrict the state 4.1. Experiments on large-scale instances (with one side constraint)
space in a DP context by using a feasibility check before extending
the labels, but without the look-ahead component of our lower We used the exact same benchmark problems used by Santos
bounds. et al. [29]; a testbed that comprises 180 instances organized in
The strength of this pruning strategy relies on the fact that it three groups. The first group contains the smaller instances with
prunes infeasible partial paths at early stages of the exploration. 9N 9 ¼ 10,000 and 9A9 ¼ 15,000, 25,000, 50,000, 100,000, 150,000,
For problems where the resource constraint is very tight, the and 200,000; the second group contains the middle-sized
algorithm takes advantage of this (early) aggressive pruning to instances with 9N 9 ¼ 20,000 and 9A9 ¼ 30,000, 50,000, 100,000,
explore the network fast, yet ensuring that the optimal solution is 200,000, 300,000, and 400,000; and the third group contains the
never pruned (only infeasible paths). larger instances with 9N 9 ¼ 40,000 and 9A9 ¼ 60,000, 100,000,
200,000, 400,000, 600,000, and 800,000. For each 9N 99A9 com-
3.3. Pruning by bounds bination in each group, there are 10 different instances. Further-
more, Santos et al. [29] define the maximum resource
This pruning strategy is inspired on the power of primal and consumption T based on the tightness of the constraint denoted
dual bounds in branch and bound. These bounds store the best by p and defined as follows:
objective function found (primal bound) and the most optimistic
TtðP nt Þ
information made available at the unexplored nodes (dual p¼
tðP nc ÞtðP nt Þ
bound).
Every time a pulse reaches the last node ve, a feasible path is where P nc and P nt are the paths that minimize cost and resource
generated and the primal bound c may be updated with the best consumption, respectively; and tðPÞ is the resource consumption
 L. Lozano, A.L. Medaglia / Computers & Operations Research 40 (2013) 378–384 381

Table 1
Computational results for the Santos et al. [29] instances.

Nodes Arcs p ¼ 0.1 p ¼0.2 p ¼ 0.4 p ¼0.6 p¼ 0.8

Pulse LRA Speedup Pulse LRA Speedup Pulse LRA Speedup Pulse LRA Speedup Pulse LRA Speedup

10,000 15,000 0.01 0.60 60.00 0.01 0.60 60.00 0.01 0.60 60.00 0.01 0.60 60.00 0.01 0.60 60.00
10,000 25,000 0.02 0.70 35.00 0.02 0.70 35.00 0.02 0.70 35.00 0.02 0.70 35.00 0.02 0.70 35.00
10,000 50,000 0.02 1.10 55.00 0.02 1.10 55.00 0.02 1.10 55.00 0.02 1.10 55.00 0.02 1.10 55.00
10,000 100,000 0.04 1.80 45.00 0.05 1.80 36.00 0.04 1.80 45.00 0.05 1.80 36.00 0.04 1.80 45.00
10,000 150,000 0.07 2.50 35.71 0.06 2.50 41.67 0.07 2.50 35.71 0.07 2.60 37.14 0.06 2.80 46.67
10,000 200,000 0.07 2.70 38.57 0.14 2.70 19.29 0.16 2.70 16.88 0.10 2.80 28.00 0.08 2.80 35.00

20,000 30,000 0.02 1.20 60.00 0.02 1.20 60.00 0.02 1.20 60.00 0.02 1.20 60.00 0.02 1.20 60.00
20,000 50,000 0.03 1.50 50.00 0.03 1.50 50.00 0.03 1.50 50.00 0.03 1.50 50.00 0.03 1.50 50.00
20,000 100,000 0.06 2.30 38.33 0.06 2.30 38.33 0.06 2.30 38.33 0.06 2.30 38.33 0.06 2.40 40.00
20,000 200,000 0.11 3.70 33.64 0.09 3.70 41.11 0.09 3.70 41.11 0.09 3.80 42.22 0.12 3.90 32.50
20,000 300,000 0.17 4.80 28.24 0.14 4.80 34.29 0.15 4.90 32.67 0.15 4.90 32.67 0.19 5.00 26.32
20,000 400,000 0.17 5.90 34.71 0.22 5.90 26.82 0.28 5.90 21.07 0.23 6.00 26.09 0.18 6.20 34.44

40,000 60,000 0.05 2.40 48.00 0.05 2.40 48.00 0.05 2.40 48.00 0.05 2.40 48.00 0.05 2.40 48.00
40,000 100,000 0.08 3.00 37.50 0.08 3.00 37.50 0.08 3.00 37.50 0.08 3.10 38.75 0.08 3.10 38.75
40,000 200,000 0.12 4.90 40.83 0.12 4.90 40.83 0.12 4.90 40.83 0.12 5.00 41.67 0.12 5.10 42.50
40,000 400,000 0.21 7.70 36.67 0.22 7.70 35.00 0.26 7.80 30.00 0.22 8.00 36.36 0.21 8.30 39.52
40,000 600,000 0.29 10.70 36.90 0.34 10.70 31.47 0.31 10.90 35.16 0.31 11.20 36.13 0.31 11.70 37.74
40,000 800,000 0.39 13.00 33.33 0.47 13.10 27.87 0.52 13.40 25.77 0.55 13.80 25.09 0.43 14.40 33.49

Geometric mean 40.32 38.41 37.26 41.22 41.10

for a given path P. Each one of the 180 instances was solved the first set consider one resource constraint while those in the
setting p ¼0.1, 0.2, 0.4, 0.6, and 0.8, thus leading to a total of 900 second set consider 10 resource constraints. For each problem in
runs of the pulse algorithm. Note that p A ½0,1; and low (high) Beasley and Christofides [1], Zhu and Wilhelm [35] build 100
values of p mean that the resource consumption constraint is random instances with different arc costs, for a total of 600
tight (loose). instances per set. Zhu and Wilhelm [35] extended the label-setting
Table 1 compares the results of the pulse algorithm against the algorithm (LSA) of Dumitrescu and Boland [13] to handle multiple
Lagrangian relaxation algorithm (LRA) proposed by Santos et al. resources and used it as benchmark algorithm for their own three-
[29] over different values for the tightness of the resource stage approach (TSA). Given that TSA is intended for solving
constraint (namely, parameter p). The instances are organized multiple subproblems under a column generation scheme, Zhu
by rows in three categories (i.e., small, medium, and large-sized and Wilhelm [35] define a threshold as the number of instances
instances) based on the number of nodes and arcs (columns 1 and (i.e., subproblems) solved with the benchmark algorithm required
2). Columns 3–7 show, for the different values of p, the average to match the total time spent by TSA, including its preprocessing
execution time (in seconds) for the pulse algorithm, the average procedures. If under a column generation scheme, it happens that
execution time as reported by Santos et al. [29], and the speedup the number of subproblems to be solved is larger than the
calculated as the ratio between the execution times. Finally, the threshold, it means that it is worth using TSA over the benchmark
last row reports the geometric mean of the speedups. Note that algorithm. Zhu and Wilhelm [35] define the threshold as follows:
the geometric mean, opposed to the arithmetic mean, avoids 8& 1 ’
being overly optimistic with large ratios obtained on few instance >
> t þ t2 3
< TSA 3 TSA if t  4t TSA
sizes (node-arc combination), thus it provides a fairer comparison tðTSA,Þ ¼ t  t TSA
>
>
of speedups [3]. For the sake of fairness we scaled our computa- :1 if t rt
3
 TSA
tion times using the LINPACK benchmark [12]. According to this
benchmark, our computer is approximately 2.1 times faster than where t 1TSA
and t 2TSA
are the times it takes TSA to complete stages
the one used by Santos et al. [29], so all the reported times were 1 and 2 of the preprocessing procedure at the beginning of the
3
scaled by this factor. column generation scheme; and t TSA and t  are the average times it
The results presented in Table 1 show that the pulse algorithm takes a single run of the stage 3 of TSA and the benchmark
consistently outperforms LRA in every instance, regardless of the algorithm (). Note that if the average time of the benchmark
constraint tightness. Speedups range from 16 to 60 while the algorithm () is less than the average time of stage 3 of TSA, it is
geometric mean shows that the pulse algorithm is consistently said that the benchmark algorithm outperforms TSA regardless of
roughly 40 times faster than LRA. In terms of absolute computing the number of subproblems to be solved, thus the threshold
time (not scaled), the pulse algorithm solved all the small net- tðTSA,Þ 91.
works in less than 0.1 s, the middle-sized networks in less than Table 2 compares the pulse algorithm with LSA and TSA on the
0.15 s, and the large networks in less than 0.3 s. Profiling our exact same set of instances by Zhu and Wilhelm [35]. We
algorithm we found that most of the computational time for these extended the pulse algorithm to handle multiple resources but
instances is spent at the initialization. to keep the computational budget under control, we obtained
lower bounds with the one-to-all shortest path for the cost and
4.2. Experiments on small and mid-sized instances (with multiple just one resource, while for the other resources the lower bounds
side constraints) were set to zero. Columns 1–4 show the instance, number of
resource constraints, number of nodes, and number of arcs,
Zhu and Wilhelm [35] present two sets of instances adapted respectively. Columns 5 and 6 present the average execution
from Beasley and Christofides [1] that range in size from 100 time to solve each subset of 100 instances for LSA and for the
nodes and 959 arcs to 500 nodes and 4978 arcs. The problems in pulse algorithm. Column 7 exhibits the speedup achieved with
382 L. Lozano, A.L. Medaglia / Computers & Operations Research 40 (2013) 378–384

Table 2
Computational results for the Zhu and Wilhelm [35] instances.

Instance Resources Nodes Arcs LSA (s) Pulse (s) Speedup TSA (s) sðTSA,PulseÞ sðTSA,LSAÞ

t 1TSA þ t 2TSA 3
t TSA

rcsp3 1 100 959 0.005 0.00067 7.5 0.140 0.00016 274.5 30

rcsp4 1 100 959 0.005 0.00069 7.3 0.125 0.00016 236.5 25
rcsp11 1 200 1971 0.004 0.00073 5.5 0.625 0.00093 1 222
rcsp12 1 200 1971 0.006 0.00072 8.3 0.516 0.00094 1 100
rcsp19 1 500 4978 0.006 0.00112 5.3 1.094 0.00297 1 411
rcsp20 1 500 4978 0.006 0.00119 5.1 0.891 0.00187 1 203

rcsp7 10 100 999 0.198 0.00443 44.7 0.906 0.00000 204.6 5

rcsp8 10 100 999 0.085 0.00418 20.4 0.188 0.00016 46.8 2
rcsp15 10 200 1960 0.513 0.00654 78.5 1.594 0.00000 243.9 3
rcsp16 10 200 1960 0.192 0.00459 41.8 0.313 0.00000 68.2 2
rcsp23 10 500 4868 6.224 0.00670 928.9 99.658 0.00109 17,763.1 16
rcsp24 10 500 4868 2.975 0.01361 218.6 7.859 0.00000 577.4 3

the pulse algorithm, that is, how many times faster is the pulse Table 3
over LSA. Columns 8 and 9 present the time used in preprocessing Comparing the execution time (s) for the single-thread and the multi-thread
versions of the pulse algorithm on instances with 40,000 nodes and 800,000 arcs.
and stage 3 for TSA. Column 10 shows the threshold that
compares the pulse algorithm with TSA, namely, the number of Instance Threads Speedup
runs (subproblems) that could be solved with the pulse, before it
starts to pay off to use TSA. Finally, column 11 shows the 1 2 4 8
threshold that compares LSA with TSA. Zhu and Wilhelm [35]
1 0.384 0.279 0.208 0.194 1.98
coded their algorithms in C/Cþþ and executed their experiments 2 0.237 0.204 0.195 0.149 1.59
on an Intel 3.2 GHz dual core CPU with 2 GB of RAM. Given that 3 0.221 0.185 0.161 0.143 1.55
their computer is better than ours, we assumed this handicap by 4 0.192 0.202 0.138 0.170 1.39
not scaling our times. 5 0.182 0.163 0.161 0.168 1.13
6 0.196 0.146 0.183 0.190 1.34
For the first set of instances (with one resource constraint),
7 0.170 0.171 0.152 0.182 1.12
Table 2 shows that the pulse algorithm consistently outperforms 8 0.208 0.224 0.148 0.140 1.48
LSA in all the experiments, with speedups ranging from 5.1 to 8.3. 9 0.228 0.208 0.142 0.137 1.66
The pulse algorithm dominated TSA in 4 out of 6 experiments 10 0.284 0.255 0.170 0.166 1.71
(denoted by the threshold of 1), meaning that regardless of the
number of subproblems solved, the total execution time of the
pulse algorithm will always be less than that of TSA. For those two processors Intel Xeon X5450 running at 3 GHz, 8 GB of RAM,
instances without proof of dominance by the pulse, the threshold and a 64-bit Windows Vista Ultimate operating system.
value gives the pulse a comfortable margin of at least 236 Table 3 shows that the parallel version of the pulse algorithm
subproblems. In terms of absolute computing time, all the is faster than the single-thread version for all the tested instances.
experiments were solved under 0.0012 s in average. For the The speedups range from 1.12 to 1.98. These promising speedups
second set of instances (with 10 resource constraints) the pulse might be the result of obtaining tighter primal bounds earlier by
algorithm significantly outperforms LSA reaching speedups of up exploring different regions of the network at the same time in
to 900 times in the largest instances. It also compares favorably different threads. However, further experimentation must be
with TSA, achieving threshold values of up to 17,700 and solving made in order to fully exploit and understand the advantages
all the experiments under 0.014 s. Over both sets of instances, the and tradeoffs while using the parallel version of the algorithm.
threshold value obtained with the pulse is several orders of
magnitude larger than the one obtained with LSA.
5. Embedding the pulse algorithm within column generation
4.3. Parallelizing the pulse algorithm schemes

Given that the recursive pulse function explores one node at a Our experience has shown that an effective and practical use
time until it reaches the end node, it is possible to call the of the pulse algorithm arises while solving auxiliary problems in
function over different nodes in parallel threads taking special column generation procedures. Note, however, that it is often the
care on how to update the global information. We parallelized the case that these auxiliary problems contain negative costs gener-
execution of the pulse algorithm by triggering at node vs different ated from the dual variables of the master problem. We success-
threads that begin the exploration starting from the outgoing fully adapted the pulse algorithm to solve the auxiliary problems
arcs. Two different threads cannot execute the same functions for a Multi-Activity Shift Scheduling Problem (MASSP) and a Bus
over the same node at the same time, if this happens, one thread Rapid Transit Route Design Problem (BRTRDP). The computational
must wait for the other to finish. In order to test our approach, we experiments in this section were executed on a Dell Precision
used the hardest instance (40,000 nodes and 800,000 arcs) from 7400 with 8 GB of RAM, two processors Intel Xeon X5450 running
Santos et al. [29] and derived 10 new instances by randomly at 3 GHz, on a 64-bit Windows Vista Ultimate operating system.
replacing the end node. Table 3 shows the execution time in
seconds required to solve the instances with 1, 2, 4, and 8 threads. 5.1. Multi-activity shift scheduling problem
The last column shows the speedup calculated as the ratio
between the best multi-thread approach and the single-thread The Shift Scheduling Problem (SSP) deals with the selection of
approach. For this experiment we used a machine with two a set of shifts to satisfy a daily demand for staff requirements. The
 L. Lozano, A.L. Medaglia / Computers & Operations Research 40 (2013) 378–384 383

Multi-Activity Shift Scheduling Problem (MASSP) is an extension Table 5

of the SSP where the staff may perform several work activities in Speeding up the solution of the BRTRDP with the
pulse algorithm.
the same shift and the assignment of activities to a shift is
restricted by work rules [8,7,10]. Depending on the problem at Instance Time (s)
hand, these work rules include activity length, number of breaks,
shift starting time, shift length, and maximum number of activ- IPa Pulse
ities per shift, among others. Typically, in the MASSP the schedule
BRT-S8 9 2
spans over a multi-day planning horizon. BRT-S10 39 7
Under a column generation scheme for the MASSP we modeled BRT-S15 2566 63
the auxiliary problem as a Shortest Path Problem with Resource BRT-S19 140,369 370
Constraints (SPPRC) with local time windows over each node and a
Xpress-MP optimizer version 19.00.00.
global resource constraints [28]. This variant of the SPPRC is
defined over an acyclic graph with negative costs, several pro-
blem-specific side constraints, and up to 2000 nodes and 50,000 to go backwards, meaning that buses stop in stations in a given
arcs. To solve the SPPRC, we modified the pulse algorithm as direction (between terminal stations). To address this auxiliary
follows: being the graph acyclic, we did not have to check for problem, we modified the pulse algorithm as follows: the dom-
cycles; the dominance pruning strategy cannot prune partial inance pruning strategy remains unchanged; the infeasibility
paths that have not met the minimum requirements for all pruning strategy discards any path that visits a forbidden station;
resources; the infeasibility pruning strategy was extended to and the bounds pruning strategy is no longer used because of the
prune using multiple resource constraints as in Section 4.2; and negative cost cycles. Table 5 compares two column generation
the bounds pruning strategy still applies. Table 4 compares the schemes for the BRTRDP: one that solves the auxiliary problem as
performance on three real-world MASSP instances solved using an integer program using Xpress-MP optimizer version 19.00.00;
traditional column generation, where the auxiliary problem is and a second one that uses the pulse algorithm. The first column
solved using integer programming, against a pulse-embedded shows the instance name where S is the prefix for the number of
column generation. The first column shows the instance name stations in the system. Columns 2 and 3 show the total computa-
where A, TIL and D are the prefixes for the number of activities, tional time needed for both approaches.
the time interval length in minutes, and the planning horizon in The computational speedup is up to 379 times due to the pulse
days. Columns 2 and 3 show the total computational time needed and it is noteworthy that larger instances show the most
using integer programming (IP) and the pulse algorithm for the significant improvement. Although the bound pruning strategy
auxiliary problems. All the integer programs were solved using was not used, the remaining strategies were strong enough to
Xpress-MP optimizer version 19.00.00. achieve these results. The faster computational times allowed us
Adapting the pulse algorithm to address this SPPRC (arising to solve instances that the approach with integer programming
from the MASSP) was a relatively easy task given that the search was simply unable to solve.
logic of the algorithm remains unchanged and the multiple work
rules for the shifts made the infeasibility pruning strategy very
strong. The faster computational times allowed us to solve harder 6. Conclusions and future work
and larger instances that were simply not tractable with one of
the best branch and bound implementations found in commercial We introduced an exact algorithm for the CSP that solves large
optimizers (i.e., Xpress-MP). networks with up to 40,000 nodes and 800,000 in a reasonable
amount of time (under 0.3 s). The algorithm, based on the idea of
5.2. Bus rapid transit route design problem propagating pulses through a network, performs better than the
algorithm by Santos et al. [29] over a large testbed from the
The Bus Rapid Transit Route Design Problem (BRTRDP) is the literature, achieving speedups of up to 60 times. It also outper-
problem of finding a set of routes and frequencies that minimizes forms the label-setting algorithm by Dumitrescu and Boland [13]
the operational and the passenger costs (travel time) while reaching speedups of up to 900 times on the proposed instances
simultaneously satisfying the system’s technical constraints, such and compares favorably with the up-to-date approach by Zhu and
as meeting the demands for trips, bus frequencies, and lane Wilhelm [35].
capacities. Additionally, this algorithm can be seen as a flexible frame-
To solve the BRTRDP we used a decomposition approach work able to handle difficult constraints modeled as pruning
where new routes are added dynamically to a master problem strategies. The idea of the algorithm is easy to understand and to
[14,17]. Similar to Section 5.1, an auxiliary problem searches a implement, and the pruning strategies can be used as modules
network for suitable routes to be added to the master problem. In and adapted to problem specific conditions. Also, the only para-
contrast to the SPPRC (in the MASSP) this network contains meters to be defined by the user are the maximum number of
negative costs and possibly negative cost cycles. Given the labels stored at each node (i.e., Q), and the label storage rule.
physical infrastructure, there is a constraint that forbids a route We show that the algorithm is easy to parallelize, and
illustrate this fact with a first parallel version of the algorithm
Table 4 that exhibits promising results when compared with the single-
Speeding up the solution of real-world MASSP instances with the pulse algorithm. thread version of the algorithm on large-scale networks.
From a practical perspective, the pulse algorithm has been
Instance Time (s)
successfully extended for solving auxiliary problems under a
IPa Pulse column generation scheme arising in two different contexts: a
bus rapid transit route design problem and a multi-activity shift
A5-TIL30-D7 3549 17 scheduling problem. A very similar auxiliary problem arises in the
A7-TIL30-D7 71,216 83
A16-TIL30-D7 848,402 629
context of the integrated vehicle- and crew-scheduling problem
[31], where the pulse algorithm is likely to accelerate the column
a
Xpress-MP optimizer version 19.00.00. generation process.
384 L. Lozano, A.L. Medaglia / Computers & Operations Research 40 (2013) 378–384

Work currently underway includes: extending the algorithm [12] Dongarra JJ. Performance of various computers using standard linear equa-
to tackle the Elementary Shortest Path Problem with Resource tions software. Technical report CS-89-85, University of Tennessee, USA;
2009.
Constraints (ESPPRC), a problem that arises while solving the [13] Dumitrescu I, Boland N. Improved preprocessing, labeling and scaling
Vehicle Routing Problem with Time Windows (VRPTW) under a algorithms for the weight-constrained shortest path problem. Networks
column-generation procedure; applying the algorithm to a real 2003;42(3):135–53.
[14] Feillet D, Gendreau M, Medaglia AL, Walteros JL. A note on branch-and-cut-
medical task scheduling problem; extending its functionality to and-price. Operations Research Letters 2010;38(5):346–53.
tackle a subproblem in the public transport line planning problem [15] Gamache M, Soumis F, Marquis M, Desrosiers J. A column generation
[4]; extending the pulse algorithm to the biobjective shortest path approach for large-scale aircrew rostering problems. Operations Research
1999;47(2):247–63.
problem; and exploring further the parallelization of the algo- [16] Garey MR, Johnson DS. Computers and intractability: a guide to the theory of
rithm to take full advantage of multiple core processors or NP-completeness. New York: W.H. Freeman; 1979.
graphics processing units. [17] González JE, Lozano L, Walteros JL, Feillet D, Medaglia AL. Solving the bus
rapid transit route design problem with general topologies via simultaneous
column and cut generation. Technical report COPA-2012-10, Universidad de
los Andes; 2012.
Acknowledgments [18] Graves GW, McBride RD, Gershkoff I, Anderson D, Mahidhara D. Flight crew
scheduling. Management Science 1993;39(6):657–82.
[19] Gromicho J, van Hoorn JJ, Kok AL, Schutten JMJ. Restricted dynamic program-
The authors would like to thank Professors Joa~ o Coutinho- ming: a flexible framework for solving realistic VRPs. Computers & Opera-
Rodrigues at Universidade de Coimbra (Portugal), Xiaoyan Zhu at tions Research 2011. http://dx.doi.org/10.1016/j.cor.2011.07.002.
[20] Grönkvist M. Accelerating column generation for aircraft scheduling using
the University of Tennessee (U.S.A), and Wilbert E. Wilhelm at
constraint propagation. Computers & Operations Research 2006;33:2918–34.
Texas A&M University (U.S.A) for their generosity in sharing with [21] Handler GY, Zang I. A dual algorithm for the constrained shortest path
us their testbeds for the CSP. Professors Louis-Martin Rousseau at problem. Networks 1980;10:293–310.
École Polytechnique de Montréal and Juan G. Villegas at Uni- [22] Joksch H. The shortest route problem with constraints. Journal of Mathema-
tical Analysis and Applications 1966;14(1):191–7.
versidad de Antioquia for their insightful comments. Last, but not [23] Lavoie S, Minoux M, Odier E. A new approach for crew pairing problems by
least, we want to thank Daniel Duque, our MS student at column generation with an application to air transportation. European
Universidad de los Andes, who developed the double-bucket Journal of Operational Research 1988;35:45–58.
[24] Lawler EL. Combinatorial optimization: networks and matroids. New York:
shortest path implementation and coded an animation of the Dover Publications; 2001.
pulse algorithm that facilitates its understanding. [25] Muter _I, Birbil S- _I, Bülbül K, S- ahin G, Yenigün H, Tas- D, et al. Solving a robust
airline crew pairing problem with column generation. Computers & Opera-
tions Research 2010. http://dx.doi.org/10.1016/j.cor.2010.11.005.
References [26] Giacomo Nannicini, Daniel Delling, Dominik Schultes, Leo Liberti. Bidirec-
tional An search on time-dependent road networks. Networks 2012;59(2):
240–51.
[1] Beasley JE, Christofides N. An algorithm for the resource constrained shortest [27] Reinhardt LB, Pisinger D. Multi-objective and multi-constrained non-additive
path problem. Networks 1989;19:379–94. shortest path problems. Computers & Operations Research 2011;38:605–16.
[2] Bertsekas DP, Tsitsiklis JN. An analysis of stochastic shortest path problems. [28] Restrepo MI, Lozano L, Medaglia AL. Constrained network-based column
Mathematics of Operations Research 1991;16(3):580–95. generation for the multi-activity shift scheduling problem. International
[3] Bixby RE. Solving real-world linear programs: a decade and more of progress. Journal of Production Economics. http://dx.doi.org/10.1016/j.ijpe.2012.06.
Operations Research 2002;50(1):3–15. 030; 2012.
[4] Borndörfer R, Grötschel M, Pfetsch ME. A column-generation approach to line [29] Santos L, Coutinho-Rodrigues J, Current JR. An improved solution algorithm
planning in public transport. Transportation Science 2007;41(1):123–32. for the constrained shortest path problem. Transportation Research Part B
[5] Carlyle WM, Royset JO, Wood RK. Lagrangian relaxation and enumeration for 2007;41:756–71.
solving constrained shortest-path problems. Networks 2008;52(4):256–70. [30] Smith OJ, Boland N, Waterer H. Solving shortest path problems with a weight
[6] Cherkassky BV, Goldbberg AV, Radzik T. Shortest path algorithms: theory and constraint and replenishment arcs. Computers & Operations Research
experimental evaluation. Mathematical Programming 1996;73:129–74. 2012;39(5):964–84.
[7] Côté M-C, Gendron B, Quimper C-G, Rousseau L-M. Formal languages for [31] Steinzen I, Gintner V, Suhl L, Kliewer N. A time-space network approach for
integer programming modeling of shift scheduling problems. Constraints the integrated vehicle- and crew-scheduling problem with multiple depots.
2009;16(1):54–76. Transportation Science 2010;44(3):367–82.
[8] Côté M-C, Gendron B, Rousseau L-M. Grammar-based integer programming [32] Stojković M, Soumis F, Desrosiers J. The operational airline crew scheduling
models for multi-activity shift scheduling. Management Science 2011;57(1): problem. Transportation Science 1998;32(3):232–45.
1–13. [33] Villeneuve D, Desaulniers G. The shortest path problem with forbidden paths.
[9] Deb K, Pratap A, Agarwal S, Meyarivan T. A fast and elitist multi-objective European Journal of Operational Research 2005;165:97–107.
genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation [34] Zhu X, Wilhelm WE. Three-stage approaches for optimizing some variations
2002;6(2):182–97. of the resource constrained shortest-path sub-problem in a column genera-
[10] Demassey S, Pesant G, Rousseau L-M. A cost-regular based hybrid column tion context. European Journal of Operational Research 2007;183:564–77.
generation approach. Constraints 2006;11(4):315–33. [35] Zhu X, Wilhelm WE. A three-stage approach for the resource-constrained
[11] Derigs U, Friederichs S, Schäfer S. A new approach for air cargo network shortest path as a sub-problem in column generation. Computers & Opera-
planning. Transportation Science 2009;43(3):370–80. tions Research 2012;39:164–78.