Youcef Amarouche
email: youcef.amarouche@hds.utc.fr

Rym Nesrine Guibadj
email: rym.guibadj@univ-littoral.fr

Elhadja Chaalal
email: de_chaalal@esi.dz

Aziz Moukrim
email: aziz.moukrim@hds.utc.fr

Effective Neighborhood Search with Optimal Splitting and Adaptive Memory for the Team Orienteering Problem with Time Windows

Keywords: Routing, Team Orienteering Problem, Time windows, Adaptive memory, Splitting algorithm, Heuristics

The Team Orienteering Problem with Time Windows (TOPTW) is an extension of the well-known Orienteering Problem. Given a set of locations, each one associated with a profit, a service time and a time window, the objective of the TOPTW is to plan a set of routes, over a subset of locations, that maximizes the total collected profit while satisfying travel time limitations and time window constraints. Within this paper, we present an effective neighborhood search for the TOPTW based on (1) the alternation between two different search spaces, a giant tour search space and a route search space, using a powerful splitting algorithm, and

(2) the use of a long term memory mechanism to keep high quality routes encountered in elite solutions.

We conduct extensive computational experiments to investigate the contribution of these components, and measure the performance of our method on literature benchmarks. Our approach outperforms state-of-theart algorithms in terms of overall solution quality and computational time. It finds the current best known solutions, or better ones, for 89% of the literature instances within reasonable runtimes. Moreover, it is able to achieve better average deviation than state-of-the-art algorithms within shorter computation times.

Moreover, new improvements for 57 benchmark instances were found.

Introduction

In the Team Orienteering Problem (TOP), we are given a transportation network in which a starting and an ending point are specified. The network connects a set of points that correspond to customer locations.

Each one of them is associated with a profit and a service time. For each pair of locations, a travel time is specified. The aim of the problem is to find a fixed number of disjoint paths from the starting point to the final destination through a subset of locations, each not exceeding a given time limit, that maximizes the total profit collected from visiting customers. Fig. 1 shows an example of a TOP instance with two available vehicles, and its solution. Each customer is represented with a circle area proportional to its profit. Given two available vehicles, the selected customers and their service order are shown in Fig. 1b. In this paper, we consider the Team Orienteering Problem with Time Windows (TOPTW), a natural extension of the TOP motivated by different practical situations. Possible applications of the problem range from logistics [START_REF] Golden | The orienteering problem[END_REF][START_REF] Tang | A tabu search heuristic for the team orienteering problem[END_REF] to leisure related applications like tourism (Vansteenwegen et al., 2009). In the TOPTW, each customer must be visited within a predefined time interval, specified by an earliest and a latest time, into which the service must start. We assume that the time windows are hard constraints. This means that early arrivals to a location are permitted, but the agent must wait for it to be "open" before the service can start. Late arrivals, however, are not allowed.

Herein, we propose an effective approach that follows the basic structure of a Multi-Start Iterated Local

Search (MS-ILS) to solve the TOPTW. We design a rather straightforward method that is able to effectively explore the search space to achieve high-quality solutions within very short computational times.

• First, we investigate the alternation between two search spaces: a giant tour search space and a route search space. In the route space, solutions are represented as genuine TOPTW solutions i.e., a set of feasible routes, one for each vehicle in use, while in the giant tour space, they are represented with an ordered list of customers with no route delimiters. The transition from the giant tour search space to the route search space is achieved using a powerful split algorithm. This idea is a follow up on preliminary work carried out by [START_REF] Guibadj | Memetic algorithm with an efficient split procedure for the team orienteering problem with time windows[END_REF]. They proposed a Memetic Algorithm for the TOPTW, which uses a giant tour representation for encoding individuals. MS-ILS, however, uses the giant tour representation more efficiently within an algorithm that is conceptually simpler and much faster than the MA.

• Second, we integrate an adaptive memory mechanism to overcome the drawbacks due to pure multistart heuristics being memoryless. The adaptive memory is used to store individual routes extracted from diverse high quality solutions. These routes are then combined to construct promising new starting solutions at each iteration of the multi-start algorithm.

• Third, we conduct extensive computational experiments to investigate the contribution of these components to the search performance of a basic MS-ILS, and measure the performance of our approach on literature benchmarks. The obtained results show that our MS-ILS outperforms state-of-the-art algorithms in terms of solution quality and computation time. It is able to find the current best-known solutions, or better ones, for 78% of the available benchmark instances. It achieves an overall average relative gap of 0.30% and 0.26% on the two benchmarks of the literature, respectively, while being faster than most state-of-the-art algorithms. Additionally, it was able to improve the solutions of 57 instances for which no optimal solutions have yet been found. In comparison, the previous best performing approach in the literature finds 65% of the current best-known solutions, and achieves a relative gap of 0.80% and 0.34% respectively, on the two benchmarks.

• Finally, we show that our algorithm can be tuned to either favor solution quality at the cost of more computational effort, or to considerably reduce computation times while maintaining good solution quality. As such, it can serve as a good basis for future developments on more complex variants of selective vehicle routing problems.

The remainder of this paper is organized as follows. Section 2 provides a brief overview of the literature related to the TOPTW, and Section 3 gives a formal description of the problem. Section 4 gives a detailed description of our proposed algorithm. First, the general framework of the method is introduced, and then each of its aspects is described in detail, including the solution representation, the optimal split procedure, and other components and parameters. The effectiveness of our approach is shown in Section 5. Finally, Section 6 concludes this paper and provides possible directions for future research.

Literature review

Orienteering problems are well-known NP-hard optimization problems, and solving them within a reasonable amount of time may prove to be rather difficult.

Exact solution methods

There are relatively few studies that deal with exact methods for the OPTW and TOPTW, and most of them are designed for the OPTW. [START_REF] Righini | Decremental state space relaxation strategies and initialization heuristics for solving the orienteering problem with time windows with dynamic programming[END_REF] solve the OPTW using a dynamic programming (DP) algorithm with Decremental State Space Relaxation (DSSR) [START_REF] Righini | New dynamic programming algorithms for the resource constrained elementary shortest path problem[END_REF], a relaxation of the problem so that customers can be visited more than once. The solution process consists in solving the DSSR, and then incrementally tightening the relaxed problem by disallowing multiple visits to a critical customer set of increasing size. [START_REF] Duque | Solving the orienteering problem with time windows via the pulse framework[END_REF] extend the pulse algorithm [START_REF] Lozano | On an exact method for the constrained shortest path problem[END_REF][START_REF] Cabrera | An exact bidirectional pulse algorithm for the constrained shortest path[END_REF] to solve the OPTW. The algorithm is presented as a general-purpose framework for hard shortest path problems, and includes two novel pruning strategies to discard sub-optimal and infeasible solutions. The proposed algorithm performs better than the dynamic programming approach of [START_REF] Righini | Decremental state space relaxation strategies and initialization heuristics for solving the orienteering problem with time windows with dynamic programming[END_REF]. The authors, however, only reported results on a subset of the benchmark instances.

The first exact method for the TOPTW was presented by [START_REF] Tae | A branch-and-price approach for the team orienteering problem with time windows[END_REF]. They proposed a Branch-&-Price (B&P) algorithm where the pricing problem is expressed as a Resource Constrained Elementary Shortest Path Problem (RCESPP). They tested their approach on the benchmark instances of [START_REF] Righini | New dynamic programming algorithms for the resource constrained elementary shortest path problem[END_REF], but omitted those of [START_REF] Montemanni | Ant colony system for team orienteering problems with time windows[END_REF], stating their difficulty, and those of Vansteenwegen et al. (2009) since their optimal solutions are already known. More recently, Gedik et al.

(2017) introduced an exact approach for the TOPTW based on a constraint programming (CP) formulation of the problem. Compared to the B&P, the CP-optimizer was able to solve benchmark instances from Montemanni and[START_REF] Montemanni | Ant colony system for team orienteering problems with time windows[END_REF]Vansteenwegen et al. (2009) which, were omitted by [START_REF] Tae | A branch-and-price approach for the team orienteering problem with time windows[END_REF], and provides good upper bounds within smaller runtimes, for the instance of [START_REF] Righini | New dynamic programming algorithms for the resource constrained elementary shortest path problem[END_REF]. However, the B&P is able to provide optimal solutions for more instances from the latter set.

Heuristic solution methods

Most of the research on orienteering problems with time windows focuses on the design of heuristic approaches. The fact that profits and travel times are independent, and that a good solution with respect to one criterion is often unsatisfactory with respect to the other, make it more difficult to devise consistently good heuristics to solve orienteering problems, despite them being seemingly simple [START_REF] Gendreau | A tabu search heuristic for the undirected selective travelling salesman problem[END_REF].

The presence of time window constraints increases this difficulty since it becomes harder to pinpoint the customers that should be included in the solution. Since the OPTW is a special case of the TOPTW where the fleet consists of a single vehicle, recent contributions tend to not make the distinction between single and multi-vehicle variants. Apart from the algorithm of Gunawan et al. (2015a), the heuristic methods that were introduced during the last decade were designed to solve both the OPTW and the TOPTW. [START_REF] Montemanni | Ant colony system for team orienteering problems with time windows[END_REF] proposed a heuristic method to solve the TOPTW using an Ant

Colony System (ACS) algorithm. The ACS was later improved by [START_REF] Montemanni | An enhanced ant colony system for the team orienteering problem with time windows[END_REF] The MuPOPTW is a generalization of the TOPTW where visits to customers span a period of time (days), and where multiple time windows are allowed. As such, each customer may be visited on a different day, and may have several time windows for each given day. To solve the MuPOPTW, the authors designed a Variable Neighborhood Search (VNS) algorithm, which also provides good solutions for the TOPTW.

However, their algorithm requires relatively large computation times. An algorithm that combines a Greedy Randomized Adaptive Search Procedure (GRASP) and an Evolutionary Local Search algorithm (ELS) was introduced by [START_REF] Labadie | Hybridized evolutionary local search algorithm for the team orienteering problem with time windows[END_REF]. Their method uses simple constructive heuristics inside the GRASP to generate distinct initial solutions, which are then improved using the ELS algorithm. In another study, [START_REF] Labadie | The team orienteering problem with time windows: An lp-based granular variable neighborhood search[END_REF] applied the concept of granularity to a VNS approach in order to reduce the size of the local search neighborhoods to moves that are more likely to lead to good quality solutions. Using the dual optimal solutions of an LP-problem, they partition the arc set into intervals of granularity from which they choose those with the most promising arcs. Two Simulated Annealing based algorithms for the TOPTW were presented by [START_REF] Lin | A simulated annealing heuristic for the team orienteering problem with time windows[END_REF]. The first one is a Fast Simulated Annealing (FSA) that is directed towards applications that require quick responses, while the second one is a Slow Simulated Annealing (SSA) that focuses on solution quality. [START_REF] Souffriau | The multiconstraint team orienteering problem with multiple time windows[END_REF] introduced the Multi-Constraint Team Orienteering

Problem with Multiple Time Windows (MC-TOP-MTW), a variant of the TOPTW that includes additional knapsack constraints to limit node selection. Such constraints can be used, for example, to model entrance fees in the case of tourist trip planning applications. The authors proposed a hybridization of GRASP and ILS, and evaluated the performance of their algorithm on benchmarks for the related problems, including the TOPTW. An iterative framework comprising three components, namely I3CH, was proposed by [START_REF] Hu | An iterative three-component heuristic for the team orienteering problem with time windows[END_REF]. The first two components are a Local Search (LS) and a Simulated Annealing (SA) that are used to explore the solution space and to discover a set of routes that are stored in a pool. The last component recombines the stored routes using a Set Packing formulation to produce good quality solutions. [START_REF] Cura | An artificial bee colony algorithm approach for the team orienteering problem with time windows[END_REF] solves the TOPTW using an Artificial Bee Colony (ABC) approach that mimics the foraging behavior of honey bees. He introduced a new food source acceptance criterion based on SA and a new scout bee search behavior based on a local search procedure. [START_REF] Guibadj | Memetic algorithm with an efficient split procedure for the team orienteering problem with time windows[END_REF] designed a Memetic Algorithm (MA) that uses the route-first cluster-second approach to solve the TOPTW. Solutions are represented as permutations of a subset of reachable customers called giant tours, and a splitting procedure is applied to extract the optimal set of feasible routes with respect to the order of visits in the giant tour. In [START_REF] Gunawan | Well-tuned ils for extended team orienteering problem with time windows[END_REF], the authors extend the ILS algorithm of Gunawan et al. (2015a) to solve the TOPTW. The algorithm generates an initial solution using a greedy insertion heuristic that chooses the customers to be inserted based on roulette-wheel selection. The initial solution is then improved using an ILS that uses different local search operators and a combination of acceptance criteria and perturbation mechanisms to balance between diversification and intensification in the search process. In a later work, Gunawan et al.

(2015b) presented a hybridization that embeds an ILS into a SA algorithm to address the drawback that is the early termination of the ILS. [START_REF] Schmid | An effective large neighborhood search for the team orienteering problem with time windows[END_REF] proposed an "Effective Large Neighborhood Search" (ELNS) approach to solve the TOPTW. The algorithm starts by generating an initial solution. During each iteration of the algorithm, the solution is destroyed and repaired through operators specifically designed for the TOPTW. To the extent of our knowledge, the ELNS is currently the approach that provides the largest proportion of best-known solutions for the TOPTW benchmark instances.

Finally, for a more detailed overview of the literature about orienteering problems, we invite the reader to refer to Vansteenwegen et al. (2011) and [START_REF] Gunawan | Orienteering problem: A survey of recent variants, solution approaches and applications[END_REF] who provide reviews on several relevant variants of the orienteering problem, including the TOPTW, and discuss applications and solution methods, both exact and heuristic, for OP variants.

Problem definition

The TOPTW is defined on a complete directed graph G = (V, A), where V = {0, 1, 2, ..., n} is the vertex set, and A = {(i, j) : i = j, i, j ∈ V } the arc set. Vertex 0 represents the depot, which corresponds to the starting and ending points, and each vertex i ∈ V \ {0} represents a customer associated with a non-negative profit (score) p i , a non-negative service time σ i , and a predefined time window [e i , l i ]. The profit of each customer i ∈ V \ {0} can be collected at most once by a vehicle within its associated time window, i.e., a vehicle cannot visit the customer i if it arrives later than l i , and in the event that it arrives earlier than e i , it has to wait until e i before the service can start. A time window [e 0 , l 0 ] is associated with the depot, where e 0 = 0 refers to the earliest departure time and l 0 refers to the latest possible arrival time at the depot. A non-negative travel time c i,j is associated with each arc (i, j) ∈ A. Travel times are assumed to satisfy the triangle inequality.

A fleet of m identical vehicles is available at the depot, and each of them performs at most one route.

A route is an ordered subset of customers that are visited by the same vehicle. Routes must start and end at the depot. We note the total profit collected during a route r, P (r) = i=size(r) i=1

p r[i] , and its total duration (travel cost)

C(r) = c 0,r[1] + i=size(r) i=2 (c r[i-1],r[i] + W r[i] + σ r[i]
) + c r[size(r)],0 , where r[i] denotes the i th customer of the route and W u the vehicle's waiting time at customer u. Note that this travel cost corresponds to the arrival time at the depot. The total duration of each route is constrained within a predefined time limit L max . We assume that L max equals the latest possible arrival time at the depot, i.e., L max = l 0 . A route is considered feasible if and only if C(r) ≤ L max and if each customer is visited within its time window. Thus, a feasible TOPTW solution S consists of at most m feasible routes in which each customer is visited at most once. The objective is to find a solution S that maximizes the total collected profit, i.e., that maximizes r∈S P (r). For mixed integer linear programming formulations of TOPTW see [START_REF] Montemanni | Ant colony system for team orienteering problems with time windows[END_REF]Vansteenwegen et al., 2009).

Solution approach

Method overview

In this paper, we propose a randomized Multi-Start Iterated Local Search procedure (MS-ILS) enhanced by an adaptive memory mechanism. Originally, MS procedures are simple memoryless algorithms that sample the solution space by applying a local or neighborhood search from multiple randomly generated initial solutions. If an initial solution falls inside the attraction basin of a global optimum, the local search will pull it to this global optimum. Due to their simplicity, more often than not, MS procedures alone have difficulties to compete with more aggressive metaheuristics and must be strengthened by complementary diversification techniques to help surmount local optima.

Our approach employs an Iterated Local Search (ILS) procedure as the local search step of the MS metaheuristic. The ILS is a simple metaheuristic whose principle is to build a sequence of improved local optima to explore the search space. More specifically, starting from a solution S, at each iteration, the ILS generates a new solution S by perturbation of S, then improves it using a local search to obtain S . If S satisfies an acceptance criterion, it becomes starting solution for the next iteration; otherwise the algorithm goes back to S.

The efficiency of our algorithm stems from two key aspects. The first one is to alternate between two search spaces, each one using a different solution representation. This idea has proven to be very effective for various vehicle routing problems [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF][START_REF] Prins | Tour splitting algorithms for vehicle routing problems[END_REF][START_REF] Mendoza | A multi-space sampling heuristic for the vehicle routing problem with stochastic demands[END_REF][START_REF] Montoya | A multi-space sampling heuristic for the green vehicle routing problem[END_REF], including the Team Orienteering Problem (TOP) [START_REF] Bouly | A memetic algorithm for the team orienteering problem[END_REF][START_REF] Dang | An effective pso-inspired algorithm for the team orienteering problem[END_REF]. In our case, the two solution representations are the route representation and the giant tour representation. The route representation is a genuine TOPTW solution, i.e., a set of feasible routes, one for each vehicle in use. To ensure fast feasibility checks upon insertion of a customer, we record, for each route, additional information on the arrival time, the waiting time and the maximum delay allowed for the service at each customer. On the other hand, the giant tour representation consists of an ordered list T of all the accessible customers in V with no route delimiters. It is a permutation T = (T [1], T [2], ..., T [n]) of n customers that ignores route length constraints and time windows. The interesting part about this indirect solution representation is (1) the fact that one giant tour corresponds to multiple TOPTW solutions, and (2) that it is possible to retrieve the optimal solution with respect to the order of customers in polynomial time using a splitting procedure.

Thus, it is possible to search the space of giant tours instead of the TOPTW solution set, without loss of information. The reverse transformation from the route representation to the giant tour representation is achieved through a simple concatenation procedure.

The second key component is the integration of a long-term memory mechanism in order to improve the performance of the algorithm by learning from local optima found during previous iterations. One important drawback of a pure multi-start procedure is being a memoryless method: each iteration is independent of the previous one and no information about the solutions is passed from one iteration to the other. In our approach, we integrate an adaptive memory mechanism to retain information about interesting customer sequences observed in high-quality solutions. The use of a memory mechanism is inspired from the probabilistic tabu search of [START_REF] Rochat | Probabilistic diversification and intensification in local search for vehicle routing[END_REF]. It consists in storing routes extracted from high-quality solutions inside a pool of elite but diverse routes, and then using them to generate new starting solutions. In our approach, the stored routes are used to construct giant tours using a probabilistic function that is biased towards the selection of elements that appear more frequently in high-quality solutions.

Algorithm 1 describes the general structure of our MS-ILS approach. At the beginning, an iterative Algorithm 1: MS-ILS algorithm for the TOPTW Data: TOPTW instance; M: adaptive memory;

iter ILS : number of iterations of the ILS;

iter max : maximum number of restarts.

Result: S best overall best solution for the TOPTW instance.

1 begin 2 iter ← 0 3 (M, S best ) ← initializeM emory() // see Algorithm 2 4 repeat // Multi-start loop 5 T ← constructGiantT our(M) 6 S ← split(T ) // see Section 4.5 7 for j ← 1 to iter ILS do // ILS loop 8 T ← randomRotation(T ) // Perturbation 9 S ← split(T ) S ← localSearch(S ) // see Section 4.6 updateM emory(M, S ) // see Section 4.3 if f (S ) ≥ f (S) then T ← concat(S ) S ← S if f (S) ≥ f (S best ) then S best ← S if new routes have been added into M then // see Section 4.3 iter ← 0 else iter ← iter + 1
until iter = iter max return S best heuristic (described in Section 4.2) is called to initialize the adaptive memory M and the overall best solution S best . At each iteration of the main loop, the algorithm constructs a giant tour T using the routes stored inside the adaptive memory M and extracts the solution S associated with tour T using the split procedure described in Section 4.5. It then performs an ILS (ILS loop) with T as the starting point. The ILS loop is executed until the maximum number of iterations iter ILS is reached. At each iteration of the ILS loop, T is perturbed by applying a random rotation, and a new giant tour T is derived. The perturbation procedure considers T as a circular array and shifts its elements by a random number of positions. Afterwards, the associated solution S is extracted from T , improved using the local search procedure described in Section 4.6, and the obtained routes are inserted into M if they are of good enough quality; otherwise they are ignored.

If S is better than S, it is concatenated into a new giant tour which replaces T in the next ILS iteration, and S is updated. After the ILS loop, S best is updated and the number of main loop iterations iter is reset to 0. Otherwise, iter is incremented and the main loop restarts. The algorithms stops when it fails to insert new routes in the adaptive memory during iter max consecutive iterations.

Memory initialization

The initial set of routes that compose the adaptive memory is generated using a fast heuristic procedure that relies on a Best Insertion Algorithm (BIA) to build feasible solutions. At each iteration, the BIA evaluates all the feasible insertions of unrouted customers then selects and carries out the best one.

The insertion of a customer u between two successive nodes i and j is feasible if and only if u is visited within its time window, the visits to j and to the subsequent nodes remain within their respective time windows, and the length of the route remains less than L max . Note that i and j can be either customers or the depot. Similar to Vansteenwegen et al. (2009), in order to ensure that the feasibility of an insertion is checked in O(1), we record, for each customer included in the solution, its waiting time W i and the maximum duration by which its service can be delayed without rendering the route infeasible M axShif t i . These two measures are given by the following formulas:

W i = max{0, e i -a i } M axShif t i = min{l i -(a i + W i ), W j + M axShif t j }
where a i is the arrival time at node i, [e i , l i ] the time window associated with node i, and j is the node visited after i in the solution. The insertion of u between the nodes i and j delays the arrival at node j by Shif t u,i = (c i,u +W u +σ u +c u,j -c i,j ). The insertion of u is feasible if a u ≤ l u and Shif t u,i ≤ W j +M axShif t j .

Fig. 2, illustrates a route containing four customers. The maximum length of the route is L max = 40.

For each visited customer, we display information on its time window, arrival time, waiting time, and the maximum duration by which its service can be delayed. Customer u can be inserted between customers 2 and 3 since Shif t u,2 = (c 2,u

+ W u + σ u + c u,3 -c 2,3 ) = 6 ≤ 12 = W 3 + M axShif t 3 .
For each feasible insertion, the following criterion is calculated:

z u,i = Shif t u,i /(p u ) α 0 1 2 3 4 u 0 5 6 2 3 4 3 2 [0,10] [10,25] [5,25] [30,40] [0,40] a 1 = 5 W 1 = 0 MaxShift 1 = 5 a 2 = 10 W 2 = 0 MaxShift 2 = 12 a 3 = 13 W 3 = 0 MaxShift 3 = 12 a 4 = 22 W 4 = 8 MaxShift 4 = 4 a 0 = 36 MaxShift 0 = 4 a u = 14 W u = 1 Shift u,2 = 6 [15,30] σ 1 = 2 σ u = 2 σ 3 = 3 σ 4 = 2 σ 2 = 1 a 0 = 36 MaxShift 0 = 4 Figure 2: Feasibility check.
where α is a control parameter. The best insertion is determined by the couple (u, i) for which the criterion z u,i is minimized. The aim of this criterion is to favor the insertions which offer a good trade-off between the increase of the tour length and the collected profit. The parameter α is used to prioritize profit over the time consumption of an insertion. The value of α is discussed in Section 5.2. After the insertion of u between i and j, the information associated with the successors of u in the route are updated using the following formulas:

W * j = max{0, W j -Shif t u,i } Delay j = max{0, Shif t u,i -W j } M axShif t * j = M axShif t j -Delay j
where Delay j represents the amount of time by which the service of customer j is delayed. When updating the data of a visit after j, Shif t u,i is replaced with the delay of the visit that precedes it.

A pseudo-code of the initialization heuristic is given in Algorithm 2. First, a feasible solution is constructed from scratch using the BIA described above. After that, at each iteration of the heuristic, a sequence of consecutive customers is removed from each route of the solution. The solution is then rebuilt by inserting unrouted customers using the BIA. The sequences to remove are identified by a starting position start and a length of q customers. The length of the removed sequences q is first set to 1, and at the end of each iteration, start is moved by q positions and q is incremented by one. If the solution is improved after repair, q is reset to its initial value. Destruction is applied in a circular manner, i.e., when the end of a route is reached, the removal continues from its beginning. This part of the algorithm is similar to the parallel remove-and-repair operator described in Section 4.6. In order to keep q from becoming too big and causing a large part of the solution to be destroyed at each iteration, it is reset to 1 each time its value reaches half the number of routed customers of the largest route (r max ) in the solution, in terms of the number of visited customers.

The destruction and construction phases are repeated until the algorithm reaches iter init iterations without improving the best solution.

q ← q + 1 19 if q ≥ size(r max )/2 then q ← 1 20 return (M, S best )

Memory update

For performance reasons, we limit the size of the adaptive memory to Size M . When it becomes full, we need to remove some routes to be able to insert new ones. In order to keep promising routes, we apply the following strategy. Before its insertion, each route is labeled with the total collected profit and the total duration of the solution to which it belongs. The adaptive memory is then sorted in decreasing order of profits using the duration to resolve any equality. When the adaptive memory is full, the route to be inserted is first compared to the weakest route of the memory. If the new route's label is weaker, i.e., it has a lower total collected profit, the route is ignored. Otherwise, it is added to the memory and the weakest route is deleted.

The idea is that routes that belong to solutions with higher scores are more likely to contain elements of optimal or sub-optimal solutions. To speed up the initialization of the adaptive memory, the routes of each solution generated by the BIA are considered for inclusion in the memory, so the heuristic is only used once.

Giant tour construction

The construction of giant tours plays an important role in our algorithm since it allows us to use information collected through past iterations for diversification and intensification purposes. To construct a giant tour, we select m routes from the adaptive memory and combine them as follows. Let M be a copy of the current memory. First, we extract a route r from M in a probabilistic way, and discard from M all the routes r that share at least one common customer with r. This process is then repeated until m routes are extracted, or until M becomes empty. When selecting routes, we favor those extracted from solutions with higher quality. To that effect, we use the roulette wheel selection mechanism. After that, the depot is removed from each of the selected routes. The resulting routes are then concatenated into a random order to form a partial giant tour, which we complete by randomly spreading out the unrouted customers over the tour. These unrouted customers are randomly inserted at the beginning of the tour, at the end, or between the routes in such a way that the m selected routes remain untouched, as can be seen in the example depicted in Fig. 3.

Splitting algorithm

In this section, we present the optimal splitting algorithm used to extract a TOPTW solution from a giant tour. Splitting a giant-tour T consists of extracting m distinct sub-sequences from T such that, (1) each extracted sub-sequence corresponds to a feasible TOPTW route, (2) there are no shared customers between the extracted sub-sequences, and (3) the total profit of the extracted sub-sequences is maximized. In the remainder of this section, we will refer to the above sub-sequences of T by extracted routes or simply routes, for convenience.

When splitting a giant-tour T , not all the feasible sub-sequences of T need to be considered to select m routes that maximize the total profit. [START_REF] Bouly | A memetic algorithm for the team orienteering problem[END_REF] introduced the notion of saturated routes and showed that solutions containing only saturated routes are dominant when splitting giant-tours for the TOP.

In the following, we provide an adapted definition of "saturated routes", and show that the same dominance property holds when splitting giant-tours for the TOPTW.

Definition 1. Given a giant tour

T = (T [1], T [2], . . . T [n]
), a "saturated route" in T is a sub-sequence

π i = (T [i], T [i + 1], . . . T [j]), 1 ≤ i ≤ j ≤ n, such that the route r = (0, T [i], T [i + 1], . . . , T [j], 0) is feasible
and either j = n, or the route r = (0,

T [i], T [i + 1], . . . , T [j], T [j + 1], 0) is infeasible.
Proposition 1. For any instance of the TOPTW where m is the number of available vehicles, for any gianttour T of this instance, there exists an optimal solution to the splitting problem of T where all the extracted routes are saturated.

Proof. Let T be a giant-tour visiting all the customers of a TOPTW instance in any order, and let S be an optimal solution obtained when T is split into m routes. S contains at most m extracted routes copy of the adaptive memory

π i k = (T [i k ], T [i k + 1], . . . , T [j k ]), 1 ≤ k ≤ m, 1 ≤ i k ≤ j k ≤ n which, can be re-ordered in such a way that 1 ≤ i ≤ j 1 < i 2 ≤ j 2 < • • • < i m ≤ j m ≤ n.
Let π i k be the k-th route in S. Two cases are identified:

-Case 1 : route π i k is such that j k + 1 < i k+1 (i.e. there are unrouted customers between the end of π i k and the start of π i k+1 ). If π i k were not saturated, then route

π i k = (T [i k ], T [i k + 1], . . . , T [j k ], T [j k + 1])
would be feasible, would have a greater profit, and would not conflict with the other extracted routes.

As such, it would be possible to improve S by replacing π i k with π i k . This contradicts the hypothesis that S is optimal. Hence, π i k is saturated.

-Case 2 : route π i k is such that j k + 1 = i k+1 (i.e. route π i k+1 starts immediately after route π i k ends).

If

π i k is not saturated, then route π i k = (T [i k ], T [i k + 1], . . . , T [j k ], T [i k+1 ]
) would be feasible, and route

π i k+1 = (T [i k+1 + 1], T [i k+1 + 2], . . . , T [j k+1 ]
) would also be feasible, due to the triangle inequalities.

Furthermore, π i k and π i k+1 do not share customers with one another, and the sum of their profits is equal to that of π i k and π i k+1 . Thus, by replacing π i k and π i k+1 with π i k and π i k+1 , respectively, in S, we obtain a new optimal S with one less unsaturated route. Hence, by reiterating the process on every unsaturated route of S, we can produce an optimal solution where all the routes are saturated. [START_REF] Dang | An effective pso-inspired algorithm for the team orienteering problem[END_REF] presented an evaluation procedure that efficiently uses the limited number of saturated routes to reduce the complexity of the splitting process. They reduced the splitting problem to a knapsack problem with conflicts (KPC) [START_REF] Takeo | Heuristic and exact algorithms for the disjunctively constrained knapsack problem[END_REF], and used a dynamic programming algorithm to determine the optimal splitting in polynomial time. In this paper, we extend this splitting procedure to tackle time window constraints.

Extraction of saturated routes

When splitting a giant-tour, we extract one saturated route starting from every customer in the tour. Let Π T = {π 1 , π 2 , ..., π n } denote the set of the n possible saturated routes extracted from the giant tour T . When extracting these routes, it must be ensured that each customer is visited within its time window and that the route is completed within the given time limit L max , i.e., that C(π i ) ≤ L max ∀i ∈ {1, ..., n}. Starting from customer T [i], we initialize the route π i = (0, T [i], 0). We then extend it by consecutively including the following customers T [j], j ≥ i as long as they are visited within their time windows. If the vehicle arrives at a customer T [j] before the beginning of its time window, a waiting time is added. If the vehicle arrives too late at customer T [j], if the inclusion of T [j] violates the time limit constraint, or if the process reaches the end of the giant tour, the extension of π i is stopped and the route is considered saturated. In the split algorithm for the basic TOP, all the saturated routes can be extracted in O(n) because the cost of any sub-sequence (i + 1, ..., j, j + 1) can be deduced in O(1) from that of (i, ..., j). However, this property is not satisfied by TOPTW sub-sequences. In our best implementation, route extraction runs in O(n 2 ).

Fig. 4 shows a giant-tour for a TOPTW instance with 8 customers, and all the saturated routes that can be extracted from it. The maximal length of a route is L max = 60. The giant-tour is represented in Fig. 4a.

It is a permutation of the 8 customers. The figure displays the travel time between every two successive nodes, and gives for each node, its service time σ i , its time widow (between brackets), and its travel time to and from the depot. The saturated routes are presented in Fig. 4b. Route π 1 , for example, corresponds to the route (0, 5, 7, 6, 0) and has a length C(π 1 ) = 5 + 5 + 10 + 10 + 15 + 5 + 10 = 60. Customer 1 cannot be included in this route, since its inclusion makes the total travel time greater than L max . In the case of route π 3 , customer 3 cannot be included in the route, not because of the time limit, but because the vehicle reaches it at instant 35, when its time window has already closed. Route π 8 consists only of customer 8, because it is the last customer of T .

Selection of saturated routes

Once all the saturated routes have been extracted from T , the problem of splitting T into a feasible TOPTW solution becomes that of choosing m routes that do not share customers and that maximize the total collected profit. This can be formulated as a Knapsack Problem with Conflicts (KPCG). In a KPCG, we consider a set of items to be put into a knapsack of limited volume. Each item is associated with a value and a volume, and some items are in conflict with each other. Conflicting items cannot be put in the knapsack together. The objective of the KPCG is to find a subset of non-conflicting items that fits into the knapsack and maximizes the total packed value. In our splitting procedure, the volume of the knapsack is equal to m, the maximal number of available vehicles. Each saturated route π i corresponds to an item i, and two items i and j are in conflict if their corresponding routes π i and π j share customers.

Usually, conflicts between items in a KPCG are modeled using a graph, called a conflict graph. In our particular case, the conflict graph for the splitting problem is an interval graph (Tarjan, 1975). In short, a graph H = (X, U ) is an interval graph if there is a mapping I between the vertex set X and sets of consecutive integers (called intervals) such that two vertices in X are adjacent if and only if their respective corresponding intervals intersect i.e. ∀i ∈ X, ∀j ∈ X, (i, j) ∈ U ⇐⇒ I(i) ∩ I(j) = ∅. To solve the KPCG, we use the algorithm proposed by [START_REF] Sadykov | Bin packing with conflicts: A generic branch-and-price algorithm[END_REF]. They devised a pseudo-polynomial algorithm to solve a KPCG with interval conflict graphs, defined with n items and a knapsack capacity equal to W , in O(n • W ). Hence, the following result for the TOPTW emerges.

Proposition 2. Given a TOPTW instance with m available vehicles and n saturated tours extracted from a giant tour T , the splitting of T can be done optimally in O(m • n) time and space.

Proof. Each saturated route extracted from a giant-tour T can be represented by the set of positions it covers, i.e., each route

π i = (T [i], T [i + 1], . . . , T [j]), 1 ≤ i ≤ j ≤ n is represented by the set I(π i ) = {i, i + 1, . . . , j}
of consecutive integers. Hence, the set of saturated routes can be mapped to the set of vertices of an interval graph H. Moreover, a non-empty intersection between two intervals I(π i ) and I(π i ) indicates the presence of shared customers between the two routes π i and π i , and corresponds to an edge (i, i ) in the graph H.

As outlined above, once the saturated routes are extracted, the problem of splitting a giant-tour becomes that of choosing m saturated routes that do not share customers and such that the total collected profit is maximized. This translates into solving a KPCG where the volume of the knapsack is equal to m, all the items have unitary weights and their values are equal to the profit of the corresponding route. The conflict graph H associated to this KPCG is an interval graph. Given the capacity m of the knapsack and the number of saturated routes n, and based on the work of [START_REF] Sadykov | Bin packing with conflicts: A generic branch-and-price algorithm[END_REF], we conclude that the splitting can be done in O(m • n) time and space.

Figure 5 illustrates the splitting process of the giant-tour presented in Fig. 4. We want to split T into a solution with two routes (m = 2) with L max = 60. Fig. 5b shows the saturated routes that were extracted from T and their collected profit P (π i ). Observe that each saturated route corresponds to a sub-interval of the interval [1,8]. For example, route π 1 corresponds to the interval [1, 3] (as it covers the first three nodes of T ), route π 3 corresponds to the interval [3,4], and route π 8 corresponds to [8,8]. Furthermore, if two routes share at least one customer, then the intersection of their intervals is not empty, and vice versa. Fig. 5c represents the interval graph that corresponds to the conflict graph of the saturated routes. Every vertex of the graph is mapped to the interval of positions of the corresponding saturated route, and two vertices are adjacent if and only if their respective intervals intersect with each other. Finally, the optimal solution obtained after solving the KPCG is shown in Fig. 5d. It is composed of the two routes π 3 and π 6 and its total score is 155.

Local search

The local search (LS) procedure is a key component of our algorithm since it serves as an intensification mechanism to improve the quality of the solutions constructed at each iteration. As previously mentioned, our algorithm uses two different solution representations, each with its own characteristics. Hence, we divide the neighborhood operators that compose the local search procedure into two sets: the first set is applied on giant tours, while the second set is performed on routes.

Giant tour neighborhoods

The LS includes two neighborhoods for giant tours:

• shift operator : moves a randomly chosen customer from its current position in the giant tour to a different one;

• swap operator : chooses a random customer from the giant tour and exchanges its position with that of another randomly chosen customer;

Because the routes that compose a solution are not explicitly defined in the giant tour representation, the shift and swap operators do not need to check if feasibility is maintained, since every move is considered feasible in a giant tour representation. All that is needed to evaluate a move is to extract the new routes by splitting the giant tour. For performance purposes, we use the greedy splitting algorithm described in greater detail below.

Route neighborhoods

For route improvement, the LS procedure uses the classical routing neighborhoods Or-Opt and 2-Opt* to reduce the total travel time of each route, and three other neighborhoods to fill the resulting time saved in order to increase the total profit. The 2-Opt* operator replaces arcs (i, i + 1) from route u and (j, j + 1) from route v with arcs (i, j + 1) and (j, i + 1), thus interchanging the two sub-paths without altering the order of visits. The Or-Opt operator relocates a sequence of one or two successive customers in the same route, without altering their order of visit. The remaining operators follow the same remove-and-repair principle of the initialization heuristic described in Section 4.2 but differ in the way they select the customers to be removed and in how they repair the solution:

• Random remove-and-repair : removes d randomly chosen customers from the solution and, then, using the Best Insertion Algorithm (BIA), inserts currently unrouted customers. At each iteration, d is randomly generated in the interval [1, D max ].

• Sequential remove-and-repair : this operator focuses on improving one single route at a time. This operator removes from each route r a sequence of q consecutive customers starting from position start mod size(r) and uses BIA to insert new customers into the destroyed route. At each iteration, start is moved by q customers, and q is increased by one. When all the customers in a route have been removed at least once, it skips to the next route in the solution. This operator stops if an improvement is found, or once all the routes of the solution have been tested. Fig. 6 shows an iteration of this operator. • Parallel remove-and-repair : Similar to the previous operator, this operator removes from each route r a sequence of q consecutive customers starting from position start mod size(r). However, instead of focusing on one route, it is applied in parallel on all the routes of the solution, and then inserts unrouted customers using the BIA. An example of this move is shown in Fig. 7. The search in this neighborhood stops when an improvement is found, or when each customer in the solution has been removed at least once. Removing parallel sequences from different routes creates free time slots across the whole solution and gives the opportunity to BIA to move customers between routes, and introduce more profitable ones in the solution.

Local search algorithm

The local search procedure works as follows. At each iteration, the LS randomly selects one of the previous neighborhoods and explores it. All the neighborhoods have the same probability of being chosen whether they consider routes or giant tours. The search in a given neighborhood stops as soon as an improvement is found 18 or when no improving move can be found. When a neighborhood operator fails to find an improvement, it is discarded until the solution is improved by another operator. This process is repeated until all neighborhoods fail to find an improvement to the current solution.

Greedy split and concatenation

Because the LS procedure uses two different sets of neighborhoods, and because it can randomly switch from an operator of one set to that of the other, the number of times a solution shifts from one representation to the other can grow large very quickly. This negatively impacts computation times, mainly due to the optimal splitting procedure. To address this aspect of the LS, we devise a fast heuristic splitting method and a giant tour construction process to go with it.

Inside the LS algorithm, the construction of a giant tour T is achieved by concatenation of the routes that compose the solution one after the other, and then by appending the remaining unrouted customers at the end of the tour.

On the same principle as the above fast concatenation method, we describe a greedy splitting heuristic to evaluate a giant tour

T = (T [1], T [2], ..., T [n]). Consider θ k = (T [i], T [i + 1], ..., T [n]) a sub-sequence of T . It is possible to divide θ k into two parts by extracting the saturated route π i k = (T [i k ], T [i k + 1], ..., T [j k ]) and considering the second part θ k+1 = (T [j k + 1], T [j k + 2], ..., T [n]
) as a new sequence of unrouted customers.

Our greedy splitting heuristic consists in performing the previous bisection m times on consecutive sequences θ k , starting with θ 1 = T . Fig. 8 shows an example of splitting using the heuristic procedure. The heuristic splitting is used as a quick way to evaluate moves in the swap and shift neighborhoods. The optimal splitting algorithm is applied at the end of the LS procedure to extract the best routes from the giant tour.

Computational experiments

Our algorithm was coded in C++ using the Standard Template Library (STL) for data structures, and compiled using the GNU GCC compiler in a Linux environment. In the following, we present the computational experiments that were carried out in order to tune the algorithm, to evaluate its performance relative to state-of-the-art algorithms, and to assess the contribution of each of its components. All the experiments were conducted on a single thread on an Intel Xeon X7542 CPU at 2.67GHz processor.

Benchmark instances

Our algorithm was tested on the various TOPTW benchmark instances. All the literature benchmarks are available online at https://www.mech.kuleuven.be/en/cib/op.

TOPTW benchmark instances are derived from Solomon et al.'s (Solomon, 1987) instances for the Vehicle Routing Problem with Time Windows (VRPTW), and Cordeau's [START_REF] Cordeau | A tabu search heuristic for periodic and multi-depot vehicle routing problems[END_REF] instances for the Multi Depot Periodic VRPTW (MDPVRPTW). The instances were adapted by considering the customer demands in the original data sets as node profits in the TOPTW instances. Travel time between two vertices is assumed to be equal to the euclidean distance. It is rounded down to the first decimal for Solomon's instances and to the second decimal for Cordeau's instances. Solomon's instances are organized into six sets: C100, C200, R100, R200, RC100 and RC200, divided according to the distribution of vertices on the plane (clustered, random, random-clustered) and the width of the time windows (narrow, wide). Each instance contains either 50 or 100 customers. Since Cordeau's instances were originally meant for a periodic variant of the VRP, when using them for the TOPTW, we consider that all customers are available on the same day.

Instances in this data set contain a number of customers ranging from 48 to 288. [START_REF] Righini | Decremental state space relaxation strategies and initialization heuristics for solving the orienteering problem with time windows with dynamic programming[END_REF] constructed the first OPTW instances using 29 of Solomon's instances from sets C100, R100 and RC100, and ten instances from Cordeau's (pr1 to pr10). The first benchmark instances for the TOPTW were introduced by [START_REF] Montemanni | Ant colony system for team orienteering problems with time windows[END_REF]. They used the earlier OPTW instances for the case where m = 1 and extended them by increasing m to two, three, and four available vehicles. They also proposed another 27 instances from sets C200, R200, and RC200 from Solomon's instances, and ten others from Cordeau's instances (pr11 to pr20). Furthermore, Vansteenwegen et al. (2009) used the original instances of Solomon and Cordeau to introduce new benchmark instances where the number of vehicles considered for each instance makes it possible to visit all the customers. Since it is possible to visit all customers, the optimal solution for these instances is known: it is equal to the sum of all customer profits.

Like [START_REF] Hu | An iterative three-component heuristic for the team orienteering problem with time windows[END_REF], we refer to these instances as the "OPT" data set. As Vansteenwegen et al. (2009) point out, algorithms to solve the TOPTW are generally not designed expecting the possibility of including all the customers, which makes solving some of the "OPT" instances to the optimum rather difficult.

Parameter tuning

As described in Section 4, our proposed MS-ILS has five sets of parameters:

• α, the control parameter of the BIA;

• D max , the maximum number of customers removed by the random removal operator;

• Size M , the size of the adaptive memory;

• iter init , the number of iterations of the initialization heuristic;

• iter max and iter ILS , the number of iterations of the main algorithm.

In order to calibrate these parameters, we carried out preliminary experiments on a subset of randomly selected instances from the TOPTW instances described in Section 5.1. We started from base parameter values that we identified during the development of our MS-ILS and then tuned each parameter, one after the other, by varying the base value. For each parameter, we kept the best setting found before proceeding to tune the next parameter. Given that parameters iter max , iter ILS , and Size M have a similar effect on the performance of the algorithm, we chose to tune them in a joint fashion. More details, on the tuning experiments can be found in Section A of the supplementary material associated with this paper.

The final calibration results are displayed in

Computational comparisons

In order to evaluate the performance of our proposed MS-ILS for TOPTW, we compared its results with those of the following state-of-the-art algorithms:

• Variable Neighborhood Search (VNS) proposed by Tricoire et al. ( 2010),

• Greedy Randomized Adaptive Search procedure with Evolutionary Local Search algorithm (GRASP-ELS) of [START_REF] Labadie | Hybridized evolutionary local search algorithm for the team orienteering problem with time windows[END_REF],

• LP-Based Granular Variable Neighborhood Search (GVNS) of [START_REF] Labadie | The team orienteering problem with time windows: An lp-based granular variable neighborhood search[END_REF],

• Iterative Three-Component Heuristic (I3CH) of [START_REF] Hu | An iterative three-component heuristic for the team orienteering problem with time windows[END_REF].

• Large Neighborhood Search (ELNS) of [START_REF] Schmid | An effective large neighborhood search for the team orienteering problem with time windows[END_REF].

The results of GVNS, GRASP-ELS, I3CH, ELNS, and our MS-ILS were all obtained with five runs of the algorithm on each instance, while VNS was run ten times on each instance. Note that some results of VNS that were not originally reported by Tricoire et al. (2010) were made available on the authors' website http : //prolog.univie.ac.at/research/op/. Since they are better results than those originally published, we used them instead for our comparisons.

For the sake of fair comparison between the algorithms in terms of computational effort, the reported running times of each algorithm were adjusted to account for the speed difference between the different computation setups. Similar to [START_REF] Hu | An iterative three-component heuristic for the team orienteering problem with time windows[END_REF], we used the Super Pi benchmark for this purpose. Super

Pi is a single-threaded program that computes the digits of π up to a specified number and is commonly used as an estimate of CPU speed. Table 2 indicates the CPU used by each algorithm, its Super Pi score, which corresponds to the number of seconds needed to compute the first one million digits of π, and the associated time scaling factor. To obtain this factor, we estimated the performance of each processor by considering the performance of our machine to be 1. For GRASP-ELS, GVNS and I3CH, we opted to use the Super Pi scores reported by [START_REF] Hu | An iterative three-component heuristic for the team orienteering problem with time windows[END_REF]. Unfortunately, the scores of VNS and ELNS are not available.

Furthermore, we only had limited information on the experimental setup used by each of them to be able to estimate their Super Pi scores. For the VNS, given that the paper was published in 2010, we decided to use the same scaling factor as the one used for the GVNS (0.32), because of the proximity of their publication dates. As for ELNS [START_REF] Schmid | An effective large neighborhood search for the team orienteering problem with time windows[END_REF], the authors did not report comparisons based on CPU times, and we were not able to estimate their Super Pi score by only relying on the family of their CPU and its clock rate. Therefore, we assumed similar performance to our setup. In the remainder of this section, the time values reported in previous works are adjusted by the associated factors, in order to account for CPU differences. For example, the adjusted computational time of a solution obtained by I3CH can be obtained by multiplying its CPU time by 0.95. to the best-known solutions in the literature. These results are presented in tables, of which, each consists of two identically structured parts. Each part contains the name of the instance, the best-known solution (BKS) to the instance, including the ones found by our method and the optimal solutions of Duque et al.

(2015) and [START_REF] Tae | A branch-and-price approach for the team orienteering problem with time windows[END_REF], the maximum score (Z max ) obtained by our algorithm, the relative error (rpe), the average score (Z avg ), the average error (arpe), and the average computational time in seconds (cpu avg ).

Table 3 reports the results obtained by all the state-of-the-art algorithms and the MS-ILS using the previous configuration on the standard benchmark. As can be seen in the table, our MS-ILS achieves an average relative gap lower than those achieved by the other algorithms for every set of instances. As for the relative gap, it achieves better results than the literature on most instances, apart from instances rc200 with m = 3, and rc100 with m = 4 where ELNS does slightly better, but at the cost of higher computational effort.

On Cordeau's instances, our method requires a little more computation time to find high-quality solutions compared to the likes of the GRASP-ELS and the GVNS, but is still much faster than I3CH and the VNS, and the quality of the solutions it obtains is much better than that of the other algorithms. Furthermore, note that Solomon's instances with wider time windows become easier to solve as m becomes larger.

The performance of our algorithm on the "OPT" data set is shown in Tables 4 and5. The column #OP T indicates the number of optimal solutions found by the MS-ILS. Note that, in the case of column #OP T , the value reported at the bottom of the two tables corresponds to the sum of the values of the column. Note also that the authors of GRASP-ELS [START_REF] Labadie | Hybridized evolutionary local search algorithm for the team orienteering problem with time windows[END_REF] and GVNS [START_REF] Labadie | The team orienteering problem with time windows: An lp-based granular variable neighborhood search[END_REF] Tables 6 and7 display the results of two runs of the MS-ILS using two different settings of iter max , iter ils and Size M . The first setting is the same as the one used in the above comparisons with the literature.

The second one is a slower version of the algorithm obtained by setting iter max = 30n/m, iter ils = 10 and Size M = 100m. The main observation here is that, if allowed to run longer, the MS-ILS is able to further improve solution quality. It is able to find the best known solutions for all of Solomon's instances at least one time out of five. On the "OPT" benchmark, it can obtain the optimal solutions for 59 out of 66 instances.

Overall, it is able to find the best known solutions for 93% of the instances available in the literature.

During our experiments, including the tuning experiments, the MS-ILS found 57 new best-known solution the literature finds 65% of the current best-known solutions, and achieves an arpe of 0.80% on the standard benchmark and 0.34% on the "OPT" benchmark, and a rpe of 0.51% and 0.27%, on the two benchmarks, respectively. Finally, our experiments show that the MS-ILS can be tuned to either favor speed or solution quality and remain relatively better than other approaches.

Performance analysis

In the following, we discuss the results of the experiments conducted to evaluate the impact of each of the All of the above configurations were run ten times on each instance of both the standard and the "OPT" benchmark using the parameter values given in Section 5.2. Table 9 shows the average gap relative to the best know solutions, and the average CPU time achieved by each of the configurations from (a) to (d).

Comparisons between the configurations (a) to (c) highlight the impact of alternating between route and giant tour search spaces instead of only using one of them. Disabling either one of the search spaces translates into a decrease of solution quality compared to configuration (a), but removing the route operators has a bigger impact on solution quality than removing giant tour operators. We still decided to keep both spaces in the design of the algorithm because both of them help improve solution quality; besides, the giant tour search spaces is still needed to find good solutions for some instances.

Comparisons between configurations (a) and (d) show the contribution of including the adaptive memory into the multi-start framework. As we can see, adding the adaptive memory significantly improves the average gap on the two benchmarks, with a reasonable impact on CPU time. In terms of solution quality, the use of the adaptive memory has a more significant impact when solving instances of the standard benchmark with m ≥ 2 than when solving instances with m = 1. This is because when m = 1, the algorithm chooses a previous local optimum and tries to improve it, whereas when m ≥ 2, the algorithm constructs new solutions using several from previous local optima.

Most of the time, disabling one component results in a decrease in computation times. However, depending on the instance, the opposite can happen, and disabling one component may result in larger computation algorithms and is able to outperform them in terms of overall solution quality and computation times. In particular, the MS-ILS is able to find the current best-known solutions, or better ones, for 78% of the benchmark instances within reasonable runtimes, and achieves an overall average relative gap of 0.30% and 0.26% on the two benchmarks of the literature, respectively. If given more time, the MS-ILS finds the current best-known solutions for 89% of the literature instances. Furthermore, our approach was also able to find new best solutions for 57 instances for which no optimal solution has yet been found.

Finally, the method proposed herein is flexible in the sense that it can be easily adapted to accommodate new constraints or to address other variants of the problem. One potential direction for future research would be to extend our algorithm to more realistic versions of orienteering problems with time windows.

  Example of a TOP instance. Depot Customer (b) Example of a TOP solution.

Figure 1 :

 1 Figure 1: Team Orienteering Problem.

  who identified some of its drawbacks and provided ways to overcome them by only considering the best solution found during the construction phase, and only applying the local search procedure on solutions on which it has not yet been applied. Vansteenwegen et al. (2009) introduced a fast and simple Iterated Local Search (ILS) for the problem. Their algorithm was made fast and simple by only keeping the insert and shake steps. The authors also introduced a new data set with more difficult instances constructed from the instances of Solomon (1987) and Cordeau et al. (1997). Tricoire et al. (2010) defined the Multi-Period Orienteering Problem with Multiple Time Windows (MuPOPTW) as a means for scheduling sale representatives to visit customers.

  Choose one route and discard all that share at least one customer with it from the adaptive memory. Choose a second route the same way as in Fig.3b. Arrange the chosen routes and complete the giant tour by inserting the remaining customers at random.

Figure 3 :

 3 Figure 3: Illustration of the giant tour construction process.

  Saturated routes extracted from T .

Figure 4 :

 4 Figure 4: Extraction of saturated tours.

  Saturated routes with the corresponding intervals.

  Conflict graph for the KPCG.

Figure 5 :

 5 Figure 5: Example of the splitting procedure for a TOPTW instance with m = 2, n = 8, and Lmax = 60.

Figure 6 :

 6 Figure 6: Sequential remove-and-repair.

Figure 7 :

 7 Figure 7: Parallel remove-and-repair.

Figure 8 :

 8 Figure 8: Heuristic splitting procedure when m = 2.

Figure 9 :Figure 10 :

 910 Figure 9: Comparison of different MS-ILS settings with the literature based on cpu and rpe.

  key components of our algorithm, namely the combination of two different search spaces and the adaptive memory. To do this, we derive alternative versions of the MS-ILS by disabling search spaces one at a time, or by disabling the adaptive memory. The considered configurations are as follows: a. Standard: the standard algorithm described in Section 4. b. No Routes: an alternate version of MS-ILS where the route search space operators are disabled. c. No Tours: an alternate version of MS-ILS where the giant tour search space operators and the splitting procedures are disabled. To construct solutions from the memory, we simply select m routes that do not share customers. d. No Memory: version of the standard algorithm where the adaptive memory and giant tour construction are disabled. Each iteration of the main loop starts from a randomly generated solution.

Table 1 .

 1 The final parameter values were chosen to provide a good trade-off between solution quality and runtime from the authors' perspective.

	Parameter Description

Table 1 :

 1 Final parameter values.

Table 2 :

 2 Estimation of single-thread performance.standard benchmark, and in terms of average relative percentage error (arpe) for the "OPT" benchmark, but we chose to include comparisons on the basis of rpe and arpe for both benchmarks. These two metrics are computed as: rpe = (BKS-Zmax) BKS * 100% and arpe = (BKS-Zavg) BKS * 100%, where Z max denotes the best score obtained over different runs and Z avg the average score. Column cpu avg of each table reports the average computational time of the above mentioned algorithms in seconds. The detailed results obtained by our MS-ILS are presented in the supplementary material associated with this paper, where they are compared

	Algorithm	CPU	Super Pi Estimate Factor
	VNS	2.4 GHz CPU (reference unknown)	Unknown	≤ 0.32
	GRASP-ELS	Intel Pentium 4 processor, 3.00 GHz	44.3	0.32
	GVNS	Intel Pentium (R) IV, 3 GHz CPU	44.3	0.32
	I3CH	Intel Xeon E5430 CPU clocked at 2.66 GHz	14.7	0.95
	ELNS	Intel Xeon 3.1 GHz (reference unknown)	Unknown	≈ 1
	MS-ILS	Intel Xeon X7542 CPU at 2.67GHz	14.1	1

Tables 3 to 5 summarize the results achieved by each algorithm. Comparison of solution quality is usually drawn in terms of relative percentage error (rpe) with respect to the best-known solution (BKS) for the

Table 6 :

 6 Comparison of two different settings of MS-ILS on the standard benchmark.

	Instances	MS-ILS (2,4,80)	MS-ILS (30,10,100)
		rpe%	arpe%	cpu	rpe%	arpe%	cpu
	m=1
	c100	0.00	0.00	2.0	0.00	0.00	32.5
	r100	0.00	0.02	2.1	0.00	0.00	25.9
	rc100	0.00	0.02	1.1	0.00	0.00	20.0
	c200	0.00	0.08	9.9	0.00	0.00	80.4
	r200	0.03	0.10	22.9	0.00	0.03	167.1
	rc200	0.00	0.12	13.4	0.00	0.05	131.8
	pr01-pr10	0.02	0.18	13.0	0.00	0.01	131.2
	pr11-pr20	0.85	1.15	21.5	0.85	0.97	205.0
	m=2
	c100	0.00	0.03	3.6	0.00	0.00	30.2
	r100	0.00	0.04	2.7	0.00	0.01	23.8
	rc100	0.00	0.05	2.9	0.00	0.01	20.8
	c200	0.00	0.14	16.9	0.00	0.07	85.0
	r200	0.08	0.26	30.9	0.01	0.19	187.9
	rc200	0.15	0.39	25.9	0.02	0.22	171.9
	pr01-pr10	0.10	0.43	18.1	0.03	0.20	124.4
	pr11-pr20	0.84	1.29	32.1	0.69	0.94	254.1
	m=3
	c100	0.00	0.09	5.1	0.00	0.04	39.8
	r100	0.00	0.06	3.5	0.00	0.00	27.5
	rc100	0.00	0.25	2.9	0.00	0.09	26.1
	c200	0.09	0.18	5.2	0.03	0.11	32.1
	r200	0.03	0.06	3.2	0.00	0.05	18.7
	rc200	0.07	0.11	8.0	0.00	0.09	36.5
	pr01-pr10	0.34	0.74	25.7	0.12	0.41	168.7
	pr11-pr20	0.48	0.94	43.3	0.12	0.64	349.4
	m=4
	c100	0.10	0.31	7.6	0.10	0.22	44.7
	r100	0.10	0.18	4.4	0.01	0.09	34.9
	rc100	0.11	0.33	3.5	0.01	0.08	31.8
	c200	0.00	0.00	0.1	0.00	0.00	0.1
	r200	0.00	0.00	0.1	0.00	0.00	0.1
	rc200	0.00	0.00	0.1	0.00	0.00	0.1
	pr01-pr10	0.42	0.85	34.4	0.07	0.58	241.2
	pr11-pr20	0.48	1.08	59.1	0.17	0.73	415.7
	Average	0.13	0.30	13.3	0.07	0.18	98.7

Acknowledgments. The authors would like to thank the reviewers for their valuable comments and suggestions which improved the quality of the paper. This work is supported by the French Environment and Energy Management Agency (ADEME), the Hauts-de-France region and the European Regional Development Fund (ERDF). It was carried out within the framework of the Labex MS2T, funded by the French Government (Reference ANR-11-IDEX-0004-02). It is also partially supported by the TCDU project (Collaborative Transportation in Urban Distribution, ANR-14-CE22-0017).

Table 5: Performance comparison based on rpe average for "OPT" data set.

while Table 5 shows the comparisons based on the rpe. Even though the "OPT" data set is known to be difficult to solve, MS-ILS achieves both the smallest average relative gap and the second smallest relative gap, and is able to obtain 59 out of 66 optimal solutions in a reasonable amount of time. The I3CH is the only other algorithm that finds slightly better solutions than the MS-ILS, but it requires significantly higher computational times to do so. In some cases, it is up to twenty times slower than the MS-ILS.

For a more thorough comparison, we conducted additional experiments using other combinations of values for parameters iter max , iter ils and Size M in order to achieve smaller and larger computation times than those presented above. We then compared the performance of every combination with those of the stateof-the-art algorithms, namely VNS, GRASP-ELS, GVNS, I3CH and ELNS in terms of computation time (CPU) and in terms of either relative error(rpe) or average relative error (arpe). The combinations of values used for these experiments were chosen based on the results of the Pareto test described in Section 5.2.

Figures 10 and9 show the performance of the MS-ILS using different combinations of parameters compared to state-of-the-art algorithms; the first in terms of computation time and relative error (rpe), the second in terms of cpu and average relative error (arpe). As can be seen in both figures, the MS-ILS achieves better results compared to the remaining algorithms: not only does it achieve smaller deviations when computation times are equivalent, but each combination of parameter values we tested resulted in smaller deviations Altogether, our MS-ILS is very competitive. It is able to find the current best-known solutions for 78% of the literature instances and was able to improve the best-known solutions for many of them. It achieves an average relative error (arpe) of only 0.30% on the standard benchmark and 0.26% on the "OPT" benchmark, and a relative percentage error (rpe) of 0.13% and 0.20%, on the two benchmarks, respectively, while still being faster than most state-of-the-art algorithms. In comparison, the previous best performing approach in Table 9: Performance analysis of MS-ILS components.

times. For example, in the standard benchmark on Solomon's instances with m = 1; disabling the route space local search operators hinders the progress of the algorithm towards good solutions, which explains the increase in CPU time.

Conclusion

In this paper, we introduced a very effective Multi-Start Iterated Local Search (MS-ILS) for the Team Orienteering Problem with Time Windows. Our algorithm is based on a local search that alternates between two different search spaces: the route search space that corresponds to actual solutions to the TOPTW, and the giant tour search space that makes it possible to explore the solution space without being limited by time windows and length constraints. In order to improve solution quality, several local search operators are included to be used in each search space. Additionally, the algorithm integrates an adaptive memory mechanism to further improve performance by making use of previous local optima to build better solutions.

Computational results have shown that our approach performs very well compared to state-of-the-art