Introducing heterogeneous users and vehicles into models and algorithms for the dial-a-ride problem.

Dial-a-ride problems deal with the transportation of people between pickup and delivery locations. Given the fact that people are subject to transportation, constraints related to quality of service are usually present, such as time windows and maximum user ride time limits. In many real world applications, different types of users exist. In the field of patient and disabled people transportation, up to four different transportation modes can be distinguished. In this article we consider staff seats, patient seats, stretchers and wheelchair places. Furthermore, most companies involved in the transportation of the disabled or ill dispose of different types of vehicles. We introduce both aspects into state-of-the-art formulations and branch-and-cut algorithms for the standard dial-a-ride problem. Also a recent metaheuristic method is adapted to this new problem. In addition, a further service quality related issue is analyzed: vehicle waiting time with passengers aboard. Instances with up to 40 requests are solved to optimality. High quality solutions are obtained with the heuristic method.

Introduction

Dial-a-ride problems deal with the transportation of people between pre-specified pickup and delivery locations. Usually, users specify time windows for either the pickup or the drop-off location, depending on the type of request. In the case of an outbound request (e.g., from home to a hospital), a time window is defined for the delivery location (destination). In the case of an inbound request (e.g., from the hospital back home), a time window for the pickup location (origin) is given. In addition, user ride time as well as maximum route duration limits have to be respected. The objective is to generate a transportation plan serving all patient transportation requests at minimum routing costs. In this article we introduce heterogeneous users and vehicles into the standard dial-a-ride problem (DARP), as defined, e.g., by Cordeau (2006). The resulting problem will be denoted as heterogeneous DARP (HDARP).

The introduction of heterogeneous users and vehicles is motivated by observations made at the Austrian Red Cross (ARC) in the field of patient transportation. The ARC distinguishes three patient types. A patient may demand to be transported seated, on a stretcher, or in a wheelchair. In addition, an accompanying person may be present. The ARC disposes of two different types of vehicles. Each type provides different capacities for four modes of transportation (staff seat, patient seat, stretcher, wheelchair place). In the following, these are referred to as resources. Resource 0 refers to staff seats, resource 1 to patient seats, resource 2 to stretchers, and resource 3 to wheelchair places. Each passenger can only be transported by a vehicle that provides the appropriate resource. Usually, accompanying persons use staff seats; seated patients take patient seats; a patient demanding a stretcher is transported on a stretcher; and a patient in a wheelchair occupies a wheelchair place. However, certain so-called upgrading conditions apply: an accompanying person may use the patient seat, in case there is no vacant staff seat available, or, in case neither of the two is currently available, he/she can sit on the stretcher. According to Austrian law, a seated patient cannot use a staff seat, he/she may, however, sit on the stretcher, in case there is no additional patient seat available. Patients that demand a stretcher can only be transported on a stretcher; patients in wheelchairs can only be transported by vehicles providing space for a wheelchair.

Another aspect that has only rarely been considered in the literature is user waiting time while aboard a vehicle. In general, one can assume that some waiting time within the defined time window is acceptable for users of dial-a-ride systems. However, when already aboard a vehicle, waiting times should be avoided in order to keep user inconvenience at a reasonably low level. Psychologically, waiting for a vehicle is perceived in a different way than waiting time on board a vehicle (except for loading or unloading operations). The latter decreases the perceived service quality more than the former. Therefore, in the following, we will also integrate a term into the objective function that penalizes vehicle waiting time with passengers aboard.

Previous publications, considering heterogeneous versions of the DARP, mostly employ metaheuristic methods in order to solve the problem at hand. Toth and Vigo (1997), e.g., devise a parallel insertion heuristic and a tabu thresholding algorithm for a DARP considering two modes of transportation (seated and in a wheelchair) and several types of vehicles. Another heterogeneous version of the DARP is described in Melachrinoudis et al. (2007). Different types of vehicles in terms of capacity limits but only one mode of transportation are considered. The solution method developed is a tabu search algorithm. Heterogeneous vehicles in terms of different capacity limits for one mode of transportation in the context of the DARP are also considered in the work of Rekiek et al. (2006). The proposed problem is solved by a grouping genetic algorithm.

In a dynamic environment, Beaudry et al. (2009) adapt the tabu search heuristic developed by Cordeau and Laporte (2003) to solve a heterogeneous DARP that arises in large hospitals. It involves transportation requests demanding three different modes of transportation (seated, on a bed, or in a wheelchair) and several different types of vehicles. Hanne et al. (2009) develop a computer based planning system for a dynamic problem in a large German hospital, considering hospital specific constraints such as multi-dimensional capacities.

In terms of exact algorithms for DARPs, early publications such as Psaraftis (1980, 1983) and Desrosiers et al. (1986) focus on dynamic programming algorithms. More recently, Cordeau, 2006 proposed a branch-and-cut algorithm for the DARP. The algorithm is based on a 3-index mixed-integer programming formulation. New valid inequalities as well as previously developed ones for the pickup and delivery problem and the vehicle routing problem are employed. The largest instance solved to optimality comprises 36 requests. Two new branch-and-cut algorithms followed in Ropke et al. (2007). Instead of a 3-index formulation, more efficient 2-index problem formulations and additional valid inequalities are proposed. In this case, a 96-request instance is the largest instance solved to optimality. The algorithms of Cordeau (2006) and Ropke et al. (2007) are thus state-of-the-art in terms of exact algorithms for the DARP. The work presented in this article is based on their findings.

In a related paper, Parragh et al. (2009a) incorporate driver related constraints as well as heterogeneous users and vehicles into a 3-index and a set partitioning formulation of the DARP. These involve, e.g., the assignment of drivers and additional personnel to vehicles and the scheduling of lunch breaks. In these formulations, maximum ride times are only implicitly taken into account by constructing time windows at both the pickup and the drop-off location. This modification makes the application of column generation possible. The authors introduce a collaborative framework combining variable neighborhood search and column generation. They report results on instances with up to 54 requests.

For further information on the DARP, we refer the reader to Cordeau and Laporte (2007). For overviews of the more general pickup and delivery problem see Berbeglia et al. (2007), Parragh et al. (2008a,b).

The contributions of this article are threefold. First, two problem formulations, a 3-index and a 2-index formulation, for the HDARP are introduced (Section 2). A major new aspect in the 2-index formulation refers to how the differing capacity limits, depending on which vehicle travels along which arc, are integrated. Second, all valid inequalities proposed for the standard DARP are adapted to the new problem (Section 3) and employed in a 3-index and in a 2-index-based branch-and-cut algorithm (Section 4). Third, a state-of-the-art variable neighborhood search (VNS) heuristic (Mladenovic and Hansen, 1997), initially designed for the standard DARP (Parragh, 2009, Chapter 3), is adapted and applied to the HDARP (Section 4.3). Computational results for three sets of instances are reported and the impact of penalizing user waiting time while aboard a vehicle is analyzed. The article ends with a conclusion.

Problem definition

The HDARP is modeled on a complete directed graph G = (V, A) where V is the set of all vertices and A the set of all arcs. A set K of m heterogeneous vehicles has to serve all n transportation requests. Each vehicle k ∈ K is associated with a vector C that gives the amount of resource r available on vehicle k. The ARC disposes of two basic vehicle types. Type 1 (T1) provides 1 staff seat, 6 patient seats, and 1 wheelchair place. Type 2 (T2) provides two staff seats, one patient seat, one stretcher, and one wheelchair place (see also Fig. 1). Each vehicle k has to start its route at the origin depot 0 and end at the destination depot 2n + 1, respecting a route duration limit T. For each arc (i, j) and each vehicle k ∈ K a non-negative travel cost and a non-negative travel time is considered. Each transportation request consists of a pickup and delivery vertex pair {i, n + i}. The set of pickup vertices is given by P = {1, … , n}, the set of delivery vertices by D = {n + 1, … , 2n}. At every pickup vertex one patient waits to be transported. This patient may demand one of three different modes of transportation. Passengers may have to be transported seated (), on a stretcher (), or in a wheelchair (); each patient may be accompanied by a friend, relative or nurse (). The demand at every delivery vertex is equal to for all r ∈ R = {0, 1, 2, 3}. Every user either specifies a time window [e, l] for the pickup (origin) or the drop-off (destination) location and beginning of service has to start within this time window. In case a vehicle arrives too early, it has to wait until service is possible. A maximum passenger ride time is also considered, in order to keep quality of service at a reasonably high level. At each vertex loading or unloading operations last for a given service time d. Thus, the set of all vertices is given by V = P ∪ D ∪ {0, 2n + 1}, and the set of all arcs by A = {(i, j)∣i ∈ V⧹{2n + 1}, j ∈ V⧹{0}, i ≠ j}.

Fig. 1

Vehicle types at the ARC.

As mentioned above, so-called upgrading conditions apply. Patients demanding to be transported seated may use a patient seat or the stretcher. Patients demanding a stretcher can only be transported on a stretcher. The same applies to wheelchair passengers. Accompanying persons, however, may use a staff seat, a patient seat, or the stretcher, if no other transportation mode is available. Table 1 gives on overview of the different upgrading options.

Table 1

Transportation mode upgrading options.

Passenger type	Staff seat	Patient seat	Stretcher	Wheelchair place
Accompanying person	x	x	x
Seated patient		x	x
Patient on stretcher			x
Patient in wheelchair				x

The model can be formulated with the following decision variables:

A 3-index formulation

The mathematical program, covering all of the above described real world conditions, based on the DARP formulation by Cordeau (2006), is as follows:s.t.:

The objective function (1) minimizes total routing costs and penalizes waiting time when passengers are aboard the vehicle. The penalty term ρ will be defined in such a way that waiting time while on board a vehicle is avoided. Constraints (2) and (3) guarantee that each request is served exactly once and that each origin-destination pair is visited by the same vehicle. Equalities (4)–(6) ensure that each vehicle starts at and returns to the depot at the end of its route.

Consistency with respect to resource and load variables is guaranteed by constraints (7)–(9). Upgrading constraints for resources 0, 1, and 2 are given in (8). They guarantee that capacity restrictions regarding these resources are not violated and that each patient demanding resource 0, 1 or 2 can only be loaded if there is either enough capacity of the resource demanded or another one with a higher number (0 = staff seat, 1 = patient seat, 2 = stretcher). Constraints (9) guarantee that patients demanding resource 3 (wheelchair place) can only be transported if there is enough capacity of resource 3.

Equalities (12) define the arrival times for each vertex. Constraints (13) guarantee that beginning of service only starts after having arrived at vertex i. Waiting times are set in (14). Waiting is only penalized if there is someone aboard the vehicle when arriving at vertex i. This is ensured by the way are set in (11). Note that constraints (12) and (13) also take care of subtour elimination given that (t + d) > 0 for all i, j ∈ V, i ≠ j.

Equalities (15) define the ride time of each user. Total route duration is limited by (16), time window and maximum ride time compliance is ensured by (17) and (18).

Note that the capacity restrictions given by (8) and (9) are tailored to the problem situation of the ARC, they can, however, be rewritten in a more generic way as follows:where the set R simply contains the different resources that can be used for the transportation of passengers demanding resource r. In our case, these would be R0 = {0, 1, 2}, R1 = {1, 2}, R2 = {2}, and R3 = {3}. Using this generic formulation, all different types of up-or down-grading could easily be accommodated with any number of resources.

Non-linear constraints

The above problem formulation contains a number of non-linear constraints, (7), (12) and (14). These can be reformulated in a linear way using so-called “big M” terms. All transformations described below are based on Cordeau (2006). The first set of non-linear constraints are the load propagation inequalities given in (7). Using “big M” terms tailored to each constraint, they can be rewritten as follows:where and . The “big M” terms employed and denoted by represent the largest value the term may take. Furthermore, as shown by Desrochers and Laporte (1991) they can then be lifted intofor the DARP by taking the reverse arc into account. Time consistency constraints given in (12) are transformed as follows:

Here, the original equality constraints are replaced by two sets of constraints, providing an upper and a lower bound on A. In the case where arc (i, j) is not used by vehicle k, these bounds should not be binding. Therefore, if arc (i, k) is not used by vehicle k, the upper bound is increased by l2 (the later time window at the end depot) and the lower bound is decreased by . Following a similar idea, inequalities (14), responsible for defining vehicle waiting time at each vertex i, are replaced by the following linear constraints,

Variable aggregation

The complexity of a problem formulation can be substantially decreased if all those variables that are homogeneous across the different vehicles are aggregated. In our case this is possible for all time related variables as these do not differ from one vehicle to the other. This implies that each vehicle drives at the same speed which is a reasonable assumption especially in rural areas. Therefore, can be replaced by t and aggregate time variables A, B, L, and can be used at all vertices except the two depots 0 and 2n + 1. Furthermore, the binary variables indicating whether a vehicle arrives full or empty at a vertex i, can be substituted by y. Here, the reasoning is that every vertex (except the depots) can only be visited exactly once. Therefore, some vehicle will at most once arrive at each vertex ∈P ∪ D either empty or with passengers aboard. All constraints on these variables have to be reformulated as follows. Equalities (12), defining the arrival time at each vertex, are now split into three separate sets of equalities:

Inequalities (13), ensuring that the beginning of service variables are correctly set, are also transformed into three separate sets of constraints:

Constraints (14) and (11), responsible for defining a correct lower bound on waiting time with passengers aboard, are replaced byand time window and ride time constraints given in (15), (17) and (18) are replaced by

Let us further assume that the per vehicle arc costs are homogeneous across the different vehicle types, i.e. can be replaced by c. Eventually, the objective function can be rewritten as follows:

A 2-index formulation

When comparing the 3-index-based branch-and-cut algorithm for the DARP of Cordeau (2006) to the 2-index-based branch-and-cut algorithms proposed in Ropke et al. (2007), the 2-index versions obviously outperform the earlier 3-index-based method. Given the additional complexity of our problem, it can be assumed that a branch-and-cut algorithm based on a 2-index formulation will also perform better than a branch-and-cut algorithm based on a 3-index program. We note that, in order to reformulate the HDARP as a 2-index program, we have to adhere to the basic assumption that per vehicle arc costs c and per vehicle arc travel times t are homogeneous across the different vehicle types considered. In the transportation of the ill or disabled this assumption is reasonable; usually vehicles do not differ much in their dimensions or their weight but in terms of their on board equipment, expressed in the different transportation mode options. However, the reformulation of the HDARP as a 2-index program is still not as straightforward as one may assume, given the heterogeneous vehicle fleet and the resulting differing capacity limits, depending on which vehicle travels along which arc. How this issue can be resolved is described in the following.

In order to account for the heterogeneous nature of the vehicle fleet, artificial origin and destination depots for each vehicle are introduced. A similar approach is followed, e.g., in Baldacci et al. (2009) in the context of the heterogeneous fleet size and mix vehicle routing problem. Let D = {2n + 1, … , 2n + m} denote the set of vehicle origin depots and D = {2n + m + 1, … , 2n + 2m} the set of vehicle destination depots, the vertex set V is redefined as V = P ∪ D ∪ D ∪ D and the arc set as A = {(i, j):i ∈ V⧹D, j ∈ V⧹D, i ≠ j}. The demand/supply at the artificial depots is set to zero, for all k ∈ K, r ∈ R; and the total maximum vehicle capacity shall be defined as .

To impose request wise precedence and pairing constraints, let denote the set of all vertex subsets S ⊆ V, such that for all k ∈ K the respective origin depot 2n + k ∈ S and the destination depot 2n + m + k ∉ S, and there is at least one request i for which i ∉ S and n + i ∈ S. For the correct pairing of the artificial depots, define as the set of all vertex subsets U ⊆ V, such that for exactly one k ∈ K the origin depot 2n + k ∈ U and the destination depot 2n + m + k ∉ U, while for all other l ∈ K⧹{k} origin depots 2n + l ∉ U and destination depots 2n + m + l ∈ U.

Furthermore, to incorporate the corresponding vehicle capacity constraints, let denote the set of infeasible paths with respect to load violations H = {j1, … , j}, such that j1 ∈ D and let A(H) denote the arc set of H. Omitting superscript k from all decision variables, defining Q as the aggregated load when leaving vertex i, and as the aggregated demand/supply at vertex i, the HDARP can be formulated as the following 2-index program (inspired by Ropke et al., 2007):s.t.:

The objective function (41) minimizes total routing costs and penalizes waiting time with passengers aboard a vehicle. Constraints (42) and (43) guarantee that every vertex i ∈ P ∪ D is entered and left and that exactly one arc leaves (enters) every origin (destination) depot, i.e. i ∈ D (i ∈ D). Precedence and pairing restrictions are guaranteed by (44) and (45). As shown by Ropke et al. (2007) inequalities (44), initially designed for the single vehicle case, also apply in a multi-vehicle context. Since one depot per vehicle is considered in our mathematical model, the original right hand side ∣S∣ − 2 has to be replaced by ∣S∣ − (m + 1). Suppose the set is S = {2n + 1, … , 2n + m, n + i}. From each start depot one arc has to leave S (i.e. m arcs). One of these arcs will be part of a path leading to i ∉ S. If pairing and precedence constraints are respected, this path has to contain one arc leading back to n + i ∈ S. Finally, it has to arrive at one of the end depots, demanding another arc to be traversed leaving S and connecting to some end depot 2n + m + k ∉ S. Thus, m + 1 arcs have to leave S, reducing the number of arcs that can be used within S to ∣S∣ − (m + 1). Constraints (45) are depicted in Fig. 2. Here, d1+,d2+, and d3+ denote the origin depots of vehicles 1, 2, and 3, respectively, while d1−, d2−, and d3− denote the corresponding destination depots. At least one arc (represented as an arrow) has to leave set U (here U consists of d1+,d2−, and d3− and some other vertices) and enter set V⧹U, in order to ensure the correct pairing of the depots (here d1+ and d1−). Aggregate loading variables are introduced and set by constraints (46). These are needed for the definition of y, defined in (48). The actual loading restrictions are taken care of by inequalities (49); on every path , each starting at one of the depots, a loading restriction is violated. Constraints of this type are known as tournament inequalities (Ascheuer et al., 2000). Arrival times are set in (50), beginning of service and waiting times in (51) and (52). Equalities (53) determine each user’s ride time. Compliance with maximum route duration restrictions, time windows, and user ride time limits is guaranteed by (54)–(56).

Fig. 2

Pairing of artificial origin and destination depots.

All non-linear inequalities are reformulated as shown for the 3-index model. Thus, load propagation constraints (46) are substituted bywith and lifted intoby taking the reverse arc into account. Finally, arrival time inequalities (50) are replaced by,and vehicle waiting time inequalities (52) by

Valid inequalities

Cordeau (2006) and Ropke et al. (2007) introduced new and adapted a number of existing valid inequalities for the DARP. In the following section we will review and adapt all these families of inequalities to the HDARP. Most valid inequalities can be used to strengthen both formulations, some are only valid for either of the two. If not stated otherwise, in the 3-index formulation. Furthermore, for ease of exposition, let denote the number of arcs traversed in a set of vertices S.

Strengthened bounds on time and load variables

In both the 3-index and the 2-index formulation, bounds on time variables, also denoted as time windows, can be strengthened as follows (Cordeau, 2006; Desrochers and Laporte, 1991):

In the case of the 3-index formulation bounds on load variables can also be strengthened (based on Cordeau, 2006). In order to do so, let orig(i) be the set that contains the origin of i, if i is a destination, and otherwise the empty set. Thus, lower bounds can be strengthened as follows:

The intuition is that if arc (j, i) is used and , given that j is not the origin of i, the load at i hast to be greater or equal than (plus , if i is an origin).

In case of upper bounds, we have to distinguish between origins and destinations. If i is an origin, i.e. i ∈ P = {1, … , n},where and , applies. The intuition behind this strengthening is that if i is visited directly after the depot 0, the load can be at most . In addition, if the vertex visited directly after i is an origin and , then the load at i can be at most . If i is a destination, i.e. i ∈ D = {n + 1, … , 2n},is used. The last term used in (68) cannot be used to strengthen (69). Consider , if and then would have to be ⩽−1 which is clearly not valid, although visiting j after i is feasible with respect to capacity limits.

Subtour elimination constraints

Standard subtour elimination constraints are given by x(S) ⩽ ∣S∣ − 1 for S ⊂ P ∪ D. They ensure that a given subset of vertices S cannot contain a cycle through all the vertices forming S. As shown by Cordeau (2006) these constraints can be lifted in several ways, exploiting the fact that an origin i must be visited before its destination n + i. Let S = {i1,i2, … , i} ⊆ P ∪ D and . The predecessors of S are denoted as π(S) = {i ∈ P∣n + i ∈ S}, its successors as σ(S) = {n + i ∈ D∣i ∈ S}. The following four sets of inequalities are valid for the DARP (Cordeau, 2006) and are also applicable in the context of the HDARP:

All four sets of inequalities can be used to strengthen both the 3-index and the 2-index formulation. The first two sets of inequalities were originally introduced for the precedence constrained asymmetric TSP by Balas et al. (1995). They are also referred to as successor and predecessor inequalities, respectively. The ordering of the vertices in S does not play a role here. In the second two sets of inequalities the ordering of the vertices is important. They are based on those of Grötschel and Padberg (1985) for the asymmetric TSP.

These four sets of inequalities are further illustrated in Fig. 3. Dotted and dashed arcs represent liftings. In Fig. 3a an example for successor inequality (70) is given. The set S = {i, j} consists of two origin vertices. The successors of these two are their destinations σ(S) = {n + i, n + j}. Clearly, at most one (∣S∣ − 1 = 1) of the four arcs shown in the figure can be used in a feasible solution. The same is true for the example of (71) in Fig. 3b. For inequalities (72), (73) the reasoning is a bit different. Instead of sets, sequences of vertices are considered. Such that a sequence does not result in a cycle, at most h − 1 arcs can be used, given that the sequence consists of h arcs. If a reverse arc is used between two vertices, the vertices connected by this arc should no longer be connected to the rest of the sequence. This reduces the number of arcs that can be used in the sequence by two. For this reason a reverse arc (dashed in the figure) is counted twice. Furthermore, in case of inequality (72), illustrated in Fig. 3c, if a destination (successor) ∉S of some origin ∈S is connected to the first vertex ∈S, the first vertex cannot be connected to the origin of this destination anymore, reducing the number of vertices that can be used in the sequence by one. The same reasoning, considering predecessors instead of successors of the vertices ∈S, leads to inequalities (73). An according example is given in Fig. 3d.

Fig. 3

Lifted subtour elimination constraints (adapted from Cordeau, 2006).

Generalized order constraints

Generalized order constraints form another set of inequalities that can be used to strengthen the above formulations. They were originally proposed by Ruland and Rodin (1997). In order to properly define them, let U1, … , U be mutually disjoint subsets and let i1, … , i ∈ P be users such that either 0, 2n + 1 ∉ U (3-index formulation) or 2n + k, 2n + m + k ∉ U for all k ∈ K (2-index formulation) and i, n + i ∈ U for l = 1, … , h (i + 1 = i1). Then, the following inequalities are valid for the (H)DARP,

For the (H)DARP these can be lifted in two ways, as shown by Cordeau (2006),

Infeasible path constraints

Another popular way to further strengthen problem formulations in the field of vehicle routing, are constraints based on infeasible paths. We consider two families of valid inequalities that are based on the notion of infeasible paths. The first are referred to as strengthened infeasible path constraints, the second have been denoted as fork constraints.

In order to define strengthened infeasible path constraints, let denote the set of infeasible paths and for each set , let A(F) be the arc set of F and V(F) the vertex set. Then the following inequalities are valid for the (H)DARP (Ropke et al., 2007),

If F = (j1, … , j) denotes an infeasible path, they can be strengthened into so-called tournament constraints (see Ascheuer et al., 2000; Ropke et al., 2007),

If F links a vertex pair {i, n + i}, and F = (i, j1, … , j,n + i) is infeasible due to time window or ride time constraints, (77) can be lifted in the following way (Cordeau, 2006),

If both path F = (j1,…,j) and the reverse path F′ = (j, … , j1) are infeasible, the following inequality can be applied (Ropke et al., 2007),

These inequalities are generated for each vehicle in turn in case of the 3-index formulation. In case of the 2-index formulation has to be replaced by x and only time related infeasibilities are considered since these restrictions are homogeneous across all vehicles.

The second family of inequalities that is based on the notion of infeasible paths has been proposed by Ropke et al. (2007) in the context of the PDPTW and the DARP. In this case we only consider paths that are infeasible regarding time infeasibilities. They are eliminated by taking a whole bundle of such infeasible paths into account, sharing a number of common arcs. Let F = (j1, … , j) be a feasible path. By adding a vertex at the beginning and at the end of this path an infeasible path is generated. A bundle of infeasible paths is defined by (i, F, l) being infeasible for every i ∈ S and l ∈ T, S, T ⊂ V, resulting in the so-called fork inequality,

These inequalities can be lifted in two ways. In the first, an additional set of vertices T is appended to each vertex j ∈ {j1, … , j} resulting in sequences (k, j1, … , j, l) that are infeasible for every k ∈ S,l ∈ T, 1 ⩽ i ⩽ h. This so-called outfork inequality is given by

In the second, a set of vertices S is added before each vertex j ∈ {j2, … , j}, leading to sequences (k, j, … , j, l) that are infeasible for every k ∈ S, l ∈ T, 1 ⩽ i ⩽ h. This so-called infork inequality is given by

Capacity inequalities

For the 3-index formulation additional valid inequalities that are based on capacity issues can be appended to the model. We distinguish two families of inequalities, the adapted rounded capacity inequalities and the strengthened capacity inequalities. In order to define the former, let denote the minimum number of vehicles necessary to visit all nodes in S ⊆ P ∪ D, then the constraint is a valid inequality for the DARP (Cordeau, 2006), . A lower approximation of is given by ⌈q(S)/C⌉, where .

This so-called rounded capacity inequality can be adapted to the HDARP as follows. Let denote the cumulative capacity limit for resource r on vehicle k and be , then the constraintis a valid rounded capacity inequality for the HDARP.

Introduced in Ropke et al. (2007), the strengthened capacity cuts for the DARP are defined as follows. Let S, T ⊆ P ∪ D denote two disjoint subsets such that q(S) > 0. Furthermore, let , the strengthened capacity inequality for the DARP is given bywith . In case of the HDARP q(·) has to be replaced by q(·) as defined above, and instead of has to be used.

Reachability constraints

In case of the 2-index formulation another group of inequalities, introduced in Lysgaard (2006) and used in the context of the DARP by Ropke et al. (2007), can be used to strengthen the model. It is based on the notion of conflicting vertices. Vertices are said to be conflicting if they cannot be served by the same vehicle. For every vertex i ∈ P ∪ D let define the minimum arc set such that any feasible sequence from an origin depot 2n + k (k ∈ K) to vertex i can be constructed only using arcs . Furthermore, let denote the minimum arc set such that any feasible sequence from i to a destination depot 2n + m + k (k ∈ K) only traverses arcs . Finally, let be a set of conflicting vertices. defines the reaching arc set of and the reachable arc set of . The following two inequalities can be defined,with and .

Branch-and-cut algorithms

Based on each of the above introduced formulations, we have implemented a branch-and-cut algorithm. The algorithm based on the 3-index program will be denoted as 3indexBC in the following, the algorithm based on the 2-index formulation as 2indexBC. Branch-and-cut algorithms combine the branch-and-bound and the cutting-plane idea. Every mixed integer program (MIP) can be reformulated as a linear program (LP) by dropping all integrality constraints. In our case, these are (19) and (20) in the 3-index, and (57) and (58) in the 2-index formulation. The optimal solution to the LP relaxation will yield a lower bound for the original MIP.

Branch-and-cut algorithms depart from the solution of the LP relaxation, considering only a reasonable subset of the original constraints. Typically, all constraint families of exponential size are not included. In case of 2indexBC, these are the pairing constraints (44) and (45), and the infeasible path constraints (49). Also all families of valid inequalities that are only used to strengthen the model are added in a cutting-plane fashion. Separation algorithms will check the current solution for violations of the omitted constraints and the above defined valid inequalities. In case of omitted constraints, the separation procedures have to correspond to exact methods such that, if one of the omitted constraints is violated, we can guarantee that it will be identified. In case of additional valid inequalities, which are not needed to ensure feasibility, heuristic separation procedures can be employed.

In the case where at least one violated constraint is detected by the separation procedures, it is added to the current LP in the form of a cut and the updated LP is solved again. This process is repeated until the separation procedures fail to detect additional violated constraints. In this case, the optimal solution to the LP relaxation of the original MIP has been found. If this solution is integer, the optimal solution to the original MIP has been identified. Otherwise, the problem is decomposed into two new problems. As in branch-and-bound, by branching on a variable that is associated with a fractional value in the current solution. The following branching rules are applied. In case of 3indexBC, as in Cordeau (2006), additional artificial variables are introduced that handle the assignment of requests to vehicles. Branching priority is given to these variables. Both 3indexBC and 2indexBC branch on the variable furthest away from the nearest integer. Let us assume that this variable is x1,3 (e.g., x1,3 = 0.4) in 2indexBC. Branching refers to the generation of two child nodes. At one of the child nodes an LP is built with an additional constraint setting a lower bound on the chosen variable. This lower bound is equal to the fractional value of the variable ceiled to the next integer, e.g. x1,3 ⩾ 1. At the second child node another LP is built with an additional constraint setting an upper bound on the chosen variable. This upper bound is equal to the fractional value of the chosen variable floored to the next integer, e.g. x1,3 ⩽ 0. Then, each of the two LPs is solved in the same way as the first LP relaxation; cuts are added until no more violated inequalities can be detected. In case the optimal solution is not integer, the child node serves as a new parent node for two new child nodes in the tree. The optimal solution to the original MIP is the best of the first two thus recursively solved problems. For further details and additional references on the branch-and-cut method in general we refer the interested reader to Naddef and Rinaldi (2002).

To accelerate the solution process several preprocessing steps can be performed prior to starting the optimization procedure. These steps refer to graph pruning and time window tightening techniques. Furthermore, some “easy” cuts can be generated in advance to strengthen the formulation. In our case these are all strengthened bounds on time and load variables; and, in order to break symmetry, in case of the 3-index formulation, variable fixing techniques can be employed; e.g. in the homogeneous fleet case, those requests that cannot be served by the same vehicle are fixed to different vehicles. Finally, before applying the branch-and-cut algorithms upper bounds are computed by means of an adapted version of a state-of-the-art heuristic method for the standard DARP. These bounds serve to prune the branch-and-bound tree; whenever the solution value at some node exceeds that of the upper bound, this node and all its child nodes can be excluded from the search. Our branch-and-cut implementations are widely based on the work of Cordeau (2006) and Ropke et al. (2007). The employed preprocessing steps, the different separation procedures, and the computation of upper bounds are briefly described in the following.

Initialization

In a first step, preprocessing procedures, using graph pruning and time window tightening techniques as described by Cordeau (2006), adapted to the HDARP, are applied. Then an initial pool of inequalities is generated. All inequalities part of this pool are enumerated exhaustively at every node of the branch-and-bound tree. It contains all cuts part of the initial pool of inequalities described in Cordeau (2006): subtour elimination constraints for ∣S∣ = 2 (see Section 3.2), the following four particular cases of precedence constraints,generalized order constraints for h = 2 and ∣U1∣ = ∣U2∣ = 2 (see Section 3.3), and infeasible path constraints of the form x + x + x ⩽ 1 if this path violates the ride time constraint of request i. Furthermore, four inequalities for each pair of incompatible users i,j ∈ P and every node l ∈ P ∪ D are added to the pool:

Finally, in case of 3indexBC, variable fixing methods are used. If the vehicle fleet of an instance is homogeneous, the same variable fixing procedures as in Cordeau (2006) are applied. If the fleet is heterogeneous, variables can be fixed as follows. If only one vehicle with resource 2 or 3 (stretcher or wheelchair) is available, all users demanding these transportation modes can be fixed to the corresponding vehicle. If some vehicle does not provide these resources all users demanding them can be forbidden on this vehicle. Furthermore, if two users are identified as incompatible for one vehicle, a constraint guaranteeing that only one of the two can use this vehicle is appended to the model.

Separation heuristics

In both algorithms, violated subtour elimination (see Section 3.2) and generalized order constraints (see Section 3.3) are separated by means of several (meta)heuristics as described in Cordeau (2006). For the separation of strengthened infeasible path inequalities (see Section 3.4), an enumerative procedure as described in Ropke et al. (2007) is used. Because of the heterogeneous fleet requirements, these constraints are checked and generated for each vehicle in turn in case of the 3-index formulation. The same applies to violated fork constraints which are also determined by means of enumeration procedures (Ropke et al., 2007). In case of the 2-index formulation, only time related infeasibilities are considered since these are homogeneous across all vehicles.

In 3indexBC, as in Cordeau (2006), for the separation of rounded capacity inequalities (see Section 3.5), a tabu search algorithm is employed. To detect violated strengthened capacity cuts, the heuristics as described by Ropke et al. (2007) are used, i.e. a randomized construction heuristic and another tabu search algorithm. Note that, in case of 3indexBC, the separation procedures briefly described above are only evoked if all assignment variables are integer, i.e. no vertex is partially assigned to more than one vehicle.

For 2indexBC additional separation procedures need to be devised. Precedence and pairing constraints, given in (44) are separated similar to Ropke et al. (2007) (the implementation of Goldberg’s algorithm for solving the individual max-flow problems is used). The pairing constraints in the 2-index model for the introduced artificial depots are separated by a similar procedure. Every origin depot together with all other destination depots are in turn used as the super source node and its corresponding destination depot together with all other origin depots as the super sink node. If the max flow on the network (arc capacities are defined by the current values of the x variables) is smaller than one, the according depot pairing constraint, given in (45), is generated.

For the detection of infeasible paths regarding heterogeneous vehicles and resources, violating a capacity restriction in the 2-index formulation, another enumerative procedure is employed: every path starting at one of the start depots is extended on all arcs with some flow (x > 0) and checked for capacity violations. If the capacity is exceeded at some vertex i, the corresponding inequality (49) is generated. Path extension ends as soon as one of the destination depots has been reached.

As in Ropke et al. (2007), to separate reachability cuts in 2indexBC, in a first step, for each vertex i ∈ P ∪ D the arc sets and are computed. Based on these arc sets conflicting vertex sets are identified (if there is no path that contains a certain vertex j leading from one of the start depots to i, or from i to one of the destination depots, i and j are conflicting vertices). After having determined all conflicting pairs, conflicting sets of larger cardinality are generated on the basis of conflicting pairs. Obviously, at most sets with a cardinality of m are considered. If sets of higher cardinality existed the respective problem instance would be infeasible. At each fractional solution encountered at some node of the branch and bound tree, each conflicting set is considered in turn. In order to identify violated inequalities (86), a maximum flow problem between D and is solved using only arcs . In order to do so, an artificial source node that is connected to each vertex ∈D with arcs of infinite capacity, and an artificial sink node connecting all vertices with arcs of infinite capacity to this sink have to be introduced. If the minimum cut’s capacity obtained is smaller than , the according reachability cut is generated. Similarly, to identify violated inequalities (87), a max-flow problem between and D is solved using only arcs .

Heuristic upper bounds

In order to accelerate the optimization process initial upper bounds are calculated by means of an adapted version of a state-of-the-art VNS (Mladenovic and Hansen, 1997) designed for the standard DARP in Parragh (2009, Chapter 3). The VNS developed in Parragh (2009, Chapter 3) departs from a possibly infeasible initial solution which is constructed taking time window and spatial closeness considerations into account. The thus constructed initial solution constitutes the first incumbent solution s. Then, in every iteration a random solution s′ is generated in the current neighborhood of s (shaking). In the shaking phase a neighborhood operator belonging to one of four different neighborhood classes is applied. The neighborhood classes are denoted as swap (S), move (M), chain (C), and zero split (Z). Throughout the search they are traversed in alternating order and increasing size. Given that s is a promising solution, s undergoes local search-based improvement (see Parragh, 2009; Parragh et al., 2010). If s′ meets the acceptance criteria, it replaces s and the search continues with the first (smallest) neighborhood. If s′′ does not meet the acceptance criteria, s is not replaced and the next (larger) neighborhood is used in the subsequent iteration (move or not). The acceptance criteria employed refer to simulated annealing type (Kirkpatrick et al., 1983) rules. This means that also deteriorating solutions may become incumbent solutions with a certain probability. Whenever the last (largest) neighborhood is attained, the search continues with the first neighborhood. This is repeated until some stopping criterion is met. Furthermore, intermediate infeasible solutions are allowed and infeasibilities are penalized. At the end, the best feasible solution s encountered during the search is returned. Finally, in order to reduce the search space, the graph pruning and time window tightening techniques, described by Cordeau, 2006, are applied prior to starting the optimization procedure.

We adapted this state-of-the-art method in the following ways. First, a new evaluation function that accommodates the characteristics of the HDARP has been devised. Besides routing costs and constraint violations, also the total waiting time with passenger aboard has to be considered:

The routing costs of a solution s are given by ĉ(s), the total vehicle waiting time with passengers aboard by , weighted by ρ. The terms refer to the different resource violations r ∈ R. Each of these is penalized by a separate penalization parameter . Furthermore, duration violation , time window violation , and ride time violation are penalized by the self-adjusting terms , and , respectively (for details see Parragh et al., 2010). A solution s can only become a new best solution s if for all r ∈ R.

Second, a regret insertion procedure for the construction of the initial solution has been designed. The first m requests are each assigned to a different vehicle. All subsequent requests are inserted as follows. A regret value for every request is calculated: the best insertion positions for each request on each route are determined and the corresponding evaluation function values are calculated; the regret value for a given request corresponds to the difference in terms of evaluation function value between its two best insertion positions on two different routes. In every iteration the request associated with the largest regret value is inserted at its best position.

Third, only 13 different neighborhoods are considered (S1–M1–C1–S2–M2–C2–S3–M3–C3–S4–M4–C4–Z), where the number in addition to the neighborhood abbreviation indicates the neighborhood size, see also Parragh (2009, Chapter 3)). Initial tests with 19 neighborhoods showed that, as in the context of the multi-objective DARP (Parragh et al., 2009b), 13 neighborhoods are sufficient.

Further, the above mentioned graph pruning techniques are only applied after the first feasible solution has been found. And last but not least, as stopping criterion a limit of 5 × 105 iterations is employed.

Computational experiments

All programs are implemented in C++. In the branch-and-cut algorithms CPLEX 11.0 together with Concert Technology 2.5 are employed. All experiments are carried out on a 3.2 GHz Pentium D computer with a memory of 4 GB. All solution procedures are tested on three artificial data sets. These are based on an existing data set from the literature, enriched with the different real world characteristics described above. In the following, first, the characteristics of the generated test instances are described. Then, the obtained results are discussed. First the results obtained by means of the two branch-and-cut algorithms are summarized; then, the solutions obtained by means of the adapted VNS are presented. Both procedures are first tested with ρ = 0 (waiting is not penalized) and then with ρ = 100, as suggested by the ARC, such that waiting on board a vehicle is avoided.

Test instances

The test instances used are based on the 12 “A” instances proposed by Cordeau (2006) for the standard DARP, containing 12 instances with 2–4 vehicles and 16–48 requests. In all instances time window length l − e = 15 min, maximum user ride time  min, and service time d = 3 min for all users. For each instance three instances with different degrees of heterogeneity are generated.

Heterogeneous users are introduced based on the probabilities given in Table 2. Every transportation request consists of at most one patient. In data set “U” only seated passengers and no accompanying persons are considered. In data set “E” half of the patients are considered to be seated, 25% to be on a stretcher, and 25% to be in a wheelchair; 10% are assumed to be accompanied by someone. Eventually, in data set “I” users are transformed into seated patients, patients on a stretcher, and wheelchair passengers according to the true distribution across all static transports carried out by the ARC in the city of Graz. Furthermore, half of them are assumed to be accompanied by someone.

Table 2

Test instances.

Data set	Probability for patient to be			Probability
	Seated	On stretcher	In wheelchair	For AP	Fleet
U	1.00	0.00	0.00	0.00	hom. (T0)
E	0.50	0.25	0.25	0.10	hom. (T2)
I	0.83	0.11	0.06	0.50	het. (T1, T2)

AP = accompanying person, hom. = homogeneous, het. = heterogeneous.

T0: 3 patient seats.

T1: 1 staff seat, 6 patient seats, 1 wheelchair place.

T2: 2 staff seats, 1 patient seat, 1 stretcher, 1 wheelchair place.

In each data set a different fleet configuration is used. In data set “U” a homogeneous fleet setting with only vehicles of type T0 is considered. Vehicle type T0 provides space for three seated passengers. In data set “E” a homogeneous vehicle fleet consisting of T2 vehicles, providing 2 staff seats, 1 patient seat, 1 stretcher, and 1 wheelchair place, is used; for data set “I” a heterogeneous vehicle fleet is generated. Here the original number of vehicles is randomly divided into T1 and T2 vehicles such that at least one vehicle of each type is available. Vehicle types T1 and T2 are derived from data provided by the ARC regarding their vehicle fleet.

Branch-and-cut results

Table 3 provides the results obtained by 3indexBC and 2indexBC, ignoring the penalization option regarding waiting with passengers aboard (ρ = 0). It contains the following information. Columns “z” give the proved optimal objective value (marked by an asterisk) or the best feasible solution found. Furthermore, the best lower bound identified within the time limit (bestLB), run times in seconds (CPU), the number of nodes visited as well as the number of cuts added are provided. The number of cuts refers to the number of user cuts generated during the optimization process. Cuts generated by CPLEX are not counted. In the last column (heurUB) the employed initial upper bounds, computed by means of one run of the adapted VNS, can be found. For all experiments reported in this section a run time limit of 4 h was used.

Table 3

3indexBC vs. 2indexBC (ρ = 0).

	3indexBC					2indexBC
	z	bestLB	CPUa	Nodes	Cuts	z	bestLB	CPUa	Nodes	Cuts	heurUB
U
a2-16	294.25∗	294.25	44.05	58	152	294.25∗	294.25	1.13	0	44	294.25
a2-20	344.83∗	344.80	1674.30	6751	2531	344.83∗	344.83	2.59	0	112	344.83
a2-24		415.17	14492.70	13632	2391	431.12∗	431.12	8.54	0	161	431.12
a3-18	300.48∗	300.45	4376.70	16350	933	300.48∗	300.48	4.55	0	249	300.48
a3-24		334	14516.40	7380	595	344.83∗	344.83	7.62	0	422	347.42
a3-30		467.25	14490	3953	433	494.85∗	494.85	9.83	0	575	494.85
a3-36		553.05	14570.80	3917	351	583.19∗	583.19	105.05	0	1099	584.44
a4-16	282.68∗	282.67	319.42	284	241	282.68∗	282.68	5.61	0	430	282.68
a4-24	375.02∗	375.02	709.89	187	240	375.02∗	375.02	5.60	0	347	378.13
a4-32		432.78	14520.60	2084	291	485.50∗	485.50	30.67	0	1056	487.81
a4-40		502.74	14468.60	88	339	557.69∗	557.63	8328.46	9723	6293	582.26
a4-48		571.27	14523.20	48	340	668.82	664.64	14542.60	5076	5164	709.47
U¯		406.12	9058.89	4561	736		429.92	1921.02	1233	1329	436.48
E
a2-16	331.16∗	331.16	37.54	14	101	331.16∗	331.13	284.17	4908	1461	331.16
a2-20	347.03∗	347.03	441.26	1523	817	347.03∗	347.03	8.06	0	150	347.03
a2-24		421.63	14537.70	7198	1255	450.25∗	450.21	891.24	5743	1298	450.25
a3-18	300.63∗	300.63	4412.52	13053	612	300.63∗	300.63	4.28	0	304	300.63
a3-24		340.65	14528.70	6721	505	344.91∗	344.91	10.15	0	505	346.22
a3-30		484.26	14514.80	3199	445	500.58∗	500.53	1608.63	4898	2600	500.58
a3-36		565.55	14554.10	2165	312	583.19∗	583.19	101.40	0	1133	585.94
a4-16	285.99∗	285.99	194.69	31	155	285.99∗	285.99	759.27	7664	1832	285.99
a4-24	383.84∗	383.84	1321.24	112	383	383.84	380.48	14471.60	19808	22083	390.87
a4-32		459.87	14499.90	870	562		488.14	14494.10	4733	12610	508.51
a4-40		534.35	14716.30	43	516		558.09	14518.50	2390	17918	609.08
a4-48		593.36	14625.30	0	635		665.02	14539.50	1476	11919	704.07
E¯		420.69	9032.00	2911	525		436.28	5140.08	4302	6151	446.69
I
a2-16	294.25∗	294.25	92.17	59	148	294.25∗	294.25	0.88	0	25	294.25
a2-20	355.74∗	355.74	10142.80	42432	3728	355.74∗	355.71	115.44	1103	428	355.74
a2-24		421.75	14553.30	13825	1324	431.12∗	431.12	7.06	0	236	431.12
a3-18		297.02	14538.70	17172	651	302.17∗	302.15	30.78	187	515	302.17
a3-24		332.49	14528.50	5570	354	344.83∗	344.83	6.78	0	436	344.83
a3-30		472.35	14490.80	2896	481	494.85∗	494.85	9.50	0	574	494.85
a3-36		566.23	14503.70	1173	263	618.63	592.69	14500.10	8742	17455	625.90
a4-16	299.05∗	299.05	360.75	804	296	299.05∗	299.05	36.51	35	606	299.05
a4-24	375.07	372.07	14499.20	4221	244	375.02∗	375.02	5.23	0	391	375.07
a4-32		419.60	14497.50	1073	225	486.93∗	486.93	40.63	0	1093	498.48
a4-40		499.03	14525.60	23	188	557.69∗	557.63	2600.02	2093	4692	583.53
a4-48		567.51	14630.20	32	276	678.59	663.30	14498.00	1820	9526	708.68
I¯		408.09	11780.27	7440	682		433.13	2654.25	1165	2998	442.81
UEI¯		411.64	9957.05	4971	648		433.11	3238.45	2233	3493	441.99

Run times in seconds.

When comparing the results obtained by the two branch-and-cut algorithms, 2indexBC clearly outperforms 3indexBC. 2indexBC finds the optimal solution to 29 out of 36 test instances within the 4-h time limit, while 3indexBC only finds 13. Only in one case (instance a4-24 of data set “E”), 3indexBC finds the optimal solution and proves it to be optimal, while 2indexBC does not. In 2indexBC, on average five times more cuts are generated but only half of the number of nodes are explored. Taking a closer look at the number of nodes explored in the enumeration tree, in 2indexBC, in many cases the optimal solution is already found at the root node; this is never the case for 3indexBC.

The number of cuts added during the optimization process is further analyzed. In Table 4 the percentage share of each class of cuts, as presented above, is given. While the 3-index formulation could be solved by a commercial MIP solver such as CPLEX directly, in case of the 2-index formulation constraint sets (44), (45) and (49) have to be added in a cutting plane fashion. Their percentage share is given in italic letters in Table 4. Almost 70% of the cuts generated in case of 2indexBC correspond to precedence and pairing cuts (44) and (45). Also parts of the infeasible path cuts are needed to ensure feasibility, i.e. those related to loading restrictions (49). These amount to about 2% of the total number of cuts added. Hence, only approximately 30% of the number of cuts given in Table 3 are based on additional valid inequalities.

Table 4

Percentage share of each class of cuts on the total number of cuts (ρ = 0).

	3indexBC				2indexBC
	Subtour elim.	Gen. order	Infeasible path	Capacity inequ.	subtour elim.	Gen. order	Infeasible path	Reachability	Pairing/precedence	Infeasible path (load)
U¯	11.40	6.12	57.18	17.64	0.93	0.19	6.14	18.27	73.61	0.03
E¯	15.64	8.57	31.64	35.27	0.28	0.06	28.83	7.43	59.73	3.40
I¯	13.18	6.40	68.70	3.84	0.27	0.10	18.71	15.48	62.97	2.03
UEI¯	13.4	7.03	52.51	18.91	0.49	0.12	17.89	13.73	65.44	1.82

Elim. = elimination, gen. = generalized, inequ. = inequalities.

Comparing those families of inequalities that are used in both 3indexBC and 2indexBC, infeasible path cuts make up the largest portion of cuts. In case of 3indexBC also a considerable number of subtour elimination constraints and capacity inequalities are generated. In 2indexBC, excluding those cuts needed to ensure feasibility, reachability cuts form the second largest portion of cuts.

Analyzing the number of cuts added on a data set level, in 3indexBC, the percentage share of the capacity cuts varies considerably. Data set “U” corresponds to the original homogeneous data set; only one type of passenger and only one type of vehicle are considered. Here quite a lot of capacity cuts are identified. In case of data set “E”, even more capacity cuts are separated; heterogeneous users but only one type of vehicle are considered. Fewer capacity cuts are generated in case of data set “I”. The reason is that, in addition to heterogeneous users, also heterogeneous vehicles are considered and thus the per resource capacity upper bound across all vehicles deteriorates.

The conclusion that can be drawn from these results is that, as expected, the 2-index-based formulation leads to a more efficient branch-and-cut algorithm. Therefore, in the following experiment only 2indexBC is considered.

Table 5 provides the results obtained with 2indexBC and ρ = 100. This means that waiting with passengers aboard is penalized rather severely. All results thus obtained do not contain any waiting time with passengers aboard. Penalizing waiting with passengers aboard makes the problem slightly more difficult to solve; instead of 29 out of 36 instances, 28 out of 36 instances can be solved to optimality within the 4-h time limit.

Table 5

2indexBC (ρ = 100).

	2indexBC
	z	bestLB	CPUa	Nodes	Cuts	heurUB
U
a2-16	300.17∗	300.17	0.92	0	63	300.17
a2-20	345.74∗	345.74	2.43	0	117	345.74
a2-24	448.07∗	448.07	15.67	0	147	448.07
a3-18	316.78∗	316.78	3.59	0	204	316.78
a3-24	346.97∗	346.97	11.55	0	516	347.37
a3-30	501.68∗	501.68	23.46	0	549	502.39
a3-36	598.53∗	598.53	21.55	0	646	599.55
a4-16	282.68∗	282.68	5.54	0	357	282.68
a4-24	386.38∗	386.38	5.73	0	398	386.38
a4-32	493.15∗	493.15	83.48	0	1200	493.15
a4-40	557.94∗	557.90	98.35	0	1204	560.31
a4-48	707.89	697.88	14467.20	4514	4132	711.12
U¯		439.66	1228.29	376	794	441.14
E
a2-16	331.16∗	331.13	92.46	1406	657	331.16
a2-20	347.89∗	347.89	4.76	0	77	347.89
a2-24	461.93∗	461.89	231.05	1310	582	462.82
a3-18	316.78∗	316.78	5.99	0	224	316.78
a3-24	347.05∗	347.05	17.32	0	637	347.45
a3-30	501.68∗	501.68	27.00	0	617	504.15
a3-36	604.35∗	604.30	390.94	447	1820	606.08
a4-16	291.55	288.70	14479.30	76702	10520	291.55
a4-24	386.38∗	386.38	5.66	0	336	386.38
a4-32		491.66	14490.20	4836	9864	507.73
a4-40		560.00	14507.70	2387	21330	590.19
a4-48		697.28	14498.10	2355	7316	717.50
E¯		444.56	4895.87	7454	4498	450.81
I
a2-16	300.17∗	300.17	0.99	0	39	300.17
a2-20	356.64∗	356.61	28.06	42	205	356.64
a2-24	450.37∗	450.37	30.04	21	251	450.37
a3-18	318.47∗	318.47	17.17	0	311	318.47
a3-24	346.97∗	346.97	16.57	0	644	347.37
a3-30	501.68∗	501.68	10.45	0	478	501.68
a3-36	631.12	612.35	14511.20	9386	18371	633.29
a4-16	301.81∗	301.78	198.41	1564	1244	302.23
a4-24	386.38∗	386.38	7.35	0	396	386.38
a4-32	494.59	493.58	14488.90	18973	9994	499.06
a4-40	559.45∗	559.42	223.67	51	1589	561.35
a4-48		698.39	14516.60	2125	8138	721.40
I¯		443.85	3670.78	2680	3472	448.20
UEI¯		442.69	3264.98	3503	2921	446.72

Run times in seconds.

Heuristic results

Table 6 gives the results obtained by means of the adapted VNS for five random runs with ρ = 0 and ρ = 100. For each data set the average value (avg.), the best value (best), the deviation from the optimal solution (%), where known, the run time in seconds (CPU), and the total waiting time with passengers aboard (waiting), averaged over five runs, are given. The adapted VNS is able to cope with both settings. The largest deviation from the optimal solution in case of ρ = 0 is 1.41%, taking average values over five runs. When considering the best values out of these five runs, the maximum gap is equal to 0.74%. In most cases the optimal solution is found. The largest deviation in case of ρ = 100, taking average values over five runs, is 1.06%. The largest gap, taking the best solution values out of these five runs, is equal to 0.27%. This indicates that for the adapted VNS, penalizing waiting with passengers aboard makes the problem slightly easier to solve. In both versions run times are low, less than 2 min on average. When comparing the amount of waiting time contained in the average solutions, setting ρ = 0 results in, on average, 18.98 min of waiting, while setting ρ = 100 results in no waiting for the passengers. Only for one instance (a4-32 data set “I”) one of the five random runs yielded a solution with passenger waiting time (0.08 min). Only slight adaptations were thus necessary in order to adapt the previously developed VNS for the standard DARP to the new problem version and to yield high quality solutions.

Table 6

VNS (five runs).

	ρ = 0						ρ = 100
	Avg.	%	Best	%	CPUa	Waiting	Avg.	%	Best	%	CPUa	Waiting
U
a2-16	294.25	0.00	294.25	0.00	68.20	0.57	300.17	0.00	300.17	0.00	70.20	0.00
a2-20	344.83	0.00	344.83	0.00	133.80	13.25	345.74	0.00	345.74	0.00	128.80	0.00
a2-24	431.12	0.00	431.12	0.00	187.80	22.21	448.07	0.00	448.07	0.00	202.40	0.00
a3-18	300.48	0.00	300.48	0.00	45.40	12.19	316.78	0.00	316.78	0.00	47.40	0.00
a3-24	345.26	0.12	344.83	0.00	86.80	14.47	347.37	0.11	347.37	0.11	86.60	0.00
a3-30	495.11	0.05	494.85	0.00	105.60	15.67	502.46	0.16	501.68	0.00	111.20	0.00
a3-36	583.89	0.12	583.30	0.02	162.60	55.06	602.23	0.62	599.02	0.08	175.20	0.00
a4-16	282.68	0.00	282.68	0.00	26.00	0.00	282.68	0.00	282.68	0.00	26.80	0.00
a4-24	375.04	0.00	375.02	0.00	50.80	30.90	386.92	0.14	386.38	0.00	51.20	0.00
a4-32	488.27	0.57	486.88	0.28	86.00	11.79	496.05	0.59	493.15	0.00	84.60	0.00
a4-40	565.58	1.41	561.80	0.74	130.60	0.75	563.87	1.06	559.45	0.27	130.00	0.00
a4-48	680.98		673.64		253.80	50.57	717.58		711.12		0.00	0.00
U¯	432.29		431.14		111.45	18.95	442.49		440.97		92.87	0.00
E
a2-16	331.16	0.00	331.16	0.00	65.60	0.00	331.16	0.00	331.16	0.00	67.20	0.00
a2-20	347.03	0.00	347.03	0.00	120.00	13.25	347.89	0.00	347.89	0.00	111.80	0.00
a2-24	450.25	0.00	450.25	0.00	160.40	16.03	462.57	0.14	461.93	0.00	172.60	0.00
a3-18	300.63	0.00	300.63	0.00	47.60	9.27	316.78	0.00	316.78	0.00	49.00	0.00
a3-24	345.59	0.20	344.91	0.00	76.20	14.47	347.37	0.09	347.05	0.00	80.60	0.00
a3-30	501.41	0.17	500.58	0.00	107.60	7.24	503.70	0.40	501.68	0.00	106.40	0.00
a3-36	583.79	0.10	583.19	0.00	161.60	58.43	607.10	0.46	605.46	0.18	161.80	0.00
a4-16	285.99	0.00	285.99	0.00	25.00	20.00	291.55	0.00	291.55	0.00	24.40	0.00
a4-24	384.03	0.05	383.84	0.00	52.60	10.88	386.64	0.07	386.38	0.00	51.00	0.00
a4-32	504.79		502.52		83.00	15.88	509.76		507.72		81.00	0.00
a4-40	588.40		585.64		121.00	12.94	595.81		590.19		121.00	0.00
a4-48	681.80		675.37		252.20	49.48	717.95		715.62		285.00	0.00
E¯	442.07		440.93		106.07	18.99	451.52		450.28		109.32	0.00
I
a2-16	294.25	0.00	294.25	0.00	68.40	0.57	300.36	0.06	300.17	0.00	71.40	0.00
a2-20	355.74	0.00	355.74	0.00	141.80	13.25	356.64	0.00	356.64	0.00	137.40	0.00
a2-24	431.12	0.00	431.12	0.00	211.00	22.21	450.37	0.00	450.37	0.00	207.20	0.00
a3-18	302.17	0.00	302.17	0.00	47.20	12.19	318.62	0.05	318.47	0.00	50.20	0.00
a3-24	344.99	0.05	344.83	0.00	83.60	14.47	347.55	0.17	347.37	0.11	87.20	0.00
a3-30	495.13	0.06	494.85	0.00	106.80	15.67	501.68	0.00	501.68	0.00	111.60	0.00
a3-36	619.64		618.58		170.60	43.11	630.61		627.39		207.80	0.00
a4-16	299.05	0.00	299.05	0.00	27.00	20.00	302.06	0.08	301.81	0.00	26.40	0.00
a4-24	376.19	0.31	375.07	0.01	51.60	26.81	386.89	0.13	386.38	0.00	51.40	0.00
a4-32	488.64	0.35	486.93	0.00	88.00	10.08	498.99		497.07		84.80	0.02
a4-40	563.34	1.01	561.35	0.66	132.20	2.36	563.07	0.65	561.35	0.34	132.20	0.00
a4-48	687.44		680.43		262.40	47.17	717.98		713.21		303.60	0.00
I¯	438.14		437.03		115.88	18.99	447.90		446.83		122.60	0.00
UEI¯	437.50		436.37		111.13	18.98	447.31		446.03		108.26	0.00

Run times in seconds.

When comparing average routing costs for ρ = 0 and ρ = 100 a cost increase of about 2.5% can be observed. When comparing best values out of five runs, the largest cost increase amounts to almost 6% (instance a4-48 of data set “E”). The implications from a company perspective are that only moderate cost increases will lead to improved solutions from a customer perspective. What remains to be seen is how this relationship translates into larger problem instances. Possibly, the difference between the two extremes, only minimizing costs and avoiding user waiting time when aboard a vehicle, will increase.

Conclusions

In this article heterogeneous vehicles and users have been introduced into state-of-the-art models and algorithms for the standard dial-a-ride problem. Two problem formulations have been proposed: a 3-index and a more sophisticated 2-index formulation that integrates the given heterogeneous fleet requirements in a clever way. When used in a branch-and-cut framework, the 2-index formulation clearly outperforms the 3-index-based one in all three data sets considered. A variable neighborhood search heuristic, originally developed for the standard dial-a-ride problem, has also been adapted to this new problem version. High quality solutions are computed within short computation times. These results suggest that the application of the proposed method is also suitable for larger (real world) instances. Furthermore, the impact of penalizing vehicle waiting time with passengers aboard has been investigated. Heuristic results for the two extreme settings have been compared: on the one hand, only considering the cost part in the objective function, on the other hand, setting the penalization term to a rather high value. The latter setting results in solutions without user waiting time when aboard a vehicle. When comparing the results, the latter setting yields average routing costs that are only 2.5% higher than in the minimum cost setting. The maximum increase amounts to about 6%. This indicates that from a company perspective, avoiding waiting time with users aboard a vehicle does not lead to a significant cost increase.