Solving the clustered traveling salesman problem via traveling salesman problem methods

The Clustered Traveling Salesman Problem (CTSP) is a variant of the popular Traveling Salesman Problem (TSP) arising from a number of real-life applications. In this work, we explore a transformation approach that solves the CTSP by converting it to the well-studied TSP. For this purpose, we first investigate a technique to convert a CTSP instance to a TSP and then apply powerful TSP solvers (including exact and heuristic solvers) to solve the resulting TSP instance. We want to answer the following questions: How do state-of-the-art TSP solvers perform on clustered instances converted from the CTSP? Do state-of-the-art TSP solvers compete well with the best performing methods specifically designed for the CTSP? For this purpose, we present intensive computational experiments on various benchmark instances to draw conclusions.


INTRODUCTION
The Clustered Traveling Salesman Problem (CTSP), originally proposed by Chisman (1975), is an extension of the classic Traveling Salesman Problem (TSP) where the cities are grouped into clusters and the cities of each cluster must be visited contiguously. Formally, the problem is defined on a symmetric complete weighted graph G = (V, E) with a set of vertices V = {1,2,…,n} and a set of edges E = {(i, j):i, j ∈ V, i ≠ j}. The vertex set V is partitioned into disjoint clusters V 1 ,V 2 ,…,V m (V 1 ∪ V 2 ∪…∪ V m = V). Let C be an n × n symmetric distance matrix such that c ij (i, j = 1,2…,n, i ≠ j) represents the travel cost between two corresponding vertices i and j, and satisfies the triangle inequality rule. The objective of the CTSP is to find a minimum cost Hamiltonian circuit over all the vertices, where the vertices of each cluster must be visited consecutively. Figure 1 shows a feasible solution for a CTSP instance, where the solution corresponds to a Hamiltonian cycle such that the vertices of each cluster are visited contiguously.
The CTSP can be formally modelled as the following integer programming model described in Chisman (1975) where, without loss of generality, the salesman is assumed to leave origin city 1 and return to 1.
u i À u j þ ðn À 1Þx ij n À 2 2 i 6 ¼ j n (4) X i2V k X j2V k x ij ¼ jV k j À 1 8V k & V; jV k j ! 1; k ¼ 1; 2; …; m (5) x ij 2 f0; 1g 8i; j 2 V (6) u i ! 0 2 i n In this model, the binary variable x ij = 1 if city j is visited immediately after city i; x ij = 0 otherwise. Objective function (1) seeks to minimize the total distance traveled by the salesman. Constraints (2) and (3) ensure that each city is visited exactly once. Constraints (4) eliminate subtours, while constraints (5) guarantee that the cities of each cluster are visited contiguously. The remaining constraints are related to the decision variables.
The above subtour elimination constraints (4) are called MTZ formulation (Miller, Tucker & Zemlin, 1960). Although MTZ is simple to implement, it provides a very poor linear relaxation (Campuzano, Obreque & Aguayo, 2020). Many compact formulations have been proposed to replace Constraints (4). According to the literature, a multi-commodity flow formulation (Wong, 1980;Claus, 1984) was proven to provide a strong linear relaxation, without compromising its simplicity. In the multi-commodity flow formulation, let k = 2,3,…,n be n − 1 commodities, and let y k ij be a nonnegative decision variable which represents the flow on the arc (i, j)∈ E for the commodity k from city 1 to city k. Then, another alternative mathematical model for the CTSP is constituted of the objective function (1) and the constraints (2), (3), (5), (6) along with the following subtour elimination constraints: Constraints (8) only allow flow in an arc (i, j) if and only if it is traversed by the salesman (i.e., x ij = 1). Constraints (9) ensure that city 1 is the source of one unit of each commodity k ∈ V \{1} and Constraints (10) avoid that the flow of each commodity k ∈ V \{1} returns to city 1. Constraints (11) and (12) guarantees that one flow unit of commodity k enters to city k and it does not leave the city k. Constraints (13) ensure flow conservation at each city, apart from city 1 and for commodity k at city k.
One notices that the CTSP is equivalent to the TSP when there is a single cluster or when each cluster contains exactly one vertex. Therefore, the CTSP is NP-hard, and thus computationally challenging in the general case. From a practical perspective, the CTSP is a versatile modeling tool for several operational research applications arising in a wide variety of areas, including automated warehouse routing (Chisman, 1975), emergency vehicle dispatching (Weintraub et al., 1999), production planning (Lokin, 1979), disk defragmentation (Laporte & Palekar, 2002), and commercial transactions with supermarkets, shops and grocery suppliers (Ghaziri & Osman, 2003). As a result, effective solution methods for the CTSP can help to solve these practical problems. Indeed, the computational challenge and the wide range of applications of the problem have motivated a variety of approaches that are reviewed in the "Literature Review on Existing Solution Methods" section. However, unlike the classic TSP problem for which many powerful methods have been introduced in the past decades, studies on the CTSP are still quite limited.
Moreover, the CTSP belongs to the large class of traveling salesman problems. Among the TSP variants, the generalized traveling salesman problem (GTSP) (Srivastava et al., 1969;Cosma, Pop & Cosma, 2021) and the family traveling salesman problem (FTSP) (Morán-Mirabal, Velarde & Resende, 2014;Pop, Matei & Pintea, 2018) share similarities with the CTSP. In the GTSP, the set of vertices is divided into clusters and the objective is to find a minimum-cost tour passing through one vertex from each cluster. In the FTSP, the set of vertices is also divided into clusters (called families) and the objective is to visit a predefined number of vertices in each family at a minimum cost.
In this work, we investigate the problem transformation approach proposed in Chisman (1975), which converts the CTSP to the TSP and assess the interest of popular modern TSP solvers for solving the resulting TSP instances. To our knowledge, this is the first large computational study testing modern TSP solvers on solving the CTSP. The work is motivated by the following considerations. First, intensive researches on the TSP have led to the development of many very powerful solvers. Thus, it is interesting to know whether we can take advantage of these solvers to effectively solve the CTSP. Second, the TSP instances converted from the CTSP are characterized by their cluster structures. These instances constitute interesting test cases for existing TSP solvers. This work aims thus to answer the following questions.
1. How do state-of-the-art exact TSP solvers perform on clustered instances converted from the CTSP?
2. How do state-of-the-art inexact (heuristic) TSP solvers perform on clustered instances converted from the CTSP?
3. Do state-of-the-art TSP solvers compete well with the best performing methods specifically designed for the CTSP?
To our knowledge, Questions 1 and 3 have never been investigated in the literature. Regarding Question 2, two previous studies (Neto, 1999;Helsgaun, 2014) are of interest. However, they are limited because they only concern one TSP algorithm, i.e., the local search based LKH solver (Helsgaun, 2009), while ignoring other powerful TSP solvers like GA-EAX (Nagata & Kobayashi, 1997) and Concorde . Answering these questions helps to enrich the state-of-the-art of solving the CTSP and gain novel knowledge on using modern TSP methods to solve new problems. Finally, we mention that the transformation approach was also tested in Lokin (1979) and Jongens & Volgenant (1985). However, these studies are clearly outdated and don't provide useful information as to the questions we want to investigate.
The remainder of this paper is organized as follows. "Literature Review on Existing Solution Methods" reviews existing solution methods for the CTSP. "Solving the CTSP via TSP Methods" presents the transformation of the CTSP to the TSP and three powerful TSP methods (solvers). "Computational Experiments" shows computational studies of the TSP solvers applied to the clustered instances and comparisons with existing algorithms dedicated to the CTSP. "Discussion" provides additional explanations regarding the behaviors of the three TSP solvers. Finally, concluding remarks are provided in "Conclusion".

LITERATURE REVIEW ON EXISTING SOLUTION METHODS
There are several dedicated solution algorithms for solving the CTSP that are based on exact, approximation, and metaheuristic approaches.
Along with the introduction of the CTSP, Chisman (1975) proposed a branch-andbound algorithm to solve the integer programming model presented in the Introduction section. Jongens & Volgenant (1985) developed an algorithm based on the 1-tree relaxation to provide lower bounds as well as a heuristic to find satisfactory upper bounds. Mestria, Ochi & de Lima Martins (2013) used the mathematical formulation of Chisman (1975) and IBM Parallel CPLEX solver (version 11.2) to obtain lower bounds for medium CTSP instances (|V| ≤ 1,000).
Various a-approximation algorithms (Anily, Bramel & Hertz, 1999;Gendreau, Laporte & Hertz, 1997;Guttmann-Beck et al., 2000) have been developed for the CTSP. These approximation algorithms require either the starting and ending vertices in each cluster or a prespecified order of visiting the clusters in the tour as inputs, and solve the inter-cluster and intra-cluster problems independently. Bao & Liu (2012) presented a new 2.17approximation algorithm where no starting and ending vertices were specified. Later, Bao et al. (2017) provided a 2.5-approximation algorithm for another version of the CTSP where the starting vertex of each cluster is given while the ending vertex is not specified. Recently, Kawasaki & Takazawa (2020) improved the approximation ratio for the CTSP by incorporating a recent approximation algorithm for the TSP by Zenklusen (2019).
Given that the CTSP is a NP-hard problem, a number of heuristic and metaheuristic algorithms have also been investigated, which aim to provide high-quality solutions in acceptable computation time, but without provable optimal guarantee of the attained solutions. For example, Laporte, Potvin & Quilleret (1997) presented a tabu search algorithm to solve a particular case of the CTSP, where the clusters are visited in a prespecified order. Potvin & Guertin (1996) developed a genetic algorithm for the CTSP that finds inter-cluster paths and then intra-cluster paths. Later, Ding, Cheng & He (2007) proposed a two-level genetic algorithm for the CTSP. In the first level, a genetic algorithm is used to find the shortest Hamiltonian cycle for each cluster. In the second level, a modified genetic algorithm is applied to merge the Hamiltonian cycles of all the clusters into a complete tour.
In addition to these early heuristic algorithms, Mestria, Ochi & de Lima Martins (2013) investigated GRASP (Greedy Randomized Adaptive Search Procedure) with pathrelinking. Among the six proposed heuristics, one heuristic corresponds to the traditional GRASP procedure whereas the other heuristics include different path relinking procedures. Mestria (2016) studied a hybrid heuristic, which is based on a combination of GRASP, Iterated Local Search (ILS) and Variable Neighborhood Descent (VND). Recently, Mestria (2018) presented another complex hybrid algorithm (VNRDGILS) which mixes GRASP, ILS, and Variable Neighborhood Random Descent to explore several neighborhoods. According to the computational results reported in Mestria, Ochi & de Lima Martins (2013) and Mestria (2016Mestria ( , 2018, these GRASP-based algorithms are among the best performing heuristics specially designed for the CTSP currently available in the literature. In addition, Hà et al. (2022) proposed a metaheuristic method based on the ILS framework with problem-tailored operators for a version of the CTSP where the order of visiting the clusters is prespecified.
Existing studies have significantly contributed to better solving the CTSP. According to the computational results reported in the literature, due to the NP-hardness of the problem, only small CTSP instances were able to be solved to optimality with the exact algorithms. The approximation algorithms provide solutions for the CTSP within a given approximation factor. However, due to the high approximation factors involved (e.g., 5/3 (Anily, Bramel & Hertz, 1999), 3/2 Gendreau, Laporte & Hertz (1997), 2.17 (Bao & Liu, 2012), and 2.5 (Bao et al., 2017)), these approximation algorithms are not practical for solving large instances. To deal with large CTSP instances, heuristic and metaheuristic algorithms are often preferred to find sub-optimal solutions within an acceptable computation time.

SOLVING THE CTSP VIA TSP METHODS
Transformation of the CTSP to the TSP As the literature review shows, a number of dedicated solution approaches have been developed to solve the CTSP. However, one observes that these approaches have difficulty producing robustly and consistently high-quality solutions for large-scale CTSP instances with tens of thousands of vertices. Moreover, the best performing CTSP methods (e.g., VNRDGILS (Mestria, 2018), HHGILS (Mestria, 2016), and GPR1R2 (Mestria, Ochi & de Lima Martins, 2013)) are computationally expensive (e.g., requiring 1,080 s to find good solutions for instances with 1,173 ≤ n ≤ 2,000).
On the other hand, problem transformation has been highly successful in solving several difficult optimization problems such as the latin square completion problem via graph coloring (Jin & Hao, 2019) and the winner determination problem via weighted maximum cliques (Wu & Hao, 2015). It is known that the CTSP can be transformed to the conventional TSP (Chisman, 1975). Therefore, in principle, the CTSP can be solved by any TSP algorithm. However, to our knowledge, no computational study on using problem transformation to solve the CTSP has been presented in the literature. This work fills the gap by exploring the problem transformation approach of Chisman (1975) and testing three representative state-of-the-art TSP solvers including both exact and inexact solution approaches.
The basic idea of transforming the CTSP to the TSP is to add a large artificial cost M to all inter-cluster edges in order to force the salesman to visit all the cities within each cluster before leaving it.
Given a CTSP instance G = (V, E) with distance matrix C, we define a TSP instance G′ = (V′, E′) with distance matrix C′ as follow.
Define V = V′ and E = E′.
Define the travel distance c′ ij in G′ by Obviously, if the value of M is sufficiently large, then the best Hamiltonian cycle in G′ is a feasible CTSP solution in G, in which the vertices of each cluster are visited contiguously.
Property. An optimal solution to the TSP instance corresponds to an optimal solution to the original CTSP instance.
Proof. Let S′ and S be the optimal solutions of the TSP instance G′ and the original CTSP instance G, respectively. Let m be the number of clusters of G. To minimize the total travel cost, there are only m inter-cluster edges in S′. Therefore, S′ is a feasible CTSP solution for G and satisfies the following relation: Obviously, S′ corresponds to S by subtracting the constant m × M.

Solution methods for the TSP
There are numerous solution methods for the TSP. In this work, we adopt three very powerful TSP solvers whose codes are publicly available, including one exact solver (Concorde ) and two inexact (heuristic) solvers (LHK-2 (Helsgaun, 2009) and GA-EAX (Nagata & Kobayashi, 2013)).
Notice that the TSP instance converted from a CTSP instance has a particular feature that the vertices are grouped into clusters and the distance between each pair of vertices within a same cluster is in general small, while this distance is large for two vertices from different clusters. Along with the presentation of the TSP solvers, we discuss their suitability for solving such clustered instances each time this is appropriate.

Exact Concorde solver
Concorde is an advanced exact TSP solver for the symmetric TSP based on Branch-and-Bound and problem specific cutting plane methods . It makes use of a specifically designed QSopt linear programming solver. According to Hoos & Stützle (2014), Concorde is the best performing exact algorithm for the TSP. As shown in , Concorde can solve benchmark instances from TSPLIB with up to 1,000 vertices to optimality within a reasonable computation time and it also solves large TSP instances at the cost of a long computation time.
The run time behavior of Concorde has been investigated essentially on random uniform instances. For instance,  investigated the run time required by Concorde for solving random uniform instances and indicated that the run time increases as an exponential function of instance size |V|. Hoos & Stützle (2014) further demonstrated that the median run time required by Concorde scales with instance size |V| of the form ab ffiffiffiffi ffi jVj p (a ≈ 0.21, b ≈ 1.24) on the widely studied class of uniform random TSP instances. To our knowledge, no study has been reported concerning the behavior of Concorde on sharply clustered instances. As a result, the current study will provide useful information on this issue.

Lin-Kernighan based heuristic solver
According to the TSP literature, a majority of the best performing TSP heuristic algorithms is based on the Lin-Kernighan (LK) heuristic (Lin & Kernighan, 1973) and its extensions. The LK heuristic is a variable-depth k-opt local search procedure, where the k-opt neighborhood is partially searched with a smart pruning strategy. LK explores the most promising neighbors within the k-opt neighborhood, that is, the set of feasible tours obtained by removing k edges and adding other k edges such that the resulting tour is feasible. Several improved versions of the basic LK heuristic have been introduced within the iterated local search framework (e.g., Applegate, Cook & Rohe, 2003;Helsgaun, 2000;Helsgaun, 2009;Martin, Otto & Felten, 1991).
Among these iterated LK algorithms, Helsgaun's LKH (Helsgaun, 2000(Helsgaun, , 2009) is the uncontested state-of-the-art heuristic TSP solver. Helsgaun (2000) developed an iterated version of LK together with an efficient implementation of the LK algorithm, known as the Lin-Kernighan-Helsgaun (LKH-1) heuristic, where a 5-opt move is used as the basic move to broaden the search and an α-measure method based on sensitivity analysis of minimum spanning trees is used to restrict the search to relative few of the α-nearest neighbors of a vertex to speed up the search process. Later, Helsgaun (2009) further extended LKH-1 by developing a highly effective implementation of the k-opt procedure (called LKH-2), which eliminated many of the limitations and shortcomings of LKH-1. Furthermore, LKH-2 specially extended the data structures of LKH-1 to solve very large TSP instances. The main features of LKH-2 include (1) using sequential and nonsequential k-opt moves, (2) using several partitioning procedures to partition a large TSP instance into smaller subproblems, (3) using a tour merging procedure to generate a better solution from two or more local optimum solutions, and (4) applying a backboneguided search to guide the local search to make biased local perturbations. LKH-2 is considered to be one of most effective heuristic methods for finding very high-quality solutions for various large TSP instances (Dubois-Lacoste, Hoos & Stützle, 2015).
However, the LK algorithm and any LK-based algorithms require much longer running times on clustered instances of the TSP than on uniformly distributed instances (Neto, 1999). The main reason why the LK heuristic stumbles on clustered instances is that relatively large inter-cluster edges serve as bait edges. During the LK search, when removing such a bait edge, the LK heuristic is tricked into long and often fruitless searches. More precisely, each time an edge bridging two clusters is removed, the cumulative gain rises enormously, and the procedure is encouraged to perform very deep searches. To alleviate the problem, a cluster compensation technique was proposed in Neto (1999) for the Lin-Kernighan heuristic to limit its performance degradation. Helsgaun (2009) showed that the LKH-2 algorithm performs significantly worse on sharply clustered instances than on uniform random instances. To remedy this difficulty, Helsgaun (2014) considered the unusual structure of clustered instances, and adjusted the parameter settings of LKH-2 to better solve the clustered instances. The resulting solver is named CLKH, which is used in this study.

Edge assembly crossover based genetic algorithm
Population-based evolutionary algorithms are another well-known approach for the TSP. A popular example is the powerful genetic algorithm introduced by Nagata & Kobayashi (2013). This algorithm (called GA-EAX, see Algorithm 1) is characterized by its powerful edge assembly crossover (EAX) operator introduced in Nagata & Kobayashi (1997) and Nagata & Soler (2012) with an efficient implementation and a cost-effective selection strategy for maintaining population diversity.
The key EAX operator generates, from two high-quality tours (parents), one offspring tour by first inheriting the edges from the parents to construct disjoint subtours and then connecting the subtours with new edges in a greedy fashion (similar to building a minimal spanning tree). Let S A and S B be the parents, EAX operates as follows (see Fig. 2 for an example): the sets of edges of parents S A and S B , respectively.
2. Extract all AB-cycles from G AB . An AB-cycle is defined as a cycle in G AB , such that edges of E A and edges of E B are alternately linked.
3. Construct an E-set by selecting AB-cycles according to a given selection strategy (e.g., single, k-multiple, block and block2 (Nagata & Kobayashi, 2013)), where an E-set is a set of AB-cycles. Note that different versions of EAX can be developed by using different selection strategies of AB-cycles for constructing E-sets. The GA-EAX algorithm employs the single and block2 strategies to generates offspring solutions from parent solutions. To maintain a healthy population diversity, GA-EAX also uses an edge entropy measure to select the solution to be used to replace a parent in the population.

Copy parent S
Other studies (e.g., Hains, Whitley & Howe, 2012) also indicated the usefulness of edge-assembly-like crossovers for solving clustered instances of the TSP. As shown in the next section, the EAX-based genetic algorithm performs remarkably well on all the clustered instances transformed from the CTSP.

COMPUTATIONAL EXPERIMENTS
In this section, we evaluate the capacity of the TSP solvers presented in "Solution Methods for the TSP" to solve the CTSP via its transformation to the TSP. For this purpose, we examine their qualitative performances and run time efficiencies on various benchmark instances and make comparisons with the best dedicated CTSP algorithms in the literature.
Set 3 (38 instances): These large instances have 1,000 to 24,978 vertices and come from GTSPLIB for the generalized traveling salesman problem (GTSP). They were generated from TSP instances by using Fischetti, Salazar González & Toth (1997) clustering algorithm and tested in Helsgaun (2014) by considering them as CTSP instances. These instances are available at http://www.ruc.dk/~keld/research/CLKH. In Helsgaun (2014), six very large instances with 31,623 to 85,900 vertices were also tested. We ignore these instances, because they are too large for the exact Concorde solver and the GA-EAX solver stops abnormally when solving these instances.

TSP solvers and experimental protocol
For our study, we employed three representative TSP solvers presented in "Solution Methods for the TSP", which are among the most powerful methods for the TSP in the literature.
Exact Concorde TSP solver (http://www.math.uwaterloo.ca/tsp/concorde/index.html): We used version Concorde-03.12.19 and ran the solver with its default parameter setting with a cutoff time of 24 CPU hours per instance. Inexact CLKH solver (http://www.ruc.dk/~keld/research/CLKH): We used the version CLKH-1.0 which is based on the latest version 2.0.9 (http://akira.ruc.dk/~keld/ research/LKH/) of LKH-2. The default parameter setting given in Helsgaun (2014) was adopted to run CLKH. Notice that to reduce its run time, the maximum number of trials (iterations) is set to 1,000 in CLKH, while this number is set to n (instance size) by default in LKH-2. Inexact GA-EAX TSP solver (https://github.com/sugia/GA-for-TSP): We used GA-EAX with its default parameter setting given in Nagata & Kobayashi (2013): p = 300, r = 30 and GA-EAX terminates if the difference between the average tour length and the shortest tour length in the population is less than 0.001. Following Kerschke et al. (2018) and Kotthoff et al. (2015), we reset the random seed for GA-EAX for each run (which was set to a fixed value in the official implementation).
The experiments were carried out on a computer running Linux operating system with an Intel E5-2670 processor (2.8 GHz and 4G RAM). Given the stochastic nature of CLKH and GA-EAX, we ran each algorithm 10 times for each instances while the deterministic Concorde TSP solver was run one time to solve each instance.

Computational results and comparison of popular TSP solvers
Our computational studies aim to answer the following questions: How do state-of-the-art exact TSP solvers perform on clustered instances converted from the CTSP? How do state-of-the-art inexact (heuristic) TSP solvers perform on clustered instances converted from the CTSP?
The results of the three TSP solvers (Concorde, CLKH, GA-EAX) on the 20 medium and 15 large CTSP benchmark instances are summarized in Table 1 (Set 1) and Table 2 (Set 2). Columns 1 to 3 show the basic information of each instance: the instance name (Instance), the number of vertices (|V|) and the number of clusters (m). Column 4 gives the optimal objective value reported by the exact Concorde TSP solver, followed by the required run time in seconds. For both the CLKH and GA-EAX solvers, we show the best (Gap best ) and average (Gap avg ) results over 10 independent runs in the form of the percentage gap to the optimal solution, as well as the average run time in seconds. If the best solution over 10 independent runs equals the optimal solution obtained with the exact Concorde TSP solver, the corresponding cell in column Gap best shows '=' along with the number of runs that succeeded in finding the optimal solution. Finally, row 'Avg.' provides the average run time in seconds for each approach, and the average gap between the average objective values obtained with CLKH/GA-EAX and the optimal values obtained with the Concorde TSP solver. From Tables 1 to 2, we can make the following observations: First, the exact Concorde TSP solver performs very well on these 35 instances and is able to solve all of them exactly. Specifically, the 20 medium instances can be solved easily in a short run time (an average of about 14 s). The 15 large instances are more difficult and the run time needed to solve these instances increases considerably (an average of 1,133.8 s, reaching 7,214.3 s for the most difficult instance).
Second, the CLKH solver performs globally very well on these 35 instances. For the 20 medium instances, CLKH attains all the optimal solutions but one with an average run time of 47.8. For the 15 large instances, CLKH reaches the optimal solutions for 13 instances with an average run time of 257.3 s.
Third, the GA-EAX solver performs remarkably well by attaining the optimal values for all 35 instances but one. For the 20 medium instances, GA-EAX consistently hits the optimal solutions for each of its 10 run (except for one instance for which it has a hit of 9 out of 10). The average run time is only 5.7 s for the medium instances and 33.6 s for the large instances. Compared to Concorde and CLKH, GA-EAX is thus extremely time efficient. Moreover, contrary to the Concorde and CLKH solvers, the computation time required by GA-EAX remains very stable across the instances of the same set, indicating a high robustness and scalability of this solver. Table 3 presents the results of the three TSP solvers on the 38 large GTSP instances of Set 3. Notice that the Concorde solver failed to exactly solve 17 instances in 24 h, the corresponding cell (in parentheses) in column 'Optimum' indicates the best tour length (best upper bound) found by CLKH and GA-EAX. In this case, the percentage gaps (Gap best and Gap avg ) are calculated by using the best bound, and column Gap best shows '=' the number of runs for an algorithm to find the best bound.
From Table 3, we can make the following observations. First, Concorde manages to optimally solve 21 large GTSP instances with up to 3,162 vertices with a run time ranging from 17.4 s to 45,008.4 s while its solution time is not completely consistent with the size of the problem instances. For the 21 instances that can be solved exactly by Concorde, CLKH attains 15 best upper bounds, while GA-EAX reaches all best upper bounds in less computing time. Second, for most of the instances with |V| < 10,000, compared with CLKH, GA-EAX has a better performance both in terms of solution quality and computation time. For the instances with 10,000 ≤ |V| ≤ 24,978, the solution quality of GA-EAX is better than that of CLKH in most cases, while requiring more computation time.
To sum up, the exact Concorde solver is very efficient for the instances with up to 1,000 vertices (order of seconds) and can even find optimal solutions for instances with up to some 3,000 vertices at a price of more run time (order of minutes to hours). For larger instances, both inexact solvers (CLKH and GA-EAX) are reliable alternatives to find optimal or sub-optimal solutions with some advantages for GA-EAX. These heuristic solvers also perform very well on smaller instances.
To deepen our computational study, we call upon to the performance profile, an analytic tool for evaluating the performances of multiple compared optimization algorithms (Dolan & Moré, 2002). The performance profile uses a cumulative distribution function for a performance metric, such as run time, objective function values, number of iterations, and so on. Precisely, let S be a set of algorithms and P be a set of problem instances. For a given performance metric f s,p (that is the performance of algorithm s ∈ S solving instance p ∈ P), the performance ratio is defined by r s;p ¼ f s;p minff a;p : a 2 Sg . Then, for each algorithm s ∈ S, the performance function is given by q s ðsÞ ¼ jfp 2 P : r s;p sgj jPj .
Thus, the value of ρ s (1) corresponds to the fraction of problem instances that algorithm s can achieve many times the performance of the best algorithm (meaning the probability that the algorithm s will win over the rest of the compared algorithms). For a large value τ, the value of ρ s (τ) indicates a high robustness of algorithm s.
To make a fair and meaningful comparison with this tool, we focus on the two inexact solvers CLKH and GA-EAX and run each solver 10 times on each of the 73 instances. We use the software 'perprof-py' (Siqueira, da Silva & Santos, 2016) to draw the performance profiles (see Fig. 3) where the quality of the solution is measured by the average objective value and average run time. These performance profiles tend to show an advantage of GA-EAX over CLKH for solving these clustered instances with up to 24,978 vertices.

TSP solvers v.s. state-of-the-art CTSP heuristics
In "Computational Results and Comparison of Popular TSP Solvers", we observed that the exact Concord TSP solver and the inexact CLKH and GA-EAX TSP solvers are powerful tools for solving clustered TSP instances converted from the CTSP. We now answer the following question: Do these general TSP solvers compete well with state-of-the-art CTSP heuristics specially designed for the problem? For this purpose, we adopt GA-EAX as our representative TSP solver and compare it with three best performing CTSP heuristics in the literature: VNRDGILS (Mestria, 2018), HHGILS (Mestria, 2016), and GPR1R2 (Mestria, Ochi & de Lima Martins, 2013). Indeed, according to the experimental studies reported in Mestria, Ochi & de Lima Martins (2013) and Mestria (2016, 2018, these three heuristics perform the best among the recent CTSP heuristics available in the literature (see Table 4). This study is based on the 35 medium and large instances of Sets 1 and 2 (no results for the three CTSP heuristics are available on the large GTSP instances of Set 3). Table 5 provides the comparative results of the GA-EAX TSP solver along with the results reported by the three CTSP algorithms on the medium and large instances. For each instance and algorithm, columns 'f best ', 'f avg ' and 't(s)' show respectively the best objective value over 10 independent runs, the average objective value and the average run time in seconds. Furthermore, the row 'Avg.' shows the average performances for each compared algorithm, including the average percentage gap of the best/average result to the optimal result obtained with the Concorde TSP solver and the average run time in seconds.
To determine whether there exists a statistically significant difference in performance between the GA-EAX TSP solver and each CTSP algorithm in terms of best and average results, the p-values from the Wilcoxon signed-rank tests are given in the last row of Table 5. Entries with "-" mean that the corresponding results are not available in the literature. The best objective values obtained by the compared algorithms are indicated in bold if they attain the optimal solution. Notice that the results of the CTSP algorithms (VNRDGILS, HHGILS and GPR1R2) correspond to 10 executions per instance on a computer with 2.83 GHz Intel Core 2 CPU and 8 GB RAM and the time limit per run was set to 720 s for medium instances and 1,080 s for large instances. Table 6 summarizes the statistical results for each compared algorithm on the two sets of medium and large instances. The first row indicates the number of optimal solutions found by each approach. The average percentage gap of the best/average result from the optimal result is provided in row 'Average Gap best /Gap avg '. Finally, row 'Average time (s)' provides the average run time in seconds for each algorithm.
From Tables 5 and 6, we observe that the GA-EAX solver significantly outperforms the three CTSP algorithms on the medium and large instances in terms of both the best and the average results. For the large instance set, the improvement gaps between the results of GA-EAX and those of the CTSP methods are very high, ranging from 10.39% to 15.49%. Furthermore, in terms of the average run time, GA-EAX is about 30 to 130 times faster than the CTSP algorithms. The above results thus indicate that the GA-EAX TSP solver has a strong dominance over current best performing CTSP approaches in the literature. In addition, the small p-values (<0.05) from the Wilcoxon signed-rank tests further confirm the statistically significant difference of the compared results.
To have a finer analysis of the compared algorithms, Fig. 4 provides boxplot graphs to compare the distribution and range of the average results for each compared algorithm, except GPR1R2 for the medium instances since its results on several medium instances are not available. In this figure, the average objective value f avg of a given algorithm is  (2007) A two-level genetic algorithm normalized according to the relation y = 100 Ã (f avg − f opt )/f opt , where f opt is the optimal value. The plots in Fig. 4 show clear differences in the distributions of the average results between GA-EAX and each compared CTSP heuristic, which further confirms the efficiency of the GA-EAX TSP solver with respect to these dedicated CTSP heuristics. Finally, considering the results of the Concorde solver and the CLKH solver reported in "Computational Results and Comparison of Popular TSP Solvers", we conclude that these TSP solvers also dominate the current best CTSP algorithms in the literature.

DISCUSSION
We now provide additional explanations regarding the behaviors of the three TSP solvers. First, given the NP-hard nature of the CTSP and the exponential time complexity of the exact Concorde solver, it is expected that the exact Concorde solver reaches its limit when the instance to be solved reaches some size (about 3,000 vertices for the studies instances). Indeed, when the search space becomes extremely large, the exact Branch-and-  Bound search even equipped with the best problem specific cutting plane methods cannot effectively enumerate all candidate solutions. In fact, such a behavior has already been observed in previous studies on Concorde applied to classical TSP instances Hoos & Stützle, 2014). Second, regarding the two heuristic solvers CLKH and GA-EAX, the CLKH solver exhibits a worse performance compared to GA-EAX. As discussed in "Lin-Kernighan Based Heuristic Solver", the underlying LK heuristic stumbles on clustered instances because relatively large intercluster edges serve as bait edges. With the presence of these bait edges, the LK heuristic may be tricked into long and often fruitless search trajectories. Third, the GA-EAX solver performs its search mainly with its edge assembly crossover, which inherits the edges of the parents to construct disjoint subtours and then connect the subtours. This crossover proves to be meaningful and helps the algorithm avoid local optimal traps. Once again, the excellent behavior of GA-EAX on the CTSP instances is consistent with its performance on conventional TSP instances as shown in Nagata & Kobayashi (2013).

CONCLUSION
This work presents the first extensive computational study on the transformation approach of solving the Clustered Traveling Salesman Problem with general TSP solvers. Based on the results from the exact Concorde solver and the heuristic CLKH and GA-EAX solvers on 20 medium (101 ≤ |V| ≤ 1,000) and 15 large (1,173 ≤ |V| ≤ 2,000) CTSP benchmark instances and 38 large GTSP benchmark instances (with up to 24,978 vertices) available in the literature, we can draw the following conclusions.
The exact Concorde solver can optimally solve all medium and large CTSP instances. It also solves exactly large GTSP instances with up to 3,162 vertices in a reasonable time, but fails to solve larger GTSP instances in 24 h. Its solution time is not completely consistent with the size of the problem instances. The heuristic CLKH and GA-EAX solvers perform very well both in terms of solution quality and computational efficiency. Both solvers have a good scalability, making them particularly suitable for solving very large instances with at least several thousands of vertices. For the tested instances with up to some 24,978 vertices, GA-EAX exhibits a better performance than CLKH. The general TSP solvers significantly dominate, both in terms of solution quality and computational efficiency, the current best performing CTSP heuristics specially designed for the problem. In particular, the TSP heuristics are several orders of faster than the state-of-the-art CTSP heuristics to find much better results.
This study indicates that the existing CTSP benchmark instances in the literature are not challenging for modern TSP solvers even if they remain difficult for the existing CTSP algorithms.
Finally, given the findings of this study, it would be interesting to investigate the problem transformation approach for solving other TSP variants that can be converted to the TSP or to a TSP variant for which effective algorithms are available.