A Scatter Search Approach for the Minimum Sum-of- Squares Clustering Problem Abstract

A Scatter Search Approach for the Minimum Sum-of-Squares Clustering ProblemJoaquín A. PachecoDepartment of Applied Economics, University of Burgos,Plaza Infanta Elena s/n BURGOS 09001, SPAINTf: +34 947 259021; fx +34 947 258956; email: jpacheco@ubu.esLatest version: September 19, 2003AbstractA metaheuristic procedure based on the Scatter Search approach is proposed for the non-hierarchical clustering problem under the criterion of minimum Sum-of-Squares Clustering. This algorithm incorporates procedures based on different strategies, such as Local Search, GRASP, Tabu Search or Path Relinking. The aim is to obtain quality solutions with short computation times. A series of computational experiments has been performed. The proposed algorithm obtains better results than previously reported methods, especially with small numbers of clusters.Keywords: Clusterization, Metaheuristics, Scatter Search, Local Search, GRASP,Tabu Search, Path Relinking1. IntroductionConsider a set X = {x 1, x 2, ..., x N } of N points in R q and let m be a predetermined positive integer. The Minimum Sum-of-Squares Clustering (MSSC) problem is to find a partition of X into m disjoint subsets (clusters) so that the sum of squared distances from each point to the centroid of its cluster is minimum. Specifically, let P m denote the set of all the partitions of X in m sets, where each partition PA ∈ P m is defined as PA = {C 1,C 2, ..., C m } and where C i denotes each of the clusters that forms PA . Thus, the problem can be expressed as:∑∑=∈∈−m i C x i l PA i l m x x 12 min P ,where the centroid i x is defined asl C x i i x n x i l ∑∈=1, with n i = |C i | .The problem can be written as∑=−N l k l lx x 12min ,where k l is the cluster to which point x l belongs.The design of clusters is a well known exploratory Data Analysis issue called Pattern Recognition. The aim is to find whether a given set of cases X has some structure and,in if so, to display it in the form of a partition. This problem belongs to the area of Non-Hierarchical cluster design, which has many applications in economics, social and natural sciences. It is known to be NP-Hard [5].Various exact methods for MSSC can be found in the literature (see, for example [17] and [7]), some of them, such as the method proposed by du Merle et al. [8], have succeeded in resolving problems with up to 150 points. For larger-sized problems the use of heuristic algorithms is still necessary. The most popular are those based on Local Search methods, such as the well-known K-Means [15] and H-Means [14] procedures. In a recent work, Hansen and Mladenovic [13] propose a new Local Search procedure, J-Means, along with variants H-Means+ or HK-Means. In recent years algorithms using Metaheuristic strategies have been designed, such as Simulated Annealing [16], Tabu Search [2], Genetic Algorithms [3] or most recently Variable Neighborhood Search or VNS [8], [13] and Memetic Algorithms [19].An algorithm that is able to obtain good solutions in short times is proposed for this problem. This method is based in the a recent Metaheuristic strategy named Scatter Search (SS). This method also incorporates others procedures based in others methods, such as Local Search, Tabu Search, GRASP and Path Relinking. This Scatter Search approach is analyzed and compared with other recents techniques. In all cases, our proposed technique gives adequate solutions, compared with others recent techniques, in reasonable time, especially with small values of m.2. Solution ApproachThe solution approach that we have developed for the MSSC problem consists of an adaptation of scatter search. Scatter Search is an instance of the so-called evolutionary methods, which is not based solely on randomization as the main mechanism for searching. SS has been successfully implemented in a variety of settings including combinatorial optimization and nonlinear optimization in continuous variables. Scatter search embodies principles and strategies that are still not emulated by other evolutionary methods, and that prove advantageous for solving a variety of complex optimization problems. More about the origin and multiple applications of scatter search can be found in Glover [11], Glover, Laguna and Martí [12] and Laguna [18].SS uses a set of solutions named Reference Set, (RefSet). RefSet is composed by the best solutions by quality (subset RefSet1) and diversity (subset RefSet2). In each iteration new solutions are generated from those of Reference Set, and then Reference Set is updated whit these new solutions. A ‘static’ version is designed for this problem. An outline of this implementation is shown below.Procedure Static_Scatter_SearchStep 1. Generate an initial set of solutions P by using a Diversification- Generation MethodStep 2. Improve these solutions by an Improvement MethodStep 3. With these solutions build an initial RefSetStep 4. Repeat4.1. Obtain all subsets of pairs from RefSet4.2. Combine these subsets and obtain new solutions4.3. Improve these solutions by the Improvement Method4.4. Update RefSet with these news solutionsuntil RefSet is stable (i.e. no new solutions have been included) {final 4.}Step 5. If max_iter iterations (Steps 1-4) elapse without improvement stop else return Step 1Figure 1. SS implementationLet’s denote by n_pob the size of P (Step 1). Also let’s denote by b1 and b2 the sizes of RefSet1 and RefSet2.For building the initial RefSet (Step 3) first the best solutions (by quality) are taken from P. Then, the next solutions are added to RefSet by diversity. For that, the following measure of diversity is used. Let λ be a solution in P \ RefSet, we defineδmin(λ) = min {dif(λ,λ’) / λ’∈RefSet};where dif(λ,λ’) = number of assignments in λ that are different from λ’. We then select the candidate solution λ that maximizes δmin(λ).The updating of the RefSet is based only on quality. That is, only the new solutions that improve the quality of the worst solution in RefSet are added, (Step 4.4). We now provide descriptions of the diversification, improvement and combination methods.2.1 Diversification MethodOur diversification method is based on GRASP constructions. GRASP, or greedy randomized adaptive search procedure, is a heuristic that constructs solutions with controlled randomization and a greedy function. Most GRASP implementations also include a local search that is used to improve upon the solutions generated with the randomized greedy function. GRASP was originally proposed in the context of a set covering problem [9]. Details of the methodology and a survey of applications can be found in Feo and Resende [10] and Pitsoulis and Resende [20].The method proposed in this paper consists of two stages.In the first, a subset S with points from X and sufficiently far from each other (seed-points) is built.This is denoted as S = {x s1, x s2, …, x sm}.We also use m one-point clusters, corresponding to the seed points, i.e., C i = {x si}, i=1,…,m. In the second stage, the remaining points are assigned to different clusters according to the distance to their centroids.The set S is determined as follows:-Determine x j*, the farthest point from the centroid of X, and do S = {x j*}.-While |S| < m do:1 Calculate ∆j= min{||x j- x l|| : x l∈S}, ∀x j∉S.2. Calculate ∆max = max { ∆j : x j∉S }3. Build L = {x j / ∆j≥α∆max }4. Choose x j*∈L randomly and do S = S∪ {x j*}.The α parameter (0 ≤α≤ 1) controls the level of randomization for the greedy selections. Randomization decreases as the value of α increases. This controlled randomization results in a sampling procedure where the best solution found is typicallybetter than the one found by setting α = 1. A judicious selection of the value of αprovides a balance between diversification and solution quality.The first time that the diversification method is employed (Step 1 in Figure 1), there is no history associated with the number of times each point has been selected as seed point. However, this information is valuable when the method is applied to rebuild the reference. The information is stored in the following array:freq (j ) =number of times that point x j has been selected to belongto S in the execution of Step 0 and previous executions ofStep 8 of Figure 1.The information accumulated in freq (j ) is used to modify the ∆j values in the application of the first phase of diversification method. The modified evaluation is:()maxmax freq j freq j j ∆−∆=∆′βwhere freq max = max { freq (j ) : ∀ i }. The modified j ∆′ values are used to calculate ∆’max and execute the diversification method. Large values of β encourage the selection of seed points that have not been frequently made. The use of frequency information within a diversification method is inspired by Campos, et al. [6].In phase 2 of our diversification method, the points of X \ S are assigned as follows:1. Set A = X \ S , (A = Set of unassigned points)2. For each point x j ∈ A and for each cluster C i , i = 1, …, m , calculate the increase in objective function ∆ij that results from assigning point x j to cluster C i . This value Γij is calculated as follows:21j i i i ij x x n n −⋅+=Γwhere i x is the centroid of C i and n i = |C i |.3. Calculate Γi*j* = min { Γij : x j ∉ A , i = 1, …, m }4. Assign x j * to C i * and set A = A – { x j* }5. If A ≠ ∅ return to 2, else stop.2.3 Improvement MethodAs an improvement method the H-Means+ algorithm [13] is used together with a simple Tabu Search procedure. The H-means+ algorithm is a variant of the well-known H-Means algorithm [14].The Tabu Search procedure has been proposed by Pacheco and Valencia [19] and is described shortly below. This Tabu Search procedure uses the neighboring moves employed in K-Means. These moves consist at each step in the movement of an entity from a cluster to a different one. In order to avoid repetitive cycling when a move which consists in moving point x j from cluster C l to cluster C i is performed, point x j isprevented from returning to the cluster C l for a certain number of iterations.Specifically, defineMatrix_tabu (l , j ) = the number of the iteration in which point x j leaves cluster C l .The Tabu Search method is described below, where P denotes a initial solution with a value f . The parameter τ indicates the number of iterations during which a point is not allowed to return to the leaving cluster. The parameter κ indicates the maximum number of unimproved iterations. In our case we use the stop criterion κ = 100. After different tests, τ was set as m .Tabu Search1 Do Matrix_tabu(i,j) = – τ, i =1,...,m, j = 1,…,N2 Do δ = 0 and P* = P, f* = f and η =0;3 Repeat(3.1) δ = δ + 1(3.2) Determine v i*j* = min {v ij / i =1,...,m; j = 1,...,N ; x j ∉ C i verifyingniter > Matrix_tabu (i,j) + τ orf + v ij < f* (‘aspiration criterion’)}(3.3) Reassign x j* to C i* ;(3.4) Do Matrix_tabu (l*,j*) = δ (l* being the previous cluster of x j*);(3.5) If f (the value of the current solution P) < f* then do: P* = P, f* = f and η = δ ;until (δ – η > κ) or another termination criterionWhere v ij is the change in the value of the objective function when x j is reassigned to C i . The following formula is obtained from Späth [22] to simplify the calculations in K-Means. Let C l be the cluster to which x j belongs, then the value of v ij is calculated as follows:2211j l l l j i i i ij x x n n x x n n v −⋅−−−⋅+=.2.4 Combination MethodNew solutions are generated from combining pairs of reference set solutions. The number of solutions generated from each combination depends on the relative quality of the solutions being combined. Let λi and λj be the reference solutions being combined,where i < j . Assume, as before, that the reference set is ordered in a way that λ1 is the best solution and λb is the worst. Then, the number of solutions generated from each combination is: 3 if i ≤ b 1 and j ≤ b 1, 2 if i ≤ b 1 and j > b 1 and 1 if i > b 1 and j > b 1.Each subset with two elements of RefSet is used to generate new solutions. To do so, we use a strategy called Path Relinking. The basic idea is as follows: a path to join the two initial solutions is built. A number of intermediate points (solutions) from the ‘path' or ‘chain’ are selected as new solutions. The aim is for the intermediate points or solutions to be as equidistant as possible from each other. The improvement method described in Sect. 2.2 is then applied to these intermediate solutions. Figure 2 depicts this idea .λiλjλ*λ**Figure 2.- Generation of New Solutions by using Path RelinkingFrom every pair of solutions of the reference set, in figure λi and λj , a path to join them is built. Solutions in these intermediate preselected positions within the path are selected and improved. In this way new solutions are generated, (in figure λ* y λ**).In order to build the path that joins λi and λj , points with different clusters are selected in the different steps and the shift is carried out. At each step the best possible shift is selected. In this way, the intermediate solutions in each step have another element in common with λj .Path Relinking is a strategy traditionally associated with the intensification phase of the Tabu Search. The underlying idea is that in the path between two good solutions there should be solutions of similar quality (in some cases, even better solutions). See Glover, Laguna and Martí [12] for more details.3. Parameter Fine TuningOne of the most time consuming tasks in the development of metaheuristic procedures for optimization is the tuning of search parameters. Adenso-Diaz and Laguna [1] propose an automated parameter tuning system called CALIBRA that employs statistical analysis techniques and a local search procedure to create a systematic way of fine-tuning algorithms. The goal of CALIBRA is to provide an automated system to fine-tune algorithms, where a user needs only to specify a training set of instances and a measure of performance.CALIBRA works with a set of problem instances and a range of values for each search parameter. The procedure utilizes design of experiments and heuristics to search for the best parameter values. The quality of a set of values is tested on the specified set of problem instances. CALIBRA is available at opalo.etsiig.uniovi.es/~adenso/file_d.html ,where a user’s manual can also be found.Our complete set of test problems consists of 15 instances, as we point out in Section4. We selected a relatively small training set of 3 problems, because we determined that this sample was representative and that the parameter values found with CALIBRA would also perform well when applied to the entire test set. The parameters to be adjusted were α and β in the (0.1, 0.9) range and b 1 and b 2 in the (3, 7) range. The values for MaxIter and PSize were set to 2 and 20, respectively. CALIBRA obtainedthe following values, after running for approximately 5 hours on a Pentium III machine at 600 MHz: α = 0.8, β = 0.5, b1 = 5 and b2 = 5.Although a computational time of 5 hours may seems excessive, it actually represents a very reasonable effort when compared to manual parameter fine-tuning. We use these parameter values for all of the experiments reported in Section 4.4. Computational ExperimentsTo compare the efficiencies of our Scatter Search algorithm and those of other recent strategies, a series of tests is performed. Next the results of this set of computational experiments using these algorithms are shown. The following algorithms are tested: Memetic-HK:Memetic Algorithm proposed in Pacheco and Valencia [19],HybMem:Hybrid Algorithm proposed in Pacheco and Valencia [19],VNS-HK:the VNS Algorithm, proposed in Hansen and Mladenovic [13],with HK-Means as the Local Search methodVNS-J:the VNS Algorithm, with J-Means as the Local Search method.SS:Scatter Search Algorithm proposed in Section 2.For this work we have used our own implementation of every algorithm. The HybMem algorithm uses an initial population generated by the greedy-random method described in Beltrán and Pacheco (2001). The VNS-HK and VNS-J algorithms use the best of these solutions as the initial solution.The sets of problems instances used in testing are:(i)the 575(ii)1060 and(iii)3038 points of the plane taken from TSPLIB [21] data base.All the tests in the current work are performed on a personal computer with a Pentium III 600 MHz processor. All the algorithms have been implemented in Pascal using Borland Delphi 5.0.In (i) the running time is limited to 300 seconds for all the algorithms. We set m = 5, 10, 15, 20, 30, 40, … and 100. In tables 1 the solutions obtained for each algorithm are presented.Table 1. Results for 575 TSPLIBm Memetic HybMem VNS-HK VNS-J SS52497887,42497887,42497887,42497887,42497887,4101110033,791110033,791110033,791110033,791110033,7915728530,353728530,353728530,353728530,353728530,35320531939,613531913,582531913,582532170,717531913,58230348880,815349160,661348823,56351480,12348563,42940256905,514256005,663255201,708255576,465254968,53150198428,047196869,891197830,014198256,332196087,80160158411,762157935,431156756,025157871,603156611,97970132176,153130547,63129276,624131955,541128863,02180111539,837111539,837109963,018111344,912110066,169097267,41896557,81394828,3896896,90694858,46810085107,01784710,19282735,24984502,78182551,242In (ii), we set m = 10, 20, 30, ..., 150. These test data were previously used in Hansen and Mladenovic [13], where the best solution known for every value of m, (except m = 40) is reported. These were obtained on a SUN Ultra I System workstation with 10 minutes computation time. The running time is limited to 600 seconds for all the algorithms. In table 2 the solutions obtained for each algorithm and the percent deviation with respect to the best-known solution (dvt) are presented.Table 2. Results and percent desviation respect to the previous best know solution.m Memetic HybMem VNS-HK VNS-J SS1017548402141754840214175808949717580894971754840214dvt : 000,1850,1850 2079179459679179459679179459679179459679179459600000482123632481369435489705165489705165481251643 300,1810,0251,7571,7570347420992343999852341675604352788332341342886 40-----260361695259291682257371638261266903256426137 501,8991,480,7292,2540,359202854402202211137201649998202070992197376472 602,8292,5032,2192,4320,052163880189161259691159593209160880193158450591 703,4271,7730,7211,5340134082965133270032130401981134420821129315712 804,0293,3981,1734,2910,330114081006114035202110417656113014913110705227 903,3183,27702,3530,26199766477,899581387,896671404,598925896,796860491,7 1003,5163,3240,3042,6440,50088519719,887146237,28502619288392643,785304193,3 1104,332,7110,2134,180,54078641016,678641016,67614941377889227,476554412,7 1204,0234,0230,7273,0281,263705148777051487768576048,770217658,168466017,6 1304,3834,3831,5133,9431,35064483367,164230385,562148721,962668365,561623521,4 1405,55,0861,682,5310,82158823525,558344170,256130419,657963193,956211593,8 1505,1824,3250,3663,6430,511 In (iii), we set m = 10, 20, …, 50, 100, 150, ..., 400. These test data were previously used in Hansen and Mladenovic [13], where the best solution known for every value of m, is reported. These were obtained on a SUN Ultra I System workstation with 50 minutes computation time. The running time is limited to 600 seconds for all the algorithms. In table 9 the solutions obtained for each algorithm and the percent deviation with respect to the best-known solution (dvt) are presented.Table 3. Results and percent desviation respect to the previous best know solution.m Memetic HybMem VNS-HK VNS-J SS1056025119156025119156025119156025119156025119100000m Memetic HybMem VNS-HK VNS-J SS266832945266812482266812482266855872266812482 200,008000,0160175789668176114040176843237176960722175598606 300,1090,2940,7090,7760126321403126106594126542821126689259124961056 400,1700,3460,462-0,90998943049,59886223998876217,79893161498340138,6 50-0,001-0,083-0,068-0,012-0,6148962285,24872596348166635,648943300,847899112,7 1002,6042,1090,9372,5640,37631371538,931371538,930685399,731333000,830709692,4 1502,6652,6650,4192,5390,49922690957,922690957,922005595,922646784,822311420,1 2003,5243,5240,3973,3221,79217288153,917288153,916852031,417234369,517114643,3 2503,7693,7691,1513,4462,72713820368,713820368,713461917,413747375,313657457,4 3003,4923,4920,8082,9462,27211570674,911557842,111251477,711536793,311407670,3 3504,264,1441,3843,9552,7919787059,349787059,349577893,19765280,479710487,94 4003,9883,9881,7663,7573,1758491442,198491442,198355290,558488212,658373889,37 4503,2223,2221,5673,1821,7937408097,67408097,67311187,827377619,977338106,5 5002,3922,3921,0521,971,424Tables 1-3 yields the follow observations:-The HybMem and Memetic-HK algorithms yield good solutions when the number of clusters, m, is low.However, the quality of their solutions deteriorates as the value of m increases.-The VNS-J, and specially the VNS-HK algorithms become more ‘competitive’ than the other algorithms as the value of m increase. For high values of m (m≥ 50 in (i), m≥ 90 in (ii), m≥ 200 in (iii)), the VNS-HK algorithm yields the best solution out of all the strategies proposed in 14 cases.-In 26 out of 41 cases, the SS algorithm yields the best solution.Furthermore, it is especially efficient with lower m values (m < 50 in (i), m < 90 in (ii), m < 200 in (iii)).In such cases, it always yields the best solution (except m = 150 in (iii)).Besides in (ii) SS equals the earlier best known solution for m =10, 20, 30, and 70 (for m=40 no solution was known).For m = 50, 60, and 80 the deviations from the earlier best known solution are very small. In (iii) SS equals the earlier best known solution for m=10, 20, 30 and improves that solution for m = 40 and 50.-For higher values of m, SS still proves to be very competitive, i.e., in six out of twenty cases it yields the best solution (m = 130 and 140).In the other fifteen instances, it holds second place, after the VNS-HK algorithm.Deviations in relation to the earlier best solutions are still low.-Globally speaking, the VNS-HK and the SS algorithms are the best strategies in these tests.SS seems to be better, as it yields the best solution in most cases (26 compared to 22).- Another interesting aspect is the robustness that SS shows in comparison to otherstrategies. In (ii) no case is the best known solution yielded by SS further than1.350%. VNS-HK goes up to2.219%. The other strategies show greater distances:up to 5.5% for Memetic-HK, 5.086% for HybMem, and 4.291% for VNS-J.Figure 3 shows the evolution of the solution found for the different algorithms,according to the computational time for N = 1060 and m= 80.1290000001300000001310000001320000001330000001340000001350000001360000001370000000200400Computational Time in SecondsFigure 3.- Evolution of the solution value found for the different strategies according tothe computational timeAs shown in figure 3, VNS-HK, and especially SS, have a very greedy character.After a few seconds, they yield solutions that are better than those produced by other strategies after using up the maximum computational time (600 seconds).5.- ConclusionsThe contribution of our work is the development of a specialized and sophisticated scatter search procedure for the solution of the cluster design problem. This contribution is important, because our diversification, improvement and combination methods introduce novel features that can be adapted to other situations. Our method brings together several search strategies and mechanisms, such as GRASP constructions, local search, tabu search and Path Relinking within the framework of scatter search.Using instances from the literature, we were able to show the merit of our SS design.In particular, our experiments show that our method obtains the best overall solutions against others recent methods. Our method is specially effective for small values of m .With higher values of m it is only surpassed by the VNS algorithm proposed by Hansen and Mladenovic (2001), with HK-Means as the Local Search Procedure. However, the following is worth noting: in their work, VNS is used with J-Means and J-Means+ as a local search procedure, whereas in our work the use of J-Means in VNS leads to worse solutions than VNS with HK-Means and also to worse solutions than the methods we propose.Finally our scatter search approach shows two important features: it is a ‘greedy’method (it found good solutions in very short calculation time) and is able to ‘evolve’(improve these solutions with more calculation time).References[1] Adenso-Díaz B and. Laguna M. Automated Fine Tuning of Algorithms with Taguchi FractionalExperimental Designs and Local Search. University of Colorado at Boulder; 2001.[2] Al-Sultan KH. A Tabu Search Approach to the Clustering Problem. Pattern Recognition 1995;28:1443-1451.[3] Babu GP and Murty MN. A Near-Optimal Initial Seed Value Selection in K-means Algorithm usingGenetic Algorithms. Pattern Recognition Letters 1993;14:763-769.[4] Beltrán M and Pacheco J. Nuevos métodos para el diseño de cluster no jerárquicos. Una aplicación alos municipios de Castilla y León. Estadística Española. Instituto Nacional de Estadística 2001;43, 148:209-224. (available in www.ine.es)[5] Brucker P. On the Complexity of Clustering Problems. Lecture Notes in Economics andMathematical Systems 1978;157:45-54.[6] Campos V, Glover F, Laguna M and Martí R. An Experimental Evaluation of a Scatter Search for theLinear Ordering Problem. Journal of Global Optimization 2001; 21:397-414.[7] Diehr G. Evaluation of a Branch and Bound Algorithm for Clustering. SIAM put.1985; 6:268-284.[8] du Merle, O., Hansen, P., Jaumard, B. and Mladenovic, N. An Interior Point Algorithm for MinimumSum of Squares Clustering. SIAM Journal on Scientific Computing 2000; 21(4):1485-1505.[9] Feo, T.A. and Resende, M.G.C. A Probabilistic heuristic for a computationally difficult Set CoveringProblem. Operations Research Letters 1989; 8:67-71.[10] Feo, T.A. and Resende, M.G.C. Greedy Randomized Adaptive Search Procedures. Journal of GlobalOptimization 1995;2:1-27.[11] Glover, F. A Template for Scatter Search and Path Relinking. in Artificial Evolution, Lecture Notesin Computer Science, 1363, J.-K. Hao, E. Lutton, E. Ronald , M. Schoenauer and D. Snyers (Eds.) Springer 1998, pp.13-54.[12] Glover, F., Laguna, M. and Martí, R. Fundamentals of Scatter Search and Path Relinking. Controland Cybernetics 2000; 39,3:653-684.[13] Hansen, P. and Mladenovic, N. J-Means: A new Local Search Heuristic for Minimum Sum-of-Squares Clustering. Pattern Recognition 2001; 34(2):405-413.[14] Howard, R. Classifying a Population into Homogeneous Groups. In Lawrence, J.R. (eds.),Operational Research in the Social Sciences. Tavistock Publ., London 1966.[15] Jancey, R.C. Multidimensional Group Analysis. Australian J. Botany 1966;14:127-130.[16] Klein, R.W. and Dubes, R.C. Experiments in Projection and Clustering by Simulated Annealing.Pattern Recognition 1989; 22: 213-220.[17] Koontz, W.L.G., Narendra, P.M. and Fukunuga, K. A Branch and Bound Clustering Algorithm. IEEETransactions on Computers 1975; C-24:908-915[18] Laguna, M. Scatter Search, in Handbook of Applied Optimization, P. M. Pardalos and M. G. C.Resende (Eds.), Oxford University Press, New York 2002, pp. 183-193.[19] Pacheco, J. and Valencia, O. Design of Hibrids for the Minimun Sum-of-Squares Clustering Problem.Computational Statistics and Data Analysis 2003; 43,2:235-248.[20] Pitsoulis, L.S. and Resende, M.G.C. Greedy Randomized Adaptive Search Procedures in Handbookof Applied Optimization, P. M. Pardalos and M. G. C. Resende (Eds.), Oxford University Press 2002, pp. 168-182.[21] Reinelt, G. TSPLIB: A Travelling Salesman Problem Library. ORSA Journal on Computing 1991; 3:376-384.[22] Späth, H. Cluster Analysis Algorithms for Data Reduction and Classification of Objects. EllisHorwood, Chichester 1980.11。