Combining ant colony optimization algorithm and dynamic programming technique for solving the covering salesman problem

Computers & Industrial Engineering 83 (2015) 244–251

Contents lists available at ScienceDirect

Computers & Industrial Engineering journal homepage: www.elsevier.com/locate/caie

Combining ant colony optimization algorithm and dynamic programming technique for solving the covering salesman problem q Majid Salari a,⇑, Mohammad Reihaneh b, Mohammad S. Sabbagh c a

Department of Industrial Engineering, Ferdowsi University of Mashhad, 91775-1111, Mashhad, Iran Isenberg School of Management, 121 Presidents Drive, University of Massachusetts Amherst, MA 01003, United States c Department of Industrial and Systems Engineering, Isfahan University of Technology, 84156-83111, Isfahan, Iran b

a r t i c l e

i n f o

Article history: Received 5 February 2014 Received in revised form 8 November 2014 Accepted 23 February 2015 Available online 5 March 2015 Keywords: Traveling salesman problem Covering salesman problem Ant colony optimization Dynamic programming Heuristics

a b s t r a c t The covering salesman problem (CSP) is an extension of the well-known traveling salesman problem in which we are allowed to leave some vertices unvisited. The goal of the CSP is to construct a minimum length Hamiltonian cycle over a subset of vertices where those vertices not visited by the tour need to be within a pre-determined distance from at least one visited vertex. In this paper, we propose a mathematical formulation and a hybrid heuristic algorithm by combining ant colony optimization algorithm and dynamic programming technique to obtain high quality solutions. Comparing the results of the proposed algorithm with available methods in the literature clearly indicates the effectiveness of our proposed heuristic algorithm. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction The covering salesman problem (CSP) is a generalization of the traveling salesman problem (TSP) in which some vertices could remain unvisited (Current & Schilling, 1989). In particular, each vertex, say i, can cover a subset of vertices located within its prespecified distance ri . The goal of solving the problem is to construct a minimum length Hamiltonian cycle over a subset of vertices in which unvisited vertices need to be within the coverage radius of at least one visited vertex by the tour. The authors referred to potential applications of this problem arising in rural healthcare delivery problem. Essentially, the goal is to visit a given number of demand zones such that the people located in unvisited zones have to travel to their nearest visited zone to receive humanitarian functions (Current & Schilling, 1989). Current and Schilling (1989) developed a simple heuristic algorithm for this problem. Their algorithm is composed of two steps. In the first step, a set covering problem (SCP) is solved to find the minimum number of vertices that cover all the vertices of the problem. In the second step, the optimal TSP tour is calculated

q

This manuscript was processed by Area Editor I-Lin Wang.

⇑ Corresponding author at: Department of Industrial Engineering, Ferdowsi University of Mashhad, P.O. Box 91775-1111, Mashhad, Iran. Tel.: +98 0511 8805113. E-mail addresses: [email protected] (M. Salari), [email protected] (M. Reihaneh), [email protected] (M.S. Sabbagh). http://dx.doi.org/10.1016/j.cie.2015.02.019 0360-8352/Ó 2015 Elsevier Ltd. All rights reserved.

over the set of vertices obtained from first step. In case of having multiple optimal solutions for the SCP, all of them are used in the second step and the least costly TSP tour is selected. Golden, Naji-Azimi, Raghavan, Salari, and Toth (2012) proposed two local search algorithms for the CSP. Their proposed algorithms take advantage of several moves like swap, removal–reinsertion, and perturbation to escape local optimal solutions. The first algorithm, denoted as LS1, improves the quality of the solutions by performing a set of removal and re-assignment moves. Essentially, during the first step of the LS1 algorithm, some visited vertices are removed from the tour and in the second step, a new feasible solution is found by substituting the removed vertices with a better combination of unvisited vertices (if any). The second algorithm, denoted as LS2, applies more sophisticated moves to improve the solutions’ quality. Particularly, removal–reassignment and perturbation procedures are applied to achieve this goal. Moreover, the Lin–Kernighan heuristic (Lin & Kernighan, 1973) is also used during this algorithm to find a better rearrangement of the visited vertices. An integer linear programming (ILP) based heuristic algorithm was developed for the CSP (Salari & Naji-Azimi, 2012). Given an initial feasible solution, the algorithm follows a general destroyand-repair paradigm to improve the tour length. To do this, some vertices are removed from the tour and reassigned, forming a new feasible solution resulting from solving an ILP-based model to optimality. The results show that ILP-based method is more efficient than other available CSP algorithms.

M. Salari et al. / Computers & Industrial Engineering 83 (2015) 244–251

Some researchers worked on a geometric variation of the CSP. In this version, each neighborhood is a compact set in the plane in a way that by intersecting this set, all the demands located at its neighborhood will be covered. Several approximation algorithms have also been proposed for this problem (see, e.g., Arkin & Hassin, 1994; Gulczynski, Heath, & Price, 2006 and Shuttleworth, Golden, Smith, & Wasil, 2008). Several other variants of the CSP have been studied. The covering tour problem (CTP) was defined by Gendreau, Laporte, and Semet (1997), where the set of vertices is divided into three groups, namely those vertices that must be visited by the tour (S1), those that must be covered (S2) and vertices that can be visited or covered (S3). The other assumptions of the CTP are the same as those explained for the CSP. The goal of CTP is to cover all the vertices in S2 by constructing a tour over all the vertices in S1 and possibly a subset of vertices in S3. Gendreau et al. (1997) proposed a heuristic algorithm and an exact branch-and-cut method. A multi-vehicle generalization of the CTP is introduced by Hachicha, Hodgson, Laporte, and Semet (2000), in which we are given a set of vehicles located at a central depot and, the goal is to visit or cover the demand of all available vertices at the minimum cost. In this problem, the classification of the vertices is the same as that defined for the CTP, while the only difference is that we have to cover the demand using more than one vehicle. Golden et al. (2012) proposed some natural generalizations of the CSP, where the assumption of visiting or covering a vertex for only one time is relaxed. They referred to some applications arising in rural healthcare delivery, for which we have to visit or cover a vertex for more than once to satisfy its corresponding demand. The ring star problem (RSP) is another related problem that has applications in telecommunications networks (Labbé, Laporte, Martín, & Salazar González, 2004). In this problem, each vertex is either visited by the tour (incurring a routing cost) or is assigned to the tour (incurring an assignment cost). The goal is to minimize the summation of routing and assignment cost. There is also a variation of the RSP called the median cycle problem, in which there is an upper bound on the maximum allocation cost (Labbé, Laporte, Martín, & Salazar González, 2005). Another problem that is very closely related to CSP is the generalized traveling salesman problem (Fischetti, Salazar González, & Toth, 1997; Karapetyan & Reihaneh, 2012). In this problem, the set of vertices is partitioned into several disjoint clusters. The goal is to find the shortest Hamiltonian tour that visits exactly one vertex from each cluster. In fact, the generalized traveling salesman problem is a special case of the CSP, in which each vertex can be covered only by the vertices of the cluster to which it belongs. In this paper, we propose a new and effective hybrid heuristic algorithm for the CSP. The proposed method combines Ant Colony Optimization (ACO) (Dorigo et al., 1992; Dorigo, Maniezzo, & Colorni, 1996) and dynamic programming technique to improve the solutions’ quality. In addition, we propose a polynomial-size mathematical formulation for the studied problem. The rest of the paper is organized as follows. Section 2 provides a formal description of the problem. The details of the hybrid heuristic algorithm are provided in Section 3. In Section 4, the efficiency of the proposed algorithm is investigated through extensive computational experiments. Finally, Section 5 concludes the paper.

2. Problem statement In this section, we propose a polynomial-size mathematical formulation for the CSP. We are given a directed graph G ¼ ðV; AÞ where V ¼ f1; . . . ; ng is the set of vertices that must be either

245

visited or covered by the tour. In addition, A ¼ fði; jÞji; j 2 Vg is the set of arcs. To each arc ði; jÞ 2 A, a cost cij is associated. For each i 2 V, we denote J i as the set of vertices that can be covered by i and also J 0i as the set of vertices that can cover vertex i. Let  0 j ¼ arg mini2V fjJ 0i jg. We create a copy i of each i 2 J 0j and add it to the set of dummy vertices D. Then the augmented graph G ¼ ðV  ¼ V [ D; A ¼ A [ A0 Þ is constructed. Here, A0 is the set of arcs between D and V, and contains the following arcs. For each 0 0 dummy vertex i 2 D, we introduce only one arc ði ; iÞ with ci0 i ¼ 0 0 0 where i is the original version of i (i is a copy of i). Also, for each 0 0 i 2 V and j 2 D we add an arc ði; j Þ to A0 with cij0 ¼ cij . The variables to model the problem are as follows:

xij ¼

yij ¼





1

if arc ði; jÞ is traversed by the tour;

0

otherwise:

1

if vertex i is assigned to vertex j;

0

otherwise:

For each arc ði; jÞ 2 A let uij be a continuous variable representing the load of the tour after leaving vertex i and before visiting vertex j when the tour passes along arc ði; jÞ. The model reads as follows:

Minimize

X

cij xij

ð1Þ

X yii ¼ 1

ð2Þ

ði;jÞ2A

subject to

i2D

X yij ¼ 1 8i 2 V X

xij þ 

ði;jÞ2A

X ðj;iÞ2A

X

ðj;iÞ2A

xji ¼ 2yii 8i 2 V 

ð4Þ

uij ¼ nyii 8i 2 D

ð5Þ

X yji 8i 2 V

ð6Þ



X

ði;jÞ2A0

uji 

X

ð3Þ

j2J 0i

uij ¼

ði;jÞ2A

j2J i

XX uij ¼ 0

ð7Þ

i2V j2D

uij 6 nxij 8ði;jÞ 2 A yki P ykj  ð1  yii Þ 8k 2 V; i;j 2 J 0k : i < j and i – j – k xij ;yij 2 f0;1g 8ði;jÞ 2 A uij P 0 8ði;jÞ 2 A :

ð8Þ ð9Þ ð10Þ ð11Þ

The objective function (1) is to minimize the tour length. The first constraint (2) is used to assure that the salesman has to visit a subset of the vertices by starting from one of the dummy vertices. Constraint (3) ensures that each vertex is assigned to exactly one routed vertex. For each i 2 V  , constraint (4) represents the relation between routing (x) and allocation (y) variables. In particular, if vertex i is visited by the tour (i.e. yii ¼ 1), there must be two edges incident to vertex i. Sub-tour elimination constraints are represented using constraints (5)–(7). Essentially, constraint (5) imposes that exactly n units of flow has to leave the dummy vertex that is visited by the tour. Constraint (6) stipulates that the load difference of the tour before and after visiting a vertex, say i, is exactly equal to the demand allocated to vertex i. Moreover, since the demand of all the n available customers has to be provided, the flow returned to the set of dummy vertices is equal to zero (given by (7)). The relation between the routing and flow variables is given in (8), in which we are not allowed to send flow along an unvisited arc. Constraints (9) are symmetry defeating constraints and are only introduced to enhance the model by breaking the symmetry caused by y. In fact, each vertex k 2 V can be covered by several vertices, and hence assigning k to any of them will create the same solution. This can

246

M. Salari et al. / Computers & Industrial Engineering 83 (2015) 244–251

lead to a situation where we have several different branch-and-cut nodes that actually represent the same solution. Constraint (9) avoids this situation by enforcing each vertex k to be either covered by itself (ykk ¼ 1) or by the smallest vertex of J0k that is also on the tour. For more information on symmetry defeating, we refer the interested reader to Ghoniem and Sherali (2011) and Jans and Desrosiers (2013). Finally, (10) and (11) represent the decision variables.

3. Solution method We develop a hybrid heuristic algorithm denoted as ‘‘hybrid ACO’’ for the CSP. The proposed algorithm combines ACO and dynamic programming technique to reach high quality results. Moreover, different local search procedures like ‘‘vertex-removal’’, ‘‘vertex-addition’’ and ‘‘3-opt’’ are also incorporated to this purpose. The general framework of the hybrid ACO is given in Algorithm 1. Essentially, the algorithm starts by constructing an initial set of feasible solutions. Following this step, the quality of the solutions is improved by applying different improvement procedures within the general framework of the ACO metaheuristic algorithm. In the following, the details of each procedure are described.

3.2. Ants construction and local pheromone update Let K ¼ f1; 2; . . . ; kg be the set of ants and T l be an ordered set of vertices visited by ant l 2 K. The goal is to construct a set of k ants where each ant is an initial feasible solution for the CSP. To do so, for each l 2 K; T l is initialized with a randomly selected vertex v 0 2 V, i.e. T l ¼ fv 0 g. As long as the solution is not feasible, the algorithm repeats the following procedure in order to construct an initial feasible tour. Assume that u is the last vertex visited by T l and denote AllowedT l as the set of all vertices that are able to cover at least one uncovered vertex of V. The next vertex to be added to T l is selected from AllowedT l . Essentially, for each vertex v 2 AllowedT l the selection of vertices is done by defining an attractiveness score, av , which is calculated as following:

av ¼ suv  ð1=cuv Þb

v 2 AllowedT

l

ð12Þ

where cuv is the length of the arc ðu; v Þ and b is an input parameter. Upon calculating this score for each v 2 AllowedT l , with the probability of q0 , ant l 2 K selects the vertex with the highest attractiveness score as the next vertex to be visited by T l , where, 0 6 q0 6 1 is an algorithm parameter. With a probability of ð1  q0 Þ, the next vertex is randomly selected, from AllowedT l , where the probability of selecting each vertex v 2 AllowedT l is calculated as follows:

pv ¼ P

av

w2AllowedT

Algorithm 1. The general framework of the proposed hybrid ACO algorithm.

8

l

aw

8v 2 AllowedT l

ð13Þ

Letting v to be the selected vertex, the algorithm adds v to T l , i.e. T l ¼ fv 0 ; . . . ; u; v g, and performs the local pheromone update as follows:

suv ¼ ð1  nÞsuv þ ns0

ð14Þ

where n is an algorithm parameter, selected in the range 0 6 n 6 1. For each ant l 2 K and at each step, just one vertex is added to T l and this is iterated until AllowedT l – /. The algorithm stops when we have a set of k initial feasible solutions. 3.3. Improvement procedure Local search procedures that can be used to improve the cost of a given CSP tour, say T, are categorized into three classes as follows:  Moves that change the visiting order of vertices on T.  Moves that replace the position of visited and covered customers.  Moves that change both visited vertices and their visiting order.

3.1. Pheromone initialization We denote as suv the corresponding pheromone of arc ðu; v Þ 2 A. At the beginning of this procedure, for each arc 1 , ðu; v Þ 2 A, we set the initial amount of pheromone as s0 ¼ ncostðT 0Þ where T 0 is an initial feasible tour and costðT 0 Þ is the corresponding objective value of T 0 . To construct T 0 , an adaptation of the nearest neighbor algorithm is used. Essentially, the tour is initialized by v 0 , i.e., T 0 ¼ fv 0 g, where v 0 is a vertex randomly selected from V. Let v t represent the last vertex visited by T 0 ; the next vertex to be visited by T 0 will be the nearest vertex to v t , which is able to cover at least one uncovered vertex. The algorithm terminates when there is no uncovered vertex in the solution. Finally, the initial feasible solution is obtained by connecting the last visited vertex by T 0 to the starting vertex v 0 .

In the hybrid ACO algorithm, one local search procedure from each of the aforementioned classes is used, which are detailed bellow. 3.3.1. Vertex-addition Let T and N ðT Þ denote a feasible CSP tour and the set of all vertices not visited by T, respectively. At first, a subset of vertices, denoted as N 0 with size D is randomly selected from NðTÞ. For all vertices in N 0 , in a random order and once at a time a vertex is selected from the set N 0 and is added into its best position in T, i.e. the position with the cheapest insertion cost. 3.3.2. Vertex-removal A ‘‘redundant vertex’’ is a vertex that if removed from a feasible CSP tour, say T, does not change the feasibility of tour T. Moreover, for a given redundant vertex v 2 T, let g v be the cost decrement achieved by removing v from T. In other words, letting v 0 and v 00

247

M. Salari et al. / Computers & Industrial Engineering 83 (2015) 244–251

be the end points of the edges visiting vertex v in T, we have g v ¼ cv 0 ;v þ cv ;v 00  cv 0 ;v 00 . During the vertex-removal procedure, until there is at least one redundant vertex, say v with g v > 0, the algorithm removes the vertex leading to the maximum decrement in the cost of T. In other words, we remove w ¼ arg max fg v : v is redundantg. 3.3.3. 3-opt heuristic For a given tour T, ‘‘3-opt’’ algorithm iteratively removes three arcs of T and then reconnects the resulting paths by trying all other possible ways to obtain a better tour length (Lin & Kernighan, 1973). 3.3.4. DP-heuristic The goal of this procedure is to improve the length of a given tour by applying dynamic programming technique. Suppose t 2 V is a vertex visited by a given CSP tour, T, and let Cov ðtÞ denote the vertices that become uncovered when t is removed from T. z Moreover, we denote as f ðtÞ ¼ fv 1 ; v 2 ; . . . ; v z g the set that consists of z exclusive vertices whose coverage has the maximum intersecz tion with Cov ðtÞ. Essentially, the elements of f ðtÞ are obtained as follows where ties broken arbitrarily:

Based on the latter recursive formula, Pmþ1 gives the length of 1 the shortest feasible CSP tour that visits each of the levels exactly once in the chosen order. Note that, in case Sðt1 Þ has more than one sequence, we can repeat the aforementioned algorithm for jSðt1 Þj times. Essentially, each time we have to consider one of the sequences of Sðt1 Þ as the only member of L1 and Lmþ1 . The shortest tour among all of the resulting tours will be the best solution of the DP-heuristic. The general framework of the ‘‘improvement procedure’’ is provided in Algorithm 2. The algorithm starts by applying the DPheuristic, vertex-removal and 3-opt procedures, respectively. Then, for a given number of iterations, I, the algorithm tries to further improve the cost of T by applying the designed moves within the improvement procedure. At the beginning of each iteration, vertex-addition and vertex-removal procedures are run in turn. If these two procedures are able to improve the cost of T, 3-opt, DP-heuristic and vertex-removal procedures are run to try to further improve the cost of T. Algorithm 2. The general framework of the improvement procedure.

v 1 ¼ argmaxfjcov ðv Þ \ cov ðtÞj; v 2 V g v 2 ¼ argmaxfjcov ðv Þ \ cov ðtÞj; v 2 V n fv 1 gg .. .

v z ¼ argmaxfjcov ðv Þ \ cov ðtÞj; v 2 V n fv 1 ; v 2 ; . . . ; v z1 gg Let SðtÞ be the set of all sequences with cardinality 1 and 2 (crez ated by combining the members of f ðtÞ) which are able to cover all the vertices of Cov ðtÞ. Each sequence s 2 SðtÞ is denoted as ðus ; v s Þ, where us and v s are the starting and ending vertices of s, respectively. In case of us ¼ v s , the sequence is a singleton. Now assume T ¼ ðt1 ; t2 ; . . . ; tm ; t 1 Þ is a feasible CSP tour and for each 1 6 i 6 m, let Sðt i Þ ¼ fsi1 ; si2 ; . . . ; sihi g. Using dynamic programming technique, the DP-heuristic procedure tries to improve the length of T by replacing each vertex ti 2 T with a sequence of Sðti Þ. The proposed DP-heuristic applies the dynamic programming technique to improve the tour length through finding the best re-arrangement of sequences generated for each ti 2 T. To do so, the algorithm creates m þ 1 levels, represented by L1 to Lmþ1 , where L1 and Lmþ1 represent the corresponding levels for the starting and ending vertex t1 . Moreover, L2 to Lm are those levels corresponding to the vertices t 2 to t m . For the sake of simplicity, we assume that there is only one sequence in Sðt 1 Þ, i.e. Sðt 1 Þ ¼ fs11 g. Consider s11 as the only member of L1 and Lmþ1 . Now, m þ 1 levels are created as in Fig. 1. For each two consecutive levels Li and Liþ1 ; 1 6 i 6 m, the dis  and tance between each two sequences sip ¼ uip ; v ip   iþ1 is defined as follows: ¼ uiþ1 siþ1 q q ; vq

  distance sip ;siþ1 ¼ q

8 <0 :c



if uip ¼ uqiþ1

uip v ip

þ cv i uiþ1 p q



and



 ; v ip ¼ v iþ1 q

otherwise: ð15Þ

where cij is the distance between i and j. Let P kj denote the length of the shortest path linking the sequence s11 (belonging to level L1 ) to the sequence skj (belonging to level Lk ) while passing through only one sequence of each level L2 ; L3 ; . . ., and Lk1 . By induction and assuming P jk1 to be the shortest path from L1 to sjk1 (belonging to Lk1 ) P kj can be calculated as h i P kj ¼ min16i6hk1 P ik1 þ csk1 sk . i

j

3.4. Global pheromone update Let T best be the tour with the minimum objective function. Within this procedure, both evaporation and pheromone deposit are applied to T best as follows:

suv ¼ ð1  qÞsuv þ

q costðT best Þ

8ðu; v Þ 2 T best

ð16Þ

where 0 6 q 6 1 is an algorithm parameter called evaporation rate (Dorigo et al., 1996). 4. Computational results To run the experiments, we used the benchmark instances available in the literature. The instances are generated using the TSPLIB instances (Reinelt, 1991) which are divided into three small, medium, and large size instances. Essentially, the small and medium size instances contain data from n = 51 to 200 vertices. In these instances, each vertex can cover its NC = 7, 9 or 11 nearest vertices. The large size instances contain data from n = 532 to 783 vertices and NC = 3, 5, and 7. Totally, there are 36, 12, and 27 instances, in the small, medium, and large groups, respectively. Four heuristic algorithms are proposed for the CSP, namely the algorithm proposed by Current and Schilling (1989), two methods (i.e. LS1 and LS2) developed by Golden et al. (2012) and an ILPbased improvement procedure denoted as SN proposed by Salari and Naji-Azimi (2012). Among these methods, LS2 and SN are the best methods that outperform the other two approaches.

248

M. Salari et al. / Computers & Industrial Engineering 83 (2015) 244–251

Fig. 1. An illustrative example of the DP-heuristic procedure.

Therefore, to compare our proposed hybrid ACO algorithm, we just consider these two available algorithms. Both SN and LS2 have been implemented in C while the hybrid ACO algorithm has been implemented in C#. Moreover, all algorithms have been tested on a PC with a 1.66 GHz Intel Core Due. Extensive computational tests have been performed to reach a good combination of parameters that are involved in the proposed algorithm. For the small and medium size instances, the termination criterion is 50 iterations of the main loop of the algorithm (see Algorithm 1), without any improvement. Moreover, for the large size instances, both LS2 and SN have used a time threshold of 60 s, as the termination criterion. Therefore, we also run our algorithm for 60 s. The rest of the involved parameters in the proposed hybrid ACO algorithm are given in Table 1. In order to investigate the efficiency of the proposed heuristic, we also solved the MIP model using CPLEX 12.6. Table 2 summarizes the results (columns 3–6). Here, the ‘‘Best’’ column reports the objective function of the best feasible solution obtained by CPLEX. Rows represented by ‘‘—’’ are related to the instances for which CPLEX was not able to achieve a feasible solution within the given time limit. ‘‘Rel-gap’’ is the relative MIP gap reported by CPLEX, which is the gap between upper bound and the smallest lower bound among all the open nodes of the search tree. In addition, ‘‘Gap’’ is the gap between the upper bound, found by CPLEX, and the best solution of all the heuristic algorithms. The results indicate that even the small size instances are very difficult for the solver. In fact, out of 36 instances, the solver was able to solve 9 to optimality while it failed to even find a feasible solution for 4 of the instances within one hour time limit. The overall average of optimality gap is 10.6%, which indicates poor quality of upper bounds found by CPLEX comparing to heuristic algorithms. The impact of symmetry defeating constraints on the overall performance of the solver was also investigated by solving the MIP model once without constraint (9). The results clearly indicated the advantage of using the symmetry defeating constraints. Basically, symmetry defeating constraints enabled the solver to solve three more instances to optimality and reduced the CPU time by 76% for those instances solved by both algorithms to optimality. Concerning the quality of the upper bounds, adding the symmetry defeating constraints led to a 8.8% decrease in ‘‘Rel-gap’’ and a 5.6 % decrease in ‘‘Gap’’. Tables 2–4 report the performance of the hybrid ACO algorithm in comparison with SN and LS2 algorithms for the small, medium and large size instances, respectively. For each instance, five Table 1 Selected values for the hybrid ACO parameters. Parameter

jKj

b

q

n

q0

z

D

J

I

Small and medium instances Large instances

10 10

3 3

0.5 0.5

0.05 0.05

0 0

11 7

40 60

50 50

10 10

different runs of the proposed algorithms are performed and the tables’ headers have the following meanings:  Instance: name of the TSP instance from which the CSP instance is created, followed by some digits that represent the total number of vertices.  NC: number of nearest vertices to each vertex i that can be covered by i.  Best cost: cost of the best solution over five independent runs.  Avg. cost: average solution cost over five different runs.  Deviation: the deviation of the Avg. cost from the best objective function found by all heuristics. The deviation is calculated as Av g:costbest best

 100, where ‘‘best’’ is the cost of the best solution found by all the heuristic algorithms.  TT: the average total computing time over five different runs (see Tables 2 and 3).  TB: the average computing time to find the best solution over five different runs (see Table 4). The last two rows of each table report the average values of each column (Avg.) and the number of the best solutions obtained by each algorithm (NB), respectively. Finally, the best results are given in bold. Table 2 shows that in the small size instances, our proposed algorithm has found the best solutions for 34 out of 36 instances, while LS2 and SN have found the best solutions for 36 and 35 instances, respectively. Taking into account the average performance of the algorithms (Avg.), obtained over five different runs of each method, on average the proposed hybrid ACO has been able to reach the best performance. Moreover, the average deviation of the hybrid ACO, for all of the small size instances is 0.02%, while this value is 0.03% and 0.04% for LS2 and SN algorithms, respectively. Finally, regarding the CPU time, the hybrid ACO is slower than LS2 and is comparable with SN. Table 3 reports the results of different algorithms for the medium size instances. Taking into account the overall average of the best performance (Best cost) of the heuristics over five different runs, SN and hybrid ACO algorithms have almost the same performance, while the hybrid ACO outperforms the other algorithms when considering the average performance (Avg. cost) of different methods. The results of different heuristic algorithms for the large size instances are given in Table 4. Essentially, as the size of instances increases, the hybrid ACO algorithm achieves a larger number of the best solutions. For large size instances, the hybrid ACO has found the best solutions for 23 out of 27 instances, while both LS2 and SN algorithms have found only 2 of the best solutions. Considering the average performance (Avg. cost) of the three different algorithms, our proposed method outperforms the SN and LS2 algorithms by obtaining the best values for 26 instances out of 27, while for only 1 instance does SN have the best average performance over five different runs. Finally, the average deviation of

Table 2 Comparison of different algorithms for the small size instances. Instance

NC

CPLEX Best

LS2

SN

Hybrid ACO

Rel-gap

Time

Gap

Best cost

Avg. cost

Deviation

TT

Best cost

Avg. cost

Deviation

TT

Avg. cost

Deviation

TT

7 9 11

164 159 147

0.00 0.00 0.00

149 220 681

0.00 0.00 0.00

164 159 147

164 159 147

0.00 0.00 0.00

1 1 1

164 159 147

164 159 147

0.00 0.00 0.00

3 2 2

164 159 147

164 159 147

0.00 0.00 0.00

2 2 3

berlin52

7 9 11

3887 3430 3262

0.00 0.00 0.00

140 212 255

0.00 0.00 0.00

3887 3430 3262

3887 3430 3262

0.00 0.00 0.00

1 1 1

3887 3430 3262

3887 3430 3262

0.00 0.00 0.00

2 2 2

3887 3430 3262

3887.4 3430 3262

0.16 0.00 0.00

2 2 2

st70

7 9 11

288 259 250

0.00 0.00 14.51

490 1391 3600

0.00 0.00 1.21

288 259 247

288 259 247

0.00 0.00 0.00

1 2 2

288 259 247

288 259 247

0.00 0.00 0.00

4 4 4

288 259 247

288 259 247

0.00 0.00 0.00

2 6 4

eil76

7 9 11

219 198 177

13.68 22.71 21.90

3600 3600 3600

5.80 7.03 4.12

207 186 170

207 186 170

0.00 0.54 0.00

1 1 1

207 185 170

207 185 170

0.00 0.00 0.00

4 4 4

207 186 170

207 186 170

0.00 0.54 0.00

2 3 2

pr76

7 9 11

50,275 45,387 44,060

0.00 5.71 12.91

2488 3600 3600

0.00 0.09 2.40

50,275 45,348 43,028

50,275 45462.2 43,028

0.00 0.25 0.00

2 2 2

50,275 45,348 43,028

50,275 45,348 43,028

0.00 0.00 0.00

4 4 4

50,275 45,348 43,028

50,275 45,348 43,028

0.00 0.00 0.00

10 6 12

rat99

7 9 11

508 530 473

17.29 40.63 35.09

3600 3600 3600

4.53 16.48 6.53

486 455 444

486 455 444

0.00 0.00 0.00

2 2 2

486 455 444

486 455 444

0.00 0.00 0.00

7 7 7

486 455 444

486 455 444

0.00 0.00 0.00

6 6 6

kroA100

7 9 11

12,762 10,130 12,843

39.07 27.61 49.46

3600 3600 3600

31.92 10.60 44.29

9674 9159 8901

9674 9159 8901

0.00 0.00 0.00

3 2 2

9674 9159 8901

9674 9159 8901

0.00 0.00 0.00

7 7 7

9674 9159 8901

9674 9159 8901

0.00 0.00 0.00

7 9 10

kroB100

7 9 11

– 10,517 –

– 36.97 –

3600 3600 3600

– 13.82 –

9537 9240 8842

9537 9240 8842

0.00 0.00 0.00

3 3 3

9537 9240 8842

9537 9240 8842

0.00 0.00 0.00

7 7 7

9537 9240 8842

9537 9240 8842

0.00 0.00 0.00

4 9 8

kroC100

7 9 11

10,477 10,020 11,226

22.34 30.76 47.91

3600 3600 3600

3.22 11.44 30.05

9723 9171 8632

9723 9171 8632

0.00 0.00 0.00

3 2 2

9723 9171 8632

9723 9171 8632

0.00 0.00 0.00

7 7 7

9723 9171 8632

9724 9171 8632

0.01 0.00 0.00

5 10 5

kroD100

7 9 11

10,316 – –

18.25 – –

3600 3600 3600

7.17 – –

9626 8885 8725

9626 8885 8725

0.00 0.00 0.00

2 3 3

9626 8885 8725

9626 8885 8725

0.00 0.00 0.00

6 7 7

9626 8885 8725

9626 8885 8725

0.00 0.00 0.00

4 5 11

kroE100

7 9 11

10,609 10,719 11,879

14.39 26.10 53.50

3600 3600 3600

9.11 16.88 40.58

10,150 8991 8450

10,150 8991 8450

0.00 0.00 0.00

2 2 2

10,150 8991 8450

10,150 8991 8450

0.00 0.00 0.00

6 7 7

10,150 8992 8450

10,150 8992 8451.2

0.00 0.01 0.01

7 5 10

rd100

7 9 11

3836 4049 3946

24.00 52.00 49.00

3600 3600 3600

10.84 26.77 35.04

3461 3194 2922

3485.6 3194 2922

0.71 0.00 0.00

2 2 2

3461 3194 2922

3493.8 3194 2922

0.95 0.00 0.00

6 6 6

3461 3194 2922

3461 3194 2922

0.00 0.00 0.00

4 6 6

Avg

–

21.12

2867.40

0.04

1.92

8325.67

8326.58

0.03

5.31

8325.72

8325.96

0.02

5.71

NB

9

36

35

34

31

10.06

8325.69

8329.55

9

35

33

M. Salari et al. / Computers & Industrial Engineering 83 (2015) 244–251

Best cost

eil51

249

250

M. Salari et al. / Computers & Industrial Engineering 83 (2015) 244–251

Table 3 Comparison of different algorithms for the medium size instances. Instance

NC

LS2

SN

Hybrid ACO

Best cost

Avg. cost

Deviation

TT

Best cost

Avg. cost

Deviation

TT

Best cost

Avg. cost

Deviation

TT

kroA150

7 9 11

11,423 10,056 9439

11800.0 10062.4 9439.0

3.30 0.06 0.00

4 3 3

11,423 10,056 9439

11423.0 10057.6 9439.0

0.00 0.02 0.00

11 10 9

11,423 10,056 9439

11423.0 10056.0 9439.0

0.00 0.00 0.00

6 10 8

kroB150

7 9 11

11,457 10,121 9611

11491.2 10121.0 9611.0

0.30 0.00 0.00

4 4 4

11,457 10,121 9611

11457.0 10121.0 9611.0

0.00 0.00 0.00

10 10 10

11,457 10,121 9611

11457.0 10121.0 9611.0

0.00 0.00 0.00

9 6 11

kroA200

7 9 11

13,285 11,708 10,748

13666.4 11716.8 10848.6

2.87 0.08 0.94

6 5 5

13,285 11,708 10,748

13327.0 11731.6 10865.6

0.32 0.20 1.09

15 14 13

13,286 11,710 10,760

13286.0 11710.0 10764.2

0.02 0.04 0.15

12 12 15

kroB200

7 9 11

13,100 11,900 10,676

13511.6 11964.8 10809.6

3.53 0.85 1.56

5 5 5

13,051 11,864 10,644

13181.2 11878.4 10656.8

1.00 0.12 0.12

15 14 13

13,051 11,864 10,644

13061.0 11871.2 10644.0

0.08 0.06 0.00

8 12 9

Avg

11127.00

11253.53

1.12

4.42

11117.25

11145.77

0.24

12.00

11118.50

11120.73

0.03

10.29

NB

9

3

12

5

9

11

Table 4 Comparison of different algorithms for the large size instances. Instance

NC

LS2

SN

Hybrid ACO

Best cost

Avg. cost

Deviation

TB

Best cost

Avg. cost

Deviation

TB

Best cost

Avg. cost

Deviation

TB

att532

3 5 7

52,399 42,634 38,186

52760.4 43280.0 38551.8

2.22 2.53 2.15

44 34 38

52,412 42,387 38,016

52794.4 43047.2 38244.0

2.28 1.98 1.33

35 39 31

51,616 42,212 37,741

51841.2 42687.4 38088.4

0.44 1.13 0.92

37 41 44

ali535

3 5 7

1370 1184 1094

1387.0 1201.2 1103.6

1.46 1.45 1.90

25 35 30

1368 1206 1086

1381.4 1210.2 1093.4

1.05 2.21 0.96

11 21 40

1367 1185 1083

1370 1188.4 1084.2

0.22 0.37 0.11

28 30 33

u574

3 5 7

23,330 19,417 16,940

23454.6 19578.6 17090.8

2.02 2.97 2.98

33 42 34

23,286 19,213 16,880

23457.0 19392.8 17056.0

2.03 1.99 2.77

37 28 44

22,990 19,014 16,597

23101.4 19113.6 16789.0

0.48 0.52 1.16

35 46 33

rat575

3 5 7

3755 3062 2667

3783.6 3088.0 2698.6

2.07 1.78 2.41

34 44 35

3754 3104 2645

3771.6 3116.2 2678.8

1.74 2.71 1.66

40 30 24

3707 3034 2635

3726.2 3053.2 2651.0

0.52 0.63 0.61

44 46 45

p654

3 5 7

25,158 23,226 22,121

25206.6 23258.4 22233.4

0.21 0.20 0.51

34 25 26

25,155 23,211 22,126

25206.0 23224.8 22138.6

0.20 0.06 0.08

36 35 18

25,166 23,242 22,125

25182.8 23289.4 22130.8

0.11 0.34 0.04

53 49 47

d657

3 5 7

30,715 26,589 23,429

30987.2 26738.2 23545.0

2.34 3.28 1.50

39 28 43

30,653 26,385 23,320

30878.0 26663.8 23551.6

1.98 2.99 1.52

36 50 41

30,278 25,890 23,198

30498.4 26132.6 23426.8

0.73 0.94 0.99

56 50 54

gr666

3 5 7

2005 1690 1485

2009.0 1695.8 1490.4

1.41 2.47 1.39

24 29 25

2006 1672 1474

2011.4 1685.4 1481.2

1.53 1.84 0.76

23 19 34

1981 1655 1470

1993.6 1665.2 1474.8

0.64 0.62 0.33

44 46 40

u724

3 5 7

25,045 20,563 18,015

25266.8 20799.0 18043.0

2.94 2.43 1.73

52 45 51

25,048 20,575 17,916

25262.4 20844.6 18230.4

2.92 2.66 2.78

51 42 44

24,545 20,305 17,737

24746.0 20543.2 17943.4

0.82 1.17 1.16

52 51 49

rat783

3 5 7

5093 4216 3651

5141.0 4237.2 3678.8

2.39 2.92 1.76

32 46 50

5076 4164 3692

5119.4 4221.6 3700.8

1.96 2.54 2.37

44 45 36

5021 4117 3615

5045.0 4172.0 3638.2

0.48 1.34 0.64

55 52 52

1.98

36.19

16215.93

16350.48

1.81

34.59

16056.52

16169.49

0.65

44.48

2

1

23

26

Avg

16260.70

16381.78

NB

2

0

the hybrid ACO algorithm for all large size instances is 0.65% while this is 1.81% and 1.98% for SN and LS2 algorithms, respectively, which clearly indicates the superiority of the hybrid ACO algorithm. 5. Conclusion We proposed an ant colony optimization algorithm for the covering salesman problem. The proposed algorithm incorporates local search procedures and dynamic programming technique to improve the quality of the solutions. Comparing the results of

the proposed algorithm with the best results available in the literature, clearly indicates the effectiveness of our algorithm, especially in large size instances with more than 500 vertices. For large size instances, the proposed algorithm outperforms the other algorithms by obtaining lower best and lower average costs in a multi-start implementation with five runs. Also with regard to the stability of the solution, the proposed method has less deviation from the best solution. In other words the relative difference between the best and the average CSP tour lengths attained by the developed algorithm turns out to be smaller than that of the other algorithms.

M. Salari et al. / Computers & Industrial Engineering 83 (2015) 244–251

Acknowledgements This work was supported in part by: Research Deputy of Ferdowsi University of Mashhad, under Grant No. 28038 (dated August 20, 2013). This support is gratefully acknowledged. References Arkin, E. M., & Hassin, R. (1994). Approximation algorithms for the geometric covering salesman problem. Discrete Applied Mathematics, 55(3), 197–218. Current, J. R., & Schilling, D. A. (1989). The covering salesman problem. Transportation Science, 23(3), 208–213. Dorigo, M. (1992). Optimization, learning and natural algorithms. PhD thesis, Politecnico di Milano, Italy. Dorigo, M., Maniezzo, V., & Colorni, A. (1996). Ant system: Optimization by a colony of cooperating agents. IEEE Transactions on Systems Man, and Cybernetics – Part B, 26(1), 29–41. Fischetti, M., Salazar González, J. J., & Toth, P. (1997). A branch and cut algorithm for the symmetric generalized traveling salesman problem. Operations Research, 45(3), 378–394. Gendreau, M., Laporte, G., & Semet, F. (1997). The covering tour problem. Operations Research, 45(4), 568–576. Ghoniem, A., & Sherali, H. D. (2011). Defeating symmetry in combinatorial optimization via objective perturbations and hierarchical constraints. IIE Transactions, 43, 575–588. Golden, B. L., Naji-Azimi, Z., Raghavan, S., Salari, M., & Toth, P. (2012). The generalized covering salesman problem. INFORMS Journal on Computing, 24(4), 534–553.

251

Gulczynski, D. J., Heath, J. W., & Price, C. C. (2006). The close enough traveling salesman problem: A discussion of several heuristics. In: F. B. Alt, M. C. Fu, & B. L. Golden (Eds.), Perspectives in operations research. Operations research/computer science interfaces (Vol. 36, pp. 271–283). Hachicha, M., Hodgson, M. J., Laporte, G., & Semet, F. (2000). Heuristics for the multi-vehicle covering tour problem. Computers & Operations Research, 27, 29–42. Jans, R., & Desrosiers, J. (2013). Efficient symmetry breaking formulations for the job grouping problem. Computers & Operations Research, 40, 1132–1142. Karapetyan, D., & Reihaneh, M. (2012). An efficient hybrid ant colony system for the generalized traveling salesman problem. Algorithmic Operations Research, 7, 22–29. Labbé, M., Laporte, G., Martín, I. R., & Salazar González, J. J. (2004). The ring star problem: Polyhedral analysis and exact algorithm. Networks, 43(3), 177–189. Labbé, M., Laporte, G., Martín, I. R., & Salazar González, J. J. (2005). Locating median cycles in networks. European Journal of Operational Research, 160, 457–470. Lin, S., & Kernighan, B. W. (1973). An effective heuristic algorithm for the traveling salesman problem. Operations Research, 21, 498–516. Reinelt, G. (1991). A traveling salesman problem library. ORSA Journal on Computing, 3(4), 376–384. Salari, M., & Naji-Azimi, Z. (2012). An integer programming-based local search for the covering salesman problem. Computers & Operations Research, 39, 2594–2602. Shuttleworth, R., Golden, B. L., Smith, S., & Wasil, E. (2008). Advances in meter reading: Heuristic solution of the close enough traveling salesman problem over a street network. In: B. L. Golden, S. Raghavan, & E. A. Wasil (Eds.), The vehicle routing problem: Latest advances and new challenges. Operations research/computer science interfaces (Vol. 43, pp. 487–501).