A stepped tabu search method for the clique partitioning problem

Given an undirected graph, a clique is a subset of vertices in which the induced subgraph is complete; that is, all pairs of vertices of this subset are adjacent. Clique problems in graphs are very important due to their numerous applications. One of these problems is the clique partitioning problem (CPP), which consists of dividing the set of vertices of a graph into the smallest number of cliques possible. The CPP is an NP-hard problem with many application fields (timetabling, manufacturing, scheduling, telecommunications, etc.). Despite its great applicability, few recent studies have focused on proposing specific resolution methods for the CPP. This article presents a resolution method that combines multistart strategies with tabu search. The most novel characteristic of our method is that it allows unfeasible solutions to be visited, which facilitates exploration of the solution space. The computational tests show that our method performs better than previous methods proposed for this problem. In fact, our method strictly improves the results of these methods in most of the instances considered while requiring less computation time.


Introduction
Given an undirected graph G = (V, E), where V represents the set of vertices and E represents the set of edges, a clique is a subset of V constituted by vertices that are pairwise adjacentthat is, for every pair of vertices i and j in such a subset, (i, j) ∈ E. A clique partition of graph G is a partition of V such that every subset in the partition is a clique. This work focuses on the problem of finding a clique partition of a graph G with the minimum cardinality. This problem is known as the clique partitioning problem (CPP) or the minimum clique partition problem, and it is an NP-hard problem [19]. Figure 1 of Section 2 shows an example of an undirected graph with 9 vertices, and Fig. 6 presents an optimal solution to this problem with three cliques.
The problem described above should be distinguished from another problem of the same name, which was addressed in Lü et al. [12] and Hu et al. [8] among other references. This problem is defined on complete graphs with weights on the edges. The aim is to divide the set of vertices into subsets such that the sum of the weights of the induced subgraphs is as small as possible.
Clique problems are very popular in the literature on graph problems. Perhaps the most widely known is the maximum clique problem (MCP), which consists of searching in G for a maximum cliquethat is, a clique of maximum cardinality. Another commonly used concept is the maximal clique, which is a clique that is not contained in another clique. A maximum clique must be maximal, but a maximal clique does not have to be maximum. Pardalos and Xue [15] performed a thorough analysis of the maximum clique problem (considering formulations, complexity, algorithms, etc.), and Wu and Hao [21] conducted a very complete survey of this problem. Other important works are those of San Segundo and Artieda [16] that shows a practical application related to unmanned autonomous vehicles and San Segundo et al. [17], a methodological paper proposing improvements in the exact algorithms for this problem. Other known problems regarding cliques are the maximal clique problem [9], the clique enumeration problem [5] and the clique coloring problem [11].
Two of these applications are described below. The work of Allignol et al. [1], addresses the problem of aviation authorities assigning flight levels (FLs) to different flights in a given planning period. Flight levels are the heights at which aircraft fly after take-off and before landing (known as "cruising altitudes"). Flight levels have a margin of 1000 ft (the difference in altitude between two consecutive levels). The problem is to assign different FLs to "incompatible" flights, i.e., flights whose trajectories may be less than 10 nautical miles apart at or near the same point in time. Assigning different FLs to incompatible flights involves establishing safety margins and avoiding accidents. The problem is stated as a graph in which each flight corresponds to a vertex and each edge joins two vertices (flights) that are compatible and can therefore fly without any type of conflict between them at the same FL. The idea is to divide the flights into groups such that all flights in a group are compatible. Therefore, each of these groups corresponds to a clique in the network. Each group of compatible flights (clique) is assigned a different FL. The objective is to minimize the number of groups and therefore the number of FLs to be used, or at least to ensure that this number does not exceed a previously established maximum number of FLs. Thus, the problem can be considered a real application of the CPP. However, other aspects are also considered, such as the costs of flying at each FL for each flight. Therefore, airlines request an ideal FL for each of their flights (normally the FL such that the flight would have the lowest fuel cost), and the allocation of FLs to flights must take these requests into account (by minimizing the number of allocations that are different from the ideal FLs, balancing the number of different allocations between airlines, etc.).
The work of Casado et al. [4] addresses the improvement of steel coil production by grouping processes. In the production of steel coils, raw materials (mainly iron and carbon but also manganese, silicon and niobium, among others) are first extracted. These ores are combined to produce steel in the form of slabs. A rolling mill then converts the slabs into coils. The characteristics of the coil determine the characteristics of the slabs. The problem concerns the production of a series of orders of steel coils requested by different customers. The coils in each order are the same as each other and different from those in other orders. Traditionally, it was ensured that there would be a one-to-one correspondence between the different orders of coils and the types of slabs being produced. However, the possibility of grouping the orders into compatible groups is now considered so that orders in the same group can be manufactured from the same type of slab. As an example, two orders of coils are compatible with respect to the type of steel if the types of steel required are of the same family (steels are divided into families, and within each family, they are ordered from lower quality to higher quality). In this case, the coils could be produced from slabs of the higher-quality steel type between the two initially needed. To check the compatibility between orders, other aspects, such as coil weights and widths, must be considered in addition to the type of steel. Grouping reduces the number of slab types that have to be manufactured: 1 type is needed per group or "cluster" instead of 1 per order. Thus, grouping has advantages such as greater continuity in the production process, cost reduction (in industrial processes with large equipment, the cost of stopping and restarting is very high), less chance of accidents and breakdowns, and improvement in storage and inventory management. This problem can be modeled as a graph in which the nodes are the orders and the edges link compatible nodes (orders). Cliques represent groups of compatible orders (and can therefore be produced by identical slabs). The objective is to minimize the number of groups and thus the number of types of slabs to be produced. It is also necessary to consider the costs that the grouping incurs (for example, by considering the cost of producing coils with steel of the same family but of a higher quality than that initially required). In this way, between two solutions with the same number of cliques, we can choose the one with the lowest cost.
The literature regarding theoretical and methodological aspects related to the CPP is less abundant than that related to other clique problems (especially the MCP). For example, Bhasker and Samad [2] formulated a new upper bound for the number of cliques, and they demonstrated that there is an optimal partitioning that includes a maximal clique (although not necessarily a maximum clique). The more recent publication of Sundar and Singh [19] developed two metaheuristic techniques based on evolutionary computation to solve the CPP. The proposed approaches were tested on 37 publicly available DI-MACS graph instances. More recently, Casado et al. [4] designed a method based on tabu search for an extension of the CPP to the field of manufacturing. This method was adapted to the CPP and tested on the same 37 DI-MACS instances. The main contribution of this work is the design of a new heuristic method for the CPP. This method combines multistart strategies with tabu search. The most interesting characteristic of our method is the type of solution exploration that is used in the tabu search procedure. As described in more detail in the next section, the method allows infeasible solutions to be visited, which will in turn make the search more flexible. In this way, the method will be able to reach feasible solutions with fewer cliques and shorter computation times than previous methods, as shown in Section 3.
The rest of the article is organized as follows: Section 2 shows a mathematical formulation of the CPP, Section 3 describes the resolution method in detail, Section 4 presents different computational tests and Section 5 details the conclusions.

Mathematical formulation
Let G = (V, E) be a graph with n vertices (i.e., V = {1, 2, …, n}). Let k max be an upper bound of the optimum number of cliques, and let C k , k = 1…k max denote each of the cliques. The problem can be formulated as a mathematical program with the following three sets of variables: x ik A binary variable that equals one if product i belongs to clique C k. y k A binary variable that equals one if clique C k is not empty n k The number of vertices assigned to clique C k The problem can be formulated as follows: Subject to In this formulation, each variable y k indicates whether clique k has an element or is empty; each variable x ik indicates whether vertex i is assigned to cluster k; and each variable n k indicates the number of elements of clique k. Target function (1) indicates the number of nonempty cliques. Restrictions (2) force y k = 1 if any vertex i is assigned to clique k, restrictions (3) force the assignment of each vertex to a clique, restrictions (4) are cardinality restrictions on each clique, and restrictions (5) force each clique to be constituted by adjacent vertices. Last, an obvious value for k max would be to take k max = n. However, to avoid excessively long formulations of the instances, the value obtained by the Constructive procedure (described in Pseudocode 2 with α = 1) was used as the value of k max .

Resolution method
Let S represent a generic solution of the corresponding CPP, K be the number of nonempty cliques that constitute it, and C 1 , C 2 …C K represent these cliques (in some cases, as an abbreviation, S = {C 1 , C 2 …C K }). Therefore, K is the target function to minimize.
The proposed method is a multistart procedure. Thus, it is an iterative process in which a solution is created in each iteration, which is then improved by a tabu search (TS) procedure. The novelty of this work is this TS procedure, which allows unfeasible solutions to be visited, thereby making the search more flexible. As described in Section 4, this strategy will allow the results obtained by previous procedures for this problem to be improved considerably. This TS procedure is described in detail below, as well as the entire proposed multistart method.
Specifically, subsection 3.1 explains the basic idea of the tabu search procedure and, more specifically, the moves it uses. Subsection 3.2 shows the flowchart of the tabu search and how moves are chosen. In subsection 3.3, some aspects of the tabu search procedure are explained in more detail and described in pseudocode. In subsection 3.4, the multistart procedure into which the tabu search procedure is integrated is described. Finally, in subsection 3.5, the differences from other heuristic methods for the CPP are explained.

Basic idea of the tabu search method: Movements
We begin with an initial feasible solution S, constituted by K nonempty cliques C 1 , C 2 …C K . K groups Gr 1 , Gr 2 …Gr K are formed, which are initially identified with each clique; i.e., Gr k = C k , k = 1, . . , K. A group is "removed" or "deactivated", and its components are reassigned to other groups. This can lead to an unfeasible solution, since some of the groups may not be cliques; that is, there may be groups with unlinked pairs of vertices. Hereafter, for simplicity, each pair of unlinked vertices in the same group will be called an "incompatibility".
In the following steps, modifications or movements are made in the solution with the aim of reducing the number of incompatibilities. If the number of incompatibilities is reduced to zero, a feasible solution will be obtained (all groups are cliques) with one fewer clique (K − 1). The movement considered in each step is to change a vertex of one group to another group (considering only active groups).
In the next step (Fig. 3), group Gr 1 is removed ("deactivated"), and its elements are reassigned among the other groups. Specifically, vertex 1 is reassigned to group Gr 2 . As shown in Fig. 4, there is an incompatibility in this group corresponding to vertices 1 and 3 (marked in red). Then, the movement that most reduces (or least increases) the number of incompatibilities is searched for and executed. Specifically, vertex 3 is moved to group Gr 3 . This eliminates the incompatibility between vertices 1 and 3, but a new incompatibility appears between vertices 3 and 6 in group Gr 3 (Fig. 5). Last, vertex 6 is moved to group Gr 4 , which removes all incompatibilities (Fig. 6 With this new feasible solution, this sequence or block of steps is repeated (removing a group, reassigning its elements, and performing movements to eliminate incompatibilities). Therefore, this is an iterative process in which the number of cliques is reduced by one unit in each block of steps. The entire process ends when it is not possible to remove the incompatibilities in one block.

Flowchart of the tabu search method
As mentioned above, this method consists of executing a sequence of blocks, where in each block, the number of clicks is reduced by one unit. The process ends when it is not possible to reduce the number of cliques in a block. The process for each block is to choose a group and "empty" it, i.e., assign the elements to other groups. This can lead to incompatibilities. In the next steps, moves are made to reduce these incompatibilities until they are eliminated, if possible. These moves consist of changing a vertex from one group to another. A block ends when a stop criterion is reached. This criterion is reached when all incompatibilities have been eliminated or when a certain number of iterations have elapsed without reducing the incompatibilities. If all incompatibilities have been eliminated, a feasible solution is found with one fewer clique (all groups are cliques), and a new block is started; otherwise, the process ends with the output being the best feasible solution found. Figure 7 shows the flowchart of this process.
In the flowchart of Fig. 7, C k , k = 1, . . , K, are the cliques that define the initial solution; C k * , k = 1, . . , K, are the cliques  in which the best solution is found and therefore in which the output of the method is stored. Note that some of these cliques can be empty. In this case, it is necessary to eliminate them and renumber the rest. NInc is the variable in which the number of incompatibilities is stored. The set of criteria or method by which to select a move in each step follows a basic tabu search strategy [7]. Initially, the best move is chosen, i.e., the one that reduces NInc the most or increases it the least; additionally, to prevent the algorithm from cycling, some moves are declared "tabu" and are not considered. In this case, returning a vertex to a group from which it was moved in recent iterations is declared tabu. The tabu status can be ignored if it results in a lower NInc value than that found in previous iterations of the block. The next subsection will explain the details of the tabu search method (which we name the SteppedTabu procedure) and describe its full pseudocode.

The SteppedTabu procedure: Details and pseudocode
Next, some aspects of the previous process are explained in more detail. The first aspect to explain is which criterion is used to select the group that is initially removed and the groups to which its vertices are reassigned. The idea is to select the group and perform the corresponding reassignments so that the number of "incompatibilities" increases as little as possible.
Let S be a solution (feasible or not) constituted by K groups, Gr 1 , Gr 2 …Gr K ; we define: Similarly, we define: Therefore, km(i) indicates, for each vertex i ∈ V, the index of the "active" group (without considering the group it belongs to) with fewer vertices that are "incompatible" with i (i.e., vertices j such that (i, j) ∉ E). Thus, it is the index of the group that would produce the smallest number of incompatibilities if i were "reassigned" to that group. SInc(k) indicates the number of incompatibilities if group Gr k is removed and each of its members i is reassigned to group Gr kmin(i) . Last, k1 is the index of the group that would produce the smallest number of incompatibilities if it were removed and each of its elements i were reassigned to the corresponding groups Gr km(i) . Consequently, group Gr k1 is eliminated.
As an example, the following table shows the initial values of Inc(i, j), which correspond to the initial solution of Fig. 2. The lines correspond to the vertices, and the columns correspond to the groups. For each vertex i, the group that corresponds to km(i) is highlighted in gray (Table 1). Table 2 shows the values of SInc. Therefore, the "best" group to be removed (j1) could be group. Gr 1 , Gr 2 or Gr 3 . The draws are resolved in lexicographic order, although other criteria can be considered. Thus, in this case, we set k1 = Gr 1 ; we "deactivate" this group and reassign its only element, vertex 1, to group Gr 2 . In the subsequent steps, group Gr 1 is no longer considered in the possible movements.
Next, for the following steps, we search for the "best" movement of a vertex i to a different groupthat is, the movement that most reduces the number of incompatibilities in the solution. Specifically, ∀i ∈ V, ∀ k = 1. . K i ∉ Gr k . This increase dif (if it is positive) or reduction (if it is negative) is calculated as: (This is the increase in incompatibilities in group k minus the decrease in incompatibilities in group Grp(i) if the movement were to be executed.) Consequently, we continue to make use of auxiliary variables Inc(i, j), which are conveniently updated. Table 3 shows the values of Inc and dif for the solution of Fig. 4.
The lowest value of dif is dif = 0, and the corresponding movement is to move vertex 3 to group Gr 3 (in this case, we use the lexicographic order of the draws). We execute this movement, obtaining the solution shown in Fig. 5. For the next step, Table 4 presents the values of Inc and dif for the solution of Fig. 5.  The lowest value of dif (−1) corresponds to moving vertex 6 to group Gr 4 . The execution of this movement removes all incompatibilities, resulting in a feasible solution (Fig. 6).
As observed in the previous example, in some cases, it is convenient to admit movements that do not improve the solution (i.e., movements that do not reduce the number of incompatibilities) to prevent the blockage and termination of the process. Thus, as shown in Table 3, in the solution presented in Fig. 4, it is not possible to find any movement that reduces the number of incompatibilities. However, we do not end the process at this point, and we perform the best possible move (moving vertex 3 to groupGr 3 ), although it does not improve the solution. Then, vertex 6 of group Gr 3 is moved to Gr 4 , obtaining a solution without incompatibilities. On the other hand, if we admit movements that do not improve the solution (and even movements that may worsen it), the process must be equipped with mechanisms that prevent it from cycling. Specifically, if a vertex has just left a group, the aim is to prevent it from returning to the group in subsequent iterations. To this end, we define tabumatrix(i, k) ∀ i ∈ V, ∀ k = 1. . K, as follows: tabumatrix(i, k) Number of iterations or steps (within each block) in which vertex i exits group Gr k Therefore, the movement defined by vertex i and group Gr k is declared "tabu" (and its execution is prevented) if: where iter is the number of the current iteration within each block. The parameter tenure indicates the number of iterations in which the return of i to group Gr k is declared tabu after the moment at which it leaves it. Excessively high values of tenure may lead the algorithm to be greatly hindered, and very low values can make it cycle. Consequently, the selection of this parameter is critical for the satisfactory functioning of the process. On the other hand, the tabu state of a movement can be ignored if this movement produces a solution with fewer incompatibilities than any solution visited during that block ("aspiration criterion").
In summary, the entire process is a sequence of blocks. Each block searches for a solution with one fewer clique than the previous block. Therefore, this is a descending "stepped" process (each block is a "step" with one fewer clique than the previous block). Each block follows a tabu search strategy to obtain a feasible solution (without incompatibilities). We call this process or algorithm TabuStepped, and it is thoroughly described in Pseudocode 1. In this pseudocode, the main variables are: -S: a feasible initial solution, where K is the number of cliques that constitute it. These cliques are represented as C 1 , C 2 …C K . -S * : the final solution, that is, the best feasible solution found, where K * is the number of cliques that constitute it (represented as C * 1 ; C * 2 …C * K * ). -Gr k , k = 1, . . , K: the groups that make up the solutions (feasible or not) visited during the process. -Kf: number of active groups at each moment.
In Pseudocode 1, steps 1-4 represent the beginning of the entire process (initiate S * , Gr, Inc, Kf). The repeat-until loop (steps 5-17) represents the execution of each block. Thus, steps 5-9 determine which group k1 is to be deactivated, reassign its elements to other groups and determine the number of incompatibilities (NInc). Steps 10-12 initiate the variables that will be used in the NInc reduction phase (iter, iterbest, tabumatrix, Niterbest). The do-while loop (steps [13][14][15][16] represents each of the iterations of this phase. Thus, step 15 explores the entire set of movements, and step 16 executes the best movement and updates the different variables. Once this do-while loop ends, step 17 verifies whether the number of incompatibilities was reduced to 0. If this is the case, the obtained solution is saved in S * , and a new block is executed. Otherwise, the algorithm ends. The auxiliary variable NIncbest indicates the minimum number of incompatibilities during the NInc reduction phase; NIncbest is used in the "aspiration criterion" (step 15b, which verifies whether (NInc + dif < NIncbest)). The variable iterbest indicates the iteration in which NIncbest was modified (step 16b). maxiter is the parameter of the algorithm that indicates the maximum number of iterations of the do-while loop after the modification of NIncbest. Therefore, maxiter and NIncbest determine the stopping criterion for this loop.
In the search for the best movement (step 15), a movement is considered if it improves Nincbest (NInc + dif < NIncbest) or if it is not tabu (iter > tabumatrix(i, k) + tenure). The values of the best movement found are saved in the variables difb (variation in the number of incompatibilities), ib (vertex) and kb (group).

MultiStartStepped (procedure)
This SteppedTabu procedure is inserted into the multistart procedure, as previously stated. The initial solutions that are subsequently improved by the SteppedTabu procedure are built with the constructive method proposed in Casado et al. [4], which is briefly described in Pseudocode 2.
The "guide" function g(i) measures the suitability of each candidate vertex to be chosen. The parameter α takes values between 0 and 1 and regulates the degree of randomness of the method. If α = 1, the process always produces the same solution (except in cases in which there is an iteration with more than one vertex i ∈ U ′ corresponding to gmax). If α = 0, the values of g(i) are irrelevant, and the process is totally random. A more detailed explanation of this method is provided in Casado et al. [4]. The multistart procedure (which we name MultiStartStepped or MSS) is described in Pseudocode 3.

Pseudocode 3 MultiStartStepped (procedure)
As can be observed, the final solution obtained is S best , with the corresponding number of cliques K best . Therefore, this procedure depends on the following parameters: α (from the constructive procedure), tenure, maxiter (from the SteppedTabu procedure) and maxiterMSS. The next section analyses the different computational tests used to measure the performance of this procedure.

Differences with other recent heuristics
To the best of our knowledge, there are 3 recent heuristic methods for the CPP: two methods proposed by Sundar and Singh [19] and one method proposed in Casado et al. [4]. The two methods proposed by Sundar and Singh [19] are evolutionary methods that work on a population of solutions that can interact with each other. These two evolutionary methods are based on the genetic algorithm and an artificial bee colony. In the case of the genetic algorithm, the Selection, Crossover, Repair, Mutation and Replacement operators/procedures are described in detail. In the case of the artificial bee colony method, the procedure for building the neighboring solutions, selecting solutions for an onlooker bee, etc., are described. In both methods, the same procedure is used to build the initial population, and only feasible solutions are considered.
Our MSS method, as explained above, is a multistart method in which at each iteration, a solution is constructed that is improved by a tabu search procedure. At each step, it operates on a single solution that continually changes. Therefore, there is no interaction between different solutions, and the abovementioned operators are not used. It also allows infeasible solutions to be visited.
The method proposed in Casado et al. [4] has a similar strategy to ours: it is a multistart method that uses a tabu search procedure in the improvement phase. The difference between this tabu search procedure and the one proposed in our work is that the procedure of Casado et al. [4] does not allow moves to infeasible solutions. Neighboring solutions must be feasible; i.e., the sets must be cliques (no "incompatibilities" are allowed). To determine which solutions are better, two hierarchical criteria are used: the number of cliques and the "imbalance" in the number of elements in the cliques. As a measure of "imbalance", the sum of the squares of the number of elements in each clique is taken. Solutions with fewer cliques are preferred, and in the case of equality in the number of cliques, the solution with the highest "imbalance" is chosen.
Thus, among the 3 solutions S 1 ¼ 1; 2; 3; 4; 5; 6 f g ; 7; 8; 9; 10 f g f g S 2 ¼ 1; 2; 3; 4; 5; 6 f g ; 7; 8 f g; 9; 10 f g f g S 3 ¼ 1; 2; 3; 4 f g ; 5; 6; 7 f g; 8; 9; 10 f g f g solution S 1 is preferred over S 2 and S 3 since S 1 consists of two cliques while S 2 and S 3 consist of 3. On the other hand, S 2 is preferred over S 3 since the "imbalance" of S 2 is 44 (36 + 4 + 4) and that of S 2 is 34 (16 + 9 + 9). The idea of using the imbalance criterion is to force empty clusters with fewer elements. Our tabu search procedure is more aggressive: starting from a feasible solution, a clique is eliminated, and its elements are distributed among the other cliques, which produces "incompatibilities", i.e., infeasible solutions. In the next steps, the procedure searches for solutions with fewer incompatibilities until it reaches a solution without incompatibilities, that is, a feasible solution with one fewer clique than the previous feasible solution. This process is repeated with the new feasible solution. Allowing infeasible solutions to be visited gives more flexibility to the search, and as will be seen below, the resulting MSS method produces better solutions than the 3 methods above.

Computational tests
This section describes a series of computational tests used to compare the performance of our multistart algorithm. The algorithm was implemented in Delphi (Object Pascal) with the development environments Rad Studio 10.3 and 11. The tests with this code were conducted using a workstation with an AMD 3990X 2.9 GHz processor with 256 Gb RAM. We also used CPLEX 20.1 in our tests. This section is divided into the following parts: Subsection 4.1 describes the parameter fitting carried out in our method, Subsection 4.2 compares our method with the commercial software CPLEX, and Subsection 4.3 compares the results with other recent heuristic results for this problem.

Parameter fine-tuning
Parameter fitting was carried out using a set of 9 instances of the DIMACS library. Table 5 shows this set and its characteristics. Specifically, it indicates the names, sizes (n), and densities (percentages of links in the network over the total possible number of links) of the instances. The choice of instances has been made in such a way as to combine instances of different size and density. Note that there are 3 instances of small size (C125.9, brock200_2 and p_hat300-1), 3 instances of medium size (DSJC500_5, p_hat700-1 and C1000.9) and 3 instances of large size (p_hat1500-1, C2000.5 and keller6). On the other hand, there are 3 low density instances (p_hat300-1, p_hat700-1 and p_hat1500-1), 3 medium density instances (C2000.5, brock200_2 and DSJC500_5) and 3 high density instances (C125.9, C1000.9 and keller6).
As previously stated, our MSS method has 4 parameters: α, tenure, maxiter and maxiterMSS. The parameter α indicates the degree of randomness of the constructive procedure, and the parameter tenure regulates the number of "tabu" movements in the SteppedTabu procedure. The parameters maxiter and maxiterMSS are used as stopping criteria in the tabu procedure and in the general MSS procedure. To fit the other two parameters (α and tenure), we fixed the values maxiter = 10 · n and maxiterMSS = 20. For α, the values considered were α = 0, 0.1, 0.5, 0.9, 0.99 and 1. The values of tenure considered were tenure ¼ n = 2 ; n; 2 Á n and 5 · n. Of all the combinations, the one that yielded the best results was α = 0.99 and tenure ¼ n = 2 . With these values, we analyzed the parameters maxiter and maxiterMSS. It was determined that there were no significant improvements with values above maxiter = 10 · n and maxiterMSS = 100. Therefore, these values were used in the rest of the tests. Finally, it should be noted that the stopping criterion used by the MultiStartStepped procedure is to achieve maxiterMSS consecutive starts without improvement in solution quality subject to a maximum of 3600 seconds.

Tests with CPLEX
Since the CPP is an NP-hard problem, it is reasonable to expect an excessive computation time to be needed to solve large instances exactly. However, it would be interesting to determine the evolution of these computation times as a function of the size of the instances and determine the maximum size of the instances that can be solved exactly within a reasonable time. If this size is small, it would justify the use of heuristic strategies, such as the one proposed in this work. It would also be interesting to know whether our strategy can find the optimum solution in instances in which such an optimal solution is known and, if so, determine the deviation from this optimal solution. In this section, the commercial software CPLEX is used to solve small instances exactly. The results obtained by CPLEX are compared with those obtained by our MSS. Specifically, instances of sizes n = 10, 20, 30, …, 120 were randomly generated. For each value of n, three values of different densities were considered: low (approximately 0.3), medium (approximately 0.5) and high (approximately 0.7). To avoid excessive computation times, the time used per instance and method was limited to 3600 seconds. Table 6 shows the obtained results. For each method, the value of the best solution found (val) is indicated, and the computation time is shown in seconds (C.T.). In the case of CPLEX, the value of the lower bound obtained is added. If this lower bound coincides with val, then the solution obtained by CPLEX is optimal, and the process ends (if the process did not end within 3600 seconds, this is indicated by "> 3600" in the C.T. column). For each instance, the value of the best solution obtained is indicated in bold. In Table 6, the following can be observed: -CPLEX was able to finish, thereby obtaining the optimal solution, for instances up to size n = 30. In the instances of sizes n = 40, 50, 60 and 70, CPLEX was able to finish in some cases (those of greater density) and not in others. After n = 80, CPLEX was not able to finish in any instance. -It seems that CPLEX behaves better for instances of greater density than those of lower density: among the small instances (n ≤ 30), the computation time was shorter for instances of greater density; among the medium-sized instances (40 ≤ n ≤ 70), CPLEX was able to finish for instances of greater density; and among the large instances (n ≥ 80), the gap between the lower bound and the value obtained was smaller for the instances of greater density. -MSS was able to obtain the best result in all instances: in 20 instances, MSS and CPLEX obtained solutions with the same value, and in 16 instances (in most of the larger instances), MSS obtained strictly better solutions. Moreover, the computation time was very short (it exceeded one second in only one instance). Therefore, in the instances in which the optimal solution was known (those where CPLEX finished), MSS reached the optimal solution with very short calculation times.
In summary, only small instances were resolved exactly (n ≤ 30, and in some instances of 40 ≤ n ≤ 70). Larger instances seem to require the use of heuristic techniques, such as the MSS method proposed. In all the instances in which the optimal solution is known, our MSS method reached the optimal solution.

Computational tests against recent heuristics
As shown in the previous subsection, this problem can only be solved exactly in small instances. This justifies the development of heuristic methods such as our MSS method. In this subsection, we compare our method with other recent heuristics in the literature for this problem. Specifically, we consider the two methods proposed in Sundar and Singh [19]: the method based on genetic algorithms (SSGGA) and the method based on the artificial bee colony (GABC). To check the performance of these methods, the authors used a set of 37 instances from the well-known DIMACS library (http:// dimacs.rutgers.edu/programs/challenge/). These instances were selected due to their difficulty (the presence of overlapping cliques). Then, Casado et al. [4] designed a method based on tabu search (MSTS) for an extension of the CPP to a manufacturing problem. This method was adapted to the CPP and tested in these 37 DIMACS instances. To compare our MSS method with these heuristics, we executed it in these same instances.
In the first set of tests, MSS was run 10 times (as were SSGGA and GABC), and the results obtained are shown in Table 7 together with those of SSGGA and GABC reported in Sundar and Singh [19]. Specifically, Table 7 shows the name, size (n) and density ("den.") of each instance; for each method and instance, the best result of the 10 runs ("Best"), the average of the results of the 10 runs ("Avg.") and the average execution time over the 10 runs ("Time") are shown. For each instance, the best solution is marked in bold.
Regarding the best values obtained ("Best" columns) by each method, our MSS method finds the best solution in 32 out of 37 instances, the SSGGA method finds the best solution in 11 instances, and the GABC method finds the best solution in 6 instances. On the other hand, the MSS method is strictly better than SSGGA in 25 out of 37 instances and strictly worse in 5. The MSS method is strictly better than GABC in 29 out of 37 instances and strictly worse in 1. Regarding the mean values ("Avg." columns), our MSS finds the best mean value in 32 out of 37 instances, the SSGGA method finds the best mean value in 9 instances, and the GABC method finds the best mean value in 5 instances. On the other hand, the MSS method obtains strictly better mean values than SSGGA in 28 out of 37 instances and strictly worse in 5. The MSS method obtains strictly better mean values than GABC in 31 out of 37 instances and strictly worse in 1.
To reinforce the conclusions that can be drawn from Table 7, the following Wilcoxon signed rank tests with the corresponding null and alternative hypotheses are proposed: These two tests are performed both for the best values ("Best" columns in Table 7) and for the average values ("Avg." columns). Therefore, 4 tests are performed. The results of these tests are shown in Table 8. Each table of this test shows the number of instances in which there is no tie (n * ); the sum of ranks with a positive difference (i.e., where the values of the solutions of SSGGA and GABC are higher than those of MSST), denoted by W + ; the sum of ranks with a negative difference (W − ); the lowest values of W + and W − (minW) and the corresponding probability tail (p-tail); and the z value obtained (Z score) and the corresponding probability tail (p-tail z).
As seen in Table 8, considering both the minW values and the Z score values, the corresponding probability tails are insignificant: the null hypothesis is rejected in all 4 tests, and therefore, it can be concluded that MSST obtains significantly better values than SSGGA and GABC (both in terms of the best solution found and in terms of mean values).
In a second set of tests, our MSS is run once, and its results are shown in Table 9 together with those obtained from MSTS (reported in Casado et al. [4]). Table 9 shows the data for each instance (as in Table 7), and for each method and instance, the value of the solution obtained ("Val") and the elapsed computation time to obtain that solution ("Time-Best") are shown. The best solution of each instance is marked in bold.
As seen in Table 9, in the comparison with MSTS, our MSS obtains the best solution in all instances (37), while MSTS obtains the best solution in 7. Therefore, MSS strictly outperforms MSTS in 30 instances, while in the other 7 instances, both methods obtain solutions of equal quality. Applying the Wilcoxon signed rank test, there are 30 instances with strictly positive differences and no negative differences. Therefore, following the notation of Table 8, n * = 30, W + = 465, W − = 0 and minW=0. Therefore, our MSS is significantly better than MSTS.
Therefore, our MSS method considerably improves the results of the previous methods. Regarding the computation times, in some instances where MSS exceeds the other methods, the computation time of MSS is much longer than those of the other methods. This happens in instances keller4, keller5, p_hat300-2, p_hat300-3, p_hat700-1, p_hat700-2, p_hat700-3, p_hat1500-1, p_hat1500-2 and p_hat1500-3. However, observing the evolution of the results with respect to the computation times in these instances, it can be concluded that our method achieves better results than the previous methods in a similar or shorter computation time.  Tables 10 and 11 show the evolution of the value of the best solution found by MSS (column "Val.") and the computation time used to obtain it (column "Time") for each of these instances. These values were obtained in the execution corresponding to Table 9. Specifically, Table 10 shows the evolution corresponding to the instances keller4, keller5, p_hat300-2, p_hat300-3 and p_hat700-1. Table 11 shows the evolution corresponding to the instances p_hat700-2 p_hat700-3, p_hat1500-1, p_hat1500-2 and p_hat1500-3.
As an example, Figs. 8 and 9 show the graphs corresponding to the evolution of the keller5 and p_hat700-2 instances.
In Tables 10 and 11 Table 11 Evolution of the best solution obtained by MSS (in the execution corresponding to Table 9) in the instances p_hat700-2, p_hat700-3, p_ hat1500-1, p_hat1500-2 and p_hat1500-3 p_hat700-2 p _ h a t 7 0 0 -3 p_hat1500-1 p_hat1500-2 p_hat1500-3  It should be noted, however, that the SSGGA and GABC times refer to the average times of the 10 complete runs corresponding to Table 7, and the solution value refers to the best solution obtained in these 10 runs. For MSTS, the solution value refers to the best solution obtained in the run corresponding to Table 9, and the time refers to the time needed to reach this solution.
In summary, from the above, it can be concluded that a) in most of the cases, our method obtains similar or strictly better solutions than the previous methods; b) in most of the cases, our method can reach solutions of the same quality as the final solutions obtained by the previous methods while using less total time than those methods (in some cases, two orders of magnitude less); and c) our method shows that it has the capacity to evolve, as it not only obtains good solutions quickly but is also able to improve them during their execution.
All the solutions obtained by MSS corresponding to Tables 7 and 9 can be found at the following link: www. ubu.es/metaheuristicos-grinubumet/ejemplos-y-datos-deproblemas. However, for a more rigorous comparison of the computation times, we provide Table 12 showing the CPUs used to run SSGGA, GABC, MSTS and our MSS, as well as their relative speeds. Also, we include some benchmarks that can be found at cpu.userbenchmark.com, (specifically those corresponding to the Gaming, Desktop and Workstation tests, as well as the average result).

Conclusions
Clique problems in graphs are interesting because they are attractive from a theoretical perspective and they have a wide range of applications. One of these problems is the so-called Clique Partitioning Problem, which has great applicability in different areas, such as timetabling, manufacturing, scheduling, and telecommunications. This is an NP-hard problem, and as described in this work, it can only be solved exactly in relatively small instances. Despite its practical importance, few resolution methods have been proposed for this problem in the recent literature. This work proposes a heuristic method that uses a tabu search procedure within a multistart strategy. An interesting characteristic of our method is that it allows unfeasible solutions to be visited. This gives flexibility to the exploration of the solution space, and in this way, the method achieves a dramatic improvement in both quality and computational time over the results of previous methods.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This work was partially supported by FEDER funds and the Spanish State Research Agency (Projects PID2019-104263RB-C44 and PDC2021-121021-C22); the Regional Government of "Castilla y León", Spain (Project BU071G19); the Regional Government of "Castilla y León"; and FEDER funds (Project BU056P20).

Declarations
Ethics approval This article does not contain any studies with human participants or animals performed by any of the authors.