Adaptive guided ejection search for pickup and delivery with time windows

Abstract

The pickup and delivery problem with time windows is an NP-hard discrete optimization problem with two objectives—to minimize the fleet serving transportation requests, and to minimize the distance traveled during this service. Although there exist exact algorithms for tackling this problem, they are still difficult to apply in massively large practical scheduling scenarios due to their time complexities. Hence, the approximate methods became the main stream of research in this field. In this paper, we propose an adaptive guided ejection search algorithm for solving the pickup and delivery with time windows. The pivotal part of this technique is the pre-processing step, in which the instance characteristics concerning its underlying structure are extracted in the clustering and histogram-based analyses. Then, the k-nearest neighbor algorithm is applied to classify the instance to an appropriate class. Finally, the most suitable variant of our enhanced guided ejection search algorithm is adaptively chosen for solving this instance based on the classification outcome. An extensive experimental study performed on the full Li and Lim’s benchmark (encompassing 354 problem instances belonging to 6 classes) revealed that our pre-processing allows for achieving very high classification accuracy, thus for selecting the best variant of the enhanced guided ejection search.

Keywords

Adaptation clustering guided ejection search PDPTW

1 Introduction

Route scheduling is an important problem due to its wide practical applicability. The feasibility of routing schedules is affected by numerous factors, including the time slots in which customers should be served, the maximum truck capacities, the precedence constraints imposed on the order of visiting travel points, potential dynamic changes of the route network, and many others [1]. There exist a plethora of different variants of routing problems which reflect these scheduling circumstances. In the multiple traveling salesman problem (mTSP), multiple salesmen can serve customers. In the capacitated vehicle routing problem (CVRP), the customers define their non-negative demands which need to be satisfied using vehicles of a given maximum capacity. The vehicle routing problem with time windows (VRPTW) defines the time slots, within which the service of each customer should be started (however, it may be finished after closing this time window). Other practical VRP variants encompass heterogeneous fleet VRP, the VRP with stochastic demands, the periodic VRP, and more [2].

Solving the pickup and delivery problem with time windows (PDPTW) consists in finding a feasible schedule for serving a set of transportation requests (combinations of the pickup and delivery requests) using homogeneous trucks [3]. In each request (which must be served within exactly one route), goods are picked up from one location, and delivered to another one. These travel points are scattered around the map and must be visited within the corresponding time windows (the service may finish after closing this time slot). The pickup location must be visited before the delivery one by the serving truck, and the capacity of all vehicles must not be exceeded. All vehicles start and finish their routes in a single depot which also specifies its own time window.

The PDPTW is a hierarchical-objective optimization problem [4]: the main objective is to minimize the fleet size, whereas the second one is to optimize the travel distance. Due to numerous applications of the VRPs (e.g., rail distribution, food delivery or industrial waste collection), they attracted research attention [5]. Since the problem is NP-hard, exact approaches are applicable for relatively small-scale scenarios because of their unacceptable execution times [3].

Approximate algorithms, which do not ensure retrieving exact solutions but are capable of obtaining nearly-optimal schedules in short time, are intensively being developed to tackle the PDPTW. They include genetic algorithms, various neighborhood searches, simulated annealing, and many others [6]. Guided ejection searches (GESes) were proven very efficient (and outperforming the state of the art) in this task [7]. In our very recent work [8], we proposed several improvements and showed that they significantly affect the final routing schedules retrieved using the GES, as well as its behavior and capabilities. Although this enhanced GES was able to elaborate very high-quality solutions, it was unclear how to select an appropriate algorithm variant which suits the instance to be solved as best as possible. This issue is tackled in this work.

The paper is organized as follows. Section 1.1 gives the contribution of this paper. The PDPTW is formulated in Section 2. Section 3 reviews the current state of the art in solving the PDPTW. In Section 4, we present our adaptive guided ejection search algorithm for tackling the PDPTW, with an emphasis put on the pre-processing procedure which is the main contribution of this work. Section 5 discusses the experimental results. The paper is concluded in Section 6, which also gives an outlook to the future work.

1.1 Contribution

As already mentioned, selecting an appropriate GES variant for solving a given PDPTW instance is a challenging and demanding task. In [8], we showed that applying certain algorithm components for some classes of problems (e.g., encompassing clustered customers) is very beneficial and helps boost the convergence capabilities of the algorithm (the differences were shown to be statistically important). However, the characteristics of a test case are not known in advance and should be carefully extracted before the execution.

This paper extends our previous work presented at the Asian Conference on Intelligent Information and Database Systems (ACIIDS 2016) [8], held in Da Nang, Vietnam (on March 14-16, 2016). Here, we propose the pre-processing procedure aimed at extracting discriminative features of test instances, and classifying them into appropriate classes. In general, six classes of PDPTW tests are distinguished. They include problems with:

clustered customers with tight time windows and small vehicle capacities (c1 class),

clustered customers with wide time windows and large vehicle capacities (c2),

randomly scattered customers with tight time windows and small vehicle capacities (r1),

randomly scattered customers with wide time windows and large vehicle capacities (r2),

mixed (clustered and randomized) customers with tight time windows and small vehicle capacities (rc1), and

mixed (clustered and randomized) customers with wide time windows and large vehicle capacities (rc2).

The pre-processing encompasses the analysis of the customer locations, and the expected minimum numbers of trucks serving the requests. First, we calculate 27 various clustering validity indices (using four different aggregation methods), and average their responses for determining the number of clusters which is the most appropriate for the underlying instance structure. This procedure is complemented with the affinity propagation clustering [9]—it is used to determine the number of clusters independently. To fully exploit the information about the geometrical structure of the instance, we perform the histogram-based analysis (at various scales), and quantify the two-dimensional histograms of customer locations. Finally, we investigate the theoretically minimum numbers of trucks which have to be utilized in a solution to serve all requests, without exceeding the truck capacities.

All of the aforementioned extracted features are combined into the feature vectors, which are classified using the k-nearest neighbor algorithm. An extensive experimental study, performed using the full Li and Lim’s benchmark set containing 354 test instances, revealed that our pre-processing step allows for achieving very high classification accuracy. Therefore, it allows for retrieving the best-fitted guided ejection search variant before its execution—the algorithm is adaptively tailored for the problem at hand.

2 Problem formulation

The PDPTW is a problem of serving N requests by K vehicles. More formally, the PDPTW is defined on a directed graph G = (V, E), with a set V of C + 1 vertices. The vertices v_i, i ∈ {1, …, C}, represent the customers, whereas v₀ denotes a single depot. The edges E = {(v_i, v_i+1) |v_i, v_i+1 ∈ V, v_i ≠ v_i+1} represent the travel connections. The travel costs c_i,j, i, j ∈ {0, 1, …, C}, i ≠ j, are equal to the distances (in the Euclidean metric) between the travel points.

Each request h_i, i ∈ {0, 1, …, N}, where N = C/2, is a pair of the pickup (origin, P) and delivery (destination, D) operations (also referred to as the jobs [3]), where P∩ D = ∅, and P ∪ D = V \ {v₀}. The pickup occurs before the corresponding delivery, and both must be served in the same route. For each request h_i, the amounts of delivered (q^d (h_i)) and picked up (q^p (h_i)) goods are defined, where q^d (h_i) = - q^p (h_i). A customer v_i defines its own demand (either delivery or pickup), its service time s_i (where s₀ = 0), and time window [e_i, l_i]. The service of each customer should be started within its time window. The fleet is homogenous—the capacity of each truck is Q.

The PDPTW solution σ is a set of routes. Each route r is given as $r = < v_{0}, v_{r (1)}, v_{r (2)}, \dots, v_{n + 1} >$ , and it starts and finishes in the depot (thus v₀ = v_n+1).

2.1 Objectives

The PDPTW is a two-objective discrete optimization problem. The primary objective is to minimize the fleet size K, K ≥ K_min (, and $D = \sum_{i = 1}^{N} q^{d} (h_{i})$ is the total delivery demands).

The second objective is to minimize the total distance (T) traveled by all the vehicles: $T = \sum_{i = 0}^{C} \sum_{j = 0}^{C} \sum_{k = 1}^{K} x_{(i, j, k)} c_{(i, j)} .$ (1) The problem is characterized by 3 decision variables: $x_{(i, j, k)} (i, j \in {0, 1, \dots, C}, i \neq j, \forall k \in K),$ (2) $a_{i} (i \in {0, 1, \dots, C}), and$ (3) $w_{i} (i \in {1, 2, \dots, C}) .$ (4) If the k-th vehicle travels from v_i to v_j (where i ≠ j), then x_(i,j,k) = 1 (0 otherwise). Two other decision variables indicate the arrival and the waiting times at v_i (a_i and w_i, respectively, and a₀ = e₀).

Let σ_A and σ_B be two feasible PDPTW solutions. Then, σ_A is of a higher quality than σ_B, if (K (σ_A) < K (σ_B)) or (K (σ_A) = K (σ_B) and T (σ_A) < T (σ_B)).

2.2 Constraints

The PDPTW constraints may be formally expressed using the following equations: $\begin{matrix} \sum_{k = 1}^{K} \sum_{j = 0, j \neq i}^{C} x_{(i, j, k)} & = & \sum_{k = 1}^{K} \sum_{i = 0, j \neq i}^{C} x_{(i, j, k)} = 1 \\ (\forall i, j \in V), \end{matrix}$ (5) $\begin{matrix} \sum_{j = 1}^{C} x_{(i, j, k)} & = & \sum_{j = 1}^{C} x_{(j, i, k)} = 1 \\ (i = 0, \forall k \in K), \end{matrix}$ (6) $\sum_{j = 1}^{C} \sum_{k = 1}^{K} x_{(i, j, k)} = K (i = 0),$ (7) $\sum_{i = 1}^{C} q^{d} (h_{i}) \sum_{j = 0, j \neq i}^{C} x_{(i, j, k)} \leq Q (\forall k \in K),$ (8) $e_{i} \leq a_{i} + w_{i} \leq l_{i} (\forall i \in V), and$ (9) $max {a_{i} + s_{i} + c_{i, j}, e_{j}} \leq l_{j},$ (10) for all pairs (i, j) such that i ≠ j, and x_(i,j,k) = 1, for a given k. Every customer is visited exactly once (Equation 5), all routes start and finish in the depot (Equation 6), and the fleet size is equal to K (Equation 7). Also, if a truck arrives at a certain travel point, then it departs from the same point. The capacity (Equation 8) and the time window constraints (Equations 9–10) must hold for each route as well. Therefore, the total amount of goods being delivered cannot exceed Q, and the service of each customer must be started before its time window closes.

3 Related literature

State-of-the-art techniques for tackling the PDPTW (and numerous of its variants addressing a plethora of real-life scheduling cases [10]) are divided into exact and approximate approaches. The former algorithms deliver the exact solutions, whereas the latter are aimed at obtaining very high-quality (quite often nearly-optimal) solutions in acceptable time. The exact algorithms were devised for solving relatively small problem instances (up to 30 requests [11]) due to their computation time. Albeit they are still being actively developed, this execution time becomes their bottleneck and it is still difficult to apply them for tackling massively large real-life scheduling scenarios. The exact techniques encompass dynamic programming and column generation methods, branch-and-cut, branch-and-price solvers, and more [12, 13]. The set-partitioning-like integer formulation of the problem was presented by Baldacci et al. [14]—two dual ascent heuristics were combined with the cut-and-column generation to find a dual solution of the linear programming relaxation of this formulation. The exact algorithms were summarized in several surveys [15, 16].

Approximate algorithms for minimizing the number of trucks in the PDPTW include construction and improvement heuristics. The construction (insertion-based) heuristics create solutions from scratch by inserting consecutive requests iteratively into the partial solution (in which a subset of requests is already handled) according to certain criteria, e.g., the maximum cost savings, the minimum additional travel distance, or the cost of reducing the time window slacks [17, 18]. These insertions should not violate the solution feasibility, therefore all of the constraints must be satisfied after inserting each request into the partial solution.

Improvement heuristics modify a low-quality solution by executing refinement procedures in search for better neighboring solutions. The metaheuristics usually embed construction and improvement heuristics [19]. Temporarily infeasible solutions, along with the quality deterioration during the optimization are very often allowed in metaheuristic algorithms [6].

Due to the hierarchical objective of the PDPTW, most of the state-of-the-art heuristic algorithms incorporate a two-stage approach—the number of routes is minimized in the first stage, whereas the travel distance is optimized afterwards [7]. Therefore, efficient techniques for both stages may be designed independently, perhaps without affecting each other. Heuristic algorithms comprise tabu [11] and neighborhood searches [20], agent-based approaches [4], guided ejection searches (GESes), enhanced in our very recent work [8], evolutionary algorithms [21, 22], and many more [7]. Parallel algorithms were intensively explored for solving rich routing problems [23], including the PDPTW [24, 25]. A thorough survey on approximate approaches for the PDPTW was presented by Parragh et al. [6, 26].

Selecting appropriate parameter values of a given algorithm became a challenging task, especially in the case of approaches in which this selection can easily jeopardize the optimization process [5 , 28]. The algorithm parameters can be either retrieved before the execution (in a time-consuming tuning process), or adaptively updated on the fly (using adaptation schemes). However, if certain algorithm components should be applied or not for a given problem instance, then the tuning process appears very difficult, since these components are most likely not independent from each other. Extracting the characteristics of a problem which is to be solved is therefore extremely helpful, if the algorithm behavior is known for certain classes of instances. In this paper, we address this issue and propose a pre-processing procedure which exploits the characteristics of PDPTW instances (importantly, this pre-processing does not have any parameters itself). This allows for an adaptive selection of the best-fitted variant of our enhanced GES [8] for solving the problem at hand as efficiently as possible.

4 Adaptive guided ejection search

4.1 Pre-processing

The pre-processing step aims at determining the appropriate PDPTW problem instance class. In this procedure, we build the feature vectors at first, and then apply the k-nearest neighbor (k-NN) algorithm [29] for classifying the corresponding instance to one class, out of six classes—they are summarized in Table 1 (the corresponding characteristics have been averaged across all problem sizes). The truck capacities are much larger, and the time windows are notably wider for the second-type tests (c2, r2, and rc2) compared with the first-type ones (c1, r1, rc1). Intuitively, a smaller numbers of trucks should be necessary to feasibly serve jobs in the second-type problems (note that the service times are the same for both types of instances of the same structure—see c1 compared with c2, r1 with r2, and rc1 with rc2). Hence, different real-life scenarios are reflected in this dataset. The following subsections discuss in detail the process of extracting features, and classifying PDPTW instances.

4.1.1 Extraction of features

Extracting an appropriate (i.e., discriminative) set of features is a crucial task in all classification problems. Since this step affects the classification performance, it should be executed with care.

In order to exploit the underlying information concerning the structure of the PDPTW instance (being a set of travel points scattered in the two-dimensional space), we incorporate various techniques to determine the number of clusters which fits best to the characteristics of this instance (i.e., the most “appropriate” number of clusters for this dataset). In all of these methods, we utilize the Euclidean metric to find the distance between two travel points (see Section 2). The following aggregation (linkage) techniques are exploited:

Ward—the total within-cluster variance is minimized. At each step, two clusters with the minimum cluster distance are merged. This operation should lead to the minimum increase in the total within-cluster variance. The dissimilarities are squared before the update, as proposed in [30].

McQuitty—the distance D_(i,j) between two clusters C_i and C_j is the weighted mean of the between-cluster dissimilarities [31]: $D_{(i, j)} = \frac{D_{(i, k)} + D_{(i, l)}}{2},$ (11) where the cluster C_j is formed by aggregating the clusters C_k and C_l. The average of the distances between these clusters are considered [32].

Centroid—the distance between two clusters is the squared Euclidean distance between the mean vectors of these clusters: $D_{(i, j)} = {| | \bar{c_{i}} - \bar{c_{j}} | |}^{2},$ (12) where $\bar{c_{i}}$ and $\bar{c_{j}}$ denote the mean vectors of the corresponding clusters. This aggregation appears most robust in the case of isolated points [33].

K-means—this is a reallocation-based aggregation technique. First, the initial (random) cluster centers are selected, and the other vectors are associated with the closest centers to form initial clusters. Then, the centers of the newly created clusters are calculated. The other vectors (which are not the cluster centers) are relocated between clusters, if the new center of a different cluster is closer than the current center of the cluster which the considered vector currently belongs to.

It is worth mentioning that we utilize both hierarchical and non-hierarchical (k-means) clustering.

To find the number of clusters “hidden” in the PDPTW instance, we exploit 27 clustering validity indices (Table 2). The indices combine the information about the inter-cluster isolation along with the information on the intra-cluster compactness in search for the best number of clusters fitted to the underlying structure of the dataset. All of these indices are calculated for all of the aforementioned aggregation techniques. For each aggregation, we calculate the average value of the clusters ( $C$ ) retrieved using these indices: $C = \frac{\sum_{i} C_{i}}{27},$ (13) where $C_{i}$ is the number of clusters suggested by the i-th clustering validity index. Since there are four aggregation methods exploited, four $C$ values (for each technique) are appended to the feature vector.

The analysis of the validity indices is coupled with the clustering by passing messages between data points [9], which was shown to be very efficient (and extremely fast) in clustering of various types of data, especially in computer vision and computational biology. In this technique, all data points are considered as potential cluster centers (they are called exemplars). The messages are then transferred between the data points—the magnitude of each message reflects the current affinity that one data point has for selecting another point as its exemplar (this method is often referred to as the affinity propagation). The messages (being either the responsibility messages r_(i,j), sent from the i-th to the j-th point, and showing how suitable is the j-th point to serve as the exemplar for the i-th point, or the availability messages a_(i,j), sent from the j-th to the i-th point, reflecting how suitable is the j-th point to serve as the exemplar for the i-th point, considering other points for which the j-th point should serve as the exemplar) are exchanged between the exemplars until a high-quality set of clusters is found. The availabilities and responsibilities are finally combined to identify the most appropriate exemplars (clusters). The final number of clusters retrieved using this algorithm is appended to the feature vector.

The clustering-based analysis of the instance is complemented with the histogram-based investigation of the customer locations. We create two-dimensional histograms which show the number of travel locations in various areas of the map (it is divided into n × n cells, and the number of points in each cell is counted). These histograms are extracted at various scales—the number of bins varies, and $2 \leq n \leq ⌈ \sqrt{C} ⌉$ , where C is the number of travel points. Hence—theoretically—each travel point could fall into a separate bin if they were scattered evenly. Finally, we create the histogram counting the bins of each size in these two-dimensional histograms (for each scale). This histogram is quantified using several measures—skewness (γ₁), kurtosis (γ₂), standard deviation (SD), and entropy (h): $γ_{1} = \frac{μ_{3}}{μ_{2}^{3 / 2}},$ (14) $γ_{2} = \frac{μ_{4}}{μ_{2}^{2}} - 3,$ (15) $SD = \sqrt{\frac{1}{(n - 1)} \sum_{i} (b_{i} - \bar{b})^{2}},$ (16) $h = - \sum_{i} (p (b_{i}) {log}_{2} p (b_{i})),$ (17) where n is the number of bins, $\bar{b}$ is the mean histogram value, p (b_i) is the probability mass function of b_i, and μ₂, μ₃, and μ₄ denote the second, third, and fourth moments, respectively. The k-th moment is given as: $μ_{k} = \frac{1}{n} \sum_{i = 1}^{n} {(b_{i} - \bar{b})}^{k} .$ (18)

Skewness is a measure of asymmetry of a histogram—the negative skewness means that its left tail is longer than the right one (it is opposite if skewness is positive). Kurtosis reflects the “peakedness” of a histogram—the positive kurtosis shows that the histogram is more peaked than the Gaussian distribution (it is flatter otherwise). Standard deviation quantifies the variation from the mean value (the higher values, the larger variations). Finally, entropy may be considered as a measure of “uncertainty” in the data distribution [53]. Also, the minimum (b_min), mean ( $\bar{b}$ ), and maximum (b_max) values in the histogram are utilized. All of these seven measures are averaged for all scales, and they are appended to the feature vector.

An exemplary process of extracting these measures (for the instance visualized in Fig. 1) is rendered in Table 3. This analysis is performed for six scales (n = {5, 10, 15, 20, 25, 30}). The two-dimensional histograms are given in Table 3(a), whereas Table 3(b) shows the histograms quantifying the number of bins encompassing a given number of travel points in the corresponding two-dimensional histograms. The extracted measures are presented in Table 4. The averaged values would be appended to the feature vector representing this PDPTW instance. It is worth mentioning that this process may be considered as the analysis of texture of the two-dimensional histograms [53].

As mentioned in Section 2, the transportation requests are served using homogeneous vehicles in the PDPTW. Although the capacity Q of these vehicles is known in advance, it is unclear if this capacity can be considered “large” or “small”. To verify if the minimum number of trucks which can feasibly serve the requests correspond to the number of clusters found in this instance, we calculate the following ratio: $C_{Q} = \frac{D}{Q \cdot \bar{C}},$ (19) where D is the total demand, and $\bar{C}$ denotes the average number of clusters retrieved using all of the aforementioned five clustering techniques. If $C_{Q} > 1$ , then the number of clusters is larger than the minimum number of feasible routes, and the truck capacities can be considered large (i.e., it may be possible to serve more than one cluster of travel points using a single vehicle). This ratio is appended to the feature vector.

Therefore, each feature vector contains 13 entries:

4 average numbers of clusters (the analysis of the clustering validity indices),

1 number of clusters extracted using the affinity propagation technique,

7 characteristics averaged for various scales in the two-dimensional histogram based analysis,

1 ratio quantifying how “large” is the maximum truck capacity.

4.1.2 Classification

In this work, we utilize the k-NN algorithm for classifying the incoming instance to be solved to one (out of six) class (see Table 1). The nearest-neighbor decision rule assigns the appropriate class to the new (unseen) feature vector based on the classification outcome of k already-classified vectors nearest (in the Euclidean distance) to the incoming one. If k is odd, then the classification is decided by the majority vote. On the other hand, if k is even and the voting process leads to the tie, then the class may be either assigned randomly (for the tied classes), or selected based on quantifying the “importance” of separate votes.

Although the k-NN algorithm is a relatively basic classification engine, it was shown to be effective in many pattern recognition tasks [53]. Also, it is easily applicable to multi-class classification problems.

Algorithm 1 Enhanced guided ejection search.

1: σ_B ← initial solution; finished ← false;

2: while notfinisheddo

3: Initialize EP with requests from random r;

4: p[i] ←1; (i = 1, 2, …, N); σ_I ← σ;

5: while (EP ≠ ∅) and (notstuck) do

6: Select and remove request h_in from EP;

7: if $S_{in}^{fe} (h_{in}, σ) \neq \emptyset$ then

8: σ ← $rand (σ^{'}) \in S_{in}^{fe} (h_{in}, σ)$ ;

9: elseσ ← Squeeze(h_in, σ);

10: end if

11: ifh_in ∉ σthen

12: p [h_in] ← p [h_in] +1;

13: fork ← 1 tok_mdo

14: Generate $S_{ej}^{fe} (h_{in}, σ)$ ;

15: if $S_{ej}^{fe} (h_{in}, σ) \notin \emptyset$ then

16: σ← $rand (σ^{'}) \in S_{ej}^{fe} (h_{in}, σ)$ ;

17: Add { $h_{out}^{1}, \dots, h_{out}^{k}$ } to EP;

18: break;

19: end if

20: end for

21: end if

22: Perturb σ; c ← c + 1;

23: Verify termination conditions;

24: Update finished and stuck;

25: end while

26: if EP ≠ ∅ σ ← σ_I;

27: else σ_B ← σ;

28: end if

29: end while

30: returnσ_B;

4.2 Enhanced guided ejection search

Our enhanced guided ejection search is an improvement heuristics which boosts the quality of the initial (feasible—all constraints discussed in Section 2.2 are satisfied) PDPTW schedule, in which every request is served in a separate route (Algorithm 1). At each step, a random route r is drawn, its requests are inserted into the ejection pool (EP) (line 3), and the attempts to re-insert them back into the partial solution are undertaken to minimize K (see Section 2.1). The penalty counters p indicate the re-insertion difficulty of the corresponding request (the higher p value, the more difficult is to include the request back in the solution). A request h_in is popped from the EP (line 6). If there are several feasible positions to insert h_in into σ, then a random one is chosen (line 8)— $S_{in}^{fe} (h_{in}, σ)$ is the set of feasible insertions of h_in. This request is otherwise inserted so as it violates the constraints (see Section 2.2), and the solution is squeezed to restore its feasibility (line 9)—this is the local-search based procedure [8]. If it fails, then its p is increased (line 12), and at most k_m requests (with minimal sum of their p’s) are ejected to insert h_in feasibly (lines 13–20). Finally, σ is perturbed with local moves (out-relocate and out-exchange) to diversify the search (line 22). This process continues until the maximum time has elapsed, or the minimum theoretical number of routes (K_min) has been retrieved. For more details on our GES,see [8].

In our recent paper [8], we proposed three significant improvements to the GES which clearly affected the quality of the retrieved PDPTW routing schedules:

Squeezing of infeasible solutions (S),

Randomizing the EP (R),

Imposing the maximum iteration limit on the analysis of “difficult” requests (M).

All of these enhancements (highlighted in Algorithm 1) notably affected the GES capabilities. We showed that they should be applied for different classes of instances (the guidelines are presented in Table 5—the letters in the acronym correspond to the enhancements, which are switched on or off). The selection of an appropriate GES variant (adapting its important components) before the execution was unclear—this issue is tackled by the adaptation procedure.

4.3 Adaptation

In the adaptation process, which is performed before the GES run, we select the appropriate GES variant based on the PDPTW instance class. First, we classify the incoming instance into one (out of six) classes (see Section 4.1). Then, Table 5 is used to retrieve the best GES variant to tackle the problem at hand. This variant does not change during the execution.

5 Experimental validation

The pre-processing procedure was implemented in R, and it was run on an Intel Xeon 3.2 GHz computer (16 GB RAM). In this study, we focus on the analysis of the efficacy of this algorithm step, since based on its outcome, the appropriate GES variant is selected (the results of our previous efforts showed how the GES enhancements affect its capabilities) [8].

To verify the classification performance, we conducted 10-fold cross-validation on the full Li and Lim’s benchmark set of 354 problem instances belonging to 6 classes (see Table 1). In Fig. 2, we visualize some exemplary PDPTW tests from this benchmark—it is worth noting that assigning a proper class is certainly not obvious (see e.g., the clustered test for 400 customers). In each class, there are problems with 100, 200, 400, 600, 800, and 1000 customers (the numbers of different-size problems are roughly the same) 1 .

Extraction of histogram-based features is extremely fast, and took less than 1 second for all PDPTW instances (of all sizes). Similarly, the affinity propagation time is negligible compared with the time of the clustering validity indices analysis. The latter feature extraction approach was the most time-consuming, and required less than 3 seconds for 100-customer tests, up to approximately 300 seconds for 1000-customer instances containing clustered travel points. However, the pre-processing is executed once for each PDPTW instance (before the GES execution), hence this initial analysis time is acceptable even for large instances.

Table 6 reports the classification accuracy (in %) along with the standard error of the k-NN algorithm obtained for various k’s. To verify how the extracted features affect the scores, we investigated four types of feature vectors, encompassing: (a) 5 clustering-based features (the Clustering variant), (b) 7 histogram-based features (Histogram), (c) both clustering and histogram-based ones (Clustering+histogram), and (d) all of 13 features discussed in Section 4.1 (All). The results revealed that the selection of k has a significant impact on the classification scores for all variants of the feature vectors (the best performance was retrieved using k ≤ 5). Similarly, the standard errors are notably smaller for lower k’s in most cases (e.g., this error is more than 8.3× larger for k = 15 compared with k = 3 for All). This indicates that the classifiers with smaller k’s generalize well to unseen data.

Utilizing the feature vectors containing all extracted features allows for obtaining the best classification performance. Although both Clustering and Histogram variants gave very high classification scores (88.98% and 92.4% for k = 2, respectively), combining these features resulted in the classification accuracy close to 99% (see Clustering+histogram). Finally, the perfect classification (the accuracy equal to 100%) was elaborated using All, in which all features are exploited (see All for k = 2). Importantly, there were no classification errors in this case—the set of the extracted features is quite discriminative and can be effectively used in this classification task (we deal with the six-class classification problem, therefore it is not trivial).

In order to investigate the statistical significance of the classification results, we conducted the two-tailed Wilcoxon tests. We verify the null hypothesis saying that “using different feature vectors (e.g., including certain subsets of all features) leads to obtaining the classification of the same quality”. The levels of the statistical significance are gathered in Table 7. The results confirm that the differences are important for all of the investigated variants (the p-values are smaller than 0.0001, therefore the null hypothesis can be safely rejected). Hence, all of the extracted features appear important and should be utilized to classify the incoming PDPTW instances, since exploiting all features resulted in the highest classification quality (see Table 6).

6 Conclusions and future work

In this paper, we proposed an adaptive guided ejection search for the PDPTW. The pre-processing step, which is the pivotal part of our algorithm, classifies the PDPTW instance to be solved to one (out of six) class, reflecting its most important characteristics. We extract 13 discriminative features (including clustering-based features, and the histogram-based ones—extracted at various scales). In our previous work [8], we showed that selecting an appropriate variant of the GES is crucial and affects its convergence capabilities. It is addressed in the adaptation procedure, which—based on the classification outcome—retrieves the best GES variant before its execution. The extensive experimental study revealed that the classification performance is very high, and applying our pre-processing procedure allows for retrieving the most appropriate GES variant efficiently. The two-tailed Wilcoxon tests proved that the differences in the classification accuracy obtained using various feature vectors (with and without certain features appended) are statistically important, thus utilizing all of the extracted features is beneficial.

Our ongoing research is focused on coupling our GES with the evolutionary algorithm to minimize the travel distance. This will result in creating a full optimization framework for the PDPTW (and potentially other challenging VRPs [54])—we plan to compare our adaptive GES with other state-of-the-art algorithms. We will investigate the performance of our pre-processing on massive problems, and those elaborated using our benchmark generator [55]. Finally, we work on the adaptive algorithm which will handle dynamic scenarios (e.g., changes in the transportation network due to the traffic congestion) on the fly.

Footnotes

For details, see: . Last access date: May 4, 2016.

Acknowledgments

This research was supported by the National Science Centre under research Grant No. DEC-2013/09/N/ST6/03461, and performed using the Intel CPU and Xeon Phi platforms provided by the MICLAB project No. POIG.02.03.00.24-093/13.

References

Pillac

, Gendreau

, Guéret

and Medaglia

A.L.

, A review of dynamic vehicle routing problems, European Journal of Operational Research225(1) (2013), 1–11.

Gulczynski

, Golden

and Wasil

, The period vehicle routing problem: New heuristics and real-world variants, Transportation Research Part E: Logistics and Transportation Review47(5) (2011), 648–668.

Lau

H.C.

and Liang

, Pickup and delivery with time windows: Algorithms and test case generation, International Journal on Artificial Intelligence Tools11(03) (2002), 455–472.

Kalina

and Vokrinek

, Parallel solver for vehicle routing and pickup and delivery problems with time windows based on agent negotiation, in Proc IEEE SMC, 2012, 1558–1563.

Nalepa

and Blocho

, Adaptive memetic algorithm for minimizing distance in the vehicle routing problem with time windows, Soft Comp20(6) (2016), 2309–2327.

Parragh

S.N.

, Doerner

K.F.

and Hartl

R.F.

, A survey on pickup and delivery problems, Journal fur Betriebswirtschaft58(1) (2008), 21–51.

Nagata

and Kobayashi

, Guided ejection search for the pickup and delivery problem with time windows, in Proc EvoCOP, vol. 6022 of LNCS, Springer, 2010, pp. 202–213.

Nalepa

and Blocho

, Intelligent Information and Database Systems: Proc 8th Asian Conference, ACIIDS 2016, ch. Enhanced Guided Ejection Search for the Pickup and Delivery Problem with Time Windows, Springer, 2016, pp. 388–398.

Frey

B.J.

and Dueck

, Clustering by passing messages between data points, Science315(5814) (2007), 972–976.

10.

Grandinetti

, Guerriero

, Pezzella

and Pisacane

, The multi-objective multi-vehicle pickup and delivery problem with time windows, Social and Beh Sc111 (2014), 203–212.

11.

Nanry

W.P.

and Barnes

J.W.

, Solving the pickup and delivery problem with time windows using reactive tabu search, Transportation Research34(2) (2000), 107–121.

12.

Cordeau

J.-F.

, A branch-and-cut algorithm for the dial-a-ride problem, Op Res54(3) (2006), 573–586.

13.

Bettinelli

, Ceselli

and Righini

, A branch-and-price algorithm for the multi-depot heterogeneous-fleet pickup and delivery problem with soft time windows, Mathematical Programming Computation6(2) (2014), 171–197.

14.

Baldacci

, Bartolini

and Mingozzi

, An exact algorithm for the pickup and delivery problem with time windows, Opererations Research59(2) (2011), 414–426.

15.

Cordeau

J.-F.

, Laporte

and Ropke

, The Vehicle Routing Problem: Latest Advances and New Challenges, ch. Recent Models and Algorithms for One-to-One Pickup and Delivery Problems, Boston, MA: Springer, 2008, pp. 327–357.

16.

Baldacci

, Mingozzi

and Roberti

, Recent exact algorithms for solving the vehicle routing problem under capacity and time window constraints, European Journal of Operational Research218(1) (2012), 1–6.

17.

and Dessouky

M.M.

, A new insertion-based construction heuristic for solving the pickup and delivery problem with time windows, European Journal of Operational Research175(2) (2006), 672–687.

18.

Zhou

, Tan

, Liao

and Liu

, Solving the multi-vehicle pick-up and delivery with time widows by new construction heuristic, in Proc CISDA, vol. 2, 2006, pp. 1035–1042.

19.

and Lim

, A metaheuristic for the pickup and delivery with time windows, in Proc ICTAI, 2001, pp. 160–167.

20.

Ropke

and Pisinger

, An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows, Transp Science40(4) (2006), 455–472.

21.

Pankratz

, A grouping genetic algorithm for the pickup and delivery problem with time windows, OR Spectrum27(1) (2005), 21–41.

22.

Nagata

and Kobayashi

, Proc PPSN XI, ch. A Memetic Algorithm for the Pickup and Delivery Problem with Time Windows Using Selective Route Exchange Crossover, Springer, Heidelberg, 2010, pp. 536–545.

23.

Nalepa

and Blocho

, Co-operation in the parallel memetic algorithm, International Journal of Parallel Programming43(5) (2015), 812–839.

24.

Blocho

and Nalepa

, A parallel algorithm for minimizing the fleet size in the pickup and delivery problem with time windows, in Proc of 22nd European MPI Users’ Group MeetingEuroMPI ’15, New York, USA, ACM2015, pp. 15:1–15:2.

25.

Nalepa

and Blocho

, A parallel algorithm with the search space partition for the pickup and delivery with time windows, in Proc 3PGCIC, 2015, pp. 92–99.

26.

Cherkesly

, Desaulniers

and Laporte

, A populationbased metaheuristic for the pickup and delivery problem with time windows and LIFO loading, Computers & Operations Research62 (2015), 23–35.

27.

Nalepa

and Kawulok

, A memetic algorithm to select training data for support vector machines, in Proc of the 2014 Annual Conference on Genetic and Evolutionary Computation, GECCO ’14, New York, USA, ACM2014, pp. 573–580.

28.

Nalepa

and Kawulok

, Adaptive memetic algorithm enhanced with data geometry analysis to select training data for SVMs, Neurocomputing185 (2016), 113–132.

29.

Cover

and Hart

, Nearest neighbor pattern classification, IEEE Trans on Inf Theory13(1) (1967), 21–27.

30.

Murtagh

and Legendre

, Ward’s hierarchical agglomerative clustering method, Journal of Classification31(3) (2014), 274–295.

31.

McQuitty

L.L.

, Similarity analysis by reciprocal pairs for discrete and continuous data, Educational and Psychological Measurement26(4) (1966), 825–831.

32.

Clifford

, Wessely

, Pendurthi

and Emes

R.D.

, Comon of clustering methods for investigation of genomewide methylation data, Front in Gen2(88) (2011)–paris.

33.

Charrad

, Ghazzali

, Boiteau

and Niknafs

, NbClust: An R package for determining the relevant number of clusters in a data set, Journal of Statistical Software61(1) (2014), 1–36.

34.

Calinski

and Harabasz

, A dendrite method for cluster analysis, Comm in Statistics3(1) (1974), 1–27.

35.

Duda

R.O.

and Hart

P.E.

, Pattern Classification and Scene Analysis. New York: John Willey & Sons, 1973.

36.

Hubert

L.J.

and Levin

J.R.

, A general statistical framework for assessing categorical clustering in free recall, ERIC Clearinghouse (1975).

37.

Frank

L.J.H.

and Baker

, Measuring the power of hierarchical cluster analysis, Journal of the American Statistical Association70(349) (1975), 31–38.

38.

Milligan

G.W.

, A monte carlo study of thirty internal criterion measures for cluster analysis, Psychometrika46(2) (1981), 187–199.

39.

Davies

D.L.

and Bouldin

D.W.

, A cluster separation measure, IEEE Transactions on Pattern Analysis and Machine IntelligencePAMI-1(2) (1979), 224–227.

40.

Hartigan

J.A.

, Clustering Algorithms, New York, NY, USA: John Wiley & Sons, Inc., 99th ed., 1975.

41.

Ratkowsky

D.A.

and Lance

G.N.

, A criterion for determining the number of groups in a classification, Australian Computer Journal10 (1978), 115–117.

42.

Scott

M.J.S.A.J.

, Clustering methods based on likelihood ratio criteria, Biometrics27(2) (1971), 387–397.

43.

Marriott

F.H.C.

, Practical problems in a method of cluster analysis, Biometrics27(3) (1971), 501–514.

44.

Ball

and Hall

, Isodata: A novel method of data analysis and pattern classification, tech. rep., Stanford Research Institute, Menlo Park, 1965.

45.

Milligan

G.W.

and Cooper

M.C.

, An examination of procedures for determining the number of clusters in a data set, Psychometrika50(2) (1985), 159–179.

46.

Friedman

J.R.H.P.

, On some invariant criteria for grouping data, Journal of the American Statistical Association62(320) (1967), 1159–1178.

47.

John

V.R.R.

and McClain

, CLUSTISZ: A program to test for the quality of clustering of a set of objects, Journal of Marketing Research12(4) (1975), 456–460.

48.

Krzanowski

Y.T.L.W.J.

, A criterion for determining the number of groups in a data set using sum-of-squares clustering, Biometrics44(1) (1988), 23–34.

49.

Rousseeuw

P.J.

, Silhouettes: A graphical aid to the interpretation and validation of cluster analysis, Journal of Computational and Applied Mathematics20 (1987), 53–65.

50.

Dunn

J.C.

, Well-separated clusters and optimal fuzzy partitions, Journal of Cybernetics4(1) (1974), 95–104.

51.

Halkidi

, Vazirgiannis

and Batistakis

, Proc PKDD, ch. Quality Scheme Assessment in the Clustering Process, Springer, Heidelberg, 2000, pp. 265–276.

52.

Halkidi

and Vazirgiannis

, Clustering validity assessment: Finding the optimal partitioning of a data set, in Proc IEEE ICDM, 2001, pp. 187–194.

53.

Nalepa

, Szymanek

, Hayball

M.P.

, Brown

S.J.

, Ganeshan

and Miles

, Computer Vision and Graphics: Proc ICCVG 2014, ch. Texture Analysis for Identifying Heterogeneity in Medical Images, Springer, 2014, pp. 446–453.

54.

Trinh

, Dang

and Dinh

, New Challenges for Intelligent Information and Database Systems, ch. An Approximation Approach for a Real–World Variant of Vehicle Routing Problem, Springer, Heidelberg, 2011, pp. 87–96.

55.

Cwiek

, Nalepa

and Dublanski

, Intelligent Information and Database Systems: Proc 8th Asian Conference, ACIIDS 2016, ch. How to Generate Benchmarks for Rich Routing Problems?Springer, Heidelberg, 2016, pp. 399–409.