A Branch-and-Cut Algorithm for the
Dial-a-Ride Problem
JEAN-FRANC?OIS CORDEAU
Canada Research Chair in Distribution Management,HEC Montr′e al 3000,chemin de la C?o te-Sainte-Catherine
Montr′e al,Canada H3T2A7
October28,2003
Submitted to Operations Research
Abstract
In the dial-a-ride problem,users formulate requests for transportation from a speci?c origin to a speci?c 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 pickup and delivery and the vehicle routing https://www.sodocs.net/doc/e216072203.html,putational experiments performed on randomly generated instances show that the proposed approach can be used to solve small to medium size instances.
Subject classi?cations:Transportation:vehicle routing.Programming:cutting plane. Area of review:Transportation.
1Introduction
In the Dial-a-Ride Problem(DARP),users formulate requests for transportation from a speci?c origin(or pick-up point)to a speci?c destination(or drop-o?point).Transportation 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 presence of two con?icting objectives: minimizing operating costs and minimizing user inconvenience.Operating costs are mostly related to?eet size and distance traveled while user inconvenience is often measured in terms of deviations from desired pick-up and drop-o?times and in terms of excess ride time(i.e., the di?erence 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 speci?es 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 pre-speci?ed 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 with15-minute time windows and a maximum ride time of60minutes, a user wishing to arrive at destination at9h00would be picked-up no earlier than7h45and dropped-o?between8h45and9h00.It should be emphasized that imposing time windows is not su?cient to accurately impose a maximum ride time.In the above example,a pick-up at7h45and a drop-o?at9h would exceed the maximum ride time of60minutes.However, constraining the pick-up to take place no earlier than8h would be too restrictive because it would,in particular,prohibit a pick-up at7h45and a drop-o?at8h45.
Dial-a-ride services may operate according to a static or to a dynamic mode.In the?rst 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 since 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 di?cult problem that generalizes the Vehicle Routing Problem with Pick-up and Delivery (VRPPD).Finding a feasible solution for the DARP is itself NP-hard since it also generalizes 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.This algorithm uses adaptations of known valid inequalities for the VRP and VRPPD as well as new inequalities that take advantage of the special structure of the problem.
The remainder of the article is organized as follows.The next section brie?y reviews relevant work on the DARP and closely related problems.Section3de?nes the DARP formally and introduces a mixed-integer formulation.Section4introduces several families of valid inequalities used in the branch-and-cut algorithm which is then described in Section5,along with preprocessing https://www.sodocs.net/doc/e216072203.html,putational experiments are reported in Section6and the conclusion follows in Section7.
2Literature Review
Early research on the DARP was carried out by Psaraftis(1980,1983)who developed dynamic programming algorithms 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 to40users.Ruland(1995)and Ruland and Rodin(1997)proposed a branch-and-cut algorithm for a special case of the DARP 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-o?node in the vehicle route.These authors provided an integer programming formulation of the problem based on an undirected graph and described four families of valid inequalities:subtour elimination constraints,precedence constraints,generalized order constraints and order matching constraints.These inequalities are also valid for the more general problem studied in this paper and they will be discussed in more depth in Section4.
Most of the algorithms for the multiple-vehicle case are heuristics or meta-heuristics.An insertion method capable 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 heuristic for the dynamic case while Toth and Vigo(1996,1997)used local search and a tabu-thresholding procedure 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. 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).
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.
3Formulation
Let n denote the number of users(or requests)to be served.The DARP may be de?ned 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-o?nodes,respectively, while nodes0and2n+1represent 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 Q k and the total duration of its route cannot exceed T k. With each node i∈N are associated a load q i and a non-negative service duration d i such that q0=q2n+1=0,q i=?q n+i(i=1,...,n)and d0=d2n+1=0.A time window[e i,l i] is also associated with node i∈N where e i and l i represent the earliest and latest time, respectively,at which service may begin at node i.With each arc(i,j)∈A are associated a routing cost c ij and a travel time t ij.Finally,denote by L the maximum ride time of a user. For each arc(i,j)∈A and each vehicle k∈K,let x k ij=1if vehicle k travels from node i to node j.For each node i∈N and each vehicle k∈K,let B k i be the time at which vehicle k begins service at node i,and Q k i be the load of vehicle k after visiting node i.Finally,for each user i,let L k i be the ride time of user i on vehicle k.The formulation is as follows:
Min k∈K i∈N j∈N c k ij x k ij(1)
subject to
k∈K j∈N x k ij=1?i∈P(2)
j∈N x k ij? j∈N x k n+i,j=0?i∈P,k∈K(3)
j∈N x k0j=1?k∈K(4) j∈N x k ji? j∈N x k ij=0?i∈P∪D,k∈K(5)
i∈N x k i,2n+1=1?k∈K(6)
B k j≥(B k i+d i+t ij)x k ij?i∈N,j∈N,k∈K(7)
Q k j≥(Q k i+q j)x k ij?i∈N,j∈N,k∈K(8)
L k i=B k n+i?(B k i+d i)?i∈P,k∈K(9)
B k2n+1?B k0≤T k?k∈K(10)
e i≤B k i≤l i?i∈N,k∈K(11)
t i,n+i≤L k i≤L?i∈P,k∈K(12) max{0,q i}≤Q k i≤min{Q k,Q k+q i}?i∈N,k∈K(13)
x k ij∈{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 origin depot and ends at the destination depot.Consistence of the time and load variables is ensured by constraints(7)and(8).Equalities(9)de?ne the ride time of each user which is bounded by constraints(12).It is worth mentioning that the latter also act as precedence constraints because the non-negativity of the L k i 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 non-linear because of constraints(7)and(8).Introducing constants M k ij and W k ij,these constraints can,however,be linearized as follows:
B k j≥B k i+d i+t ij?M k ij(1?x k ij)?i∈N,j∈N,k∈K(15)
Q k j≥Q k i+q j?W k ij(1?x k ij)?i∈N,j∈N,k∈K.(16)
The validity of these constraints is ensured by setting M k ij≥max{0,l i+d i+t ij?e j}and W k ij≥min{Q k,Q k+q i}.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 B i at every node except the origin depot0and the destination depot2n+1. In this case,one replaces(7)and(9)with the constraints
B j≥(B k0+d0+t0j)x k0j?j∈N,k∈K(17)
B j≥(B i+d i+t ij) k∈K x k ij?i∈N,j∈N(18)
B k2n+1≥(B i+d i+t i,2n+1)x k i,2n+1?i∈N,k∈K(19)
L i=B n+i?(B i+d i)?i∈P,(20) to which the same linearization process can be applied.Along the same lines,if the?eet of vehicles is homogeneous in the sense that Q k=Q for every k∈K,one can replace(8)with
Q j≥(Q k0+q j)x k0j?j∈N,k∈K(21)
Q j≥(Q i+q j) k∈K x k ij?i∈N,j∈N(22) Q k2n+1≥(Q i+q2n+1)x k i,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:
Q j≥Q i+q j?W ij(1? k∈K x k ij)+(W ij?q i?q j) k∈K x k ji?i∈N,j∈N,k∈K.(24)
Observe that in our case an equivalent lifting cannot be performed with constraints(18) because waiting will sometimes take place after the beginning of a time window(i.e.,B i>e i) in order to reduce the ride time of a user.
4Valid 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 de?ne facets of the DARP polytope appears to be a rather challenging task.However,their usefulness 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,de?neδ(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 x ij= k∈K x k ij denote the total?ow on arc(i,j)and de?ne x(S)= i,j∈S x ij.Similarly,let x(A′)= (i,j)∈A′x ij for any arc set A′?A.
4.1Bounds on time and load variables
As?rst suggested by Desrochers and Laporte(1991)in the context of the TSP with time windows,bounds on the time variables B i can be strengthened as follows:
B i≥e i+ j∈N\{i}max{0,e j?e i+d j+t ij}x ji(25)
B i≤l i? j∈N\{i}max{0,l i?l j+d i+t ij}x ij.(26)
These inequalities were used,for example,for solving the Asymmetric Travelling Salesman Problem with Time Windows by branch-and-cut(Ascheuer et al.,2001).
Similarly,bounds on load variables Q i can also be strengthened as follows:
Q i≥max{0,q i}+ j∈N\{i}max{0,q j}x ji(27)
{q j}?q i)x0i? j∈N\{i}max{0,q j}x ij.(28) Q i≤min{Q,Q+q i}?(Q?max
j∈N\{i}
4.2Subtour elimination constraints
Consider the simple subtour elimination constraint x(S)≤|S|?1for S?P∪D.In the case of the DARP,this inequality can be lifted in two di?erent ways by taking into account the fact that for each user i,node i must be visited before node n+i.
Proposition1.Let S?P∪D and P+(S)={i|i∈S∩P and n+i∈S}.The following inequality is valid for the DARP:
x(S)+ i∈P+(S) j∈S x n+i,j≤|S|?1.(29)
Proof.Because of precedence relationships,set S must be entered at least once before using any of the lifted arcs.As a result,x(δ?(S))≥1+ i∈P+(S) j∈S x n+i,j.In addition, x(δ+(S))=x(δ?(S))implies that x(δ+(S))≥1+ i∈P+(S) j∈S x n+i,j.Since2x(S)+ x(δ+(S))+x(δ?(S))=2|S|,one obtains2x(S)+2+2 i∈P+(S) j∈S x n+i,j≤2|S|and,?nally,x(S)+ i∈P+(S) j∈S x n+i,j≤|S|?1.
Example.Consider the node set S={i,j}?P.The resulting lifted subtour elimination constraint is x ij+x ji+x n+i,j+x n+j,i≤1.This example is illustrated in Figure1.Observe that variables x n+i,i and x n+j,j need not be introduced in the lifted inequality since the
graph.
corresponding arcs can be trivially removed from the
Figure1:Lifted subtour elimination constraint for S={i,j}?P
Proposition2.Let S?P∪D and D?(S)={i|n+i∈S∩D and i∈S}.The following inequality is valid for the DARP:
x(S)+ i∈S j∈D?(S)x ij≤|S|?1.(30)
Proof.The proof is similar to that of Proposition1by observing that because of precedence relationships,set S must be exited at least once after using any of the lifted arcs and thus x(δ+(S))≥1+ i∈S j∈D?(S)x ij.
Example.Consider the node set S={n+i,n+j}?D.The resulting lifted subtour elimination constraint is x n+i,n+j+x n+j,n+i+x n+i,j+x n+j,i≤1.This is illustrated in 2.
Figure
In the case of a directed formulation,one can also lift subtour elimination constraints by taking into account the orientation of the arcs.For a set S={i1,i2,...,i h}?N with h≥3nodes,Gr¨o tschel and Padberg(1995)proposed the following inequalities for the asymmetric TSP:
h?1
j=1x i j,i j+1+x i h,i1+2h?1 j=2x i j,i1+h?1 j=3j?1 l=2x i j,i l≤h?1(31)
h?1
j=1x i j,i j+1+x i h,i1+2h j=3x i1,i j+h j=4j?1 l=3x i j,i l≤h?1.(32)
Of course,di?erent liftings are obtained by considering di?erent 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.
Proposition3.Let S={i1,i2,...,i h}?P∪D and P+(S)={i|i∈S∩P and n+i∈S}. The following inequality is valid for the DARP:
h?1
j=1x i j,i j+1+x i h,i1+2h?1 j=2x i j,i1+h?1 j=3j?1 l=2x i j,i l+ i p∈P+(S)x n+i p,i1≤h?1.(33)
Proof.Suppose that one arc of the form(n+i p,i1)with i p∈P+(S)is part of the solution. Then all arcs of the form(i j,i1)with2≤j≤h?1cannot 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.But because set S contains the origin node i p,this subpath together with the arc(n+i p,i1)would violate the precedence constraint for user i p. 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 x ij+x jk+x ki+2x ji+x n+j,i+
3.
x n+k,i≤2.This is illustrated in Figure
Figure3:Lifted directed subtour elimination constraint for S={i,j,k}?P
Proposition4.Let S={i1,i2,...,i h}?P∪D and D?(S)={i|n+i∈S∩D and i∈S}. The following inequality is valid for the DARP:
h?1
j=1x i j,i j+1+x i h,i1+2h j=3x i1,i j+h j=4j?1 l=3x i j,i l+ i p∈D?(S)x i1,i p≤h?1.(34)
Proof.The proof is similar to that of Proposition3by observing that if one arc of the form(i1,i p)with i p∈D?(S)is part of the solution,then all arcs of the form(i1,i j)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 elimination constraint(obtained with i1=n+i,i2=n+j,i3=n+k)is
4.
x n+i,n+j+x n+j,n+k+x n+k,n+i+2x n+i,n+k+x n+i,j+x n+i,k≤2.This is illustrated in Figure
Figure4:Lifted directed subtour elimination constraint for S={n+i,n+j,n+k}?D
4.3Capacity constraints
For any subset S?P∪D,let R(S)denote the minimum number of vehicles needed to visit all nodes in S.The constraint x(δ(S))≥2R(S)is then a valid inequality.Although computing R(S)is di?cult,a lower approximation is provided by?q(S)/Q?where q(S)= i∈S q i.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 positive and negative q i values.In particular, q(P∪D)=0by de?nition and,in the absence of time windows,all nodes can be visited by the same https://www.sodocs.net/doc/e216072203.html,eful inequalities can,however,be obtained by restricting S to be a subset of either P or D and setting q(S)=| i∈S q i|.It should be emphasized that in our context,the quantity?q(S)/Q?estimates the required number of vehicle passages in set S rather than the number of vehicles per se.Indeed,by leaving and entering set S more than once,the same vehicle may be able to visit all nodes in the set even though q(S)>Q.
4.4Precedence constraints
Consider a node set S such that0∈S,i∈S,n+i∈S and2n+1∈S for at least one user i.Because node i must be visited before node n+i,x(S)≤|S|?2and,equivalently, x(δ(S))≥4.The same applies to node sets S such that0∈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 pickup and delivery problem.
4.5Generalized order constraints
Let U1,...,U m?N be mutually disjoint subsets and let i1,...,i m∈P be users such that 0,2n+1∈U l and i l,n+i l+1∈U l for l=1,...,m(where i m+1=i1).The following inequality,introduced by Ruland and Rodin(1997),is also valid for the DARP:
m
l=1x(U l)≤m l=1|U l|?m?1.(35)
In the directed case,these inequalities can be lifted in two di?erent ways,as shown by the following propositions.
Proposition5.The following inequality is valid for the DARP:
m
l=1x(U l)+m?1 l=2x i1,i l+m l=3x i1,n+i l≤m l=1|U l|?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 since they all belong toδ+(i1).To demonstrate the validity of the lifting,we show that if any of the lifted arcs is part of the solution then there are at least two subsets U l for which x(U l)≤|U l|?2.Since x(U l)≤|U l|?1for all other subsets,the validity of the inequality follows directly.The details of the proof are given in the appendix.
Remark1.The numbering of the sets U1,...,U m is arbitrary and leads,by symmetry,to m!possibly di?erent liftings.
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 x i,n+j+x n+j,i+x j,n+k+
5.
x n+k,j+x k,n+i+x n+i,k+x ij+x i,n+k≤2.This is illustrated in Figure
Figure5:Lifted generalized order constraint with m=3(?rst lifting) Proposition6.The following inequality is valid for the DARP:
m
l=1x(U l)+m?2 l=2x n+i1,i l+m?1 l=2x n+i1,n+i l≤m l=1|U l|?m?1.(37)
Proof.The reasoning is similar to that of Proposition5and 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 x i,n+j+x n+j,i+x j,n+k+ x n+k,j+x k,n+i+x n+i,k+x n+i,n+j≤2.This is illustrated in Figure6.Observe that in inequality(37),the value of m must be larger than or equal to4for the second term of the inequality to have any
e?ect.
Figure6:Lifted generalized order constraint with m=3(second lifting)
4.6Order 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?ows 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.
Proposition7.Let i1,...,i m be m users and let H?P∪D and T h?P∪D,h=1,...,m be node sets such that{i h,n+i h}?T h and H∩T h={i h}.The following inequality is valid for the DARP:
x(H)+
m
h=1x(T h)≤|H|+m h=1|T h|?2m.(39)
Proof.First observe that x(H)≤|H|?1and x(T h)≤|T h|?1for h=1,...,m.If x(T h)=|T h|?1for any given subset T h,then there exists a path connecting all nodes in T h,
including i h and n+i h.However,this path cannot?nish at node i h because of the precedence constraint for user i h.Letαbe the number of sets T h for which x(T h)=|T h|?1.Then, x(δ+(H))≥α.Since x(δ+(H))=x(δ?(H))and2x(H)+x(δ+(H))+x(δ?(H))=2|H|,one obtains that2x(H)≤2|H|?2αand thus x(H)≤|H|?α.Finally,because x(T h)≤|T h|?2 for the remaining m?αsets,one may conclude that x(H)+ m h=1x(T h)≤|H|?α+ m h=1(|T h|?1)?(m?α),which simpli?es to expression(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)≤4and is illustrated in Figure
7.
Figure7:Generalized order matching constraint with m=2
Remark2.Inequality(39)is stronger than the corresponding TSP comb constraint,de?ned for m≥3and odd,and which is obtained from(39)by replacing the right-hand side with |H|+ m h=1|T h|?(3m+1)/2.
4.7Infeasible path inequalities
Ride time constraints may give rise to paths that are infeasible in an integer solution but nonetheless feasible in a fractional solution.Forbidding such paths can be accomplished through the valid inequality introduced in the next proposition.
Proposition8.For any directed path P={i,k1,k2,...,k p,n+i}such that t i,k
1+d k
1
+
t k
1,k2+d k
2
+···+t k
p,n+i
>L the following inequality is valid for the DARP:
x i,k
1
+
p?1
h=1x k h,k h+1+x k p,n+i≤p?1.(40)
Proof.Suppose that the value of the left-hand side of inequality(40)is equal to p in a feasible integer solution.Then there is a single arc from path P not belonging to the solution because the path contains p+1arcs.Since 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. But if the triangle inequality holds for travel times,the path from i to n+i has a larger duration than that of P.As a result,its duration must exceed L,which contradicts the assumption that the solution is feasible.
Remark3.As will be shown in the next section,the case p=1can be handled directly through the simple elimination of arcs in a preprocessing step.
5Branch-and-Cut Algorithm
Branch-and-cut is a popular approach for solving combinatorial problems.It has been applied successfully to several routing problems(see,e.g.,Ascheuer et al.,2001;Gendreau et al., 1998;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 Section5.1,the algorithm?rst solves the LP relaxation of the problem.If the solution to the LP relaxation is integer,an optimal solution has been identi?ed.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 Section5.2.
To control the branch-and-cut process,additional variables y k i= j∈N x k ij are introduced in the formulation.At each node of the search tree,if the y k i variables are all integer in the current relaxation but there is at least one fractional x k ij variable,the separation heuristics are executed in the hope of identifying violated valid inequalities.If at least one of the heuristics succeeds in?nding one or more violated inequalities,the relaxation is solved with all identi?ed cuts and the heuristics are executed again.The cut generation process at a node terminates when all heuristics fail to?nd any violated inequality.If the solution to the relaxation is still fractional after the generation of cuts,branching is performed on a fractional y k i variable,if there is any,or on a fractional x k ij variable,otherwise.The variable selected for branching is that whose value is the farthest from the nearest integer.
Prior to the application of the branch-and-cut algorithm,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.
Finally,an initial set of valid inequalities is added prior to solving the initial LP relaxation:all bounds on time and load variables,lifted2-node subtour elimination constraints,generalized order constraints with m=2and|U l|=2,and precedence constraints with|S|=3are introduced in a pool of cuts whose violations are checked at each node of the branch-and-bound tree,including those where not all y variables take integer values.
5.1Preprocessing
This section describes the time window tightening,arc elimination and variable?xing steps that can be applied prior to the branch-and-cut algorithm.
5.1.1Time 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 e i=max{0,e n+i?L?d i}and l i= min{l n+i?t i,n+i?d i,T}.In the case of an inbound user,the time window at the destination node can be tightened by setting e n+i=max{0,e i+d i+t i,n+i}and l n+i=min{l i+d i+L,T}. The time window on nodes0and2n+1can also be tightened by setting e0=e2n+1= min i∈P∪D{e i?t0i}and l0=l2n+1=max i∈P∪D{l i+d i+t i,2n+1}.
5.1.2Arc elimination
Formulation(1)-(14)is de?ned 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 e i+d i+t ij>l j;
?arcs(i,j)and(j,n+i)with i∈P,j∈N are both infeasible if t ij+d j+t j,n+i>L. As?rst 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 P={j,i,n+j,n+i}is infeasible;
?arc(n+i,j)is infeasible if path P={i,n+i,j,n+j}is infeasible;
?arc(i,j)is infeasible if paths P1={i,j,n+i,n+j}and P2={i,j,n+j,n+i}are both infeasible;
?arc(n+i,n+j)is infeasible if paths P1={i,j,n+i,n+j}and P2={j,i,n+i,n+j} are both infeasible.
In the presence of ride time constraints,further 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 windows and ride time constraints.Incompatible user pairs {i,j}can be identi?ed 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},{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 B i=e i 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 F i at node i in a path{i,i+1,...,q}can be computed as
follows:
F i=min
i≤j≤q i
where W i denotes the waiting time at node i,P i denotes the ride time of the user whose destination node is i if n+1≤i≤2n,and P i=?∞,otherwise.This de?nition of the forward time slack generalizes that of Savelsbergh(1992)for the TSP with time windows.
5.1.3Variable?xing
The identi?cation of incompatible user pairs can also be used to permanently assign users to speci?c vehicles.If the?eet of vehicles is homogeneous,one can 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 di?erent 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 vehicle k of users that are incompatible with i.Finding a clique of maximum cardinality in G′may be very time consuming when n is large.In our implementation,we thus use a greedy heuristic described by Johnson(1974).
5.2Separation Heuristics
This section describes the separation heuristics used to identify violated inequalities.When all y k i variables are integer but there is at least one fractional x k ij variable,the following heuristics are executed sequentially.
5.2.1Subtour elimination constraints
The identi?cation of violated inequalities of the form x(S)≤|S|?1can be achieved by solving a series of maximum?ow 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 heuristic inspired from that proposed by Augerat et al.(1999).
Using the fact that2x(S)+x(δ(S))=2|S|in a feasible integer solution,violations of(29) can be identi?ed by?nding node sets S such that
x(δ(S))?2 i∈P+(S) j∈S x n+i,j<2.(42)
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(42).Whenever a node is removed from set S,its reinsertion is declared tabu forθiterations.The heuristic runs for a preset number of iterations(100in our implementation)and may identify several violated inequalities during a single execution.At each iteration,the current set S is also checked for possible violations of inequality(33).To this purpose,the node with the largest outgoing ?ow is numbered as i1and 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?ow is numbered as i1.
5.2.2Capacity 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 subset S?P or with a random subset S?D.At each iteration,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.Again,the heuristic runs for100iterations and multiple violated inequalities may be identi?ed during a single execution.
5.2.3Precedence constraints
As explained by Ruland(1995),identifying violated precedence constraints can be done by solving a multi-terminal maximum?ow problem for each user.To check whether the precedence constraint for user i is violated,one can compute the maximum?ow from the sources0and n+i to the sinks i and2n+1.If the value of this?ow is less than2,then a precedence constraint x(S)≤|S|?2is violated for a set S such that0,n+i∈S and i,2n+1∈S.The set S corresponds to one of the shores of the corresponding minimum cut.We have implemented this approach and use the maximum?ow algorithm provided in the GTL library(see http://infosun.fmi.uni-passau.de/GTL).
5.2.4Generalized order constraints
Ruland(1995)proposed an approach to identify violated generalized order constraints with m=2.This approach requires the computation of O(|N|2)maximum?ows.Here,we use instead three simple heuristics,the second and third of which take advantage of the lifted forms(36)and(37).
The?rst heuristic attempts to identify violations of inequalities(35)for the special case where m=2and|U1|=|U2|=3.For each pair of users i,j∈P,we?rst form the subsets U1={i,n+j}and U2={j,n+i}.We then identify nodes k1and k2such 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?nding violations of inequalities(36)and(37) for the special case where m=3and|U1|=|U2|=|U3|=2.The second heuristic?nds,for
each user i,a user j that maximizes x i,n+j+x n+j,i+x ij.It then?nds a user k such that the left-hand side of(36)is maximized.Similarly,the third heuristic?nds,for each user i,a user j that maximizes x i,n+j+x n+j,i+x n+i,n+j and then?nds a user k that maximizes the left-hand side of(37).
5.2.5Order matching constraints
To identify violated order matching constraints,Ruland(1995)has proposed an algorithm requiring the solution of O(|N|2)maximum?ow problems.
We again resort to a simpler heuristic that aims to identify violations of(39)for the special case where m=2,|H|=3and|T1|=|T2|=3.For each pair of users i,j∈P,we?rst form the subsets T1={i,n+i}and T2={j,n+j}.We then select the nodes k1and k2such that x(T1)and x(T2)are maximized.Finally,we select the node k such x(H)is maximized and check for a violation of the resulting inequality.
5.2.6Infeasible path inequalities
Violated path inequalities are identi?ed by means of a path-construction heuristic applied to each user i∈P.The heuristic starts with node i and gradually constructs a path P i={i,k1,k2,...}by iteratively moving from the current node k j to the node k l that
maximizes the value of x k
j k l .The process stops if a cycle is formed or if the heuristic reaches
either node n+i or node2n+1.If the path stops at node n+i,it is checked for a violation of inequality(40).
6Computational Experiments
The branch-and-cut algorithm just described was implemented in C++by using ILOG Concert1.3and CPLEX8.1.It was run on a2.5GHz Pentium4computer with512Mb of memory.
The algorithm was applied to two sets of randomly generated instances comprising up to32 users.In an instance with n users,where n is even,users1,...,n/2are assumed to formulate outbound requests while users n/2+1,...,n formulate inbound requests.In all instances, the coordinates of pick-up and drop-o?nodes are randomly and independently chosen in the square[?10,10]×[?10,10]according to a uniform distribution,and the depot is located at the center of this square.For every arc(v i,v j)∈A,the routing cost c ij and travel time t ij are equal to the Euclidean distance between the two nodes.
A time window[e i,l i]is also associated to each node.For an outbound user i,a time window was generated by?rst choosing a number e n+i in the interval[420,1080](i.e.,between7h00 and18h00)and then setting l n+i=e n+i+15.For an inbound user,the value of e i was chosen in the interval[360,1020]and the value of l i was set equal to e i+15.Time windows on the
origin nodes of outbound requests and on the destination nodes of inbound requests are set as explained in Section5.1.1.
In the?rst set of instances,Q=3for every vehicle,q i=1for every user and the maximum ride time L is equal to30minutes.In the second set,Q=6for every vehicle and the values of q i are chosen randomly according to a uniform distribution on the set{1,...,Q}.The maximum ride time is set equal to60minutes.In all cases,we assume that service time is proportional to the number of passengers and d i=q i.The maximum route duration is set to720minutes.The?rst set of instances represents the context where cars are used for the transportation of individuals whereas the second set re?ects the situation of a transporter using mini-busses for the transportation of individuals or groups of individuals.
Tables1and2provide the characteristics of these instances,the number of constraints and variables in the resulting integer programming formulation,the value of the LP relaxation and the cost of a heuristic solution.This solution was obtained by running the tabu search algorithm of Cordeau and Laporte(2003b)for1000iterations.
In Tables3and4we indicate the best lower bounds reached by running the branch-and-cut algorithm on each instance with a maximum CPU time limit of60minutes.Column CPLEX indicates the bound obtained by running CPLEX itself with no user-cut generation but after the application of the preprocessing steps of Section5.1.The next columns report the bounds obtained by generating the lifted bounds on time and load variables(Section4.1) with only one of the following:subtour eliminations constraints(SEC),capacity constraints (CC),generalized order constraints(GOC),order matching constraints(OMC)and infeasible path constraints(IPC).Average values reported on the last line of each table show that the strength of the various inequalities varies from one set of instances to the other.Furthermore, for the?rst set of instances,using any kind of inequality alone did not yield signi?cantly Table1:Characteristics of the?rst set of instances
p2-16216330855733270.25298.00
p2-2022033011341050300.55345.14
p2-2422433018951679304.10375.15
p2-2822833026102344337.60417.12
p2-3223233035363098350.15466.53
p3-163163309201020197.39258.10
p3-2032033015201776238.88282.81
p3-2432433019572311312.90374.18
p3-2832833026093059368.43478.58
p3-3233233036764593291.22423.08
p4-164163309601299229.12294.82
p4-2042033015022104256.25343.37
p4-2442433021713188254.20334.41
p4-2842833029114479267.28395.90
p4-3243233035755564376.65491.82
Table2:Characteristics of the second set of instances g2-16216660666600223.14269.98 g2-2022066012321176280.41374.77 g2-2422466013799379307.64429.90 g2-2822866024482163336.35427.06 g2-3223266027472395400.81544.95 g3-16316660782917207.28256.70 g3-2032066012001459251.37301.41 g3-2432466017962045319.89449.88 g3-2832866024012807325.37451.09 g3-3233266026253293391.38562.73 g4-164166608751220217.09273.60 g4-2042066011461642257.69363.75 g4-2442466017532639300.70423.10 g4-2842866024843399261.27363.70 g4-3243266028804221351.21516.00 Instance CPLEX SEC CC GOC OMC IPC
Avg.344.98347.50343.56338.91344.30344.67
Table4:Best lower bounds-second set of instances
g2-16268.36268.36268.36268.36268.36268.36
g2-20370.77370.77370.77370.77370.77370.77
g2-24397.60429.27426.05429.27415.52413.46
g2-28409.30418.95424.27416.53418.02417.29
g2-32502.84528.20523.01516.85508.62510.84
g3-16255.18255.18255.18255.18255.18255.18
g3-20300.45300.45300.45300.45300.45300.45
g3-24382.67402.90411.80399.11395.84393.96
g3-28375.30383.36382.60381.95377.08382.27
g3-32443.54453.94462.18456.98452.35451.84
g4-16270.24270.24270.24270.24270.24270.24
g4-20346.42361.56361.56361.35361.56361.56
g4-24349.38360.96359.82357.07358.03356.16
g4-28285.90293.76295.17295.42291.11294.37
g4-32390.11405.24414.47411.68406.33405.04
Combining the di?erent types of inequalities does,however,yield important improvements over the basic version of CPLEX.Tables5and6reports the results obtained by using the basic version of CPLEX and by using the?ve types of inequalities concurrently.In both cases,a maximum of12hours of CPU time was allowed for the solution of each instance.We Table5:Comparisons with CPLEX-?rst set of instances
p2-16*298.000.0154*298.000.033627
p2-20*345.140.11628*345.140.1216077
p2-24*375.15 5.9621,844*375.15 3.125,368203
p2-28*417.05188.57472,114*417.058.5912,516359
p2-32448.70720.00909,500*466.53142.51148,501545
p3-16*258.100.05288*258.100.0721836
p3-20*282.81 1.375,118*282.81 1.333,42085
p3-24*374.07 2.046,042*374.07 1.121,785108
p3-28*462.80379.36802,594*462.80153.56240,402269
p3-32354.45720.00446,400354.71720.00377,030802
p4-16*294.82 3.4522,233*294.82 1.355,76474
p4-20*341.93524.591,403,353*341.93308.97765,146181
p4-24309.73720.00711,000*332.96685.54451,9151203
p4-28337.34720.00584,900341.00720.00433,158390
p4-32412.11720.00263,100418.71720.00267,975471