supply chain

2006_Original.pdf

Home >Business & Finance homework help >Operations Management homework help >supply chain

This article was downloaded by: [107.133.16.252] On: 28 October 2021, At: 19:43 Publisher: Institute for Operations Research and the Management Sciences (INFORMS) INFORMS is located in Maryland, USA

Operations Research

Publication details, including instructions for authors and subscription information: http://pubsonline.informs.org

A Branch-and-Cut Algorithm for the Dial-a-Ride Problem Jean-François Cordeau,

To cite this article: Jean-François Cordeau, (2006) A Branch-and-Cut Algorithm for the Dial-a-Ride Problem. Operations Research 54(3):573-586. https://doi.org/10.1287/opre.1060.0283

Full terms and conditions of use: https://pubsonline.informs.org/Publications/Librarians-Portal/PubsOnLine-Terms-and- Conditions

This article may be used only for the purposes of research, teaching, and/or private study. Commercial use or systematic downloading (by robots or other automatic processes) is prohibited without explicit Publisher approval, unless otherwise noted. For more information, contact [email protected].

The Publisher does not warrant or guarantee the article’s accuracy, completeness, merchantability, fitness for a particular purpose, or non-infringement. Descriptions of, or references to, products or publications, or inclusion of an advertisement in this article, neither constitutes nor implies a guarantee, endorsement, or support of claims made of that product, publication, or service.

Please scroll down for article—it is on subsequent pages

With 12,500 members from nearly 90 countries, INFORMS is the largest international association of operations research (O.R.) and analytics professionals and students. INFORMS provides unique networking and learning opportunities for individual professionals, and organizations of all types and sizes, to better understand and use O.R. and analytics tools and methods to transform strategic visions and achieve better outcomes. For more information on INFORMS, its publications, membership, or meetings visit http://www.informs.org

http://pubsonline.informs.org

https://doi.org/10.1287/opre.1060.0283

https://pubsonline.informs.org/Publications/Librarians-Portal/PubsOnLine-Terms-and-Conditions

http://www.informs.org

OPERATIONSRESEARCH Vol. 54, No. 3, May–June 2006, pp. 573–586 issn0030-364X�eissn1526-5463�06�5403�0573

informs ® doi10.1287/opre.1060.0283

A Branch-and-Cut Algorithm for the Dial-a-Ride Problem

Jean-François Cordeau HEC Montréal, 3000 chemin de la Côte-Sainte-Catherine, Montréal, Quebec, Canada H3T 2A7, [email protected]

In the dial-a-ride problem, users formulate requests for transportation from a specific origin to a specific destination. Transportation is carried out by vehicles providing a shared service. The problem consists of designing a set of minimum- cost vehicle routes satisfying capacity, duration, time window, pairing, precedence, and ride-time constraints. This paper introduces a mixed-integer programming formulation of the problem and a branch-and-cut algorithm. The algorithm uses new valid inequalities for the dial-a-ride problem as well as known valid inequalities for the traveling salesman, the vehicle routing, and the pick-up and delivery problems. Computational experiments performed on randomly generated instances show that the proposed approach can be used to solve small to medium-size instances.

Subject classifications: transportation: vehicle routing; programming: cutting plane. Area of review: Transportation. History: Received September 2003; revision received August 2004; accepted March 2005.

1. Introduction In the Dial-a-Ride Problem (DARP), users formulate re- quests for transportation from a specific origin (or pick-up point) to a specific destination (or drop-off point). Trans- portation is carried out by vehicles that provide a shared service in the sense that several users may be in a vehicle at the same time. The aim is to design a minimum-cost set of vehicle routes accommodating all requests under a number of side constraints. A common DARP application arises in door-to-door transportation services for the elderly and the disabled. In this context, users often formulate two requests per day: an outbound request from home to a destination, and an inbound request for the return trip. Most dial-a-ride services are characterized by the pres-

ence of two conflicting objectives: minimizing operating costs and minimizing user inconvenience. Operating costs are mostly related to fleet size and distance traveled, while user inconvenience is often measured in terms of devia- tions from desired pick-up and drop-off times and in terms of excess ride time (i.e., the difference between the actual ride time of a user and the minimum possible ride time). One way to achieve a balance between these objectives is to treat cost minimization as the primary objective and to impose service quality constraints. As in the work of Jaw et al. (1986) and Cordeau and

Laporte (2003b), we assume here that the user specifies either a desired arrival time at destination (in the case of an outbound request) or a desired departure time from the origin (in the case of an inbound request). In both cases, a time window of a prespecified width is constructed around the desired time. In addition, an upper bound is imposed on the ride time of the user. This approach seems to be in line with the current practice of several North American

transporters. For example, in a system with 15-minute time windows and a maximum ride time of 45 minutes, a user wishing to arrive at destination at 9:00 would be picked up no earlier than 8:00 and dropped off between 8:45 and 9:00. It should be emphasized that imposing time windows is not sufficient to accurately impose a maximum ride time. In the above example, a pickup at 8:00 and a drop-off at 9:00 would exceed the maximum ride time of 45 minutes. How- ever, constraining the pickup to take place no earlier than 8:15 would be too restrictive because it would, in particular, prohibit a pickup at 8:00 and a drop-off at 8:45. Dial-a-ride services may operate according to a static or

to a dynamic mode. In the first case, all requests are known beforehand, while in the second case requests are gradually received throughout the day and vehicle routes are adjusted in real time to meet demand. In practice, pure dynamic problems rarely exist because a large subset of requests is often known in advance. This paper deals with the static variant of the problem. Because it incorporates time windows and maximum

ride-time constraints, the DARP is a difficult problem that generalizes the Vehicle Routing Problem with Pickup and Delivery (VRPPD). Finding a feasible solution for the DARP is itself NP-hard because it also general- izes the Traveling Salesman Problem with Time Windows (TSPTW) (see, e.g., Savelsbergh 1985). As a result, most previous work has concentrated on the development of heuristics. Nevertheless, when the problem is moderately constrained, exact solution approaches can be devised to solve small to medium-size instances. Our aim is to describe valid inequalities and a branch-

and-cut algorithm for the DARP. Branch-and-cut has proved to be an effective technique for solving the TSP, the

573

VRP, and several other routing problems (see, e.g., Naddef and Rinaldi 2002). The algorithm proposed here uses adap- tations of known valid inequalities for the precedence- constrained TSP, the VRP, and the VRPPD, as well as new inequalities that take advantage of the special structure of the problem. The remainder of this paper is organized as follows. The

next section briefly reviews relevant work on the DARP and closely related problems. Section 3 defines the DARP formally and introduces a mixed-integer formulation. Sec- tion 4 introduces several families of valid inequalities used in the branch-and-cut algorithm, which is then described in §5, along with separation heuristics and preprocessing techniques. Computational experiments are reported in §6 and the conclusion follows in §7.

2. Literature Review Early research on the DARP was carried out by Psaraftis (1980, 1983), who developed dynamic programming algo- rithms for the single-vehicle case. An improved labeling scheme was later proposed by Desrosiers et al. (1986), who were able to solve instances with up to 40 users. Exact algorithms have also been devised for closely

related single-vehicle problems. Ruland (1995) and Ruland and Rodin (1997) proposed a branch-and-cut algorithm for a special case of the pick-up and delivery problem in which there are no capacity or time window constraints. In this case, if distances satisfy the triangle inequality, a single vehicle will serve all users. As a result, the problem reduces to a TSP with the additional constraint that the pick-up node of each user must precede his drop-off node in the vehicle route. These authors provided an integer program- ming formulation of the problem based on an undirected graph and described four families of valid inequalities: sub- tour elimination constraints, precedence constraints, gen- eralized order constraints, and order-matching constraints. These inequalities are also valid for the more general prob- lem studied in this paper, and they will be discussed in more depth in §4. The more general precedence-constrained asymmetric

TSP (in which each node may have multiple predecessors and successors) has been studied, among others, by Balas et al. (1995) and Ascheuer et al. (2000), who proposed several families of valid inequalities based, in large part, on the strengthening of existing inequalities for the asym- metric TSP. This problem has applications, for example, in the scheduling of flexible manufacturing systems (see, e.g., Ascheuer 1996). Because the TSP with pickup and delivery is a particular case of the precedence-constrained TSP, some of the proposed inequalities are also useful for solving the DARP. In particular, predecessor and successor inequalities proposed by Balas et al. will be used in our approach. Most of the algorithms for the multiple-vehicle DARP

are heuristics or metaheuristics. An insertion method capa- ble of handling large-scale instances was described by

Jaw et al. (1986), while Bodin and Sexton (1986) developed an exchange heuristic that uses Benders decomposition to optimize the individual vehicle routes. Dumas et al. (1989) later presented a clustering and column generation method that can handle instances with several thousand users. More recently, Madsen et al. (1995) proposed an insertion heuris- tic for the dynamic case, while Toth and Vigo (1996, 1997) used local search and a tabu-thresholding proce- dure to address a real-life problem arising in Bologna. Finally, Cordeau and Laporte (2003b) described a tabu search heuristic for the variant of the problem addressed in this paper. This heuristic will be used in §6 to compute upper bounds. Relevant work was also performed on the closely related

VRPPD with time windows, arising in contexts such as urban courier services. For this problem, Dumas et al. (1991) presented an exact column generation method using the dynamic programming algorithm of Desrosiers et al. (1986) to solve the pricing subproblem. A similar approach was also described by Savelsbergh and Sol (1998), who proposed several refinements to improve performance. For recent overviews of the DARP and VRPPD, the

reader is referred to the surveys of Cordeau and Laporte (2003a) and Desaulniers et al. (2002), respectively.

3. Formulation Let n denote the number of users (or requests) to be served. The DARP may be defined on a complete directed graph G = �N�A�, where N = P ∪ D ∪ 0�2n + 1�, P = 1��n�, and D = n+1��2n�. Subsets P and D contain pick-up and drop-off nodes, respectively, while nodes 0 and 2n + 1 represent the origin and destination depots. With each user i are thus associated an origin node i and a destination node n+ i. Let K be the set of vehicles. Each vehicle k ∈ K has

a capacity Qk, and the total duration of its route can- not exceed Tk. With each node i ∈ N are associated a load qi and a nonnegative service duration di such that q0 = q2n+1 = 0, qi = −qn+i �i = 1��n�, and d0 = d2n+1 = 0. A time window �ei�li� is also associated with node i ∈ N , where ei and li represent the earliest and latest time, respec- tively, at which service may begin at node i. With each arc �i�j� ∈ A are associated a routing cost cij and a travel time tij. Finally, denote by L the maximum ride time of a user. As mentioned in the introduction, we assume here that a

time window is specified either for the origin node or for the destination node of a request, but not both. However, a valid time window can always be determined for the other node by taking the direct and the maximum ride times into account. We will explain in §5.1.1 how these time windows can be made as tight as possible. For each arc �i�j� ∈ A and each vehicle k ∈ K, let xkij = 1

if vehicle k travels from node i to node j. For each node i ∈ N and each vehicle k ∈ K, let Bki be the time at which

vehicle k begins service at node i, and Qki be the load of vehicle k after visiting node i. Finally, for each user i, let Lki be the ride time of user i on vehicle k. The DARP can be formulated as the following mixed-

integer program:

Min ∑ k∈K

∑ i∈N

∑ j∈N ckijx

k ij (1)

subject to

∑ k∈K

∑ j∈N xkij = 1 ∀i ∈ P� (2)

∑ j∈N xkij −

∑ j∈N xkn+i�j = 0 ∀i ∈ P� k ∈ K� (3)

∑ j∈N xk0j = 1 ∀k ∈ K� (4)

∑ j∈N xkji −

∑ j∈N xkij = 0 ∀i ∈ P ∪D� k ∈ K� (5)

∑ i∈N xki�2n+1 = 1 ∀k ∈ K� (6)

Bkj � �B k i +di + tij�xkij ∀i ∈ N� j ∈ N� k ∈ K� (7)

Qkj � �Q k i +qj�xkij ∀i ∈ N� j ∈ N� k ∈ K� (8)

Lki = Bkn+i −�Bki +di� ∀i ∈ P� k ∈ K� (9) Bk2n+1 −Bk0 � Tk ∀k ∈ K� (10) ei � B

k i � li ∀i ∈ N� k ∈ K� (11)

ti�n+i � L k i � L ∀i ∈ P� k ∈ K� (12)

max 0�qi��Q k i �min Qk�Qk+qi� ∀i∈N�k∈K� (13)

xkij ∈ 0�1� ∀i ∈ N� j ∈ N� k ∈ K� (14)

The objective function (1) minimizes the total routing cost. Constraints (2) and (3) ensure that each request is served exactly once and that the origin and destination nodes are visited by the same vehicle. Constraints (4)–(6) guarantee that the route of each vehicle k starts at the ori- gin depot and ends at the destination depot. Consistency of the time and load variables is ensured by constraints (7) and (8). Equalities (9) define the ride time of each user, which is bounded by constraints (12). It is worth men- tioning that the latter also act as precedence constraints because the nonnegativity of the Lki variables ensures that node i will be visited before node n + i for every user i. Finally, inequalities (10) bound the duration of each route, while (11) and (13) impose time windows and capacity constraints, respectively. This formulation is nonlinear because of constraints (7)

and (8). Introducing constants Mkij and W k ij , these constraints

can, however, be linearized as follows:

Bkj �B k i +di+tij −Mkij�1−xkij� ∀i∈N�j∈N�k∈K� (15)

Qkj � Q k i +qj −Wkij�1−xkij� ∀i ∈ N� j ∈ N� k ∈ K� (16)

The validity of these constraints is ensured by setting Mkij � max 0�li +di +tij −ej� and Wkij � min Qk�Qk +qi�. These constraints are very similar to the Miller-Tucker-Zemlin subtour elimination constraints for the TSP (Miller et al. 1960). One can reduce the number of variables and constraints

in Model (1)–(14) by using aggregate time variables Bi at every node except the origin depot 0 and the destination depot 2n+1. In this case, one replaces (7) and (9) with the constraints

Bj � �B k 0 +d0 + t0j�xk0j ∀j ∈ N� k ∈ K� (17)

Bj � �Bi +di + tij� ∑ k∈K xkij ∀i ∈ N� j ∈ N� (18)

Bk2n+1 � �Bi +di + ti�2n+1�xki�2n+1 ∀i ∈ N� k ∈ K� (19) Li = Bn+i −�Bi +di� ∀i ∈ P� (20) to which the same linearization process can be applied. Along the same lines, if the fleet of vehicles is homoge- neous in the sense that Qk = Q for every k ∈ K, one can replace (8) with

Qj � �Q k 0 +qj�xk0j ∀j ∈ N� k ∈ K� (21)

Qj � �Qi +qj� ∑ k∈K xkij ∀i ∈ N� j ∈ N� (22)

Qk2n+1 � �Qi +q2n+1�xki�2n+1 ∀i ∈ N� k ∈ K� (23) Finally, as shown by Desrochers and Laporte (1991), the

linearized form of constraints (22) can be lifted as follows by taking the reverse arc �j�i� into account:

Qj � Qi +qj −Wij ( 1− ∑

k∈K xkij

) +�Wij −qi −qj�

∑ k∈K xkji

∀i ∈ N� j ∈ N� k ∈ K� (24) Observe that in our case an equivalent lifting cannot be performed with constraints (18) because waiting will some- times take place after the beginning of a time window (i.e., Bi >ei) to reduce the ride time of a user.

4. Valid Inequalities We now describe several families of valid inequalities for the DARP. All of these inequalities are of course redundant for Model (1)–(14) but can strengthen its LP-relaxation. Because of the structure of the model, analyzing whether these inequalities define facets of the DARP polytope appears to be a rather challenging task. However, their use- fulness is demonstrated through computational experiments in the last section. The following additional notation will be used to

describe the valid inequalities. Given a node set S ⊆ N , define S̄ = i ∈ N � i S� and "�S� = "+�S� ∪ "−�S�, where "+�S� = �i�j� ∈ A � i ∈ S� j S� and "−�S� = �i�j� ∈ A � i S� j ∈ S�. For notational convenience, let xij =

∑ k∈K x

k ij denote the total flow on arc �i�j� and define

x�S� = ∑i�j∈S xij. Similarly, let x�A′� = ∑�i�j�∈A′ xij for any arc set A′ ⊆ A.

4.1. Bounds on Time and Load Variables

As first suggested by Desrochers and Laporte (1991) in the context of the TSP with time windows, bounds on the time variables Bi can be strengthened as follows:

Bi � ei + ∑

j∈N\ i� max 0�ej −ei +dj + tij�xji� (25)

Bi � li − ∑

j∈N\ i� max 0�li − lj +di + tij�xij� (26)

These inequalities were used, for example, for solving the asymmetric TSP with time windows by branch and cut (Ascheuer et al. 2001). Similarly, bounds on load variables Qi can also be

strengthened as follows:

Qi � max 0�qi�+ ∑

j∈N\ i� max 0�qj�xji� (27)

Qi � min Q�Q+qi�− ( Q− max

j∈N\ i� qj�−qi

) x0i

− ∑ j∈N\ i�

max 0�qj�xij� (28)

4.2. Subtour Elimination Constraints

Consider the simple subtour elimination constraint x�S� � �S�−1 for S ⊆ P ∪D. In the case of the DARP, this inequal- ity can be lifted in many different ways by taking into account the fact that for each user i, node i must be visited before node n+ i. For any set S ⊆ P ∪D, let #�S� = i ∈ P � n+i ∈ S� and

$�S� = n + i ∈ D � i ∈ S� denote the sets of predecessors and successors of S, respectively. Balas et al. (1995) have proposed two families of inequalities for the precedence- constrained asymmetric TSP that also apply to the DARP by observing that each node i ∈ P ∪D is either the prede- cessor or the successor of exactly one other node. For S ⊆ P ∪D, the following inequality, called a succes-

sor inequality (or $-inequality) is valid for the DARP:

x�S�+ ∑ i∈S̄∩$�S�

∑ j∈S xij +

∑ i∈S̄\$�S�

∑ j∈S∩$�S�

xij � �S� −1� (29)

Example. Consider the node set S = i�j� ⊆ P. The result- ing inequality is xij +xji +xn+i�j +xn+j�i � 1. This example is illustrated in Figure 1, where dotted arcs correspond to lifted arcs (i.e., arcs that are not part of the classical sub- tour elimination constraint). Observe that variables xn+i�i and xn+j�j are not to be introduced in the lifted inequal- ity because the corresponding arcs can be trivially removed from the graph. Similarly, for any set S ⊆ P ∪D, the following predeces-

sor inequality (or #-inequality) is valid for the DARP:

x�S�+ ∑ i∈S

∑ j∈S̄∩#�S�

xij + ∑

i∈S∩#�S�

∑ j∈S̄\#�S�

xij � �S� −1� (30)

Figure 1. Successor inequality for S = i�j� ⊆ P.

j n + j

n + i

Example. Consider the node set S = n + i�n + j� ⊆ D. The resulting predecessor inequality is xn+i�n+j +xn+j�n+i + xn+i�j +xn+j�i � 1. This is illustrated in Figure 2. In the case of a directed formulation, one can also

lift subtour elimination constraints by taking into account the orientation of the arcs. For an ordered set S = i1�i2��ih� ⊆ N with h � 3 nodes, Grötschel and Padberg (1985) proposed the following inequalities, called D+k and D

− k , respectively, for the asymmetric TSP:

h−1∑ j=1 xij�ij+1 +xih�i1 +2

h−1∑ j=2 xij�i1 +

h−1∑ j=3

j−1∑ l=2 xij�il � h−1� (31)

h−1∑ j=1 xij�ij+1 +xih�i1 +2

h∑ j=3 xi1�ij +

h∑ j=4

j−1∑ l=3 xij�il � h−1� (32)

Of course, different liftings are obtained by considering different orderings of the nodes in set S. In the case of the DARP, these inequalities can be further lifted by taking precedence relationships into account. This results in the following two propositions.

Proposition 1. Let S = i1�i2��ih� ⊆ P ∪ D. The fol- lowing inequality is valid for the DARP:

h−1∑ j=1 xij�ij+1 +xih�i1 +2

h−1∑ j=2 xij�i1 +

h−1∑ j=3

j−1∑ l=2 xij�il

+ ∑ n+ip∈S̄∩$�S�

xn+ip�i1 � h−1� (33)

Figure 2. Predecessor inequality for S = n+ i�n+ j� ⊆ D.

in + i

n + j

Figure 3. Lifted directed subtour elimination constraint for S = i�j�k� ⊆ P.

n + j

kn + k

Proof. Suppose that one arc of the form �n+ ip�i1� with n+ ip ∈ S̄ ∩$�S� is part of the solution. Then, all arcs of the form �ij�i1� with 2 � j � h − 1 cannot belong to the solution. As a result, if the left-hand side of inequality (33) is larger than h−1, then there exists a subpath linking the h nodes in set S. However, because set S contains the origin node ip, this subpath, together with the arc �n + ip�i1�, would violate the precedence constraint for user ip. �

Example. Consider the node set S = i�j�k� ⊆ P. One possible lifted directed subtour elimination constraint (obtained with i1 = i, i2 = j, i3 = k) is xij + xjk + xki + 2xji +xn+j�i +xn+k�i � 2. This is illustrated in Figure 3. Proposition 2. Let S = i1�i2��ih� ⊆ P ∪ D. The fol- lowing inequality is valid for the DARP:

h−1∑ j=1 xij�ij+1 +xih�i1 +2

h∑ j=3 xi1�ij +

h∑ j=4

j−1∑ l=3 xij�il

+ ∑ ip∈S̄∩#�S�

xi1�ip � h−1� (34)

Proof. The proof is similar to that of Proposition 3 by observing that if one arc of the form �i1�ip� with ip ∈ S̄ ∩ #�S� is part of the solution, then all arcs of the form �i1�ij� with 3 � j � h cannot belong to the solution. �

Example. Consider the node set S = n + i�n + j� n + k� ⊆ D. One possible lifted directed subtour elimina- tion constraint (obtained with i1 = n + i, i2 = n + j, i3 = n+k) is xn+i�n+j +xn+j�n+k +xn+k�n+i +2xn+i�n+k +xn+i�j + xn+i�k � 2. This is illustrated in Figure 4. For the case h = 3, Ascheuer (1996) proposed slightly

stronger liftings. These do not, however, generalize to the case h>3.

4.3. Capacity Constraints

For any subset S ⊆ P ∪ D, let R�S� denote the minimum number of vehicles needed to visit all nodes in S. The con- straint x�"�S�� 2R�S� is then a valid inequality. Although

Figure 4. Lifted directed subtour elimination constraint for S = n+ i�n+j�n+k� ⊆ D.

n + i

n + j

n + k 2

computing R�S� is difficult, a lower approximation is pro- vided by �q�S�/Q�, where q�S� = ∑i∈S qi. The resulting inequality is called a rounded capacity inequality in the context of the VRP. While capacity inequalities play an important role in

most branch-and-cut algorithms for the VRP (see, e.g., Naddef and Rinaldi 2002), they are not likely to be very strong for the DARP because of the presence of both pos- itive and negative qi values. In particular, q�P ∪D� = 0 by definition and, in the absence of time windows, all nodes can be visited by the same vehicle. Useful inequalities can, however, be obtained by restricting S to be a subset of either P or D and setting q�S� = �∑i∈S qi�. It should be emphasized that in our context, the quantity �q�S�/Q� esti- mates the required number of vehicle passages in set S rather than the number of vehicles per se. Indeed, by leav- ing and entering set S more than once, the same vehi- cle may be able to visit all nodes in the set even though q�S�>Q.

4.4. Precedence Constraints

Consider a node set S such that 0 ∈ S, i S, n+ i ∈ S, and 2n+1 S for at least one user i. Because node i must be visited before node n+ i, x�S� � �S� −2 and, equivalently, x�"�S�� 3. The same applies to node sets S such that 0 S, i ∈ S, n+i S, and 2n+1 ∈ S for at least one user i. These inequalities, which are also valid for the DARP, were introduced by Ruland and Rodin (1997) for the pick-up and delivery problem.

4.5. Generalized Order Constraints

Let U1��Um ⊂ N be mutually disjoint subsets and let i1��im ∈ P be users such that 0�2n + 1 Ul and il, n+ il+1 ∈ Ul for l = 1��m (where im+1 = i1). The fol- lowing inequality, introduced by Ruland and Rodin (1997), is also valid for the DARP:

m∑ l=1 x�Ul� �

m∑ l=1

�Ul� −m−1� (35)

Similar inequalities, called precedence cycle breaking in- equalities, have also been proposed by Balas et al. (1995). In the directed case, generalized order constraints can

be lifted in two different ways, as shown by the following propositions.

Proposition 3. The following inequality is valid for the DARP:

m∑ l=1 x�Ul�+

m−1∑ l=2 xi1�il +

m∑ l=3 xi1�n+il �

m∑ l=1

�Ul� −m−1� (36)

Proof. First observe that in any feasible integer solution, at most one of the lifted arcs may be part of the solution because they all belong to "+�i1�. To demonstrate the valid- ity of the lifting, we show that if any of the lifted arcs are part of the solution, then there are at least two subsets Ul for which x�Ul� � �Ul�−2. Because x�Ul� � �Ul�−1 for all other subsets, the validity of the inequality follows directly. The details of the proof are given in the appendix. �

Remark 1. The numbering of the sets U1��Um is arbi- trary and leads, by symmetry, to m! possibly different liftings.

Proposition 4. The following inequality is valid for the DARP:

m∑ l=1 x�Ul�+

m−2∑ l=2 xn+i1�il +

m−1∑ l=2 xn+i1�n+il �

m∑ l=1

�Ul� −m−1� (37)

Proof. The reasoning is similar to that of Proposition 5, and the details of the proof are given in the appendix. �

Example. Consider the subsets U1 = i�n + j�, U2 = j�n+k�, and U3 = k�n+i� with i1 = i, i2 = j, and i3 = k. The lifted generalized order constraint is xi�n+j + xn+j�i + xj�n+k +xn+k�j +xk�n+i +xn+i�k +xn+i�n+j � 2. This is illus- trated in Figure 6. Observe that in inequality (37), the value of m must be larger than or equal to 4 for the second term of the inequality to have any effect.

Figure 5. Lifted generalized order constraint with m = 3 (first lifting).

n + j n + k

j k

n + i

Figure 6. Lifted generalized order constraint with m = 3 (second lifting).

n + j n + k

j k

n + i

4.6. Order-Matching Constraints

Consider two nodes i�j ∈ P and a subset H such that i�j� ⊆ H ⊆ N\ 0�n + i�n + j�2n + 1�. The following inequality, introduced by Ruland and Rodin (1997), is also valid for the DARP:

x�H�+x� i�n+ i��+x� j�n+j�� H�� (38) As suggested by Ruland (1995), this inequality can be

lifted by introducing the term x� h�n+h�� for every user h such that h ∈ H and n+h H. In the directed case, order- matching constraints are in fact redundant because any arc of the form �n+h�h� can be removed from the graph. As a result, all arcs in the left-hand side of (38) have their origin node in set H. Hence, the sum of the flows on these arcs cannot exceed �H�, even in a fractional solution. Order-matching constraints can, however, be generalized

by replacing the arcs �h�n+h� by node sets. This leads to the following proposition.

Proposition 5. Let i1��im be m users and let H ⊂ P ∪ D and Th ⊂ P ∪ D, h = 1��m be node sets such that ih�n+ih� ⊆ Th and H∩Th = ih�. The following inequality is valid for the DARP:

x�H�+ m∑ h=1 x�Th�+ � �H� +

m∑ h=1

�Th� −2m� (39)

Proof. First observe that x�H� � �H� − 1 and x�Th� � �Th� −1 for h = 1��m. If x�Th� = �Th� −1 for any given subset Th, then there exists a path connecting all nodes in Th, including ih and n + ih. However, this path cannot finish at node ih because of the precedence constraint for user ih. Let , be the number of sets Th for which x�Th� = �Th� − 1. Then, x�"+�H�� ,. Because x�"+�H�� = x�"−�H�� and 2x�H�+x�"+�H��+x�"−�H�� = 2�H�, one obtains that 2x�H� � 2�H� −2, and thus x�H� � �H� −,. Finally, because x�Th� � �Th� − 2 for the remaining m−, sets, one may conclude that x�H� + ∑mh=1 x�Th� � �H� − ,+ ∑mh=1��Th� −1�−�m−,�, which simplifies to expres- sion (39). �

Example. Consider the sets H = i�j�, T1 = i�n + i�k� and T2 = j�n + j�l� with i1 = i and i2 = j. The resulting inequality is x�H�+x�T1�+x�T2� � 4 and is illustrated in Figure 7.

Figure 7. Generalized order-matching constraint with m = 2.

n + i

n + j

Remark 2. Inequality (39) is stronger than the correspond- ing TSP comb constraint, valid for m � 3 and odd, which is defined as x�H� + ∑mh=1 x�Th� � �H� + ∑mh=1 �Th� − �3m+1�/2. Inequalities (39) can be further lifted by introducing the

arcs from any node in H to any node outside of H. Nev- ertheless, in the presence of lifted subtour elimination con- straints (30), these inequalities are redundant, as will be shown in the following proposition.

Proposition 6. Define �H = ⋃mh=1 ih�. The following inequality, which is a strengthening of (39), is redundant in the presence of inequalities (30):

x�H�+ m∑ h=1 x�Th�+

∑ i∈ �H

∑ j∈ �H\Ti

xij + ∑ i∈H\ �H

∑ j∈ �H xij

� �H� + m∑ h=1

�Th� −2m� (40)

Proof. Inequality (40) can be written equivalently as

∑ i∈H

∑ j∈N xij +

m∑ h=1

∑ i∈Th\ ih�

∑ j∈Th xij � �H� +

m∑ h=1

�Th� −2m�

For every subset Th, h = 1��m, one has that x�Th\ ih�� +

∑ i∈Th\ ih�xi�ih � �Th� − 2 because of inequal-

ity (30) defined for subset S = Th\ ih�, whose cardinality is �Th� − 1. As a result,

∑ i∈Th\ ih�

∑ j∈Th xij � �Th� − 2.

In addition, one has that ∑ i∈H

∑ j∈N xij = �H� by con-

straints (2) and (3). The result follows directly from these two observations. �

4.7. Infeasible Path Constraints

Ride-time constraints may give rise to paths that are infea- sible in an integer solution, but nonetheless feasible in a fractional solution. Forbidding such paths can be accom- plished through the valid inequality introduced in the next proposition.

Proposition 7. Assume that the triangle inequality holds for travel times. For any directed path � = i�k1�k2�� kp�n+i� such that ti�k1 +dk1 +tk1�k2 +dk2 +···+tkp�n+i > L, the following inequality is valid for the DARP:

xi�k1 + p−1∑ h=1 xkh�kh+1 +xkp�n+i � p−1� (41)

Proof. Suppose that the value of the left-hand side of inequality (41) is equal to p in a feasible integer solution. Then, there is a single arc from path � not belonging to the solution because the path contains p+1 arcs. Because the solution is feasible, it must contain a path from i to n+i in which the missing arc has been replaced by at least two other arcs. However, because the triangle inequality holds for travel times, the path from i to n+ i has a larger duration than that of �. As a result, its duration must exceed L, which contradicts the assumption that the solu- tion is feasible. �

Remark 3. As will be shown in the next section, the case p = 1 can be handled directly through the simple elimina- tion of arcs in a preprocessing step.

5. Branch-and-Cut Algorithm Branch and cut is a popular approach for solving combina- torial problems. It has been applied successfully to several routing problems (see, e.g., Ascheuer et al. 2001, Gendreau et al. 1998, and Naddef and Rinaldi 2002). This section describes the branch-and-cut algorithm used for the DARP. It focuses on the preprocessing techniques developed to reduce problem size and on the separation heuristics used to identify violated inequalities. After applying the preprocessing steps presented in §5.1

and generating an initial pool of inequalities, the algorithm first solves the LP relaxation of the problem. If the solution to the LP relaxation is integer, an optimal solution has been identified. Otherwise, an enumeration tree is constructed and violated valid inequalities are generated at some nodes of this tree by means of the separation heuristics described in §5.3. Because all inequalities described in §4 are valid for the original formulation, the inequalities added at any node of the tree are also valid for all other nodes. Hence, whenever the LP bound is evaluated at a given node of the tree, the linear program incorporates all cuts generated so far. To control the branch-and-cut process, additional aggre-

gate variables yki = ∑ j∈N x

k ij� i ∈ P, are introduced in the

formulation to reflect the assignment of requests to vehi- cles. At a given node of the search tree, if the yki variables are all integer in the current relaxation but there is at least one fractional xkij variable, the separation heuristics are exe- cuted in the hope of identifying violated valid inequalities (the heuristics are also executed at the root node). If at least one of the heuristics succeeds in finding one or more vio- lated inequalities, the relaxation is solved with all identified

cuts and the heuristics are executed again. The cut genera- tion process at a node terminates when all heuristics fail to find any violated inequality. If the solution to the relaxation is still fractional after the generation of cuts, branching is performed on a fractional yki variable, if there is any, or on a fractional xkij variable, otherwise. The variable selected for branching is that whose value is the farthest from the nearest integer. After a node has been processed, the next node selected for cutting and branching is that with the best bound. In addition to the application of the separation heuris-

tics at some nodes of the branch-and-bound tree, a pool of inequalities is checked exhaustively for violations at each node of the tree, including those where not all yki variables take integer values. The inequalities that belong to this pool are described in §5.2. Finally, prior to solving the problem by branch and cut, an upper bound is computed by using the tabu search heuristic of Cordeau and Laporte (2003b). This upper bound is then used to prune the enumeration tree whenever the solution value at a given node exceeds that of the upper bound.

5.1. Preprocessing

This section describes the time-window tightening, arc elimination, and variable-fixing steps that are performed prior to applying the branch-and-cut algorithm.

5.1.1. Time-Window Tightening. Let T denote the end of the planning horizon. In the case of an outbound user, the time window at the origin node can be tightened by setting ei = max 0�en+i − L − di� and li = min ln+i − ti�n+i − di�T�. In the case of an inbound user, the time window at the destination node can be tightened by set- ting en+i = max 0�ei +di +ti�n+i� and ln+i = min li +di + L�T�. The time window on nodes 0 and 2n + 1 can also be

tightened by setting e0 = e2n+1 = mini∈P∪D ei−t0i� and l0 = l2n+1 = maxi∈P∪D li +di + ti�2n+1�. 5.1.2. Arc Elimination. Formulation (1)–(14) is

defined on a complete graph G. However, because of time windows, pairing, and ride-time constraints, several arcs can in fact be removed from the graph as they cannot belong to a feasible solution. A simple analysis leads to the following observations: • arcs �0�n+i�, �i�2n+1�, and �n+i�i� are infeasible

for i ∈ P; • arc �i�j� with i�j ∈ N is infeasible if ei +di +tij >lj; • arcs �i�j� and �j�n + i� with i ∈ P, j ∈ N are both

infeasible if tij +dj + tj�n+i >L. As first proposed by Dumas et al. (1991), combining

time windows and pairing constraints leads to even stronger elimination rules: • arc �i�n+j� is infeasible if path � = j�i�n+j�n+i�

is infeasible; • arc �n+i�j� is infeasible if path � = i�n+i�j�n+j�

is infeasible;

• arc �i�j� is infeasible if paths �1 = i�j�n+ i�n+j� and �2 = i�j�n+j�n+ i� are both infeasible;

• arc �n+i�n+j� is infeasible if paths �1 = i�j�n+i� n+j� and �2 = j�i�n+ i�n+j� are both infeasible. In the presence of ride-time and capacity constraints, fur-

ther elimination can be performed by identifying pairs of users that are incompatible (i.e., that cannot be assigned to the same vehicle) because of the interaction between time- window, capacity, and ride-time constraints. Incompatible user pairs i�j� can be identified by checking the feasibility of the following paths: i�j�n+i�n+j�, i�j�n+j�n+i�, j�i�n+i�n+j�, j�i�n+j�n+i�, i�n+i�j�n+j�, and j�n+j�i�n+ i�. If none of these six paths is feasible, then all eight arcs

between i�n+i� and j�n+j� can be eliminated. Because of ride-time constraints, checking the feasibility of a path is not always straightforward. In the case of the path i�j� n + i�n + j�, for example, setting Bi = ei may lead to the violation of the ride-time constraint for user i in the event that unnecessary waiting time then occurs at node j. For this reason, the forward time slack should be computed at node i so as to delay the beginning of service as much as possible without violating any of the time windows. The same can be said about node j where the beginning of service should be delayed as much as possible by taking into account the additional constraint that the maximum ride time for user i should not be exceeded. In general, the forward time slack Fi at node i in a path i�i + 1��q� can be computed as follows:

Fi = min i�j�q

{ ∑ i<p�j

Wp +min lj −Bj�L−Pj� } � (42)

where Wi denotes the waiting time at node i, Pi denotes the ride time of the user whose destination node is i if n + 1 � i � 2n, and Pi = −�, otherwise. This definition of the forward time slack generalizes that of Savelsbergh (1992) for the TSP with time windows.

5.1.3. Variable Fixing. The identification of incom- patible user pairs can also be used to permanently assign users to specific vehicles. If the fleet of vehicles is homo- geneous, one can create a graph G′ = �N ′�E′�, where N ′ = 1��n� and E′ contains an edge �i�j� if i and j are incompatible users. Given a clique in G′, each user in the clique can be assigned to a different vehicle. In addition, if user i is assigned to vehicle k, constraints can be added to the formulation so as to forbid the assignment to vehi- cle k of users that are incompatible with i. Finding a clique of maximum cardinality in G′ may be very time consum- ing when n is large. In our implementation, we thus use a greedy heuristic described by Johnson (1974).

5.2. Initial Pool of Inequalities

This section describes the inequalities that are enumerated exhaustively and whose violations are checked individu- ally at every node of the branch-and-bound tree. These

inequalities are chosen to ensure that the size of the pool remains reasonable and the time spent processing a node does not overly increase. Bounds on time and load variables. For each node i ∈

P ∪D, the four bounding inequalities described in §4.1 are generated and added to the pool. Subtour elimination constraints. For each pair of users

i�j ∈ P, inequality (29) yields the following three particular cases: • S = i�j� ⇒ xij +xji +xn+i�j +xn+j�i � 1, • S = i�n+j� ⇒ xi�n+j +xn+j�i +xn+i�n+j � 1, • S = i�n+ i�j� ⇒ xij +xji +xi�n+i +xj�n+i +xn+i�j +

xn+j�i +xn+j�n+i � 2. Similarly, inequality (30) yields the following cases: • S = n + i�n + j� ⇒ xn+i�n+j + xn+j�n+i + xn+i�j +

xn+j�i � 1, • S = i�n+j� ⇒ xi�n+j +xn+j�i +xi�j � 1, • S = i�n + i�n + j� ⇒ xi�n+j + xn+j�i + xi�n+i +

xn+j�n+i +xn+i�j+j +xij +xn+i�j � 2, while inequalities (33) and (34) provide the following two cases: • S = n+ i�j�n+ i� ⇒ xn+i�j +2xj�n+i +xji +xi�n+i +

xn+j�n+i � 2� • S = i�n + i�n + j� ⇒ xi�n+i + xn+i�n+j + xn+j�i +

2xi�n+j +xij � 2. Finally, the four-node subtour inequality x�S� � 3 is gen- erated for every user pair i�j ∈ P yielding the subset S = i�j�n+ i�n+j�. Precedence constraints. For every pair of users i�j ∈ P,

the following four particular cases of precedence con- straints are generated: • S = 0�i�n+j� ⇒ x0i +xi�n+j +xn+j�i � 1, • S = i�n+j�2n+1� ⇒ xi�n+j +xn+j�i +xn+j�2n+1 � 1, • S = 0�i�n+i�n+j� ⇒ x0i +xi�n+i +xi�n+j +xn+j�i +

xn+i�n+j +xn+j�n+i � 2, • S = i�j�n+j�2n+1� ⇒ xij +xji +xi�n+j +xn+j�i +

xj�n+j +xn+j�2n+1 � 2. Generalized order constraints. For every pair of requests

i�j ∈ P, the following inequality is generated:

xi�n+j +xn+j�i +xn+i�j +xj�n+i � 1�

Infeasible path constraints. Finally, the following in- equality is generated for every pair of requests i�j ∈ P such that tij +dj + tj�n+jdn+j + tn+j�n+i >L:

xij +xj�n+j +xn+j�n+i � 1�

In addition, for each pair of incompatible users i�j ∈ P and every node l ∈ P ∪ D, the following inequalities are generated: • xil +xli +xlj +xjl � 1, • xil +xli +xl�n+j +xn+j�l � 1, • xn+i�l +xl�n+i +xlj +xjl � 1, • xn+i�l +xl�n+i +xl�n+j +xn+j�l � 1.

5.3. Separation Heuristics

We now describe the separation heuristics used to identify additional violated inequalities. At the root node as well as when all yki variables are integer but there is at least one fractional xkij variable at a given node, the following heuristics are executed sequentially.

5.3.1. Subtour Elimination Constraints. The identi- fication of violated inequalities of the form x�S� � �S� −1 can be achieved by solving a series of maximum-flow problems between any node i and all other nodes j ∈ N\ i�. However, in addition to being time consuming, this approach does not take the possible liftings into account. For these reasons, we resort to a simple tabu search heuris- tic inspired from that proposed by Augerat et al. (1999). Using the fact that 2x�S� + x�"�S�� = 2�S� in a feasi-

ble integer solution, violations of (29) can be identified by finding node sets S such that

x�"�S��−2 ∑ i∈S̄∩$�S�

∑ j∈S xij −2

∑ i∈S̄\$�S�

∑ j∈S∩$�S�

xij <2� (43)

The heuristic starts with an empty set S. At each iteration, it either adds or removes a node from the set S so as to minimize the left-hand side of (43). Whenever a node is removed from set S, its reinsertion is declared tabu for 1 iterations. The heuristic runs for a preset number of iter- ations (25 in our implementation) and may identify sev- eral violated inequalities during a single execution. At each iteration, the current set S is also checked for possible vio- lations of inequality (33). To this purpose, the node with the largest outgoing flow is numbered as i1, and the other nodes are numbered at random. A similar heuristic is used to identify violations of

inequalities (30) and (34). In the latter case, the node with the largest incoming flow is numbered as i1.

5.3.2. Capacity Constraints. Again, we use a tabu search heuristic to identify sets S such that q�S�>Q and x�"�S��<4. The heuristic starts either with a random sub- set S ⊆ P or with a random subset S ⊆ D. At each itera- tion, a node is either removed or added to the set S so as to minimize the value of x�"�S��, with the constraint that q�S� > Q. The heuristic runs for 25 iterations and multi- ple violated inequalities may be identified during a single execution.

5.3.3. GeneralizedOrderConstraints. Ruland(1995) proposed an approach to identify violated generalized order constraints with m = 2. This approach requires the com- putation of O��N�2� maximum flows. Here, we use instead three simple heuristics, the second and third of which take advantage of the lifted forms (36) and (37). The first heuristic attempts to identify violations of

inequalities (35) for the special case where m = 2 and �U1� = �U2� = 3. For each pair of users i�j ∈ P, we first form the subsets U1 = i�n + j� and U2 = j�n + i�. We

then identify nodes k1 and k2 such that x�U1� and x�U2� are maximized, and check for a violation of the resulting inequality. The second and third heuristics are aimed at finding vio-

lations of inequalities (36) and (37) for the special case where m = 3 and �U1� = �U2� = �U3� = 2. The second heuristic finds, for each user i, a user j that maximizes xi�n+j + xn+j�i + xij. It then finds a user k such that the left-hand side of (36) is maximized. Similarly, the third heuristic finds, for each user i, a user j that maximizes xi�n+j +xn+j�i +xn+i�n+j, and then finds a user k that max- imizes the left-hand side of (37).

5.3.4. Infeasible Path Inequalities. Violated path in- equalities are identified by means of a path-construction heuristic applied to each user i ∈ P. The heuristic starts with node i and gradually constructs a path �i = i�k1� k2�� by iteratively moving from the current node kj to the node kl that maximizes the value of xkjkl. The process stops if a cycle is formed or if the heuristic reaches either node n+ i or node 2n+1. If the path stops at node n+ i, it is checked for a violation of inequality (41).

6. Computational Experiments The branch-and-cut algorithm was implemented in C++ by using ILOG Concert 1.3 and CPLEX 8.1. It was run on a 2.5 GHz Pentium 4 computer with 512 Mb of memory. The algorithm was applied to two sets of randomly gen-

erated instances comprising up to 48 users. In an instance with n users, where n is even, users 1��n/2 are assumed to formulate outbound requests while users n/2+1��n formulate inbound requests. In all instances, the coordinates of pick-up and drop-off nodes are randomly and indepen- dently chosen in the square �−10�10�× �−10�10� accord- ing to a uniform distribution, and the depot is located at the center of this square. For every arc �vi�vj� ∈ A, the rout- ing cost cij and travel time tij are equal to the Euclidean distance between the two nodes.

Table 1. Characteristics of the first set of instances.

Without reductions With reductions Upper

Instance �K� n T Q L Cons Vars LP Cons Vars LP bound a2-16 2 16 480 3 30 1�289 1�211 219.60 893 665 269.54 294.25 a2-20 2 20 600 3 30 1�858 1�767 251.06 1�392 1�153 301.66 344.83 a2-24 2 24 720 3 30 2�714 2�585 316.38 1�997 1�628 364.22 431.12 a3-18 3 18 360 3 30 1�773 2�281 193.05 1�296 1�258 256.28 300.48 a3-24 3 24 480 3 30 2�898 3�880 221.63 2�244 2�446 260.56 344.83 a3-30 3 30 600 3 30 4�369 5�938 338.36 2�998 3�282 383.10 496.52 a3-36 3 36 720 3 30 6�077 8�329 360.84 3�919 4�171 484.57 600.75 a4-16 4 16 240 3 30 1�657 2�364 175.61 1�193 1�173 207.41 282.68 a4-24 4 24 360 3 30 3�279 5�181 235.68 2�190 2�616 293.07 375.07 a4-32 4 32 480 3 30 5�225 8�977 311.20 3�828 4�941 361.01 486.57 a4-40 4 40 600 3 30 8�104 13�829 366.04 5�535 8�200 408.13 571.07 a4-48 4 48 720 3 30 10�884 19�745 366.04 8�171 12�824 428.60 680.99

Avg. 279.62 334.85 434.10

A time window �ei�li� is also associated with each node. For an outbound user i, a time window was generated by first choosing a number ln+i in the interval �60�T� and then setting en+i = ln+i − 15, where T denotes the length of the planning horizon. For an inbound user, the value of ei was chosen in the interval �0�T − 60� and the value of li was set equal to ei +15. Time windows on the origin nodes of outbound requests and on the destination nodes of inbound requests are set as explained in §5.1.1. In the first set of instances, Q = 3 for every vehicle,

qi = 1 and di = 3 for every user, and the maximum ride time L is equal to 30 minutes. In the second set, Q = 6 for every vehicle and the values of qi are chosen randomly according to a uniform distribution on the set 1��Q�. The maximum ride time is set equal to 45 minutes. In this case, we assume that service time is proportional to the number of passengers and di = qi. The maximum route duration is equal to T . The first set of instances represents the context where cars are used for the transportation of individuals whereas the second set reflects the situation of a transporter using mini-busses for the transportation of indi- viduals or groups of individuals. All instances used in computational experiments are

available on the website http://www.hec.ca/chairedistributique/data/darp.

Tables 1 and 2 provide the characteristics of these instances, the number of constraints and variables in the resulting mixed-integer programming formulation, and the value of the LP relaxation. The corresponding values are computed with and without the application of the pre- processing techniques described in §5.1. In both cases, the problem reduction step of CPLEX is nevertheless per- formed. These results show that despite the strength of the CPLEX preprocessing routines, a considerable reduction in problem size (and improvement in the LP bound at the root node) is achieved through the application of the specialized arc elimination and variable-fixing techniques of §5.1. The improvement in the LP lower bound is 19.75% for the first set of instances and 15.17% for the second set. We also

Table 2. Characteristics of the second set of instances.

Without reductions With reductions Upper

Instance �K� n T Q L Cons Vars LP Cons Vars LP bound b2-16 2 16 480 6 45 1�030 993 218.80 834 653 250.55 309.61 b2-20 2 20 600 6 45 1�175 1�241 290.33 768 612 314.05 334.93 b2-24 2 24 720 6 45 1�922 1�873 311.71 1�436 1�189 372.80 445.11 b3-18 3 18 360 6 45 1�334 1�754 218.26 1�068 1�074 256.15 301.80 b3-24 3 24 480 6 45 2�230 2�915 284.78 1�625 1�669 329.14 394.57 b3-30 3 30 600 6 45 3�410 4�457 344.68 2�383 2�626 448.31 536.04 b3-36 3 36 720 6 45 4�200 6�035 437.37 2�984 3�403 493.79 611.79 b4-16 4 16 240 6 45 1�233 1�771 209.25 919 943 223.11 299.07 b4-24 4 24 360 6 45 2�426 4�129 252.87 1�850 2�262 312.28 380.27 b4-32 4 32 480 6 45 3�928 6�405 341.76 2�778 3�719 409.25 500.92 b4-40 4 40 600 6 45 5�361 9�481 477.90 3�797 5�368 525.28 662.91 b4-48 4 48 720 6 45 8�854 15�611 496.37 6�404 9�608 538.43 685.46

Avg. 323.67 372.76 455.21

indicate in the last columns of Tables 1 and 2 the value of an upper bound obtained by running the tabu search heuris- tic of Cordeau and Laporte (2003b) for 1,000 iterations. In Tables 3 and 4, we indicate the lower bound obtained

at the root node (after applying the problem reduction steps) under seven different scenarios. Column LP indicates the bound obtained without any kind of cut generation. The next four columns indicate the bound obtained by generat- ing violated inequalities of one of the following families: capacity constraints (CC), generalized order constraints (GOC), infeasible path constraints (IPC) and subtour elim- ination constraint (SEC). Column CPX then reports the bound obtained when activating the automatic cut genera- tion of CPLEX (but not the generation of user cuts). In this case, several types of inequalities (e.g., implied bound cuts, flow cuts, Gomory fractional cuts) are automatically gener- ated by CPLEX at the root node. Finally, column Full indi- cates the bound obtained when applying both the automatic cut generation of CPLEX and the four separation heuristics of §5.3 (as well as the inequality pool of §5.2). The last col- umn shows the improvement (in percentage) in the lower bound obtained when generating user cuts together with the

Table 3. Initial lower bounds—first set of instances.

Instance LP CC GOC IPC SEC CPX Full Imp. (%)

a2-16 269.54 269.54 270.54 278.88 269.66 287.38 290.18 0.97 a2-20 301.66 302.21 304.30 310.73 308.51 321.70 326.13 1.38 a2-24 364.22 365.55 373.30 383.94 365.45 377.20 395.57 4.87 a3-18 256.28 256.28 258.71 260.08 262.43 259.94 269.15 3.54 a3-24 260.56 260.56 260.97 268.45 261.98 276.98 291.00 5.06 a3-30 383.10 386.94 383.26 409.99 384.31 386.62 421.60 9.05 a3-36 484.57 484.57 484.91 516.69 484.87 502.74 529.47 5.32 a4-16 207.41 209.00 207.41 214.45 209.13 212.48 229.19 7.86 a4-24 293.07 293.08 299.58 303.84 293.07 299.52 318.93 6.48 a4-32 361.01 361.01 361.39 368.13 361.24 375.18 390.02 3.96 a4-40 408.13 418.23 412.42 431.30 414.43 421.52 450.43 6.86 a4-48 428.60 431.33 428.61 437.82 429.57 432.09 457.25 5.82

Avg. 334.85 336.53 337.12 348.69 337.05 346.11 364.08 5.10

CPLEX cuts instead of relying only on the latter class (i.e., relative difference between the last two scenarios). These results show that using any of the four types

of inequalities yields a small improvement in the average lower bound obtained. For both sets of instances, the largest improvement is obtained with the generation of infeasible path inequalities, while subtour elimination, capacity, and generalized order constraints have a limited impact. Simi- lar improvements are also obtained by enabling the auto- matic cut generation of CPLEX. Nevertheless, combining these cuts with our own valid inequalities yields significant improvements. The average improvement exceeds 5% for both sets of instances. It is worth mentioning that similar conclusions are reached when comparing the best bound obtained after running the branch-and-cut algorithm for a certain number of nodes. Finally, Tables 5 and 6 report the results obtained by

using the branch-and-cut algorithm of CPLEX alone and with the concurrent generation of all families of inequali- ties. In both cases, a maximum of four hours of CPU time were allowed for the solution of each instance. We indicate the best bound obtained, the CPU time (in minutes) used,

Table 4. Initial lower bounds—second set of instances.

Instance LP CC GOC IPC SEC CPX Full Imp. (%)

b2-16 250.55 250.55 264.13 252.67 251.35 268.08 289.66 8.05 b2-20 314.05 321.24 314.05 314.05 315.15 327.46 332.61 1.57 b2-24 372.80 377.07 380.54 377.29 380.02 389.81 404.80 3.85 b3-18 256.15 256.15 260.86 260.15 256.15 261.82 268.67 2.62 b3-24 329.14 329.14 331.34 338.37 329.38 332.98 344.85 3.56 b3-30 448.31 448.93 453.21 482.79 448.31 465.74 495.56 6.40 b3-36 493.79 495.04 506.92 502.21 493.79 520.59 547.83 5.23 b4-16 223.11 223.11 230.35 228.07 224.17 225.86 247.97 9.79 b4-24 312.28 312.28 315.36 320.69 312.28 320.03 336.40 5.12 b4-32 409.25 409.25 413.98 420.59 412.90 412.77 439.98 6.59 b4-40 525.28 525.28 531.44 544.31 525.28 531.17 561.36 5.68 b4-48 538.43 538.43 541.94 553.12 539.05 547.85 567.06 3.51

Avg. 372.76 373.87 378.68 382.86 373.99 383.68 403.06 5.16

and the number of nodes explored in the branch-and-bound tree. For the branch-and-cut algorithm with user cut gen- eration, we also indicate the total number of cuts gener- ated. An asterisk preceding a lower bound indicates that the instance was solved to optimality. One can observe that the generation of user cuts sig-

nificantly improves the performance of the branch-and-cut algorithm. For the 14 instances solved to optimality by both approaches within the time limit, the average CPU time was 25.93 minutes for CPLEX alone compared with 1.81 minutes for CPLEX with user cut generation. The average number of nodes explored went from 95,466 to 3,248. For the five instances that could not be solved optimally by any of the two approaches, the average lower bound reached by CPLEX alone was 500.91 compared with 517.48 for the branch-and-cut with user cut generation. Finally, the latter algorithm solved five more instances to optimality. The small number of cuts generated with respect to the

total number of nodes explored in the branch-and-bound tree is explained by the fact that the separation procedures are executed only when all yki variables are integer, but at least one xkij variable is fractional. For most instances, this occurs at approximately 5% of the nodes. Experiments performed with the generation of cuts at every node of

Table 5. Comparisons with CPLEX—first set of instances.

CPLEX B and C B and C with user cuts

Instance Bound CPU Nodes Bound CPU Nodes Cuts

a2-16 ∗294�25 0�01 23 ∗294�25 0�01 4 46 a2-20 ∗344�83 0�05 313 ∗344�83 0�08 181 93 a2-24 ∗431�12 1�42 10�868 ∗431�12 0�27 444 140 a3-18 ∗300�48 0�41 3�596 ∗300�48 0�20 565 172 a3-24 ∗344�83 76�59 310�667 ∗344�83 4�17 6�863 324 a3-30 472�17 240�00 515�931 ∗494�85 42�39 43�632 425 a3-36 570�26 240�00 504�553 ∗583�19 12�46 7�451 423 a4-16 ∗282�68 21�49 145�680 ∗282�68 2�54 6�615 404 a4-24 359�52 240�00 442�000 ∗375�02 25�19 28�939 351 a4-32 427�65 240�00 189�900 447�66 240�00 105�500 806 a4-40 462�21 240�00 65�000 475�05 240�00 32�200 716 a4-48 466�70 240�00 40�400 486�03 240�00 21�400 1�019

the tree have on average produced slightly worse results because of the extra time spent by the separation proce- dures. It is likely that a larger number of cuts could be generated by using exact or more sophisticated heuristic separation procedures. The development of such proce- dures is by itself an important area of research that will be addressed in future work. All instances used for testing are symmetric (i.e., cij = cji

for every pair of nodes i, j). In additional experiments, we have also tried solving asymmetric instances. Unlike the asymmetric TSP, which is usually easier to solve than its symmetric counterpart, it appears that the DARP is unaffected by the asymmetry of the cost matrix. This is explained by the fact that a while a TSP tour is feasible in both directions, thus leading to highly fractional solutions, a DARP route is feasible in only one direction because of precedence constraints. The computational results presented in this section show

that instances with up to 30 users can be solved optimally in reasonable CPU times. Using column generation to address the VRPPD with time windows, Savelsbergh and Sol (1998) were also able to consistently solve instances with up to 30 requests in short computation times. Experiments performed with the branch-and-cut algorithm on similar

Table 6. Comparisons with CPLEX—second set of instances.

CPLEX B and C B and C with user cuts

Instance Bound CPU Nodes Bound CPU Nodes Cuts

b2-16 ∗309�41 0�21 5�815 ∗309�41 0�15 487 130 b2-20 ∗332�64 0�01 26 ∗332�64 0�01 2 30 b2-24 ∗444�71 2�76 32�399 ∗444�71 0�13 146 103 b3-18 ∗301�64 1�29 12�223 ∗301�64 0�70 2�458 254 b3-24 ∗394�51 7�27 42�950 ∗394�51 3�62 8�147 252 b3-30 ∗531�44 189�74 574�281 ∗531�44 6�82 9�862 274 b3-36 588�44 240�00 447�474 ∗603�79 62�07 66�079 291 b4-16 ∗296�96 2�44 20�189 ∗296�96 0�79 2�579 196 b4-24 ∗371�41 59�31 175�495 ∗371�41 5�85 7�119 255 b4-32 460�78 240�00 222�600 ∗494�82 176�82 126�332 375 b4-40 570�37 240�00 153�400 591�76 240�00 81�100 632 b4-48 577�64 240�00 50�000 586�91 240�00 21�000 753

instances (obtained by relaxing the ride-time constraints) have, however, shown a deterioration of the algorithm’s per- formance. This is partly explained by the fact that infeasible path inequalities as well as most preprocessing techniques then become irrelevant. In addition, column generation usu- ally provides tighter formulations, leading to stronger LP relaxation values. Nevertheless, it seems that handling ride- time constraints in a column generation approach is far from trivial. Indeed, the states of the dynamic program- ming algorithm used for pricing would have to take into account the individual ride time of each user. Label elimi- nation based on dominance is then likely to be very weak.

7. Conclusion We have introduced several new valid inequalities and a branch-and-cut algorithm for the dial-a-ride problem. Although no polyhedral analysis was carried out, the use- fulness of these inequalities has been demonstrated through computational experiments. A comparison with version 8.1 of CPLEX shows that the branch-and-cut algorithm pro- posed here reduces both the CPU time and the number of nodes explored in the branch-and-bound tree. The method- ology proposed in this paper obviously cannot be used to solve large-scale instances containing hundreds or thou- sands of users, as is sometimes the case in large cities. It is, however, fast enough to be used as a postprocessor to optimize individual routes or small subsets of routes pro- duced by a metaheuristic. Future work will concentrate on the development of more refined separation procedures for the various types of inequalities introduced in this paper.

Appendix Details of the Proof of Proposition 3. There are two cases to distinguish. Case 1. An arc �i1�il� with 2 � l � m−1 is part of the

solution. First, we show that there is at least one subset Uk

with 1 � k � l − 1 for which x�Uk� � �Uk� − 2. Suppose this is not true. Then, for every k, there is a path in Uk

connecting ik and n + ik+1. Because the arc �i1�il� is in the solution, the path for U1 must visit node n+ i2 before visiting node i1. In addition, this path must appear imme- diately before the arc �i1�il� in the solution. Now, for k = 2��l − 1, the path for Uk must appear before that for Uk−1 in the solution to avoid violating the precedence con- straint for user ik. However, for k = l−1, this implies that node n+il appears before node il, which is a contradiction. Next, we show that there is at least one subset Uk with

l � k � m for which x�Uk� � �Uk� −2. Suppose this is not true. Then, for every k, there is a path in Uk connecting ik and n+ ik+1. Because the arc �i1�il� is in the solution, the path for Ul must visit node il before visiting node n+ il+1. In addition, this path must appear immediately after the arc �i1�il� in the solution. Now, for k = l +1��m, the path for Uk must appear before that for Uk−1 in the solution to avoid violating the precedence constraint for user ik. How- ever, for k = m, this implies that node n+i1 appears before node i1, which is a contradiction. Case 2. An arc �i1�n+ il� with 3 � l � m is part of the

solution. First, we show that there is at least one subset Uk with

1 � k � l − 1 for which x�Uk� � �Uk� − 2. Suppose this is not true. Because the arc �i1�n+ il� is in the solution, the path for U1 must visit node n+ i2 before visiting node i1. In addition, this path must appear immediately before the arc �i1�n + il� in the solution. Now, for k = 2��l − 1, the path for Uk must appear before that for Uk−1 in the solution to avoid violating the precedence constraint for user ik. However, for k = l−1, this implies that node n+il appears before node i1, which is a contradiction because the arc �i1�n+ il� is part of the solution. Next, we show that for there is at least one subset Uk

with l � k � m for which x�Uk� � �Uk�−2. Suppose this is not true. Because the arc �i1�n+il� is in the solution, then for k = l��m, the path for Uk must appear before the arc �i1�n+il� in the solution to avoid violating the precedence constraint for user ik. However, for k = m, this implies that the path connecting im and n + i1 appears before node i1, which is a contradiction. �

Details of the Proof of Proposition 4. Again, there are two cases to distinguish. Case 1. An arc �n+ i1�il� with 2 � l � m−2 is part of

the solution. For k = l−1��1, the path connecting ik and n+ ik+1

must appear after the arc �n+i1�il� in the solution to avoid violating the precedence constraint for user ik+1. However, for k = 1 this implies that the path connecting i1 and n+i2 appears after the arc �n+ i1�il�, which is a contradiction. Because the arc �n + i1�il� is in the solution, the sub-

path for Ul must visit node il before visiting node n+ il+1. For k = l + 1��m, the path for Uk must appear before that for Uk−1 to avoid violating the precedence constraint for user ik. However, for k = m, this implies that the path connecting im and n+i1 appears before the path for at least one other subset Uh with l+1 � h � m−1, and thus node n+ i1 appears in two different places in the solution. Case 2. An arc �n+i1�n+il� with 2 � l � m−1 is part

of the solution. For k = 1��l−1, the path connecting ik and n+ ik+1

must appear after the arc �n+i1�n+il�, but for k = 1, this implies that node i1 appears after node n + i1, which is a contradiction. The path for Ul must appear before the arc �n+i1�n+il�

and, for k = l +1��m−1, the path for Uk must appear before that for Uk−1. However, for k = m, this implies that the path connecting im and n + i1 appears before the path for at least one other subset Uh with l � h � m − 1, and thus node n + i1 appears in two different places in the solution. �

Acknowledgments This work was supported by the Natural Sciences and Engineering Research Council of Canada under grant 227837-00. This support is gratefully acknowledged. Thanks are due to Luigi Moccia for comments on an ear- lier version of this paper. The author is also thankful to two anonymous referees.

References Ascheuer, N. 1996. Hamiltonian path problems in the on-line optimiza-

tion of flexible manufacturing systems. Ph.D. thesis, Konrad-Zuse- Zentrum für Informationstechnik Berlin, Berlin, Germany.

Ascheuer, N., M. Fischetti, M. Grötschel. 2001. Solving the asymmetric travelling salesman problem with time windows by branch-and-cut. Math. Programming 90 475–506.

Ascheuer, N., M. Jünger, G. Reinelt. 2000. A branch & cut algorithm for the asymmetric traveling salesman problem with precedence con- straints. Comput. Optim. Appl. 17 61–84.

Augerat, P., J. M. Belenguer, E. Benavent, A. Corberán, D. Naddef. 1999. Separating capacity inequalities in the CVRP using tabu search. Eur. J. Oper. Res. 106 546–557.

Balas, E., M. Fischetti, W. R. Pulleyblank. 1995. The precedence-con- strained asymmetric traveling salesman polytope. Math. Programming 68 241–265.

Bodin, L. D., T. Sexton. 1986. The multi-vehicle subscriber dial-a-ride problem. TIMS Stud. Management Sci. 26 73–86.

Cordeau, J.-F., G. Laporte. 2003a. The dial-a-ride problem (DARP): Variants, modeling issues and algorithms. 4OR—Quart. J. Belgian, French Italian Oper. Res. Soc. 1 89–101.

Cordeau, J.-F., G. Laporte. 2003b. A tabu search heuristic for the static multi-vehicle dial-a-ride problem. Transportation Res. B 37 579–594.

Desaulniers, G., J. Desrosiers, A. Erdmann, M. M. Solomon, F. Soumis. 2002. VRP with pickup and delivery. P. Toth, D. Vigo, eds. The Vehi- cle Routing Problem. SIAM Monographs on Discrete Mathematics and Applications. Philadelphia, PA, 225–242.

Desrochers, M., G. Laporte. 1991. Improvements and extensions to the Miller-Tucker-Zemlin subtour elimination constraints. Oper. Res. Lett. 10 27–36.

Desrosiers, J., Y. Dumas, F. Soumis. 1986. A dynamic programming solu- tion of the large-scale single-vehicle dial-a-ride problem with time windows. Amer. J. Math. Management Sci. 6 301–325.

Dumas, Y., J. Desrosiers, F. Soumis. 1989. Large scale multi-vehicle dial- a-ride problems. Technical Report G-89-30, GERAD, HEC Montréal, Montréal, Quebec, Canada.

Dumas, Y., J. Desrosiers, F. Soumis. 1991. The pickup and delivery prob- lem with time windows. Eur. J. Oper. Res. 54 7–22.

Gendreau, M., G. Laporte, F. Semet. 1998. A branch-and-cut algorithm for the undirected selective traveling salesman problem. Networks 32 263–273.

Grötschel, M., M. W. Padberg. 1985. Polyhedral theory. E. L. Lawler, J. K. Lenstra, A. H. G. Rinnooy Kan, D. B. Shmoys, eds. The Trav- eling Salesman Problem. Wiley, New York, 251–305.

Jaw, J., A. R. Odoni, H. N. Psaraftis, N. H. M. Wilson. 1986. A heuristic algorithm for the multi-vehicle advance-request dial-a-ride problem with time windows. Transportation Res. B 20 243–257.

Johnson, D. S. 1974. Approximation algorithms for combinatorial opti- mization problems. J. Comput. System Sci. 9 256–278.

Madsen, O. B. G., H. F. Ravn, J. M. Rygaard. 1995. A heuristic algorithm for a dial-a-ride problem with time windows, multiple capacities, and multiple objectives. Ann. Oper. Res. 60 193–208.

Miller, C. E., A. W. Tucker, R. A. Zemlin. 1960. Integer programming for- mulations and traveling salesman problems. J. Association Comput. Machinery 7 326–329.

Naddef, D., G. Rinaldi. 2002. Branch-and-cut algorithms for the capac- itated VRP. P. Toth, D. Vigo, eds. The Vehicle Routing Prob- lem. SIAM Monographs on Discrete Mathematics and Applications. Philadelphia, PA, 53–84.

Psaraftis, H. N. 1980. A dynamic programming approach to the single- vehicle, many-to-many immediate request dial-a-ride problem. Trans- portation Sci. 14 130–154.

Psaraftis, H. N. 1983. An exact algorithm for the single-vehicle many-to- many dial-a-ride problem with time windows. Transportation Sci. 17 351–357.

Ruland, K. S. 1995. Polyhedral solution to the pickup and delivery problem. Ph.D. thesis, Sever Institute of Technology, Washington University, St. Louis, MO.

Ruland, K. S., E. Y. Rodin. 1997. The pickup and delivery problem: Faces and branch-and-cut algorithm. Comput. Math. Appl. 33 1–13.

Savelsbergh, M. W. P. 1985. Local search in routing problems with time windows. Ann. Oper. Res. 4 285–305.

Savelsbergh, M. W. P. 1992. The vehicle routing problem with time win- dows: Minimizing route duration. ORSA J. Comput. 4 146–154.

Savelsbergh, M. W. P., M. Sol. 1998. Drive: Dynamic routing of indepen- dent vehicles. Oper. Res. 46 474–490.

Toth, P., D. Vigo. 1996. Fast local search algorithms for the handi- capped persons transportation problem. I. H. Osman, J. P. Kelly, eds. Meta-Heuristics: Theory and Applications. Kluwer, Boston, MA, 677–690.

Toth, P., D. Vigo. 1997. Heuristic algorithms for the handicapped persons transportation problem. Transportation Sci. 31 60–71.