Pengfei He
email: pengfeihe606@gmail.com

Jin-Kao Hao
email: jin-kao.hao@univ-angers.fr

Qinghua Wu
email: qinghuawu1005@gmail.com

Grouping memetic search for the colored traveling salesmen problem

Keywords: Memetic algorithm, Local search, Colored traveling salesmen problem, TSP and Multiple TSP

The colored traveling salesmen problem is a node routing problem with multiple salesmen, where the cities are divided into m exclusive city sets and one shared city set. The objective is to minimize the total traveling distance of m Hamiltonian circuits (routes) under the following constraints: each exclusive city is to be visited by the corresponding salesman, while each shared city can be visited by any salesman. In this work, we present the rst grouping memetic algorithm for solving this challenging problem. The algorithm includes three main components: (i) a greedy randomized heuristic for population initialization; (ii) a dedicated local search procedure for local optima exploration; (iii) a backbone-based crossover operator for solution recombination. We show computational results on three sets of 65 popular benchmark instances to demonstrate the competitiveness of our algorithm. We especially report improved upper bounds for 38 instances (for more than 58% cases). We also present rst computational results with the general CPLEX solver, including 10 proven optimal solutions. Finally, we shed lights on the impacts of the key components of the algorithm. We make the code of the algorithm publicly available.

1 Introduction

The colored traveling salesmen problem (CTSP) can be stated as follows [START_REF] Li | Colored traveling salesman problem[END_REF]. [START_REF] Arnold | Knowledge-guided local search for the vehicle routing problem[END_REF] Let G=(V, A) be a complete undirected graph, where V = {0, 1, 2, • • • , n} is 3 the set of nodes (or cities) and A = {{i, j} : i, j ∈ V, i = j} is the set of edges. Each edge {i, j} ∈ A has a non-negative weight c ij representing the traveling distance between cities i and j. All cities are divided into m + 1 disjoint sets: m exclusive city sets CTSP generalizes a variant of the classical traveling salesman problem (TSP), known as the multiple traveling salesmen problem (MTSP) where all cities are shared [START_REF] Bektas | The multiple traveling salesman problem: an overview of formulations and solution procedures[END_REF][START_REF] Carter | A new approach to solving the multiple traveling salesperson problem using genetic algorithms[END_REF]. Besides, if there is only one salesman (m = 1), CTSP becomes TSP [START_REF] Applegate | The Traveling Salesman Problem: a Computational Study[END_REF]. Given that CTSP generalizes these NP-hard problems, solving CTSP is computationally challenging.

CTSP is a useful model for a number of practical problems [START_REF] Li | Colored traveling salesman problem[END_REF]. For instance, Li et al. [START_REF] Li | Collision-free scheduling of multi-bridge machining systems: a colored traveling salesman problem-based approach[END_REF] presented typical applications concerning multi-bridge machining systems, i.e., a dual-bridge waterjet machining center and a dualmanipulator hull welding system. In these systems, there are two independent individual machines (cutting machines and manipulators) with their individual workspaces (exclusive cities) and a shared workspace (shared cities).

Finding a collision-free scheduling of these machines corresponds to solving a CTSP with two salesmen. The rice harvesters problem studied in He et al. [START_REF] He | Fields distinguished by edges and middles visited by heterogeneous vehicles to minimize non-working distances[END_REF][START_REF] He | Optimisation of the harvesting time of rice in moist and non-moist dispersed elds[END_REF] can be considered from the CTSP perspective. Rices from moist elds and non-moist elds need to be harvested by harvesters (salesmen) such that moist elds (exclusive cities) can only be visited by crawler-harvesters, while non-moist elds (shared cities) can be visited by both crawler-harvesters and wheelers-harvesters. Scheduling the harvesters comes down to solving the CTSP problem.

Heuristics and metaheuristics are natural approaches for nding sub-optimal solutions of dicult optimization problems that cannot be solved exactly. Unlike the related MTSP for which numerous heuristics are available [START_REF] Soylu | A general variable neighborhood search heuristic for multiple traveling salesmen problem[END_REF][START_REF] Yuan | A new crossover approach for solving the multiple travelling salesmen problem using genetic algorithms[END_REF][START_REF] Carter | A new approach to solving the multiple traveling salesperson problem using genetic algorithms[END_REF][START_REF] Singh | A new grouping genetic algorithm approach to the multiple traveling salesperson problem[END_REF][START_REF] Wang | Memetic algorithm based on sequential variable neighborhood descent for the minmax multiple traveling salesman problem[END_REF],

only six metaheuristic-based algorithms were studied for CTSP: genetic algorithm [START_REF] Li | Colored traveling salesman problem[END_REF], variable neighborhood search [START_REF] Meng | Variable neighborhood search for a colored traveling salesman problem[END_REF], articial bee colony (ABC) [START_REF] Pandiri | A swarm intelligence approach for the colored traveling salesman problem[END_REF][START_REF] Dong | Articial bee colony algorithm with generating neighbourhood solution for large scale coloured traveling salesman problem[END_REF], ant colony optimization [START_REF] Dong | Ant colony optimisation for coloured travelling salesman problem by multi-task learning[END_REF], and iterated two phase local search [START_REF] He | Iterated two-phase local search for the colored traveling salesmen problem[END_REF].

These algorithms have reported valuable computational results on various benchmark instances. Meanwhile, it is observed that they lack robustness and stability in particular when they are applied to solve large scale instances.

In this work, we investigate for the rst time the powerful memetic algo-rithm (MA) framework for solving CTSP and present a competitive grouping memetic algorithm (GMA) dedicated to the problem. Indeed, eective MAs have been proposed to solve the related MTSP [START_REF] Wang | Memetic algorithm based on sequential variable neighborhood descent for the minmax multiple traveling salesman problem[END_REF][START_REF] Kashan | Grouping evolution strategies: An eective approach for grouping problems[END_REF][START_REF] Kóczy | Enhanced discrete bacterial memetic evolutionary algorithm-an ecacious metaheuristic for the traveling salesman optimization[END_REF][START_REF] Lu | Hybrid evolutionary search for the traveling repairman problem with prots[END_REF] and several vehicle routing problems [START_REF] Nagata | Edge assembly-based memetic algorithm for the capacitated vehicle routing problem[END_REF][START_REF] Vidal | A hybrid genetic algorithm for multidepot and periodic vehicle routing problems[END_REF]5,[START_REF] Prins | A simple and eective evolutionary algorithm for the vehicle routing problem[END_REF]. However, most of these MAs are based on the so-called giant tours and split algorithms, which are not suitable for CTSP due to the presence of exclusive cities. We consider CTSP from the perspective of grouping problems [START_REF] Falkenauer | Genetic Algorithms and Grouping Problems[END_REF] and introduce an eective grouping memetic algorithm (GMA). The proposed algorithm integrates two complementary key components: an original local optima exploration procedure (to nd high quality local optima, Section 3.3) and a dedicated backbone-based crossover operator (to generate promising new ospring, Section 3.4). As demonstrated by the computational results shown in Section 4, the proposed algorithm competes very favorably with the state-of-the-art CTSP algorithms on three sets of benchmark instances.

The rest of this paper is organized as follows. Section 2 presents a literature review and related works. The proposed grouping memetic algorithm is presented in Section 3. Computational results and comparisons with state-ofthe-art algorithms are presented in Section 4. In Section 5, the impacts of key components of the algorithm are discussed. Section 6 presents conclusions and future research directions.

2 Literature review and related works Given the theoretical and practical signicance, CTSP has attracted considerable attention in recent years and several heuristic methods have been presented. In this section, we review the existing solution approaches for CTSP and related works.

In 2014, Li et al. [START_REF] Li | Colored traveling salesman problem[END_REF] introduced the colored traveling salesmen problem to optimize routes of a dual-bridge waterjet cutting machine tool. As solution methods, they presented four genetic algorithms (basic GA, GA with greedy initialization, hill-climbing GA and simulated annealing GA), where the dualchromosome encoding was used to represent the candidate solutions. The rst chromosome is a permutation of all cities except depot 0, while the second chromosome assigns a salesman to each of the shared and exclusive cities in the corresponding position of the rst chromosome. They presented computational results on 20 small scale benchmarks created from existing symmetric TSP instances (with up to 101 vertices). They showed that the hybrid algorithm combining simulated annealing and GA dominated the three other algorithms and their algorithms performed better than the general mixed integer programming tool Lingo.

Then, in 2017, Meng et al. [START_REF] Meng | Variable neighborhood search for a colored traveling salesman problem[END_REF] proposed a variable neighborhood search (VNS) which employs a direct-route encoding to represent the solutions. VNS consists of two phases. The rst phase perturbs the current solution by two shaking operations (Interchange and Relocation), while the second phase improves the perturbed solution by applying a local search based on two search operations (neighborhood change and 2-opt). Compared with the four GAs [START_REF] Li | Colored traveling salesman problem[END_REF], VNS showed its competitiveness on the 20 benchmark instances.

Later, in 2018, Pandiri and Singh [START_REF] Pandiri | A swarm intelligence approach for the colored traveling salesman problem[END_REF] presented an articial bee colony (ABC)

based on the m-tour encoding. This encoding uses m arrays, and each array includes all the cities visited by the corresponding salesman. They provided a proof that the size of the solution space of CTSP with the m-tour encoding is smaller than that of the dual-chromosome encoding. They showed that ABC could match or update the best results reported in [START_REF] Li | Colored traveling salesman problem[END_REF][START_REF] Meng | Variable neighborhood search for a colored traveling salesman problem[END_REF] on the 20 small scale benchmark instances with very short cuto times. Besides, they introduced 8 new medium scale instances (with 229-666 cities) and reported computational results.

Also in 2018, Dong et al. [START_REF] Dong | Ant colony optimisation for coloured travelling salesman problem by multi-task learning[END_REF] employed an ant colony optimization (ACO) with multi-tasks learning for solving CTSP. The multi-task cooperative learning was proposed to improve the eciency of ACO. To assess their algorithm, they introduced 6 medium (with 202-431 cities) and 5 large instances (with 1002 cities) and showed the competitiveness of ACO compared with the four GAs [START_REF] Li | Colored traveling salesman problem[END_REF]. Nevertheless, this algorithm did not compete well with VNS [START_REF] Meng | Variable neighborhood search for a colored traveling salesman problem[END_REF] on the 20 small scale instances in terms of the best and average results.

In 2019, Dong et al. [START_REF] Dong | Articial bee colony algorithm with generating neighbourhood solution for large scale coloured traveling salesman problem[END_REF] presented another articial bee colony algorithm (ABC)

and reported computational results on 26 new large instances (with 2461-7397 cities). These new large scale instances could be used by subsequent studies to evaluate their algorithms. However, this ABC algorithm performed worse than the ABC algorithm of [START_REF] Pandiri | A swarm intelligence approach for the colored traveling salesman problem[END_REF] on the 20 small scale instances.

Finally, He and Hao [START_REF] He | Iterated two-phase local search for the colored traveling salesmen problem[END_REF] proposed an iterated two-phase local search (ITPLS), which is based on a new adjacency representation of the candidate solutions.

This representation relies on an array A[m, n + 1] such that for each route r (r = 1, . . . , m), A[r, i] = j (i, j = 0, . . . , n, i = j) if and only if the route goes from city i to city j. ITPLS applies jointly inter-route optimization and intra-route optimization for solution improvement, reinforced by a probabilistic greedy perturbation strategy to diversify the search. Extensive computational results were reported on all the benchmark instances available in the literature (a total of 65 instances), including 29 improved best known results.

According to the computational results reported in the literature, we identify ABC [START_REF] Pandiri | A swarm intelligence approach for the colored traveling salesman problem[END_REF] and ITPLS [START_REF] He | Iterated two-phase local search for the colored traveling salesmen problem[END_REF] as the current state-of-the-art CTSP algorithms.

CTSP generalizes the popular multiple traveling salesmen problem (MTSP), which has attracted much interest in the last decades. For instance, Wang et al. [START_REF] Wang | Memetic algorithm based on sequential variable neighborhood descent for the minmax multiple traveling salesman problem[END_REF] introduced a memetic algorithm for solving MTSP, which includes a variable neighborhood descent to search local optima. Another evolutionary algorithm was proposed by Kashan et al. [START_REF] Kashan | Grouping evolution strategies: An eective approach for grouping problems[END_REF] for solving MTSP from the perspective of grouping problems. Other representative studies were reported in [START_REF] Soylu | A general variable neighborhood search heuristic for multiple traveling salesmen problem[END_REF][START_REF] Yuan | A new crossover approach for solving the multiple travelling salesmen problem using genetic algorithms[END_REF][START_REF] Carter | A new approach to solving the multiple traveling salesperson problem using genetic algorithms[END_REF][START_REF] Singh | A new grouping genetic algorithm approach to the multiple traveling salesperson problem[END_REF]. However, these methods are not suitable for CTSP, because of the presence of exclusive cities.

In this work, we are interested in designing a practically eective algorithm for solving CTSP with the memetic framework. This is motivated by two considerations. First, one notices that the route of each salesman can be considered as a TSP solution. Therefore, the optimization of each individual route can naturally benet from existing powerful TSP methods. Second, we can consider CTSP from the perspective of grouping problems in the sense that the shared cities are to be dispatched into m groups (m being the number of salesmen). As such, the population-based memetic framework with a meaningful crossover represents an attractive approach given that it has been applied with great success to several dicult grouping problems (e.g., [START_REF] Falkenauer | Genetic Algorithms and Grouping Problems[END_REF][START_REF] Galinier | An ecient memetic algorithm for the graph partitioning problem[END_REF][START_REF] Zhou | Memetic search for identifying critical nodes in sparse graphs[END_REF]).

3 Grouping memetic algorithm for CTSP Given a CTSP instance, the search space explored by CTSP is a multi-route problem whose candidate solutions consist of m tours where the k-th tour includes city 0, the exclusive cities of C k and some shared cities of S.

In this section, we present the grouping memetic algorithm (GMA) for solving CTSP. For a CTSP instance, GMA explores a search space Ω composed of all candidate feasible solutions, where a candidate solution ϕ consists of m tours {r 1 , r 2 , . . . , r m } with r k (k = 1, 2, . . . , m) being the k-th route visited by the k-th salesman. Given a solution ϕ ∈ Ω, its objective value f (ϕ) is given by the total distance of its m routes. The goal of GMA is thus to nd a solution with the smallest objective value as possible.

In the literature, three common methods were used to represent solutions of CTSP: dual chromosome encoding [START_REF] Li | Colored traveling salesman problem[END_REF], m-tour encoding [START_REF] Pandiri | A swarm intelligence approach for the colored traveling salesman problem[END_REF] and adjacency representation encoding [START_REF] He | Iterated two-phase local search for the colored traveling salesmen problem[END_REF]. In this work, we adopt the adjacency representation encoding, which has the advantage of encoding each route (group of cities) independently to facilitate inter-routes operations. The interested reader is referred to [START_REF] He | Iterated two-phase local search for the colored traveling salesmen problem[END_REF] for more details and an illustrative example.

General scheme

The proposed GMA algorithm consists of four main components: population initialization, local optima exploration, backbone-based crossover and population updating. GMA starts with an initial population P of p solutions generated by the population initialization procedure (Section 3.2). It then repeats a number of generations during which new candidate solutions are sampled.

At each generation, the backbone-based crossover combines two randomly and uniformly selected parent solutions to generate a promising ospring solution (Section 3.4). The local optima exploration (LOE) is then applied to improve the ospring solution (Section 3.3), followed by population update (Section 3.5). This evolution process is terminated when a predened stopping condition (e.g., an allowed number of generations, an allotted cuto time limit) is reached. In this work, we use a cuto time limit. The pseudo-code of GMA is shown in Algorithm 1.

Algorithm to the following steps: 1) build a partial route for each of the m salesmen by using the corresponding exclusive cities; 2) dispatch the shared cities among the m partial routes to obtain a complete solution. The pseudo-code of the greedy randomized heuristic is shown in Algorithm 2. During the rst step (lines 4-12), to create the k-th partial route r k , one rst initiates the route with the city 0. Then, the exclusive cities in C k are selected randomly and uniformly, and inserted one by one into r k such that the insertion gives the smallest increase of the route distance. When all exclusive cities of every salesman are inserted into the corresponding route, the rst step stops, leading to a partial solution ϕ composed of m partial routes. During the second step (lines 13-18), the shared cities j from S \ {0} are processed randomly and uniformly, and inserted one by one into a route of the partial solution ϕ such that its total distance increase is minimal. When all shared cities are inserted, an initial solution is obtained. The rst step has a time complexity of

O(|C m | 2 × m),
while the second step is bounded by O(|S|×n). Therefore, the time complexity of the greedy randomized heuristic is O(|S| × n).

Local optima exploration

Local optimization plays a key role in a memetic algorithm and constitutes one of the driving forces for nding solutions of increasing quality. For an effective examination of local optima, GMA employes a specic strategy that combines an inter-route optimization and an intra-route optimization procedure heuristic. Specically, our local optima exploration procedure (LOE) iterates two complementary search components: the constrained cross-exchange operator (CCE) (Section 3.

Constrained cross-exchange

The conventional cross-exchange was initially designed for vehicle routing problems [START_REF] Arnold | Knowledge-guided local search for the vehicle routing problem[END_REF][START_REF] Taillard | A tabu search heuristic for the vehicle routing problem with soft time windows[END_REF][START_REF] Chen | Iterated variable neighborhood descent algorithm for the capacitated vehicle routing problem[END_REF]. It is a generic local search operator which performs exchanges of two consecutive substrings (sub-routes) ri and rj from two dierent routes r i and r j . However, given the particularity of exclusive cities in CTSP, the cross-exchange cannot be used directly in our context. For this reason, we propose a constrained cross-exchange (CCE) in this work. Moreover, it is known that the cross-exchange has a high time complexity [START_REF] Arnold | Knowledge-guided local search for the vehicle routing problem[END_REF][START_REF] Taillard | A tabu search heuristic for the vehicle routing problem with soft time windows[END_REF]. CCE uses a suitable pruning technique to reduce this complexity.

The evaluation of a CCE move for CTSP is summarized in two steps. The rst step is to determine the start of two substrings and the second step is to identify the suitable length of both substrings ( r k 1 and r k 2 ). For the start of substring r k 1 , we rst need to nd an edge which will break route r k 1 . Suppose the edge is {I 1 , I 2 }. Then, a suitable new neighbor of city I 1 needs to be determined. To limit the number of candidate moves, CCE uses the following heuristic pruning technique that only considers the neighbors among the N n nearest cities. Suppose that city J 3 is such a neighbor, and city J 3 belongs to route r k 2 . If edge {I 1 , J 3 } is added as a new edge, edge {J 2 , J 3 } or edge {J 3 , J 4 } should be removed. Once the starts of two substrings (I 2 and J 3 ) are determined, we need to identify the length of each substring. It is worth noting that each substring should not include any exclusive cities because these cities are only visited by the corresponding salesman.

Because the number of cities of each substring can vary from 0 to α (a parameter), all feasible combinations of the two substrings with their given starts can be listed, and the move gain δ for each combination can also be calculated.

There are at most (α + 1) 2 combinations of two substrings. When a substring is empty and the other is non-empty ( r k 1 = ∅ or r k 2 = ∅), these two cases are Or-opt [START_REF] Or | Traveling salesman-type combinatorial problems and their relation to the logistics of blood banking[END_REF][START_REF] Taillard | A tabu search heuristic for the vehicle routing problem with soft time windows[END_REF]. However, both substrings cannot be empty simultaneously. Therefore, at most (α + 1) 2 -1 combinations of two substrings could be listed for two given starts. Then, we need to identify the best move (i.e., with the largest gain δ l ) in these combinations. So far, a CCE move < r k 1 , r k 2 > is acquired and the lengths of two substrings are determined. For all combinations of the two starts, the global minimal move gain δ b can be identied. If δ b < 0, For example, Fig. 1 illustrates two cases of determining the starts of two substrings. Then two complete CCE moves ( r k 1 = {I 2 } and r k 2 = {J 3 , J 4 } or r k 2 = {J 3 , J 2 }) are illustrated in Fig. 2, where cities {I 2 , J 2 , J 3 } are shared. If edge {J 2 , J 3 } is broken in the rst step, the substring r k 2 = {J 3 , J 4 } is serial and in order. However, if edge {J 3 , J 4 } is broken in the rst step, the substring r k 2 = {J 3 , J 2 } is serial and reverse. One may note the following dierences between CCE and cross-exchange [START_REF] Arnold | Knowledge-guided local search for the vehicle routing problem[END_REF].

F lag ← true, RT [k 1 ] ← 1 and RT [k 2 ] ← 1; then,
First, the cross-exchange operator used in [START_REF] Arnold | Knowledge-guided local search for the vehicle routing problem[END_REF] does not limit the length of the substrings to be exchanged; however, in CCE, the length of the substrings must be less than or equal to the value xed by the parameter α. Second, in CCE, exclusive cities are constrained to stay in a route and cannot be moved to other routes. Therefore, the two substrings to be exchanged should not include any exclusive cities. Finally, unlike vehicle routing for which cross-exchange was designed, there is no capacity limitation for each salesman in CTSP. So CCE does not consider this capacity constraint.

Edge assembly crossover (EAX) for TSP

For the optimization of each individual route, the constraint of exclusive cities can be ignored. Thus optimizing each route corresponds to solving a TSP.

There are several sophisticated and powerful heuristics designed for solving TSP. For example, the well-known fast 2-opt heuristic or LK algorithm can be used to improve each route [START_REF] Helsgaun | An eective implementation of the linkernighan traveling salesman heuristic[END_REF][START_REF] Arnold | Knowledge-guided local search for the vehicle routing problem[END_REF][START_REF] Soylu | A general variable neighborhood search heuristic for multiple traveling salesmen problem[END_REF]. In this work, we adopt the EAX heuristic [START_REF] Nagata | A powerful genetic algorithm using edge assembly crossover for the traveling salesman problem[END_REF] 1 , which is among the best TSP heuristics. In our case, EAX helps to keep each route to being optimal or near-optimal in the iterative process of LOE. 1 The code of EAX is available at: https://github.com/sugia/GA-for-TSP 10

Backbone-based crossover

Crossover is another important ingredient of a memetic algorithm and should be designed with care in order to favor transmissions of useful information from parents to ospring [START_REF] Hao | Memetic algorithms in discrete optimization[END_REF], while respecting the problem specic structure. One popular way of designing meaningful crossover for grouping problems such as CTSP is to explore the so-called backbone information, which typically corresponds to solution attributes shared by elite solutions [START_REF] Sun | Memetic search for the equitable coloring problem[END_REF][START_REF] Galinier | Hybrid evolutionary algorithms for graph coloring[END_REF][START_REF] Galinier | An ecient memetic algorithm for the graph partitioning problem[END_REF][START_REF] Zhou | Memetic search for identifying critical nodes in sparse graphs[END_REF]. In this work, we follow this idea and design a dedicated backbone crossover for CTSP.

Let ϕ F and ϕ M be two parent solutions in the population. Based on ϕ F and ϕ M , we divide the set of shared cities except the depot (S \ {0}) into two categories, i.e., common elements and non-common elements.

(i.e., i ∈ r F k ∩ r M k ). If i appears in r F k and r M l (k = l), city i is a non-common element.
Then, an ospring solution ϕ O is constructed in two steps. In the rst step, a donor parent is rst chosen randomly and uniformly between ϕ F and ϕ M .

A partial ospring solution ϕ O is then created by inheriting all m routes of the donor parent without the shared cities. In the second step, for each city

C 1 = {1, 2, 3}, C 2 = {4,

Experimental results and comparisons

This section presents a performance assessment of the GMA algorithm. We show computational studies on well-known benchmark instances from the literature, and comparisons with existing state-of-the-art algorithms for CTSP.

Benchmark and experimental protocol

We employ three sets of 65 benchmark instances, which were commonly used in previous studies on CTSP. The rst set (Set I) contains 20 small instances which were introduced in [START_REF] Li | Colored traveling salesman problem[END_REF], and the number of cities is between 21 to 101 while the number of salesmen m is between 2 and 7. The second set (Set II), introduced in [START_REF] Pandiri | A swarm intelligence approach for the colored traveling salesman problem[END_REF][START_REF] Dong | Ant colony optimisation for coloured travelling salesman problem by multi-task learning[END_REF], contains 14 medium size instances with 202, 229, 431, 666 cities, and 10 -40 salesmen. The last set (Set III), presented in [START_REF] Dong | Articial bee colony algorithm with generating neighbourhood solution for large scale coloured traveling salesman problem[END_REF][START_REF] Dong | Ant colony optimisation for coloured travelling salesman problem by multi-task learning[END_REF], contains 31 large instances with 1002 -7397 cities and 3 -60 salesmen.

GMA was coded in C++ and complied with a g++ compiler with the -O3 option 2 . All experiments were conducted on a computer with an AMD-6134 processor (2.3 GHz and 2 GB RAM) running Linux.

To assess the performance of GMA, we show comparisons with the following algorithms: articial bee colony (ABC) [START_REF] Pandiri | A swarm intelligence approach for the colored traveling salesman problem[END_REF] (2018), ant colony optimization (ACO) [START_REF] Dong | Ant colony optimisation for coloured travelling salesman problem by multi-task learning[END_REF] (2018) and iterated two phase local search (ITPLS) [START_REF] He | Iterated two-phase local search for the colored traveling salesmen problem[END_REF] (2021). Indeed, computational results reported in the literature indicate that these three algorithms represent the state-of-the-art of solving the above benchmark instances, while ABC [START_REF] Pandiri | A swarm intelligence approach for the colored traveling salesman problem[END_REF] and ITPLS [START_REF] He | Iterated two-phase local search for the colored traveling salesmen problem[END_REF] are clearly two dominating algorithms.

So we use ABC (source code unavailable) and ITPLS (source code available)

as the main reference algorithms and cite ACO (source code unavailable) when it is appropriate.

To make the comparisons as fair as possible, we faithfully re-implemented the best ABC algorithm of [START_REF] Pandiri | A swarm intelligence approach for the colored traveling salesman problem[END_REF] 3 . We veried that our implementation (denoted as re-ABC) was able to reproduce the results reported in [START_REF] Pandiri | A swarm intelligence approach for the colored traveling salesman problem[END_REF] (and in fact, our ABC implementation even obtained some better results than those reported in [START_REF] Pandiri | A swarm intelligence approach for the colored traveling salesman problem[END_REF]). To ensure a fair comparison, we ran the compared algorithms on our computer under the same cuto limits: 1, 10 and 60 minutes for Sets I, II and III, with the exception of 240 minutes for the largest instances with at least 7000 cities of Set III. Given the stochastic nature of the compared algorithms, we ran each algorithm 20 times independently to solve each instance with the above time limits.

In order to assess the gaps between the heuristic solutions (from GMA and the reference algorithms) and the optimal solutions, we also investigated the general mixed integer programming solver CPLEX (version 12.7) based on the binary linear programming model from [START_REF] Li | Colored traveling salesman problem[END_REF] (see Appendix A). Our experiment indicated that CPLEX with this model can only solve optimally 10 smallest instances of Set I within 7200 seconds, but fails to solve any instance of Sets II and III due to memory overow.

Parameter tuning

GMA requires 3 parameters: population size p, number of nearest cities N n and parameter α. In order to identify a set of suitable parameters, we used the popular 'IRACE' package [START_REF] López-Ibáñez | The irace package: Iterated racing for automatic algorithm conguration[END_REF] for automatic parameters tuning. The tuning was performed on 8 benchmark instances with 202 to 5397 cities. For the experiment, the tuning budget was set to 500 runs, each with a time limit of half of the cuto time. The studied and nal values (suggested by IRACE) of these parameters are shown in Table 1. Finally, to assess the statistically signicant dierence between GMA and each main compared algorithm, Table 5 shows the p-values from the Wilcoxon signed-rank test. With a condence level of 95%, a p-value smaller than 0.05 indicates a statistically signicant dierence between the pair of compared results.

From Table 2 on the 20 small instances of Set I, the following observations can be made. First, CPLEX is able to solve optimally the 10 smallest instances with 21 -51 cities and 2 -4 salesmen. For the remaining instances, the gap between UP and LP remains reasonable and tends to increase with the size of the instance. For the three heuristic algorithms, they perform equally well in terms of solution quality by reaching their best solutions consistently including the 10 optimal values. The geometric means indicates that the three heuristic algorithms can reach the same results in terms of both the best and average results. Meanwhile, the heuristic algorithms have smaller geometric means compared with CPLEX and thus perform better for this set of instances.

In terms of computational eciency, GMA and re-ABC perform better than ITPLS since they require signicantly less computation times to reach the same results. It is worth mentioning that none of the other algorithms in the literature, such as GA [START_REF] Li | Colored traveling salesman problem[END_REF] (2014), VNS [START_REF] Meng | Variable neighborhood search for a colored traveling salesman problem[END_REF] 4 (2017), ACO [START_REF] Dong | Ant colony optimisation for coloured travelling salesman problem by multi-task learning[END_REF] (2018), and ABC [START_REF] Dong | Articial bee colony algorithm with generating neighbourhood solution for large scale coloured traveling salesman problem[END_REF] (2019) can reach such a performance (they report worse results for some instances or their best results cannot be reached consistently). Table 3 indicates that GMA nds better results for 7 out of the 14 instances, and matches the best results of the reference algorithms for 3 other instances.

The Wilcoxon signed-rank test on the f best and f avg values in Table 5 also conrms that GMA signicantly outperforms the two main reference algorithms.

We do not insist on computation time because the main compared algorithms report solutions of dierent quality. Nevertheless, the three main compared algorithms (re-ABC, ITPLS and GMA) require comparable computation times to reach their best solutions. Note that the results of ACO [START_REF] Dong | Ant colony optimisation for coloured travelling salesman problem by multi-task learning[END_REF] are somewhat inconsistent. Among the six instances tested by ACO, even if it reports three better results than the other algorithms (indicated with a star), its results for the three other instances as well as for most of the 20 small instances of Set I are considerably worse than algorithms, such as ABC [START_REF] Pandiri | A swarm intelligence approach for the colored traveling salesman problem[END_REF] and ITPLS [START_REF] He | Iterated two-phase local search for the colored traveling salesmen problem[END_REF].

Table 4 presents the computational results of the compared algorithms for the re-ABC [START_REF] He | Iterated two-phase local search for the colored traveling salesmen problem[END_REF] ITPLS [START_REF] He | Iterated two-phase local search for the colored traveling salesmen problem[END_REF] GMA (this work) 3. These results clearly show the dominance of the proposed GMA algorithm over the reference algorithms for these large instances, by systematically reporting better results in terms of the best and the average objective values. Moreover, GMA requires the shortest computation times to reach its solutions for this set of large instances, demonstrating its remarkable search capacity and high computational eciency. According to the p-values (less than 0.05) from the Wilcoxon signed-rank test shown in Table 5, the dierence between GMA and each compared algorithm is statistically signicant.

Tables 2 For a more intuitive illustration of the performance assessment of the algorithms, we use the performance prole [START_REF] Dolan | Benchmarking optimization software with performance proles[END_REF], which is a popular benchmarking tool for rigorous comparison of dierent algorithms. In general, performance proles adopt a specic performance metric (in our case, we use f best and f avg ) on all sets of instances. To compare a set of algorithms X over a set of problems Q, the performance ratio is dened by r x,q = fx,q min{f x,q:x∈X,q∈Q } . If an algorithm

x does not report result for a problem q, r x,q = +∞. The performance function of an algorithm x is computed by Q

x (τ ) = |q∈Q|rx,q≤τ | |Q| . The value Q x (τ )
computes the fraction of problems that the algorithm x can solve with at most τ many times the cost of the best algorithm. For example, Q x (1) equals the number of problems that algorithm x solved better than, or as good as the other algorithms in Q. Similarly, the value Q x (r f ) is the maximum number of problems that algorithm x solved. Therefore, Q x (1) and Q x (r f ) represent the eciency and robustness of algorithm x, respectively. Fig. 4 shows the performance proles of GMA, ITPLS and re-ABC. We observe that GMA dominates the reference algorithms in terms of the best and average values. Indeed, GMA has a much higher Q x (1) value compared with the reference algorithms, indicating that GMA can nd better or equal results for all instances. Furthermore, GMA reaches Q x (r f ) rstly, indicating a high robustness of our algorithm. Finally, Table 5 summarizes the results reported by the compared algorithms on the three sets of 65 instances. Column 2 gives the set name and the number of instances in the set. Column 3 shows the quality indicators (f best and f avg ).

Columns 4-6 count the number of instances for which GMA achieves a better, equal or worse value compared with each reference algorithm. The last column presents the p-values from the Wilcoxon signed-rank test. Table 5 reveals large performance gaps between GMA and each reference algorithm on Sets II and III. We conclude that GMA is very competitive for solving CTSP and this is particularly true for large instances.

5 Discussion and analysis

Benet of the key components

In this section, we justify the design choices behind the proposed GMA algorithm. For this, we investigate the impacts of its key components: Constrained Cross-exchange, EAX as well as backbone-based crossover. For our experiments, we used the 45 instances of Sets II and III and ignored the instances of Set I. Indeed, for the instances of Set I, their best-known results can be consistently reached by the state-of-the-art algorithms including ABC, ITPLS and GMA within a very short time. As such, these instances are too easy to be used to compare algorithm variants. shares the other ingredients of GMA.

From the results in Tables 6 and7, we observe that GMA signicantly outperforms GMA 1 on all instances, except gr202-25. For gr202-25, the best result of GMA 1 is slightly better than GMA. Furthermore, the small p-value (less than 0.05) in Table 8 from the Wilcoxon singed-rank test attests the signicant dierence between GMA and GMA 1 . Moreover, the performance proles of To study the eectiveness of the backbone-based crossover, we compared GMA with a third variant GMA 2 . In GMA 2 , the backbone-based crossover is re-placed by a crossover proposed by Singh and Baghel [START_REF] Singh | A new grouping genetic algorithm approach to the multiple traveling salesperson problem[END_REF], which was designed for the related MTSP problem. This crossover selects one of the two parents uniformly at random and copies, from the parent to the ospring, the most promising route (i.e., the route having the smallest ratio of route length to the number of cities in that route). Then all the cities belonging to the route are deleted from both parents by connecting the predecessor of each city to its successor, and the length of the route is updated accordingly.

From the comparative results (f best and f avg ) of GMA and GMA 2 in Tables 6 and7, we observe that in terms of f best , GMA dominates GMA 2 by acquiring 25 better results, 10 equal results and 10 worse results. The Wilcoxon signedrank test, shown in Table 8, also conrms that GMA outperforms signicantly GMA 2 on the large instances (set III). This experiment demonstrates that the backbone-based crossover operator contributes positively to the performance of GMA, in particular for solving large instances.

Finally, Fig. 5 provides other useful information for the importance of each ingredient of GMA. For example, GMA 1 performs the worst because it has the worst (smallest) Q x (1) value and reaches Q x (r f ) lastly. Therefore, we can summarize that EAX is the most important component of GMA, followed by CCE, nally the backbone-based crossover operator.

Inuences of selection, pool updating and mutation

In addition to the local optimization and crossover components, the performance of a memetic algorithm such as GMA could be inuenced by other factors such as parent selection, pool updating and mutation. According to our experiments with the roulette-wheel selection strategy and the rank-pool updating strategy [START_REF] Zhou | Memetic search for identifying critical nodes in sparse graphs[END_REF], no signicant changes were observed regarding the performance of the GMA algorithm. In this subsection, we focus on studying the inuence of mutation. Specically, when the best solution ϕ * is not improved for maxN oImpor consecutive iterations (we empirically set maxN oImpor = 50), the search is judged to be stagnating. Then a mutation operator is triggered to modify one third of the solutions in the population (i.e., each solution is mutated with an equal probability of 1/3). The mutation consists of displacing a certain number of randomly and uniformly chosen shared cities. To be specic, for a solution to be mutated, |S| × 0.3 shared cities are rst removed, leading to a partial solution. Then these removed shared cities are inserted into the partial solution one by one, using the second step of the greedy randomized heuristic to minimize the distance increase. After that, each mutated solution is optimized by the local optima exploration procedure of Section 3.3. Comparative results of GMA and the GMA variant extended with the mutation (called GMA 3 ) are shown in Tables 9 and10. One notices rst that the curves of the population-based GMA and GMA 3

do not start at time 0. This is because that these algorithms spent a nonnegligible portion of the time on generating the initial population (around 60 and 100 seconds for gr431-25 and gr666-30, respectively). From Fig. 6, one observes that re-ABC and ITPLS improve their solution quality quickly at the beginning of the search, and slow down or even stagnate as the time going.

For GMA and GMA 3 , the population initialization step allowed them to start the search with high-quality solutions. The best solution in the population is continually updated when the time goes on, implying that GMA and GMA 3 can better benet from the allowed time to improve their solutions.

Conclusions

The colored traveling salesmen problem is a relevant variant of the popular traveling salesmen problem and generalizes the well-known multiple traveling salesmen problem. In this work, we presented the rst grouping memetic al- indicate that our algorithm is very competitive compared with existing leading algorithms. In particular, it reports 38 new upper bounds while matching 24 best-known results. We also investigated the interest of CPLEX for solving CTSP and reported 10 proven optimal solutions for the rst time.

For future work, there are several perspectives. First, the cross-exchange operator used for inter-route optimization has a high time complexity. This implies that the local optimization component of the proposed algorithm is time consuming. As such, for a given time unit (e.g., a short cuto time), the algorithm will not be able to sample many candidate solutions, limiting thus its performance. To cope with this problem, one possible way is to reduce the number of candidate solutions considered by CCE. To this end, it is interesting to investigate the idea of neighborhood pruning that proves to be successful for vehicle routing [START_REF] Arnold | Knowledge-guided local search for the vehicle routing problem[END_REF][START_REF] Toth | The granular tabu search and its application to the vehiclerouting problem[END_REF] and TSP [START_REF] Helsgaun | An eective implementation of the linkernighan traveling salesman heuristic[END_REF]. Second, recent research on using learning technique to reinforce optimization algorithm showed interesting results (e.g. [START_REF] Zhou | Improving probability learning based local search for graph coloring[END_REF]). As such, it would be useful to study in depth hybrid approaches that combine learning strategies and GMA. Third, CTSP is strongly related to MTSP and TSP, for which powerful algorithms exist. Ideas of these algorithms could be useful for solving CTSP. Finally, to the best of our knowledge, no dedicated exact algorithm exists for CTSP. Eorts are needed to ll the gap.

  {C 1 , C 2 , • • • , C m }, and one shared city set S such that ∪ m k=1 C k ∪ S = V and ∩ m k=1 C i ∩ S = ∅. The cities of an exclusive set C k (k = 1, 2, • • • , m) are to be visited by salesman k only, while the shared cities can be visited by any of the m salesmen. Besides, city 0 (the depot) belongs to the shared city set S and is visited by all salesmen. CTSP is to determine m shortest Hamiltonian tours (routes) starting from the depot and ending at the depot such that each exclusive city in C k is visited exactly once by salesman k and each shared city is visited exactly once by one of the m salesmen. A mathematical model of CTSP is provided in Appendix A.

  3.1) and a TSP heuristic called Edge Assembly Crossover (EAX) [26] (Section 3.3.2). CCE improves solutions by exchanging two substrings (sub-routes) from two routes. The routes modied by CCE are indicated by a binary vector RT of length m (RT [i] = 1 if route i is changed by CCE, RT [i] = 0 otherwise). Then each modied route is further optimized by EAX. CCE and EAX are repeated until the current solution ϕ cannot be further improved. Algorithm 3 shows the pseudo-code of the local optima exploration procedure integrating the CCE operator and the EAX heuristic. Algorithm 3: Pseudo-code of local optima exploration Input: Solution ϕ, number of the nearest cities N n , parameter α Output: Improved solution ϕ B begin ϕ B ← ϕ; F lag ← true; RT [k] ← false ∀k ∈ {1, . . . , m}; /* RT is a binary vector, indicating the routes modified by CCE */ while F lag do < ϕ, F lag, RT >← CCE(ϕ, N n , α); /* Solution improvement by CCE, Section 3.3.1 */ for k = 1, . . . , m do if RT [k] = true then ϕ ← ϕ \ {r k }; 10 r k ← EAX(r k ); /* Intra-route improvement by EAX, Section 3

  solution ϕ is updated by swapping two substrings ( r k 1 and r k 2 ); otherwise, solution ϕ, F lag and matrix RT are returned, because the stopping condition (δ b ≥ 0) of CCE is met. As for the time complexity of CCE, there are O(|S| × (α + 1)) ways to select the rst substring in any given route, while O(N n × (α + 1)) ways exist to select the second substring in another route. Therefore, the time complexity of CCE is O(|S| × N n × ((α + 1) 2 -1)).

1

 1

4 JFig. 1 . 3 IFig. 2 .

 4132 Fig. 1. Illustrative example of starts for a CCE move.

iFig. 3 .

 3 Fig. 3. Illustrative example of the backbone-based crossover Fig. 3 shows an example of the crossover operator for a CTSP instance with 11 cities {0,1,. . . ,10} and m = 2 salesmen with their sets of exclusive cites

  Thus a negativeImp(%) value indicates that GMA improved the best results of the reference algorithms. For Set II, we ignored the geometric means given that the lower bounds needed for their calculations are unavailable. In fact, we tried to obtain LB for these instances by solving, with CPLEX, the linear programming relaxation of the model presented in the Appendix. But CPLEX terminates abnormally due to memory overow without proving any results or bounds.

  of Set III (1002 -7397 cities and 3 -60 salesmen) with the same information as in Table

Fig. 4 .

 4 Fig. 4. Performance proles of GMA and two reference algorithms on the 65 instances of sets I, II, and III. The left part corresponds to the best results while the right part is for the average results.

Fig. 5 Fig. 5 .

 55 Fig.5indicate that GMA surpasses GMA 1 in terms of f best and f avg . Indeed, GMA arrives at Q x (r f ) rstly, much earlier than GMA 1 . These observations illustrate the benet of EAX in GMA.

Fig. 6 .

 6 Fig.6. Convergence charts (running proles) of re-ABC, ITPLS, GMA and GMA 3 for solving two representative instances (gr431-25 and gr666-30). The results were obtained from 20 independent executions of each compared algorithms

  Pseudo-code of the greedy randomized heuristic Input: Instance I (exclusive city sets {C 1 , C 2 , • • • , C m }, shared city set

		Algorithm 2: S) and distance matrix
		Output: A feasible solution ϕ
		begin
		ϕ ← ∅; /* First step: build m partial routes with
		exclusive cities	*/
		for k = 1 to m do
		r k ← {0}; /* Initiate the route with the city 0 1: Pseudo-code of GMA algorithm Input: Instance I, population size p, number of the nearest cities N n , while C k = ∅ do	*/
		parameter α
		Output: The best solution ϕ * found
		begin
		P = {ϕ 1 , ϕ 2 , • • • , ϕ p } ← PopInitilize (I, p); /* Build an initial
		population of p elite solutions, Section 3.2	*/
		ϕ * ← arg min {f (ϕ i ) : i = 1, 2, • • • , p};
		while Stopping condition is not met do
		randomly and uniformly select two parents ϕ F and ϕ M from P ;
		ϕ O ← Backbone_Crossover(ϕ F , ϕ M ); /* Generate an
		offspring solution by backbone-based crossover,
		Section 3.4	*/
		ϕ O ← LOE (ϕ O , N n , α); /* Improve the new solution by
		local optima exploration, Section 3.3	*/
		if f (ϕ O ) < f (ϕ * ) then
		ϕ * ← ϕ O ;
	10	end
	11	P ← P oolU pdating(P, ϕ O ); /* Update the population with
		the new solution, Section 3.5	*/
	12	end
	13	return ϕ * ;
	14 end
		3.2 Population initialization
		The GMA starts its search with an initial population P of p high-quality (elite)
		solutions. To construct a population, we use a greedy randomized heuristic to
		generate a feasible solution, which is further improved by LOE described in
		Section 3.3. The improved solution is then inserted into P if the solution is
		dierent from any existing solution in P ; otherwise, this solution is discarded.
		This process is repeated until p dierent solutions are generated. Thanks to
		the greedy randomized heuristic and subsequent LOE improvement step, we
		obtain a diverse and high-quality population.
		A feasible solution is constructed by the greedy randomized heuristic according

6 Select randomly and uniformly a city i from set C k ; 7 Insert city i in route r k such that the route distance increase is minimal; 8 Remove city i from set C k ; end 10 ϕ ← ϕ ∪ {r k }; 11 end /* Second step: dispatch the shared cities S \ {0} among m partial routes */ 12 S ← S \ {0}; 13 while S = ∅ do 14 Select randomly and uniformly a city j from set S ; 15 Insert city j into a route of ϕ such that the total distance increase is minimal; 16 Remove city j from set S ; 17 end 18 return ϕ; 19 end

Denition 1 :

 1 Given two parent solutions ϕ F = {r F \ {0} with respect to ϕ F and ϕ M is a common element if there exists a k ∈ {1, . . . , m} such that i appears in both r F k and r M k

	1 , r F 2 , . . . , r F m } and ϕ M =
	{r M 1 , r M 2 , . . . , r M

m }, a city i ∈ S

  5, 6}, and the set of shared cities S \{0} = {7, 8, 9, 10} (marked in red color). Let ϕ F = {r F 1 = {0, 1, 2, 3, 8, 9}; r F 2 = {0, 4, 5, 6, 7, 10}}and ϕ M = {r M 1 = {0,2, 1, 3, 9, 10}; r M 2 = {0, 6, 4, 5, 7, 8}} be the parent solutions. By Denition 1, cities 7 and 9 are common elements, while 8 and 9 are non-common elements. First, suppose that ϕ M is the donor parent. Then ospring ϕ O inherits the routes r M 1 and r M 1 by deleting the four shared cities, leading to ϕ O ← {r O 1 = {0, 2, 1, 3}; r O 2 = {0, 6, 4, 5}}. Then the shared cities {7, 8, 9, 10} are successively considered until they are all inserted. City 7 is a common element of the second routes of the parent solutions, it is thus greedily inserted into the partial route r O 2 , supposing this is the cheapest insertion that increases the least the total distance. City 8 is a non-common element, it is greedily inserted into the partial route r O 1 or r O 2 with equal probability. Suppose that route r O 1 is selected, and city 8 is inserted into route r O 1 at the cheapest place leading to the smallest increase of the route distance. Cites 9 and 10 are processed in the same way. When all shared cities {7, 8, 9, 10} are inserted into ϕ O , a feasible ospring solution is constructed successfully, which is then submitted to LOE for further improvement.The time complexity of the crossover can be estimated as follows. The rst step needs to scan all the cities of the donor parent to allow its m routes to be partially inherited. This is achieved in O(n) time. In the second step, the shared cities in S \ {0} are inserted into the partial ospring at the most suitable places, while the time complexity of evaluating each move gain is

O

[START_REF] Applegate | The Traveling Salesman Problem: a Computational Study[END_REF]

. The second step can be performed in O(|S| × n) time. As the result, the time complexity of the crossover is O(|S| × n).

3.5 Pool updating strategy

For each new ospring solution ϕ O improved by LOE in Section 3.3, the pool updating strategy uses ϕ O to update the population P as follows. If the ospring ϕ O is dierent from any existing solutions and better than the worst solution in P , ϕ O replaces the worst solution; otherwise ϕ O is discarded.

Table 1

 1

	Parameters tuning results
	Parameters	Section	Description	Considered values	Final value
	p	3.1	population size	{10,15,20,25,30}	20
	Nn	3.3.1	number of nearest cities	{30,40,50,60,70,80,90}	50
	α	3.3.1	maximum length of substring	{1,2,3,4,5,6,7}	7
	4.3 Computational results and comparisons with existing algorithms
	Computational results of GMA and the reference algorithms on set I are shown
	in

Table 2

 2

	h i=1	U P i LB i )

. For CPLEX, we report for each instance the upper bound (UB), the lower bound (LB) and the Gap given by (U B -LB)/LB ×100. So Gap = 0 implies that an optimal solution is found. Columns 6 -17 report respectively the results of re-ABC, ITPLS and GMA in terms of the best objective value f best (over 20 runs), the average objective value f avg , standard deviation σ and the average time in seconds to reach the best objective value (Time(s)).

For the f best and f avg indicators, equally best values are shown in italic font.

Given that both the upper bounds and lower bounds are available for these instances, we include the geometric mean of each algorithm for a global assessment (row Geomean). For CPLEX, the geometric mean is calculated with the gaps between UB and LB by ( 1 h where U P i and LB i are the upper and lower bound of the ith instance, respectively. Similarly, for the other algorithms (re-ABC, ITPLS, and GMA), we calculate their geometric means for the best and average objective values by replacing U P i with the f best and f avg values, respectively.

Table 3

 3 presents the results of the compared algorithms (re-ABC, ITPLS and GMA) on the 14 medium instances of Set II with 202 -666 cities and 10 -40 salesmen. In addition to the main reference algorithms re-ABC and ITPLS, we also include in this comparison ACO[START_REF] Dong | Ant colony optimisation for coloured travelling salesman problem by multi-task learning[END_REF] for indicative purposes, which only reported results for six instances. For each algorithm except ACO, we show the best and average objective values (f best and f avg ), the standard deviation (σ) and the average time to reach the best objective value (Time(s)). Equally best values are indicated in italic font, while strictly best values are highlighted in boldface. Moreover, the last column Imp(%) provides the percentage improvement of GMA's best result f best over the best objective value f bk of the reference algorithms, computed as (f best -f bk )/f bk × 100.

Table 2 .

 2 Comparative results of GMA and reference algorithms on Set I. The equally best values are indicated in italic.

		Time(s)	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.40	1.00	1.00	1.60	1.00	1.00	1.00	1.00	1.70	1.80	-
	GMA (this work)	f best favg σ	144.92 144.92 0.0	157.48 157.48 0.0	259.36 259.36 0.0	295.31 295.31 0.0	315.97 315.97 0.0	346.24 346.24 0.0	367.84 367.84 0.0	392.14 392.14 0.0	478.08 478.08 0.0	469.50 469.50 0.0	489.99 489.99 0.0	525.98 525.98 0.0	593.28 593.28 0.0	603.79 603.79 0.0	651.99 651.99 0.0	672.73 672.73 0.0	726.82 726.82 0.0	779.15 779.15 0.0	759.55 759.55 0.0	798.85 798.85 0.0	1.0235 1.0235 -
		Time(s)	18.32	13.15	12.70	12.85	16.90	14.45	22.05	11.55	21.55	20.40	28.55	14.35	16.40	17.75	6.75	31.40	18.45	10.80	12.75	12.56	-
	ITPLS	f best favg σ	144.92 144.92 0.0	157.48 157.48 0.0	259.36 259.36 0.0	295.31 295.31 0.0	315.97 315.97 0.0	346.24 346.24 0.0	367.84 367.84 0.0	392.14 392.14 0.0	478.08 478.08 0.0	469.50 469.50 0.0	489.99 489.99 0.0	525.98 525.98 0.0	593.28 593.28 0.0	603.79 603.79 0.0	651.99 651.99 0.0	672.73 672.73 0.0	726.82 726.82 0.0	779.15 779.15 0.0	759.55 759.55 0.0	798.85 798.85 0.0	1.0235 1.0235 -
		Time(s)	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.05	1.00	1.00	1.10	1.00	1.50	1.00	2.10	1.30	1.05	1.25	1.20	-
	re-ABC	f best favg σ	144.92 144.92 0.00	157.48 157.48 0.00	259.36 259.36 0.00	295.31 295.31 0.00	315.97 315.97 0.00	346.24 346.24 0.00	367.84 367.84 0.00	392.14 392.14 0.00	478.08 478.08 0.00	469.50 469.50 0.00	489.99 489.99 0.00	525.98 525.98 0.00	593.28 593.28 0.00	603.79 603.79 0.00	651.99 651.99 0.00	672.73 672.73 0.00	726.82 726.82 0.00	779.15 779.15 0.00	759.55 759.55 0.00	798.85 798.85 0.00	1.0235 1.0235 -
		Gap(%)	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.85	4.39	2.19	3.09	5.85	10.12	7.04	13.76	10.84	15.16	-
	CPLEX	UB LB t(s)	144.92 144.92 1	157.48 157.48 1	259.36 259.36 2	295.31 295.31 20	315.97 315.97 61	346.24 346.24 7	367.84 367.84 46	392.14 392.14 120	478.08 478.08 126	469.50 469.50 773	489.99 485.88 7201	525.98 503.84 7212	596.18 583.41 7211	603.79 585.69 7202	656.56 620.25 7206	687.43 624.25 7202	746.93 697.83 7204	854.23 750.92 7203	783.08 706.47 7209	840.60 729.92 7201	1.0335 --
		Instance	eil21-2	eil21-3	eil31-2	eil31-3	eil31-4	eil41-2	eil41-3	eil41-4	eil51-2	eil51-3	eil51-4	eil51-5	eil76-3	eil76-4	eil76-5	eil76-6	eil101-4	eil101-5	eil101-6	eil101-7	Geomean

Table 3 .

 3 Comparative results of GMA and reference algorithms on Set II. Equally best values are indicated in italic. The strictly best values are indicated in boldface.

		Imp(%)	39.00	74.00	97.00	0.00	0.00	0.00	0.01	-0.23	-0.21	-0.21	-0.52	-0.40	-0.39	-0.26
		Time(s)	396.65	211.80	93.85	201.90	154.85	83.20	73.10	29.50	394.95	159.45	473.55	506.20	512.25	535.70
	GMA (this work)	f best favg σ	99871.00 100162.50 185.46	173477.00 173594.65 75.72	233871.00 234003.35 72.81	222167.00 222173.75 30.19	264146.00 264146.00 0.00	319669.00 319880.15 547.77	406701.00 407389.75 279.37	248447.00 248447.00 0.00	347335.00 347559.80 420.13	415314.00 415387.45 88.31	387562.00 389594.80 3417.34	446475.00 447123.60 328.49	519121.00 519773.45 397.47	650116.00 650974.90 417.87
		Time(s)	93.55	141.25	156.85	226.70	67.00	128.20	186.20	221.45	253.70	296.90	485.90	223.60	249.50	255.05
	ITPLS[14]	favg σ	100009.50 112.58	173523.80 46.77	233857.80 73.17	222347.65 103.50	264146.00 0.00	319671.90 12.97	406884.00 225.72	250036.95 613.23	349238.10 417.38	417963.75 958.14	396841.55 2716.00	449635.25 800.17	522650.90 1006.57	653318.10 927.19
		f best	99871.00	173418.00	233749.00	222167.00	264146.00	319669.00	406664.00	249421.00	348181.00	416552.00	389583.00	448257.00	521149.00	651801.00
		Time(s)	329.45	348.85	316.15	244.60	69.60	303.05	301.35	300.25	333.30	355.20	159.45	248.15	177.55	224.80
	re-ABC [14]	ACO [8] f best favg σ	71924.00* 99871.00 100033.20 110.54	99606.00* 173547.00 173596.80 54.01	118495.00* 233749.00 233817.85 70.16	-222167.00 222354.85 164.08	-264146.00 264146.00 0.00	-319669.00 319669.00 0.00	-406664.00 407194.85 375.21	330554.00 249031.00 249682.25 293.07	464298.00 348056.00 348431.10 203.82	483977.00 416189.00 416758.40 249.58	-390188.00 392234.00 971.38	-448604.00 449997.35 716.97	-522157.00 523583.15 937.90	-652587.00 654001.50 633.57
		Instance	gr202-12	gr202-25	gr202-35	gr229-10	gr229-15	gr229-20	gr229-30	gr431-12	gr431-25	gr431-40	gr666-10	gr666-15	gr666-20	gr666-30

Table 4 .

 4 Comparative results of GMA and reference algorithms on Set III. The strictly best values are indicated in boldface.

Table 5

 5 Summary of comparative results between GMA and two reference algorithms

	Algorithm pair	Set/Instance	Indicator	Better	Equal	Worse	p -value
	GMA vs. ITPLS	I/20	f best	0	20	0	0.00E+00
			favg	0	20	0	0.00E+00
		II/14	f best	7	4	3	2.44E-04
			favg	8	1	5	4.80E-02
		III/31	f best	31	0	0	1.17E-06
			favg	31	0	0	1.17E-06
	GMA vs. re-ABC	I/20	f best	0	20	0	0.00E+00
			favg	0	20	0	0.00E+00
		II/14	f best	8	4	2	1.37E-02
			favg	9	1	4	4.79E-02
		III/31	f best	31	0	0	1.17E-06
			favg	31	0	0	1.17E-06

Table 6

 6 Comparative results on Set II between GMA and its three variants. Strictly best values are shown in boldface.

			GMA		GMA 0		GMA 1		GMA 2
	Instance	f best	favg	f best	favg	f best	favg	f best	favg
	gr202-12	99871.00	100162.50	100292.00	100722.85	100196.00	100573.25

Table 7

 7 Comparative results on Set III between GMA and its three variants. Strictly best values are indicated in boldface. In other words, only EAX is employed in GMA 0 in the local optima exploration component.Computational results of GMA and GMA 0 are shown in Tables6 and 7and summarized in Table8and Fig.5. The results indicate that GMA performs signicantly better than GMA 0 in terms of f best and f avg . For f best , GMA dominates GMA 0 by getting 42 better results out of the 45 tested instances and reporting only one worse result. Furthermore, the statistically signicant dierence between GMA and GMA 0 is veried by the Wilcoxon singed-rank test with a 95% level of condence in Table8. Therefore, this experiment conrms the usefulness of CCE for the GMA algorithm.5.1.2 Benet of EAXTo assess the benet of EAX in LOE, we created another variant GMA 1 in which EAX is replaced by 2-opt[START_REF] Helsgaun | An eective implementation of the linkernighan traveling salesman heuristic[END_REF] for individual route optimization. GMA 1

		GMA		GMA 0		GMA 1	GMA 2
	Instance	f best	favg	f best	favg	f best	favg	f best	favg
	pr1002-5	313885.00	314083.20	314065.00	314495.30	324126.00	327673.20	313867.00	313946.15
	pr1002-10	379846.00

Table 8

 8 Summary of comparative results between GMA and and its three variants

	Algorithm pair	Set/Instance	Indicator	Better	Equal	Worse	p -value
	GMA vs. GMA 0	II/14	f best	11	2	1	9.77E-04
			favg	14	0	0	1.22E-04
		III/31	f best	31	0	0	1.17E-06
			favg	31	0	0	1.17E-06
	GMA vs. GMA 1	II/14	f best	12	1	1	7.32E-04
			favg	14	0	0	1.22E-04
		III/31	f best	31	0	0	1.17E-06
			favg	31	0	0	1.17E-06
	GMA vs. GMA 2	II/14	f best	6	4	4	3.34E-01
			favg	9	1	4	9.42E-02
		III/31	f best	19	6	6	1.60E-03
			favg	23	0	8	9.94E-04

Table 9

 9 Comparative results on Set II between GMA and GMA 3 (with mutation). Strictly best values are indicated in boldface.

			GMA			GMA 3
	Instance	f best	favg	σ	f best	favg	σ
	gr202-12	99871.00	100162.50	185.46	99871.00	100136.95	201.03

Table 10

 10 Comparative results on Set III between GMA and GMA 3 (with mutation). Strictly best values are indicated in boldface.

			GMA			GMA 3
	Instance	f best	favg	σ	f best	favg	σ

The code of our algorithm will be available at http://www.info.univangers.fr/pub/hao/CTSP.html

[START_REF] Bektas | The multiple traveling salesman problem: an overview of formulations and solution procedures[END_REF] Our implementation of ABC[START_REF] Pandiri | A swarm intelligence approach for the colored traveling salesman problem[END_REF] is available from the link given in footnote 2.

VNS reports a wrong result of 465.28 for eil51-2 because this result is smaller than the proven optimal value of 478.08 from CPLEX.

Acknowledgments

We are grateful to the reviewers for their useful comments and suggestions which helped us to signicantly improve the paper. We thank Prof. Jun Li

The results indicate that in terms of f best , GMA 3 outperforms GMA by obtaining 13 better results, 28 equal results and 4 worse results. However, the Wilcoxon signed-rank test shows that there is no statistically signicant dierence. On the contrary, in terms of f avg , GMA 3 dominates GMA by acquiring 34 better results, 2 equal results and 9 worse results. The Wilcoxon signed-rank test conrms that GMA 3 signicantly outperforms GMA. This experiment indicates that the mutation strategy can indeed improve the performance of GMA. Especially, it signicantly reinforces the stability of the algorithm.

Convergence analysis

Finally, we investigate the convergence behaviors of GMA (and the GMA 3 variant with mutation) and two key reference algorithms (re-ABC and ITPLS).

For this study, we acquired the running proles of these algorithms on two representative instances of Set II (gr431-25, gr666-30). We ran each algorithm 20 times with the cuto time of 600 seconds per run and recorded the best objective values during the process. The results of this experiment are shown in Fig. 6.

and Prof. Alok Singh for providing their test instances. Support from the China Scholarship Council (CSC, No. 201906850087) for the rst author is also acknowledged.

A Appendix

In this appendix, we show the mathematical model that we used to report the results of the general ILP solver CPLEX in Section 4. The model is based on [START_REF] He | Iterated two-phase local search for the colored traveling salesmen problem[END_REF] and [START_REF] Li | Colored traveling salesman problem[END_REF].

The binary variable x ijk = 1 indicates that the k-th salesman passes through edge {i, j}, and otherwise x ijk = 0. u ik is the number of cities visited on the This ILP model has two variables for every edge (for both directions) and uses the Miller-Tucker-Zemlin subtour elimination constraints. Given that this model includes more variables than necessary and the Miller-Tucker-Zemlin subtour elimination constraints are known to be very slow for TSP from practical point of view, it would be interesting to investigate other formulations. A possible way would be to dene an undirected model using, e.g., the subtour elimination constraints of Dantzig-Fulkerson-Johnson. Even if such a formulation makes the use of ILP solvers impossible, it can form a basis for designing dedicated branch-and-cut algorithms, which can be expected to outperform the approach based on general ILP solvers.