Integer programming formulations for the elementary shortest path

Integer programming formulations for the
elementary shortest path problem
Leonardo Taccari
Dipartimento di Elettronica, Informazione e Bioingegneria,
Politecnico di Milano, Italy
[email protected]
Abstract
Given a directed graph G = (V, A) with arbitrary arc costs, the Elementary Shortest Path Problem (ESPP) consists of finding a minimumcost path between two nodes s and t such that each node of G is visited
at most once. If the costs induce negative cycles on G, the problem is
N P-hard. In this paper, several integer programming formulations for
the ESPP are compared. We present analytical results based on a polyhedral study of the formulations, and computational experiments where we
compare their linear programming relaxation bounds and their behavior
within a branch-and-cut framework. The computational results show that
a formulation with dynamically generated cutset inequalities is the most
effective.
1
Introduction
Given a directed graph G = (V, A) and arc costs cij for each (i, j) ∈ A, the
shortest path problem consists of finding a minimum-cost path between two
nodes s and t.
Often, an implicit assumption is that such path has to be elementary. A
path is defined to be elementary if it does not visit any node more than once,
i.e., if it does not contain subtours. When the costs cij induce no negative
cycles on G, the problem can be solved efficiently via ad-hoc polynomial time
algorithms like Bellman-Ford’s or Dijikstra’s (if cij ≥ 0). However, if negative
cycles do arise, subtours must be explicitly prevented, leading to the so-called
Elementary Shortest Path Problem (ESPP).
The ESPP is clearly N P-hard due to a simple reduction from the Hamiltonian path problem. The variant with only negative arc costs has been vastly
discussed in the literature (regarded as a maximization problem), and is usually referred to as the Longest Path Problem (LPP). The authors in [BHK03]
prove that the LPP is hard to approximate on directed graphs within a n1−
for any , unless P = N P. For undirected graphs, approximation algorithms
are described in, e.g., [KMR97] and [Gab07].
1
Related to the ESPP are variants of the Travelling Salesman Problem (TSP),
such as the Prize Collecting TSP [Bal89] and the Orienteering Problem [VSO11],
that involve determining optimal tours where each node is visited at most once.
While an interesting problem in its own right, the ESPP may also arise in the
pricing subproblems of branch-and-price algorithms [ACT14]. More often, the
pricing phase in Vehicle Routing Problems (VRP) involve resource-constrained
variants of the elementary shortest path problem (ESPPRC). The problems
are usually solved with fast dynamic programming-based labelling algorithms,
e.g., see [RS06], [FDGG04] and [BDD06]. It is possible to adapt this kind
of approach to the unconstrained ESPP by considering an artificial resource
for each node, and imposing that less that one unit of each resource is used, as
already proposed in [BC89]. However, this is a very weak constraint, so that the
typical approaches for the ESPPRC do not seem to carry over effectively to the
ESPP [DI14]. Even when resource constraints are present, labelling algorithms
are sometimes not applicable or inefficient, and in those cases branch-and-cut
algorithms might be a better choice [JPS08]. In the context of a branch-andprice algorithm, where the ESPP is solved repeatedly in the pricing phase,
another desirable feature of an integer programming approach is its flexibility,
that allows to easily incorporate branching decisions. These reasons motivate
the study of integer programming techniques for the ESPP.
In the literature, not much previous work has appeared on integer programming approaches tailored specifically for the ESPP. Ibrahim et al. in [IMM09]
provide computational results on the linear programming (LP) bounds of a
flow-based extended formulation, but give no details on its behavior in a branchand-bound algorithm. In [HMM13], another extended formulation is proposed.
Drexl and Irnich in [DI14] describe a branch-and-cut approach and compare its
efficiency with the extended formulation in [IMM09], while in [Dre13] the efficient separation of subtour elimination constraints for the ESPP is discussed.
In this article we present a comparison between several formulations for the
ESPP, which are described in detail in Section 2. We include integer programming formulations with exponentially many subtour elimination constraints, and
mixed-integer programming extended formulations with a polynomial number
of variables and constraints. In Section 3 we provide some analytical results,
including a proof of equivalence between the polyhedra described by the two
strongest formulations. Section 4 reports computational experiments where we
compare the LP relaxation bounds and branch-and-cut results. Finally, in Section 5, we report conclusions and discuss further research topics.
2
Integer programming formulations
Let us consider a directed graph G = (V, A) with the set of nodes V and the set
of arcs A. Let n and m be the cardinality of, respectively, V and A. A path is
a sequence of nodes v1 , . . . , vk , and is said to be elementary if no node appears
in the path more than once. A cycle or subtour is a path with v1 = vk . With
δ + (i) and δ − (i) we denote the set of outgoing and incoming arcs with respect
2
to node i, with δ + (S) and δ − (S) the arcs leaving/entering the set S ⊆ V ,
and with A(S) the set of arcsPwith both ends in S ⊆ V . Let us also define
Vi := V \ {i}, and x(S) :=
i∈S xi . In all the formulations, it is assumed
w.l.o.g. that |δ − (s)| = |δ + (t)| = 0.
A standard integer programming formulation to determine the shortest path
from node s to node t is the following:
X
min
cij xij
(1)
(i,j)∈A


1
xij −
xji = −1


(i,j)∈δ + (i) (j,i)∈δ − (i)
0
X
xij ≤ 1
X
if i = s
if i = t
else
X
i∈V
(2)
i∈V
(3)
(i, j) ∈ A,
(4)
(i,j)∈δ + (i)
xij ∈ {0, 1}
where cij ∈ R are the arc costs, and xij are binary arc variables that take
value 1 if the arc (i, j) belongs to the optimal path. Constraints (2) are balance
equations imposing there is a path connecting s and t, while Constraints (3)
ensure that the outgoing degree of each node is at most one. When the costs cij
induce negative cycles on G, i.e., there is a subtour such that the total cost of
its arcs is negative, this system of inequalities is not sufficient to guarantee the
elementarity of the path. Thus, additional constraints (and, possibly, variables)
are necessary to prevent subtours ensuring that no node is visited more than
once.
Notice that the crucial difference with respect to problems where the path
has to be Hamiltonian the absence, in the ESPP, of the strict degree constraints:
X
X
xij =
xji = 1
i ∈ V.
(i,j)∈δ + (i)
(j,i)∈δ − (i)
We now describe different set of constraints and variables that can be added
to the Formulation (1)–(4) so as to obtain a valid integer programming formulation for the ESPP.
2.1
Dantzig-Fulkerson-Johnson (DFJ)
For the TSP, success has been achieved with strong formulations with exponentially many constraints. It is possible to write a formulation for the
ESPP based on the classical Dantzig-Fulkerson-Johnson subtour elimination
constraints [DFJ54], adding to the basic formulation (1)–(4) the following inequalities:
X
xij ≤ |S| − 1
S ⊆ Vst , |S| ≥ 2.
(5)
(i,j)∈A(S)
3
In each subset S, subtours are prevented ensuring that the number of arcs in S
which are selected is smaller than the number of nodes in S. This formulation
includes O(m) variables and O(2n ) constraints.
Observation For the Asymmetric TSP (ATSP), due to the degree constraints,
the DFJ subtour eliminations constraints can be written in the equivalent form:
X
xij ≥ 1
S ⊂ V, |S| ≥ 2.
(6)
(i,j)∈δ + (S)
For the ESPP, Constraints (6) are valid only for subsets S with s ∈ S and
t∈
/ S. Moreover, they are not equivalent to (5), and they are not sufficient to
prevent subtours.
Example 1. Consider the solution reported in Figure 1, assuming it is a complete graph and that only the arcs with xij = 1 are drawn. The solution does
not violate any inequality (6) for any set S containing s, since x(δ + (S)) = 1
for any such S, although it contains a (disconnected) subtour. On the other
hand, notice that the inequality (6) is not valid, e.g., for the set S 0 , not containing s and containing only nodes that are not part of the s-t path, where
x(δ + (S 0 )) = 0.
S
S0
s
t
Figure 1: An example where Constraints (6) are not sufficient to prevent subtours.
2.2
Generalized cutset inequalities (GCS)
DFJ Constraints (6) can be adapted to the ESPP by replacing the constant
right-hand side with a variable expression. This approach was used for a symmetric version of ESPPRC in [JPS08] and applied to the asymmetric ESPP
in [DI14]. Similar subtour elimination constraints have also been used in branchand-cut algorithms for the VRP (see e.g. [NR02]). We refer to them as generalized cutset inequalities (GCS):
X
X
xij ≥
xkj
k ∈ S ⊆ Vst , |S| ≥ 2.
(7)
(i,j)∈δ + (S)
(k,j)∈δ + (k)
Constraints (7) prevent subtours by ensuring that, for each subset S, the number
of selected arcs leaving S is greater than the number of selected arcs outgoing
4
from any node in S. In an integer solution, this means that the cut induced
by S must containt at least one arc if at least one node in S belongs to the s-t
path, while, if S does not contain any node in the s-t path, the constraint is the
trivial inequality 0 ≥ 0. The number of variables in the formulation is O(m),
while the number of constraints is O(n2n ).
Note that Constraints (7) can be shown to be equivalent to a strengthened
version of (5).
Proposition 2. Constraints (7) are equivalent to:
X
X
xij ≤
(i,j)∈A(S)
X
xij
i∈S\{k} (i,j)∈δ + (i)
∀k ∈ S ⊆ Vst ,
|S| ≥ 2.
(8)
Proof.
For all k ∈ S, x(δ + (k)) ≤ x(δ + (S)) = x(δ + (S)) + x(A(S)) − x(A(S)) =
P
+
i∈S x(δ (i)) − x(A(S)), where the inequality follows from (7).
These inequalities can be interpreted as imposing that the number of selected
arcs in a subset S is strictly smaller than the number of nodes in S that belong
to the s-t path.
2.3
Sequential formulations (MTZ)
To derive an extended formulation à la Miller, Tucker and Zemlin [MTZ60]
(hereafter MTZ) it is enough to introduce, for each node, an auxiliary variable
that can be viewed as the order in which the node is visited and a constraint
for each arc:
tj ≥ ti + 1 + (n − 1)(xij − 1)
(i, j) ∈ A
i 6= s, j 6= t.
(9)
For the ATSP, this formulation is well-known to give poor linear relaxation
bounds. However, it is very compact, as it requires only O(m) additional constraints and O(n) auxiliary variables.
2.4
Reformulation-linearization based (RLT)
From the following nonlinear version of the MTZ formulation:
tj xij = (ti + 1)xij
(i, j) ∈ A
(10)
+
(11)
(s, j) ∈ δ (s),
tj xsj = xsj
5
the authors in [HMM13] use the Sherali-Adams reformulation-linearization technique to obtain the following stronger formulation for the ESPP:
αij = βij + xij
X
X
xsj +
αij −
βji = 0
−
(i,j)∈δ (j)
i6=s
X
αij −
(i,j)∈δ − (j)
(j,i)∈δ + (j)
X
βji = 0
(j,i)∈δ + (j)
xij ≤ αij
αij ≤ (n − 1)xij
(i, j) ∈ A
(s, j) ∈ δ (s),
j 6= t
(13)
(s, j) ∈
/ δ + (s)
j 6= s, j 6= t
(14)
(i, j) ∈ A,
i 6= s, j 6= s
(15)
(i, j) ∈ A, i 6= s
(i, j) ∈ A,
i 6= s, j 6= s
xij ≤ βij
βij ≤ (n − 1)xij
βij ≤ αij
(12)
+
(16)
(17)
(i, j) ∈ A : i 6= s
(18)
(i, j) ∈ A.
(19)
This extended formulation requires O(m) constraints and O(m) auxiliary variables.
2.5
Single-flow formulation (SF)
A formulation which is close in spirit to the single-flow ATSP formulation of
[GG78], can be obtained introducing an auxiliary flow q to be delivered to the
nodes belonging to the s-t path:
qij ≤ (n − 1)xij
X
X
qsj =
zk
(s,j)∈δ + (s)
X
(i,k)∈δ − (k)
X
(i, j) ∈ A
(20)
(21)
k∈V
X
qik −
qkj = zk
k ∈ Vs
(22)
k ∈ Vs
(23)
(i, j) ∈ A.
(24)
(k,j)∈δ + (k)
xik = zk
(i,k)∈δ − (k)
qij ≥ 0
Constraints (20) impose that the auxiliary flow is positive only over the arcs
where xij = 1. The auxiliary flow leaving from the node s has value equal
to the number of nodes that are reached by the s-t path. Constraints (22)
ensure that the balance of the auxiliary flow on each node is equivalent to zk ,
which, according to Constraint (23), is either 1, if node k is in the s-t path, or
0 otherwise.
6
2.6
Multicommodity-flow formulation (MCF)
An extension of the single-flow formulation is obtained by disaggregating the
auxiliary flow into n − 1 unitary flows. Subtours are prevented by enforcing, in
the support graph, one unit of a distinct auxiliary flow from s to each node that
belongs to the s-t path.
k
qij
≤ xij


zk
k
k
qij −
qji = −zk


(i,j)∈δ + (i) (j,i)∈δ − (i)
0
X
xik = zk
X
X
if i = s
if i = k
else
k ∈ Vs ,
(i, j) ∈ A
(25)
i ∈ V,
k ∈ Vs
(26)
k ∈ Vst
(27)
k ∈ Vs
(28)
(i,k)∈δ + (k)
X
xik = zk
(i,k)∈δ − (k)
X
xsj = 1
(29)
xit = 1
(30)
(s,j)∈δ + (s)
X
(i,t)∈δ − (t)
(i, j) ∈ A,
k ∈ Vs .
k
qij
≥ 0, zk ∈ {0, 1}
(31)
The formulation includes O(nm) additional variables and constraints. This
extended formulation is introduced, for the ESPP, in [IMM09], and it is very
similar to the classic multi-commodity flow formulations for the ATSP proposed
by Wong [Won80] and Claus [Cla84].
In Table 1 we summarize the presented formulations. To the best of our
knowledge, formulations MTZ and SF have not been previously considered for
the ESPP.
DFJ
GCS
MTZ
RLT
SF
MCF
vars
m
m
m
m
m
nm
cons
2n
n2n
n
m
m
nm
description
Dantzig-Fulkerson-Johnson (5)
generalized cutsets (7)
Miller-Tucker-Zemlin (9)
reformulation-linearization (12)–(19)
single-flow (20)-(24)
multi-commodity flow (25)-(30)
Table 1: A summary of the considered formulations.
3
Polyhedral results
Let us describe some analytical results for the considered ESPP formulations.
7
Proposition 3. Formulation MCF is stronger than formulation SF.
Proof. Constraints (20)–(22) can be obtained by MCF P
simply aggregating Conk
straints (25)–(26) over k ∈ Vs , and then substituting k∈Vs qij
with qij . The
example in Figure 2 shows that the inclusion is strict.
Proposition 4. Formulation GCS is stronger than formulation DFJ.
Proof. The result follows by considering GCS as stated in (8), whose righthand side is obviously smaller or equal to |S| − 1, right-hand side in (5), and
the inclusion is strict by the example in Figure 2.
t
−20
1
c
−10
s
b
a
1
−10
−10
Figure 2: Example proving strict inclusion for Proposition 3 and 4. With GCS
and MCF, the LP optimal solution is the one with xst = 1 and optimal value
−20. With SF, the optimal solution has value −26, with xst = xca = 14 ,
xsa = xcta = 34 and xab = xbc = 1. With DFJ, the solution has value −40, with
xst = 1 and a disconnected subtour with xab = xbc = xca = 32 .
Showing that formulation MCF is as tight as GCS requires to calculate the
projection of the MCF extended formulation into the space of the x variables.
We will make use of two strong results presented in [PS91], and follow a similar
approach to the equivalence proofs therein.
Theorem 5. The projection of the MCF-polytope onto the x-space is equivalent
to the GCS-polytope.
Proof. Let us rewrite formulation MCF by projecting out the zk variables. Constraints (34) and (35) follow from Constraint (28), while Constraints (27) and
(28). Assume again that the graph G does not contain the arcs δ − (s) and δ + (t).
8
The compact formulation reads:
k
qij
≤ xij
X
X
k
qij
−
k
qji
=0
(j,i)∈δ − (i)
(i,j)∈δ + (i)
X
X
k
qsj
=
(j,k)∈δ − (k)
(33)
k ∈ Vs
(34)
xik
k ∈ Vs
(35)
xik
k ∈ Vst
(36)
(i,k)∈δ − (k)
X
xik =
(i,k)∈δ − (k)
(i,k)∈δ + (k)
X
i ∈ Vs , i 6= k,
k ∈ Vs
xik
X
k
qjk
=
X
(32)
(i,k)∈δ − (k)
(s,j)∈δ + (s)
X
k ∈ Vs ,
(i, j) ∈ A
xsj = 1
(37)
xit = 1
(38)
(s,j)∈δ + (s)
X
(i,t)∈δ − (t)
k
qij
≥ 0, xij ≥ 0
(i, j) ∈ A, k ∈ Vs .
(39)
In order to compare MCF and GCS, we need to project out also the q-variables
of the M CF formulation. Let us define the sets:
X = {x ∈ Rm | x satisfies (36), (37) and (38)},
PGCS = {x ∈ X | x satisfies (7)},
PM CF = {(x, q) ∈ Rmn | (x, q) satisfies (32)–(39)},
P rojx (PM CF ) = {x ∈ X | ∃ q s.t. (x, q) ∈ PM CF },
where PGCS is the GCS-polytope, PM CF is the MCF -polytope and P rojx (PM CF )
is its projection onto the x-space. It is convenient to rewrite Constraints (32)–
(35) in matrix form as follows:
Bx + M q = 0
(40)
− Dx + Iq ≤ 0
(41)
x, q ≥ 0.
(42)
Equation (40) corresponds to (33)–(35), while (41) corresponds to (32). One
can easily notice that, since the n − 1 flows are independent, the system can be
decomposed according to the index k. The matrix M is, in fact, block-diagonal,
and can be decomposed into n − 1 blocks Mk for all k ∈ Vs . Each block Mk
represents the node-arc incidence matrix of the graph G. I is the identity matrix,
that can be decomposed into submatrices Ik of dimension m × m. In a similar
fashion, we also decompose B and D into submatrices Bk and Dk . The rows and
columns of Bk correspond, respectively, to the nodes and arcs of the graph G.
9
A submatrix Bk has zeros everywhere, except for the row corresponding to node
s, with entries of value −1 for each arc in δ − (k), and the row corresponding to
node k, with +1 entries for each arc in δ − (k). Each row of Dk corresponds to a
k
variable qij
and has zeros everywhere, except for a +1 in the column associated
with variable xij .
We now make use of Theorem 2 in [PS91], by which the projection onto the
x-space of the polytope PM CF defined by (40)–(41) can be obtained as:
P rojx (PM CF ) = {x ∈ X | (uB − vD − w)x ≤ 0
∀(u, v, w) ∈ C},
(43)
where C is the cone defined as:
C = {(u, v, w) | uM + vI ≥ 0, v ≥ 0, w ≥ 0}.
The result allows us to carry out the comparison between PGCS and P rojx (PM CF )
simply by finding a system of generators for the cone C.
From the inequalities w ≥ 0 we obtain extreme rays of the form u = 0,
v = 0, w = ei , where ei is the i-th standard basis vector of Rm , that yield the
nonnegativity constraints
xij ≥ 0
∀ (i, j) ∈ A.
(44)
This allows us to restrict our following study to the cone C 0 defined as:
C 0 = {(u, v) | uM + vI ≥ 0, v ≥ 0}.
Exploiting the decomposition of M , we can work on the even smaller cones:
Ck = {(uk , v k ) | uk Mk + v k ≥ 0, v k ≥ 0}.
(45)
Due to (43), once we have the system of generators (uk , v k ) for each cone Ck ,
the constraints in the x-space are obtained by calculating (uk Bk − v k Dk )x ≤ 0
for each k ∈ Vs .
According to Proposition 6 in [PS91], a full system of generators of a cone
Ck defined as in (45), where Mk is a node-arc incidence matrix of a digraph, is
given by:
- a basis of its lineality space, of the form:
uk = ±e, v k = 0,
where e is the all-ones vector, that in our case translate to the trivial equality
0 = 0, and
- the extreme rays, given by all the positive multiples of the vector (uk , v k ) such
10
that:
(i)
uki
k
vij
= 0 ∀i ∈ V,
(
1 ∀ i ∈ S,
(ii)
=
0 otherwise
(
−1 ∀ i ∈ S,
(iii) uki =
0 otherwise
uki
(
1 for one (i, j) ∈ A
=
0 otherwise
(
¯ j ∈ S,
1 ∀ i ∈ S,
k
vij =
0 otherwise
(
¯
1 ∀ i ∈ S, j ∈ S,
k
vij
=
0 otherwise
for any S ⊆ V , where S¯ = V \ S. The extreme rays of the form (i) give rise to
nonnegativity constraints.
From the extreme rays given by (ii) and (iii) we obtain the inequalities:
− x(δ − (S)) + uk x(δ − (k)) − us x(δ − (k)) ≤ 0
−
+
−
− x(δ (S)) − uk x(δ (k)) + us x(δ (k)) ≤ 0
(46)
(47)
where ui = 1 if i ∈ S, and 0 otherwise. For both (46) and (47), we can
distinguish four cases depending on whether s and k are in S, thus whether us , uk
are 0 or 1. If both s and k are in S, or neither of them is, the inequality is implied
by the nonnegativity constraints (44). If only the coefficient with negative sign
is nonzero, the corresponding inequality is, again, redundant. Therefore, the
only meaningful cases are the following:
x(δ − (S)) ≥ x(δ − (k))
+
−
x(δ (S)) ≥ x(δ (k))
∀S ⊆ V, s ∈
/ S, k ∈ S
(48)
∀S ⊆ V, s ∈ S, k ∈
/ S.
(49)
We have established so far that P rojx (PM CF ) is fully described by the nonnegativity constraints and Constraints (48)–(49). This set of inequalities can be
shown to be equivalent to:
x(δ + (S)) ≥ x(δ + (k))
∀S ⊆ Vst , k ∈ S.
(50)
Constraint (48) and (49) are equivalent, due to the fact that x(δ + (S)) =
¯
x(δ (S)).
Let us then consider only (48). For k = t, the inequality is trivially
satisfied by all x ∈ X, thus redundant. For k 6= t and t ∈
/ S, we obtain exactly
the inequalities in (50), since by (36)–(38), we have that x(δ − (k)) = x(δ + (k))
and x(δ − (S)) = x(δ + (S)) for any S containing neither s nor t. If k 6= t and
t ∈ S, it suffices to observe that, since δ + (t) = 0, the inequality x(δ − (S)) ≥
x(δ − (k)) is implied by x(δ − (S \{t})) ≥ x(δ − (k)), which, again, can be rewritten
in the form (50).
Hence, the projection of PM CF onto the x-space is given by
−
P rojx (PM CF ) = {x ∈ X | x satisfies (36)–(38) and (50)},
and it follows that P rojx (PM CF ) = PGCS .
From this result and Proposition 4, it also follows that formulation MCF is
stronger than formulation DFJ.
11
4
Computational comparison
We compare the described formulations with respect to their LP relaxation
bounds and their behavior within an exact branch-and-bound framework.
Three types of instances are considered for the tests. The first pool consists
of instances from the pricing phase of the unsplittable flow problem in [ACT14]
on small-sized networks from the SNDlib [OWPT10], namely, the topologies
abilene, atlanta, france, geant, germany, nobel-us, polska.
The second pool is a set of small to medium-sized random-cost graphs, either
sparse (rnd-s) or dense (rnd-d). The graphs for the instances in rnd-s are
generated by building a connected component including all the nodes, and then
randomly adding arcs until the desired sparsity is reached. The instances in
rnd-d are the dense instances in [Dre13], with random arc costs on a complete
graph.
The third pool (prc) contains the pricing instances in [Dre13]. It consists
of small and medium-size pricing problems from a column generation algorithm
for the asymmetric m-salesmen TSP at the first (f), penultimate (p) and last
(l) pricing iterations. Table 2 summarizes the features of the test instances.
abilene
polska
nobel-us
atlanta
geant
france
germany
rnd-s
rnd-s
rnd-d
prc-f
prc-p
prc-l
n
12
12
14
15
22
25
50
50/100/200
500/1000
25/50/100
27/52/102
27/52/102
27/52/102
m
30
36
42
44
72
90
176
164/660/2654
16634/66601
600/2450/9900
702/2652/10302
702/2652/10302
702/2652/10302
range
[−10000,10000]
[−10000,10000]
[−10000,10000]
[−10000,10000]
[−10000,10000]
[−10000,10000]
[−10000,10000]
[−1000,1000]
[−1000,1000]
[−1000,1000]
[−108 ,−9.5 · 107 ]
[−4 · 104 , 5.2 · 106 ]
[−4 · 104 , 5.2 · 106 ]
#
30
30
30
30
30
30
30
90
60
90
90
90
90
Table 2: Description of the instances.
4.1
Linear programming relaxation bounds
The LP relaxation bounds are computed constructing the complete model for the
extended formulations and for formulation DFJ, since we have not implemented
a separation procedure for it. For formulation GCS, we use a Min Cut-based
separation procedure (its implementation details are left to the next section).
The tests are carried out with IBM Ilog Cplex 12.5 on an Intel Xeon E5645
@2.40GHz.
Table 3 reports the average gap of the LP relaxation bounds with respect to
LP −Opt|
. The results confirm
the optimal integer values, computed as 100 |Bound
|BoundLP |
that the LP bounds of MCF and GCS are equivalent, and show that they
are by far the tightest formulations. Remarkably, MCF and GCS close the
12
gap on a good fraction of instances (roughly 50%), while all the remaining
formulations have a relaxation with 0 gap only in less than 10% of the considered
instances. However, on the largest instances, the LP relaxation of the MCF
formulation could not be solved within the time limit of 1200 seconds. Observe
the considerable difference between GCS and classic DFJ, while formulations
RLT and SF provide similar bounds.
For the prc-f instances, all the extended formulations provide good bounds.
The reason is that, on the majority of those instances, all the arcs have very
similar costs, hence the integrality gap is inherently small for all the formulations. It is worthwile to point out that, even when the bounds are very good,
weaker formulations find solutions with many more fractional values.
GCS DFJ MTZ RLT SF
abilene 10.00 93.37 93.97 92.85 93.23
polska 37.87 67.85 69.57 61.73 63.35
nobel-us 26.48 66.60 71.96 69.07 69.72
atlanta 9.91 35.25 42.80 36.46 37.53
geant 5.52
- 32.02 30.90 31.06
france 4.33
- 64.50 59.88 60.37
germany 2.03
- 31.47 31.29 31.31
rnd-50-s 0.43
- 11.38 10.91 10.94
rnd-100-s 0.22
- 1.57 1.52 1.53
rnd-200-s 0.02
- 0.42 0.42 0.42
rnd-500-s 0.00
- 0.06 0.06 0.06
rnd-25-d 0.19
- 0.56 0.46 0.47
rnd-50-d 0.04
- 0.14 0.12 0.12
rnd-100-d 0.01
- 0.03 0.03 0.03
prc-25-f 0.00
- 0.03 0.02 0.02
prc-50-f 0.00
- 0.01 0.01 0.01
prc-100-f 0.00
- 0.00 0.00 0.00
prc-25-p 9.72
- 91.68 87.41 87.65
prc-50-p 7.92
- 80.50 77.40 77.49
prc-100-p 1.96
- 60.70 42.66 42.93
prc-25-l 1.34
- 86.73 81.06 81.39
- 54.75 48.99 49.17
prc-50-l 2.18
prc-100-l 2.44
- 60.87 42.55 42.83
MCF
10.00
37.87
26.48
9.91
5.52
4.33
2.03
0.43
0.22
0.19
0.04
0.00
0.00
9.72
7.92
1.34
2.18
-
Table 3: Average LP relaxation gaps (%) w.r.t. the optimal integer solution.
Missing values are due to time (MCF ) or memory limits (DFJ ).
4.2
Branch-and-cut
Since we aim at integer solutions, we compare the behavior of a state-of-theart MIP solver with the considered formulations. The formulations and the
separation procedures were implemented in C++ with IBM Ilog Cplex/Concert
12.5, using default settings for the branch-and-cut.
For the polynomial-size extended formulations, the full model is built. Formulation DFJ is not included in these tests. For GCS, we report results obtained
with two different separation routines. For GCS-StrongComp, the separation is
carried out for fractional and integer solutions, identifying the strongly con13
nected components in the support graph induced by the variables xij . This can
be done in a O(n + m) running time with Tarjan’s algorithm. Once a strong
component S has been found, it is enough to check if Constraint (7) is violated
for any of the nodes in S. This separation procedure is efficient, but not guaranteed to find all the violated inequalities on fractional solutions. Correctness
is preserved by the fact that the procedure is exact for integer solutions. For
GCS-MinCut, the separation is carried out on fractional solutions by solving a
sequence of Min Cut (or Max Flow) problems
between each node and t, with an
√
overall worst-case complexity of O(n3 n) using the highest-label preflow-push
algorithm described, e.g., in [AMO93]. All violated inequalities are identified,
although with higher computational cost. On the incumbents, the faster strong
components-based procedure is used.
For the Min Cut problems and the identification of strongly connected components, we use the efficient implementations in the open-source LEMON Graph
Library 1.3 [DJK11]. We refer the interested reader to [Dre13] for additional considerations on the separation of subtour elimination constraints for the ESPP.
In Table 4 we summarize the results over the test instances. Column “opt”
reports the fraction of instances that are solved to optimality within the time
limit of 1200 seconds. Column “time” reports the average computing time.
Column “nodes” reports the average number of explored nodes in the branchand-bound tree. Column “cuts” reports the average number of GCS inequalities
added by the separation procedures.
Both GCS variants solve to optimality all the instances, except the ones in
rnd-1000-s, and prove to be by far the best choice for the largest instances.
Formulation GCS-MinCut usually explores fewer nodes in the search tree, while
the number of cuts is similar. However, with respect to the average computing
time, GCS-StrongComp has an advantage over GCS-MinCut, especially on the
largest instances, where solving a sequence of Min Cut problems is computationally heavy. GCS-StrongComp is the only formulation able to solve one of
the rnd-1000-s instances, with 1000 nodes and 66601 arcs, within the time
limit.
Compared to the extended formulations, GCS-StrongComp is a clear winner
on the rnd and prc instances, where it is often more than one order of magnitude
faster. Despite the better linear relaxation bounds, formulation RLT is usually
not faster than SF, and explores a larger number of B&B nodes. Formulation SF
is able to solve a significantly larger number of the largest-size pricing instances
than RLT.
Interestingly, while on the prc instances it is rather ineffective, formulation
MTZ yields good results on the rnd instances, where it is in most cases more
efficient than all the other extended formulations. The number of explored nodes
is, nevertheless, larger. The reasons for this behavior might lie in the fact that,
on the rnd pool, all the extended formulations have similarily good LP bounds.
Such small difference might not be relevant for a modern MIP solver such as
Cplex, which is able to generate effective cutting planes, and it probably pays
off to have an LP relaxation of small size.
Formulation MCF, despite the good LP bounds, appears to be too heavy to
14
be of practical interest. On small-sized instances, very few b&b nodes are necessary thanks to the strength of the relaxation. However, this is not enough to
overcome the computational load required by the linear relaxation, as reflected
by the results on graphs with 100 or more nodes, where MCF can solve only a
small fraction of the instances within the time limit.
Figure 3 summarizes the computational experiments. On the y-axis, we
report the fraction of all the instances that are solved to optimality within the
time on the x-axis. Note that the left part of the x-axis is in linear scale, while
the right part is logarithmic.
GCS-StrongComp is the topmost curve, solving more than 80% of the instances within 20 seconds, and 90% in 50. GCS-MinCut is not far behind,
although it is generally slower. Both solve over 95% of the instances within
1200 seconds. Formulations RLT, SF and MTZ have similar results on the easiest instances, although, overall, SF is able to solve around 85% of the instances
within the time limit, while RLT and MTZ solve, respectively, less then 78%
and 72% of them. Formulation MCF (bottom curve) is already significantly
slower on the easiest instances, and solves the smallest fraction of the instances
(around 65%).
100
100
90
90
80
80
70
70
60
60
GCS-SC
GCS-MC
SF
RLT
MTZ
MCF
50
40
0
10
20
30
40
50
60
70
80
10
50
40
3
Figure 3: Fraction of instances solved to optimality within a given time. The
x-axis is in linear scale before the break (up to 85), logarithmic scale after the
break (90 to 1200).
5
Conclusions
We have analyzed integer programming formulations for the ESPP, including
formulations not yet appeared in the literature.
The polyhedral results in Section 3 provide a partial hierarchy among the
ESPP formulations, and prove that the strong extended formulation MCF has
a projection on the space of the arc variables which is equivalent to the polytope of the GCS formulation, that has exponentially many subtour elimination
15
constraints.
It is also important to understand how effective the formulations are from a
computational point of view. In this regard, we report a set of extensive computational experiments, suggesting that the extended formulations are inferior
for all practical purposes, and the dynamic separation of subtour elimination
constraints appears to be the best option when tackling the ESPP as an integer
program. The GCS formulation with the strong component-based separation
procedure is able to solve small-sized instances in a few seconds and mediumsized instances, with up to 500 nodes, within a few minutes, but it is still not
sufficient to solve large-scale problems. When the implementation of a separation procedure is not possible or convenient, formulation SF is probably the
best choice.
It seems likely that the development of good primal heuristics and strong
valid inequalities, possibly extended from of the ATSP (e.g., comb inequalities),
might further speed up the computing times and allow the solution of largesized instances. It might also be useful to borrow techniques from the typical
approaches used for the ESPPRC, such as a preprocessing phase with the aim
of reducing the search space.
References
[ACT14]
Edoardo Amaldi, Stefano Coniglio, and Leonardo Taccari, Maximum throughput network routing subject to fair flow allocation,
ISCO, Lecture Notes in Computer Science, vol. 8596, Springer, 2014,
pp. 1–12.
[AMO93]
Ravindra K Ahuja, Thomas L Magnanti, and James B Orlin, Network flows: theory, algorithms, and applications, Prentice hall,
1993.
[Bal89]
Egon Balas, The prize collecting traveling salesman problem, Networks 19 (1989), no. 6, 621–636.
[BC89]
JE Beasley and Nicos Christofides, An algorithm for the resource
constrained shortest path problem, Networks 19 (1989), no. 4, 379–
394.
[BDD06]
Natashia Boland, John Dethridge, and Irina Dumitrescu, Accelerated label setting algorithms for the elementary resource constrained
shortest path problem, Operations Research Letters 34 (2006), no. 1,
58–68.
[BHK03]
Andreas Björklund, Thore Husfeldt, and Sanjeev Khanna, Approximating longest directed path, Electronic Colloq. on Comp. Complexity, Rept, no. 32, 2003.
16
[Cla84]
A Claus, A new formulation for the travelling salesman problem,
SIAM Journal on Algebraic Discrete Methods 5 (1984), no. 1, 21–
25.
[DFJ54]
George Dantzig, Ray Fulkerson, and Selmer Johnson, Solution of
a large-scale traveling-salesman problem, Operations Research 2
(1954), no. 4, 393–410.
[DI14]
Michael Drexl and Stefan Irnich, Solving elementary shortest-path
problems as mixed-integer programs, OR spectrum 36 (2014), no. 2,
281–296.
[DJK11]
Balázs Dezső, Alpár Jüttner, and Péter Kovács, LEMON – an open
source C++ graph template library, Electronic Notes in Theoretical
Computer Science 264 (2011), no. 5, 23–45.
[Dre13]
Michael Drexl, A note on the separation of subtour elimination constraints in elementary shortest path problems, European Journal of
Operational Research 229 (2013), no. 3, 595–598.
[FDGG04] Dominique Feillet, Pierre Dejax, Michel Gendreau, and Cyrille
Gueguen, An exact algorithm for the elementary shortest path problem with resource constraints: Application to some vehicle routing
problems, Networks 44 (2004), no. 3, 216–229.
[Gab07]
Harold N Gabow, Finding paths and cycles of superpolylogarithmic
length, SIAM Journal on Computing 36 (2007), no. 6, 1648–1671.
[GG78]
Bezalel Gavish and Stephen C Graves, The travelling salesman problem and related problems, Working Paper GR-078-78, Massachusetts
Institute of Technology (1978).
[HMM13]
Mohamed Haouari, Nelson Maculan, and Mehdi Mrad, Enhanced
compact models for the connected subgraph problem and for the
shortest path problem in digraphs with negative cycles, Computers
& Operations Research 40 (2013), no. 10, 2485–2492.
[IMM09]
M.S. Ibrahim, N. Maculan, and M. Minoux, A strong flow-based
formulation for the shortest path problem in digraphs with negative cycles, International Transactions in Operational Research 16
(2009), no. 3, 361–369.
[JPS08]
Mads Kehlet Jepsen, Bjørn Petersen, and Simon Spoorendonk, A
branch-and-cut algorithm for the elementary shortest path problem
with a capacity constraint, Tech. report, Department of Computer
Science, University of Copenhagen, 2008.
[KMR97]
David Karger, Rajeev Motwani, and GDS Ramkumar, On approximating the longest path in a graph, Algorithmica 18 (1997), no. 1,
82–98.
17
[MTZ60]
Clair E Miller, Albert W Tucker, and Richard A Zemlin, Integer
programming formulation of traveling salesman problems, Journal
of the ACM (JACM) 7 (1960), no. 4, 326–329.
[NR02]
Denis Naddef and Giovanni Rinaldi, Branch-and-cut algorithms for
the capacitated VRP, ch. 3, pp. 53–84, SIAM Monographs on Discrete Mathematics and Applications Philadelphia, PA, 2002.
[OWPT10] S. Orlowski, R. Wessäly, M. Pióro, and A. Tomaszewski, SNDlib
1.0 - survivable network design library, Networks 55 (2010), no. 3,
276–286.
[PS91]
Manfred Padberg and Ting-Yi Sung, An analytical comparison of
different formulations of the travelling salesman problem, Mathematical Programming 52 (1991), no. 1-3, 315–357.
[RS06]
Giovanni Righini and Matteo Salani, Symmetry helps: bounded bidirectional dynamic programming for the elementary shortest path
problem with resource constraints, Discrete Optimization 3 (2006),
no. 3, 255–273.
[VSO11]
Pieter Vansteenwegen, Wouter Souffriau, and Dirk Van Oudheusden, The orienteering problem: A survey, European Journal of Operational Research 209 (2011), no. 1, 1–10.
[Won80]
R.T. Wong, Integer programming formulations of the travelling
salesman problem, Proceedings IEEE Conference on Circuits and
Computers, 1980, pp. 149–152.
18
19
abilene
polska
nobel-us
atlanta
geant
france
germany
rnd-50-s
rnd-100-s
rnd-200-s
rnd-500-s
rnd-1000-s
rnd-25-d
rnd-50-d
rnd-100-d
prc-25-f
prc-50-f
prc-100-f
prc-25-p
prc-50-p
prc-100-p
prc-25-l
prc-50-l
prc-100-l
opt
100
100
100
100
100
100
100
100
30
0
0
0
100
86.7
0
100
70
0
100
93.3
0
100
86.7
0
opt
100
100
100
100
100
100
100
100
100
100
100
0.0
100
100
100
100
100
100
100
100
100
100
100
100
GCS-MinCut
RLT
SF
MTZ
time nodes
cuts opt
time
nodes opt
time
nodes opt
time
nodes
0
0.6
16.1 100
0
0.2 100
0
0 100
0
1.4
0
2.6
20.2 100
0
7.2 100
0
0.3 100
0
32.2
0
2.1
21.8 100
0
8.1 100
0
2.9 100
0
23.7
0
0.8
22 100
0
1.7 100
0
0.5 100
0
17.9
0.1
2.6
41.1 100
0.1
6.7 100
0.1
1.3 100
0
22.1
0.1
1.1
46 100
0.1
30.8 100
0.2
8.2 100
0.1
279.8
1.4
21.7 392.6 100
0.5
177.5 100
0.7
207.3 100
4.8
8189.5
0.3
12.1 144.7 100
0.3
119.1 100
0.3
48.9 100
0.2
118.7
2.2
33.9
239 100
5
529.2 100
3.7
326.9 100
2.4
652.5
14.7
38.7 367.5 100
55.7
930.8 100
10.4
103.8 100
11.1
655.2
238.1
52.3
707 5.6 1173.1
326.8 72.2 615.8
428.9 100
355
2606
1200
0.0 1275.5 0.0
1200
0.0 0.0
1200
0.0 0.0
1200
46.1
0.1
3.7
39.6 100
0.2
9.9 96.7
0.1
5.3 100
0.1
15.1
0.4
11.4
78.1 100
1.5
51 100
1.1
50.1 100
2.4
686.5
3.6
20.5
147 100
37.2
219.8 100
27.8
188 100
11.4
147.7
0.1
12
46.2 100
9.8
1348.5 100
0.4
17.1 100
9
13785
1.3
20.9 123.5 96.7
100
5314.6 100
23.2
1601.5
80 339.8 1.6 · 105
100.6 2436.6 415.3
10 1156.8 16036.4
50
839
11586 6.7
1148 98779.5
0.2
5.1
63.3 100
1.8
3444.7 100
1
358.9 53.3 668.7 1.1 · 106
1.7
7.2 137.1 100
24.5 23169.4 100
44.7
6617.4
40 937.1 7.4 · 105
115.8 730.6 727.8
0
1200 2.8 · 105 13.3 1155.8 2.3 · 105
0
1200 2.5 · 105
0.3
1.4
80.6 100
51.6 25478.5 100
1.5
119.7 33.3 820.8 1.4 · 106
2.6
12 184.1
50 790.4 79178.6 96.7 217.1 24884.5
10 1085.1 6.8 · 105
130.6 1489.1 746.9
0
1200 2.9 · 105 6.7 1184.5 2.3 · 105
0
1200 2.6 · 105
Table 4: Branch-and-cut results.
MCF
GCS-StrongComp
time nodes opt
time nodes
cuts
0
0 100
0
0.8
16.7
0
0 100
0
3.5
20.7
0.1
0 100
0
2.8
21.8
0.1
0 100
0
1.3
22.3
0.2
0 100
0.1
3.6
45.8
0.4
0 100
0.1
1.7
44.5
8.7
7.9 100
0.9
58.7 379.6
4.7
9 100
0.2
40.2 152.3
886 17.2 100
1.1
67.2
221
1200
0 100
6.6
46.1 351.4
1200
0 100 156.3
79.3 777.7
1200
0 3.3 1192.0
6.4 1417.6
3.4
3 100
0.1
5.1
40.2
479.1
7.7 100
0.4
12.1
85.8
1200
0 100
3.1
26.2
159
6.1
5.1 100
0.1
10
42.8
748.9 13.3 100
0.6
37.5 100.2
1200
0 100
13.7 2061.3 276.7
6
3.6 100
0.2
14
64.5
406.1
3.8 100
1
29.4 128.6
1200
0 100
34.2 2242.9 549.2
6.8
0.4 100
0.3
16.6
77.9
415.7
9.6 100
1.8
85.7 224.4
1200
0 100
38.2 3349.4 508.1

Download Report