An iterative three-component heuristic for the team orienteering problem with time windows

Accepted Manuscript Discrete Optimization An Iterative Three-Component Heuristic for the Team Orienteering Problem with Time Windows Qian Hu, Andrew Lim PII: DOI: Reference:

S0377-2217(13)00480-3 http://dx.doi.org/10.1016/j.ejor.2013.06.011 EOR 11728

To appear in:

European Journal of Operational Research

Received Date: Accepted Date:

15 May 2012 8 June 2013

Please cite this article as: Hu, Q., Lim, A., An Iterative Three-Component Heuristic for the Team Orienteering Problem with Time Windows, European Journal of Operational Research (2013), doi: http://dx.doi.org/10.1016/ j.ejor.2013.06.011

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

An Iterative Three-Component Heuristic for the Team Orienteering Problem with Time Windows Qian Hua , Andrew Lima,∗ a

Department of Management Sciences, City University of Hong Kong, Tat Chee Ave, Kowloon Tong, Hong Kong

Abstract This paper studies the team orienteering problem with time windows, the aim of which is to maximize the total profit collected by visiting a set of customers with a limited number of vehicles. Each customer has a profit, a service time and a time window. A service provided to any customer must begin in his or her time window. We propose an iterative framework incorporating three components to solve this problem. The first two components are a local search procedure and a simulated annealing procedure. They explore the solution space and discover a set of routes. The third component recombines the routes to identify high quality solutions. Our Computational results indicate that this heuristic outperforms the existing approaches in the literature in average performance by at least 0.41%. In addition, 35 new best solutions are found. Keywords: Routing, Team orienteering problem with time windows, Heuristic

1. Introduction The aim of the team orienteering problem with time windows (TOPTW) is to maximize the total scores collected by visiting a set of locations, each of which has a score, a service time and a time window. The number of routes is limited, and each location can be visited at most once. A route is feasible if no time window for a visited location is violated, and if it begins and ends at the same given location. The TOPTW has numerous real-life applications (Tsiligirides, 1984; Golden et al., 1987; Souffriau et al., 2008). In this paper, we consider the vehicle routing application in ∗

Corresponding author Email addresses: [email protected] (Qian Hu), [email protected] (Andrew Lim )

Preprint submitted to European Journal of Operational Research

June 12, 2013

which a shipper sends out a fixed number of vehicles from a depot to serve some of its customers with the aim of maximizing the total profit gained from successful service. The main contribution of this paper is an iterative three-component heuristic (I3CH) that is straightforward but every effective. The heuristic employs a local search procedure and a simulated annealing procedure as its first two components to explore the solution space. They store the discovered routes into a pool for future use. The third component then recombines these routes by solving a set packing formulation to obtain a high quality feasible solution. Embedded in an iterative structure, these three components cooperate and perform effectively. We conduct computational experiments on the TOPTW test instances in the literature, and find our heuristic to outperform the existing approaches in terms of average performance. Before these experiments, we classify the test instances with known optimal solutions into an “OPT” category and the remainder into an “INST-M” category. The average performance improvement for the “INST-M” instances is at least 0.49%, and 35 new best solutions are found. For the 66 instances in the “OPT category”, our heuristic found 55 optimal solutions, 16 more than the previous best approach in the literature. In addition, our approach runs efficiently. The remainder of this paper is organized as follows. Section 2, provides a clear definition of the TOPTW and a brief overview of the research literature addressing it. Section 3 describes the three-component heuristic and its components, and Section 4 extends the heuristic by constructing an iterative framework. Section 5 evaluates the effectiveness of our approach and compares it to the known approaches in the literature. Finally, Section 6 concludes that paper with closing remarks and possible directions for future research. The detailed solution values for each test instance are relegated to the appendix. 2. Problem definition and literature review The TOPTW examined in this paper is defined as follows. On a given network graph G = (V, A), there are n + 1 different locations denoted by the set V = {0, 1, 2, . . . , n} and a set of arcs connecting these locations which is denoted as A = {(i, j) : i
j ∈ V}. The travel time, ti j , from location i to j is equal to the Euclidean distance, di j , from i to j. Let location 0 be the depot and each of the remaining location corresponds to one customer. The number of available vehicles is 2

fixed at m. Each vehicle must begin and end its route at the depot within the depot’s time window [O0 , C0 ]. Each customer i = 1, . . . , n with a profit pi , a service time T i and a time window [Oi , Ci ] can be visited at most once. The service delivered to the customer is successful if it begins within his or her time window. In the case of an earlier arrival, the vehicle has to wait until the time window begins. Profit is gained upon successful service. Note that owing to the limited number of vehicles, some customers may not be visited in a feasible solution. The objective of the TOPTW is to identify a feasible routing plan that ensures the maximization of total profit. The TOPTW is one of many variants of the widely studied Orienteering Problem (OP) (Schilde et al., 2009; Chao et al., 1996a; Fischetti et al., 1998). Vansteenwegen et al. (2011) published an excellent review of the OP and, its applications and variants, including the team orienteering problem (TOP) (Chao et al., 1996b; Tang and Miller-Hooks, 2005; Vansteenwegen et al., 009a), the OP with time windows (OPTW) (Kantor and Rosenwein, 1992; Righini and Salani, 2009) and the TOPTW. Given a set of locations with a score, OP seeks to determine one path that visits some vertices and maximizes the total score. The TOP extends the OP to identify multiple paths that maximize the total score. The OPTW and TOPTW extend the OP and TOP, respectively, by incorporating time windows constraints. All these problems have many applications (Tsiligirides, 1984; Golden et al., 1987; Souffriau et al., 2008), and are challenging because of their complexity. As the OP is NP-hard (Golden et al., 1987), the TOPTW must also be NP-hard. Therefore, solving the highly constrained TOPTW optimally in polynomial time is highly unlikely. Recent research shows heuristics and meta-heuristics to be the most popular techniques to solve this problem. The TOPTW has only recently become popular, and several state-of-the-art papers on the problem can be found in the literature. For example, Montemanni and Gambardella, in 2009, published a paper entitled “The team orienteering problem with time windows” (Montemanni and Gambardella, 2009). They designed an ant colony system (ACO) to solve the problem and achieved the best results on OPTW instances. In addition, they contributed the earliest test instances for this benchmark research. More recently, Montemanni et al. (2011) and Gambardella et al. (2012), extended the ACO and reported new best solutions. Vansteenwegen et al. (009b) subsequently proposed an iterated local search (ILS) algorithm for the TOPTW, an approach that runs significantly faster than the ACO requiring only a few seconds of computational time to solve 3

an instance while maintaining competitive solution quality. In addition, the authors constructed a new set of test instances with known optimal solutions. Tricoire et al. (2010) solved an even more complicated problem, the multi-period orienteering problem with multiple time windows, which is a generalization of the TOPTW. They obtained good-quality solutions when they tested their variable neighbor search (VNS) meta-heuristic on the TOPTW instances. Souffriau et al. (2011) studied another generalization, the multi-constraint team orienteering problem with multiple time windows, and found its solution method also to perform well on TOPTW test instances. Concentrating on the TOPTW, Labadi et al. (2011) developed a hybrid meta-heuristic that combines the greedy randomized adaptive search procedure (GRASP) with the evolutionary local search (ELS) approach. This method improved several of the best known results on available benchmark instances at that time. In 2012, Lin and Yu developed two versions of simulated annealing algorithm for the TOPTW (Lin and Yu, 2012). The fast version computes a solution within only several seconds, while the slow version achieves better solutions at the expense of more computational time. Some current best-known solutions are obtained by the slow version (SSA). Most recently, Labadie et al. (2012) provided new results for the TOPTW. They developed a VNS algorithm that explores granular neighborhoods (GVNS) based on linear programming. 3. Three-component heuristic for the TOPTW As we have seen, most of the existing approaches are neighborhood search approaches. When a neighborhood search approach solves a TOPTW instance, not only are new solutions are explored, but a set of routes is discovered. A feasible solution can be defined as a combination comprising m routes that satisfy route feasibility and the condition that each customer can be served at most once. Given a set of discovered routes, there exists a best combination that is always no worse than the explored solutions. We derive the two following straightforward “no-worse-than” propositions. The proofs are omitted. Proposition 1. Suppose that a neighborhood search algorithm A solves the TOPTW. It searches N neighbors and obtains a solution S 1 . At the same time, it also discovers a set of mN routes. The best combination S 2 of these routes is no worse than S 1 . 4

Proposition 2. Suppose that two neighborhood search approaches, A1 and A2 , both solve the TOPTW. A1 discovers a set of routes R1 and obtains a solution S 1 . A2 discovers a set of routes R2 and obtains a solution S 2 . The best combination, S 3 , over R1 ∪ R2 is no worse than S 1 or S 2 . On the basis of these two propositions, we propose a heuristic that first employs two different neighborhood search approaches to search the solution space, then keeps the discovered routes into a pool, and finally recombines them to obtain the best combination. We call it a three-component heuristic. In our design, the first two components are local search and simulated annealing. They explore the solution space and store routes. The third component is route recombination which recombines the routes to produce the best combination. Figure 1 illustrates the heuristic’s structure.

Route pool Local search

Simulated annealing

Route recombination

Best solution

Figure 1: The structure of the three-component heuristic.

In the following subsections, we first introduce the representation of a solution, and then the neighborhood operator that will be used in the local search and simulated annealing. Finally, we demonstrate the three components. 3.1. Solution representation In the solution encoding, we use m lists to represent the m routes. Each list starts and ends with 0, which is the depot. Denote the routes by ri , i = 1, . . . , m. Customers who are not visited by any 5

route are stored in another list u. For example, given V = {0, 1, 2, . . . , 10} and m = 2, a potential solution is represented by r1 = {0, 1, 3, 5, 0}, r2 = {0, 2, 4, 6, 8, 0} and u = {7, 9, 10} (see Figure 2).

1 10 3 9

Solution representation

0

r1 = { 0, 1, 3, 5, 0 } r2 = { 0, 2, 4, 6, 8, 0 }

2

5

u = { 7, 9 10 } 8

4

Route

7

Depot Customer

6

Figure 2: Example of a solution representation.

3.2. The neighborhood operator: Eliminator Both the local search and simulated annealing procedures are neighborhood search approaches. We propose a neighborhood operator called eliminator to produce neighborhood solutions. The eliminator first removes some customers from some routes, and then inserts the unvisited customers from u. If the newly added customers are more profitable, solution quality is improved. Such newly generated solutions are then further improved through a post-processing procedure. We say a solution is complete if no more customers from u can be inserted into any route, otherwise the solution is partial. A partial solution can be improved by inserting some customers from u while keeping the routes feasible. Given a solution A, the eliminator randomly removes some customers from some routes to obtain a partial solution A . Suppose that ri , i = 1, . . . , m are routes of A and that u is the list of unvisited customers. We can then improve A to obtain a complete solution B by iteratively removing a customer c0 from the head of u and attempting to insert c0 into any route ri . The insertion should be made at the first feasible insert position. If no feasible insert position exists, c0 6

is moved to the end of u. The procedure is terminated when no more customers can be inserted. The resulting solution B is one of the many neighbors of A. Since the eliminator is stochastic, a set of solution neighbors can be produced by applying it to A multiple times. To enhance the eliminator, we set rules on how to remove customers. Let p¯ be the average profit over all customers on m routes. Customer c j is eliminated with a probability Ph if p j ≥ p, ¯ and a probability Pl if p j < p. ¯ Considering Ph < Pl , it indicates that eliminator prefers to retain customers that represent a higher profit. Therefore, with employment of the eliminator, the local search and simulated annealing are to find better neighbor solutions. We next apply the post-processing procedure to improve the newly generated solution B. This post-processing procedure iteratively invokes seven operators until further improvement is impossible. The seven operators are divided into three types: relocate, exchange and 2-opt. Relocate operators. Relocate operators insert a customer c j at a feasible insertion position on a route ri . As c j may come from different routes or from u, relocate operators have three variants: 0-relocate, 1-relocate and 2-relocate. 0-relocate functions when c j comes from u. For a successful operation, 0-relocate improves the solution by additional profit p j . If the relocation takes place on the same route, it is a 1-relocate operation. Finally, a 2-relocate operation is to relocate a customer from one route to another. Neither 1-relocate nor 2-relocate can increase the total profit, but they can reduce the travel distance. We believe that a solution with shorter travel distance has a more desirable structure, because it may provide greater opportunity for the other operators to improve the solution quality. Exchange operators. Exchange operators swap two customers in the solution. Similar to relocate operators, they also have three variants: 0-exchange, 1-exchange and 2-exchange. Of the three, 0-exchange is the only exchange operator to improve the objective value. The solution quality is improved when 0-exchange replaces an existing customer on one route with a more profitable customer from u. 1-exchange is an intra-route operator while 2-exchange is an interroute operator. Both of them are used to improve the solution structure by reducing the travel distance. 2-opt operator. 2-opt operator behaves in the same way that it does when used to solve the vehicle routing problem with time windows (VRPTW). It first selects two edges from two different 7

routes. These two edges separate the two routes into four parts. After a crossover, the 2-opt might recombine them into two new feasible routes reducing the total travel distance.

0-relcoate

0-exchange

1-relocate

1-exchange

2-relocate

2-opt

2-exchange

Figure 3: Illustration of operators in the post-processing procedure

The seven operators are illustrated in Figure 3. Only 0-relocate and 0-exchange contribute to improving solution quality. The remaining operators restructure a TOPTW solution. The restructured solutions can be desirable as it may provide opportunities for 0-relocate and 0-exchange operators to make improvements. Algorithm 1 illustrates the post-processing procedure. 8

Algorithm 1 Post-processing procedure (PP) Require: A solution S ; 1: Set impr ← 1; 2: while impr > 0 do 3:

Invoke the 2-relocate, 2-opt and 2-exchange operators on S and obtain S  ;

4:

Apply the 1-relocate and 1-exchange operators on S  and obtain S  ;

5:

Run the 0-relocate and 0-exchange operators on S  and obtain S  ;

6:

Set impr ← 1 if there is an improvement, otherwise set impr ← 0;

7:

S ← S  ;

8: end whilethe solution of S ;

3.3. Local Search and Simulated Annealing The first two components of our heuristics are Local Search and Simulated Annealing. As previously noted, their objective is to search for a better solution in the neighborhood and then store the discovered routes during the exploration. The Local Search procedure (LS) iteratively explores the neighborhood of a starting solution X. In each iteration, LS generates N neighbors using the neighborhood operator eliminator. The solution of each neighbor contains a list of routes. LS caches these routes into a pool POOL. At the end of each iteration, LS updates X with the best neighbor if there is one. Otherwise, X is replaced by the last explored neighbor. Let XLS be the best solution found. If X is better than XLS , there is an improvement and we update XLS with X. Let Ils no impr be the number of current consecutive iterations with no improvement. Once there is an improvement, set Ils no impr = 0. Otherwise, increase Ils no impr by 1. LS searches N neighbors and then decides which one to move to, whereas the Simulated Annealing procedure (SA) behaves in a different way: each time it produces only one neighbor and accepts it with a computed probability. Simulated annealing has been widely studied in the literature and employed in many applications. It is an attractive technique because of its ease of implementation, convergence properties and ability to escape from local optima. As one of the most popular approaches to solve the routing problems, it has been used to solve the traveling 9

salesman problem (Malek et al., 1989), the vehicle routing problem (Van Breedam, 1995), the vehicle routing problem with time windows (Chiang and Russell, 1996) and other routing problems (Lin et al., 2009; Yu et al., 2010; Lin and Yu, 2012). For a more in depth discussion of simulated annealing, we recommend Eglese (1990), Koulamas et al. (1994), Suman and Kumar (2005) and Gendreau and Potvin (2010). We define three parameters T 0 , α, and Isa no impr in SA to denote the initial temperature, cooling speed and number of current consecutive iterations without improvement respectively. Let Y be the starting solution and YSA the best solution found by SA. For each step, SA invokes the eliminator to obtain one neighbor of Y, denoted by Y  . If Y  is better than Y, SA will move to Y  directly. Otherwise, SA will accept this neighbor with a probability P sa , which is computed with Equation (1). At the end of each step, we update T = α ∗ T . To be consistent with LS on the term of “iteration”, we redefine an iteration in SA as the whole of N consecutive steps: that is, SA runs N steps in an iteration, where N is exactly the same as in LS. We can then update Isa no impr in a similar way. Once there is an improvement in an iteration, we set Isa no impr = 0. Otherwise, we increase Isa no impr by 1. SA reports YSA at the end. It also caches all routes of the visited neighbors into POOL. P sa = exp{

1 Y − Y }, T YSA

(1)

3.4. Route recombination The heuristic’s third component, Route recombination (RR), receives the routes from POOL and recombines them to obtain high-quality solutions. After taking over POOL, RR solves a set packing formulation over the routes to produce the best combination, which is a feasible solution of the TOPTW. Providing a set of routes in POOL = {r1 , r2 , . . . , rS pool }, where S pool is the size of POOL. For j ∈ V\{0} and k ∈ {1, 2 . . . , S pool }, let a jk = 1 if customer j is visited by route k, and a jk = 0 otherwise. The total profit of a route rk is
computed as qk = j∈V\{0} p j a jk . We define xk = 1 if route rk is selected in the final solution, and xk = 0 otherwise. Therefore, we derive the set packing formulation as follows.

10

RR :

S pool

Maximize k=1 S pool

Subject to k=1 S pool

qk xk a jk xk ≤ 1,

(2) ∀ j ∈ V\{0}

xk ≤ m

(3)

(4)

k=1

xk ∈ {0, 1},

∀ k ∈ {1, 2 . . . , S pool }

(5)

The objective (2) maximizes the total collected profit. Constraints (3) ensure that each customer can only be visited at most once. Constraint (4) guarantees that no more than m routes be combined into a new solution. This inequality equipped with “≤” is reasonable as m can be so large for some test instances that all customers may already have been visited by fewer than m vehicles. The aim of the operation is to pack more profitable routes into the generating solution to maximize the total profit. If POOL is completed by including all feasible routes, then RR will compute an optimal solution for the TOPTW. However, the complete set has exponential elements which make it difficult to solve RR in practice. Our design allows us to set an upper limit on S pool . Following this capacity setting, RR loads only partial routes from and obtains an Integer Programming (IP) solution that may not be optimal for the TOPTW. Such a strategy constitutes a trade-off between solution quality and computational efforts. With a suitable S pool , commercial MIP solvers, such as ILOG CPLEX 12.2, are capable of solving RR within a reasonable amount of computing time. To further accelerate the RR component, we can make use of the upper bound obtained from its LP relaxation. Because POOL includes only a limited number of routes, there is a possibility that the solution to be obtained by RR is no better than the best solution we already achieved. If that is the case, then the computational efforts of the MIP solver would be in vain. Thus, we first solve the LP relaxation of RR to obtain an upper bound. If this upper bound is no better than the best achieved solution, then the time-consuming computation for the IP model is skipped and RR returns an empty solution which consists of m empty routes and a random ordered list of all unvisited customers. This results in considerable computational efforts. 11

When RR is solved by an MIP solver, the best combination of the routes in POOL is obtained. We then construct a complete solution based this combination by attempting to insert all unvisited customers. In the implementation, the customers in u are sorted in non-ascending order of profit. Let ZRR be the final solution. ZRR is usually of high quality and can be used to assist LS and SA in escaping from a local optimum in the iterative framework introduced in the following section. 4. An iterative framework To ensure better cooperation among the three components, we introduce an iterative framework as the entire algorithm described in Algorithm 2. The whole algorithm is called the Iterative 3Component Heuristic (I3CH). Algorithm 2 Iterative 3-Component Heuristic for the TOPTW Require: maximum number of iterations Imax ; Require: integer N; Require: route pool POOL; 1: Apply the Initialization to obtain starting solution A, and cache the routes into POOL; 2: Set i ← 1; 3: while i ≤ Imax do 4:

Invoke the route recombination over POOL and obtain ZRR ;

5:

Apply the local search to explore N neighbors and cache the routes into POOL;

6:

Run the simulated annealing with N steps and cache the routes into POOL;

7:

Select the best solution B from the collection {A, ZRR , XLS , YSA } breaking ties arbitrarily, where XLS and YSA are the best solutions of LS and SA respectively;

8:

A ← B, i ← i + 1;

9:

If all customers are served in A, then set i ← Imax + 1;

10: end whilethe solution of A;

At the beginning of the proposed heuristic, we apply an initialization procedure to obtain a starting solution A. During initialization, it randomly enumerates 3 ∗ N solutions and selects the best one to add to set A. Additionally, it stores the routes of the enumerated solutions into 12

POOL for future use. To enumerate a solution, the procedure first initializes m empty routes and a randomly ordered list u containing all customers. It then invokes the eliminator which can eliminate nothing but can choose customers from u for the routes before a complete solution is obtained. Applying the post-processing procedure in the third stage produces an improved solution. The routes from such solution are cached into POOL. At the end of this procedure, we set A as the starting solutions of LS and SA. The iterative framework significantly strengthens our approach. As introduced, LS and SA deliver their routes to RR such that RR can produce a high quality solution ZRR . On the other hand, RR contributes and assists LS and SA to escape from local optima. In this framework, the starting solution of LS is to be reset to the global best solution A once the number of consecutive iterations without improvement exceeds a number Ino impr , i.e., Ils no impr > Ino impr . Similarly, the starting solution of SA is also updated by A. The difference is that the condition for SA is Isa no impr > Ino impr . As previously observed, the best global solution A is frequently updated by ZRR . In this sense, RR regularly brings improvements to LS and SA. The iterative framework highlights the management issues with POOL. On the one hand, LS and SA repeatedly push routes into POOL. On the other, the efficiency of RR depends on the size of POOL. Setting pool size S pool to a fixed number means that new routes cannot be immediately added when the capacity is full. This means that we have to remove some “old” routes out of POOL before adding new ones in an iteration. Furthermore, the routes in POOL determine the quality of the solution found by RR. Hence, it is important to devise an appropriate strategy for managing the route pool POOL. Our strategy is to adapt the least recently used (LRU) concept that discards the least recently used route in our application. By preserving the most recently used routes in POOL and exchanging the least recently used routes with the newly discovered routes, the LRU strategy increases RR’s opportunities to improve solution quality. We define an attribute apr for each route in POOL to indicate its importance and set each new route included into POOL with apr = 0. Because of POOL’s capacity, we eliminate routes with the least apr and reallocate the space to new routes. In the case of a tie, the more profitable routes are retained. Once RR is solved and returns a feasible solution, we update the attributes of the routes employed with apr = 100. The apr of the other routes in POOL is reduced 13

by 1. 5. Computational Experiments In computational experiments, our algorithm was implemented as a sequential algorithm in Java. The experimental results reported were obtained on a Linux server equipped with an Intel Xeon E5430 CPU clocked at 2.66 GHz, 8GB RAM and running CentOS 5.4. The IP solver used was ILOG CPLEX 12.2 (64-bit Linux edition). The detailed results can be found online at http://www.computational-logistics.org/orlib/toptw/. 5.1. Test instances and approach comparison Vansteenwegen et al. (2011) summarized all existing test instances of the TOPTW in their review paper. The data sets can be accessed at http://www.mech.kuleuven.be/en/cib/ op. Montemanni and Gambardella (2009) constructed the TOPTW instances that extend the OPTW instances by increasing the number of vehicles: m = 2, 3, 4. The earlier OPTW instances were designed by Righini and Salani (2009) based on the 48 instances in Solomon (1987) VRPTW data set (c* 100, r* 100 and rc* 100) and 10 instances in the multi-depot vehicle routing problems (pr01 - pr10) considered by Cordeau et al. (1997). Montemanni and Gambardella (2009) developed another 37 OPTW instances, of which 27 instances are adapted from Solomon (1987)’s data set (c* 200, r* 200 and rc* 200) and 10 instances from that of Cordeau et al. (1997) (pr11-pr20). We classify these TOPTW instances into the aforementioned “INST-M” category which contains four subsets which are “Cordeau 1-10”, “Cordeau 11-20”, “Solomon 100” and “Solomon 200”. Vansteenwegen et al. (009b) constructed a new data set of TOPTW instances which were extensions of those in Solomon (1987) and Cordeau et al. (1997). The major feature of this new data set is that the optimal solution for each instance is known. This is due to the specific setting on the number of provided vehicles, which enables all customers to be visited. Hence, the optimal objective value is equal to the total profit on the network graph. We include these instances into the “OPT” category.

14

To evaluate the performance of our algorithm, we considered the following state-of-the-art algorithms in our comparisons: • ACO: the ant colony system proposed by Montemanni and Gambardella (2009); • ILS: the iterated local search algorithm developed by Vansteenwegen et al. (009b); • VNS: the variable neighbor search approach proposed by Tricoire et al. (2010); • GRASP-ELS: a hybrid metaheuristic combining the greedy randomized adaptive search procedure with the evolutionary local search, proposed by Labadie et al. (2011); • SSA: slow version of the simulated annealing heuristic proposed by Lin and Yu (2012); • GVNS: the LP-based granular variable neighborhood search algorithm by Labadie et al. (2012). Of these algorithms, only ACO was not tested on the new data set “OPT”. Also, ACO was executed in five runs for each instance whereas the VNS and GVNS were executed in 10. We considered their average performance in the comparisons. Moreover, as previously noted, Lin and Yu (2012)’s SA has fast and slow versions. The fast simulated annealing (FSA) which is given less computational time and the slow simulated annealing (SSA) which finds better solutions at the expense of more computational time. As we are more concerned with solution quality, we used SSA rather than FSA in the comparisons. To ensure fair comparisons, the solutions for each algorithm were compared with the best known solutions (BKS) among ACO, ILS, VNS, GRASP-ELS, SSA, GVNS and the improved best solutions from Gambardella et al. (2012) with the computational time adjusted to the speed of the computers used in those solutions. We summarize the experimental environment of each algorithm and compare their CPU speed in Table 1. As all of the algorithms are single threaded, it is reasonable to compare their CPU speed using the Super Pi benchmark 1 . In Table 1, the Super Pi 1

Super Pi is a single-threaded program that computes π up to a specified number of digits using the Gauss-Legendre

algorithm. It is commonly used as a crude estimate of CPU speed. See http://www.superpi.net/

15

(s) column reports the number of seconds it took a processor to compute the first 1 million digits of π. Unfortunately, for ACO and VNS, the Super Pi scores of their processors are not available. Labadie et al. (2011) presents strong evidence to show the comparability of processors used by ACO and GRASP-ELS. As only limited information was available on the processor used by VNS method, as such we cannot estimate its speed by considering only the clock-rate. A review of all 2.40 GHz processors found most to be slower than ours, and we therefore assumed that the processor used in computing the VNS was slower than ours. Finally, we estimated the singlethread performance of each processor by supposing the performance of our machine to be 1. The computational time was adjusted by the single-thread performance. For example, we multiplied the computational time of a solution obtained by GVNS by 0.33. Table 1: Estimation of single-thread performance. Algorithm

Experimental environment

Super Pi

Estimate of single-thread performance

ACO

Dual AMD Opteron 250 2.4GHz CPU, 4GB RAM

unknown

0.33

ILS

Intel Core 2 2.5 GHz CPU, 3.45GB RAM

18.6

0.79

VNS

2.4 GHz CPU, 4 GB RAM

GRASP-ELS

unknown

<1

Intel Pentium 4 processor, 3.00 GHz, 1 GB RAM

44.3

0.33

SSA

Intel Core 2 CPU, 2.5 GHz

18.6

0.79

GVNS

Intel Pentium (R) IV, 3 GHz CPU

44.3

0.33

I3CH

Intel Xeon E5430 CPU clocked at 2.66 GHz, 8GB RAM

14.7

1

5.2. Parameter tuning The proposed approach I3CH has nine parameters: Imax , N and Ino impr in the iterative framework, Ph and Pl in the neighborhood operator eliminator, T 0 and α in the SA component, S pool in the RR component and a random seed, denoted as seed, to ensure randomness. As appropriate settings on these parameters strengthen our algorithm, we designed the parameter tuning experiments. We carried out four experiments to tune the parameters, using a predetermined subset of 10 instances. The chosen instances are “c203”, “c207”, “pr02”, “pr07”, “pr12”, “pr16”, “r102”, “r105”, “rc107” and “rc204” considering m = 4. Our approach performed relatively worse on these 10 instances in the preliminary tests in which the parameters were set to preliminary values. 16

We first performed tests to adjust the settings on Ph and Pl in the neighborhood operator eliminator. In our design, we invoked the I3CH procedure on the 10 chosen instances to test all probabilities Ph ∈ {0.0, 0.1, 0.2, 0.3, 0.5} and Pl ∈ {0.1, 0.2, 0.4, 0.6, 0.8, 1.0}. The experiment was repeated five times with different random seeds seed ∈ {3, 5, 7, 9, 11}. For the remaining parameters, we set T 0 = 1 and α = 0.92 for the SA component, pool size S pool = 300, N = 20, Ino impr = 20 and Imax = 1000. For each run, the percentage gap between the solution value achieved by the I3CH and the BKS is computed by Equation (6). gap =

BKS − I3CH × 100% BKS

(6)

The results are summarized in Table 2 in which we present the average gap over 50 runs (10 instances × 5 random seeds) for each pair of Ph and Pl . We can see that the I3CH performs worse when Ph = 0 which indicates that allowing more profitable customers to be eliminated affords greater flexibility and increases the opportunities to get a better solution. We finally decided to pick Ph = 0.1 and Pl = 0.3 as this pair resulted in the smallest average gap. Table 2: Parameter tuning on the neighborhood operator Eliminator (Avg. Gap (%)). Ph

Pl = 0.1

Pl = 0.2

Pl = 0.3

Pl = 0.4

Pl = 0.6

Pl = 0.8

Pl = 1

0

2.64

2.16

1.99

2.13

1.84

2.28

4.42

0.1

0.71

0.60

0.51

0.75

0.99

1.46

2.40

0.2

0.63

0.59

0.67

0.78

1.11

1.53

2.44

0.3

0.55

0.56

0.71

0.90

1.23

1.88

2.58

0.5

0.99

1.45

1.27

1.52

1.86

2.19

2.88

A second experiment was carried out to tune parameters T 0 and α in the SA component. Similar to the first experiment, we invoked the I3CH procedure on the 10 instances and repeated it five times with different random seeds seed ∈ {3, 5, 7, 9, 11}. We tested T 0 ∈ {0.1, 1, 10, 100}, α ∈ {0.9, 0.93, 0.95, 0.97, 0.99, 0.995, 0.999}. The settings on the remaining parameters were Ph = 0.1, Pl = 0.3, S pool = 300, N = 20, Ino impr = 20 and Imax = 1000. We computed the average percentage gap between the solution value by I3CH and BKS over five different random seeds for each test instance. In Table 3, we present the average of this gap over the 10 test instances 17

for each T 0 and α. It is shown both T 0 = 0.1, α = 0.995 and T = 100, α = 0.999 achieved the smallest gap 0.39%. Therefore, we picked the pair T 0 = 0.1, α = 0.995. Table 3: Parameter tuning on SA (Avg. Gap (%)). T0

α = 0.9

α = 0.93

α = 0.95

α = 0.97

α = 0.99

α = 0.995

α = 0.999

0.1

0.54

0.58

0.41

0.65

1

0.62

0.48

0.60

0.60

0.47

0.39

0.46

0.50

0.54

0.42

10

0.47

0.62

0.60

0.53

0.50

0.41

0.53

100

0.51

0.56

0.60

0.48

0.48

0.47

0.39

The third experiment tuned S pool in the RR component. Ph = 0.1, Pl = 0.3, T 0 = 0.1 and α = 0.995 were determined. We set seed ∈ {3, 5, 7, 9, 11}, N = 20, Ino impr = 20, and Imax = 1000 and aimed to select the most appropriate S pool ∈ {100, 200, 300, 400, 600, 800, 1000, 1500, 2000}. Table 4 reports the average performance for each setting of S pool . We chose S pool = 1000, which is a tradeoff between effectiveness and efficiency. Table 4: Parameter tuning on pool size S pool . S pool

Avg. Gap (%)

Avg. Time (s)

100

0.66

31.0

200

0.48

32.6

300

0.40

33.4

400

0.56

34.6

600

0.34

37.7

800

0.35

41.0

1000

0.27

44.5

1500

0.35

58.4

2000

0.25

79.2

The final experiment was designed to investigate the tuning of N, Ino impr and Imax . These three parameters affect the performance of the I3CH from the outer iterative framework. We considered N ∈ {10, 20, 50, 100}, Ino impr ∈ {10, 20, 50, 100} and Imax ∈ {500, 1000, 2000, 3000, 4000, 5000}. The other parameter settings had already been configured as Ph = 0.1, Pl = 0.3, T 0 = 0.1, α = 18

0.995 and S pool = 1000. We repeated the experiment five times with different seed ∈ {3, 5, 7, 9, 11} and computed their average performance. Figure 4 contains four subgraphs, each of which investigates the performance of the I3CH over different N and Imax by fixing Ino impr . In the subgraph for “Ino impr = 10”, we see that the curve marked “N = 50” outperforms the other three over all Imax . The three other subgraphs, also showed that one group of settings outperforming the others. We therefore selected the best four, the results of which are summarized in Table 5. The table shows that Ino impr = 20 and N = 50 constitute the best choice with regard to solution quality. Column Avg. Time (s) shows that we chose Imax = 3000 since 380 seconds for an instance is acceptable. This setting is also a tradeoff

0.90

0.90

0.80

0.80

0.70

0.70

Aver. Gap (%)

Aver. Gap (%)

between effectiveness and efficiency.

0.60 0.50

N = 10

0.40

N = 20 N = 50

0.30

N = 100

0.20 0.10 500

1000

2000

3000

4000

5000

0.60 0.50

N = 10

0.40

N = 20 N = 50

0.30

N = 100

0.20 0.10

Imax

500

1000

0.90

0.80

0.80

0.70

0.70

Aver. Gap (%)

Aver. Gap (%)

0.90

0.60 0.50

N = 10

0.40

N = 20 N = 50

0.30

N = 100

0.20 0.10 500

1000

2000

3000

2000

3000

4000

5000

Imax

Ino_impr = 20

Ino_impr = 10

4000

5000

0.60 0.50

N = 10

0.40

N = 20 N = 50

0.30

N = 100

0.20 0.10

Imax

500

Ino_impr = 50

1000

2000

3000

Ino_impr = 100

Figure 4: Parameter tuning on Ino impr , N and Imax .

19

4000

5000

Imax

Table 5: Parameter tuning on Ino impr , N and Imax . Avg. Gap (%)

Avg. Time (s)

Ino impr

N

10

50

20

50

0.35

0.32

0.21

0.18

0.16

0.16

67.1

131.1

258.2

384.9

512.1

638.9

50

50

0.52

0.45

0.40

0.37

0.33

0.29

68.6

133.4

262.8

392.9

521.7

650.0

100

50

0.45

0.39

0.37

0.34

0.32

0.30

68.1

131.1

257.8

384.2

510.6

637.7

Imax = 500

1000

2000

3000

4000

5000

Imax = 500

1000

2000

3000

4000

5000

0.48

0.38

0.29

0.26

0.26

0.24

68.1

133.0

261.6

391.2

520.0

650.0

5.3. Analysis of components In this subsection, we discuss the results of experimental tests carried out to investigate the behavior of the proposed approach. As our approach is essentially component-based, it is important to investigate the performance of each component and their combinations. These experiments allowed us to analyze the effects exerted by the components. In the experiments, we considered both each component alone and their various combinations, namely, LS, SA, LS+SA, LS+RR, SA+RR and LS+SA+RR. All six methods are embedded in the iterative framework. Note that LS+SA+RR is actually the I3CH. To investigate their performance, we tested each of them on all the “INST-M” instances (m = 1, 2, 3, 4). The methods share the same parameter settings which are Ph = 0.1 and Pl = 0.3 for the neighborhood operator, T 0 = 0.1 and α = 0.995 for SA, S pool = 1000 for RR; Imax = 3000, N = 50 and Ino impr = 20 for the iterative framework. The random seed is seed = 3. We summarize the average performance of each method in Table 6. The performance is evaluated on solution quality by computing the percentage gap from BKS for each instance. Each column under a method gives the average performance on each instance set and the overall average performance on all instances. From the table, LS performed the worst with an average gap 1.50% from BKS. SA performed better than LS on each instance set. This is due to the superior convergence properties of SA and its capability to escape local optima. The combination LS+SA achieved better solution quality than both LS and SA, but the improvement is slight. However, RR brought improvements. After combining RR with LS, SA and LS+SA respectively, the average performance is significantly improved. For instance, we can see that LS+RR reduced the average 20

gap obtained by LS by 0.45%. Table 6: Average performance of the combinations of the components on “INST-M” instances (%). Set

LS

SA

LS + SA

LS + RR

SA + RR

LS + SA + RR

Cordeau 1-10

2.42

1.34

1.53

1.67

0.97

0.60

Cordeau 11-20

3.64

2.94

2.79

3.06

2.17

1.71

Solomon 100

1.18

0.79

0.68

0.60

0.46

0.34

Solomon 200

0.72

0.64

0.52

0.57

0.59

0.44

All

1.50

1.10

1.01

1.05

0.80

0.59

Figure 5 shows the average computational time of each component for the six methods. Both LS and SA took 90 seconds, on average, to compute an instance. Incorporation of the RR component increased the average computational time by about 40 seconds. Given that this component also brought about an improvement in solution quality, this additional 40 seconds well worth the slightly longer computing time.

CPU Sceonds

200 150 Time on RR

100

Time on SA 50

Time on LS

0

Figure 5: Analysis of computational time on the components.

Furthermore, we also investigated the performance of the LS component by allowing it to compete with the solutions achieved by ILS and GRASP-ELS. The comparison results are presented in Table 7. Table 8 compares the performance of FSA and SSA with our SA component. 21

The columns headed AG (%) give the average percentage gap from the BKS, and those headed columns AT (s) show the average computational time in CPU seconds. The computational time of an approach was adjusted according to its computer’s speed. In Table 7, the solution quality achieved by LS is much better than ILS, but our method takes longer computational time. Compared to GRASP-ELS, LS performs worse. With the efforts from the RR component, LS+RR is able to achieve better solutions in average than GRASP-ELS. Table 7: Comparison of “LS” methods. ILS

GRASP-ELS

LS + RR

LS

Set AG (%)

AT (s)

AG (%)

AT (s)

AG (%)

AT (s)

AG (%)

AT (s)

Cordeau 1-10

6.04

5.9

1.92

9.1

2.42

143.8

1.67

172.0

Cordeau 11-20

8.38

6.0

3.01

12.0

3.64

176.0

3.06

223.6

Solomon 100

2.34

1.0

1.16

9.1

1.18

55.0

0.60

63.2

Solomon 200

1.76

1.4

0.55

3.4

0.72

88.5

0.57

121.0

All

3.41

2.4

1.29

7.5

1.50

94.5

1.05

119.2

The approaches in Table 8 all apply the simulated annealing technique. From the comparisons in this table, our SA method achieves better solutions on the instance sets “Cordeau 1-10” and “Cordeau 11-20” with a similar amount of computational time. On instance sets “Solomon 100” and “Solomon 200”, SSA outperforms our SA method. The RR component recombines the routes from the SA solution and produces high-quality solutions. As a result, SA+RR outperforms the other three approaches in Table 8 in terms of overall performance. Table 8: Comparison of “SA” methods. FSA

SSA

SA + RR

SA

Set AG (%)

AT (s)

AG (%)

AT (s)

AG (%)

AT (s)

AG (%)

AT (s)

Cordeau 1-10

5.35

5.9

1.88

145.9

1.34

147.2

0.97

173.0

Cordeau 11-20

7.61

6.0

3.43

177.7

2.94

179.7

2.17

210.8

Solomon 100

2.80

1.0

0.26

31.8

0.79

61.7

0.46

60.3

Solomon 200

1.74

1.4

0.56

42.3

0.64

79.4

0.59

113.2

All

3.39

2.4

1.00

69.8

1.10

94.8

0.80

113.7

These experiments show that the RR component plays a key role in our approach. Integrating 22

the RR component strengthens our approach considerably. In addition, the average computational time spent on the RR component is less than that for either the LS component or SA component. 5.4. Comparisons with other approaches In this subsection, we compare the results of the I3CH with those of the ACO, ILS, VNS, GRASP-ELS, SSA and GVNS on all TOPTW instances in the “INST-M” and “OPT”. The I3CH was configured with determined parameter settings: Ph = 0.1 and Pl = 0.3 for the neighborhood operator; T 0 = 0.1 and α = 0.995 for SA; S pool = 1000 for RR; and Imax = 3000, N = 50 and Ino impr = 20 for the iterative framework. Observing that the average performance varied little by random seed, we selected seed = 3 for the final run. All of the detailed solutions obtained by the I3CH are presented in the appendix. Table 9 summarizes the results for the instances in “INST-M”. It compares the average performance of I3CH with that of the ACO, ILS, VNS, GRASP-ELS, SSA and GVNS. The num column gives the number of instances in a set, and the AG(%) column reports the average percentage gap with the BKS. A negative percentage gap indicates an average improvement. The AT(s) column reports the average computational time in CPU seconds. The computational times of the approaches were adjusted according to computers’ speed, as discussed in Section 5.1. For the I3CH, we add additional column #Impr to show the number of new best solutions it found. We can see from Table 9 that, on average, the I3CH improved the solution quality of 10 out of the 16 instance sets. For the OPTW instances, the performance of I3CH is not good enough. This is because the RR component degenerates to pick the route with the highest profit when m = 1. In this case, I3CH is equivalent to LS+SA. The table also shows that our LS+SA is not as competitive as SSA or VNS. As m becomes larger, our dedicated approach, I3CH, becomes quite effective in solving the TOPTW instances. Overall, it outperforms the other five approaches by at least 0.41% in terms of average performance over all instances from “INST-M”. Further, its average computational time is 200 seconds, which is acceptable. The I3CH also discovered 35 new best solutions which are summarized in Table 10. I3CH was also tested on the instances in the “OPT” category. Table 11 reports the overall comparisons with ILS, GRASP-ELS, SSA and GVNS. In this table, the overall performance of each 23

Table 9: Overall comparison on “INST-M”.

ACO num

Set

ILS

VNS

AG (%)

AT (s)

AG (%)

AT (s)

AG (%)

GRASP-ELS

AT (s)

AG (%)

SSA

AT (s)

AG (%)

GVNS

AT (s)

AG (%)

I3CH

AT (s)

#impr

AG (%)

AT (s)

m=1 Cordeau 1-10

10

1.20

536.8

4.72

1.4

1.08

< 822.1

1.44

1.7

0.97

88.6

1.61

4.1

0

1.05

109.0

Cordeau 11-20

10

11.47

292.9

9.11

1.6

2.92

< 1045.9

2.92

2.6

3.25

128.3

4.06

10.3

3

3.79

130.2

Solomon 100

29

0.10

66.0

1.94

0.2

0.07

< 85.4

0.20

3.0

0.05

17.6

2.46

22.0

0

0.69

26.7

Solomon 200

27

0.89

283.1

2.87

1.3

0.89

< 857.8

1.49

5.5

0.85

35.3

2.87

24.9

1

1.34

132.1

m=2 Cordeau 1-10

10

3.41

623.6

6.06

3.8

3.72

< 524.8

2.04

6.4

2.29

137.4

1.62

12.9

1

0.94

247.1

Cordeau 11-20

10

6.14

787.0

7.84

4.1

3.62

< 618.8

3.56

9.5

3.87

159.3

2.17

27.2

1

2.69

304.6

Solomon 100

29

0.71

422.0

1.93

0.7

1.05

< 68.8

1.36

8.8

0.13

27.3

1.72

24.4

0

0.47

69.3

Solomon 200

27

3.14

733.5

3.06

2.1

0.97

< 813.7

0.65

6.9

0.92

60.7

1.42

6.5

4

0.43

463.8

m=3 Cordeau 1-10

10

3.99

714.1

6.57

7.3

3.62

< 473.2

1.97

13.4

2.33

155.6

1.13

28.3

1

0.35

424.0

Cordeau 11-20

10

6.50

735.3

9.08

7.7

3.24

< 517.5

3.15

14.2

3.70

198.9

1.69

49.7

3

1.00

497.0

Solomon 100

29

1.41

469.2

2.36

1.2

1.19

< 68.9

1.56

11.6

0.40

36.3

1.83

30.1

2

0.16

135.9

Solomon 200

27

0.83

443.8

1.10

1.4

0.12

< 322.3

0.06

1.2

0.46

41.3

0.25

2.4

1

-0.01

89.2

m=4 Cordeau 1-10

10

3.30

807.7

6.80

11.1

3.28

< 403.2

2.23

15.1

1.93

201.9

1.33

42.0

4

0.05

566.5

Cordeau 11-20

10

5.35

852.6

7.50

10.9

2.46

< 408

2.40

21.6

2.92

224.3

1.77

76.8

6

-0.64

728.6

Solomon 100

29

1.66

502.6

3.11

1.9

1.56

< 66.8

1.53

13.2

0.47

46.1

1.84

28.6

8

0.06

199.5

Solomon 200

27

0.02

81.0

0.00

0.8

0.00

< 141.2

0.00

0.0

0.00

32.0

0.00

0.2

0

0.00

0.2

Grand Total

304

2.17

452.1

3.41

2.4

1.33

< 375.6

1.29

7.5

1.00

69.8

1.66

21.3

35

0.59

200.9

Table 10: New best solutions discovered by I3CH. Name

m

New Best

Name

m

New Best

Name

m

New Best

Name

m

New Best

pr11

1

353

pr13

1

466

rc201

2

1384

pr03

4

1232

pr20

4

2062

c103

3

990

pr04

4

1585

r103

4

pr15

1

928

707

pr08

3

1139

pr05

4

1838

r106

4

906

r202 pr05

1

930

pr13

3

1145

pr10

4

1943

r107

4

950

2

1101

pr15

3

1654

pr13

4

1386

r110

4

915

pr13

2

832

pr17

3

841

pr15

4

2065

r111

4

952

r201

2

1254

rc104

3

834

pr17

4

934

rc103

4

974

r203

2

1416

rc205

3

1719

pr18

4

1539

rc104

4

1064

r205

2

1380

c103

4

1210

pr19

4

1750

24

approach is given by three columns. The columns #OPT show the number of optimal solutions achieved. The columns AG (%) and the columns AT (s) give the average percentage gap from the optimal and the average computational time in CPU seconds respectively. Table 11: Overall comparison on “OPT”.

ILS num

Set

GRASP-ELS

SSA

GVNS

I3CH

#OPT

AG (%)

AT (s)

#OPT

AG (%)

AT (s)

#OPT

AG (%)

AT (s)

#OPT

AG (%)

AT (s)

#OPT

AG (%)

AT (s)

Cordeau 1-10

10

1

2.32

24.0

5

0.92

23.6

4

1.04

447.1

3

1.25

16.9

6

0.78

326.7

Solomon 100

29

4

1.80

2.5

11

0.55

21.7

8

0.59

71.7

6

1.14

6.4

25

0.03

393.8

Solomon 200

27

20

0.39

1.2

23

0.02

1.2

22

0.08

38.3

21

0.12

1.1

24

0.04

127.2

Grand Total

66

25

1.30

5.2

39

0.39

13.6

34

0.45

114.9

30

0.74

5.8

55

0.15

274.5

Table 11 shows that the I3CH achieved the smallest average gap over each of the seven instance sets. The average percentage gap over all “OPT” instances achieved by the I3CH was 0.15%, a 0.24% improvement over the best exiting result. In addition, the approach obtained 55 optimal solutions, 16 more than the previous best approach in the literature. Column num shows that there are only 66 OPT instances in total. Overall, the I3CH outperformed the existing approaches on the “OPT” instances, although it required more computational time. On the whole, our I3CH approach is rather competitive. It discovered new best solutions and achieved the best average performance. Although our approach may at first glance appear less efficient than some of the other approaches, requiring more computational time, this is not actually the case. We conducted comparison experiments in which the I3CH was set to the same computational time as the other approaches which had been adjusted by their computers’ speed. These experiments were carried out on the “INST-M” instances, with our approached compared with the two best known approaches, GRASP-ELS and SSA. We examined the solution quality obtained by I3CH when the total execution time was set to the computational time required by GRASP-ELS and SSA. Tables 12(a) and 12(b) summarize the comparison results. Table 12(a) shows the I3CH approach’s execution time to be the same as GRASP-ELS. It outperformed GRASP-ELS on 9 of the 16 instance sets. GRASP-ELS performed quite well on the OPTW instances, whereas our dedicated approach achieved better performance in solving the 25

Table 12: Comparisons with the same computational time. (a) Comparisons to GRASP-ELS.

(b) Comparisons to SSA.

Avg. Gap (%) m 1

2

3

4

Set

Avg. Gap (%) Avg. Time (s)

GRASP-ELS

I3CH

m

Set

Avg. Time (s) SSA

I3CH

Cordeau1-10

1.44

2.56

1.7

Cordeau1-10

0.97

1.05

88.6

Cordeau11-20

2.92

4.21

2.6

Cordeau11-20

3.25

3.79

128.3

Solomon100

0.20

0.76

3.0

Solomon100

0.05

0.69

17.6

Solomon200

1.49

2.67

5.5

Cordeau1-10

2.04

1.70

6.4

Cordeau11-20

3.56

3.02

Solomon100

1.36

0.61

Solomon200

0.65

1.42

6.9

Cordeau1-10

1.97

0.93

13.4

Cordeau11-20

3.15

1.70

14.2

Solomon100

1.56

0.23

11.6

Solomon200

0.06

0.52

1.2

Cordeau1-10

2.23

1.00

15.1

Cordeau11-20

2.40

0.27

Solomon100

1.53

0.13

0.00

0.00

0.0

1.29

1.08

7.5

Solomon200 Overall

1

Solomon200

0.85

1.61

35.3

Cordeau1-10

2.29

0.94

137.4

9.5

Cordeau11-20

3.87

2.75

159.3

8.8

Solomon100

0.13

0.53

27.3

Solomon200

0.92

0.70

60.7

Cordeau1-10

2.33

0.45

155.6

Cordeau11-20

3.70

1.14

198.9

Solomon100

0.40

0.19

36.3

2

3

Solomon200

0.46

0.01

41.3

Cordeau1-10

1.93

0.14

201.9

21.6

Cordeau11-20

2.92

-0.54

224.3

13.2

Solomon100

0.47

0.15

46.1

0.00

0.00

32.0

1.00

0.67

69.8

4

Solomon200 Overall

TOPTW instances. Overall, the I3CH achieved an average percentage gap of 1.08%, which is 0.21% less than that of GRASP-ELS. Similarly, we set the execution time for the I3CH to the computational time for SSA for each instance. In comparing their average performance on the instance sets, we found the I3CH to display better performance on 10 of the 16 instance sets, and to largely outperform SSA on the TOPTW instances. Its overall average performance is 0.33% better than that of SSA. Our experimental results thus indicate that the I3CH approach proposed herein is both effective and efficient, and is well-suited to the TOPTW. 6. Conclusions and Future Work In this paper, we study the TOPTW and we propose an iterative framework that incorporates three components. This framework employs a local search procedure and a simulated annealing procedure as the first two components to explore the solution space. By tracing the neighborhood exploration, they discover a set of routes which are then recombined by the third component to 26

obtain a high quality solution. Embedded in an iterative framework, these three components cooperate and work as an effective heuristic. Our computational results indicate that the heuristic outperforms the existing approaches in the literature in terms of average performance. Our iterative three-component heuristic found 35 new best solutions and improved solution quality, on average, by 0.41% over the “INST-M” instances. It also achieved 55 optimal solutions over the 66 “OPT” instances, which is 16 more than the previous best approach in the literature. In addition, results indicate our approach is both efficient and effective. The new results we reported herein can serve as benchmarks for future studies. The algorithm we propose in this paper is highly effective. The framework is robust and can incorporate additional components. If the components were adapted to accommodate new constraints, it would also be applicable to other routing problems. One potential avenue for further research would be to extend our approach to a parallel version. Either classical parallel computing or recent popular GPU computing techniques would be promising candidates to speed up the solution process and improve solution quality. Acknowledgment The authors thank two anonymous referees, whose comments greatly improved the quality of the paper. T. Tsiligirides, Heuristic methods applied to orienteering, Journal of the Operational Research Society (1984) 797– 809. B. Golden, L. Levy, R. Vohra, The orienteering problem, Naval Research Logistics (NRL) 34 (3) (1987) 307–318. W. Souffriau, P. Vansteenwegen, J. Vertommen, G. Berghe, D. Oudheusden, A personalized tourist trip design algorithm for mobile tourist guides, Applied artificial intelligence 22 (10) (2008) 964–985. M. Schilde, K. Doerner, R. Hartl, G. Kiechle, Metaheuristics for the bi-objective orienteering problem, Swarm Intelligence 3 (3) (2009) 179–201. I. Chao, B. Golden, E. Wasil, et al., A fast and effective heuristic for the orienteering problem, European Journal of Operational Research 88 (3) (1996a) 475–489. M. Fischetti, J. Salazar Gonzalez, P. Toth, Solving the orienteering problem through branch-and-cut, INFORMS Journal on Computing 10 (1998) 133–148. P. Vansteenwegen, W. Souffriau, D. Oudheusden, The orienteering problem: A survey, European Journal of Operational Research 209 (1) (2011) 1–10.

27

I. Chao, B. Golden, E. Wasil, et al., The team orienteering problem, European journal of operational research 88 (3) (1996b) 464–474. H. Tang, E. Miller-Hooks, A tabu search heuristic for the team orienteering problem, Computers & operations research 32 (6) (2005) 1379–1407. P. Vansteenwegen, W. Souffriau, D. Van Oudheusden, A Detailed Analysis of Two Metaheuristics for the Team Orienteering Problem, Engineering Stochastic Local Search Algorithms. Designing, Implementing and Analyzing Effective Heuristics (2009a) 110–114. M. Kantor, M. Rosenwein, The orienteering problem with time windows, Journal of the Operational Research Society (1992) 629–635. G. Righini, M. Salani, Decremental state space relaxation strategies and initialization heuristics for solving the Orienteering Problem with Time Windows with dynamic programming, Computers & operations research 36 (4) (2009) 1191–1203. R. Montemanni, L. Gambardella, Ant colony system for team orienteering problems with time windows, Foundations of Computing and Decision Sciences 34 (2009) 287–306. R. Montemanni, D. Weyland, L. Gambardella, An Enhanced Ant Colony System for the Team Orienteering Problem with Time Windows, in: Computer Science and Society (ISCCS), 2011 International Symposium on, IEEE, 381– 384, 2011. L. Gambardella, R. Montemanni, D. Weyland, Coupling ant colony systems with strong local searches, European Journal of Operational Research 220 (3) (2012) 831 – 843, ISSN 0377-2217, doi:10.1016/j.ejor.2012.02.038. P. Vansteenwegen, W. Souffriau, G. Vanden Berghe, D. Van Oudheusden, Iterated local search for the team orienteering problem with time windows, Computers & operations research 36 (12) (2009b) 3281–3290. F. Tricoire, M. Romauch, K. Doerner, R. Hartl, Heuristics for the multi-period orienteering problem with multiple time windows, Computers & Operations Research 37 (2) (2010) 351–367. W. Souffriau, P. Vansteenwegen, G. Berghe, D. Van Oudheusden, The Multiconstraint Team Orienteering Problem with Multiple Time Windows, Transportation Science (2011) Published in Andvance. N. Labadi, J. Melechovsk`y, R. Calvo, An effective hybrid evolutionary local search for orienteering and team orienteering problems with time windows, Parallel Problem Solving from Nature–PPSN XI (2011) 219–228. S. Lin, V. Yu, A Simulated Annealing Heuristic for the Team Orienteering Problem with Time Windows, European Journal of Operational Research 217 (1) (2012) 94–107. N. Labadie, R. Mansini, J. Melechovsk, R. W. Calvo, The Team Orienteering Problem with Time Windows: An LPbased Granular Variable Neighborhood Search, European Journal of Operational Research 220 (1) (2012) 15 – 27, ISSN 0377-2217, doi:10.1016/j.ejor.2012.01.030. M. Malek, M. Guruswamy, M. Pandya, H. Owens, Serial and parallel simulated annealing and tabu search algorithms for the traveling salesman problem, Annals of Operations Research 21 (1) (1989) 59–84.

28

A. Van Breedam, Improvement heuristics for the vehicle routing problem based on simulated annealing, European Journal of Operational Research 86 (3) (1995) 480–490. W. Chiang, R. Russell, Simulated annealing metaheuristics for the vehicle routing problem with time windows, Annals of Operations Research 63 (1) (1996) 3–27. S. Lin, V. Yu, S. Chou, Solving the truck and trailer routing problem based on a simulated annealing heuristic, Computers & Operations Research 36 (5) (2009) 1683–1692. V. Yu, S. Lin, W. Lee, C. Ting, A simulated annealing heuristic for the capacitated location routing problem, Computers & Industrial Engineering 58 (2) (2010) 288–299. R. Eglese, Simulated annealing: a tool for operational research, European journal of operational research 46 (3) (1990) 271–281. C. Koulamas, S. Antony, R. Jaen, A survey of simulated annealing applications to operations research problems, Omega 22 (1) (1994) 41–56. B. Suman, P. Kumar, A survey of simulated annealing as a tool for single and multiobjective optimization, Journal of the operational research society 57 (10) (2005) 1143–1160. M. Gendreau, J. Potvin, Handbook of metaheuristics, Springer, 2010. M. Solomon, Algorithms for the vehicle routing and scheduling problems with time window constraints, Operations research (1987) 254–265. J. Cordeau, M. Gendreau, G. Laporte, A tabu search heuristic for periodic and multi-depot vehicle routing problems, Networks 30 (2) (1997) 105–119. N. Labadie, J. Melechovsk`y, R. Wolfler Calvo, Hybridized evolutionary local search algorithm for the team orienteering problem with time windows, Journal of Heuristics 17 (6) (2011) 729–753.

29

x x x x x

Developed an iterative three-component heuristic for the TOP with time windows /ŶƐƉŝƌĞĚďLJƚǁŽ͞ŶŽ-worse-ƚŚĂŶ͟ƉƌŽƉŽƐŝƚŝŽŶƐǁŚĞŶĚĞƐŝŐŶŝŶŐƚŚĞĂůŐŽƌŝƚŚŵ Combined local search, meta-heuristic and mathematical programming Employed a dedicated neighborhood operator Achieved best overall performance and new best solutions on the benchmark