Qing Zhou
email: qingzhou@hust.edu.cn

Una Benlic

Qinghua Wu
email: qinghuawu1005@gmail.com

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

T S Grasp+ts

Heuristic search to the capacitated clustering problem

Keywords: Tabu search, memetic algorithm, infeasible local search, capacitated clustering

Given a weighted graph, the capacitated clustering problem (CCP) is to partition a set of nodes into a given number of distinct clusters (or groups) with restricted capacities, while maximizing the sum of edge weights corresponding to two nodes from the same cluster. CCP is an NP-hard problem with many relevant applications. This paper proposes two effective algorithms for CCP: a Tabu Search (denoted as FITS) that alternates between exploration in feasible and infeasible search space regions, and a Memetic Algorithm (MA) that combines FITS with a dedicated cluster-based crossover. Extensive computational results on five sets of 183 benchmark instances from the literature indicate that the proposed FITS competes favorably with the state-of-the-art algorithms. Additionally, an experimental comparison between FITS and MA under an extended time limit demonstrates that further improvements in terms of the solution quality can be achieved with MA in most cases. We also analyze several essential components of the proposed algorithms to understand their importance to the success of these approaches.

Introduction

Given a weighted graph G = (V, E) where V is a set of n nodes and E is a set of edges, let w i ≥ 0 be the weight of node i ∈ V and let c ij ({i, j} ∈ E) be the edge weight between nodes i and j (c ij = 0, if {i, j} / ∈ E). The Capacitated Clustering Problem (CCP) is to partition V into a given number p (p ≤ n) of disjoint clusters or groups such that the sum of node weights in each cluster is constrained by an upper and a lower capacity limit, while maximizing the sum of edge weights whose two associated endpoints belong to the same cluster.

Formally, let the binary variable X ig take the value 1 if node i is assigned to group g (g ∈ {1, 2, ..., p}), and 0 otherwise. CCP can then be expressed as the following quadratic program [START_REF] Brimberg | Solving the capacitated clustering problem with variable neighborhood search[END_REF][START_REF] Deng | A reactive GRASP with path relinking for capacitated clustering[END_REF]:

maximize p g=1 n-1 i=1 n j=i+1 c ij X ig X jg (1)
subject to p g=1

X ig = 1, ∀i ∈ V (2)
L g ≤ n i=1 w i X ig ≤ U g , ∀g ∈ {1, 2, ..., p} (3)
X ig ∈ {0, 1}, ∀i ∈ V, g ∈ {1, 2, ..., p} (4)
Constraint [START_REF] Bard | Large-scale constrained clustering for rationalizing pickup and delivery operations[END_REF] guarantees that every node is assigned to exactly one cluster, while constraint [START_REF] Bartz-Beielstein | Experimental Methods for the Analysis of Optimization Algorithms[END_REF] ensures that the minimum capacity (L g ) and the maximum capacity (U g ) requirements of each cluster are satisfied.

Note that CCP is closely related to the Graph Partitioning Problem (GPP) [START_REF] Benlic | A multilevel memetic approach for improving graph k-partitions[END_REF][START_REF] Benlic | Hybrid Metaheuristics for the Graph Partitioning Problem[END_REF][START_REF] Galinier | An efficient memetic algorithm for the graph partitioning problem[END_REF] where the lower and the upper capacity limits of the clusters are respectively set to 0 and a predetermined imbalance parameter. Moreover, the Maximally Diverse Grouping Problem (MDGP) [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhood search[END_REF][START_REF] Gallego | Tabu search with strategic oscillation for the maximally diverse grouping problem[END_REF][START_REF] Johnes | Operational research in education[END_REF][START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF][START_REF] Palubeckis | Maximally diverse grouping: an iterated tabu search approach[END_REF][START_REF] Rodriguez | An artificial bee colony algorithm for the maximally diverse grouping problem[END_REF][START_REF] Weitz | An empirical comparison of heuristic methods for creating maximally diverse groups[END_REF] is a special case of CCP, when G is a complete graph with unit cost node weights. Consequently, CCP is an NP-hard problem as MDGP is known to be NPhard. Furthermore, CCP is equivalent to the Handover Minimization Problem (HMP) in mobility networks [START_REF] Martínez-Gavara | Tabu search and GRASP for the capacitated clustering problem[END_REF] where the objective is to minimize the sum of weights of the edges with endpoints in different clusters. In addition to its application in mobility networks, CCP arises in several different contexts including mail delivery [START_REF] Bard | Large-scale constrained clustering for rationalizing pickup and delivery operations[END_REF], VLSI design [START_REF] Weitz | An empirical comparison of heuristic methods for creating maximally diverse groups[END_REF] and vehicle routing [START_REF] Fisher | A generalized assignment heuristic for vehicle routing[END_REF][START_REF] Koskosidis | Clustering algorithms for consolidation of customer orders into vehicle shipments[END_REF].

Since it was first defined in 1984 by Mulvey and Beck [START_REF] Mulvey | Solving capacitated clustering problems[END_REF], a variety of solution approaches have been proposed for CCP given its practical importance and NP-hard nature. State-of-the-art approaches include a Greedy Randomized Adaptive Search Procedure with Path Relinking (GRASP-PR) by Deng and Bard [START_REF] Deng | A reactive GRASP with path relinking for capacitated clustering[END_REF]. In 2013, Morán-Mirabal et al. [START_REF] Morán-Mirabal | Randomized heuristics for handover minimization in mobility networks[END_REF] proposed three algorithms for the equivalent handover minimization problem: a GRASP with path-relinking (denoted as GQAP in the corresponding paper), a GRASP with evolutionary path-relinking (GevPR-HMP) and a population-based biased random-key genetic algorithm (BRKGA). According to their computational results, GevPR-HMP exhibits the best performance among those three algorithms. In 2015, Martínez-Gavara et al. [START_REF] Martínez-Gavara | Tabu search and GRASP for the capacitated clustering problem[END_REF] presented several approaches which may be considered as state-of-the-art methods for CCP, including a Greedy Randomized Adaptive Search Procedure (GRASP), a Tabu Search method (TS), a hybrid method combining GRASP with TS (GRASP+TS), and a Tabu Search with Strategic Oscillation (TS SO). In 2016, Lai and Hao [START_REF] Lai | Iterated variable neighborhood search for the capacitated clustering problem[END_REF] proposed a highly effective Iterated Variable Neighborhood Search (IVNS) for CCP. More recently, Martínez-Gavara et al. [START_REF] Martínez-Gavara | Randomized heuristics for the capacitated clustering problem[END_REF] applied several methods to CCP including a GRASP algorithm (denoted as GRASP2-1) and an Iterated Greedy (IG) algorithm, while Brimberg et al. [START_REF] Brimberg | Solving the capacitated clustering problem with variable neighborhood search[END_REF] presented two highly effective VNS-based heuristics denoted as GVNS and SGVNS. A comprehensive review on the most representative approaches for CCP prior to 2011 can be found in [START_REF] Deng | A reactive GRASP with path relinking for capacitated clustering[END_REF].

CCP is a constrained problem that imposes a lower and an upper capacity limit to the size of the clusters. One notices that most of the existing approaches for CCP restrict their search to the feasible region only, while only few approaches including GRASP-PR [START_REF] Deng | A reactive GRASP with path relinking for capacitated clustering[END_REF], GevPR-HMP [START_REF] Morán-Mirabal | Randomized heuristics for handover minimization in mobility networks[END_REF], GQAP [START_REF] Morán-Mirabal | Randomized heuristics for handover minimization in mobility networks[END_REF], and TS SO [START_REF] Martínez-Gavara | Tabu search and GRASP for the capacitated clustering problem[END_REF] are allowed to visit infeasible solutions. In this work, we are interested in search algorithms that examine both the feasible and the infeasible space in search of improved solutions. For this purpose, we introduce a highly effective tabu search (denoted as FITS) that alternates between feasible and infeasible regions, bringing more flexibility into the search process. In addition, we propose the first population-based memetic algorithm (MA) for CCP. It uses FITS as the local optimization mechanism, and incorporates a dedicated cluster-based crossover to transfer pertinent properties ("building blocks") from parents to offspring. Experimental results on five sets of 183 benchmark instances indicate a highly competitive performance of FITS with respect to the existing state-of-the-art algorithms. Given a longer time limit, a computational comparison between FITS and MA reveals that MA is able to further improve on the performance of its underlying FITS in terms of solution quality.

The remainder of the paper is organized as follows. Section 2 presents FITS, followed by a detailed description of the proposed memetic approach in Section 3. Experimental results on a widely-used benchmark are provided in Section 4. Section 5 analyzes the contribution of the key algorithmic ingredient to the performance of the proposed algorithms. Furthermore, we motivate the choice for the crossover used by MA, prior to conclusions drawn in the last section.

Main framework

Algorithm 1 Main scheme of FITS Instead of confining the search process to feasible regions, a number of studies on highly constrained problems [START_REF] Chen | A hybrid metaheuristic approach for the capacitated arc routing problem[END_REF][START_REF] Errico | A priori optimization with recourse for the vehicle routing problem with hard time windows and stochastic service times[END_REF][START_REF] Higgins | A dynamic tabu search for large-scale generalised assignment problems[END_REF][START_REF] Lapierre | Balancing assembly lines with tabu search[END_REF][START_REF] Paraskevopoulos | Resource Constrained Routing and Scheduling: Review and Research Prospects[END_REF] have shown that the consideration of infeasible solutions during the search may help to better explore the search space. Based on this observation, the proposed tabu search alternates between a feasible local search phase (FLS for short) that only examines feasible solutions, and an infeasible local search phase (InfLS for short) where the capacity constraint is relaxed in a controlled manner. The two phases play different roles in the search process -FLS ensures an intensified exploitation in a relevant search region, while InfLS is used to introduce more freedom (diversification) into the search. By alternating between these two complementary phases, FITS is expected to explore various zones of the search space without being easily trapped in a local optimum.

Require: Graph G = (V,
Algorithm 1 summarizes the general framework of FITS. Starting from a feasible solution generated with a construction procedure (Section 2.2), FITS first enters the FLS phase that is based on the best-improvement strategy with a joint use of three types of move operators (Section 2.3.2). This phase terminates as soon as the search is deemed to be trapped in a deep local optimum, i.e., if the best found solution cannot be improved for N cons consecutive iterations. The algorithm then switches to InfLS that relies on a penalty-based evaluation function to guide the search to move towards new search regions. The stopping condition is typically a time limit or a fixed number of iterations.

Initial solution

The starting point for the search is a feasible solution generated by means of two randomized construction methods similar to those used in [START_REF] Gallego | Tabu search with strategic oscillation for the maximally diverse grouping problem[END_REF][START_REF] Lai | Iterated variable neighborhood search for the capacitated clustering problem[END_REF]. The first method consists in two stages, where the first stage performs the following steps: (i) randomly select a node v from the set of unassigned nodes, and randomly choose a cluster g from the set of clusters whose lower capacity constraint is not satisfied; (ii) allocate v to cluster g. The two steps are repeated until all clusters satisfy the lower capacity constraint. Once the first stage is completed, the proposed construction method enters the second stage that: (i) randomly picks an unassigned node v and a cluster g such that size[g] + w v <= U g , where size[g] and w v represent respectively the current weight of cluster g and the weight of node v; (ii) assign v to g. The second stage of this procedure terminates as soon as all the nodes have been assigned.

As observed in our preliminary experiments, the above described method often fails to find a feasible assignment of all the nodes when the upper capacity limit of clusters is very tight. Consequently, we propose the second construction method which constitutes a slight modification of the first method. Instead of randomly choosing a node v in both stages of the first construction method, an unassigned node v is selected such that v has the largest weight (ties broken randomly). The other steps are kept unchanged.

The time complexity of the construction method is O(n * p).

Feasible local search (FLS)

FLS searches for the most promising solutions in the feasible space of candidate solutions, thus ensuring that the capacity constraint is verified. It is based on the general tabu search framework [START_REF] Laguna | Tabu Search[END_REF] and a combined use of three complementary move operators as described in the following subsections.

Feasible search space and evaluation function

A candidate solution to CCP is any partition of the node set V into p subsets C 1 , C 2 ,..., C p , also called clusters. The search space, including both feasible and infeasible solutions, is then formally defined as:

Ω = {{C 1 , C 2 , ..., C p } : ∪ p i=1 C i = V, C i ∩ C j = ∅} (5)
where i = j, 1 ≤ i, j ≤ p. Notice that an infeasible solution may contain one or more empty clusters.

The feasible search space includes the set of all the candidate solutions Ω f ⊂ Ω satisfying the capacity constraints:

Ω f = {{C 1 , C 2 , ..., C p } : L i ≤ |C i | ≤ U i , ∪ p i=1 C i = V, C i ∩ C j = ∅} (6)
where i = j, 1 ≤ i, j ≤ p and |C i | represents the total weight of the nodes in cluster i (i.e.,

|C i | = u∈C i w u ).
To evaluate the quality of each candidate solution s = {C 1 , C 2 , ..., C p } in Ω f , the evaluation function is equivalent to the objective function which sums up the edge weights associated to endpoints in the same cluster:

f (s) = p g=1 i,j∈Cg,i<j c ij (7)

Neighborhood Structures

As previously mentioned, the neighborhood exploited by FLS is defined by a joint use of three basic move operators, which have previously been employed in [START_REF] Brimberg | Solving the capacitated clustering problem with variable neighborhood search[END_REF][START_REF] Deng | A reactive GRASP with path relinking for capacitated clustering[END_REF][START_REF] Lai | Iterated variable neighborhood search for the capacitated clustering problem[END_REF][START_REF] Martínez-Gavara | Tabu search and GRASP for the capacitated clustering problem[END_REF][START_REF] Martínez-Gavara | Randomized heuristics for the capacitated clustering problem[END_REF]. These operators are briefly described as follows:

OneMove operator: Given a solution s = {C 1 , C 2 , ..., C p }, OneM ove transfers a node v from its original cluster i to another cluster j such that the capacity constraint is respected, i.e., |C i |w v ≥ L i and |C j | + w v ≤ U j . To rapidly evaluate the gain value for each candidate move, our algorithm employs a fast incremental evaluation technique similar to that used in [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhood search[END_REF][START_REF] Lai | Iterated variable neighborhood search for the capacitated clustering problem[END_REF][START_REF] Palubeckis | Maximally diverse grouping: an iterated tabu search approach[END_REF][START_REF] Rodriguez | An artificial bee colony algorithm for the maximally diverse grouping problem[END_REF]. The main idea is to maintain an incremental matrix γ, where each element γ[v][g] represents the sum of the edge weights between v and other nodes located in cluster g of the current solution, i.e., γ[v][g] = u∈Cg c uv . Let OneM ove(v, i, j) denote a move that consists in transferring a node v from cluster i to cluster j, the corresponding gain value can be conveniently calculated as:

∆ f (OneM ove(v, i, j)) = γ[v][j] -γ[v][i]
After each OneM ove operation, a subset of values in γ affected by the move is updated as follows:

γ[u][i] = γ[u][i] -c uv , γ[u][j] = γ[u][j] + c uv , ∀u ∈ V .
The complexity to update γ after a OneM ove operation is O(n).

SwapMove operator : This move operator swaps two nodes v and u from two different clusters i and j, such that the capacity constraint is maintained. Let SwapM ove(v, u) denote a swap move, the associated move gain can be efficiently obtained as:

∆ f (SwapM ove(v, u)) = (γ[v][j] -γ[v][i]) + (γ[u][i] -γ[u][j]) -2c vu
Since a SwapM ove can be decomposed into two consecutive OneM ove operations, γ is updated in two steps as for the corresponding OneM ove moves.

Clearly, updating γ after a SwapM ove operation can also be achieved in O(n) time.

2-1 Exchange operator: Let v, u and z be three nodes where v and u are located in the same cluster i, while z belongs to another cluster j. The 2-1 Exchange transfers v and u from cluster i to cluster j and simultaneously moves z from j to i, while respecting the capacity restriction of i and j.

Let Exchange(v, u, z) denote such a move, the resulting move gain can be computed as:

∆ f (Exchange(v, u, z)) = (γ[v][j] -γ[v][i]) + (γ[u][j] -γ[u][i]) + (γ[z][i] - γ[z][j]) + 2(c vu -c vz -c uz )
Since a 2-1 Exchange move can be decomposed into three consecutive OneM ove operations, the matrix γ is consecutively updated three times according to the corresponding OneM ove moves.

Exploration of the feasible search space

The general scheme of FLS is summarized in Algorithm 2. Starting from a feasible solution, FLS selects at each iteration the best non-prohibited move (i.e., non-tabu move of the highest gain) from the union of OneM ove, SwapM ove and 2-1 Exchange moves, where ties are broken at random. Obviously, such a combined neighborhood ensures an intensified examination of the feasible search space, and thus enhances the capacity of finding improved feasible solutions. To avoid short-term cycling, each time a node v is moved from its original cluster C, it is forbidden to move v back to C for the next tt iterations (tt is called the tabu tenure). Along with this rule, an aspiration criterion is applied to allow a move, regardless of its tabu status, if it leads to an improved best found solution. The exploration of the feasible search space terminates if no improvement is achieved in the consecutive N cons iterations (N cons is called the search depth). At this stage, the algorithm switches to the Infeasible Local Search (InfLS) phase that introduces a greater diversity, as the search is deemed to be trapped in a deep local optimum.

Aside from the diversification incurred during InfLS, FLS additionally incorporates a shake procedure similar to that used in [START_REF] Lai | Iterated variable neighborhood search for the capacitated clustering problem[END_REF]. The shake procedure is applied periodically, and consists in performing a random move from the combined feasible OneM ove and SwapM ove neighborhoods. This operation is repeated η times, where η is the shake strength.

Infeasible local search (InfLS)

The basic idea of InfLS is to relax the capacity constraint so as to allow the algorithm to visit some intermediate infeasible solutions

Evaluation function and neighborhood structures

To evaluate the quality of a solution s ∈ Ω during InfLS, we employ a penaltybased evaluation function f p which is a linear combination of the basic eval-10 uation function f (Equation. 7) and a penalty function associated with the degree of solution infeasibility:

f p (s) = p g=1 i,j∈Cg,i<j c ij -β × EX(s) (8)
where β is a self-adjustment penalty parameter that controls the degree of infeasibility introduced into the search. EX(s) is the total degree of infeasibility of s defined as EX(s) = p g=1 o g , where

o g =        L g -|C g |, if |C g | < L g |C g | -U g , if |C g | > U g 0, otherwise (9)
The InfLS phase employs the same three basic move operators (OneM ove, SwapM ove and 2-1 Exchange) defined in Section 2.3.2, but without any capacity restriction. We further use a fast incremental evaluation technique to effectively calculate the move gain that corresponds to the change in the penalty-based evaluation function f p . Specifically, for a given move denoted as mv (mv = OneM ove, SwapM ove or 2-1 Exchange), the move gain of mv can be defined as ∆ f (mv) = ∆ f (mv)β × ∆ EX (mv). In addition to the incremental matrix γ (see Section 2.3.2), another vector ω is maintained where each element ω g represents the total weight of nodes contained in cluster g. Since all moves induced by OneM ove, SwapM ove and 2-1 Exchange operators only involve two clusters, the weight of these two clusters after each move can be directly calculated by adding or subtracting the weight of the added or removed nodes. Furthermore, the move gain associated to this change can easily be updated as described in Equation. 9. The complexity of the above gain update procedure is O(1).

Exploration with InfLS

While both InfLS and FLS rely on the tabu search strategy, the main difference between these procedures lies in the previously described evaluation functions. The main scheme of InfLS is summarized in Algorithm 3.

During the search process, the variable β of the evaluation function (see Equation. [START_REF] Chen | A multi-facet survey on memetic computation[END_REF]) is periodically updated depending on the penalty counter (penalty count) which records the number of times a feasible solution has been found during λ consecutive iterations (λ is a parameter P OP ← P ool U pdating(S c , P ) /*Section 3.2 */ 11: end while Relying on the combined exploitation power of local optimization and exploration capacity of population-based search, Memetic Algorithm (MA) [START_REF] Moscato | A Gentle Introduction to Memetic Algorithms[END_REF] is an effective hybrid framework for tackling a variety of difficult combinatorial problems.

The main scheme of our MA for CCP is given in Algorithm 4. The algo-rithm consists of four basic components: a population initializing procedure, a crossover operator, the FITS procedure for local improvement and a population updating rule. Each solution from the initial population is obtained with the two randomized construction methods described in Section 2.2, and then further improved with our FITS method presented in Section 2. At each cycle (generation) of MA, two parent solutions are randomly selected from the population, and then recombined by means of the crossover operator to generate an offspring solution (Section 3.1). This new offspring is then improved by applying a fixed number of iterations of our FITS algorithm (Section 2.1). Finally, the population updating rule decides whether the improved offspring should be inserted into the population and which existing solution should be replaced (Section 3.2). This process is repeated until a predefined stopping condition (usually a fixed number of generations or time limit) is reached.

Cluster-based crossover

Algorithm 5 Cluster-based crossover Require: Two randomly selected parents

s 1 = {C 1 1 , C 1 2 , ..., C 1 p }, s 2 = {C 2 1 , C 2 2 , ..., C 2 p } Ensure: Offspring s o = {C o 1 , C o 2 , ..., C o p } 1: l ← 1 2: while l ≤ p do 3: if l%2 == 0 then 4:
Select a cluster C * ∈ s 1 with the maximum sum of edge weights Remove the subset C o l of nodes from s 1 and s 2 11:

l ← l + 1 12: end while 13: Assign in a greedy manner the unassigned nodes

V -{C o 1 ∪ C o 2 ∪ ... ∪ C o p }
Crossover operator is one of the key elements of a population-based algorithm.

It is well-known that a crossover's efficiency on a given optimization problem crucially depends on its ability to preserve pertinent properties ("building blocks") from parents to offspring [START_REF] Hao | Memetic Algorithms in Discrete Optimization[END_REF]. In the context of CCP, a building block may be defined as a cluster (i.e., group or subset of a graph partition), where the aim is to maximize the sum of edge weights whose two associated endpoints belong to the given cluster. As CCP can be classified as a grouping problem [START_REF] Falkenauer | Genetic algorithms and grouping problems[END_REF], it is more natural and straightforward to manipulate groups of objects (i.e., clusters) rather than individual objects. Furthermore, an analysis on a sample of locally optimal CCP solutions (see Section 5.3.1) discloses a high percentage of nodes that are always grouped together across high quality solutions, which provides a strong motivation for preserving the grouped nodes (i.e., cluster) from parent individuals to offspring solution. The proposed cluster-based crossover (CBX) is further inspired by the Greedy Partition Crossover (GPX) used for the classic graph coloring problem [START_REF] Galinier | Hybrid evolutionary algorithms for graph coloring[END_REF].

Given two parent solutions, CBX generates an offspring

s o = {C o 1 , C o 2 , .
.., C o p } in two sequential stages as summarized in Algorithm 5. The first stage performs p iterations (i.e., one iteration per cluster), where each iteration l consists in selecting a cluster C * from a reference parent such that the weighted sum of edges with both endpoints in C * is maximized. Cluster C * then becomes the l th building block of s o , followed by the removal of all the nodes contained in C * from both parent individuals. The reference parent is selected between s 1 and s 2 in an alternating manner. Note that the first crossover stage may result in a partial solution as some nodes may have been left unassigned. The second stage of the crossover process is a greedy construction method that consists in selecting an unassigned node v and inserting it into cluster g of the offspring solution, such that L g ≤ w v + |C o g | ≤ U g while maximizing the objective function value defined in Equation. 7. This process is repeated until all the nodes are assigned.

Population updating strategy

Population update strategy is another key element of a MA algorithm whose main role is to maintain a healthily diversified population throughout the search. To avoid premature convergence, we employ a quality-and-distance pool updating strategy which takes into account both the solution quality and the distance between individuals in the population to decide whether a new offspring should be introduced into the population. For this purpose, the distance Dist(S a , S b ) between two solutions S a and S b is defined as the minimum number of one-move steps required to transform S b to S a [START_REF] Lü | A hybrid metaheuristic approach to solving the ubqp problem[END_REF]. Given a population P = {S 1 , S 2 , ..., S |P | }, the distance between a solution S i (i ∈ 1, 2, ..., |P |) and P is computed as:

D S i ,P = min{Dist(S i , S j )|S j ∈ P, S j = S i } (10)
The main scheme of the population update procedure, which is similar to that used in [START_REF] Lü | A hybrid metaheuristic approach to solving the ubqp problem[END_REF][START_REF] Wu | A hybrid metaheuristic method for the maximum diversity problem[END_REF], is provided in Algorithm 6.

First, offspring S 0 is tentatively added to P resulting in P = P {S 0 }. The quality-and-distance score function R is then applied to rank each solution

S i ∈ P : R(S i , P ) = αX (f (S i )) + (1 -α)X (D S i ,P ) (11)
where α is a parameter set to 0.6 according to [START_REF] Lü | A hybrid metaheuristic approach to solving the ubqp problem[END_REF], and f (S i ) is the objective value of S i . X (.) is a function defined as:

X (y) = y-y min ymax-y min +1
where y max and y min are the maximum and the minimum possible values of y respectively. Finally, the solution S w with the smallest score is removed from the population to make space for offspring S 0 .

Computational experiments

This section provides an extensive assessment of the proposed algorithms on a well-known set of 133 benchmark instances from CCPLIB, as well as on two new groups of 50 large instances recently generated in [START_REF] Brimberg | Solving the capacitated clustering problem with variable neighborhood search[END_REF]. To evaluate the effectiveness of the proposed approaches, we perform comparisons with several state-of-the-art algorithms from the literature. For the ease of reading and for the sake of clarity, Table 1 shows the list of the proposed and the reference algorithms from the CCP literature. Calculate the distance between S i and P according to Eq.( 10)

4:

Calculate the goodness score R(S i , P ) of S i according to Eq.( 11) 5: end for 6: Identify the solution S w with the smallest goodness score in P :

S w = min {R(S i , P )|i = 0, ..., |P |} 7: if S w = S 0 then 8:
Replace S w with S 0 : P = P {S 0 } \ {S w } 9: end if

Benchmark instances

The instances from CCPLIB1 can be grouped into three sets:

• DB (10 instances): This set was introduced by Deng and Bard [START_REF] Deng | A reactive GRASP with path relinking for capacitated clustering[END_REF] for the Maximally Diverse Grouping Problem (MDGP), and adapted for CCP in [START_REF] Martínez-Gavara | Tabu search and GRASP for the capacitated clustering problem[END_REF] by generating node weights with a uniform distribution in the range [0, 10]. These instances are characterized by n = 82, p = 8, L g = 25, U g = 75.

• RanReal (40 instances): This set was originally proposed in [START_REF] Gallego | Tabu search with strategic oscillation for the maximally diverse grouping problem[END_REF] and first adapted for CCP in [START_REF] Martínez-Gavara | Tabu search and GRASP for the capacitated clustering problem[END_REF] [START_REF] Morán-Mirabal | Randomized heuristics for handover minimization in mobility networks[END_REF] with n ∈ {20, 30, 40, 100, 200, 400} and p ∈ {5, 10, 15, 25, 50}. The edge weights are real numbers, while L g and U g are respectively set to 0 and a real number that differs for each instance. The set is widely used in the literature for the handover minimization problem.

Additionally, we perform experiments on the following two groups of 50 large instances2 from [7]:

• RanReal960 (30 instances): This set consists of 3 subsets. Each subset contains 10 instances characterized as follows:

• n = 960, p = 30, L g = 120, U g = 180;

• n = 960, p = 40, L g = 90, U g = 135;

• n = 960, p = 60, L g = 60, U g = 90; • MDG (20 instances): This set includes 20 instances with n = 2000, p = 50, L g = 200, U g = 300. The edge weights and the node weights are generated using the same method as in [START_REF] Duarte | Tabu search and grasp for the maximum diversity problem[END_REF].

Experimental protocol

The proposed FITS algorithm requires eight parameters: tabu tenure (tt), search depth of FLS (N cons ), shake frequency (δ) and shake strength (η) of FLS, update frequency (λ) of the self-adjustment penalty β in InfLS, update coefficients of β (τ , µ 1 , µ 2 ) in InfLS, and the maximum number of InfLS iterations (M ). MA requires two additional parameters: the population size (|P |) and the number of FITS iterations (Iter). To determine the appropriate parameter settings for FITS and MA, we run the Iterated F-race (IFR) method [START_REF] Bartz-Beielstein | Experimental Methods for the Analysis of Optimization Algorithms[END_REF], implemented within the IRACE package [START_REF] López-Ibánez | The irace package: iterated racing for automatic algorithm configuration[END_REF], on a selection of 20 RanReal instances with n = 240 and n = 480. The tuning budgets of FITS and MA are set to 500 runs with the time limit of 1.0 * n seconds for each run. Table 2 shows the tested and the final values obtained in the tuning process. For λ, τ µ 1 and µ 2 , we simply adopt the values recommended in [START_REF] Chen | A hybrid metaheuristic approach for the capacitated arc routing problem[END_REF][START_REF] Hertz | A tabu search heuristic for the capacitated arc routing problem[END_REF].

The proposed algorithms are coded in C++, compiled with the g++ compiler using option "-O3", and executed on an Intel E5-2670 processor (2.8GHz) with 2GB RAM running under Linux. For time scaling purposes, the execution time of the DIMACS machine benchmark3 on our system is 0.19s for graph r300.5, 1.17s for graph r400.5 and 4.54s for graph r500.5.

Comparison between FITS and the state-of-the-art algorithms on the general CCP instances

To evaluate the performance of FITS on the first two sets of instances (DB and RanReal), we provide comparisons with several state-of-the-art algorithms including GRASP [START_REF] Martínez-Gavara | Tabu search and GRASP for the capacitated clustering problem[END_REF], TS [START_REF] Martínez-Gavara | Tabu search and GRASP for the capacitated clustering problem[END_REF], GRASP+TS [START_REF] Martínez-Gavara | Tabu search and GRASP for the capacitated clustering problem[END_REF], IVNS [START_REF] Lai | Iterated variable neighborhood search for the capacitated clustering problem[END_REF], GVNS [START_REF] Brimberg | Solving the capacitated clustering problem with variable neighborhood search[END_REF] and SGVNS [START_REF] Brimberg | Solving the capacitated clustering problem with variable neighborhood search[END_REF]. For a fair comparison with these state-of-the-art algorithms, we use their corresponding source codes and run them on our computing platform under the same computing conditions as described in Section 4.2. GRASP, TS, GRASP+TS, IVNS were previously re-implemented by Lai et al. [START_REF] Lai | Iterated variable neighborhood search for the capacitated clustering problem[END_REF] and their source codes are available at http://www.info.univ-angers.fr/ pub/hao/ccp.html. The source codes of GVNS and SGVNS were shared by the corresponding author and are available at http://www.mi.sanu.ac.rs/ ~nenad/ccp/. All the codes were compiled using g++ compiler with the '-O3' option. For all the algorithms in this experiment, the stopping condition is a fixed cutoff time limit t max set to 1.0 * n seconds, where n is the number of nodes in the given graph. Due to the stochastic nature of the algorithms, we perform 20 independent runs of each algorithm per instance. Except for the comparison with SGVNS, the statistical test reveals a significant difference in performance between each of the reference algorithms (p-value ≤ 0.05), demonstrating the efficiency of FITS on the DB and the RanReal instances.

Comparison between FITS and the state-of-the-art algorithms on the large CCP instances

RanReal960 and MDG benchmarks consist of large instances recently used to assess the performance of the algorithms for CCP in [START_REF] Brimberg | Solving the capacitated clustering problem with variable neighborhood search[END_REF]. A summary of the statistical results for these instances, obtained with FITS and the reference algorithms, are shown in Table 4. Detailed results are given in Tables 13 and14 of the Appendix.

Row p-value best /p-value avg reveals a statistically significant difference in performance between FITS and all the reference algorithms except GVNS. When considering the average percent deviation from the best-found solutions within the experiments, SGVNS outperforms all the algorithms with Dev best = 0.01% and Dev avg = 0.61%, while FITS exhibits a better performance than GRASP, TS, GRASP+TS and IVNS.

Comparison between FITS and the state-of-the-art algorithms on the handover minimization instances

Table 5 summarizes the statistical results reported with FITS and the six stateof-the-art algorithms on the handover minimization instances with n ≥ 100. These reference algorithms include IVNS [START_REF] Lai | Iterated variable neighborhood search for the capacitated clustering problem[END_REF], GVNS [START_REF] Brimberg | Solving the capacitated clustering problem with variable neighborhood search[END_REF], SGVNS [START_REF] Brimberg | Solving the capacitated clustering problem with variable neighborhood search[END_REF] and three algorithms proposed in [START_REF] Morán-Mirabal | Randomized heuristics for handover minimization in mobility networks[END_REF]: GevPR-HMP, GQAP and BRKGA. The handover minimization instances with n ∈ {20, 30, 40} do not appear to be challenging as all the considered algorithms are able to attain the best-known solution within a very short computing time for each case. For completeness, the computational results for these instances are provided in Appendix (Table 22).

In Table 5, the results reported with GevPR-HMP, GQAP and BRKGA are directly compiled from [START_REF] Morán-Mirabal | Randomized heuristics for handover minimization in mobility networks[END_REF], and were obtained over 5 independent runs with a cutoff time set to 24 hours. This is significantly longer than the time limit used for IVNS, GVNS, SGVNS and FITS, which is set to 1.0 * n seconds (24 hours vs. n ≤ 400 seconds). The results for IVNS, GVNS, SGVNS and FITS were obtained across 20 independent runs under the same computing conditions as described in Section 4.2. When handling the handover minimization instances, we use the results in the form of minimization for a direct comparison. The following relation is used to transform the CCP objective function into the equivalent objective for handover minimization:

f min = 2( i<j c ij -f max ),
where f min and f max represent the objective values of handover minimization and CCP respectively. The symbol "-" denotes the cases when the result is not reported in the literature. Detailed results are given in Tables 15 and16 of the Appendix.

From Table 5, we observe that GevPR-HMP, GQAP, BRKGA, IVNS, GVNS, SGVNS and FITS outperform the other reference algorithms on 9, 11, 2, 37, 23, 24 and 43 instances respectively in terms of the best objective value.

In terms of the average results, FITS achieves a better performance on 30 instances, while GevPR-HMP, GQAP, BRKGA, IVNS, GVNS and SGVNS outperform the other methods on 4, 0, 0, 18, 15, 25 instances respectively. Notice that if we sum up the number of times each algorithm performed best, we get 149 for Best and 92 Avg respectively. This is because several algorithms obtain the same best objective value and the same average objective value on some instances. In terms of the average percent deviation of the best/average results from the best solutions obtained within this experiment, FITS shows the best performance with Dev best = 0.00% and Dev avg = 0.20%. Finally, the p-values of the Wilcoxon pairwise test (row 'p-value best /p-value avg ') show a statistically significant difference in the best performance between FITS and the six reference approaches with p-value ≤ 0.05, which shows the benefit of FITS on the handover minimization instances with n ≥ 100.

Time-to-target analysis

To further compare the performance between FITS and the reference algorithms, we apply the time-to-target (TTT) analysis which identifies the empirical probability distribution of the time required to achieve a given target value [START_REF] Aiex | Ttt plots: a perl program to create time-to-target plots[END_REF]. We conduct this TTT experiment by executing 100 independent runs of GRASP, TS, GRASP+TS, IVNS, GVNS, SGVNS and FITS on each instance. For each instance/target pair, the running times are sorted in an increasing order. We associate with the i-th sorted running time t i a probability p i = (i-0.5)/100, and plot the points (t i , p i ). In this experiment, in order to allow all the algorithms to reach the target in all runs, the target value is set to be a value slightly smaller than the best obtained objective value. Fig. 1

Comparison between FITS and MA

In this section, we additionally compare FITS and MA on the first two sets (DB and RanReal) of CCP instances and the handover minimization instances with n ≥ 100. For this comparison, each algorithm is executed 20 times per Table 6 summarizes the statistical comparative results between FITS and MA, while detailed results are given in Tables 17 and18 in the Appendix. In a nutshell, MA improves on the best result reported by FITS for 28 instances, and fails to match the best solution obtained with FITS for only three instances. In terms of the average performance, MA outperforms FITS on 52 instances, and is outperformed by FITS on 15 instances. When considering the average percent deviation of the best/average results from the best solutions found within this experiment (Average Dev best /Dev avg ), MA achieves a better performance than FITS on all the three sets of instances. Finally, for all the three sets of instances, the Wilcoxon test indicates a statistically significant difference in the best and the average performances with p-value = 6.437e-4 and p-value = 0.049 respectively. Thus, we can conclude that MA should be considered over FITS given a longer time limit.

Analysis

This section evaluates the importance of the key elements of the proposed FITS and MA algorithms: (i) the joint exploitation of feasible and infeasible search space, (ii) the best improvement strategy vs. the first improvement strategy, and (iii) the cluster-based crossover operator. The experiments presented below are carried out on a set of 40 instances from the RanReal benchmark with 20 independent executions per instance under the same computing platform as described in Section 4.2.

Analysis of the combined exploitation of feasible and infeasible search space

While a number of existing heuristics for CCP including TS [START_REF] Martínez-Gavara | Tabu search and GRASP for the capacitated clustering problem[END_REF], IVNS [START_REF] Lai | Iterated variable neighborhood search for the capacitated clustering problem[END_REF], GVNS [START_REF] Brimberg | Solving the capacitated clustering problem with variable neighborhood search[END_REF] and SGVNS [START_REF] Brimberg | Solving the capacitated clustering problem with variable neighborhood search[END_REF] restrict their search to the feasible regions only, a key feature of FITS and several other heuristics like GRASP-PR [START_REF] Deng | A reactive GRASP with path relinking for capacitated clustering[END_REF], TS SO [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhood search[END_REF], GevPR-HMP [START_REF] Morán-Mirabal | Randomized heuristics for handover minimization in mobility networks[END_REF] and GQAP [START_REF] Morán-Mirabal | Randomized heuristics for handover minimization in mobility networks[END_REF] is the consideration of infeasible solutions. By alternating between feasible and infeasible local searches, FITS is able to visit various zones of the search space without being easily trapped in a local optimum. To evaluate the effectiveness of this hybrid scheme, we compare FITS with its two underlying components, namely the feasible local search (FLS) and the infeasible local search (InfLS).

Table 7 summarizes the statistical results of this analysis, while detailed results are given in the Appendix (Table 19). As observed in Table 7, it is evident that FITS outperforms both FLS and InfLS. More precisely, in terms of the best and the average performance, FITS respectively reports matching or better results than both algorithms on 35 and 36 out of the 40 instances with pvalues of 6.23e-5/4.20e-7 and 2.54e-10/2.54e-10. We further observe that InfLS alone appears to be the weakest of the three versions on all of the tested instances. To complement this comparison, Fig. 2 plots the performances of the three algorithms on these 40 instances. For each instance and each algorithm, the y-axis shows the percent deviation of the average result from the best 20 in the Appendix.

From Table 8, one notices that FITS outperforms FITS FI on all of the instances both in terms of the best and the average result (#Best/Avg.). In terms of the average percent deviation of the best/average results from the best solutions obtained within this experiment (Average Dev best /Dev avg ), FITS shows a better performance than FITS FI with p-value = 2.54e-10 and p-value = 2.54e-10 respectively. These observations confirm the usefulness of the best improvement strategy within the proposed tabu search framework. However, as expected, FITS FI shows to be faster than FITS in terms of the average time required to reach the final solution.

Analysis of the cluster-based crossover

Motivation behind the cluster-based crossover

As explained in Section 3.1, the basic idea behind the cluster-based crossover is to preserve building blocks (clusters of nodes) from parent individuals to offspring solution. Such crossovers have shown to be effective when there is a large percentage of nodes that are always grouped together between highquality local optima (including global optima which are technically speaking also local optima) [START_REF] Benlic | A multilevel memetic approach for improving graph k-partitions[END_REF].

Given two local optima

s 1 = {C 1 1 , C 1 2 , ..., C 1 p } and s 2 = {C 2 1 , C 2 2 , ..., C 2 p }, let E = {(C 1 i , C 2 j )|i ∈ {1, 2, .
., p}, j ∈ {1, 2, .., p}} denote the set of all the p × p cluster combinations of s 1 and s 2 , and let J denote the set of nodes that are grouped together in both s 1 and s 2 . Starting from J = ∅, we use an iterative procedure that determines the largest percentage of nodes grouped together in s 1 and s 2 with the following steps: (i) for each cluster combination (C

1 i , C 2 j ) ∈ E , compute the number of identical nodes k C 1 i C 2 j = C 1 i ∩ C 2 j ; (ii) select a combination (C 1 i , C 2 j ) ∈ E with the largest k C 1 i C 2
j and place the common nodes into J (i.e., J = J ∪ {C 1 i ∩ C 2 j }); (iii) remove from E all combinations associated with C 1 i and C 2 j . This process is repeated until E becomes empty. The percentage of nodes that are grouped together in both s 1 and s 2 is then expressed as 100% * |J| |V | .

For this analysis, we employ a selection of 10 hard instances from the Ran-Real set. For each instance, we collect a set S all of local optima of different qualities, obtained after 500 independent runs of MA and FITS with different time limits. We select the top 20% (100) local optima with the largest objective values from S all to form the subset S hq of 'high-quality solutions'. Similarly, we take the bottom 20% (100) with the smallest objective values to form the subset S lo of 'low-quality solutions'. Table 9 shows the percentages of common node groupings across all the local optima in S hq , S all and S lo respectively. From these results, we conclude that the percentage of nodes that are grouped together throughout each of the high-quality local optima from S hq is very large, ranging from 72% to 92%. Assuming that high-quality solutions might be close to an optimal solution or could themselves constitute optimal solutions, it is very likely that these clusters might form the building blocks of a global optimum.

Comparison with the uniform crossover

To evaluate the benefit of the cluster-based crossover that transfers pertinent properties from parents to offspring, we compare it with a standard uniform crossover where each node is randomly transferred to the same cluster as in one of the two parents. As the resulting offspring may violate the capacity constraint, the uniform crossover proceeds next by moving a randomly selected node from the highest to the lowest weight cluster until the solution feasibility is reached. For this experiment, the other MA components are left unchanged. The time limit is set to 15.0 * n seconds per execution.

Table 10 summarizes the statistical results for each MA version on the 40 RanReal instances, while detailed results are given in Table 21 in the Appendix. Although the cluster-based crossover results in longer computing time than the uniform crossover due to a higher complexity, it clearly outperforms the uniform crossover in terms of both the best and the average results with p-value = 4.46e-7 and p-value = 2.54e-10 respectively. Indeed, out of the 40 instances, the cluster-based crossover is outperformed by the uniform crossover only on two instances. These observations highlight the importance of preserving important building blocks from parents to children in case of CCP.

Conclusion

We presented two highly effective heuristics for the Capacitated Clustering Problem (CCP): a tabu search approach (denoted as FITS) that alternates between exploration in feasible and infeasible search space regions, and a Memetic Algorithm (MA) that extends FITS with a dedicated cluster-based crossover and a quality-and-distance pool updating strategy. The computational results on five sets of 183 CCP instances indicate that both FITS and MA compete favorably with the current state-of-the-art algorithms. The investigation of several essential components of the proposed algorithms sheds light on the following points. First, the consideration of both feasible and infeasible search space regions can greatly enhance the neighborhood search for CCP. Second, the best improvement strategy is able to outperform the first improvement strategy within the tabu search framework for CCP. Third, given an extended time limit, MA can further improve upon the performance of FITS which is greatly due to the cluster-based crossover that transfers pertinent properties from parents to offspring. The use of this crossover within MA was motivated by a large degree of similarity between high-quality CCP solutions. Finally, this work demonstrates the effectiveness of exploring both feasible and infeasible spaces for CCP -an idea that certainly deserves to be investigated on other highly constrained problems in the future.

7 :

 7 Select a cluster C * ∈ s 2 with the maximum sum of edge weights

Algorithm 6 1 :

 61 Population update strategy Require: Offspring S 0 , population P = S 1 , S 2 , ..., S |P | Ensure: Updated population P = S 1 , S 2 , ..., S |P | Tentatively add S 0 to P : P = P {S 0 } 2: for i = 0, 1, ..., |P | do 3:

  by generating node weights with a uniform distribution in the range [0, 10]. It includes 20 instances with n = 240, p = 12, L g = 75, U g = 125, and another 20 instances with n = 480, p = 20, L g = 100, U g = 150. The edge weights are real numbers randomly generated in the range [0, 100]. • MM (83 instances): This set was introduced by Morán-Mirabal et al.

Fig. 1 .

 1 Fig. 1. Probability distribution of the time required to achieve a target value.

Fig. 2 .

 2 Fig. 2. Percentage deviation of the average result reported with FITS, FLS and ILS from the best solutions found in this experiment for the 40 RanReal instances.

  Feasible solution s returned by FLS Ensure: Best feasible solution found during InfLS s local best , final feasible solution s f inal 1: s tmp ← s /*s tmp is a duplicate of the starting feasible solution s*/

	Algorithm 3 Infeasible Local Search (InfLS)
	Require: 2: f lag f s ← f alse /*f lag f s is a boolean variable that indicates whether a feasible solution was encountered during InfLS*/
	3: s local best ← s 4: penalty count ← 0 5: penalty f actor ← 2 6: M I ← 0 7: Initialize tabu list
	. In this way, a larger Choose the best allowed move m ∈ {OneM ove ∪ SwapM ove ∪ 2-1 8: while M I ≤ M do 9: Exchange}
	number of moves become available, enabling transitions between structurally different high-quality feasible solutions. 10: s ← s ⊕ m /*Perform the best move*/ 11: Update tabu list
	12:	if s is a feasible solution then
	13:	if f (s) > f (s local best ) then
	Algorithm 2 Feasible Local Search 14: s local best ← s 15: end if Require: Initial solution s 16: else Ensure: Final solution s, best solution s local best found during this phase 1: s local best ← s 17: penalty count ← penalty count + 1 18: end if 2: N I ← 0 /*number of consecutive iterations without improvement of s local best */ 19: M I ← M I + 1 20: if (M I + 1)%λ==0 then 3: Iter1 ← 0 /*iteration counter*/ 21: if penalty count > µ 1 then 4: Initialize tabu list 5: while N I < N cons do 6: 22: penalty f actor ← penalty f actor * τ 23: else if penalty count < µ 2 then Choose the best allowed move m ∈ {OneM ove ∪ SwapM ove ∪ 2-1 Exchange} 7: 24: penalty f actor ← penalty f actor/τ 25: end if s ← s ⊕ m /*Perform the best move*/ 8: Update tabu list 26: penalty count ← 0 27: end if
	9: 28:	if f (s) > f (s local best ) then if s is a feasible solution then
	10: 11: 12: 29: 30: 31:	s local best ← s N I ← 0 else s f inal ← s f lag f s ← true end if
	13: 32: end while N I ← N I + 1 14: end if 33: if f lag f s == f alse then
	15: 16: 17: 34: 35: 36: end if if (Iter1 + 1)%δ==0 then Shake() end if s tmp ← Shake() /*apply the shake procedure to s tmp , Section 2.3.3*/ s f inal ← s tmp
	18: 19: end while Iter1 ← Iter1 + 1 37: return (s f inal ,s local best )
	20: return (s, s local best )

  ). Recall that β controls the degree of infeasibility introduced into the search. More precisely, β is increased by a multiple of τ if penalty count > µ 1 , and is divided by τ if penalty count < µ 2 (τ , µ 1 and µ 2 are parameters). Furthermore, the best found feasible solution s local best is updated with the current solution s if an improved feasible solution has been found. InfLS terminates after M iterations (M is a parameter), followed by the FLS phase. The starting point for the next round of the FLS search is the most recently encountered feasible solution s f inal returned by InfLS.

	Finally, if no feasible solution is found during the InfLS process, the starting
	feasible solution s returned by FLS is perturbed with the shake procedure
	(Section 2.3.3) and is then returned as the output of InfLS. Note that the
	shake procedure explores feasible solutions only thus resulting in a feasible
	solution. Algorithm 3 summarizes the general procedure of the InfLS phase,
	in which we use a boolean variable f lag f s to indicate whether a feasible
	solution is encountered during the InfLS phase. It is important to mention
	that FLS and InfLS use two separate tabu lists, which are critical to the
	performance of the two local search phases as they prevent the search from
	short-term cycling. When switching back to FLS, the tabu list is re-initialized
	before entering the main loop.
	3 Memetic algorithm
	Algorithm 4 Main scheme of MA
	Require: Graph G = (V, E)
	Ensure: Best solution S * found so far
	1: Initialize population P = {S 1 , S 2 , ..., S |P | } 2: S * ← Best(P ) 3: while Time does not exceed t max do
	4: 5: 6: 7:	Randomly select two parent solutions S i ∈ P and S j ∈ P S c ← Crossover(S i , S j ) /*Section 3.1*/ S c ← F IT S (S c ) /*Section 2.1*/ if f (S c ) > f (S * ) then
	8: 9:	S * ← S c end if
	10:

Table 1

 1 List of the reference algorithms for CCP and its equivalent HMP.

	Algorithm	Reference	Search strategy
	name
	FITS	-	A tabu search alternating between exploration in feasible and infeasible search space
	MA	-	A memetic algorithm that extends FITS with a dedicated cluster-based crossover
			operator
	GRASP-PR	[10](2011)	A reactive GRASP
	GevPR-HMP	[31](2013)	A GRASP combined with an evolutionary path-relinking algorithm in which a repair
			procedure is applied to achieve feasibility
	GQAP	[31](2013)	A GRASP method with a new variant of path-relinking dealing with infeasibilities
	BRKGA	[31](2013)	A biased random-key genetic algorithm using a parameterized uniform crossover
	GRASP	[29](2015)	A simplified GRASP
	TS	[29](2015)	A tabu search algorithm exploiting the 2-1 exchange neighborhood
	GRASP+TS	[29](2015)	A hybrid combining GRASP with tabu search
	TS SO	[29](2015)	A tabu search with strategic oscillation that considers infeasible solutions
	IVNS	[24](2016)	An iterated variable neighborhood search combining an extended variable neighbor-
			hood descent with a randomized shake procedure
	GRASP2-1	[30](2017)	A new GRASP method in which the improvement procedure performs 2-1 exchanges
	IG	[30](2017)	An iterated greedy method alternating between destructive and constructive phases
	IG-GRASP	[30](2017)	A hybrid between GRASP2-1 and iterated greedy method
	GVNS	[7](2017)	A general variable neighborhood search that follows the standard VNS approach
			including more levels of shaking
	SGVNS	[7](2017)	A skewed general variable neighborhood search that allows moves to inferior solutions

Table 2

 2 Settings of the parameters.

	Parameter Section	Description	Considered values	Final
					value
	Ncons	2.3	search depth of each FLS phase	{500, 700, 1000, 1500, 2000}	1000
	M	2.4	maximum number of iterations of each InfLS phase	{100, 150, 200, 250, 300}	200
	tt	2.3,2.4	tabu tenure	{5, 7, 10, 12, 15}	10
	δ	2.3	frequency of shake	{300, 400, 500, 600, 700}	500
	η	2.3	shake strength	0.14*n} {0.06*n, 0.08*n, 0.10*n, 0.12*n,	0.10*n
	Iter	3	number of iterations of FITS in MA	{5000, 8000, 10000, 12000, 15000}	10000
	|P |	3	size of population	{5, 7, 10, 13, 15}	5
	λ	2.4	update frequency of β in InfLS	-	5
	τ	2.4	update coefficients of β in InfLS	-	2
	µ 1	2.4	update coefficients of β in InfLS	-	4
	µ 2	2.4	update coefficients of β in InfLS	-	1

Table 3

 3 Statistical results for FITS and six state-of-the-art algorithms on two sets of CCP instances: DB and RanReal. The best performance is indicated in bold.

	Instance		GRASP	TS	GRASP+TS IVNS	GVNS	SGVNS FITS
	set
	DB	&	#Best/Avg.	8/2	4/0	9/2	14/10	10/10	32/40	29/20
	RanReal
			valueavg p-value best /p-	11/1.54e-7.75e-	11/1.54e-7.75e-	9/1.54e-10 2.86e-	7/7.70e-1.25e-	9/2.54e-1.87e-	0.88/0.02
				12	12		8	10
			Average	6.66/7.56 1.33/2.42	1.35/2.00	0.26/0.40	0.22/0.63 0.20/0.36 0.19/0.33
			Dev best /Devavg (%)
			AvgTime(s)	149.25	50.68	187.26	224.24	243.04	187.33	201.35

Table 3

 3 summarizes the statistical results for each algorithm on the instances of the DB and RanReal sets. Row '#Best/Avg.' indicates the number of cases that each algorithm outperforms the remaining approaches in terms of the best and the average objective value. The average percent deviation of the best/average result from the best solution obtained within this experiment is provided in row 'Average Dev best /Dev avg ', while the average computing time (in seconds) required by each algorithm to reach its final objective value is provided in row 'AvgTime'. For each instance, we calculate the best and the average deviation (Dev best and Dev avg ) as (f *f )/f * × 100, where f is the best or the average result and f * is the best solution obtained with all the compared algorithms. Finally, to determine whether there exists a statistically significant difference in performance between FITS and the six reference algorithms, row 'p-value best /p-value avg ' provides the p-values obtained with the pairwise Wilcoxon statistical test on the best/average results. Detailed results on the DB and the RanReal benchmarks are given in the Appendix (Tables

Table 4

 4 Statistical results of FITS and six state-of-the-art algorithms on two large instance sets: RanReal960 and MDG. The best performance is given in bold.

	Instance		GRASP	TS	GRASP+TS IVNS	GVNS	SGVNS	FITS
	set
	RanReal960	#Best/Avg.	0/0	0/0	0/0	0/0	0/0	46/46	4/4
	& MDG
		p-value best /p-valueavg	1.54e-12/1.54e-	1.54e-12/1.54e-	1.54e-12/1.54e-	1.14e-11/1.54e-	0.09/0.16 2.86e-9/2.86e-
			12	12	12	12		9
		Average	16.94/17.44 6.60/7.36 3.77/4.47	1.27/1.70	0.61/0.80 0.01/0.61 0.89/1.07
		Dev best /Devavg (%)
		AvgTime(s)	694.09	1115.02	1221.46	1251.36	1202.49	1240.53	1282.67
	11, 12).
	From

Table 3 ,

 3 we observe that GRASP, TS, GRASP+TS, IVNS, GVNS, SGVNS and FITS respectively outperform the other algorithms on 8, 4, 9, 14, 10, 32 and 29 instances in terms of the best objective value. In terms of the average results, FITS achieves better performance on 20 instances, while GRASP, TS, GRASP+TS, IVNS, GVNS and SGVNS outperform the other methods on 2, 0, 2, 10, 10 and 40 instances respectively. In terms of the average percent deviation (Average Dev best /Dev avg ), FITS reports the smallest deviation from the best solutions obtained within this experiment (0.19%/0.33%).

Table 5

 5 Statistical results of FITS and the six state-of-the-art algorithms on handover minimization instances with n ≥ 100. The best performance is indicated in bold.

	Instance	GevPR-	GQAP	BRKGA	IVNS	GVNS	SGVNS	FITS
	set		HMP
	MM	#Best/Avg.	9/4	11/0	2/0	37/18	23/15	24/25	43/30
	(n 100)	≥
		p-value best /p-valueavg	1.97e-9/1.68e-	5.51e-9/1.97e-	5.47e-11/1.97e-	0.00/0.02	1.19e-5/0.08	7.44e-5/0.38
			8	11	11
		Average	1.32/1.74	8.37/10.15 3.19/3.94	0.01/0.28	0.20/1.53 0.19/0.33 0.00/0.20
		Dev best /Devavg (%)
		AvgTime(s)	-	-	-	97.82	66.18	75.77	120.48

Table 6

 6 Comparative statistical results between FITS and MA on the first two sets (DB and RanReal) of CCP instances and the handover minimization instances with n ≥ 100. The better performances are indicated in bold.

	Instance set		FITS	MA
	DB & RanReal	#Better1/Better2	3/7	25/33
		Average Dev best /Devavg (%)	0.03/0.11	0.00/0.10
	MM (n ≥ 100)	#Better1/Better2	0/8	3/19
		Average Dev best /Devavg (%)	0.00/0.15	0.00/0.05

Table 7

 7 Summary of statistical results obtained with FITS and its two underlying components FLS and InfLS for the RanReal benchmark set. The best performance is indicated in bold.

	Instance set		FITS	FLS	InfLS
	RanReal	#Best/Avg.	35/36	12/4	0/0
		p-value best /p-valueavg		6.23e-5/4.20e-7	2.54e-10/2.54e-10
		Average Dev best /Devavg (%)	0.01/0.18	0.03/0.21	0.33/0.53
	instance with a prolonged computing time of 15.0 * n seconds. The reason behind an extended cutoff time is due to our experiences and observations
	from previous studies [19] indicating a slower convergence pace of an MA
	compared to a local search algorithm.

Table 8

 8 Summary of statistical results obtained with FITS and its variant FITS FI on the 40 RanReal instances. The best performance is indicated in bold. The best improvement strategy v.s. the first improvement strategy As described in Section 2.1, FITS applies the best improvement strategy to select a solution from the neighborhoods induced by the three move operators, i.e., OneM ove, SwapM ove and 2-1 Exchange. To verify the importance of this feature, we create a variant of FITS (denoted by FITS FI) that uses the first improvement strategy, where the earliest visited neighboring solution of improved quality replaces the current solution. The summary of statistical results of this comparison on the 40 RanReal instances are shown in Table8, while detailed results are given in Table

	Instance set		FITS	FITS FI
	RanReal	#Best/Avg.	40/40	0/0
		p-value best /p-valueavg		2.54e-10/2.54e-10
		Average Dev best /Devavg (%)	0.00/0.15	0.56/0.71
		AvgTime(s)	249.39	157.44
	solutions found in this experiment. The figure further highlights the benefit of
	the combined use of the feasible and the infeasible local search.
	5.2

Table 9

 9 Percentage of nodes that are grouped together in local optima of various qualities

	Instance	S hq	S all	S lo	Instance	S hq	S all	S lo
	RanReal240 01	72	54	45	RanReal480 01	82	29	14
	RanReal240 02	85	49	39	RanReal480 02	88	26	14
	RanReal240 03	75	48	38	RanReal480 03	88	28	14
	RanReal240 04	80	48	41	RanReal480 04	92	29	14
	RanReal240 05	86	48	38	RanReal480 05	83	26	14

Table 10

 10 Comparison of the two MA versions using a uniform and a cluster-based crossover respectively on the 40 RanReal instances. The best performance is indicated in bold.

	Instance set		Uniform CX	Cluster-based CX
	RanReal	#Best/Avg.	2/0	31/40
		p-value best /p-valueavg	4.46e-7/2.54e-10
		Average Dev best /Devavg (%)	0.56/0.71	0.00/0.15
		AvgTime(s)	2898.81	3856.70

Table 12 .

 12 Comparison between FITS and six state-of-the-art algorithms on the first two sets of CCP instances (RanReal and DB).

	522218.13 521267.22
	519720.44
	521442.89
	509272.24
	508273.65
	469860.30

Table 13 .

 13 Comparison of FITS with six state-of-the-art algorithms on the large CCP instances (RanReal960 and MDG).

	best
	f

Table 14 .

 14 Comparison between FITS and the six state-of-the-art algorithms on the large CCP instances (RanReal960 and MDG).

	best /Devavg (%)
	Dev

Table 15 .

 15 Comparison of FITS with the six state-of-art algorithms on handover minimization instances with n ≥ 100, where the objective is to minimize the total cost function.

		favg
	FITS
		best
		f
		favg
	SGVNS
		f best
		favg
	GVNS
		best
		f
		favg
	IVNS	best
		f
	BRKGA	f best favg
		favg
	GQAP	f best
	GevPR-HMP	best favg
		f
	Name

Table 16 .

 16 Comparison between FITS and six state-of-the-art algorithms on the handover minimization instances with n ≥ 100.Table 17Comparison between FITS and MA on two sets of CCP instances (RanReal and DB).

	tavg	IVNS GVNS SGVNS FITS		1.03 0.24 0.28 3.84	1.54 0.39 0.49 25.06	0.08 0.06 0.10 0.14	0.30 0.13 0.11 2.03	1.62 0.24 0.36 11.35	1.54 0.64 0.82 4.38	1.62 1.33 0.66 5.57	5.42 1.40 0.86 12.78	6.18 0.13 0.21 19.05	6.89 0.15 0.23 14.61	9.19 0.45 0.48 25.75	10.47 0.87 0.71 7.41	7.61 0.68 0.59 8.70	5.63 1.26 0.05 14.18	9.45 1.42 0.43 42.71	25.19 2.23 2.77 67.35	70.30 20.16 55.81 116.10	38.27 2.08 1.85 13.10	48.99 13.94 20.44 43.07	11.33 29.97 31.91 40.97	41.78 24.60 54.73 83.74	87.71 47.07 19.05 91.89	53.89 86.94 51.56 103.30	49.99 35.98 19.00 92.27	9.64 6.86 5.17 78.63	76.15 62.88 65.42 87.34	69.51 25.08 14.92 105.89	69.87 46.94 83.88 101.50	81.66 64.90 55.80 111.74	86.92 65.05 68.95 131.04	184.94 114.10 154.72 218.52	197.33 195.42 162.95 230.06	156.99 73.55 125.33 253.96	245.49 123.96 171.49 226.90	234.97 136.54 243.64 271.41	239.15 124.71 127.27 280.77	206.04 109.35 114.10 263.82	213.34 176.03 182.30 259.13	209.56 108.17 142.02 274.96	243.63 192.24 215.43 295.26	260.81 214.55 244.55 255.22	284.52 228.68 229.22 280.60	302.17 197.58 258.58 267.40	265.09 228.35 230.35 294.04	267.93 210.97 249.88 284.11	97.82 66.18 75.77 120.48
		FITS		0.00/0.00	0.00/0.00	0.00/0.00	0.00/0.00	0.00/0.00	0.00/0.00	0.00/0.00	0.00/0.00	0.00/0.00	0.00/0.00	0.00/0.00	0.00/0.00	0.00/0.00	0.00/0.00	0.00/0.26	0.00/0.00	0.00/0.79	0.00/0.01	0.00/0.07	0.00/0.36	0.00/0.00	0.00/0.05	0.00/0.06	0.00/0.00	0.00/0.00	0.00/0.16	0.00/0.07	0.00/0.22	0.00/0.11	0.00/0.03	0.00/0.73	0.00/0.91	0.00/0.62	0.00/1.81	0.00/0.62	0.06/0.42	0.00/0.11	0.00/0.29	0.00/0.12	0.00/0.18	0.00/0.11	0.00/0.21	0.00/0.15	0.02/0.44	0.00/0.30	0.00/0.20
	best /Devavg (%)	IVNS GVNS SGVNS		0.00/0.00 0.00/0.00 0.00/0.00	0.00/0.00 0.00/0.00 0.00/0.00	0.00/0.00 4.30/4.30 4.30/4.30	0.00/0.00 0.00/0.00 0.00/0.00	0.00/0.00 0.65/0.65 0.65/0.65	0.00/0.00 0.00/0.00 0.00/0.00	0.00/0.00 0.00/0.00 0.00/0.00	0.00/0.00 0.00/0.00 0.00/0.00	0.00/0.00 0.00/0.00 0.00/0.00	0.00/0.00 1.10/1.10 1.10/1.10	0.00/0.00 0.00/0.00 0.00/0.00	0.00/0.00 0.00/0.00 0.00/0.00	0.00/0.00 0.00/0.00 0.00/0.00	0.00/0.00 0.00/0.00 0.00/0.00	0.00/0.02 0.00/0.04 0.00/0.00	0.00/0.00 0.00/0.00 0.00/0.00	0.00/1.02 0.00/0.42 0.00/0.26	0.00/0.03 0.00/0.00 0.00/0.00	0.00/0.29 0.00/0.07 0.00/0.00	0.00/0.68 0.06/0.06 0.06/0.06	0.00/0.00 0.46/0.46 0.46/0.46	0.00/0.13 0.07/0.07 0.07/0.07	0.00/0.01 0.00/21.95 0.00/0.01	0.00/0.02 0.18/0.19 0.18/0.19	0.00/0.00 0.14/0.14 0.14/0.14	0.00/0.14 0.00/0.03 0.00/0.03	0.00/0.10 0.18/0.21 0.18/0.21	0.00/0.15 0.00/0.04 0.00/0.01	0.00/0.25 0.38/30.93 0.38/0.40	0.01/0.09 0.25/0.32 0.25/0.28	0.00/1.25 0.00/0.49 0.00/0.47	0.00/1.17 0.05/1.04 0.05/0.68	0.00/1.60 0.14/0.80 0.14/0.68	0.00/2.15 0.00/1.09 0.00/1.75	0.00/1.03 0.08/1.04 0.08/0.64	0.08/0.33 0.08/0.42 0.00/0.23	0.00/0.37 0.00/0.33 0.00/0.19	0.06/0.63 0.06/0.53 0.06/0.40	0.05/0.30 0.10/0.31 0.12/0.37	0.00/0.12 0.06/0.30 0.06/0.20	0.03/0.12 0.19/0.28 0.16/0.26	0.03/0.27 0.26/0.47 0.17/0.38	0.03/0.19 0.18/0.35 0.04/0.26	0.06/0.19 0.00/0.27 0.00/0.19	0.00/0.18 0.06/0.19 0.01/0.14	0.01/0.28 0.20/1.53 0.19/0.33
	Dev	GQAP BRKGA		0.00/0.60 0.00/2.81	0.00/0.53 2.65/3.49	0.00/0.07 0.40/0.67	0.00/0.06 0.92/2.42	0.00/0.61 1.52/3.01	0.10/0.54 0.93/1.91	0.00/0.18 1.68/2.35	0.00/0.10 0.07/2.05	0.00/0.31 1.79/2.11	0.51/0.76 1.12/1.84	0.41/0.51 1.04/1.52	0.00/0.11 0.81/1.07	0.00/0.12 1.09/1.46	0.00/0.02 0.54/0.80	0.39/0.39 2.79/3.21	1.56/3.40 0.00/0.48	1.26/4.44 1.13/3.22	2.21/3.14 0.40/0.44	2.83/8.19 2.17/3.44	3.29/4.67 3.29/3.54	3.97/6.60 5.50/6.59	4.70/5.90 5.17/5.41	5.36/6.48 5.08/5.58	4.53/6.15 2.76/4.40	4.34/6.40 3.44/4.35	3.58/3.79 3.58/4.01	3.34/4.02 3.31/3.92	3.75/4.11 3.15/3.36	2.94/3.59 3.24/3.41	2.71/3.62 3.16/3.50	23.60/28.71 1.79/2.62	25.79/27.27 4.71/5.58	27.30/29.48 3.89/4.81	22.58/34.91 4.34/5.19	26.87/32.12 4.63/5.49	21.77/27.42 6.22/7.22	24.59/28.63 4.99/6.10	27.31/29.75 5.43/5.86	26.22/31.48 7.27/7.56	24.04/28.34 6.75/7.54	15.43/16.11 6.64/6.86	15.47/15.97 6.31/6.73	14.68/15.73 5.21/6.10	15.64/16.56 6.14/6.60	13.52/14.74 6.48/6.84	8.37/10.15 3.19/3.94
		GevPR-	HMP	0.92/0.92	0.00/0.00	0.00/0.00	0.31/0.31	0.64/0.64	0.00/0.70	1.39/1.53	0.86/1.26	1.01/1.36	0.59/0.60	0.00/0.25	0.04/0.07	0.04/0.20	0.00/0.02	0.01/0.33	0.00/0.44	0.36/1.63	0.00/0.00	0.00/0.00	0.33/0.85	0.00/0.38	1.69/2.79	1.94/2.33	0.24/0.98	0.38/0.98	2.86/2.99	2.56/2.80	3.33/3.71	2.70/2.98	0.87/1.32	0.99/1.73	1.20/2.03	1.73/2.15	0.77/1.49	0.41/1.54	4.70/4.90	3.05/3.69	4.45/5.62	4.22/4.79	1.45/2.41	3.24/3.64	2.82/3.15	2.22/2.83	1.91/2.36	3.18/3.47	1.32/1.74
		Name		100 15 270001	100 15 270002	100 15 270003	100 15 270004	100 15 270005	100 25 270001	100 25 270002	100 25 270003	100 25 270004	100 25 270005	100 50 270001	100 50 270002	100 50 270003	100 50 270004	100 50 270005	200 15 270001	200 15 270002	200 15 270003	200 15 270004	200 15 270005	200 25 270001	200 25 270002	200 25 270003	200 25 270004	200 25 270005	200 50 270001	200 50 270002	200 50 270003	200 50 270004	200 50 270005	400 15 270001	400 15 270002	400 15 270003	400 15 270004	400 15 270005	400 25 270001	400 25 270002	400 25 270003	400 25 270004	400 25 270005	400 50 270001	400 50 270002	400 50 270003	400 50 270004	400 50 270005	Average

Table 18

 18 Comparison between FITS and MA on the large handover minimization instances with n ≥ 100. For direct comparisons, we present the results in the minimization form.

	Name	FITS			MA
	f best	favg	tavg (s)	f best	favg	tavg (s)

Table 19

 19 Comparison between FITS and its two underlying components FLS and InfLS on the RanReal benchmark set.

	Name	FITS	FLS		InfLS
	f best	favg	f best	favg	f best	favg

Table 20

 20 Comparison between FITS and its variation FITS FI (using the first-improvement strategy) on the RanReal benchmark.

	Name	FITS			FITS FI
	f best	favg	tavg	f best	favg	tavg

Table 21

 21 Comparison of two MA versions using uniform and cluster-based crossover respectively on the RanReal benchmark.

	Name	Uniform CX			Cluster-based CX
	f best	favg	tavg (s)	f best	favg	tavg (s)

Table 22

 22 Computational results on small handover minimization instances. For direct comparison, the results are converted into the minimization form.

	Name		FITS
		f best	favg	tavg (s)
	20 5 270001	540	540.00	0.00
	20 5 270002	54	54.00	0.00
	20 5 270003	816	816.00	0.00
	20 5 270004	126	126.00	0.00
	20 5 270005	372	372.00	0.00
	20 10 270001	2148	2148.00	0.00
	20 10 270002	1426	1426.00	0.00
	20 10 270003	2458	2458.00	0.00
	20 10 270004	1570	1570.00	0.00
	30 5 270001	772	772.00	0.00
	30 5 270002	136	136.00	0.00
	30 5 270003	920	920.00	0.01
	30 5 270004	52	52.00	0.00
	30 5 270005	410	410.00	0.01
	30 10 270001	3276	3276.00	0.00
	30 10 270002	1404	1404.00	0.00
	30 10 270003	2214	2214.00	0.00
	30 10 270004	2150	2150.00	0.02
	30 10 270005	2540	2540.00	0.04
	30 15 270001	6178	6178.00	0.01
	30 15 270002	4042	4042.00	0.00
	30 15 270003	4126	4126.00	0.00
	30 15 270004	3920	3920.00	0.01
	40 5 270001	610	610.00	0.07
	40 5 270002	136	136.00	0.05
	40 5 270003	234	234.00	0.12
	40 5 270004	232	232.00	1.30
	40 5 270005	774	774.00	0.00
	40 10 270001	4544	4544.00	0.08
	40 10 270002	2068	2068.00	0.00
	40 10 270003	2090	2090.00	0.01
	40 10 270004	1650	1650.00	0.00
	40 10 270005	4316	4316.00	0.01
	40 15 270001	8646	8646.00	0.24
	40 15 270002	4586	4586.00	0.26
	40 15 270003	5396	5396.00	0.02
	40 15 270004	4800	4800.00	0.00
	40 15 270005	6272	6272.00	0.05

http://www.optsicom.es/ccp/

http://www.mi.sanu.ac.rs/ ~nenad/ccp/

ftp://dimacs.rutgers.edu/pub/dsj/clique/

Acknowledgments

We are grateful to the reviewers for their valuable comments which helped us to improve the paper. We thank the authors of [7,24] for kindly sharing with us their source codes that enabled fair comparisons with the proposed algorithms. This work is partially supported by the National Natural Science Foundation Program of China [Grant No. 71401059, 71771099, 71810107003, 71620107002, 71531009] and the Huazhong University of Science and Technology (5001300001).

Appendix

The purpose of this appendix is to show detailed computational results and comparisons between our two proposed algorithms (FITS and MA) and the state-of-the-art algorithms on the complete CCP benchmark consisting of 183 instances (Tables 111213141516171819202122). For each instance and approach, columns 'f best ', 'f avg ' and 't avg ' show respectively the best objective value, the average objective value and the average computing time in seconds required to reach the final solution (see Section 4.2 for the used experimental protocol). Column 'Dev best /Dev avg ' indicates the percent deviation between the best or the average result and the best solutions obtained within each experiment. Row '#Best' gives the number of cases when each algorithm outperforms the remaining approaches, while row 'Average' shows the average result for a given subset of instances. Finally, row 'p-value' indicates the outcome of the nonparametric Friedman tests on the results obtained with FITS and the reference algorithms. The best results are highlighted in bold.