ALGORITHMS FOR NON-LINEAR AND STOCHASTIC RESOURCE
CONSTRAINED SHORTEST PATH.

arXiv:1504.07880v2 [] 15 Feb 2017

AXEL PARMENTIER
Abstract. Resource constrained shortest path problems are usually solved thanks to a smart
enumeration of all the non-dominated paths. Recent improvements of these enumeration algorithms
rely on the use of bounds on path resources to discard partial solutions. The quality of the bounds
determines the performance of the algorithm. The main contribution of this paper is to introduce
a standard procedure to generate bounds on paths resources in a general setting which covers most
resource constrained shortest path problems, among which stochastic versions.
In that purpose, we introduce a generalization of the resource constrained shortest path problem where the resources are taken in a monoid. The resource of a path is the monoid sum of
the resources of its arcs. The problem consists in finding a path whose resource minimizes a nondecreasing cost function of the path resource among the paths that respect a given constraint.
Enumeration algorithms are generalized to this framework. We use lattice theory to provide polynomial procedures to find good quality bounds. These procedures solve a generalization of the
algebraic path problem, where arc resources belong to a lattice ordered monoid. The practical
efficiency of the approach is proved through an extensive numerical study on some deterministic
and stochastic resource constrained shortest path problems.

1. Introduction
1.1. Problem statement. Several methods have recently been developed to increase the efficiency
of the shortest path problem solvers [8]. A common feature to all these methods is the use of lower
bounds on costs of paths to discard partial paths in an enumeration of all the paths. In this paper,
we exploit the properties of lattices to extend this lower bound idea to algorithms for a generic
version of the resource constrained shortest path problem where arc resources belong to a lattice
ordered monoid. The lattice ordered monoid point of view covers a wide range of applications,
among which stochastic versions of the resource constrained shortest path problem.
Let (R, 6) and (S, ≤) be two partially ordered sets. A map ρ : R → S is isotone if x 6 y implies
ρ(x) ≤ ρ(y). Let (R, ⊕) be a set endowed with a law of composition. (R, ⊕) is a monoid if ⊕ is
associative and admits a neutral element in R. A partial order 6 is a compatible order on (R, ⊕)
if all translations y 7→ x ⊕ y and y 7→ y ⊕ x are isotone. An ordered monoid (R, ⊕, 6) is a monoid
endowed with a compatible order. A partially ordered set (R, 6) is a lattice if any pair of elements
(x, x̃) of R2 admits a greatest lower bound x ∧ x̃ or meet and a least upper bound x ∨ x̃. A lattice
ordered monoid (R, ⊕, 6) is an ordered monoid such that 6 induces a lattice structure.
Let (R, ⊕, 6) be a lattice ordered monoid. Consider the following problem.
Monoid Resource Constrained Shortest Path Problem
Input. A digraph D = (V, A), two vertices o, d ∈ V , a collection (xa ) ∈ RA , and two isotone
mappings c : R → R and ρ : R → {0, 1}.


L
L
Output. An o-d path P such that ρ
a∈P xa = 0 and with minimum c
a∈P xa .
L
R is the set of resources. The resource of a path P is a∈Pxa and we denote it xP , the cost of
L
L
P is c
a∈P xa , and P is feasible if and only if ρ
a∈P xa is a equal to 0. The function ρ is
later referred as the infeasibility function.
1

The Monoid Resource Constrained Shortest Path Problem is N P-hard as it contains
the usual Resource Constrained Shortest Path Problem, which is obtained by using the
set R = R2 , the cost c(x1 , x2 ) = x1 , and the function ρ((x1 , x2 )) equal to 1 if and only if x2 > M
for a given M ∈ R. Finally, we emphasize that ⊕ is possibly non commutative.
This lattice ordered monoid framework has two main strengths. First, it enables to deal with
stochasticity in path problems. And second, the lattice ordered monoid structure enables to define
polynomial procedure to compute lower bounds on paths resources, and to use these bounds to
speed-up solution algorithms. We now sketch the main ideas behind the treatment of stochastic
path problems and the solution schemes.
1.2. Application to stochastic path problems. Suppose that for each arc a we have a random
variable ξa . A large class of stochastic path problems can be expressed as follows. Given an origin
vertex o and a destination d, find an o-d path minimizing
!#
"
X
,
ξa
min E f
P

a∈P

under the constraint

P

X

ξ>τ

a∈P

!

≤ α,

where f is a non-decreasing function, and τ and α are constants. We use two main ideas to
model such problems within the Monoid Resource Constrained Shortest Path Problem
framework. First, a space of random variables endowed with the addition and the almost sure
order is a lattice ordered monoid. And second, many if not most probability functionals that
intervene in stochastic path problems are monotone with respect to this order and can therefore
be modeled withing P
this framework. On our example, it suffices to define xa = ξa , c(ξ) = E [f (ξ)],
and ρ(ξ) = 1(α,1] (P( a∈P ξ > τ )), where we introduce the notation 1I that we frequently use in
the paper to denote the indicator function of a set I.
Finally, the use of alternative stochastic orders enables to exploit assumptions such as the independence of the ξa to improve the performance of our solution algorithms.
1.3. Solution scheme. We propose
L solution schemes for the Monoid Resource Constrained
Shortest Path Problem when a∈C xa > 0 for each cycle C in D. These schemes are in two
steps. We start by computing for each vertex v a lower bound bv on the resource xP of all the v-d
paths. Then, we use these lower bounds to discard partial solutions in an enumeration of all the
paths.
Bounds bv are used in enumeration algorithms to compute lower bounds on the resource of any
o-d path starting by an o-v path P . Indeed, given an o-v path P and a v-d path P̃ ,
xP ⊕ bv 6 xP ⊕ xP̃ = xP +P̃ ,
where P + P̃ denotes the path composed of P followed by P̃ . The resource xP ⊕ bv is therefore a
lower bound on the resource of any o-d path starting by P . As ρ and c are isotone, ρ(xP ⊕ bv ) = 1
implies that there is no feasible o-d path starting by P , and c(xP ⊕ bv ) is a lower bound on the cost
of any o-d path P starting by P . The enumeration algorithm we propose therefore enumerates all
the paths satisfying
(1)

ρ(xP ⊕ bv ) = 0 and

B
c(xP ⊕ bv ) ≤ cU
od ,

B
where v is the destination of P and cU
od is an upper bound on the cost of an optimal solution.
The practical efficiency of the approach relies on our ability to compute in a preprocessing
V a tight
lower bound on the resource of all the v-d paths. As (R, 6) is a lattice, the meet bopt
=
v
P ∈Pvd xP ,
2

where Pvd denotes the set of v-d paths, is the tightest lower bound. Indeed, by definition of the
meet,
^
bv 6 xP for all P in Pvd ⇔ b 6 bopt
=
xP .
v
P ∈Pvd

Unfortunately, we have shown [71] that, unless P = N P, there is no polynomial algorithm that
enables to compute bopt
in polynomial time even on some simple lattice ordered monoid where
v
resources are positive. However, we provide polynomial procedures that compute lower bounds
on bopt
v . To that purpose, we show that the following equation always admits solutions, that its
solutions are lower bounds on bopt
v , and we introduce polynomial procedures that compute its
greatest solution.

 bd = 0,
^
(2)
bv = bv ∧
(x(v,u) ⊕ bu ) for all v ∈ V \{d}.

(v,u)∈δ+ (v)

Note that Equation (2) can be interpreted as a generalization of the Ford-Bellman dynamic programming equation for the usual shortest path problem, where the minimum has been replaced by
the meet operator ∧, and the sum by the operator ⊕. Our polynomial procedures are generalizations
of the Ford-Bellman and of the Dijkstra algorithm for the usual shortest path problem.

1.4. Contributions and plan. We can sum-up our contributions as follows:
• We introduce a versatile algebraic framework for constrained shortest path problems. This
framework notably enables to deal with non-linearity and stochasticity in the objective and
in the constraints.
• We provide polynomial procedures to compute the solution of the generalized dynamic
programming equation (1). The problem of solving this equation when the resource set
has the structure of idempotent semiring is known as the algebraic path problem and has
received much attention. Our procedures extend well-known algorithms of this community
to the more general setting of lattice ordered monoids. Besides, contrary to the algebraic
path problem community, we do not interpret the solution of (1) as the solution of the
problem but as lower bounds that can then be used in an enumeration algorithm. This new
interpretation together with the extension to lattice ordered monoids enable to apply these
algebraic methods to a much wider range of paths problems. The lattice ordered monoid
point of view is notably essential to model stochastic path problems.
• We generalize the usual enumeration algorithms for resource constrained shortest path
problems to our framework. The use of bounds is these algorithms is easy thanks to the
lattice ordered monoid structure.
• Concerning stochastic path problems, we show that our framework can deal with most
probability functionals of interest, among which the version independent risk measures.
Besides, we can deal with a wide range of probability distributions for the random variables
ξa , and we can solve approximated versions of problems with other distribution through
sampling.
• We show the practical efficiency of our approach through extensive numerical experiments
on deterministic and stochastic resource constrained shortest path problems. To the best
of our knowledge, our algorithms are the first practically efficient ones for paths problems
with probabilistic constraints.
• Finally, we provide strategies to improve the performance of our enumeration algorithms
on difficult problems thanks to a longer preprocessing.
After introducing some notions on digraphs and ordered algebraic structures in Section 2, we
detail the connections between our framework and algorithms and those of the literature on usual,
3

algebraic, resource constrained, and stochastic path problems in Section 3. Sections 4 and 5 introduce respectively our enumeration and bounding algorithms for the Monoid Resource Constrained Shortest Path Problem. Section 6 tests the numerical performance of our algorithms
on instances of the usual resource constrained shortest path problem. In Section 7, we explain how
to model stochastic path problems within our framework, and test numerically the performance of
our algorithm on some non-constrained and constrained stochastic path problems. We conclude
the paper with techniques to improve the performance of our algorithms on difficult instances in
Section 8.
2. Preliminaries
2.1. Digraphs. A digraph D is a pair (V, A), where V is the set of vertices and A is the set of arcs
of D. An arc a links a tail vertex to a head vertex. An arc a is incoming to (resp. outgoing from) v
if v is the head (resp. the tail) of v. The set of arcs incoming to (resp. outgoing from) v is denoted
δ− (v) (resp. δ+ (v)). A path is a sequence of arcs a1 , . . . , ak such that for each i ∈ {1, . . . , k − 1},
the head vertex of ai is the tail vertex of ai+1 . Note that with this definition, paths can contain
multiple copies of an arc or of a vertex. A path P is said to be elementary if it contains at most
one copy of each vertex. The origin of a path it the tail of its first arc and its destination is the
head of its last arc. Given two vertices o and d in V , an o-d path P is a path with origin o and
destination d. Finally, a cycle is a path whose origin is identical to its destination.
Given two paths P and Q such that P ends in the origin of Q, we denote P + Q the path made
of P followed by Q. Given two vertices u and v, we denote Puv the set of u-v paths.
2.2. Lattice ordered monoids and algebraic structures. We now introduce additional properties on lattices. When a subset S of a partially ordered set admits a greatest
Wa
V lower
V bound (resp.
least
W upper bound), we again call it its meet (resp. its join), and denote it S or x∈S x (resp. S
or x∈S x). Any finite subset of a lattice admits a meet and a join. A lattice is complete if any
subset S ⊆ R admits a meet and a join. It is conditionally complete if any bounded subset S ⊆ R
admits a meet and a join.
(3)

All the lattices we consider in this paper are conditionally complete.

If R is conditionally complete, then R ∪ {−∞, +∞} is complete, where −∞ (resp. +∞) is
smaller (resp. greater) than any element in R. The lattice R ∪ {−∞, +∞} is a completion of R.
We sometimes need our lattices to be complete to be able to define some quantities needed in the
paper. When such quantities are defined, we mention that they may belong to the completion
R ∪ {−∞, +∞} of R.
We denote 0 the neutral element of the operator ⊕ of a lattice ordered monoid (R, ⊕, 6). A
resource x is positive if x > 0. A lattice ordered monoid (R, ⊕, 6) is a lattice ordered group if
(R, ⊕) is a group. A lattice ordered group is Archimedean if, for each x, x̃ ∈ R such that x > 0,
there exists n in Z+ such that nx > x̃, where nx = x
· · ⊕ x}.
| ⊕ ·{z
n times

3. Literature review

As an algebraic framework for non-linear and stochastic resource constrained shortest path problems, our work is at the cross-road of several branches of the literature. First, our algorithms are
naturally interpreted as generalizations of the usual shortest path problem algorithms. Second,
our bounding algorithms are generalizations to lattice ordered monoids of those of the algebraic
path problem community. Third, our framework can be seen as a versatile alternative to other resource constrained shortest path framework with enhanced version of the enumeration algorithms.
4

Finally, the restriction of our algorithms to stochastic path problems compare favorably to the
existing literature on the topic.
3.1. Usual shortest path problem. Variants of the Shortest Path Problem have been thoroughly studied during the last six decades. As we already mentioned, there are two main types
of algorithms for the Shortest Path Problem. The first ones, such as Dijkstra’s algorithm or
Ford-Bellman algorithm, compute the shortest path between one vertex and all the other ones. In
that sense their output is a solution of the dynamic programming equation (2) where ⊕ is the usual
sum on R and the meet ∧ is the minimum. The standard algorithm to solve it is Dijkstra’s algorithm [27] when arcs costs are non-negative, and Ford-Bellman [10, 36] when they can be negative.
We generalize both of them to compute solutions of (2).
The second ones are enumeration algorithms and use bounds to discard paths in an enumeration
of all the paths. The typical example of enumeration algorithm is A∗ algorithm [45], that we
generalize to our setting. These algorithms are called goal directed algorithms as they compute
the shortest path only between a given pair of origin and destination vertices. When good bounds
are used, enumeration algorithms are faster than polynomial algorithms. Bast et al. [8] survey the
rich literature developed in the last few years on the topic in the context of online route planning
systems. However, the goals of these recent contributions are orthogonal to our ones: their objective
is to be able to compute quickly the solution of an easy path problem between any o-d pair, while
we want to compute the solution of a difficult problem between one given o-d pair.
3.2. Bounding algorithms and algebraic path problems. Equation (2) is not the first generalization of the usual dynamic programming equation. Indeed, when (R, ∧, ⊕) is an idempotent
semiring, the problem of solving Equation (2) is called the algebraic path problem. The literature
on the topic considers idempotent semirings of various generality [3, 7, 17, 34, 40, 59, 65, 76, 84]
and is surveyed in [34]. Our framework generalizes the algebraic path problem: the idempotent
semiring structure is stronger than the lattice ordered monoid one. Indeed, if (R, +, ×) is an idempotent semiring, then (R, ×, 6+ ) is a lattice ordered monoid, where 6+ is the idempotent semiring
canonical order: x 6+ x̃ if and only if x + x̃ = x. Its meet operator is +. In the other direction,
if (R, ⊕, 6) is a lattice ordered monoid with meet operator ∧, then (R, ∧, ⊕) is an idempotent
semiring if and only if ⊕ distributes with respect to ∧, i.e.
a ⊕ (b ∧ c) = (a ⊕ b) ∧ (a ⊕ c) and

(a ∧ b) ⊕ c = (a ⊕ b) ∧ (a ⊕ c).

If b†v is not necessarily equal to bopt
on lattice ordered monoids, these two quantities are equal on
v
idempotent semirings.
Two types of algorithms have been developed for the algebraic path problem. The first ones generalize respectively the Ford-Bellman [22, 23] and the Dijkstra [65] algorithm. Our algorithms can
be seen as generalization to the lattice ordered monoid setting of these algorithms. Our generalized
Dijkstra algorithm is in particular very similar to the one developed by Mohri [65]. Therefore,
we can use results from the algebraic path problem community to obtain stronger versions of our
convergence theorems when ⊕ distributes with respect to ∧. The second type of algorithms for the
algebraic path problem consider Equation (2) as a system of linear equations in the idempotent
semiring (R, ⊕, 6), and generalize the Gauss-Seidel and the Gauss-Jordan algorithms to that setting [40, 84]. Unfortunately, these algorithms do not generalize well to the lattice ordered monoid
setting.
Finally, we present in Section 8 a technique to improve the quality of bounds that builds a lower
envelope on the set {x ∈ R|xP 6 x for some P in Pvd }. Techniques to build such a lower envelope
have recently been proposed [4, 5] in the context of static program analysis. However, their process
for building the lower envelope is orthogonal to our one.
5

3.3. Enumeration algorithms and resource constrained shortest path problems. Irnich
& Desaulniers [48] provide a resource constrained shortest path framework and survey the exact and
heuristic approaches to resource constrained shortest path problems. Their resource constrained
shortest path framework is based on the notion of resource extension functions [47]: there is one such
function ζa for each arc a, and the resource of a path formed of the arcs a1 , . . . , ak is ζak ◦. . . ◦ζa1 (0).
The main difference between their framework and our one is the associativity of ⊕. This makes our
framework slightly less versatile, but enables us to compute lower bounds on paths resources and
to use them to speed up the resolution. Besides, most problems of the literature can be modeled
within our framework, as we argue in Chapter 3 of [71].
There are three main types of exact solution schemes to solve resource constrained shortest path
problems: constraint programming, branch and bound, and enumeration algorithms. Constraint
programming approaches [25, 31, 42, 55, 75] combine specifically designed search, domain reduction, and propagation algorithms. Concerning Branch-and-Bound algorithms, specific branching
patterns have been developed by branch and bound algorithms for resource constrained shortest
path problems: they branch on cycles, on arcs, and on resources [9, 15, 47, 73]. Finally, our work
enters in the field of enumeration algorithms, and we now detail the literature on that topic.
In their survey on resource constrained shortest path problems, Irnich & Desaulniers [48] describe
the enumeration algorithms of the literature as variants of a generic enumeration algorithm. This
generic enumeration algorithm has many similarities with the one we propose in Section 4. To
obtain a practical algorithm from the generic enumeration algorithm, one must choose define a
key, some bounds, and a dominance rule. The key defines which paths are processed first. The
bounds and the dominance rules are what enable to discard paths. The dominance rule is an order
on the set of resources that enables to discard paths whose resources are dominated. Desrochers
& Soumis [26] and [72] provide specific keys for routing problems with time windows, but these
strategies apply to many resource constrained shortest path problems. Irnich [47] provides general
techniques to define resource extension functions leading to good dominance rules, and techniques
to handle path with identical resources [49]. The remaining of the techniques are problem-specific
[9, 33, 46, 49, 56, 58]. Finally, we note that variants of the enumeration algorithm based on the
k-shortest path problem [30] have been proposed [9, 44, 80]. These variants and problem-specific
dominance rules can be used within our setting when appropriated.
Several techniques have been proposed to compute bounds when (R, ⊕, 6) is Rn endowed with its
product order and sum. Some contributions solve a usual shortest path problem for each component
of the resource [26, 28, 53, 62]. Another branch of the literature uses Lagrangian relaxation on an
integer formulation of the problem to obtain lower bounds [16, 28, 44, 80]. These methods require
the absence of negative cost cycles, in order to be able to solve the Langrangian relaxation using
a shortest path problem. The case with negative cost cycles is considered by Feillet et al. [33].
When (R, ⊕, 6) is Rn , these Lagrangian techniques provide bounds that are typically tighter than
our ones, but that require longer computations along the enumeration algorithm.
The strength of our framework lies in our bounding algorithms, that enable to use lower bounds
to discard path and good keys in non-linear and stochastic settings.
3.4. Stochastic path problems. There are two types of stochastic path problems: offline problems, where the entire path is chosen a priori, and online problems, whose solution is a policy that
updates the path used given the realization of uncertainty on the first arcs. Given that the solution
of the Monoid Resource Constrained Shortest Path Problem in an o-d path, it enables
to model offline stochastic path problems. However, interestingly, the lower bounds computed in
our framework provide an optimal policy [71] for the well-studied online stochastic on time arrival
problem [32, 35, 38, 39, 43, 68, 78, 79]. We now review the approaches to the different offline
stochastic path problems.
6

Stochastic shortest path problems without constraints have been extensively studied since the
seminal work of Frank [37]. Models differ by the type of distribution they use for arc random
variables, and by the probability functional they optimize. The objective of a first line of paper is
to find a path maximizing the probability of on time arrival, or analogously, a path with minimum
quantile of given order. Approaches have been developed for both continuous [19, 37, 69, 70]
and discrete [64] distributions. Chen et al. [20] describe an efficient labeling algorithm to deal
with normal distributions on the arcs. This algorithm is not so far from our label correcting
algorithm applied with the lattice ordered monoid presented in Section 5.1.2 of [71] when restricted
to ρ(·) = P(· > τ ). A second line of papers defines a shortest path as a path minimizing the
expectation of a cost function [61]. Dynamic programming can be used when cost functions are
affine or exponential [29]. Murthy & Sarkar [66, 67] present an efficient labeling algorithm when arc
distributions are normal and cost functions are piecewise-linear and concave. All these objective
functions are isotone with respect to the stochastic orders we use and can therefore be modeled
within the Monoid Resource Constrained Shortest Path Problem framework. Besides,
we show in Section 7 that we can also model version independent risk measures. Concerning the
distributions, our approach can handle all the distributions used in the literature, among which
discrete, normal, and gamma distributions for independent random variables, and scenario based
distributions for non-independent random variables. In this paper, we focus on discrete independent
distributions and scenario based distributions. Lattice ordered monoids for other distributions can
be found in Chapter 5 of [71]. The specificity of our solution approach is the use of lower bounds
for stochastic orders.
Our approach cannot deal with positive linear combinations of means and variances in the objective [69, 70, 81], as these objective functions are not isotone with respect to stochastic order.
The problem of finding a minimum cost path for deterministic arc costs under stochastic resource
constraints have been introduced in [57], and a solution algorithm based on linear programming is
derived. We are not aware of other approaches specifically dedicated to stochastic constraints in
path problems. However, such problems have been considered in the context of column generation.
A wide range of stochastic versions of the traveling salesman and the vehicle routing problems have
been studied in the last decades. When such problems are solved by column generation, the pricing
subproblem is a stochastic resource constrained shortest path problem. Uncertainty in customer
presence [50, 51], in demand [11, 12, 41, 82], and in travel time [2, 18, 52, 54, 60, 63, 77, 83] have
been considered. Using the modeling techniques we introduce in Chapter 3 of [71], most probability
functionals considered in this literature can be dealt with using the lattice ordered monoid and the
probability functionals presented in this chapter. However, we underline that, in the context of
vehicle routing problems, the graph is often complete. In that case, the bounds provided by our
approach are likely to be of poor quality, and thus our solution approach may not suit to the specific
structure of these problems on complete graphs.

4. Enumeration algorithms
4.1. Generic enumeration algorithms. In this section, we give three algorithms for the Monoid
Resource Constrained Shortest Path Problem: the generalized A∗ , the label dominance,
and the label correcting algorithms. These three algorithms share the same structure. They
enumerate all the paths in the graph using tests to discard partial paths. They differ only by the
tests they use to discard paths, and by the keys they use to determine in which order the paths
are processed. We therefore give a generic algorithm, and define the algorithms used in practice as
specializations of this generic algorithm. Later in this section, we sometimes call optimal path an
optimal solution of the Monoid Resource Constrained Shortest Path Problem.
7

Algorithm

Test

Key

Pre-processing Maintained
structures
structures
∗
∗
B
Generalized A (A )
(Low)
c(xP ⊕ bv )
bv
L, cU
od
U
B
Label dominance (dom.)
(Dom)
c(xP )
—
L, cod , Mv
B
Label correcting (cor.)
(Dom), (Low) c(xP ⊕ bv )
bv
L, cU
od , Mv
Table 1. Monoid Resource Constrained Shortest Path Problem algorithms.

We now describe the generic enumeration algorithm. A list L of partial paths P , and an upper
B
bound cU
od on the cost of an optimal solution are maintained. Initially, L contains the empty path
B
at the origin o, and cU
od = +∞. While L is not empty, the following operations are repeated.
(1) Extract a path P of minimum “key” from L. Let v be the destination of P .
(2) If v = d and P is feasible and better than the current solution, i.e. ρ(xP ) = 0, and c(xP ) <
UB
B
cU
od , then update cod to c(xP ).
(3) Else if “test” returns “yes”, extend P : for each arc a outgoing from v, add P + a to L.
We obtain our different algorithms by specifying the key and the test respectively in the first
and in the last step. We now introduce several keys and test.
Our first key and test rely on lower bounds bv on the resources of all the v-d paths. Section 5
provides polynomial algorithms to compute such bounds in a preprocessing. Given an o-v path P ,
the lower bound test is expressed as follows.
B
(Low) Does P satisfy ρ(xP ⊕ bv ) = 0 and c(xP ⊕ bv ) 6 cU
od ?
The isotony of ρ and c implies that a subpath P of an optimal path satisfies this test.
An o-v path P dominates an o-v path P̃ if xP 6 xP̃ . The dominance test maintains along the
algorithm a list Mv of non-dominated o-v paths for each vertex v, and is expressed as follows.
(Dom) Is P non-dominated by any path in Mv ?
If the answer is yes, then before extending P , we remove from Mv and L all the paths in Mv
dominated by P , and add P to Mv . The rationale behind this test is that, given the way Mv is
built, there is an optimal path whose subpaths are all non-dominated.
The aim of keys is to extend first the most promising paths. As an optimal solution is a feasible
solution of minimum cost, we use as keys estimations of the cost of an optimal path. When bounds
bv are computed, given an o-v path P , the quantity c(xP ⊕ bv ) provides a lower bound on the cost
of any o-d path starting by P , and is therefore a good estimator of how promising P is. When no
bounds are computed, we use c(xP ), the cost of the partial solution considered.
We can now describe our enumeration algorithms. The generalized A∗ is obtained from the
generic algorithm by using the lower bound c(xP ⊕ bv ) as key in Step 1, and extending a path P
in Step 3 only if it satisfies the lower bound test (Low). When (R, ⊕, 6) = (R, +, ≤), ρ = 0, and
c(x) = x, it corresponds to the usual A∗ algorithm. The label dominance algorithm is obtained
from the generic algorithm by using the partial path cost ρ(cP ) as key in Step 1, and extending an
o-v path P in Step 3 only if it satisfies the dominance test (Dom). The label correcting algorithm
is obtained from the generic algorithm by using the lower bounds c(xP ⊕ bv ) as key in Step 1, and
extending an o-v path P in Step 3 only if it satisfies both the lower bound test (Low) and the
dominance test (Dom).
The properties of these algorithms are summed up in Table 1. We underline fact that bounds bv
must be computed in a preprocessing when the generalized A∗ and the label correcting algorithms
are used, and that the label correcting and the label dominance algorithm need to maintain lists
of non-dominated paths Mv .
Alternative combinations of our tests and keys are possible but less interesting. Indeed, our
three algorithms reach trade-offs between the pre-processing time, and the quality of the keys and
8

tests. Once bounds bv have been computed, performing the lower bound test (Low), and computing
c(xP ⊕ bv ) require a fix number of operations ⊕, 6, and ∧ of the lattice ordered monoid. Thus, if
time has been spent computing bounds bv in a preprocessing, it is always better to use both the
lower bound test and the key c(xP ⊕ bv ). The alternative is to use the label dominance algorithm
to avoid spending time in the preprocessing.
The label dominance algorithm corresponds to the standard algorithm for the resource constrained shortest path problem [48] in the resource extension framework. The strength of the
Monoid Resource Constrained Shortest Path Problem framework is that it enables to
introduce the generalized A∗ and the label correcting algorithms, which both rely on the use of
bounds.
When it comes to practical performances, it is well known that when lower bounds bv can be
computed, the label correcting algorithm outperforms the label dominance algorithm on usual
resource constrained shortest path problems [28]. The relative performance of the generalized A∗
and of the label correcting algorithm depends on the problem considered. Indeed, the dominance
test enables to discard more paths, but its complexity is linear in the size of Mv , and it may slow
the algorithm if dominance is rare and lists Mv become large. This is notably the case of some
problems with numerous or stochastic constraints. Finally, the numerical experiments in Sections
6 and 7 show that the label dominance algorithm tends to perform well on easy problems, while
the generalized A∗ algorithm tends to perform better on difficult problems with many deterministic
constraints or one stochastic constraint.
4.2. Convergence of the algorithms. We now prove the convergence of our enumeration algorithms when xC > 0 for all the cycles C of D. When this assumption is not satisfied, these
algorithms must be adapted to discard non-elementary paths, using for instance the techniques
developed by Feillet et al. [33].
4.2.1. Generalized A∗ algorithm. The following assumptions will be used to define some settings
under which ones the generalized A∗ algorithm converges:
W
(4) For all a, b < R and x > 0 in R, there exists and n ∈ Z+ such that nq ⊕ a b.

−1
−1
There exists a feasible
_ o-d path P such that c ((−∞, c(xP )]) ∩ ρ (0) is upper-bounded by
a resource xM <
R.
W
In Assumptions (4) and (5), R may be in the completion of R. Assumption (4) is a weaker
version for ordered monoid of the Archimedean property. It is satisfied by all the lattice ordered
monoids considered in this paper. R2 endowed with its product sum and order is an example
of lattice ordered group that is not Archimedean but which satisfies Assumption (4). Finally,
Assumption (5) says that the set of resources of potentially optimal solutions is bounded. Indeed,
an optimal path must be feasible and of cost non-greater than c(xP ): the feasible resources belong
to ρ−1 (0), and c−1 ((−∞, c(xP )]) is the set of resources whose cost is non-greater than the cost of
path P .

(5)

Theorem 1. Suppose that at least one of the following conditions is satisfied.
(a) D is acyclic.
(b) Assumptions (4) and (5) are satisfied, xa is positive for each arc a, and bv > 0 for each
vertex v.
L
(c) Assumptions (4) and (5) are satisfied, ⊕ is commutative, and a∈C xa is positive for any
cycle C in D.
Then the generalized A∗ algorithm converges after a finite number of iterations, and at the end, if
B
cU
od is finite, then it is the cost of an optimal solution of the Monoid Resource Constrained
Shortest Path Problem. Otherwise, the problem admits no feasible solutions.
9

We note that in case (b), the hypothesis that bv > 0 is not restrictive, as xa > 0 implies that
xP > 0 for all paths P . Finally, Case (c) is notably satisfied when (R, ⊕, 6) is an Archimedean
lattice ordered group and cycle resources are positive. Indeed, Theorem 10.19 in [13] ensures that
any Archimedean lattice ordered group is commutative, and Assumption (4) is a consequence of
the Archimedean property in ordered groups.
Lemma 2. Let P be an o-d path satisfying ρ(xP ) = 0. Then at a given step of the generalized A∗
algorithm, at least one of the following statements is satisfied:
• there is a subpath P ′ of P in L,
B
• cU
od ≤ c(xP ).
Note that P ′ can be equal to P .
Remark 1. As c(xP ′ ⊕ bv ) ≤ c(xP ), where v is the destination of P , Lemma 2 implies that, if we
stop the algorithm before its convergence, the minimum c(xP ′ ⊕ bv ) for P ′ ∈ L provides a lower
bound on the cost of an optimal path.
Proof. We start with preliminary results. Paths are added to L only due to extension of paths in
Step 3. A path Q can therefore be in L only if its subpaths have been considered, removed from
L, and extended by the algorithm. Thus, at a given step of the generalized A∗ algorithm, for each
path Q with origin o, exactly one of the following statements is satisfied:
• Q has been considered by the generalized A∗ algorithm,
• a subpath Q′ of Q is in L,
• a strict subpath Q′ of Q has not been extended by the algorithm when considered.
Besides, if a feasible o-d path Q has already been considered, Step 2 of the algorithm implies
B
cU
od ≤ c(xQ ).
We now prove Lemma 2. Suppose that none of the statements of Lemma 2 are satisfied. As P is
a feasible o-d path, the two results above imply that a subpath P ′ of P has not been extended by
B
the algorithm when considered. Let P ′ be this subpath, and v ′ be its destination. As cU
od decreases
B
′
along the algorithm, the hypothesis implies that cU
od > c(xP ) when P is considered. As bv′ is a
′
lower bound on the resource of all v -d paths, we have xP ′ ⊕ bv′ 6 xP . By monotonicity of ρ and
B
feasibility of P , we have ρ(xP ′ ⊕ bv′ ) = 0. By monotonicity of c, we have c(xP ′ ⊕ bv′ ) ≤ c(xP ) < cU
od
when P ′ is considered. The two last inequalities imply that P ′ satisfies the lower bound test, which
contradicts the fact that P ′ has not been extended, and we obtain the lemma.

Lemma 3. Under Assumption (5), if an o-v path Q such that xQ ⊕ bv
algorithm, it does not satisfy the lower bound test (Low).

xM is considered by the

Proof. Suppose that Assumption (5) is satisfied, and that Q is considered by the algorithm. If
ρ(xQ ⊕ bv ) = 1, path Q does not satisfy the lower bound test and we obtain the result. We now
B
prove that, if ρ(xQ ⊕ bv ) = 0, we have cU
od < c(xQ ⊕ bv ) when Q is considered by the algorithm,
which then implies that Q does not satisfy the lower bound test, which gives the lemma. Suppose
that it is not the case. We place ourselves at the step when Q is considered. As a consequence, Q
minimizes c(xQ ⊕ bv ) among the paths in L. Let P and xM be as in Assumption (5). By definition
B
of P and xM , the hypothesis xQ ⊕ bv xM implies c(xP ) < c(xQ ⊕ bv ) ≤ cU
od when Q is considered.
Lemma 2 implies that there is a subpath P ′ of P in L when Q is considered. By monotonicity of c
we have c(xP ′ ⊕ bv′ ) 6 c(xP ) < c(xQ ⊕ bv ). This contradicts the fact that Q minimizes c(xQ ⊕ bv )
among the paths in L, and gives the lemma.

Lemma 4. Suppose that (a), (b), or (c) is satisfied, then there is a finite number of paths in D
that satisfy the lower-bound test.
10

Proof. In case (a), graph D is acyclic and there is a finite number of paths. We now suppose that
we are in case (b) or (c): let xM be as in one Assumption (5). Lemma 3 ensures that only o-v
paths P such that xP ⊕ bv 6 qM can satisfy the lower bound test. We show the lemma by proving
that there is a finite number of such paths.
As we are in case (b) or (c), any elementary cycle C in D satisfies xC > 0. Thus, given an
elementary cycle C, an elementary path Q, and a vertex v in P , as the resource of C is positive,
Assumption (4) implies that there exists an integer nC,Q,v such that (nC,Q,v xC ) ⊕ xQ ⊕ bv
xM .
As there is a finite number of elementary paths and a finite number of elementary cycles in D, we
can define n to be an integer such that
(6)

(nxC ) ⊕ xQ ⊕ bv

xM

for any elementary cycle C, elementary path Q, and bound bv . Let nC be the number of elementary
cycles in D.
The proof of case (b) relies on the following well-known result.
Any path in a directed graph can be decomposed in a sequence of elementary paths and elementary
cycles.
Suppose that we are in case (b), let P be a path with at least 2nnC |V | arcs and consider such
a decomposition. As an elementary path or an elementary cycle contains at most |V | arcs, this
decomposition contains at least nnC cycles, and thus at least n copies of a given cycle C0 . As
the resource of all arcs are positive by hypothesis of case (b), we have xP > nxC0 . We therefore
have xP ⊕ bv > nxC0 ⊕ bv , and by applying Equation (6) with the empty path as Q, we obtain
x P ⊕ bv
xM . As a consequence, only o-v paths P with fewer than 2nnC |V | arcs can satisfy
xP ⊕ bv 6 qm , and Lemma 3 ensures that there is a finite number of paths that satisfy the lower
bound test in case (b).
We now consider case (c). As the monoid is supposed to be commutative, the resource of a path
does not depend on the order of the sequence of its arcs, but only on its multiset of arcs. The proof
of case (c) relies on the following well-known result.
The multiset of arcs of any path in a directed graph can be decomposed in the union of the sets
of arcs of an elementary path and of several elementary cycles.
Suppose that we are in case (c), let P be a path with at least n|V |(nC + 1) arcs and consider such
a decomposition where Q denotes the elementary path. As by hypothesis of case (c), the operator
⊕ is commutative, the resource of P is entirely defined by the resource of its arcs, independently
of their order. As an elementary path or an elementary cycle contains at most |V | arcs, this
decomposition contains at least nnC cycles, and thus at least n copies of a given cycle C0 . As, by
hypothesis of case (c), all cycles are positive, we have xP ⊕ bv > nxC0 ⊕ xQ ⊕ bv . Equation 6 ensures
that only paths with less than n|V |(nC + 1) arcs can satisfy xP ⊕ bv 6 qm , and Lemma 3 ensures
that there is a finite number of paths that satisfy the lower bound test in case (c).

Proof of Theorem 1. As any path inserted in L is the extension of a previously considered path, a
given path is considered at most once by the algorithm. Thus, Lemma 4 implies the convergence
after a finite number of iterations as only paths satisfying the lower bound test can be extended by
the algorithm.
B
At the end of the algorithm, list L is empty and Lemma 2 ensures that cU
od is a lower bound on
B
the cost of any o-d path satisfying ρ(xP ) = 0. Besides, Step 2 of the algorithm ensures that if cU
od
U
B
is different from +∞, then there is a path P such that c(xP ) = cod and ρ(xP ) = 0. This concludes
the proof.

4.2.2. Label correcting and label dominance algorithms.
Theorem 5. Suppose that the resource of any cycle C in D satisfies xC > 0, then the label
B
correcting algorithm converges after a finite number of iterations, and at the end, if cU
od is finite,
11

B
then cU
od is the cost of a non-dominated optimal solution of the Monoid Resource Constrained
Shortest Path Problem. Otherwise, the problem admits no feasible solutions.

Theorem 6. Suppose that the resource of any cycle C in D satisfies xC > 0, then the label
B
dominance algorithm converges after a finite number of iterations, and at the end, if cU
od is finite,
U
B
then cod is the cost of a non-dominated optimal solution of the Monoid Resource Constrained
Shortest Path Problem. Otherwise, the problem admits no feasible solutions.
Remark that the label correcting and the label dominance algorithms converge under weaker
conditions than those required for the convergence of the generalized A∗ algorithm in Theorem 1.
Lemma 7. Suppose that the resource of any C in D satisfies xC > 0, then if a path P containing a
cycle is considered by the label dominance or the label correcting algorithms, then it does not satisfy
the dominance test (dom).
Proof. As paths in L are added by extension of paths previously in L, we only need to prove the
result for paths ending by a cycle. Let P be such a path, let Q + C be its decomposition into
a path and a cycle, and let v be the common destination vertex of P and Q. By hypothesis, we
have xC > 0. As a consequence, xP = xQ ⊕ xC > xQ . As P is processed, all its subpaths have
been extended by the algorithm, and thus path Q has necessarily been extended. This implies that
either Q or a path Q′ such that xQ′ < xQ ≤ xP is in Mv , and thus P is dominated by a path in
Mv and is therefore not extended.

Proof of Theorem 5. As any path inserted in L is the extension of a previously considered path, a
given path is considered at most once by the algorithm. Thus, as there is only a finite number of
acyclic paths in a graph, Lemma 7 ensures that the algorithm converges after a finite number of
iterations.
B
Step (b) of the algorithm ensures that cU
od is non-smaller than the cost of an optimal solution
of the Monoid Resource Constrained Shortest Path Problem. We now prove that at the
B
end of the algorithm, cU
od is equal to the cost of an optimal solution. Indeed, suppose that it is not
the case. Let P be an optimal solution. Let L be the set of all the paths that have been contained
in L along the algorithm, and for each vertex v in P , let Pov be the subpath of P starting v. Let
v be the last vertex of P such that there is an o-v path Q in L with xQ 6 xPov . It exists because
the empty path Poo is added to L at the beginning of the algorithm, and is therefore in L. Besides,
B
it is not equal to d, as otherwise we would have cU
od ≤ c(xP ). Among the o-v paths dominating
Pov in L, let Q be the first one generated by the algorithm. By definition of Q and as any path
that has been in Mv along the algorithm is in L, there is no path dominating Q in Mv when Q
B
is processed. As cU
od decreases along the algorithm and by hypothesis, when Q is processed we
B
have ρ(xQ ⊕ bv ) ≤ ρ(xPov ⊕ bv ) ≤ ρ(xP ) = 0 and c(xQ ⊕ bv ) ≤ c(xPov ⊕ bv ) ≤ c(xP ) < cU
od . As
a consequence, Q has been extended, and Q + (v, w) is in L, where w be the vertex after v in P .
Besides, we have xQ+(v,w) = xQ ⊕ x(v,w) 6 xPov ⊕ x(v,w) = xPow , which contradicts the definition of
v.

The proof of Theorem 6 is analogous and can be found in [71].
5. Bounding algorithms
We now come to the computation of lower bounds bv on the resource of any v-d path P . We
have already
mentioned in the introduction that such lower bounds bv satisfy bv 6 bopt
v , where
V
opt
opt
bv = P ∈Pvd xP , and Pvd is the set of v-d paths. The bound bv is therefore the best lower
bound. We prove in Proposition 4.11 of [71] the following complexity results on the computation
of bopt
v .
12

Proposition 8. Unless P = N P, there is no polynomial algorithm independent of R that enables
to compute bopt
even when restricted to a commutative monoid with positive resources.
v
The proof is a reduction of the usual resource constrained shortest path problem.
Computing bopt
is therefore difficult in the general case. However, Theorem 9 shows that when
v
can be computed in polynomial time. The remaining of the
⊕ distributes with respect to ∧, bopt
v
section introduces polynomial procedures to compute lower bounds on bopt
v .
5.1. Extended Ford-Bellman algorithm. Let (bnv )n be the sequence of tuples of resources defined recursively as follows.
 n

 bd = 0,

^
(7)

(x(v,u) ⊕ bnu ) for v ∈ V \{d}.
= bnv ∧
b0v = ∞ and bn+1
v


+
(v,u)∈δ (v)

As R is a complete lattice, we can define b∞
v =
(bnv )v∈V

tuple
equation.

. Recall that we have defined

bd = 0,




(8)


b = bv ∧

 v

^

b†

=

^

bnv for each vertex v. Let bn denote the

n∈Z+
†
(bv )v∈V to

be the greatest solution of the following

(x(v,u) ⊕ bu ) for all v ∈ V \{d}.

(v,u)∈δ+ (v)

The existence of a greatest solution of Equation (8) is a direct consequence of the Knaster-Tarski
fixed point theorem applied in the complete product lattice RV . This theorem states that the set of
fixed points of a monotone mapping in a complete lattice is a non-empty complete lattice. Details
on the Knaster-Tarski fixed point theorem can be found in [24]. We underline that both b†v and b∞
v
may be defined only in the completion of R.
Theorem 9. Let ℓ∗ be the length of the longest elementary v-d path. If xC > 0 for each cycle C
in D, then for each vertex and v-d path P , we have
∗

ℓ
opt
b†v 6 b∞
v 6 bv 6 bv 6 x P

(9)

If ⊕ distributes with respect to ∧, then three first inequalities are equalities.
∗

These bounds bℓv are good candidates to be used as bounds bv on the resource xP of all v-d paths P
in the enumeration algorithms of Section 4. Indeed, they can be computed in O(|A|ℓ∗ ) operations ⊕
and ∧ by computing the ℓ∗ first terms of sequence (bn )n using its definition in Equation (7). Besides,
the sequence (bn )n can be interpreted as a generalization of the Ford-Bellman algorithm. Indeed,
when (R, ⊕, 6) = (R, +, ≤), or more generally when (R, ⊕, 6) is a totally ordered group, the
meet x1 ∧ x2 of two resources x1 and x2 is the minimum of x1 and x2 . In that case, the sequence
of Equation (7) corresponds to the successive steps of the Ford-Bellman shortest path algorithm,
and for each integer k, the bound bkv is the value of a shortest o-v path with at most k arcs.
The proof of Theorem 9 relies on two lemmas. The mapping F : RV → RV defined as follows is
useful in the proof.
 ′
b = 0,


 d
^
(10)
F (b) = b′ with

b′v = bv ∧
(x(v,u) ⊕ bu ) for all v ∈ V \{d},


(v,u)∈δ+ (v)
13

where RV denotes the Cartesian product. Note that bn+1 = F (bn ) and that F is isotone by
monotonicity of the operators ⊕ and ∧.
n
Lemma 10. For each vertex v and integer n, we have b†v 6 b∞
v 6 bv .

Proof. A straightforward induction on n based on the isotony of mapping F defined in Equation (10)

gives b†v 6 bnv for all n, which implies that b†v 6 b∞
v .
Lemma 11. The resource bkv is a lower bound on the resource xP of the v-d paths P of length at
most k. When ⊕ distributes with respect to ∧, bkv is the meet of the resources of the v-d paths of
length at most k.
Proof. The result is proved by induction on k. The result for k = 0, i.e. b0v is equal to 0 if v = d and
∞ otherwise, follows from the fact that the only path of length 0 is the trivial path. Let k > 0 be
an integer and suppose the result is true up to k − 1, let v be a vertex, and let P be a v-d path with
ℓ(P ) ≤ k. If ℓ(P ) = 0 then v = d and xP > bkd = 0. Otherwise let (v, u) be the first arc of P and
Q be the subpath of P obtained by removing (v, u) from P . Then ℓ(Q) ≤ k − 1, thus buk−1 6 xQ
which implies x(v,u) ⊕ buk−1 6 x(v,u) ⊕ xQ and finally bkv 6 xP , which gives that bkv is a lower bound
on the resource xP of the v-d paths P of length at most k.
When ⊕ distributes with respect to ∧, we have x1 ⊕ (x2 ∧ x3 ) = (x1 ⊕ x2 ) ∧ (x1 ⊕ x3 ). Suppose
that buk−1 is the meet of the resource of all the v-d paths of length at most k for each vertex u.
Then x(v,u) ⊕ buk−1 is the meet of the resources xP of all the v-d paths P starting by (v, u) such


^
that ℓ(P ) ≤ k. Thus bvk−1 ∧
x(v,u) ⊕ buk−1 is the meet of all the v-d paths of length at
(v,u)∈δ+ (v)

most k.



Proof of Theorem 9. As the resource of any cycle C in D satisfies xC ≥ 0, for each o-d path P ,
is the meet of the resources
there is an elementary o-d path P ′ such that xP ′ 6 xP . Hence, bopt
v
of all the elementary o-d paths. As the length of any elementary v-d path is non greater than ℓ∗ ,
∗
∗
Lemma 11 implies that bℓv 6 xP for all elementary v-d paths P , and thus bℓv ≤ bopt
v . Lemma 10
then gives Equation (9).
Suppose now that ⊕ distributes with respect to ∧. We already mentioned that xC > 0 for all
cycles C implies that bopt
v is the meet of all the elementary o-d paths. As the length of an elementary
∗
ℓ∗
ℓ∗
path is non greater than ℓ∗ , Lemma 11 implies that bℓv = bopt
v . This implies that F (b ) = b and
∗
∗
bℓv is a solution of Equation (8). Thus, bℓv = b†v , which gives the result.

Remark 2. If there are cycles with negative resources in D, then Theorem 9 remains true provided
that, first, P is an elementary v-d path, and second, bopt
is defined as the meet of the resources of
v
all the elementary v-d paths. Besides, Lemma 11 implies that b†v = b∞
v = −∞.
Remark 3. The sequence (bn )n is the sequence used in the constructive proof by Cousot & Cousot
[23] of the Knaster-Tarski fixed point theorem for mapping the F defined in Equation (10). Given a
†
∞
topology and some
on R, it can be proved that bnv converges to b∞
v and bv = bv .
V weak assumptions
n
∞
∞
∞
The inequality (v,u)∈δ+ (v) (x(v,u) ⊕bu ) 6 bv is easy to prove: indeed, as x(v,u) ⊕bu 6 x(v,u) ⊕bu for
V
V
each arc (v, u) in δ+ (v) and for all n in Z+ , we have (v,u)∈δ+ (v) (x(v,u) ⊕ b∞
u )6
(v,u)∈δ+ (v) (x(v,u) ⊕
V
∞
n
n+1
bu ) = bv
for all n ∈ Z+ , which gives the result. The inequality (v,u)∈δ+ (v) (x(v,u) ⊕ b∞
u ) > bv
requires a transfinite induction.
5.2. Generalized Dijkstra algorithm for faster bound computations. In all our numerical
ℓ∗
experiments, we have observed that practically, b†v = b∞
v = bv . Therefore, the bounds we obtain if
∗
we compute b†v are as good as those we obtain if we compute bℓv . When ⊕ distributes with respect
14

to ∧, computing b†v amounts to solving the algebraic path problem. There are algorithms for the
algebraic path problem that are practically much more efficient than the generalized Ford-Bellman
algorithm of the previous section. We now adapt one of these algorithms, which generalizes Dijkstra
algorithm, to our lattice ordered monoid setting.
A resource b̃W
v and an integer ñv are attached to each vertex v and updated during the algorithm.
Initially, b̃v = R, which may be defined only in the completion of R, and ñv = +∞ for each vertex
v 6= d, and b̃d = 0. During the algorithm, a queue L of vertices “to be extended” is maintained.
Initially, the queue L contains only d. The algorithm ends when L is empty. While L is not empty
and the minimum ñv over all vertices v is non greater than ℓ∗ , where ℓ∗ is the maximum length of
an elementary path ending in d, the following operations are repeated:
• Extract from L a vertex v with minimum ñv .
• For each arc (u, v) in δ− (v), extend v along (u, v): if b̃u x(u,v) ⊕ b̃v , then
– Update ñu to min(ñu , 1 + ñv ).
– Update b̃u to b̃u ∧ (x(u,v) ⊕ b̃v ).
– Add u to L (if it is not already present).
• Set ñv = +∞.
Proposition 12. This algorithm terminates after at most ℓ∗ |V | iterations, where ℓ∗ is the maximum
length of an elementary path ending in d. The value bv of b̃v at the end of the algorithm is equal to
∗
bℓv for each v ∈ V . If L is empty at the end of the algorithm, then bv = b†v for all vertices v.
When ⊕ distributes with respect to ∧, at the end of the algorithm, we have L = ∅ and bv = b†v
for all vertices v.
The proof of Proposition 12 is technical but not difficult and relies on Theorem 9. The interested
reader can find it in Proposition 4.15 in [71]. If the complexity bounds we obtain are identical to
those of the algorithm of the previous section, this algorithm is practically faster in practice, and
should be preferred. We use it in our numerical experiments.
In the first step, we can extract a vertex v that minimizes a key function φ(b̃v ) instead of one
that minimizes ñv . Using the key φ(x) = x, our algorithm corresponds to Dijkstra algorithm when
(R, ⊕, 6) = (R+ , +, ≤). When an alternative φ(b̃v ) is used instead of ñv , the algorithm ends only
when L is empty, and we loose the convergence of Proposition 12. However, in practice, the list L
is always empty after a number of iterations that is not much larger than |V |. In our numerical
experiments, we use the ratio
number of vertices extended before convergence
(11)
γ=
|V |
to evaluate how fast the algorithm converges. As each vertex v such that there exists a v-d path
is extended at least once, if there is a v-d path P for each vertex v in D, the ratio γ is necessarily
non smaller than 1. With a carefully chosen φ, convergence is fast: the worst γ encountered in
the numerical experiments is 2.6. Using ñv , we have obtain ration γ up to 20 times larger on the
instance.
When ⊕ distributes with respect to ∧, we have a convergence result for arbitrary φ. Indeed, in
that case, the only difference between our algorithm and the one proposed by Mohri [65] for the
algebraic path problem is that the vertex picked-up in L at each iteration is arbitrarily chosen.
Mohri [65] shows that, after a finite but possibly exponential number of iterations, L is empty and
the algorithm terminates.
Remark 4. When digraph D is acylic, the bounds b†v can easily be obtained with an even faster
algorithm. Indeed, using a topological ordering, they can be computed iteratively directly from the
dynamic programming equation (8).
15

Name
road
square
long
wide
acyc
rand

Generator
extraction
grid
grid
grid
acyclic
random

Brief description
Road networks with a given number of vertices.
Square grid of size m × m
Long grid of size 16m × m
Wide grid of size m × 16m
Acyclic graph with n vertices v1 , . . . , vn and 5n arcs (vi , vj ) with i < j
Hamiltonian cycle with n vertices and 5n − n chords.
Table 2. Summary of the families of graphs used

6. Numerical experiments on usual resource constrained shortest path problem
In this section, we test numerically our enumeration algorithms on the usual resource constrained
shortest path problem with k = 1 and k = 10 constraints,
X
wa0 ,
min
P ∈Pod
a∈P
X
(12)
wai ≤ W i .
s.t.
a∈P

where wai ∈ R for each arc a and i ∈ {0, . . . , k}, thresholds W i ∈ R for i ∈ {1, . . . , k}, and Pod is the
set of o-d paths. We model it as a Monoid Resource Constrained Shortest Path Problem
with resources x ∈ Rk+1 ,
c(x) = w0

and ρ(x) = max(1(W i ,+∞)(wi )),
i

where x =
and 1(W i ,+∞) is the indicator function of the interval (W i , +∞). As
we focus on instances with positive resources, Theorems 1.(b), 5, and 6 ensure respectively the
convergence of the generalized A∗ , label correcting, and label dominance algorithms.
(w0 , . . . , wk ),

6.1. Instances. We use four families of graphs: road networks, acyclic graphs, grids, and random
graphs. The three last families of graphs are used by Cherkassky et al. [21] in their experimental
study of algorithms for the Shortest Path Problem. We have used adapted versions of their
generators spgrid, sprand, and spacyc to produce instances of these families. The adaptation
consists in the insertion of a destination vertex. Among others, these four family of graphs have
been used by Dumitrescu & Boland [28] to test the different Resource Constrained Shortest
Path Problem algorithms available in the literature. We now describe them.
o

o

o

d

d

d

Figure 1. From left to right, a grid graph of width 3 and depth 2, a random graph
with 5 vertices and 9 arcs, and an acyclic graph with 4 vertices and 7 arcs.
The road network graphs have been extracted from the Rome and the San Francisco Bay Area
instances of the Dimacs challenge [1], and the origin and the destination vertices have been chosen
16

far way from each other. Details on how these instances have been extracted are available in [71].
A grid graph of length ℓ and width m is composed of ℓ layers of m vertices. Each layer i is a
Hamiltonian cycle vi,1 , . . . , vi,m with arcs in both direction. Each vertex vi,j of layer i ∈ [ℓ − 1] is
connected to the corresponding vertex vi+1,j in the next layer. An origin vertex o and a destination
vertex d are added. There is an arc between the origin o and each vertex of the first layer, and
an arc between each vertex of the last layer and the destination. A random graph with n vertices
and m ≥ n arcs is composed of n vertices on a Hamiltonian cycle and randomly generated chords
on that cycle. Finally an acyclic graph with n vertices and m ≥ n arcs contains a path v1 , . . . , vn ,
and randomly generated arcs (vi , vj ) with i < j. The origin and the destination are respectively
v1 and vn . Figure 1 illustrates a grid, a random and an acyclic graph. We consider random and
acyclic graphs with n vertices and 5n arcs. We have obtained similar results are obtained when
using digraphs with hn arcs with h between 2 and 50 [71]. Table 2 provides a summary of the main
characteristics of these families of graph.
The weight wa0 of an arc a of a road network instance is its length in kilometers. All the other
weights wai of each type of graphs are randomly chosen using an uniform distribution on [1, 100].
We use the same technique as [28] to ensure the existence of an optimal
X we first search an
Xpath:
o-d path Pc of minimum cost, then find an o-d path Pw minimizing
wai , and finally set
a∈Pw i∈k

W i = (1 − λ)wPi w + λ max(wPi , wPi w ), where λ ∈ [0, 1] enables to choose the constraint strength. As
the relative performances of the algorithms do not depend much of λ, we always use λ = 0.5 in this
paper. A study of the influence of the constraint strength λ is available in [71].

6.2. Candidate paths, weaker bounds, and algorithms tested. We have tested the three
enumeration algorithms introduced in Section 4. When bounds bv are needed, they are computed
P
in a preprocessing using the generalized Dijkstra algorithm of Section 5.2, with x 7→ ki=0 wi as
key function φ of Equation (11), where x = (w0 , . . . , wk ). On acyclic instances, we use instead the
algorithm of Remark 4. We have also tested two variants of the enumeration algorithms, that we
now introduce.
Candidate paths. A standard technique to improve the performance of resource constrained shortest
path algorithms [28] relies on the use of v-d paths Qv that are likely to be subpaths of optimal v-d
paths. We find such paths Qv during the preprocessing: we use Dijkstra algorithm to find the v-d
P
P
path P that minimizes a∈P ki=0 wai . When an o-v path P satisfies the test of an enumeration
B
algorithm, before extending it, we test if P + Qv is a feasible solution of cost smaller than cU
od , and
B
update cU
od if it is. This procedure may enable the algorithm to find faster feasible o-d paths of good
B
quality. The upper bound cU
od is therefore likely to decrease faster, which enables to strengthen
the lower bound test and thus reduce the total number of paths enumerated by the algorithm.
Weaker bounds. To identify the respective contributions of the lower bound test (Low) and of the
keys c(xP ⊕ bv ) to the performance of the generalized A∗ and the label correcting algorithms, we
also provide numerical results for “downgraded” versions of the generalized A∗ and label correcting
algorithms where c(xP ) is used instead of c(xP ⊕ bv ) as key.
6.3. Experimental setting. The numerical experiments are performed on a Macbook Pro of 2012
with four 2.5 Ghz processors and 4 Gb of ram. The algorithms are not parallelized. The limiting
parameter for the algorithms is the memory available. Therefore, we stop the algorithms if the list
L of candidate paths contains more than 1e+05 elements. We also stop the label correcting and
the label dominance algorithms if the number of paths in the union of the sets of non dominated
paths Mv is larger than 1e+05. As we mention in Remark 1, the minimum c(xP + bv ) on the paths
P in L when the algorithm is stopped provides a lower bound on the cost of an optimal solution,
17

square

road
102

A∗
cor.
dom.

101
100

CPU time (s)

CPU time (s)

102

10−1
10−2
10−3

A∗
cor.
dom.

101
100
10−1
10−2
10−3

102

103
104
Vertices

105

102

103
104
Vertices

long

wide
102

A∗
cor.
dom.

101
100

CPU time (s)

CPU time (s)

102

10−1
10−2
10−3

A∗
cor.
dom.

101
100
10−1
10−2
10−3

102

103
104
Vertices

105

102

103
104
Vertices

acyc
102

A∗
cor.
dom.

101
100

105

rand
102
CPU time (s)

CPU time (s)

105

10−1
10−2
10−3

A∗
cor.
dom.

101
100
10−1
10−2
10−3

102

103
104
Vertices

105

102

103
104
Vertices

105

Figure 2. Usual RCSP with k = 1 constraint
B
and cU
od provides an upper bound. These two bounds are used to compute a gap on instances that
cannot be solved to optimality.

6.4. Numerical results. We now come to numerical results on Problem (12). Figure 2 plots the
performance of our algorithms for each family of instances. As there is only one degree of freedom
18

square
102

101

101

CPU time (s)

CPU time (s)

road
102

100
10−1
A∗
cor.
dom.

10−2
10−3
102

103
104
Vertices

100
10−1
A∗
cor.
dom.

10−2
10−3

105

102

102

101

101

100
10−1
10−2

A∗

10−3
103
104
Vertices

100
10−1
A∗
cor.
dom.

10−2
10−3

cor.
102

105

102

acyc

100
10−1
10−2
10−3

101
100
10−1
A∗
cor.
dom.

10−2
10−3

102

103
104
Vertices

105

rand

A∗
cor.
dom.

101

103
104
Vertices

102
CPU time (s)

CPU time (s)

102

105

wide

102
CPU time (s)

CPU time (s)

long

103
104
Vertices

105

102

103
104
Vertices

105

Figure 3. Usual RCSP with k = 10 constraints
in the definition of instances of each family of graphs, we can identify an instance of a given family
by its number of vertices. For each of our enumeration algorithms in Table 1, there is a curve that
gives the total CPU time as a function of the number of vertices. This total CPU time includes the
computation of bounds in the preprocessing. Curves A∗ , cor., and dom. correspond respectively
to the generalized A∗ , the label correcting and the label dominance algorithm. We only add points
19

corresponding to instance solved to optimality. Therefore, the longer and the lower a curve, the
better the performance of the corresponding algorithm. Figure 3 provides the same results for the
problem with k = 10 constraints.
Tables 3 and 4 provide detailed results on some difficult instances with respectively k = 1 and
k = 10 constraints. The first columns provide the instance considered, the number of vertices, and
the number of arcs in the instance. The next column provides the enumeration algorithm tested.
Again, A∗ corresponds to the generalized A∗ algorithm, cor. to the label correcting algorithm, and
dom. to the label dominance algorithm. The suffix CP indicates that the variant of the algorithm
with candidate paths is considered, and the suffix K that c(xP ) is used instead of c(xP ⊕ bv ).
The ratio γ of the next column is the one of Equation (11), and indicates the performance of
the generalized Dijkstra algorithm. The lower γ, the better is the performance of this algorithm.
The next column provides the percentage of the total CPU time spent in the preprocessing. For
the standard version of the label dominance algorithm, no preprocessing is needed, while for the
generalized A∗ and the label correcting algorithm, this preprocessing time indicates the time needed
to compute lower bounds. When candidate paths are computed, the percentage after the symbol
+ indicates the percentage of the total CPU time spent computing them. The two next columns
provide the number of paths extended and the number of paths cut. The next column indicates
the percentage of paths cut by the dominance test, the remaining being cut by the lower bound
test. The other algorithm are not concerned as they use a single test to cut paths. The column ℓ
provides the number of arcs in the solution returned. The next column provides the gap between
the lower and the upper bound when the algorithm is stopped before convergence, and the last
column provides the total CPU time.
Difficulty of the instances. We first, remark that the difficulty of the instances is mainly linked
to the number of arcs ℓ in an optimal solution. The smaller this number of arcs, the easier the
instances. The acyclic and the random graph instances are therefore the easiest to solve, and any
version of our algorithms can solve them well. Then comes the road network instances and the
wide grids. The most difficult instances are the square and the long grids.
Relative performance of the algorithms. The label dominance algorithm has the best performances
on the easy instances, i.e. the acyclic and the random instances. If we look at Tables 3 and 4, we
can see that on these instances, most of the computation time of the generalized A∗ and of the
label correcting algorithms is spent in the preprocessing, gathering information that is not needed
as the problem is easy to solve. On the other instances, the label correcting algorithm has the best
performance. Using lower bounds in addition to the dominance test is always interesting on difficult
instances. The relative performance of the generalized A∗ algorithm and of the label correcting
algorithm depends on the family of instances considered.
Influence of dimension and relative performance of the algorithms. Instance with ten constraints
are more difficult to solve than instances with one constraint. We also underline the fact that
the relative performance of the algorithms on instances with ten constraints is different from the
performance in the case of instances with one constraint. Indeed, the label dominance algorithm
exhibits poor performances in that setting, while the generalized A∗ performs almost as well as the
label correcting algorithm. We can also see in the tables that the percentage of paths cut by the
dominance test in the label correcting algorithm decreases with the dimension. Section 8 explains
why the performance of the dominance test decreases with the dimension, and provides techniques
to increase the performance of our algorithms on problems in large dimension.
Influence of keys. Figure 4 provides the performances of the generalized A∗ and of the label correcting algorithms with different keys to solve Problem (12) on wide grids with k = 1 constraint. Plain
lines correspond to the key c(xP ⊕ bv ) and dashed lines to the key c(xP ). On instances with k = 1
constraint, There is no-dashed line corresponding to the generalized A∗ because only one instance
was solved by the algorithm with key c(xP ). When bounds bv have been computed for the lower
20

Instance
road20

|V |
20000

|A|
55180

road50

50000

138112

square100

10002

30100

5122

15376

wide100

25602

78400

acyc10000

10000

50000

rand10000

10000

50000

long20

Alg.
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP

γ
1.6
1.6
1.6
1.6
1.6
1.6
–
–
2.7
2.7
2.7
2.7
2.7
2.7
–
–
1.6
1.6
1.6
1.6
1.6
1.6
–
–
1.6
1.6
1.6
1.6
1.6
1.6
–
–
1.5
1.5
1.5
1.5
1.5
1.5
–
–
1
1
1
1
1
1
–
–
2.6
2.6
2.6
2.6
2.6
2.6
–
–

Preproc.
15%
9%+10%
16%
64%
36%+37%
1%
–
– +0%
43%
25%+21%
42%
17%
12%+10%
1%
–
– +1%
9%
7%+6%
10%
3%
2%+2%
0%
–
– +0%
5%
3%+4%
6%
1%
1%+1%
0%
–
– +0%
97%
51%+47%
19%
97%
50%+48%
17%
–
– +10%
97%
32%+66%
96%
95%
29%+70%
95%
–
– +50%
99%
68%+32%
89%
99%
65%+35%
89%
–
– +12%

Ext.
54131
55762
71402
7899
7852
1128194
1505043
1505043
44613
44613
71420
153191
153221
1570406
1557409
1557409
58533
58533
50298
168036
168016
1488625
1493529
1493529
51961
51961
49996
186590
186590
1558600
1558884
1558884
107
104
67108
98
95
105678
148903
148564
12
7
22
12
7
22
5918
5909
41
37
2078
41
37
1957
67626
67378

Cut.
14121
17579
4
8306
8228
793042
995349
995349
3710
3710
0
80428
80488
1047263
1029848
1029848
17162
17162
690
155429
155396
981538
973040
973040
3933
3933
3
93798
93798
1025370
1024684
1024684
1809
1806
35812
1785
1782
99041
83968
83772
63
51
127
63
51
127
1413
1407
161
144
9393
161
144
7889
41871
41653

Dom.
–
–
–
44%
44%
92%
–
–
–
–
–
96%
96%
99%
–
–
–
–
–
58%
58%
99%
–
–
–
–
–
96%
96%
100%
–
–
–
–
–
0%
0%
58%
–
–
–
–
–
0%
0%
0%
–
–
–
–
–
0%
0%
2%
–
–

ℓ
–
252
–
243
243
243
–
243
–
494
–
–
478
–
–
462
–
131
–
118
118
–
–
131
–
438
–
–
419
–
–
437
21
21
–
21
21
21
21
21
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10

Gap
∞
2.5%
∞
opt
opt
opt
∞
21.0%
∞
13.6%
∞
∞
5.5%
∞
∞
431.2%
∞
30.4%
∞
opt
opt
∞
∞
25.1%
∞
58.6%
∞
∞
33.9%
∞
∞
202.3%
opt
opt
∞
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt

CPU (s)
3.08e-01
4.91e-01
3.21e-01
7.32e-02
1.30e-01
7.17e+00
1.06e+01
1.20e+01
4.42e-01
7.21e-01
4.71e-01
1.12e+00
1.60e+00
1.37e+01
1.21e+01
1.44e+01
2.69e-01
4.36e-01
2.62e-01
7.84e-01
1.01e+00
1.17e+01
1.05e+01
1.21e+01
2.66e-01
4.17e-01
2.36e-01
1.12e+00
1.45e+00
2.56e+01
2.36e+01
2.62e+01
7.23e-02
1.42e-01
4.22e-01
7.36e-02
1.40e-01
4.51e-01
4.62e-01
6.79e-01
9.20e-03
2.59e-02
9.14e-03
8.07e-03
2.61e-02
1.01e-02
1.99e-02
3.44e-02
7.48e-02
1.10e-01
1.05e-01
8.02e-02
1.13e-01
1.17e-01
2.26e-01
3.31e-01

Table 3. Standard resource constrained shortest path with k = 1 resource constraint
21

Instance
road20

|V |
20000

|A|
55180

road50

50000

138112

square50

2502

7550

long5

1282

3856

wide100

25602

78400

acyc10000

10000

50000

rand10000

10000

50000

Alg.
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP

γ
2.1
2.1
2.1
2.1
2.1
2.1
–
–
2.5
2.5
2.5
2.5
2.5
2.5
–
–
2.2
2.2
2.2
2.2
2.2
2.2
–
–
2.2
2.2
2.2
2.2
2.2
2.2
–
–
2.1
2.1
2.1
2.1
2.1
2.1
–
–
1
1
1
1
1
1
–
–
6.3
6.3
6.3
6.3
6.3
6.3
–
–

Preproc.
7%
6%+4%
5%
24%
18%+13%
25%
–
– +0%
47%
30%+19%
43%
12%
10%+7%
1%
–
– +0%
3%
2%+2%
3%
0%
0%+0%
1%
–
– +0%
2%
1%+1%
2%
0%
0%+0%
0%
–
– +0%
85%
59%+27%
83%
82%
60%+27%
83%
–
– +7%
91%
47%+51%
89%
83%
42%+56%
84%
–
– +77%
100%
88%+12%
100%
100%
88%+12%
100%
–
– +40%

Ext.
253214
253214
340302
45066
45034
45177
1103713
1103713
56159
56159
71421
64093
64093
466563
477467
477467
81247
81247
56775
141815
141815
179777
168274
168274
64642
64642
52538
101197
101197
196476
195012
195012
6702
6698
6711
6413
6409
6421
108128
108128
6
2
6
6
2
6
211
210
13
9
14
13
9
14
7183
7177

Cut.
323551
323551
556262
66374
66317
66416
603707
603707
9481
9481
1
17783
17783
303206
307503
307503
62539
62539
13594
103501
103501
93497
78864
78864
29295
29295
5087
52388
52388
98269
96678
96678
14999
14994
15017
13731
13726
13745
38179
38179
24
12
24
24
12
24
0
0
56
37
57
56
37
57
160
159

Dom.
–
–
–
21%
21%
21%
–
–
–
–
–
53%
53%
99%
–
–
–
–
–
39%
39%
89%
–
–
–
–
–
48%
48%
99%
–
–
–
–
–
3%
3%
3%
–
–
–
–
–
0%
0%
0%
–
–
–
–
–
0%
0%
0%
–
–

ℓ
–
221
–
221
221
221
–
221
–
400
–
–
400
–
–
400
–
52
–
–
52
–
–
52
–
81
–
–
81
–
–
81
17
17
17
17
17
17
–
17
5
5
5
5
5
5
5
5
7
7
7
7
7
7
7
7

Gap
∞
55.3%
∞
opt
opt
opt
∞
373.3%
∞
18.1%
∞
∞
18.0%
∞
∞
1321.4%
∞
57.5%
∞
∞
49.3%
∞
∞
376.1%
∞
67.9%
∞
∞
62.9%
∞
∞
599.7%
opt
opt
opt
opt
opt
opt
∞
314.6%
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt
opt

CPU (s)
9.75e-01
1.23e+00
1.50e+00
3.02e-01
3.82e-01
2.96e-01
3.48e+01
3.91e+01
4.59e-01
7.12e-01
5.02e-01
1.80e+00
2.10e+00
2.91e+01
2.67e+01
2.76e+01
3.46e-01
4.73e-01
2.85e-01
2.41e+00
2.64e+00
2.00e+00
1.64e+00
1.84e+00
2.79e-01
4.02e-01
2.39e-01
1.70e+00
1.90e+00
3.06e+00
2.81e+00
3.03e+00
1.22e-01
1.81e-01
1.38e-01
1.25e-01
1.74e-01
1.40e-01
4.92e-01
6.69e-01
2.79e-03
4.97e-03
2.46e-03
2.52e-03
4.93e-03
2.63e-03
1.08e-03
3.95e-03
2.41e-01
2.58e-01
2.47e-01
2.32e-01
2.67e-01
2.50e-01
4.18e-02
8.55e-02

Table 4. Standard resource constrained shortest path with k = 10 resource constraint
22

k=1

k = 10
102

A∗
cor.
dom.

101
100

CPU time (s)

CPU time (s)

102

10−1
10−2
10−3

101
100
10−1
A∗
cor.
dom.

10−2
10−3

102

103
104
Vertices

105

102

103
104
Vertices

105

Figure 4. Influence of keys on usual RCSP with k = 1 and k = 10 constraint
k=1

k = 10
102

A∗
cor.
dom.

101
100

CPU time (s)

CPU time (s)

102

10−1
10−2
10−3

101
100
10−1
A∗
cor.
dom.

10−2
10−3

102

103
104
Vertices

105

102

103
104
Vertices

105

Figure 5. Influence of candidate paths on usual RCSP with k = 1 and k = 10 constraint

bound test, is is always more interesting to use the key c(xP ⊕ bv ) than the key c(xP ). In practice,
the key has an influence only on weakly constrained problems with k = 1, where the keys identify
which partial solutions are the most promising. On more constrained problems with k = 10, we
are not able to guess which partial solution will lead to a feasible solution, and the importance of
keys decreases.
Influence of candidate paths. Figure 5 shows the influence of candidate paths on the performances
of the algorithms on the wide grid instances with k = 1 and k = 10 constraints. Plain lines correspond to algorithms without candidate paths, and dashed lines to algorithms with candidate paths.
Candidate paths are not interesting on instances that can be solved to optimality: the decrease in
B
cU
od is not large enough to enable to reduce significantly the number of paths enumerated. Therefore, the use of candidate paths only slow the algorithm by requiring a preprocessing and slowing
each step of the enumeration. On the contrary, we can see in Tables 3 and 4 that candidate paths
become interesting on difficult instances: they enable to find feasible o-d paths of small cost along
the algorithm, and thus to obtain a smaller gap.
23

7. Stochastic path problems: models and numerical experiments
We now come to the stochastic path problems mentioned P
in the introduction. We suppose to
have a random variable ξa for each arc a, and denote ξP = a∈P ξa for each path P . Given an
origin o and a destination v, a stochastic path problem typically seeks an o-d path P minimizing a
probability functional µ
!
X
ξa .
min µ
P ∈Pod

a∈P

Three types of probability functionals µ are of specific interest. First, given τ ∈ R, the functional
P(· > τ ) enables to model the probability of arriving on time at an event starting at τ . Second, utility theory considers the expectation E(f (·)) of a non-decreasing cost function f . Third, stochastic
optimization often deals with risk measures [6], i.e. probability functionals µ that are normalized,
µ(0) = 0, monotone with respect to the almost sure order, i.e. µ(ξ) ≤ µ(ξ̃) if ξ ≤ ξ˜ almost surely,
and invariant by translation µ(ξ + c) = µ(ξ) + c for all c ∈ R. Alternative definitions of risk
measures exist in the literature, but all assume the monotonicity with respect to the almost sure
order. Stochastic constraints are typically of the form
µ′ (ξP ) ≤ α,
where µ′ is a probability functional. The most common stochastic constraints are probability
constraints of the type
P(ξP > τ ) ≤ α.
We now provide techniques to model stochastic path problems within the Monoid Resource
Constrained Shortest Path Problem framework. To that purpose, with introduce several
lattice ordered monoids of random variables. On a given lattice ordered monoid (R, ⊕, 6), we can
model any stochastic path problem such that the functions
c(ξ) = µ(ξ) and

ρ(ξ) = 1(α,∞] (µ′ (ξ))

are isotone with respect to 6. On the lattice ordered monoids we consider, this property is satisfied
by any probability functional µ or µ′ that is monotone with respect to the almost sure order. This
is notably the case of all the probability functionals µ and µ′ of interest aforementioned.
7.1. General case. Given arbitrary random variables ξa , the vector space of the random variables
ξa endowed with the addition and the almost sure order is a lattice ordered monoid. Indeed, the
˜ and ξ ′ such
addition is compatible with the almost sure order, as given three random variables ξ, ξ,
′
′
˜
that ξ ≤ ξ a.s., we have ξ + ξ ≤ ξ̃ + ξ a.s.. Besides, the almost sure order induces a lattice structure
and the meet of two random variables is their essential-infimum, which is unique up to a.s. equality.
Practically, to be able to make the computations, we sample from the initial random variables
a finite number N of scenarios ω1 , . . . , ωN . The sampled distribution can therefore be encoded
as elements of RN . The a.s. order between the sampled distributions becomes the component by
component order on RN , and the essential infimum becomes the minimum:
˜ i )).
(ξ ∧ ξ̃)(ωi ) = min(ξ(ωi ), ξ(ω
The lattice ordered monoid we use is therefore (RN , +, ≤). We provide in [71] upper bounds on
the error committed due to sampling that are exponential in the number of samples.
7.2. Independent random variables. When the random variables ξa are independent, we can
work on their distributions. Indeed, the distribution of the sum of two random variables is their
convolution product ∗. In that case, we can use the usual stochastic order ≤st which is such that
ξ ≤st ξ̃

if

P(ξ ≤ t) ≥ P(ξ̃ ≤ t) for all t ∈ R.
24

The usual stochastic order is compatible with the convolution product. Practically, we work on the
set M of distributions of random variables with finite support on Z. We denote Fξ the cumulative
distribution of ξ. The usual stochastic order induces a lattice structure on M, and the meet ξ ∧st ξ̃
of two random variables ξ and ξ̃ is defined through its cumulative distribution:
Fξ∧st ξ̃ = max(Fξ , Fξ̃ ).
We can therefore work with resources in the lattice ordered monoid (M, ∗, ≤st ). A probability
functional µ is version independent if µ(ξ) depends only on the distribution of ξ. As we work on
the random variables distributions, we can use only version-independent probability functionals.
All the probability functionals of interest mentioned at the beginning of the section are version
independent. And as they are also monotone with respect to the almost sure order, the following
proposition ensures that we can deal with them using the Monoid Resource Constrained
Shortest Path Problem framework with resource in (M, ∗, ≤st ).
Proposition 13. A version independent probability functional that is monotone with respect to the
almost sure order is monotone with respect to the usual stochastic order ≤st .
The proof of this proposition is straightforward and can be found in [71]. The converse statement
holds even for probability functionals that are not version independent as the usual stochastic order
is coarser than the almost sure order, i.e.
ξ ≤ ξ˜ a.s. implies ξ ≤st ξ̃.
As ≤st is coarser than the almost sure order, the lower bounds we obtain using our algorithms
with resources in (M, ∗, ≤st ) are tighter than the one we obtain using the sampling approach of
the previous section with resources in (RN , +, ≤st ). Therefore, when the ξa are independent, it is
always more interesting from an algorithmic point of view to use resources in (M, ∗, ≤st ) than to
use the sampling approach.
We can also use parametric families of distributions that are stable by convolution product
instead of discrete distributions. We notably introduce in [71] lattice ordered monoid structures
to model the families of distributions considered in the literature: normal distributions, gamma
distributions, and Cauchy distributions.
Remark 5. (M, ∗, ≤st ) is an example of lattice ordered monoid that is not an indempotent semiring,
as ∗ does not distribute with respect to ∧st .
7.3. Numerical results. We now benchmark our algorithms on two stochastic path problems.
We consider here independent random variables ξa . We provide in [71] numerical results on the
same problems and graph instances, but with sampled distributions of non-independent random
variables ξa .
The Condition Value-at-Risk [74] of level β ∈ [0, 1) is
Z 1
1
VaRα (ξ)dα,
CVaRβ (ξ) =
1−β β

where VaRα (ξ) = inf {t|P(ξ ≤ t) ≥ α}. Intuitively, the Conditional Value at Risk of level β can be
interpreted as the expectation in the β worst case. Parameter β enables to choose a level of risk
awareness. As the conditional value at risk is one of the most popular risk measures in stochastic
optimization, the first problem on which we test our algorithms is
!
X
ξa .
(13)
min CVaRβ
P ∈Pod

a∈P

For instance, the optimal solution with β = 0.05 of this problem is the best itinerary for a commuter
going to work: it gives an itinerary that is still good the most congested day of the month.
25

The second problem we consider is a shortest path problem with deterministic cost and probability constraints:
X
wa ,
min
P ∈Pod
a∈P
!
(14)
X
ξa > τ ≤ α.
s.t. P
a∈P

where τ ∈ R and wa are weights in R. This problem can be interpreted as a truck delivery problem.
The objective is to find a path of minimum cost among those that guarantee a certain probability
of arriving on time.
On both problems, we use independent distributions ξa , and model them as elements of the
lattice ordered monoid (M, ∗, ≤st ). The resources of Problem (13) therefore belong to M, and those
of Problem (14) to R × M endowed with the product sum and order. As we focus on instances with
positive resources, Theorems 1.(b), 5, and 6 ensure respectively the convergence of the generalized
A∗ , label correcting, and label dominance algorithm.
square

road
103

A∗
cor.
dom.

102

CPU time (s)

CPU time (s)

103

101
100
10−1
10−2

102

103

104

101
100
10−1
10−2

105

A∗
cor.
dom.

102

102

Vertices

103
104
Vertices

long

wide
103

A∗
cor.
dom.

102

CPU time (s)

CPU time (s)

103

101
100
10−1
10−2

102

103
104
Vertices

105

101
100
10−1
10−2

105

A∗
cor.
dom.

102

102

103
104
Vertices

105

Figure 6. Stochastic Shortest Path Problem with generic distributions and β = 0.05
Instances. We build instances by generating resources on the road, square grid, long grid, and wide
grid digraphs introduced in Section 6. We have generated generic discrete distributions with finite
support as follows. First, the length of the arcs are rescaled to be non greater than 200. Second, their
26

support size is chosen as 10 plus a random integer between 0 and the scaled length of the arc. For
each t in the support, a weight P̃(ξ = t) is randomly chosen using a uniform distribution on [0, 1] for
each t. The distributions P̃ is then normalized to obtain a probability distribution P. The measure
is then normalized to obtain a probability distribution. The probability functionals P(· > τ )
and CVaRβ have been tested for different values of τ and β. We have chosen τ as follows: a first
parameter τ −1 is chosen between 0 and 1, then τ is chosen as the smallest t such that P(bo > t) ≤ τ −1
where bo is the bound provided by the bounding algorithm for the origin vertex o. It enables to
obtain a threshold τ such that the P(ξP > τ ) ≃ τ −1 where ξP is the resource of an optimal path.
To choose the parameter α, we use the same technique
as the one we used in Section 6 to choose
 

  P
P
P
i
ξ
>
τ
)
/2 ,
ξ
>
τ
+
P(
ξ
>
τ
,
P
the parameter ω : we set α = max P
a∈Pρ a
a∈Pc a
a∈Pc a
P
where Pc is an o-d path P with minimum
a∈P wa , and Pρ is an o-d path P with minimum
P
P
a∈P ξa > τ .
For both problems, numerical experiments on the same digraphs but with truncated and discretized lognormal distributions are available in [71]. The performances of the algorithms with
these distributions are similar.
Experimental setting. On medium and large instances, the order of magnitude of the number of
non-zero terms in the support of ξ is a few thousands. Storing the resources therefore requires more
memory than in the usual RCSP case. We have therefore set a maximum size of 1e + 04 for the list
of candidate paths L. To avoid spending too much time in the convolution products, we compute
them using a Fast Fourier Transform. To that purpose, we use the Fast Fourier Transform C++
library kissFFT [14]. Concerning the computation of lower bounds, we use the generalized Dijkstra
algorithm of Section 5.2. We use ξ 7→ E(ξ) as key function φ when solving Problem (13), and
(w, ξ) 7→ w + E(ξ) when solving Problem (14).
Non-constrained problem. Figure 6 and Table 5 provide numerical results for Problem (13) with
β = 0.05. They can be read like the figures and tables of Section 6. This non-constrained stochastic
problem is fairly easy to solve: large instances are solved in reasonable time. The main limit on
the size of the instances we can solve is the memory needed to solve the instance. On Figure 6, we
can see that our three algorithms exhibit similar performances in terms of computing time needed.
In Table 5, we can see that these similar CPU time recover totally different realities. The label
dominance algorithm enumerates hundreds of thousands of paths when the two other algorithms
enumerate at most a few thousands. We can see in the proportion of paths cut by the dominance
test in the label correcting algorithm the explanation: the lower bound test cut paths much better
than the dominance test in that setting. This enables the generalized A∗ and the label correcting
algorithm to tackle with larger instances than the label dominance algorithm. The similar CPU
time we obtain at the end come from the fact that the preprocessing time needed to compute the
lower bounds is rather long. We also note that, on this non-constrained problem, the key plays
an important role in the performance of the generalized A∗ and the label correcting algorithm.
Indeed, the versions of the algorithms with the key c(xP ) instead of c(xP ⊕ bv ), identified by the
suffix K in Table 5 exhibit much poorer performances. On difficult instances, candidate paths are
an important element of the performance of the generalized A∗ algorithms. Finally, the numerical
results confirm that the choice of the expectation of ξ as key in the bounding algorithm is relevant:
the parameter γ of Equation (11) remains small.
Constrained problem. Figure 7 and Table 6 are the analogues of Figure 6 and Table 5 for Problem
(14). This constrained problem is much more difficult than the non constrained one. The relative
behaviour of the algorithms is quite similar to the one we observe on the usual resource constrained
shortest path problem (12) with k = 10 constraints. The label correcting and the generalized A∗
algorithms behave much better than the label dominance algorithm. Using candidate paths enable
27

Instance
road50

|V |
|A|
50000 138112

Alg.
γ
A∗
1
A∗ CP
1
A∗ K
1
cor.
1
cor. CP 1
cor. K
1
dom.
–
dom. CP –
square200 40002 120200
A∗
1
∗
A CP
1
A∗ K
1
cor.
1
cor. CP 1
cor. K
1
dom.
–
dom. CP –
long50
12802 38416
A∗
1
∗
A CP
1
A∗ K
1
cor.
1
cor. CP 1
cor. K
1
dom.
–
dom. CP –
wide100
25602 78400
A∗
1
∗
A CP
1
A∗ K
1
cor.
1
cor. CP 1
cor. K
1
dom.
–
dom. CP –

Preproc.
74%
65%+34%
80%
91%
64%+35%
25%
–
– +14%
82%
60%+39%
82%
92%
60%+39%
22%
–
– +5%
36%
52%+35%
57%
34%
52%+35%
6%
–
– +4%
98%
61%+38%
76%
98%
61%+38%
39%
–
– +28%

Ext.
Cut.
Dom.
ℓ
Gap
CPU (s)
5044
0
–
–
∞ 3.15e+01
281
610
– 468
opt 3.27e+01
7330
0
–
–
∞ 2.91e+01
1648
3135
4% 468
opt 2.63e+01
281
610
0% 468
opt 3.49e+01
147208 97207 100%
–
∞ 9.69e+01
147208 97207
–
–
∞ 1.74e+01
147208 97207
– 469 140.3% 7.84e+01
2801
5800
– 261
opt 1.22e+01
194
576
– 261
opt 1.54e+01
4903
0
–
–
∞ 1.27e+01
1162
2488
1% 261
opt 1.14e+01
194
576
0% 261
opt 1.53e+01
147641 97640 100%
–
∞ 4.64e+01
147641 97640
–
–
∞ 1.09e+01
423004 280038
– 261
opt 1.18e+02
4995
0
–
–
∞ 2.07e+01
591
1172
– 1014
opt 1.28e+01
4995
0
–
–
∞ 1.32e+01
5427
288 100%
–
∞ 2.18e+01
591
1172
0% 1014
opt 1.29e+01
148493 98492 100%
–
∞ 1.34e+02
148493 98492
–
–
∞ 2.28e+01
148493 98492
– 1018 112.6% 1.21e+02
20
1638
–
20
opt 2.48e+00
1
1591
–
20
opt 3.60e+00
4203
0
–
–
∞ 3.17e+00
20
1638
0%
20
opt 2.46e+00
1
1591
0%
20
opt 3.59e+00
46003 27868 100%
–
∞ 6.26e+00
46003 27868
–
–
∞ 1.56e+00
46003 27868
–
20 33.3% 5.22e+00

Table 5. Stochastic shortest path problem with generic distributions and β = 0.05.

to speed-up the algorithms on difficult instances. The key in the enumeration algorithm is less
important than for the non-constrained problem.
8. What to do on difficult problems
We have noted in Section 6 that the performance of the algorithms decreases with the number
of constraints. To illustrate what happens, we plot on Figure 8 the resources x = (w0 , w1 ) ∈ R2 of
the usual resource constrained shortest path problem (12) with k = 1 constraint. Each symbol ×
corresponds to the the resource xP of a path P . Figure 8.(a) illustrates why the performance of
the dominance test decreases with the dimension. All it takes is one coordinate such that wPi  wP̃i
for P not to dominate P̃ . When the number of coordinates increases, it becomes rare that a
path dominates another, and the dominance test does not enable to cut well paths. Figure 8.(b)
illustrates the fact thatVthe lower bound bv on the resource of all the v-d paths P must be nongreater than the meet xP of the resources of the v-d paths,
V illustrated by the red diamond on
the figure. When the dimension increases, the gap between xP and the resources xP tends to
28

square

road
103

A∗
cor.
dom.

102

CPU time (s)

CPU time (s)

103

101
100
10−1
10−2

102

103

104

A∗
cor.
dom.

102
101
100
10−1
10−2

105

102

103
104
Vertices

Vertices
long

wide
103

A∗
cor.
dom.

102

CPU time (s)

CPU time (s)

103

101
100
10−1
10−2

102

103
104
Vertices

105

A∗
cor.
dom.

102
101
100
10−1
10−2

105

102

103
104
Vertices

105

Figure 7. Stochastic Resource Constrained Shortest Path Problem with generic distributions
increase, and the quality of the bounds bv decreases. However, a bound bv can still be computed,
which explains the reasonably good performance of the algorithms using the lower bound test (Low)
when the dimension increases.
w1

w1

(a)
xP
×

×

×

(b)
×

×xP̃

V•
xP

w0

×

×
w0

Figure 8. Dominance and lower bound tests when dimension increases. On Figure
(b), each symbol × corresponds to the resource of a v-d path.
Suppose that for each vertex v, we have a set Bv of resources in R such that, for each v-d path
P , there exists a resource b ∈ Bv satisfying b 6 xP . We can then replace the lower bound test by
the following clustered lower bounds test.
B
(Clu) Is there a bound b in Bv such that ρ(xP ⊕ b) = 0 and c(xP ⊕ b) ≤ cU
od ?
29

Instance
road20

|V |
20000

|A|
55180

Alg.
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP
square200 40002 120200
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP
long20
5122 15376
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP
wide100
25602 78400
A∗
∗
A CP
A∗ K
cor.
cor. CP
cor. K
dom.
dom. CP

γ
1.4
1.4
1.4
1.4
1.4
1.4
–
–
1.5
1.5
1.5
1.5
1.5
1.5
–
–
1.5
1.5
1.5
1.5
1.5
1.5
–
–
1
1
1
1
1
1
–
–

Preproc.
11%
9%+3%
3%
88%
66%+24%
76%
–
– +7%
20%
15%+7%
19%
40%
31%+14%
40%
–
– +72%
3%
2%+1%
2%
3%
2%+1%
3%
–
– +4%
99%
65%+34%
99%
99%
66%+34%
99%
–
– +71%

Ext.
79249
79207
314771
1544
1502
3538
113537
113537
87822
87822
95221
30302
30301
32288
37691
37691
115855
115855
128449
62412
62412
62800
140783
140783
55
47
68
45
37
54
8631
8631

Cut.
Dom.
ℓ
Gap
CPU (s)
166646
– 240
opt 7.21e+01
166575
– 240
opt 8.54e+01
669005
–
–
∞ 2.87e+02
2647
9% 240
opt 9.07e+00
2576
9% 240
opt 1.15e+01
6105
9% 240
opt 1.07e+01
63535
–
–
∞ 1.26e+01
63535
–
–
∞ 3.28e+01
165841
–
–
∞ 7.10e+01
165841
– 255 51.3% 8.91e+01
180639
–
–
∞ 7.65e+01
50233
11% 255
opt 3.42e+01
50234
11% 255
opt 4.31e+01
53435
11% 255
opt 3.65e+01
21859
–
–
∞ 1.92e+00
21859
–
–
∞ 6.17e+00
221722
–
–
∞ 8.42e+01
221722
– 395 67.0% 1.04e+02
246911
–
–
∞ 9.43e+01
99413
12%
–
∞ 7.27e+01
99413
12% 395 7.4% 9.00e+01
99675
13%
–
∞ 7.49e+01
90781
–
–
∞ 1.80e+01
90781
–
–
∞ 1.79e+01
1705
– 21
opt 1.99e+00
1692
– 21
opt 2.83e+00
1731
– 21
opt 2.02e+00
1683
0%
21
opt 1.97e+00
1670
0%
21
opt 2.83e+00
1701
0%
21
opt 2.01e+00
2939
–
–
∞ 3.34e-01
2939
–
–
∞ 1.18e+00

Table 6. Stochastic resource constrained shortest path problem with generic distributions.

The idea behind this new test is illustrated on Figure 9.(a), where the lower bounds b ∈ Bv are
represented by blue circles. If we partition the v-d paths Vinto clusters of paths with “similar”
resources, these lower bounds b ∈ Bv are much tighter than xP , and thus enable to discard more
paths. We can also replace the key c(xP ⊕ bv ) by min{c(xP ⊕ b)|b ∈ Bv , ρ(xP ⊕ b) = 0}, which is a
better approximation of the minimum cost of a feasible o-d path starting by P .
Figure 9.(b) illustrates another idea to improve the quality of the bounds. Consider the usual
resource constrained shortest path problem with k constraints of Equation (12). An o-v path P can
0 ≤ cU B − w 0 .
be a subpath of an optimal path only if there exists a feasible v-d path Q such that wQ
od
P
Thus, instead of a bound bv on the resources of all the v-d paths, we can use a bound b′v on the
0 ≤ cU B − w 0 . If cU B is not too large, we can expect the
resource of the v-d paths Q such that wQ
od
P
od
set of such paths Q to be much smaller than the complete set of v-d paths, and thus b′v to be larger
than bv .
We now formalize this idea. We assume to have a weight
morphism ω from (R, ⊕, 6) to (R, +, ≤)
P
ω(x
such that there is no cycle C of negative weight
a ). As ω is a morphism, we have
a∈C
P
ω(xP ) = a∈P ω(xa ). Moreover, we assume to have a scalar ω U B such that ω(xP ) > ω U B implies
30

w1

×
◦×

w1

(a)

b2•

◦
×
V⋄
xP

×
b1• ×

×
◦ ×

(b)
×

b3•
ω1 ω2

w0

×

×
ω3

w0

Figure 9. Lower bounds on clusters of paths (a) and conditional lower bounds (b).

that P is not an optimal feasible path. Note that such a scalar always exists if the set of optimal
paths is finite, which is always the case if we consider only elementary paths. For each vertex v in
V , let nv be an integer, and ωv1 < ωv2 < . . . < ωvnv < ω U B 6 ωvnv +1 = +∞ be a sequence of real
numbers such that ωv1 is the minimum weight of a v-d path, which is well defined because there is no
cycle of negative weight. Finally, for each vertex v, we suppose to have a set of bounds b1v , . . . , bnv v
such that biv is a lower bound on the resource of any v-d path P such that ω(xP ) < ωvi+1 . We can
now define the conditional lower bounds test.
B
Do we have ρ(xP ⊕ biv ) = 0 and c(xP ⊕ biv ) ≤ cU
od , where i is the minimum index such
(Con)
i+1
U
B
that ωv ≥ ω
− ω(xP )?
We can also replace the key c(xP ⊕bv ) by c(xP ⊕biv ), where i is defined as in the test. When the cost
B
U B . This is for instances the case of the
function c is a morphism, we can use it as ω, and cU
od as ω
usual resource constrained shortest path problem (12), and of the stochastic resource constrained
shortest path problem (14). An alternative choice for Problem (14) would be ω : (w, ξ) 7→ E(ξ).
We can adapt the proofs of Theorems 1 and 5 to show that they remain true if we replace the
lower bound test (Low) by the clustered lower bounds test (Clu) or the conditional lower bounds
test (Con) [71]. After a brief overview of the technique that enables to build the sets of bounds
required by both tests, we detail the relative advantages of each of them, and we conclude with
numerical experiments showing the gain they enable.
Remark 6. Provided that we are able to build the bounds biv , any isotone function ω : R → R can
be used in the conditional lower bounds test. However, our technique to build the bounds biv works
only when ω is a morphism.
8.1. Computing the sets of bounds. It is not required to define new bounding algorithms to
compute the lower bounds set of the clustered lower bound test. Our strategy to compute the lower
bounds Bv is to blow-up the digraph D = (V, A) in a much larger digraph D = (V, A), and to use
the algorithm of Section 5 in D. Digraphs V and V are respectively illustrated on the left and on
the right part of Figure 10. We provide here the properties of D that enable to retrieve the bounds
Bv from the blown-up digraph D. These properties are enforced when D is built. The procedures
we provide in [71] to build such a graph D are not difficult but technical: we therefore refer the
interested reader to the Chapter 6 of [71].
We denote ϑ the vertices of D. We have a surjective mapping ϕ that associates to each vertex
ϑ of V a vertex ϕ(ϑ) in V . We denote ϕ−1 (v) the set of vertices ϑ such that ϕ(ϑ) = v. There are
typically many vertices in ϕ−1 (v), the only exception being the set ϕ−1 (d) of vertices corresponding
to the destination d, which is the singleton {ϑd }. We deduce from ϕ a mapping θ from the set of
paths π in D ending in ϑd to the set of paths in P in d. We define the resources of the arcs in A
is such a way that the resource xπ of a path π is equal to the resource xP of the corresponding
path P = θ(π) in D. The most important assumption is the following one: the mapping θ induces
31

a bijection between the ϕ−1 (v)-ϕ−1 (d) paths in D and the v-d paths in D, where a ϕ−1 (v)-ϕ−1 (d)
path is a path between a vertex ϑ ∈ ϕ−1 (v) and ϑd .

D

ϑ1
ϑ2

vertices
π

ϑ3

ϕ−1 (v)

ϑd

D

ϕ
θ
paths
θ −1

v

θ(π)
d

Figure 10. Blown-up graph. ϕ is a surjective mapping between vertices, and θ a
bijective mapping between paths.
As we illustrate on Figure 10, under these assumptions, the bijection θ provides a natural partition the set v-d paths P into clusters. The cluster of a v-d path P is given by the origin vertex
ϑ of the corresponding path π = θ −1 (P ). To illustrate this partition in terms of lower bounds on
resources, the red diamond on Figure 9.(a) is the a lower bound on the resource all the v-d path
in D, and the blue circles are the lower bounds on the resource of all the ϑ-ϑd paths in D for each
ϑ in φ−1 (d). Practically, in the same way we use the algorithm of Equation (7) to compute the
∗
lower bounds bℓv or b†v on the resource of all the v-d paths in D, we can now use these algorithm
∗
in digraph D to obtain the lower bounds bℓϑ or b†ϑ on the resources of all the ϑ-ϑd paths. As θ
induces a bijection, for each v-d path P , there is a vertex ϑ ∈ ϕ−1 (v) and a ϑ-ϑd path π such that
∗
∗
b†ϑ 6 bℓϑ 6 xπ = xP . We can therefore use {bℓϑ , ϑ ∈ ϕ−1 (v)} or {b†ϑ , ϑ ∈ ϕ−1 (v)} as lower bounds
set Bv in the clustered lower bounds test.
The bounds of the conditional graph can also be computed using a similar “blown-up graph”
approach. In both case, the procedure building graph D takes as input the maximum cardinal κ of
ϕ−1 (v). Parameter κ corresponds to the maximum number of bounds we obtain for a given vertex
v: |Bv | ≤ κ for the clustered lower bounds test, and the maximum nv ≤ κ for the conditional lower
bounds test. See Chapter 6 of [71] for more details.
8.2. When to use the clustered and the conditional lower bound test. The clustered and
conditional lower bounds tests (Clu) and (Con) are stronger versions of the lower bound test (Low)
which require a longer preprocessing but enable to reduce the number of paths enumerated by
the generalized A∗ and label correcting algorithms. They are therefore not interesting on easy
instances, where even the preprocessing for the lower bound test (Low) is a waste of time. On
the contrary, they are interesting on difficult instances where the time spent in the enumeration is
much larger than the one spent in the preprocessing, and even more interesting on very difficult
instances that we cannot solve to optimality because the list L becomes to large.
Increasing the size of D increases the preprocessing time and reduces the number of paths enumerated by and the time spent in the enumeration algorithms. On instances that cannot be solved
to optimality, we advise to use the largest graph D that can be stored in the memory. On instances
that can be solved to optimality, we advise to choose the size of D in order to spend about the
same time in the preprocessing and in the enumeration algorithm.
Clustered and conditional lower bound tests have both their pros and cons. The clustered lower
bound test can in theory be applied to any instance of the Monoid Resource Constrained
Shortest Path Problem. Indeed, the procedures that build the digraph D only requires a
similarity measure that enable to compare two resources x and x̃. Besides, as we can see on Figure
9, the bounds it produces tend to be tighter than those produced by the conditional lower bound
32

test. Its first drawback comes from the fact that, each time a test (Clu) is performed at a vertex
v, up to |Bv | operations ⊕, c(·), and ρ(·) must be performed, which can be time consuming. The
other drawbacks are linked to the procedure we use to build the bounds Bv . Indeed, this procedure
calls a clustering subroutine for each vertex v of the digraph D [71], which is time consuming when
D is large. Finally, our procedure that builds D has been thought for acyclic digraphs. We have
extended it to digraphs with cycles, but the quality of the bounds returned is mitigated.
Due to these properties, the clustered lower bound approach works particularly well in the context
of column generation if the pricing subproblem is a resource constrained shortest path problem on
an acyclic digraph. This is often the case when the problem solved by column generation is a
vehicle or a crew scheduling problem. Along a column generation scheme, many instances of the
same resource constrained shortest path problem must be solved successively, the only difference
between two instances being the reduced costs. Hence, we can compute once and for all the
digraph D in a (time-consuming) pre-processing, and then use it to compute the lower bound sets
Bv and speed-up the resolutions of all the instances of the pricing subproblem. We have applied
this technique to a column generation approach to the airline crew pairing problem. As it can be
seen in Table 9.2 of [71], the use of the clustered lower bound test within a generalized A∗ or a
label correcting algorithm enables to divide by three the total time spent in the column generation
scheme on some industrial instances.
Compared to the clustered lower bound test, the main advantage of the conditional lower bound
test is that operations ⊕, c(·), and ρ(·) need to be performed only once when a test (Con) is
performed. Its main drawback is that a morphism ω is required, which reduces the number of
problems to which it can be applied.
The performance of the test depends on the choice of ω. If the set of o-d paths P such that
ω(xP ) ≤ ω U B is much smaller than the set of o-d paths, then the conditional lower bound test
tends to exhibit good performances.
8.3. Numerical results. We now test the performance of the clustered and conditional lower
bounds tests on difficult instances of the usual resource constrained shortest path problem (12)
with k = 10. We use the cost as morphism ω for the conditional lower bounds test: ω(x) = w0
where x = (w0 , . . . , wk ). When we build the sets of bounds Bv or the conditional bounds biv , the
number of bounds we can take for a given vertex v is practically limited by the memory available. We
therefore consider the instances long5 and square50 of Table 4 because they are difficult instances
of reasonable size, which enables to build large sets of bounds Bv . The numerical experiments have
been performed on the same computer as those of Section 6. Both the generalized A∗ and the
label correcting algorithms have been tested with the new tests (Clu) and (Con). We have used
algorithms with candidate paths, as these versions are the most efficient ones on difficult instances.
In Section 6, in order to have comparable results between algorithms and instances, we have set
the maximum size of L to 1e+05 because it is the maximum size that could fit in the memory for
all the instances and all the algorithms. Here, we want to obtain the best possible results with each
algorithm. For a list L of identical size, the label correcting algorithm requires more memory than
then generalized A∗ as it needs to store the lists of non-dominated paths Mv . We have therefore
set 2e+05 as the of maximum size of L for the label correcting, and 2e+06 for the generalized A∗ .
The numerical results are available in Table 7. The column “Test” provides the test used, which
can be the lower bound test (Low), the clustered lower bound test (Clu), or the conditional lower
bounds test (Con). When the clustered or the conditional lower bound test is used, we provide the
|V|
ratio |V
| of the number of vertices in the blown-up digraph D and in the digraph D: it is is the
average number of bounds used per vertex v. The percentage of time spent in the preprocessing
includes the construction of the digraph D, the bounding algorithm in D, and the computation of
candidate paths. The other columns are identical to those of Table 4.
33

Alg.
Test
|V|/|V |
cor. CP (Low)
–
(Con) 1.8e+01
1.7e+02
9.1e+02
(Clu) 1.7e+01
8.2e+01
3.8e+02
A∗ CP (Low)
–
(Con) 1.8e+01
1.7e+02
9.1e+02
(Clu) 1.7e+01
8.2e+01
3.8e+02
square50 cor. CP (Low)
–
(Con) 9.7e+00
8.6e+01
4.7e+02
(Clu) 1.6e+01
7.4e+01
3.3e+02
A∗ CP (Low)
–
(Con) 9.7e+00
8.6e+01
4.7e+02
(Clu) 1.6e+01
7.4e+01
3.3e+02
Instance
long5

γ Prep.
Ext.
Cut.
Dom. ℓ
Gap
CPU (s)
2.2
0%
230537
130784
50% 81 60.8% 7.61e+00
2.2
6%
167283
50956
82% 81 60.0% 4.91e+00
2.6
41%
137348
24950 100% 81 32.2% 8.14e+00
1
73%
133873
22655 100% 81 29.2% 2.05e+01
1.7
17%
190559
88647
52% 81 55.9% 1.17e+01
1.8
22%
187568
83110
55% 81 53.4% 2.54e+01
1.8
35%
190368
85960
55% 81 49.5% 9.85e+01
2.2
0% 1424776
849563
– 81 64.4% 9.49e+00
2.2
2% 1011810
23631
– 81 62.3% 1.41e+01
2.6
16% 1000137
285
– 81 31.0% 2.32e+01
1
37% 1000309
629
– 81 27.8% 4.07e+01
1.8
6% 1276845
553700
– 81 57.9% 3.45e+01
1.8
5% 1218043
436097
– 81 54.5% 1.09e+02
1.8
7% 1165819
331651
– 81 51.7% 4.35e+02
2.2
0%
350553
273819
42% 52 44.9% 1.15e+01
2.2
6%
217971
100571
67% 52 41.1% 5.61e+00
2.5
41%
158700
57177
53% 52 15.1% 8.85e+00
1
73%
165658
71086
42% 52 12.3% 2.24e+01
1.7
23%
299643
207262
46% 52 36.5% 1.64e+01
1.9
30%
285296
203651
41% 52 35.4% 3.32e+01
1.9
47%
307779
231886
40% 52 31.6% 1.22e+02
2.2
0% 9954128 11908301
– 52 44.9% 5.93e+01
2.2
1% 4627965 1255974
– 52 40.2% 5.96e+01
2.5
3% 6772542 5545130
– 52 7.5% 1.11e+02
1
3% 30098520 60197088
– 52
opt 4.73e+02
1.7
2% 8425382 8850810
– 52 34.9% 1.61e+02
1.8
2% 8146176 8292398
– 52 33.8% 4.73e+02
1.9
3% 9374152 10748349
– 52 27.9% 1.71e+03

Table 7. Results of clustered and conditional lower bounds tests Problem (12) with
k = 10 constraints. List L maximum size is 2e+05 for label correcting algorithm
and 2e+06 for A∗ algorithm

As most algorithms do not solve the instance to optimality, the interesting statistic to evaluate
the performance of the algorithm is the gap. With each algorithm, both the clustered lower bound
test and the conditional lower bound test enable to reduce the gap. Indeed, they enable to divide
by two the gap on the instance long5, and to solve the instance square50 to optimality. Besides, the
|V|
larger the number of bounds per vertex, i.e., the larger the ratio |V
| , the smaller is the gap at the
|V|
end. Larger |V
| means longer preprocessing. We have been limited by the memory of the computer
on the choice of the size of D. As, on each algorithm launched, either the instance is not solved to
optimality, or most of the time is spent in the enumeration algorithm, better performances would
be obtained with larger digraphs D.
On our two instances, the conditional lower bounds test performs better than the clustered lower
bounds test. This is not surprising because instances long5 and the square50 have cycle, and, as
we already mentioned, our procedure that builds the lower bounds sets Bv does not work well on
digraphs with cycles.
We finish with two statistics that can look surprising. First, the ratio γ of Equation (11) is
equal to 1 on very large graphs D for the conditional lower bound test. This comes from the fact
that, due to our building procedure, these graphs are acyclic. Second, on the long5 instance, the
34

percentage of paths cut by the dominance test increases with the size of D. This is due to the fact
that, as the key used is a good approximation, only the promising paths are considered, and thus
few paths are cut by the lower bound test at the beginning of the algorithm. The lower bound test
enable to cut paths at the end of the algorithm. As the long5 instance is difficult, the algorithm is
stopped early, and few paths are cut by this test. On the contrary, on the easier square50 instance,
the algorithm is stopped later, and the proportion of paths cut by the dominance test decreases
when the size of D increases.
Acknowledgments
I greatly thank my PhD advisor Frédéric Meunier for his numerous and deep remarks on the
mathematics and the way to write this article.
References
9th

[1] 2006.
DIMACS
Implementation
Challenge,
Shortest
Paths.
http://www.dis.uniroma1.it/challenge9/. Accessed: 2016-12-07.
[2] Adulyasak, Yossiri, & Jaillet, Patrick. 2015. Models and algorithms for stochastic and robust
vehicle routing with deadlines. Transportation Science, 50(2), 608–626.
[3] Aho, Alfred V, & Hopcroft, John E. 1974. The design and analysis of computer algorithms.
Pearson Education India.
[4] Amato, Gianluca, Scozzari, Francesca, Seidl, Helmut, Apinis, Kalmer, & Vojdani, Vesal. 2016.
Efficiently intertwining widening and narrowing. Science of Computer Programming, 120, 1–24.
[5] Apinis, Kalmer, Seidl, Helmut, & Vojdani, Vesal. 2013. How to combine widening and narrowing
for non-monotonic systems of equations. ACM SIGPLAN Notices, 48(6), 377–386.
[6] Artzner, Philippe, Delbaen, Freddy, Eber, Jean-Marc, & Heath, David. 1999. Coherent measures
of risk. Mathematical finance, 9(3), 203–228.
[7] Backhouse, Roland C, & Carré, Bernard A. 1975. Regular algebra applied to path-finding
problems. IMA Journal of Applied Mathematics, 15(2), 161–186.
[8] Bast, Hannah, Delling, Daniel, Goldberg, Andrew, Müller-Hannemann, Matthias, Pajor,
Thomas, Sanders, Peter, Wagner, Dorothea, & Werneck, Renato. 2014. Route Planning in Transportation Networks.
[9] Beasley, JE, & Christofides, Nicos. 1989. An algorithm for the resource constrained shortest
path problem. Networks, 19(4), 379–394.
[10] Bellman, Richard. 1958. On a routing problem. Quarterly of applied mathematics, 87–90.
[11] Bertsimas, Dimitris J. 1992. A vehicle routing problem with stochastic demand. Operations
Research, 40(3), 574–585.
[12] Bertsimas, Dimitris J, & Simchi-Levi, David. 1996. A new generation of vehicle routing research: robust algorithms, addressing uncertainty. Operations Research, 44(2), 286–304.
[13] Blyth, Thomas Scott. 2005. Lattices and ordered algebraic structures. Vol. 1. Springer.
[14] Bogerding, Mark. 2013. KissFFT library.
[15] Borndörfer, Ralf, Grötschel, Martin, & Löbel, Andreas. 2001. Scheduling duties by adaptive
column generation.
[16] Carlyle, W Matthew, Royset, Johannes O, & Kevin Wood, R. 2008. Lagrangian relaxation
and enumeration for solving constrained shortest-path problems. Networks, 52(4), 256–270.
[17] Carré, Bernard A. 1971. An algebra for network routing problems. IMA Journal of Applied
Mathematics, 7(3), 273–294.
[18] Chang, Tsung-Sheng, Wan, Yat-wah, & Ooi, Wei Tsang. 2009. A stochastic dynamic traveling
salesman problem with hard time windows. European Journal of Operational Research, 198(3),
748–759.
35

[19] Chen, Anthony, & Ji, Zhaowang. 2005. Path finding under uncertainty. Journal of advanced
transportation, 39(1), 19–37.
[20] Chen, Bi Yu, Lam, William HK, Sumalee, Agachai, Li, Qingquan, Shao, Hu, & Fang, Zhixiang.
2013. Finding reliable shortest paths in road networks under uncertainty. Networks and spatial
economics, 13(2), 123–148.
[21] Cherkassky, Boris V, Goldberg, Andrew V, & Radzik, Tomasz. 1996. Shortest paths algorithms:
Theory and experimental evaluation. Mathematical programming, 73(2), 129–174.
[22] Cousot, Patrick, & Cousot, Radhia. 1977. Abstract interpretation: a unified lattice model
for static analysis of programs by construction or approximation of fixpoints. Pages 238–252
of: Proceedings of the 4th ACM SIGACT-SIGPLAN symposium on Principles of programming
languages. ACM.
[23] Cousot, Patrick, & Cousot, Radhia. 1979. Constructive versions of Tarski’s fixed point theorems. Pacific Journal of Mathematics, 82(1), 43–57.
[24] Davey, Brian A, & Priestley, Hilary A. 2002. Introduction to lattices and order. Cambridge
university press.
[25] De Silva, Amal. 2001. Combining constraint programming and linear programming on an
example of bus driver scheduling. Annals of Operations Research, 108(1-4), 277–291.
[26] Desrochers, Martin, & Soumis, François. 1988. A generalized permanent labeling algorithm
for the shortest path problem with time windows. INFOR Information Systems and Operational
Research.
[27] Dijkstra, Edsger W. 1959. A note on two problems in connexion with graphs. Numerische
mathematik, 1(1), 269–271.
[28] Dumitrescu, Irina, & Boland, Natashia. 2003. Improved preprocessing, labeling and scaling
algorithms for the Weight-Constrained Shortest Path Problem. Networks, 42(3), 135–153.
[29] Eiger, Amir, Mirchandani, Pitu B, & Soroush, Hossein. 1985. Path preferences and optimal
paths in probabilistic networks. Transportation Science, 19(1), 75–84.
[30] Eppstein, David. 1998. Finding the k shortest paths. SIAM Journal on computing, 28(2),
652–673.
[31] Fahle, Torsten, Junker, Ulrich, Karisch, Stefan E, Kohl, Niklas, Sellmann, Meinolf, & Vaaben,
Bo. 2002. Constraint programming based column generation for crew assignment. Journal of
Heuristics, 8(1), 59–81.
[32] Fan, Yueyue, & Nie, Yu. 2006. Optimal routing for maximizing the travel time reliability.
Networks and Spatial Economics, 6(3-4), 333–344.
[33] Feillet, Dominique, Dejax, Pierre, Gendreau, Michel, & Gueguen, Cyrille. 2004. An exact
algorithm for the elementary shortest path problem with resource constraints: Application to
some vehicle routing problems. Networks, 44(3), 216–229.
[34] Fink, Eugene. 1992. A survey of sequential and systolic algorithms for the algebraic path
problem.
[35] Flajolet, Arthur, Blandin, Sébastien, & Jaillet, Patrick. 2014. Robust adaptive routing under
uncertainty. arXiv preprint arXiv:1408.3374.
[36] Ford, Lester R, & Fulkerson, Delbert R. 1956. Maximal flow through a network. Canadian
journal of Mathematics, 8(3), 399–404.
[37] Frank, H. 1969. Shortest paths in probabilistic graphs. Operations Research, 17(4), 583–599.
[38] Fu, Liping. 2001. An adaptive routing algorithm for in-vehicle route guidance systems with
real-time information. Transportation Research Part B: Methodological, 35(8), 749–765.
[39] Fu, Liping, & Rilett, Larry R. 1998. Expected shortest paths in dynamic and stochastic traffic
networks. Transportation Research Part B: Methodological, 32(7), 499–516.
[40] Gondran, Michel, & Minoux, Michel. 2008. Graphs, dioids and semirings: new models and
algorithms. Vol. 41. Springer Science & Business Media.
36

[41] Gounaris, Chrysanthos E, Wiesemann, Wolfram, & Floudas, Christodoulos A. 2013. The
robust capacitated vehicle routing problem under demand uncertainty. Operations Research, 61(3),
677–693.
[42] Gualandi, Stefano, & Malucelli, Federico. 2009. Constraint programming-based column generation. 4OR, 7(2), 113–137.
[43] Hall, Randolph W. 1986. The fastest path through a network with random time-dependent
travel times. Transportation science, 20(3), 182–188.
[44] Handler, Gabriel Y, & Zang, Israel. 1980. A dual algorithm for the constrained shortest path
problem. Networks, 10(4), 293–309.
[45] Hart, Peter E, Nilsson, Nils J, & Raphael, Bertram. 1968. A formal basis for the heuristic
determination of minimum cost paths. Systems Science and Cybernetics, IEEE Transactions on,
4(2), 100–107.
[46] Ioachim, Irina, Gelinas, Sylvie, Soumis, Francois, & Desrosiers, Jacques. 1998. A dynamic
programming algorithm for the shortest path problem with time windows and linear node costs.
Networks, 31(3), 193–204.
[47] Irnich, Stefan. 2008. Resource extension functions: Properties, inversion, and generalization
to segments. OR Spectrum, 30(1), 113–148.
[48] Irnich, Stefan, & Desaulniers, Guy. 2005. Shortest path problems with resource constraints.
Springer.
[49] Irnich, Stefan, & Villeneuve, Daniel. 2006. The shortest-path problem with resource constraints
and k-cycle elimination for k ≥ 3. INFORMS Journal on Computing, 18(3), 391–406.
[50] Jaillet, Patrick. 1988. A priori solution of a traveling salesman problem in which a random
subset of the customers are visited. Operations research, 36(6), 929–936.
[51] Jaillet, Patrick, & Odoni, A. 1988. The probabilistic vehicle routing problem. Vehicle routing:
methods and studies. North Holland, Amsterdam.
[52] Jaillet, Patrick, Qi, Jin, & Sim, Melvyn. 2016. Routing optimization under uncertainty. Operations Research.
[53] Joksch, Hans C. 1966. The shortest route problem with constraints. Journal of Mathematical
analysis and applications, 14(2), 191–197.
[54] Jula, Hossein, Dessouky, Maged, & Ioannou, Petros A. 2006. Truck route planning in nonstationary stochastic networks with time windows at customer locations. IEEE Transactions on
Intelligent Transportation Systems, 7(1), 51–62.
[55] Junker, Ulrich, Karisch, Stefan E, Kohl, Niklas, Vaaben, Bo, Fahle, Torsten, & Sellmann,
Meinolf. 1999. A framework for constraint programming based column generation. Pages 261–274
of: International Conference on Principles and Practice of Constraint Programming. Springer.
[56] Kohl, Niklas, Desrosiers, Jacques, Madsen, Oli BG, Solomon, Marius M, & Soumis, Francois.
1999. 2-path cuts for the vehicle routing problem with time windows. Transportation Science,
33(1), 101–116.
[57] Kosuch, Stefanie, & Lisser, Abdel. 2010. Stochastic shortest path problem with delay excess
penalty. Electronic Notes in Discrete Mathematics, 36, 511–518.
[58] Larsen, Jesper. 1999. Parallelization of the vehicle routing problem with time windows. Ph.D.
thesis, Technical University of DenmarkDanmarks Tekniske Universitet, Department of Informatics and Mathematical ModelingInstitut for Informatik og Matematisk Modellering.
[59] Lehmann, Daniel J. 1977. Algebraic structures for transitive closure. Theoretical Computer
Science, 4(1), 59–76.
[60] Li, Xiangyong, Tian, Peng, & Leung, Stephen CH. 2010. Vehicle routing problems with time
windows and stochastic travel and service times: Models and algorithm. International Journal of
Production Economics, 125(1), 137–145.
37

[61] Loui, Ronald Prescott. 1983. Optimal paths in graphs with stochastic or multidimensional
weights. Communications of the ACM, 26(9), 670–676.
[62] Lozano, Leonardo, & Medaglia, Andrés L. 2013. On an exact method for the constrained
shortest path problem. Computers & Operations Research, 40(1), 378–384.
[63] Mazmanyan, Lilit, Trietsch, Dan, & Baker, KR. 2009. Stochastic traveling salesperson models
with safety time. Tech. rept. Working paper.
[64] Mirchandani, Pitu B. 1976. Shortest distance and reliability of probabilistic networks. Computers & Operations Research, 3(4), 347–355.
[65] Mohri, Mehryar. 2002. Semiring frameworks and algorithms for shortest-distance problems.
Journal of Automata, Languages and Combinatorics, 7(3), 321–350.
[66] Murthy, Ishwar, & Sarkar, Sumit. 1996. A relaxation-based pruning technique for a class of
stochastic shortest path problems. Transportation Science, 30(3), 220–236.
[67] Murthy, Ishwar, & Sarkar, Sumit. 1998. Stochastic shortest path problems with piecewiselinear concave utility functions. Management Science, 44(11-part-2), S125–S136.
[68] Nie, Yu, & Fan, Yueyue. 2006. Arriving-on-time problem: discrete algorithm that ensures
convergence. Transportation Research Record: Journal of the Transportation Research Board,
1964(1), 193–200.
[69] Nikolova, Evdokia. 2010. High-performance heuristics for optimization in stochastic traffic
engineering problems. Pages 352–360 of: Large-Scale Scientific Computing. Springer.
[70] Nikolova, Evdokia, Kelner, Jonathan A, Brand, Matthew, & Mitzenmacher, Michael. 2006.
Stochastic shortest paths via quasi-convex maximization. Pages 552–563 of: Algorithms–ESA
2006. Springer.
[71] Parmentier, Axel. 2016. Algorithms for shortest path and airline problems. Ph.D. thesis, École
des Ponts Paristech, Université Paris Est.
[72] Powell, Warren B, & Chen, Zhi-Long. 1998. A generalized threshold algorithm for the shortest path problem with time windows. DIMACS Series in Discrete Mathematics and Theoretical
Computer Science, 40, 303–318.
[73] Righini, Giovanni, & Salani, Matteo. 2009. Decremental state space relaxation strategies and
initialization heuristics for solving the orienteering problem with time windows with dynamic
programming. Computers & Operations Research, 36(4), 1191–1203.
[74] Rockafellar, R Tyrrell, & Uryasev, Stanislav. 2000. Optimization of conditional value-at-risk.
Journal of risk, 2, 21–42.
[75] Rousseau, Louis-Martin, Gendreau, Michel, Pesant, Gilles, & Focacci, Filippo. 2004. Solving
VRPTWs with constraint programming based column generation. Annals of Operations Research,
130(1-4), 199–216.
[76] Roy, Bernard. 1959. Transitivité et connexité. Comptes Rendus Hebdomadaires Des Seances
De L Academie Des Sciences, 249(2), 216–218.
[77] Russell, RA, & Urban, TL. 2008. Vehicle routing with soft time windows and Erlang travel
times. Journal of the Operational Research Society, 59(9), 1220–1228.
[78] Sabran, Guillaume, Samaranayake, Samitha, & Bayen, Alexandre M. 2014. Precomputation
techniques for the stochastic on-time arrival problem. Pages 138–146 of: ALENEX. SIAM.
[79] Samaranayake, Samitha, Blandin, Sebastien, & Bayen, A. 2012. A tractable class of algorithms for reliable routing in stochastic networks. Transportation Research Part C: Emerging
Technologies, 20(1), 199–217.
[80] Santos, Luis, Coutinho-Rodrigues, João, & Current, John R. 2007. An improved solution algorithm for the constrained shortest path problem. Transportation Research Part B: Methodological,
41(7), 756–771.
[81] Sivakumar, Raj A, & Batta, Rajan. 1994. The variance-constrained shortest path problem.
Transportation Science, 28(4), 309–316.
38

[82] Sungur, Ilgaz, Ordónez, Fernando, & Dessouky, Maged. 2008. A robust optimization approach
for the capacitated vehicle routing problem with demand uncertainty. IIE Transactions, 40(5),
509–523.
[83] Taş, D, Gendreau, M, Dellaert, N, Van Woensel, Tom, & De Kok, AG. 2014. Vehicle routing
with soft time windows and stochastic travel times: A column generation and branch-and-price
solution approach. European Journal of Operational Research, 236(3), 789–799.
[84] Zimmermann, Uwe. 1981. Linear and Combinatorial Optimization in Ordered Algebraic Structures. Elsevier.
A. Parmentier, École Nationale des Ponts et Chaussées, CERMICS, 6-8 avenue Blaise Pascal, Cité
Descartes, 77455 Marne-la-Vallée, Cedex 2, France
E-mail address: axel.parmentier@cermics.enpc.fr

39

