Effective metaheuristic algorithms for the minimum differential dispersion problem

European Journal of Operational Research 258 (2017) 829–843

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Discrete Optimization

Effective metaheuristic algorithms for the minimum differential dispersion problem Yang Wang a, Qinghua Wu b,∗, Fred Glover c a

School of Management, Northwestern Polytechnical University, 127 Youyi West Road, 710072 Xi’an, China School of Management, Huazhong University of Science and Technology, 430074 Wuhan, China c OptTek Systems, Inc., 2241 17th Street Boulder, CO 80302, USA b

a r t i c l e

i n f o

Article history: Received 21 December 2015 Accepted 18 October 2016 Available online 22 October 2016 Keywords: Metaheuristics Combinatorial optimization Dispersion problems Tabu search Memetic search

a b s t r a c t This paper presents tabu search and memetic search algorithms for solving the minimum differential dispersion problem. The tabu search algorithm employs a neighborhood decomposition candidate list strategy and a rarely used solution-based tabu memory. Unlike the typical attribute-based tabu list, the solution-based tabu strategy leads to a more highly focused intensification process and avoids tuning the tabu tenure, while employing coordinated hash functions that accelerate the determination of tabu status. The memetic search algorithm incorporates the tabu search procedure within it and makes use of a crossover operator that generates solution assignments by an evaluation mechanism that includes both quality and distance criteria. Experimental results on a benchmark testbed of 250 problems reveal that our tabu search algorithm is capable of discovering better solutions for 179 (71.6%) of the problem instances, while our memetic search algorithm finds better solutions for 157 (62.8%) of the instances, collectively yielding better solutions for 181 (72.4%) of the test problems than recently reported state-ofthe-art algorithms. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Let a set N = {1, . . . , n} and a distance dij (di j ∈ R, di j = d ji ) between any two elements i, j ∈ N. The equitable dispersion problem consists in selecting a subset M ⊆ N of cardinality m (hence |M| = m), to optimize a given equity measure (function) z(M) of dij . Prokopyev, Kong, and Martinez-Torres (2009) proposed several equity measures that define different dispersion problems, including: Max-Mean DP that maximizes the mean dispersion of the selected elements, Generalized Max-Mean DP that extends the mean dispersion by weighting the elements, Min-DiffSum DP that minimizes the difference between the maximum aggregate dispersion and minimum aggregate dispersion, as well as Max-MinSum DP that maximizes the minimum aggregate dispersion among selected elements. In addition, Punnen, Taghipour, Karapetyan, and Bhattacharyya (2014) proposed a balanced quadratic optimization problem that minimizes the difference between the maximum dispersion and the minimum dispersion among selected elements. Equitable dispersion problems are a class of NP-hard combinatorial optimization problems (Prokopyev et al., 2009) and have ap∗

Corresponding author. E-mail addresses: [email protected] (Y. Wang), [email protected] (Q. Wu), [email protected] (F. Glover). http://dx.doi.org/10.1016/j.ejor.2016.10.035 0377-2217/© 2016 Elsevier B.V. All rights reserved.

plications in the selection of homogeneous groups (Brown, 1979a), equity based measures in network flows (Brown, 1979b), dense subgraph identification (Kortsarz & Peleg, 1993), urban public facility location (Teitz, 1968) and tour package planning (Punnen et al., 2014). In particular, equity dispersion problems include applications that are closely related to public facility locations (Barbati & Piccolo, 2015). For instance, the equity measure used in facility location problems is often designed to give members of a population a fair exposure to desirable facilities such as schools, hospitals and parks, or to undesirable facilities such as nuclear power plants and oil storage tanks. Another example is to locate franchise stores in a designated region using an equity measure that assures each franchise shop will be located roughly a similar distance from the others, as a means to provide an equitable coverage of the region. Among those problems that derive from an efficiency measure, the maximum dispersion problem and the maximum diversity problem receive the preponderance of attention in the literature (see, e.g., the survey in Martí, Gallego, Duarte, & Pardo (2013)). For the problems based on equity measures, the minimum dispersion problem (Max-Min DP) is the most widely studied equitable dispersion problem so far, addressed by a wide variety of approaches including polynomial time algorithms for special cases (Chandrasekaran & Daughety, 1981), branch and bound (Erkut, 1990), integer programming (Kuby, 1987), local search heuristics

830

Y. Wang et al. / European Journal of Operational Research 258 (2017) 829–843

(Kincaid, 1992; Porumbel, Hao, & Glover, 2011), GRASP with Path Relinking (Resende, Martí, Gallego, & Duarte, 2010) and clique based heuristics (Della Croce, Grosso, & Locatelli, 2009). The MaxMean DP has additionally been tackled via a GRASP and path relinking algorithm by Martí and Sandoya (2013). In this paper, we focus on the Min-DiffSum DP, whose equity   measure is given by z(M ) = maxi∈M j∈M di j − mini∈M j∈M di j . An equivalent binary optimization formulation arises as follows. Let xi = 1 for elements i ∈ M and xi = 0 for elements i ∈ NࢨM. Then the Min-DiffSum DP consists of assigning values 0 or 1 to each variable in x = {x1 , x2 , . . . , xn }, subject to m variables receiving value 1 and the other n − m variables receiving value 0, to minimize the func  tion f (x ) = maxi∈N:xi =1 j∈N di j x j − min j∈N:x j =1 i∈N di j xi . Alternatively, let Li and Ui represent lower and upper bounds  on the value of  j ∈ M dij , i.e. Li = j∈M min{di j , 0} and Ui =  j∈M max{di j , 0}. Then a mixed linear 0–1 programming model of the Min-DiffSum DP can be formulated as follows (Prokopyev et al., 2009).

min t = r − s 

s.t. r ≥ s≤ 



(1) di j x j − Ui (1 − xi ) + M− (1 − xi ), i ∈ N

(2)

j∈N, j=i

di j x j − Li ( 1 − xi ) + M + ( 1 − xi ), i ∈ N

(3)

j∈N, j=i

xi = m

(4)

i∈N

xi ∈ {0, 1}, i ∈ N

(5)

between a local search phase and a shaking phase. With a randomized initial solution, the local search phase employs swap moves for neighborhood exploration and investigates both the first improvement and best improvement strategies. Once the local search phase is trapped in a local optimum, the shaking phase is triggered by randomly performing a sequence of swap moves to launch the search in a new region. Although the proposed VNS algorithm is quite simple, computational assessment shows that it outperforms the more complicated hybrid metaheuristics proposed in Duarte et al. (2014). Aringhieri, Cordone, and Grosso (2014) devised a collection of innovative two-phase heuristics for the Max-MinSum DP and Min-DiffSum DP, which alternate between a constructive phase and an improvement phase. The constructive phase removes nonpromising vertices or edges to quickly identify a clique of cardinality m or terminates on determining that the construction cannot find such a clique. The improvement phase modifies the subgraph inherited from the constructive phase by means of a tabu search algorithm with dynamically adjusted tabu tenure. A diversification strategy is included to drive the search to cover a large area. Experiments on six sets of benchmarks disclose the high efficacy of the proposed heuristics. The new tabu search and memetic search algorithms We propose two new algorithms for the Min-DiffSum DP problem, a tabu search (TS) algorithm and a memetic search (MS) algorithm, with the following features. •

where M+ is an upper bound on the Ui values, M− is a lower bound on the Li values. State-of-the-art Min-DiffSum DP algorithms We briefly sketch the history of recent efforts to solve the MinDiffSum DP problem as a foundation for understanding the algorithms that constitute the current state-of-the-art for this problem class. Each of the methods that follows obtains high quality results and each improves upon the algorithms that preceded it. Prokopyev et al. (2009) proposed a greedy randomized adaptive search procedure (GRASP) which alternates between a solution construction phase and a local search phase. The solution construction phase selects an element according to a randomized greedy function to gradually enlarge the current solution until the number of elements in the solution reaches the cardinality m. The local search phase uses a descent heuristic based on swap moves. In addition, Prokopyev solved the mixed integer programming model of the problem by the general CPLEX optimization software. Experiments on small instances indicate that the proposed GRASP algorithm consumes less time than CPLEX and obtains high quality solutions. Duarte, Sánchez-Oro, Resende, Glover, and Martí (2014) present several hybrid metaheuristics combining GRASP and path relinking algorithms. Two initial solution constructions are proposed, mainly differing in the way the restricted candidate list of GRASP is constructed. Three local search procedures utilizing swap moves are developed that rely on different orders for scanning the elements. The performance of the algorithm is further enhanced by using variable neighborhood search as an improvement method. Both interior path relinking and exterior path relinking schemes are investigated, in which greedy and random moves are employed to generate solutions on the path connecting the starting solution and guiding solution. Experimental analysis indicates that the hybrid variant combining the variable neighborhood search with exterior path relinking performs the best among all the proposed variants. ´ Todosijevic, ´ and Uroševic´ (2016) present a basic Mladenovic, variable neighborhood search (VNS) algorithm, which alternates

•

•

•

The tabu search algorithm employs a neighborhood decomposition candidate list strategy to exclude the examination of unpromising moves, and to promote a more extensive neighborhood exploration in a given time period. For most instances, our candidate list strategy consumes on average one third the amount of time required by the typical approach of exploring the full neighborhood, while achieving slightly better solution quality. The tabu search method additionally employs a solution-based determination of tabu status which incorporates dedicated hash functions to avoid the need of tuning a tabu tenure. We show our TS solution-based strategy achieves an extensive neighborhood exploration and is more effective than the attribute-based tabu list implemented in Aringhieri et al. (2014). We introduce the first adaptation of memetic search for handling the Min-DiffSum DP, incorporating a distance-quality guided crossover operator to create offspring solutions with high quality and good diversity. The proposed memetic search proves particularly effective for two of the benchmark sets studied. Experimental assessment on seven sets of benchmarks with a total of 250 instances discloses that our tabu search algorithm finds better solutions for 179 (71.6%) of the instances and our memetic algorithm finds better solutions for 157 (62.8%) of the instances, performing statistically better than the state-of-theart algorithms for most of the benchmark sets.

The rest of the paper is organized as follows. Sections 2 and 3 describe the proposed tabu search and memetic search algorithms, respectively. Section 4 presents experimental results and comparisons with state-of-the-art algorithms in the literature. Section 5 analyzes the influence of the essential components of our algorithms on their performance. Concluding remarks are given in Section 6. 2. Tabu search Following the general tabu search framework (Glover & Laguna, 1997; Glover, Ye, Punnen, & Kochenberger, 2015), our tabu search

Y. Wang et al. / European Journal of Operational Research 258 (2017) 829–843

Algorithm 1 Outline of the tabu search algorithm. 1: 2: 3: 4: 5: 6:

7: 8: 9: 10: 11: 12: 13: 14:

Input: a given solution M with its solution value f (M ) Output: the locally optimal solution M∗ found by tabu search with its solution value f (M∗ )  Initialize Iter = 0, i = j∈M di j , i ∈ N repeat Identify the swap moves in the neighborhood decomposition candidate list (see Section 2.2) Determine the tabu status of each neighbor solution obtained by performing the swap move with respect to the solution M (see Section 2.3) Identify the best non-tabu solution in terms of the objective function value, say this move corresponds to (u, v ) Update the current solution by making M = M\{u} ∪ {v} and record the objective function value f (M ) Update i , i ∈ N according to Eq. (6) if f (M ) < f (M∗ ) then M ∗ = M , f (M ∗ ) = f (M ) end if Iter = Iter + 1 until Iter ≥ MaxIter

algorithm (see Algorithm 1) is specifically adapted to the MinDiffSum DP by integrating a neighborhood decomposition candidate list and a solution-based tabu strategy to record solutions encountered. The main components of our tabu procedure are detailed in the following sections. 2.1. Neighborhood evaluation The solution of the Min-DiffSum DP is represented as an m-element subset M of the set N. To transit from the current solution to a neighbor solution, we employ a customary swap move that interchanges u ∈ M and v ∈ N \ M. The resulting solution remains feasible as the size of the subset M is unchanged. In order to quickly evaluate the objective function value of each neighbor solution, we apply a streamlined technique to reduce the complexity of this calculation from O(m2 ) to O(n). For a given so lution M, let i = u∈M dui be the sum of distances between the element i ∈ N and all the elements in M, let max = maxu∈M u and min = minu∈M u . Then the objective function of the solution M, denoted by f(M), is computed as f (M ) = max − min . Once a move (u, v ) is performed that replaces the element u ∈ M by the element v ∈ N \ M, then i can be updated using the following equation:

 i − dui , if i = v if i = u i = i + d v i , i − dui + dvi , otherwise

(6)

2.2. Neighborhood decomposition candidate list A simple way to explore the search area with tabu search is to start with any initial m-subset M and successively examine all swaps of the form that interchange an element in M with another element in NࢨM in order to pick the non-tabu swap that yields the best objective function value. Such a full neighborhood examination was used in the leading state-of-the-art method (Aringhieri et al., 2014) and further explored in Duarte et al. (2014). One advantage of this swap move is that it maintains solution feasibility. However, this full neighborhood examination entails evaluating m × (n − m ) swaps. Moreover, the time complexity to calculate the objective gain of a swap move is O(m) even if the streamed technique to accelerate the evaluation procedure is employed. Hence, the time complexity for each tabu search iteration with a full

831

neighborhood examination is O(nm2 ), which is too computationally expensive. To improve the computational efficiency of our tabu search procedure, we devise a candidate list strategy based on neighborhood decomposition (see, e.g., Glover, Taillard, & de Werra, 1993) which limits the swap moves to be examined to come from two specifically identified subsets X ⊆ M and Y ⊆ NࢨM such that at least one of |X| and |Y| is somewhat smaller than its source set, and the resulting neighborhood contains the most promising solutions within the full neighborhood induced by M and NࢨM. From the computation of the objective value f(M) of the solution M, it’s quite likely that swapping an element v ∈ N \ M with v > max or v < min with another element u ∈ M tends to produce a worse solution than the one currently considered. Based on this observation, we use the following probability formula to determine if a vertex v ∈ N \ M should be included in the subset Y:

P [v is included in Y ]



=

1, e(v −min )/(0.1∗min ) , e(max −v )/(0.1∗max ) ,

if min ≤ v ≤ max , if v < min , if v > max .

The reason for specifying a small probability of permitting vertices with v > max and v < min to be involved in a swap move is to avoid the occurrence of an empty subset Y during the run, especially for some small instances with m = 5 or 10 and thus to avoid an unusual termination of the run. Moreover, this probability equation signifies that under the case v < min or v > max , as v becomes closer to max or min , the vertex v becomes increasingly preferable and thus is assigned a higher probability to be included in the subset M. Finally, in order not to make results vary too much for different instances and keep the equation simple, we use the factor 0.1 in the probability equation even if fine-tuning may produce better results. Considering that this probability mechanism to build Y is not sufficient to guarantee that Y is never empty, we examine the swap moves in the full neighborhood instead of the constrained neighborhood in case this mechanism fails. On the other hand, we simply set subset X to be equal to M. To obtain a neighboring solution M from M, we swap one element u ∈ X with another element v ∈ Y . All possible swap moves induced by X and Y define our neighborhood decomposition candidate list CN(M) which is therefore given by,

CN (M ) = {M |M = M\{u} ∪ {v}, u ∈ X, v ∈ Y }. It is noteworthy that the candidate list does not necessarily contain the best swap move, which is illustrated by a small example shown in Fig. 1. This example is a graph with n = 5 vertices A, B, C, D, E and m = 3 for the subset size. The current solution is S = {A, B, C } and its six neighbor solutions are S1 = {D, A, B}, S2 = {D, A, C }, S3 = {D, B, C }, S4 = {E, A, B}, S5 = {E, A, C } and S6 = {E, B, C } after performing swap moves. We first calculate v for each vertex v in S, giving A = dAB + dAC = 6, B = dBA + dBC = 3, C = dCA + dCB = 5. With max = 6 and min = 3, the solution S gets an objective value f (S ) = 3. Since min ≤ E = dEA + dEB + dEC = 6 ≤ max and D = dDA + dDB + dDC = 12 > max , only the vertex E will be considered for a swap move with vertices in the set S determined by the neighborhood decomposition candidate list. When swapping the vertex E with any vertex in S = {A, B, C }, the objective values of the resulting neighbor solutions are f (S4 ) = f (S5 ) = f (S6 ) = 3 and thus the solution S will not be improved. However, the best neighborhood solution provided by the full neighborhood is S2 = {D, A, B} with an objective value f (S2 ) = 2, which leads to an improvement over the solution S. It can be appropriate to employ a process that does not always perform the best swap move as a means to help diversify the

832

Y. Wang et al. / European Journal of Operational Research 258 (2017) 829–843

D 2 E

3 5

Neighbor Solutions

4 1

1 4 4

A

C

S1={D, A, B}

S4={E, A, B}

S2={D, A, C}

S5={E, A, C}

S3={D, B, C}

S6={E, B, C}

1 2

S={A, B, C} B

Fig. 1. An example of candidate list not containing the best swap move.

search, which is crucial to the search performance. Computational experience indicates that the time required per iteration with the neighborhood decomposition strategy is on average one third as long as that required by using the full swap neighborhood. A detailed analysis can be found in section 5.4.

Given a solution x = {x1 , x2 , . . . , xn } where xi = 1 if i ∈ M and xi = 0 otherwise, we associate each element i with a pre-computed weight wi = iγ  where γ is a positive real number and iγ  is the integer part of iγ . Then we propose the following type of hash functions:

2.3. Solution-based determination of tabu status

h (x ) =

n 

wi × xi

(7)

i=1

Tabu search frequently employs an attribute-based definition of tabu status to identify moves that are excluded from consideration for a duration called the tabu tenure. A solution attribute in the simplest case can be a vertex or edge or the value of a variable, which is used to perform moves that transform one solution into another. While attaching a tabu status to such a simple attribute proves effective in many applications, there also arise applications in which more complex attributes are relevant. One of the key ways to generate more complex attributes for problems where solutions can be expressed as vectors of integer variables (and particularly 0–1 variables) is to encode a solution by creating a special function which often is expressed as an integer combination of the component variables. The methods for doing this make use of socalled hash functions, and a strategy for determining tabu status by reference to hash functions is called a solution-based strategy, since the composite attribute created is a function of all the variables in the solution. (See, e.g., Section 7.4 of Glover and Laguna, 1997.) In short, where hash functions can be used expeditiously, and are able to sufficiently differentiate among different solutions, they have the advantage of freeing the algorithm designer from the onerous work of tuning tabu tenure, because they make it unnecessary to limit tabu status to a specified time duration. Done properly, a solution-based determination of tabu status can be relatively easy to implement, since the main operations consist of maintaining hashing vectors that record the values computed by hash functions, and of comparing values in these vectors to corresponding hash values of solutions that are candidates to be generated to determine whether these solutions have been previously visited. As pointed out by Woodruff and Zemel (1993), a proper hash function capable of efficiently recording a large number of solutions and checking for duplicates should comply with three goals, one of which is to yield a low probability of collision. A collision occurs when two different vectors are encountered with the same hash function value and such encounters may occur frequently due to the well known “birthday paradox” (Carlton & Barnes, 1996). To effectively reduce the collision rate, we propose in this work to use multiple hash functions to keep track of previously visited solutions. (Using our hash function definition below, we find that using two hash functions works particularly well in the present context.)

For a solution x and its hash function value h(x), we can quickly calculate the hash function value of its neighboring solution x resulting from interchanging elements u ∈ M and v ∈ N \ M (i.e. setting xu = 0 and xv = 1) as follows:

h ( x ) = h ( x ) − wu + wv

(8)

The time complexity to compute the hash function value of any neighbor solutions is O(1), which is computationally rather cheap compared to computation of the objective function (Eq. (7)) which requires a time complexity O(n). Using different weighting methods for wi , we can obtain different hash functions. In our experiment, we set w1i = i1.3  and w2i = i1.8  to obtain the two hashing functions h1 (x) and h2 (x). The rationale underlying these parameter settings is discussed in Section 5.1. Our method to determine if a solution should be classified as tabu or not then operates as follows. First, we associate each hash function with a hashing vector of length L (|L| = 107 ). Initially the hashing vectors are filled with 0’s. For each candidate solution x, the hashing values h1 (x) and h2 (x) are computed and we identify the index position i1 in the hashing vector H1 where i1 = h1 (x ) mod L and the index position i2 in the hashing vector H2 where i2 = h2 (x ) mod L. Then we look at the index position i1 in the hash list H1 . If H1 (i1 ) = 0, we readily confirm that the candidate solution has not been visited. Otherwise, we additionally look at the index position i2 in the hashing vector H2 . If now H2 (i2 ) = 1 (as well as H1 (i1 ) = 1), we assume that the candidate solution was previously visited and mark it as a tabu solution. If H2 (i1 ) = 0 (as well as H1 (i1 ) = 1), we assume that this candidate solution has not been visited. Finally, for each visited solution in the search trajectory, we fill in the corresponding positions in the hashing vectors with 1. As demonstrated in Section 5.3, our two hash functions can effectively determine tabu status of solutions and lead to a much lower collision rate than using a single hash function. 2.4. Tabu search algorithm Our proposed tabu search algorithm repeatedly carries out the following iterative procedure until reaching the pre-determined maximum number of iterations. At each iteration, the algorithm identifies a subset of swap moves by our candidate list strategy

Y. Wang et al. / European Journal of Operational Research 258 (2017) 829–843

based on neighborhood decomposition. Then solutions in the candidate list are categorized as tabu or non-tabu by the solutionbased tabu determination. The non-tabu solution with the smallest objective function value is chosen as the solution for the next iteration, breaking ties randomly. The pseudo-code of our tabu search method is shown in Algorithm 1. During each tabu search iteration, determining the candidate lists X and Y is of complexity O(n), the tabu status of each neighbor solution is of complexity O(1), the objective value of each neighbor solution is of complexity O(m), and updating the vector  is of complexity O(n). Therefore, the total time complexity of each iteration in our tabu search algorithm using the constrained neighborhood is bounded by O(n ) + O(m2 |Y | ). 3. Memetic search Memetic search (MS) is a hybrid method that couples a population based genetic algorithm with a local search procedure for solution improvement (Hao, 2012, Chap. 6; Neri, Cotta, and Moscato, 2011; Wu and Hao, 2015; Xu, Lü, Yin, Shen, and Buscher, 2014). We employ the general MS framework to design our memetic search algorithm for the Min-DiffSum DP, which is described in Algorithm 2. Algorithm 2 Outline of the memetic search algorithm. 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11:

P = {M1 , . . . , M|P| } ←− Init_Pop() Record the best solution M∗ with its objective function value f (M ∗ ) repeat M p1 , M p2 ←− Select_Parents(P ) Mc ←− Crossover_Operator(M p1 , M p2 ) Mc ←− Tabu_Search(Mc ) P = {M1 , . . . , M|P| } ←− Update_Pop(M1 , . . . , M|P| , Mc ) if f (Mc ) < f (M∗ ) then M∗ = Mc end if until the stopping criterion is satisfied

The first step of MS is to create an initial population P with high quality and good diversity. Then two solutions are chosen from P as parents to produce an offspring solution by use of a crossover operator, whereupon an improving procedure is launched to optimize the newly generated offspring solution. As in most MS algorithms, we use tabu search as the improving procedure, in the present case drawing on the TS algorithm described in the preceding section. After applying the improving method, a population updating rule decides whether the resulting (potentially) improved solution should be inserted into the population and, if so, identifies which existing individual in the population should be replaced. The procedure repeats until a given stopping criterion is met. 3.1. Population maintenance The initial population of |P| solutions is generated as follows. First, a solution M is randomly created by randomly selecting and adding m elements to M. At each step where such a solution M is generated, it is passed to tabu search for improvement. The resulting improved solution is added to the population if it is not identical to any solution currently in the population. This procedure is repeated until the population is filled up with |P| solutions (|P | = 10 in our experiments). Once a population is complete, the crossover operator is used to produce an offspring solution from selected parent solutions. We use the tournament selection strategy to choose each parent

833

solution for producing the offspring solution as follows. A subpopulation of two solutions is randomly selected from the population P and the best solution in the subpopulation becomes one of the parent solutions. In the same way a second parent solution is chosen. Note that the second subpopulation should not include the best solution in the first subpopulation to avoid identical solutions chosen from the two subpopulations. The second parent solution is then joined with the first by the crossover operator to give an offspring solution, which in turn is submitted to tabu search for solution improvement. Whether the improved offspring solution qualifies to be inserted in the population depends on the population updating mechanism. In our approach, if the improved solution has a better solution quality than the worst solution and is not a duplicate of any solution in the population, then we add it to the population to replace the worst solution. This simple population updating mechanism has been shown quite effective in Peng, Lü, and Cheng (2015); Wang, Hao, Glover, and Lü (2014); Wang, Lü, Hao, and Glover (2012).

3.2. Crossover operator In a memetic search algorithm, the determination of a crossover operator strongly influences the performance of the search. To design an effective operator, it is often preferable to integrate some pertinent properties (“building blocks”) of the problem being solved and devise a recombination mechanism to preserve these properties from parents to offspring (Hao, 2012, Chap. 6). To identify building blocks for the Min-DiffSum DP, we carried out a detailed analysis of locally optimal solutions (see Section 5.5) which disclosed that in some instances, especially those with a large m/n value, high quality solutions share a significant number of common elements (This is often a feature of high quality solutions in a number of settings). Generally speaking, if there are common elements shared by many high quality solutions, there is a reasonably strong chance that they are a part of globally optimal solution. Therefore, as a starting point for our crossover operator, we require it to preserve shared elements among the chosen parent solutions when creating an offspring. Upon first creating a partial solution by extracting common elements shared by the two selected parents, we then enlarge this solution by gradually adding elements using a quality-and-distance based construction strategy. Specifically, given two parents M1 ⊂ N and M2 ⊂ N, the partial offspring is first constructed as Mc = M1 ∩ M2 . Then, the elements in M1 ࢨMc and M2 ࢨMc are considered in turn until the size of Mc reaches m, to guarantee the offspring solution being equally far from each parent solution. The preference of adding an element u into the partial solution Mc is evaluated  based on a greedy function g(u|Mc ) = maxi∈Mc ∪{u} j∈Mc ∪{u} di j −  mini∈Mc ∪{u} j∈Mc ∪{u} di j . To be specific, we first examine all the el-

ements in M1 ࢨMc to identify the element with the smallest g(u|Mc ) value and move it from M1 to Mc . Next, we identify the element in M2 ࢨMc with the smallest g(u|Mc ) value and displace it from M2 to Mc . The indicated construction procedure finally creates an offspring solution with high quality and good diversity. Our crossover operator is especially effective for instances with a large ratio of m/n, which appears to be an indicator of problems in which high quality solutions often share a large number of elements. During each MS iteration, choosing the parent solutions is of complexity O(1), the solution recombination mechanism of complexity O(nm2 ), the tabu search of complexity O((n + m2 |Y | )MaxIter ), the population updating of complexity O(n). Hence, the total complexity for each MS iteration is O(nm2 + (n + m2 |Y | )MaxIter ).

834

Y. Wang et al. / European Journal of Operational Research 258 (2017) 829–843

Table 1 Benchmarks properties. Benchmark

No. of instances

n

m

dij

MDG-a Duarte and Martí (2007) MDG-b Duarte and Martí (2007) MDG-c Duarte and Martí (2007) SOM-b Silver, Ochi, and Martins (2004) GKD-b Glover, Kuo, and Dhir (1998) GKD-c Glover et al. (1998) APOM Aringhieri et al. (2014)

60 40 20 20 50 20 40

[50 0,20 0 0] [50 0, 20 0 0] 30 0 0 [10 0,50 0] [25,150] 500 [50, 250]

[0.1n, 0.4n] 0.1n [0.1n, 0.2n] [0.1n, 0.4n] [2,45] 50 [0.2n, 0.4n]

Random real value in [0, 10] Random real value in [0, 10 0 0] Random integer in [0, 10 0 0] Random integer in [0, 9] Euclidean Euclidean Euclidean or Random

4. Computational results and comparisons 4.1. Experimental settings In order to evaluate our proposed TS and MS algorithms, we use seven sets of benchmarks1 , 2 with a total of 250 instances, which are also used in the best previous state-of-the-art algorithms, denoted by Construction and Improvement Heuristics (CIH) (Aringhieri et al., 2014), Evolutionary Path Relinking (EPR) (Duarte et al., 2014) and Variable Neighborhood Search (VNS) (Mladenovic´ et al., 2016), for performance evaluation. The properties of the benchmarks, including the number of elements n, the cardinality of selected elements m and the distance dij between elements i and j are shown in Table 1. Note that we put the 20 DM1A instances from Aringhieri et al. (2014) into the MDG-a benchmarks given that they share similar characteristics and name them from MDG-a_41_n50 0_m20 0 to MDG-a_60_n50 0_m20 0. Our TS and MS algorithms are programmed in C++ and compiled using GNU g++ with - O3 flag on a Xeon E5440 with 2.83 gigahertz CPU and 8 gigabytes RAM. To make fair comparisons with state-of-the-art algorithms CIH, EPR and VNS, we re-implement each of them in the same programming language C++ and run them on our computing platform. For each instance, our TS algorithm along with the reference algorithms CIH, EPR and VNS are given 20 independent runs. The stopping condition of the EPR algorithm is adopted as the uniform criterion for all tested algorithms to terminate each run, which is set to be the total time of 100 GRASP iterations plus 45 Path Relinking iterations, subject to a maximum allowed computational time of n seconds (n being the number of vertices) (Duarte et al., 2014). Column 2 in Tables 2–8 shows the time limit for each instance during each run. Our MS algorithm is run just once with the time limit equaling to the total time of the 20 independent TS runs. Within the MS algorithm, the tabu search procedure for solution improvement is limited to 100 iterations for 20 MDG-a instances with n = 500, m = 200 and the SOM-b benchmarks, and to 10 0 0 iterations for the other benchmarks. Our policy of choosing different iteration limits for applying TS as an improvement method for different benchmarks is motivated by the goal of exploiting the tradeoff between search intensification and diversification. Generally speaking, a larger number of iterations helps TS to conduct a more focused and deeper intensification process in a specific search region, while frequent applications of the crossover operator facilitate the generation of more diverse solutions in a broader search area. 4.2. Experimental results Tables 2–8 respectively present results of our proposed TS and MS algorithms, as well as the algorithms CIH, EPR and VNS over each set of benchmarks. For each algorithm, we report the best solution value Bst and the average solution value Avg. In addition, 1 2

http://www.optsicom.es/mdp. http://www.di.unito.it/∼aringhie/benchmarks.html.

Column 2 shows the given time limit Time (in seconds) for each run, which is the same for each algorithm. As shown in Tables 2–8, our proposed TS and MS algorithms are collectively able to attain better solutions (marked in bold) than the compared algorithms for a total of 181 instances, accounting for 72.4% of the tested 250 instances. Fig. 2 shows the bar plots for the number of better solutions obtained by TS and MS over each set of benchmarks, respectively. From this figure, we first observe that TS by itself is able to obtain better solutions for a total of 179 instances, encompassing 52 of the 60 MDG-a instances, 37 of the 40 MDG-b instances, 20 of the 20 MDG-c instances, 14 of the 20 SOM-b instances, 13 of the 50 GKD-b instances, 20 of the 20 GKD-c instances and 23 of the 40 APOM instances. The MS algorithm is able to get better solutions for a total of 157 instances, encompassing 48 of the 60 MDG-a instances, 29 of the 40 MDG-b instances, 20 of the 20 MDG-c instances, 14 of the 20 SOM-b instances, 7 of the 50 GKDb instances, 18 of the 20 GKD-c instances and 21 of the 40 APOM instances. Moreover, we compute for each instance the deviation of the best solution values found by each algorithm from the best known results among the five algorithms, computed as devbst % = (Bst − BK R )/BK R ∗ 100% and summarize the average deviations over each set of benchmark instances in Table 9. As shown in Table 9, TS achieves the smallest average deviation of 2.16% from the best known results, slightly less than the 2.51% of MS, 3.6 times less than the 9.86% of CIH, 10.6 times less than the 25.00% of VNS and 18.5 times less than the 42.15% of EPR in terms of best solution quality. By the same token, we calculate the deviation of the average solution values found by each algorithm (except the MS algorithm) from the best known results among the five algorithms, computed as devavg % = (Avg − BK R )/BK R ∗ 100% and summarize the average deviations over each set of benchmark instances in Table 10. As can be seen from Table 10, TS achieves the smallest average deviation of 12.08% from the best known results, 0.7 times less than the 20.34% of CIH, 2.4 times less than the 40.52% of VNS and 3.9 times less than the 59.67% of EPR in terms of average solution quality. It is noteworthy that for instances with n = 20 0 0 or 30 0 0 that use a computational time of n seconds as the stopping criterion, our reported results for the EPR algorithm are much better than the results for this algorithm reported in Duarte et al. (2014). This is mainly due to the fact that the EPR program is tested on our computing platform and compiled with the optimization flag - O3, which achieves a 2.77 computational speed improvement over the computing platform used in Duarte et al. (2014), i.e. we have implicitly given a longer time limit to the EPR algorithm in our experiment. For the VNS algorithm, we report results of using the first improvement search strategy because it leads to better overall performance. For small instances, our reported results may not be as good as the previously reported results in Mladenovic´ et al. (2016) because of the shorter time limit used in our experiments. For large instances, we use the optimization compiling flag - O3 and thus obtain better results for the VNS algorithm than reported in Mladenovic´ et al. (2016) using the same time limit. For the CIH

Y. Wang et al. / European Journal of Operational Research 258 (2017) 829–843

835

Table 2 Computational comparisons on the MDG-a benchmark instances (time in seconds). Instances

MDG-a_1_n500_m50 MDG-a_2_n500_m50 MDG-a_3_n500_m50 MDG-a_4_n500_m50 MDG-a_5_n500_m50 MDG-a_6_n500_m50 MDG-a_7_n500_m50 MDG-a_8_n500_m50 MDG-a_9_n500_m50 MDG-a_10_n500_m50 MDG-a_11_n500_m50 MDG-a_12_n500_m50 MDG-a_13_n500_m50 MDG-a_14_n500_m50 MDG-a_15_n500_m50 MDG-a_16_n500_m50 MDG-a_17_n500_m50 MDG-a_18_n500_m50 MDG-a_19_n500_m50 MDG-a_20_n500_m50 MDG-a_21_n20 0 0_m20 0 MDG-a_22_n20 0 0_m20 0 MDG-a_23_n20 0 0_m20 0 MDG-a_24_n20 0 0_m20 0 MDG-a_25_n20 0 0_m20 0 MDG-a_26_n20 0 0_m20 0 MDG-a_27_n20 0 0_m20 0 MDG-a_28_n20 0 0_m20 0 MDG-a_29_n20 0 0_m20 0 MDG-a_30_n20 0 0_m20 0 MDG-a_31_n20 0 0_m20 0 MDG-a_32_n20 0 0_m20 0 MDG-a_33_n20 0 0_m20 0 MDG-a_34_n20 0 0_m20 0 MDG-a_35_n20 0 0_m20 0 MDG-a_36_n20 0 0_m20 0 MDG-a_37_n20 0 0_m20 0 MDG-a_38_n20 0 0_m20 0 MDG-a_39_n20 0 0_m20 0 MDG-a_40_n20 0 0_m20 0 MDG-a_41_n50 0m20 0 MDG-a_42_n50 0m20 0 MDG-a_43_n50 0m20 0 MDG-a_44_n50 0m20 0 MDG-a_45_n50 0m20 0 MDG-a_46_n50 0m20 0 MDG-a_47_n50 0m20 0 MDG-a_48_n50 0m20 0 MDG-a_49_n50 0m20 0 MDG-a_50_n50 0m20 0 MDG-a_51_n50 0m20 0 MDG-a_52_n50 0m20 0 MDG-a_53_n50 0m20 0 MDG-a_54_n50 0m20 0 MDG-a_55_n50 0m20 0 MDG-a_56_n50 0m20 0 MDG-a_57_n50 0m20 0 MDG-a_58_n50 0m20 0 MDG-a_59_n50 0m20 0 MDG-a_60_n50 0m20 0

Time

65 64 66 65 65 65 65 65 65 64 64 65 64 66 66 65 65 65 64 64 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 500 500 500 500 500 500 500 500 500 500 500 500 500 500 500 500 500 500 500 500

TS

MS

CIH

Bst

Avg

Bst

Bst

Avg

Bst

Avg

Bst

Avg

10.46 10.94 10.74 10.95 10.58 10.76 10.35 10.85 9.97 10.58 10.63 10.62 10.93 11.02 11 10.45 11.04 10.56 10.46 10.58 38 37 38 39 38 38 38 38 37 38 38 38 39 38 39 38 38 38 38 37 36.49 38.72 38.34 38.6 38.18 38 37.34 37.91 38.68 38.03 38.07 38.58 38.77 38.85 38.31 39.19 38.5 37.15 38.91 38.37

11.44 11.52 11.68 11.56 11.74 11.50 11.65 11.47 11.25 11.56 11.48 11.56 11.57 11.57 11.71 11.45 11.52 11.62 11.43 11.39 39.45 39.25 40.05 39.65 39.45 39.95 40.30 39.50 39.15 39.50 39.50 39.60 40.35 39.50 40.35 39.40 39.45 39.50 39.45 39.45 40.31 40.18 39.94 40.65 39.62 39.64 39.79 39.30 40.06 40.00 40.07 40.26 40.21 40.38 40.22 40.53 40.32 39.70 40.82 39.89

10.53 10.58 11.19 11.08 11.04 10.87 10.6 10.2 11.22 11.1 11.12 11.37 10.31 10.7 10.9 11.29 10.33 11.03 11.25 10.54 38 37 38 38 38 39 39 38 39 38 39 38 38 38 39 37 38 39 38 38 33.37 34.35 33.23 34.28 35.02 35.55 35.41 38.25 33.23 34.32 36.48 33.98 35.84 33.2 35.89 34.4 38.28 35.37 36.46 36.28

11.22 11.42 10.88 11.23 10.98 10.08 11.11 11.36 10.67 10.61 10.57 10.63 10.92 9.95 10.4 10.4 10.97 10.83 11.27 10.96 41 40 41 42 41 40 40 41 41 38 41 40 42 41 41 41 41 41 41 41 41.29 42.8 41.88 41.22 42.28 41.94 41.42 40.43 41.08 41.66 42.93 42.76 42.58 41.66 41.98 41.72 40.67 42.58 41.18 41.21

12.37 12.08 11.86 12.10 12.03 11.89 11.97 12.21 11.83 11.88 11.87 11.85 11.85 11.94 12.13 11.97 12.05 11.92 11.99 12.11 43.30 42.20 43.45 43.15 42.55 42.15 42.20 42.50 42.40 42.30 42.65 42.45 43.10 42.50 42.10 42.60 42.65 42.50 42.35 42.15 44.82 44.51 44.56 43.95 44.00 44.10 43.99 43.49 44.47 44.22 44.14 44.22 44.06 43.96 44.47 44.35 43.82 43.65 44.93 44.78

12.02 11.94 12.22 12.02 11.78 12.16 11.91 12.74 11.98 11.82 11.93 12.32 12.7 11.64 12.09 12.19 11.36 11.49 12.44 12.31 49 51 50 49 50 48 51 47 49 51 51 50 51 49 50 50 50 52 50 50 55.26 56.03 53.44 53.23 54.84 54.66 54.87 55.09 53.82 54.18 56.78 56.35 57.07 54.19 57.38 54.45 52.11 53.58 54.06 55.27

13.09 12.98 13.15 13.06 13.06 13.25 13.35 13.32 13.00 13.21 13.14 13.08 13.43 13.15 13.01 13.16 12.90 12.39 13.17 13.02 53.80 54.15 53.70 54.05 54.80 54.00 55.15 56.05 53.05 54.85 54.25 54.15 53.90 55.20 55.75 53.70 54.90 55.70 53.70 55.25 58.33 60.19 57.72 58.33 57.58 58.01 57.64 57.95 57.55 57.22 58.66 58.64 59.48 58.04 59.27 58.78 57.29 56.36 58.32 57.85

11.42 11.38 11.93 12.08 12.47 11.11 11.75 11.71 12.1 12.01 12.21 12.18 10.93 11.6 11.4 11.97 12.23 11.64 11.73 10.68 48 49 50 50 49 47 45 47 47 45 44 46 45 50 47 51 47 47 48 48 49.15 50.69 47.64 46.85 47.19 48.38 47.15 46.93 47.59 46.29 48.74 49.09 47.88 49.1 49.28 48.1 48.75 44.16 45.83 48.21

13.05 13.05 12.78 12.92 13.10 12.70 12.85 12.81 12.77 12.74 12.90 13.06 12.68 13.01 12.83 12.79 12.72 12.70 12.96 12.58 50.40 50.85 52.70 53.10 52.85 50.10 49.40 50.40 50.30 50.85 49.40 49.10 49.35 52.60 50.35 52.60 49.35 50.90 50.55 50.45 52.40 52.86 50.03 50.96 49.98 50.90 51.31 49.71 51.54 51.44 52.84 52.00 52.58 51.87 52.39 50.82 51.96 50.33 50.59 51.73

algorithm, we use the randomization initialization strategy because it is much faster and has the nonnegligible advantage of being simpler as pointed out in Aringhieri et al. (2014). In order to further verify whether our proposed algorithms are significantly better than the previous state-of-the-art algorithms, we conduct Wilcoxon statistical tests (http://www.vassarstats.net/ textbook/ch12a.html) in terms of the measure devbst % and devavg % for each set of benchmarks. Table 11 shows the p-values obtained from conducting these tests between TS and each compared al-

EPR

VNS

gorithm in terms of the measure devbst %. Table 11 discloses that TS performs significantly better (p-value < 0.05) than each of the algorithms CIH, EPR and VNS for each set of benchmarks. Note that the dominating algorithm is determined by further comparing the average deviations from the best known results between the two algorithms. When comparing our proposed algorithms TS and MS, we find that TS performs significantly better than MS for the MDG-b and MDG-c benchmarks but is significantly outperformed by MS for the MDG-a and GKD-b benchmarks. For the SOM-b,

836

Y. Wang et al. / European Journal of Operational Research 258 (2017) 829–843

Table 3 Computational comparisons on the MDG-b benchmark instances (time in seconds). Instances

MDG-b_1_n500_m50 MDG-b_2_n500_m50 MDG-b_3_n500_m50 MDG-b_4_n500_m50 MDG-b_5_n500_m50 MDG-b_6_n500_m50 MDG-b_7_n500_m50 MDG-b_8_n500_m50 MDG-b_9_n500_m50 MDG-b_10_n500_m50 MDG-b_11_n500_m50 MDG-b_12_n500_m50 MDG-b_13_n500_m50 MDG-b_14_n500_m50 MDG-b_15_n500_m50 MDG-b_16_n500_m50 MDG-b_17_n500_m50 MDG-b_18_n500_m50 MDG-b_19_n500_m50 MDG-b_20_n500_m50 MDG-b_21_n20 0 0_m20 0 MDG-b_22_n20 0 0_m20 0 MDG-b_23_n20 0 0_m20 0 MDG-b_24_n20 0 0_m20 0 MDG-b_25_n20 0 0_m20 0 MDG-b_26_n20 0 0_m20 0 MDG-b_27_n20 0 0_m20 0 MDG-b_28_n20 0 0_m20 0 MDG-b_29_n20 0 0_m20 0 MDG-b_30_n20 0 0_m20 0 MDG-b_31_n20 0 0_m20 0 MDG-b_32_n20 0 0_m20 0 MDG-b_33_n20 0 0_m20 0 MDG-b_34_n20 0 0_m20 0 MDG-b_35_n20 0 0_m20 0 MDG-b_36_n20 0 0_m20 0 MDG-b_37_n20 0 0_m20 0 MDG-b_38_n20 0 0_m20 0 MDG-b_39_n20 0 0_m20 0 MDG-b_40_n20 0 0_m20 0

Time

65 64 66 65 65 65 65 65 65 64 64 65 64 66 66 65 65 65 64 64 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0 20 0 0

TS

MS

CIH

Bst

Avg

Bst

Bst

Avg

EPR Bst

Avg

VNS Bst

Avg

1055.33 1088.56 1100.60 1086.20 1044.08 1062.90 1063.67 1093.98 1075.16 1069.54 1031.02 1063.76 1026.86 1018.69 1022.19 1057.20 1078.89 1082.23 1071.88 1022.66 3421.21 3420.91 3448.59 3305.12 3360.30 3534.09 3361.44 3454.52 3457.26 3373.50 3519.23 3442.42 3444.89 3454.03 3457.00 3442.17 3458.43 3390.50 3476.10 3375.62

1169.67 1151.26 1165.21 1159.17 1156.84 1166.76 1147.94 1168.29 1157.85 1151.70 1165.64 1164.92 1159.24 1175.03 1129.22 1163.72 1145.27 1166.58 1150.50 1149.51 3544.32 3564.41 3550.02 3532.08 3603.87 3630.28 3530.74 3545.25 3553.72 3547.15 3609.88 3566.98 3584.87 3578.48 3580.56 3574.16 3593.93 3572.96 3590.59 3523.60

1132.88 1044.03 1086.91 1053.40 1138.17 1087.69 1093.40 1088.63 1096.05 1102.06 1119.71 1089.34 1098.75 1097.70 1133.37 1120.43 1061.48 1041.14 1066.78 1138.11 3440.18 3389.63 3445.18 3418.32 3447.94 3342.92 3530.08 3473.98 3351.36 3457.83 3520.19 3467.52 3519.10 3498.76 3372.26 3442.49 3352.08 3582.43 3488.40 3351.17

1112.64 1038.08 1126.46 1091.76 1096.98 1070.57 1105.41 1114.92 1069.26 1111.66 1128.63 1138.73 1068.24 1084.51 1106.33 1074.01 1045.20 1151.36 1110.18 1090.95 3592.78 3610.15 3608.12 3599.84 3527.50 3644.37 3693.03 3643.33 3707.34 3678.40 3752.73 3673.65 3706.50 3773.05 3699.91 3715.52 3664.97 3661.20 3672.97 3719.84

1213.27 1214.02 1210.79 1206.73 1204.08 1194.34 1197.72 1216.55 1185.80 1209.16 1200.90 1205.03 1191.48 1190.32 1190.69 1180.83 1196.39 1221.18 1197.71 1193.60 3883.27 3879.67 3808.08 3839.34 3825.67 3880.27 3868.30 3810.18 3870.87 3797.06 3861.12 3797.78 3815.30 3894.40 3883.25 3897.08 3857.85 3803.77 3863.94 3816.35

1185.72 1174.76 1233.75 1140.59 1252.53 1195.36 1237.14 1164.18 1208.88 1240.54 1188.77 1256.29 1164.11 1197.58 1233.09 1193.14 1183.71 1206.22 1203.04 1114.82 4600.85 4333.36 4566.91 4483.36 4429.91 4523.01 4533.26 4389.26 4400.64 4349.86 4313.65 4315.46 4385.88 4632.31 4429.15 4321.26 4549.56 4476.97 4470.91 4426.71

1316.48 1304.88 1317.50 1302.29 1311.55 1309.99 1293.32 1308.19 1314.17 1316.95 1309.01 1311.03 1312.56 1315.08 1302.87 1312.48 1300.38 1307.36 1307.20 1296.75 4778.31 4661.84 4722.15 4707.11 4794.93 4730.99 4701.02 4698.69 4681.13 4764.17 4801.32 4778.58 4697.26 4791.64 4728.08 4653.35 4836.76 4685.33 4698.42 4670.78

1190.37 1239.42 1165.43 1236.85 1164.59 1190.69 1080.78 1223.58 1115.89 1182.99 1254.10 1187.77 1255.16 1214.19 1214.68 1224.80 1085.15 1212.05 1206.56 1223.63 4232.27 4280.79 4196.89 4188.47 4362.02 4145.28 4068.17 4195.74 4039.83 4270.79 4083.42 4240.51 4387.52 4113.29 4119.50 4131.32 4232.38 4295.61 4114.55 4136.50

1297.20 1303.28 1270.39 1292.06 1267.97 1304.83 1260.30 1300.79 1268.64 1279.33 1312.54 1294.75 1312.67 1297.49 1310.45 1305.00 1289.28 1288.89 1290.57 1320.37 4435.83 4520.33 4390.30 4472.02 4557.13 4391.32 4385.32 4477.90 4301.16 4420.86 4415.22 4366.35 4574.32 4529.20 4342.11 4356.16 4381.58 4405.56 4291.46 4391.52

Table 4 Computational comparisons on the MDG-c benchmark instances (time in seconds). Instances

MDG-c_1_n30 0 0_m30 0 MDG-c_2_n30 0 0_m30 0 MDG-c_3_n30 0 0_m30 0 MDG-c_4_n30 0 0_m30 0 MDG-c_5_n30 0 0_m30 0 MDG-c_6_n30 0 0_m40 0 MDG-c_7_n30 0 0_m40 0 MDG-c_8_n30 0 0_m40 0 MDG-c_9_n30 0 0_m40 0 MDG-c_10_n30 0 0_m40 0 MDG-c_11_n30 0 0_m50 0 MDG-c_12_n30 0 0_m50 0 MDG-c_13_n30 0 0_m50 0 MDG-c_14_n30 0 0_m50 0 MDG-c_15_n30 0 0_m50 0 MDG-c_16_n30 0 0_m60 0 MDG-c_17_n30 0 0_m60 0 MDG-c_18_n30 0 0_m60 0 MDG-c_19_n30 0 0_m60 0 MDG-c_20_n30 0 0_m60 0

Time

30 0 0 30 0 0 30 0 0 30 0 0 30 0 0 30 0 0 30 0 0 30 0 0 30 0 0 30 0 0 30 0 0 30 0 0 30 0 0 30 0 0 30 0 0 30 0 0 30 0 0 30 0 0 30 0 0 30 0 0

TS

MS

CIH

Bst

Avg

Bst

Bst

Avg

EPR Bst

Avg

VNS Bst

Avg

4796 4830 4913 4830 4881 6466 6480 6255 6607 6297 7793 7719 7767 7678 7659 9337 8618 9118 9387 9013

5018.60 5020.70 5107.45 4988.05 5118.75 6680.65 6855.30 6518.55 6913.70 6469.70 8064.00 8101.60 8206.10 8114.90 7991.05 9878.05 9529.30 9540.30 9696.40 9618.75

4817 4827 5033 4858 4809 6349 6334 6627 6346 6514 7884 7785 7711 7645 7793 9455 9821 9226 9537 9937

5215 5203 5174 5164 5175 6883 6916 7417 6652 6797 8477 8293 8078 8470 8536 10,066 10,091 10,451 12,313 10,284

5537.60 5393.10 5604.60 5493.75 5431.60 7599.85 7763.75 7894.35 7027.35 7188.35 9086.55 8927.50 9207.35 8859.75 9174.90 11,516.70 11,226.35 11,098.75 13,038.15 11,390.65

6661 6482 6518 6245 6500 8646 8016 8198 8321 9206 10,130 10,081 10,847 10,472 10,489 12,104 13,924 13,322 12,329 12,219

7139.85 7197.70 7294.30 7152.85 6845.75 9513.10 9273.25 9258.80 9116.20 10,022.30 11,486.05 11,965.35 12,232.10 12,394.55 11,945.55 13,846.90 14,663.65 14,411.05 14,364.90 13,966.90

6145 5975 6105 6465 6152 8313 7890 8248 8298 8514 10,236 10,428 10,318 10,327 10,320 12,007 12,083 12,538 12,216 12,231

6393.85 6378.40 6545.25 6723.30 6290.95 8714.50 8690.90 8566.05 8651.60 8912.15 10,896.90 10,735.35 10,692.20 10,885.55 11,032.65 12,406.05 12,978.90 13,077.40 12,870.45 12,707.40

Y. Wang et al. / European Journal of Operational Research 258 (2017) 829–843

837

Table 5 Computational comparisons on the SOM-b benchmark instances (time in seconds). Instances

SOM-b_1_n100_m10 SOM-b_2_n100_m20 SOM-b_3_n100_m30 SOM-b_4_n100_m40 SOM-b_5_n200_m20 SOM-b_6_n200_m40 SOM-b_7_n200_m60 SOM-b_8_n200_m80 SOM-b_9_n300_m30 SOM-b_10_n300_m60 SOM-b_11_n300_m90 SOM-b_12_n300_m120 SOM-b_13_n400_m40 SOM-b_14_n400_m80 SOM-b_15_n400_m120 SOM-b_16_n400_m160 SOM-b_17_n500_m50 SOM-b_18_n50 0_m10 0 SOM-b_19_n500_m150 SOM-b_20_n50 0_m20 0

Time

1 2 2 4 3 8 16 26 9 29 50 82 18 83 103 156 45 152 324 500

TS

MS

CIH

Bst

Avg

Bst

Bst

Avg

Bst

Avg

Bst

Avg

1 4 7 10 4 9 13 18 6 13 19 25 9 16 23 32 10 19 29 38

1.3 5.3 8 11.35 4.75 10.25 14.9 20.15 7.7 14.4 20.6 27.2 9.45 17.85 25.75 33.7 11.5 21.15 31.05 40.2

1 4 7 10 4 9 13 19 6 13 18 24 9 17 24 27 10 20 26 34

0 4 8 10 4 9 16 22 7 14 22 29 9 18 26 36 11 22 33 43

0.95 5.6 8.65 12.45 4.9 10.75 16.65 22.8 7.75 15.4 22.8 30.35 9.75 19.15 28.5 37.6 12.15 23.5 34.2 44.65

1 4 9 12 4 11 17 25 7 17 26 37 10 23 37 50 13 26 44 61

1.95 5.8 9.65 13.7 5.2 12.5 19.15 27.85 8.25 17.85 28.35 40.3 11.3 24.2 38.2 52.45 13.95 30 47.35 65.15

1 4 9 11 4 10 17 23 8 16 23 28 9 20 30 41 12 24 35 49

1.55 5.65 9.4 13.3 5.25 11.7 18.6 25.3 8.2 17.25 25.45 34.55 10.7 22.5 32.15 45.2 13.35 27.55 41.1 54.25

GKD-c and APOM benchmarks, TS and MS present no significant differences. Table 12 shows the p-values obtained from conducting Wilcoxon tests between MS and each compared algorithm in terms of the measure devbst %. As indicated in Table 12, MS is significantly better than EPR for each set of benchmark. In comparison with CIH and VNS, MS presents no significant differences for the GDK-b benchmark and is significantly better for all the other benchmarks. Table 13 shows the p-values obtained from conducting Wilcoxon tests between TS and CIH, EPR and VNS algorithms in terms of the measure devavg %. As indicated in Table 13, TS performs significantly better than these algorithms for each benchmark set. 5. Analysis In this section, we first perform parameter sensitivity analysis for weights of hash functions and time-to-target analysis to evaluate the performance of different algorithms. Then we analyze the role of the hash functions in correctly determining the tabu status of solutions, and also examine the computational gain of using the neighborhood decomposition candidate list compared to using a full neighborhood examination. Finally, we explore the reason why the memetic search algorithm works particularly well for some benchmarks. 5.1. Parameter sensitivity analysis In this section, we analyze the influence on the performance of the tabu search algorithm of the weights iγ1 and iγ2 in defining the hash functions h1 and h2. Proper settings for the parameters γ 1 and γ 2 should satisfy the following two conditions. First, the hash function values of each possible solution need to be no more than the maximum integer 232 − 1 under a 32-bit unsigned integer representation to avoid integer overflow. Second, the hash function values of different solutions should be distributed as widely as possible to reduce hashing collisions. Preliminary experiments find that setting γ 1 , γ 2 > 1.8 will encounter the error of integer overflow for instances with n ≥ 30 0 0 while setting γ 1 , γ 2 ≤ 1 will make adjacent elements i and i + 1 acquire the same weights and thus result in an increased probability of different solutions getting the same hash function value. Hence, the values in the interval (1.0, 1.8] should succeed in satisfying the desired conditions.

EPR

VNS

In order to verify if different settings in the interval (1.0, 1.8] for γ 1 and γ 2 present differences, we fetch 12 different pairs of settings from this interval and conduct the following new experiments. First, we select 13 representative instances, each of which is given 20 independent tabu search runs with each pair of chosen settings while using the same time limit as shown in Section 4.1 to terminate each run. Then, we report average solution values for each instance with each pair of settings in Table 14. Finally, we conduct Friedman statistical tests for the tabulated results. Our experimental findings disclose that different settings of the hash function parameters within the indicated range present no significant difference with the p-value of 0.633. In summary, the parameters γ 1 and γ 2 in hash functions are not sensitive in our algorithm and any values in the interval (1.0, 1.8] are workable for the tested benchmarks.

5.2. Time-to-target analysis To compare the efficiency of the proposed algorithms and other reference algorithms, we apply a time-to-target (TTT) analysis which identifies the empirical probability distribution for the time to achieve a given target value (Aiex, Resende, & Ribeiro, 2007). An alternative tool to compare two or more algorithms is a closed form index (Ribeiro, Rosseti, & Vallejos, 2012). During the TTT experiment, we perform 200 independent runs for each compared algorithm on a single instance. The time needed to attain an objective value at least as good as a given target value for each run is recorded, and these times are then sorted in increasing order so that ti represents the ith smallest time and we associate a probability pi = (i − 0.5 )/200 with each time ti . The points (ti , pi ) are finally plotted. Fig. 3 shows the experimental results of TS, MS, CIH, EPR and VNS algorithms on the instances MDG-a_10_n500_m50 (Top) and MDG-b_20_n500_m50 (Bottom). From Fig. 3, we observe that the probability of reaching a given target value increases for each algorithm as time grows. Notably, the lines of our TS and MS algorithms strictly lie above the lines of the CIH, EPR and VNS algorithms, indicating that TS and MS have a higher probability of reaching the target value within a given time period or, equivalently, requires a shorter computational time to obtain the same probability of successfully striking the target value. In addition, the lines of EPR and VNS algorithms are very close and thus present

838

Y. Wang et al. / European Journal of Operational Research 258 (2017) 829–843

Table 6 Computational comparisons on the GKD-b benchmark instances (time in seconds). Instances

GKD-b_01_n25_m2 GKD-b_02_n25_m2 GKD-b_03_n25_m2 GKD-b_04_n25_m2 GKD-b_05_n25_m2 GKD-b_06_n25_m7 GKD-b_07_n25_m7 GKD-b_08_n25_m7 GKD-b_09_n25_m7 GKD-b_10_n25_m7 GKD-b_11_n50_m5 GKD-b_12_n50_m5 GKD-b_13_n50_m5 GKD-b_14_n50_m5 GKD-b_15_n50_m5 GKD-b_16_n50_m15 GKD-b_17_n50_m15 GKD-b_18_n50_m15 GKD-b_19_n50_m15 GKD-b_20_n50_m15 GKD-b_21_n100_m10 GKD-b_22_n100_m10 GKD-b_23_n100_m10 GKD-b_24_n100_m10 GKD-b_25_n100_m10 GKD-b_26_n100_m30 GKD-b_27_n100_m30 GKD-b_28_n100_m30 GKD-b_29_n100_m30 GKD-b_30_n100_m30 GKD-b_31_n125_m12 GKD-b_32_n125_m12 GKD-b_33_n125_m12 GKD-b_34_n125_m12 GKD-b_35_n125_m12 GKD-b_36_n125_m37 GKD-b_37_n125_m37 GKD-b_38_n125_m37 GKD-b_39_n125_m37 GKD-b_40_n125_m37 GKD-b_41_n150_m15 GKD-b_42_n150_m15 GKD-b_43_n150_m15 GKD-b_44_n150_m15 GKD-b_45_n150_m15 GKD-b_46_n150_m45 GKD-b_47_n150_m45 GKD-b_48_n150_m45 GKD-b_49_n150_m45 GKD-b_50_n150_m45

Time

0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.4 0.4 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.5 4 4 4 4 4 2 1 1 1 1 8 9 8 8 8 3 3 3 3 3 20 18 16 17 17

TS

MS

CIH

Bst

Avg

Bst

Bst

Avg

Bst

Avg

Bst

Avg

0.00 0.00 0.00 0.00 0.00 12.72 14.10 16.76 17.07 23.27 1.93 2.05 2.36 1.66 2.85 42.75 48.11 43.20 46.41 47.72 9.33 8.60 6.91 7.59 8.00 159.19 124.17 106.38 135.85 127.27 11.05 13.05 9.18 11.83 9.20 125.55 194.22 184.27 155.39 172.80 18.17 12.38 13.30 16.58 16.43 207.81 211.77 177.29 197.88 220.76

0.00 0.00 0.00 0.00 0.00 12.72 14.10 16.76 17.07 23.27 1.93 2.05 2.36 1.66 2.85 42.87 50.58 43.20 46.41 50.83 10.26 11.01 9.56 11.83 10.72 162.94 128.55 124.00 136.35 129.95 13.26 15.72 12.40 16.46 11.14 128.23 196.25 195.87 155.58 182.55 21.62 17.41 15.65 19.49 20.75 221.25 212.17 180.58 233.42 230.50

0.00 0.00 0.00 0.00 0.00 12.72 14.10 16.76 17.07 23.27 1.93 2.05 2.36 1.66 2.85 42.75 48.11 43.20 46.85 47.72 9.43 10.32 8.64 8.16 9.35 159.19 124.17 141.62 136.88 127.27 11.24 13.54 9.76 11.83 9.20 125.55 194.22 184.27 155.39 174.34 19.20 17.86 11.83 17.06 19.30 207.81 212.97 210.79 205.39 243.14

0.00 0.00 0.00 0.00 0.00 12.72 14.10 16.76 17.07 23.27 1.93 2.05 2.36 1.66 2.85 42.75 48.11 43.20 46.41 47.72 9.33 9.12 7.59 8.72 10.90 159.19 124.17 106.38 135.85 127.27 11.05 11.81 10.64 13.71 10.12 130.35 194.22 184.27 159.27 161.68 16.48 15.23 14.90 17.56 17.65 207.81 212.97 177.29 197.88 220.76

0.00 0.00 0.00 0.00 0.00 12.72 14.10 16.76 17.07 23.27 1.93 2.05 2.36 1.66 2.85 42.75 48.11 43.20 46.41 47.72 10.95 11.27 10.03 11.72 12.12 166.68 125.98 106.38 137.25 130.77 11.53 15.95 13.82 17.25 12.96 134.19 196.50 196.58 164.41 179.20 21.99 20.41 18.63 22.05 21.48 219.37 215.73 189.10 205.09 232.73

0.00 0.00 0.00 0.00 0.00 12.72 14.10 16.76 17.07 23.27 1.93 2.05 2.36 1.66 2.85 42.75 48.11 43.20 46.41 47.72 9.33 10.32 9.64 8.64 6.91 159.19 124.17 106.38 135.85 127.27 11.05 11.43 13.80 12.92 13.06 135.63 194.22 184.27 165.88 172.09 19.89 18.52 19.76 20.71 22.90 223.17 214.69 200.01 205.39 236.37

0.00 0.00 0.00 0.00 0.00 12.72 14.10 16.76 17.07 23.27 1.93 2.12 2.36 1.66 2.85 42.75 48.11 43.69 46.66 47.72 12.35 13.84 12.34 13.64 12.85 168.31 130.14 107.08 137.66 129.43 11.54 18.56 16.99 18.88 15.39 154.59 201.70 191.32 174.68 180.42 25.87 25.51 24.69 23.57 25.95 233.55 227.81 217.73 227.59 250.63

0.00 0.00 0.00 0.00 0.00 12.72 14.10 16.76 17.07 23.27 1.93 2.05 2.36 1.66 2.85 42.75 48.11 43.20 46.41 47.72 9.33 10.32 7.59 8.14 8.62 159.19 124.17 106.38 135.85 127.27 11.05 13.34 12.31 12.70 11.36 129.19 194.22 184.27 159.27 177.38 19.97 18.68 16.96 19.21 17.72 207.81 212.97 180.00 197.88 234.53

0.00 0.00 0.00 0.00 0.00 12.72 14.10 16.76 17.07 23.27 1.93 2.05 2.36 1.66 2.85 42.75 48.11 43.20 46.41 47.72 11.43 12.58 11.05 11.30 13.05 166.43 127.21 106.54 139.55 131.06 11.40 17.74 15.93 18.21 15.47 147.09 200.60 186.80 169.87 187.83 25.91 25.33 23.02 25.32 24.10 228.81 219.83 198.57 215.30 246.32

similar time performance, while both are outperformed by the CIH algorithm. We take the instance MDG-a_10_n500_m50 (Top) as an example to provide further insight. First, with a time limit of 5 seconds, TS and MS are able to achieve the target value 13 with a probability of 90% and 85%, which are much better than the 40% of CIH and the 5% of EPR and VNS. Second, to make certain that each run of a given algorithm will reach the target value 13, the shortest time limits for TS, MS and CIH are 10 seconds, 13 seconds and 60 seconds, while the time limits for EPR and VNS are unknown because the limit of 63 seconds applied to the instance MDG-a_10_n500_m50 is much shorter than EPR and VNS need to solve MDG-a_10_n500_m50. Third, in order to have a probability of 40% to reach the target value 13, TS and MS consume the shortest time of around 1 second, CIH ranks second at 5 seconds, and EPR and VNS take the longest time of more than 60 seconds. In conclusion, our experiments demonstrate that TS and MS exhibit a high degree of computational efficiency.

EPR

VNS

5.3. Error rates of the hashing vectors In order to show the role of the proposed multiple hash functions for reducing hashing collisions, we compare the “error rates” caused by using two hash functions and a single hash function. An error is counted when an unvisited solution is incorrectly determined as tabu. To be specific, we perform a single tabu search run for 104 iterations on the instance MDG-a_10_n500_m50 and preserve each solution encountered on the search in a population P. Each iteration checks all the neighbor solutions of the current solution to find those containing a 1 in the corresponding positions of each hash vector. For all such neighbor solutions, an error counter is incremented by one if the population P does not contain a duplicate of this solution; otherwise the non-error counter nec is incremented by 1. Finally, the error rate is calculated as ec/(ec + nec ). Fig. 4 depicts how the error rate changes as a function of the number of tabu search iterations for the instance MDG-a_10_n500_m50. From this figure, we observe that the tabu

Y. Wang et al. / European Journal of Operational Research 258 (2017) 829–843

839

Table 7 Computational comparisons on the GKD-c benchmark instances (time in seconds). Instances

Time

GKD-c_1_n500_m50 GKD-c_2_n500_m50 GKD-c_3_n500_m50 GKD-c_4_n500_m50 GKD-c_5_n500_m50 GKD-c_6_n500_m50 GKD-c_7_n500_m50 GKD-c_8_n500_m50 GKD-c_9_n500_m50 GKD-c_10_n500_m50 GKD-c_11_n500_m50 GKD-c_12_n500_m50 GKD-c_13_n500_m50 GKD-c_14_n500_m50 GKD-c_15_n500_m50 GKD-c_16_n500_m50 GKD-c_17_n500_m50 GKD-c_18_n500_m50 GKD-c_19_n500_m50 GKD-c_20_n500_m50

70 69 70 70 68 69 68 68 67 67 70 70 68 69 68 68 69 70 70 71

TS

MS

CIH

Bst

Avg

Bst

Bst

Avg

EPR Bst

Avg

VNS Bst

Avg

7.26 7.06 6.88 6.64 7.50 6.79 7.67 7.65 6.84 7.51 6.42 7.22 6.52 6.38 7.44 7.13 6.61 6.88 7.03 7.20

7.92 8.42 8.00 7.70 8.36 8.21 8.46 9.07 7.96 8.54 7.19 8.17 8.21 7.78 8.14 7.98 7.56 8.67 7.88 8.18

6.39 6.13 6.65 6.78 7.38 7.16 6.84 7.01 8.09 7.37 7.06 6.50 6.62 7.02 6.99 6.51 6.31 6.89 6.84 6.32

8.39 8.10 7.93 7.79 7.66 8.17 8.33 8.25 7.37 8.18 6.71 8.09 8.58 7.61 8.79 8.20 7.36 8.23 8.46 8.18

8.93 9.32 8.97 8.78 9.63 8.82 9.47 9.29 9.28 9.33 8.55 9.54 9.25 8.85 9.38 8.96 8.42 9.51 9.06 9.33

14.75 16.09 17.13 15.21 16.74 16.88 17.20 17.40 16.64 15.95 14.92 16.78 15.88 15.61 15.85 16.12 15.96 16.48 14.31 15.40

17.60 17.79 19.04 18.62 18.55 18.16 18.49 19.57 18.87 18.58 17.34 18.56 17.59 18.73 18.39 18.06 18.21 18.70 17.65 18.49

10.33 10.51 10.74 11.46 10.63 10.82 9.90 10.34 10.98 11.44 10.53 10.59 10.76 11.01 10.40 11.58 10.38 10.67 10.33 11.53

12.68 12.77 13.71 13.50 14.14 12.91 13.93 13.54 13.42 13.87 13.50 13.51 13.55 13.52 13.29 14.15 12.70 13.30 13.74 13.89

Table 8 Computational comparisons on the APOM benchmark instances (time in seconds). Instances

01a050m10 02a050m20 03a100m20 04a100m40 05a150m30 06a150m60 07a200m40 08a200m80 09a250m50 10a250m100 11b050m10 12b050m20 13b100m20 14b100m40 15b150m30 16b150m60 17b200m40 18b200m80 19b250m50 20b250m100 21c050m10 22c050m20 23c100m20 24c100m40 25c150m30 26c150m60 27c200m40 28c200m80 29c250m50 30c250m100 31d050m10 32d050m20 33d100m20 34d100m40 35d150m30 36d150m60 37d200m40 38d200m80 39d250m50 40d250m100

Time

0.1 0.3 2 5 6 23 10 33 24 48 0.2 0.3 2 3 5 21 13 41 28 84 0.1 0.2 2 4 6 21 17 34 23 165 0.3 0.4 2 6 10 23 15 55 45 177

TS

MS

CIH

Bst

Avg

Bst

Bst

Avg

EPR Bst

Avg

Bst

VNS Avg

1.406 14.721 3.645 25.497 6.562 46.987 11.39 63.478 14.595 82.085 1355 5552 4045 9576 6913 13,449 9238 17,502 11,427 21,832 1124 6205 2239 11,098 3550 13,087 4865 19,393 5650 22,050 1049 4564 2374 8979 3234 12,444 4927 18,683 6032 21,001

1.576 14.879 3.89 27.83 8.398 50.814 13.202 67.6 18.96 95.63 1959.65 6162.2 5212.55 11,369.2 7983.75 14,996.5 9954.3 19,261.75 12,194.35 22,712.8 1248.95 6222.4 2872 11,099.4 4211.2 13,087 5700.8 19,494.55 7016.7 22,200.8 1193.6 4575.95 2851.55 8979 4116.15 12,444 5431.8 18,683 6395.05 21,168.35

1.406 14.721 3.645 25.497 6.562 46.987 11.39 63.478 14.557 82.085 1091 5552 4058 9540 6769 13,875 9191 17,782 11,534 21,832 1149 6205 2239 11,098 3709 13,087 5397 19,393 6434 22,050 1049 4564 2513 8979 3910 12,444 4752 18,683 5856 21,001

1.406 14.721 3.645 25.497 6.562 46.987 11.843 63.525 14.951 82.298 1091 5552 3996 10,316 7194 15,945 8197 19,740 11,877 23,935 1124 6205 2435 11,098 3849 13,101 5289 19,537 6569 22,782 1049 4564 2592 8979 3950 13,354 5293 19,175 6616 21,811

1.406 14.973 4.096 28.48 7.389 50.733 13.207 67.307 18.417 93.501 1539.25 6161.95 5471 12,378.75 8213.6 16,873.9 10,705.95 21,293.55 12,897.7 25,552.35 1130.25 6216.6 2943.65 11,838.15 4739.65 13,978.1 5980.15 20,835.8 7682.05 23,611.6 1049 4568.05 2966.15 9355.1 4408.2 14,046.45 5827.75 19,774.1 6989.15 22,473.3

1.406 14.721 3.645 25.497 7.927 46.987 13.898 67.219 19.416 90.191 2019 5805 4565 12,227 7545 17,134 10,134 23,934 12,887 28,992 1124 6205 2686 11,238 4443 13,087 6566 21,280 7917 26,321 1141 4564 2967 8979 4559 13,990 6222 20,793 7980 24,253

1.531 14.721 4.285 25.785 9.005 48.947 15.95 70.83 21.438 94.149 2538.65 6774.9 5746.7 13,050.75 8714.65 18,567.75 11,687.5 25,374.95 14,126.25 30,545.25 1281.35 6268.8 3248.55 11,992.9 4993.85 14,418.55 7011.9 22,722.25 8375.6 27,622.45 1346.25 4577.4 3325.15 9518.75 5007.7 14,663.9 6859.05 22,145.15 8344.05 25,698.65

1.406 14.721 3.918 25.497 7.248 46.987 12.192 64.289 15.165 85.99 1091 5552 5370 11,603 7702 16,843 10,609 22,499 13,161 25,650 1124 6205 2787 11,112 3937 13,087 5334 20,787 7313 24,099 1141 4564 2753 8979 3706 14,195 5762 19,674 7268 22,566

1.471 14.748 4.222 25.607 8.586 49.756 13.838 68.414 18.544 91.4 1595.35 6266.1 5863.6 12,865.9 8610.5 18,025.8 11,332.7 24,095.1 14,024.5 28,025.8 1179.55 6268.8 3249.3 12,211.4 5003.75 14,153.7 6704 21,944.2 8279.9 25,513 1201.05 4564 3252.5 9564.15 4820.6 15,005.5 6633.6 21,603.1 7903.25 24,079

840

Y. Wang et al. / European Journal of Operational Research 258 (2017) 829–843

Number of instances

40 30 20 10

30 20 10

APOM

GKD−c

GKD−b

SOM

MDG−c

MDG−a

APOM

GKD−c

GKD−b

SOM

MDG−a

MDG−c

0 MDG−b

0

40

MDG−b

Number of instances

50

Fig. 2. Number of better solutions attained by TS (Left figure) and MS (Right figure) for each set of benchmarks.

Table 9 Average deviations of the best solution values from the best known results. Benchmark

TS

MS

CIH

EPR

VNS

MDG-a MDG-b MDG-c SOM-b GKD-b GKD-c APOM Average

4.00% 0.86% 0.01% 2.58% 1.19% 5.26% 1.19% 2.16%

1.78% 2.70% 2.24% 1.07% 5.43% 2.39% 1.94% 2.51%

10.24% 6.07% 10.67% 13.25% 4.06% 20.16% 4.55% 9.86%

34.66% 22.38% 36.04% 35.14% 8.61% 140.68% 17.51% 42.15%

25.00% 18.41% 31.68% 22.07% 6.12% 61.14% 10.59% 25.00%

Table 10 Average deviations of the average solution values from the best known results. Benchmark

TS

CIH

EPR

VNS

MDG-a MDG-b MDG-c SOM-b GKD-b GKD-c APOM Average

9.15% 7.56% 5.40% 15.98% 11.57% 21.60% 13.32% 12.08%

17.09% 13.69% 19.46% 25.51% 13.77% 36.82% 16.06% 20.34%

44.22% 31.76% 52.38% 52.59% 32.67% 174.94% 29.12% 59.67%

33.79% 26.51% 38.35% 40.58% 19.65% 101.97% 22.77% 40.52%

search with the single hash function leads to a dramatic growth of the error rates from 0 to 0.84 during the first 20 0 0 iterations. The error rate of 0.84 means that a majority of neighboring solutions are mistakenly tagged as tabu, resulting in a very limited neighborhood exploration. By contrast, the error rate of tabu search with the two hash functions over the first 20 0 0 iterations remains very close to 0 and afterwards slowly grows up to the maximum value of 0.06, which is considered acceptable in view of the superior performance of the tabu search approach. In addition, the same experiments are also conducted for the other instances and the findings are similar to those shown in Fig. 4. Hence, we consider that Fig. 4 stands for a general behavior of error rates of the hash vectors. In conclusion, this experiment discloses the merit of using two hash functions for effectively determining tabu status.

5.4. Computational contribution of the neighborhood decomposition In order to assess the computational gain by using the neighborhood decomposition candidate list over using the full neighborhood, we conduct the following experiment. We select 13 representative instances, extracting two or three from each set of benchmarks, and for each instance, we run the tabu search algorithm for 2 × 104 iterations with each type of neighborhood. The results are reported in Table 15, which shows for each problem instance the best solution value Best, the deviation of Best from the best known result BKR and the computational time in seconds to reach Best. As can be seen from Table 15, the time consumed by the neighborhood decomposition is on average one third as long as consumed by the full neighborhood examination, while achieving slightly better solution quality. It is noteworthy that the neighborhood decomposition strategy does not always execute the best swap move (see Section 2.2), which to some extent enhances search diversification. In conclusion, this experiment clearly confirms the merit of the candidate list strategy based on neighborhood decomposition in accelerating the tabu search procedure. 5.5. Why MS works well on some benchmarks As indicated in the results, the proposed memetic search (MS) achieves significantly better results than the other algorithms for instances with a large m/n value, as exemplified by the DM1A benchmark instances. Considering that an effective crossover operator is generally regarded as the core ingredient to influence the performance of MS, we hypothesize that this phenomenon occurs because a large m/n value implies more elements are shared between two parent solutions, which reinforces the rational underlying the design of our crossover operator. To verify this hypothesis, we conduct an analysis of the structural similarity between local optima in terms of the number of commonly shared elements. Thus, we define the similarity of two locally optimal solutions Ms and Mt as the number of their shared elements divided by the cardinality m, i.e. Sim(Ms , Mt ) = |Ms ∩ Mt |/m. To carry out this experiment, we use the same instances as in the neighborhood analysis and run both the memetic search and the tabu search algorithms to collect high quality solutions. For

Table 11 p-values of TS versus MS, CIH, EPR and VNS in terms of best deviations (Wilcoxon’s test at the 0.05 level).

MS CIH EPR VNS

MDG-a

MDG-b

MDG-c

SOM-b

GKD-b

GKD-c

APOM

0.01686 2.987e-09 1.132e-11 1.69e-11

0.01994 5.821e-10 9.095e-13 9.095e-13

0.001662 1.907e-06 1.907e-06 1.907e-06

0.4409 0.001097 0.0 0 03204 0.0 0 04814

0.0 0 06356 0.02036 0.001979 0.0 0 0143

0.1327 1.907e-06 1.907e-06 1.907e-06

0.2066 0.0 0 08499 1.234e-06 1.421e-05

Y. Wang et al. / European Journal of Operational Research 258 (2017) 829–843

841

Table 12 p-values of MS versus CIH, EPR and VNS in terms of best deviations (Wilcoxon’s test at the 0.05 level).

CIH EPR VNS

MDG-a

MDG-b

MDG-c

SOM-b

GKD-b

GKD-c

APOM

5.399e-08 1.129e-11 1.32e-11

3.035e-05 1.819e-12 1.819e-12

1.907e-06 1.907e-06 1.907e-06

0.001097 0.0 0 03204 0.0 0 04824

0.3243 1.907e-06 0.8202

2.67e-05 1.907e-06 1.907e-06

0.0 0 03987 1.01e-06 5.209e-06

Table 13 p-values of TS versus CIH, EPR and VNS in terms of average deviations (Wilcoxon’s test at the 0.05 level).

CIH EPR VNS

MDG-a

MDG-b

MDG-c

SOM-b

GKD-b

GKD-c

APOM

1.671e-11 1.671e-11 1.671e-11

1.819e-12 1.819e-12 1.819e-12

1.907e-06 1.907e-06 1.907e-06

3.815e-06 3.815e-06 3.815e-06

0.03348 0.0 0 01124 0.002117

1.907e-06 1.907e-06 1.907e-06

0.008623 3.702e-09 2.105e-05

Table 14 Objective values of different parameter settings for γ 1 and γ 2 in the hash vectors. (γ 1 , γ 2 )

(1.1,1.6)

(1.2,1.5)

(1.3,1.6)

(1.3,1.7)

(1.3,1.8)

(1.4,1.8)

20b250m100 40d250m100 MDG-a_10_n500_m50 MDG-a_20_n500_m50 MDG-a_10_n50 0_m20 0 MDG-a_20_n50 0_m20 0 MDG-b_10_n500_m50 MDG-b_20_n500_m50 SOM-b_10_n300_m60 SOM-b_18_n50 0_m10 0 SOM-b_20_n50 0_m20 0 GKD-c_10_n500_m50 GKD-c_20_n500_m50

22,852.75 21,236.35 11.41 11.53 39.8 39.93 1153.19 1133.52 13.9 21.45 40.55 8.59 8.03

22,590.6 21,182.5 11.46 11.66 40.37 40.16 1156.53 1165.15 14.3 21.4 40.35 8.45 7.95

22,915.1 21,123.35 11.62 11.64 40.56 39.88 1157.3 1155.2 14.25 21.6 40.25 8.6 7.94

22,875.35 21,191.4 11.44 11.29 39.63 39.74 1161.96 1130.18 14.15 21.4 40.35 8.17 8.01

22,598.35 21,183.65 11.58 11.64 40.42 40.33 1167.2 1145.16 14.35 21.45 40.25 8.45 7.86

22,426.25 21,207.55 11.51 11.6 39.81 39.87 1137.24 1140.27 14.1 21.45 40.1 8.73 8.14

20b250m100 40d250m100 MDG-a_10_n500_m50 MDG-a_20_n500_m50 MDG-a_10_n50 0_m20 0 MDG-a_20_n50 0_m20 0 MDG-b_10_n500_m50 MDG-b_20_n500_m50 SOM-b_10_n300_m60 SOM-b_18_n50 0_m10 0 SOM-b_20_n50 0_m20 0 GKD-c_10_n500_m50 GKD-c_20_n500_m50

22,715.6 21,145.4 11.56 11.51 39.7 39.85 1147.55 1155.29 14.15 21.4 40.25 8.72 8

(1.5,1.6)

(1.5,1.7) 22,712.85 21,121.7 11.61 11.53 40.35 39.69 1130.78 1163.95 14.25 21.35 40.4 8.65 8.59

each chosen problem instance, we collect 600 local optima of different solution quality and measure the similarity between each pair. Table 16 shows the ratio m/n of the elements in the subset divided by the total elements, the average similarity AvgSim of each pair of collected solutions and the best results of TS and MS extracted from the previous tables. From Table 16, we confirm the hypothesis that a large m/n ratio leads to a large similarity between solutions as measured by the number of their common elements. Moreover, the proposed MS works well for instances that exhibit such a large degree of similarity. For example, the instance SOM-b_20_n50 0_m20 0 with the highest similarities of 0.56 yields a best solution value 34 using MS, which obtains the quality improvement of 10.53% over the best solution value of 38 obtained by TS. In conclusion, MS is shown as a better approach for handling instances of problems with high similarity measures. 6. Conclusion Our proposed tabu search and memetic search algorithms for solving the minimum difference dispersion problem, although departing somewhat from customary designs for TS and MS meth-

(1.6,1.4) 22,678.6 21,169.55 11.56 11.43 40.16 40.13 1130.98 1153.02 14.45 21.45 39.95 8.61 8.51

(1.6,1.8) 22,613.55 21,156.9 11.64 11.69 40.12 40.18 1150.83 1156.39 14.4 21.25 40.35 8.47 8.38

(1.7,1.1) 22,607.8 21,185.35 11.49 11.66 40.25 39.93 1147.15 1176.96 14.3 21.15 40.3 8.35 7.92

(1.7,1.2) 22,659.35 21,154.35 11.67 11.68 39.98 40.06 1157.68 1156.09 14.3 21.1 40.15 8.48 8.15

ods, exhibit remarkable efficacy. The tabu search algorithm integrates a neighborhood decomposition candidate list strategy and a solution-based determination of tabu status implemented by means of two coordinated hash functions. The memetic search algorithm uses a distance-quality based crossover operator to produce a solution with both high quality and good diversity. Extensive experiments on seven benchmark test sets with a total of 250 problem instances demonstrate that the proposed algorithms are collectively capable of finding better solutions for 181 instances, while attaining this performance with a solution effort comparable to that of previous state-of-the-art algorithms. Statistical tests also confirm that our algorithms perform significantly better than the other algorithms tested for most benchmark sets. In addition, we analyze the roles of the key algorithmic components and reveal their contribution to the high performance of our proposed methods. Our study points to an interesting avenue for future research, based on the fact that our memetic algorithm incorporates our TS algorithm and thus can be viewed in two ways – one in which a population-based approach uses tabu search to improve solutions produced by combining solutions, and the other in which a tabu search approach uses a population-based procedure as an

842

Y. Wang et al. / European Journal of Operational Research 258 (2017) 829–843

Table 15 Results of the neighborhood decomposition vs. the full neighborhood. BKR

20b250m100 40d250m100 MDG-a_10_n500_m50 MDG-a_20_n500_m50 MDG-a_30_n50 0_m20 0 MDG-a_40_n50 0_m20 0 MDG-b_10_n500_m50 MDG-b_20_n500_m50 SOM-b_10_n300_m60 SOM-b_18_n50 0_m10 0 SOM-b_20_n50 0_m20 0 GKD-c_10_n500_m50 GKD-c_20_n500_m50

Full Neighborhood

21,832 21,001 10.58 10.38 34.02 34.57 1063.64 1022.66 13 19 34 7.24 6.32

Neighborhood Decomposition

Best

Dev%

Time

Best

Dev%

Time

21,938 21,001 10.58 11.21 37.42 37.63 1063.64 1071.5 13 20 38 7.74 7.38

0.49% 0 0 8.00% 9.99% 8.85% 0 4.78% 0 5.26% 11.76% 6.91% 16.77%

258 107 143 198 1166 1228 121 119 78 620 1504 92 87

21,832 21,001 10.73 10.38 37.37 36.93 1063.64 1071.5 13 20 38 7.32 7.38

0 0 1.42% 0 9.85% 6.83% 0 4.78% 0 5.26% 11.76% 1.10% 16.77%

68 47 39 55 498 457 42 47 45 134 835 26 30

1.0

Instances

1

a single hash function two hash functions

0

Instances:MDG−a_10_n500_m50 Target:13.00

0

5

10

15

20

25

30 35 40 45 Time to target value

50

55

60

65

1 TS MS CIH VNS EPR

0.8

0.2

0.2

0.0

0.4

0.4

0.6

0.6

Error rates

Probability

0.8

0.8

TS MS CIH VNS EPR

0

2000

4000

6000

8000

10000

Probability

Iterations 0.6

Fig. 4. Error rates of the tabu search iterations. Table 16 Results of the structural similarity between local optima.

0.4

0.2 Instances:MDG−b_20_n500_m50 Target:1300 0

0

5

10

15

20

25

30 35 40 45 Time to target value

50

55

60

65

Fig. 3. Empirical probability distribution for the time to achieve a target value.

intensification/diversification strategy invoked at specified intervals (after executing a selected number of TS iterations). Given that the tabu search algorithm performs better than the memetic algorithm on most of the benchmark test problems, our findings motivate an investigation of other intensification/diversification mechanisms as a basis for achieving better performance. A variety of intensification/diversification strategies have been proposed in connection with tabu search (see, e.g., Glover & Laguna, 1997) that could be implemented and explored in this connection. Also, from a population-based standpoint, scatter search and path relinking procedures afford an alternative to the GA-based crossover method used in the memetic algorithm of our present study. Various intensification/diversification methods have likewise been proposed in

Instances

m/n

AvgSim

MS

TS

20b250m100 40d250m100 MDG-a_10_n500_m50 MDG-a_20_n500_m50 MDG-a_30_n50 0_m20 0 MDG-a_40_n50 0_m20 0 MDG-b_01_n500_m50 MDG-b_20_n500_m50 SOM-b_10_n300_m60 SOM-b_18_n50 0_m10 0 SOM-b_20_n50 0_m20 0 GKD-c_10_n500_m50 GKD-c_20_n500_m50

0.4 0.4 0.1 0.1 0.4 0.4 0.1 0.1 0.2 0.2 0.4 0.1 0.1

0.46 0.56 0.10 0.10 0.54 0.56 0.10 0.10 0.21 0.21 0.56 0.12 0.12

21,832 21,001 10.58 11.02 37.37 34.57 1086.22 1126.35 13 19 34 7.49 7.57

21,832 21,001 10.73 10.38 41.10 36.89 996.87 1071.5 13 20 38 7.24 7.38

conjunction with these procedures(see, e.g., Glover, 1998), thereby affording a source of additional mechanisms for enhancing the algorithms of this paper. Acknowledgement We are grateful to the reviewers whose comments have helped to improve our paper. We would like to thank Dr. Abraham

Y. Wang et al. / European Journal of Operational Research 258 (2017) 829–843

Duarte and Jesús Sánchez-Orz Calvo for providing the source code of the EPR algorithm, and Dr. Roberto Aringhieri and the coauthors of Aringhieri et al. (2014) for explaining the computational testing of the CIH algorithm. This work was supported by the National Natural Science Foundation of China (grant nos. 71501157, 71401059, 71320107001, 715310 09, 50 0130 0 0 01), the Natural Science Basic Research Plan in Shaanxi Province of China (grant no. 2016JQ7007) and the China Postdoctoral Science Foundation (grant no. 2015M580873). References Aiex, R. M., Resende, M. G. C., & Ribeiro, C. C. (2007). TTT plots: A perl program to create time-to-target plots. Optimization Letters, 1, 355–366. Aringhieri, R., Cordone, R., & Grosso, A. (2014). Construction and improvement algorithms for dispersion problems. European Journal of Operational Research, 242(1), 1–13. Barbati, M., & Piccolo, C. (2015). Equality measures properties for location problems. Optimization Letters. doi:10.1007/s11590-015-0968-2. Brown, J. R. (1979a). The knaspack sharing problem. Operations Research, 27(2), 341–355. Brown, J. R. (1979b). The sharing problem. Operations Research, 27(2), 324–340. Carlton, W. B., & Barnes, J. W. (1996). A note on hashing functions and tabu search algorithms. European Journal of Operational Research, 95, 237–239. Chandrasekaran, R., & Daughety, A. (1981). Location on tree networks: p-centre and n-dispersion problems. Mathematics of Operations Research, 6(1), 50–70. Della Croce, F., Grosso, A., & Locatelli, M. (2009). A heuristic approach for the max-min diversity problem based on max-clique. Computers & Operations Research, 36(8), 2429–2433. Duarte, A., & Martí, R. (2007). Tabu search and GRASP for the maximum diversity problem. European Journal of Operational Research, 178(1), 71–84. Duarte, A., Sánchez-Oro, J., Resende, M. G. C., Glover, F., & Martí, R. (2014). Greedy randomized search procedure with exterior path relinking for differential dispersion minimization. Information Science, 296(1), 46–60. Erkut, E. (1990). The discrete p-dispersion problem. European Journal of Operational Research, 46, 48–60. Glover, F. (1998). A template for scatter search and path relinking. In J. K. Hao, E. Lutton, E. Ronald, M. Schoenauer, & D. Snyers (Eds.), Lecture notes in computer science: Vol. 1363. Artificial evolution (pp. 13–54). Springer. Glover, F., Kuo, C. C., & Dhir, K. S. (1998). Heuristic algorithms for the maximum diversity problem. Journal of Information & Optimization Sciences, 19(1), 109–132. Glover, F., & Laguna, M. (1997). Tabu search. Boston: Kluwer Academic Publishers. Glover, F., Taillard, E., & de Werra, D. (1993). A user’s guide to tabu search. Annals of Operations Research, 41, 12–37. Glover, F., Ye, T., Punnen, A. P., & Kochenberger, G. (2015). Integrating tabu search and VLSN search to develop enhanced algorithms: A case study using bipartite boolean quadratic programs. European Journal of Operational Research, 241(3), 697–707.

843

Hao, J. K. (2012). Memetic algorithms in discrete optimization. In F. Neri, C. Cotta, & P. Moscato (Eds.), Handbook of memetic algorithms. In Studies in computational intelligence: Vol. 379 (pp. 73–94). Berlin, Heidelberg: Springer. Kincaid, R. K. (1992). Good solutions to discrete noxious location problems via metaheuristics. Annals of Operations Research, 40(2), 65–81. Kortsarz, G., & Peleg, D. (1993). On choosing a dense subgraph. Proceedings of the 34th annual symposium on foundations of computer science (pp. 692–701). Kuby, M. J. (1987). Programming models for facility dispersion: the p-dispersion and maximum dispersion problems. Geographical Analysis, 19(4), 315–329. Martí, R., Gallego, M., Duarte, A., & Pardo, E. G. (2013). Heuristics and metaheuristics for the maximum diversity problem. Journal of Heuristics, 19, 591–615. Martí, R., & Sandoya, F. (2013). GRASP and path relinking for the equitable dispersion problem. Computers & Operations Research, 40(12), 3091–3099. ´ N., Todosijevic, ´ R., & Uroševic, ´ D. (2016). Less is more: Basic variable Mladenovic, neighborhood search for minimum differential dispersion problem. Information Sciences, 326, 160–171. Neri, F., Cotta, C., & Moscato, P. (Eds.) (2011). Handbook of memetic algorithms. Studies in computational intelligence (Vol. 379). Berlin, Heidelberg: Springer. Peng, B., Lü, Z., & Cheng, T. C. E. (2015). A tabu search/path relinking algorithm to solve the job shop scheduling problem. Computers & Operations Research, 53, 154–164. Porumbel, D., Hao, J. K., & Glover, F. (2011). A simple and effective algorithm for the MaxMin diversity problem. Annals of Operations Research, 186, 275–293. Prokopyev, O. A., Kong, N., & Martinez-Torres, D. L. (2009). The equitable dispersion problem. European Journal of Operational Research, 197(1), 59–67. Punnen, A. P., Taghipour, S., Karapetyan, D., & Bhattacharyya, B. (2014). The quadratic balanced optimization problem. Discrete Optimization, 12, 47–60. Resende, M. G. C., Martí, R., Gallego, M., & Duarte, A. (2010). GRASP and path relinking for the max-min diversity problem. Computers & Operations Research, 37, 498–508. Ribeiro, C. C., Rosseti, I., & Vallejos, R. (2012). Exploiting run time distributions to compare sequential and parallel stochastic local search algorithms. Journal of Global Optimization, 54, 405–429. Silver, G. C., Ochi, L. S., & Martins, S. L. (2004). Experimental comparisons of greedy randomized adaptive search procedures for the maximum diversity problem. In Experimental and efficient algorithms, Lecture notes in computer science: Vol. 3059 (pp. 498–512). Berlin, Heidelberg: Springer. Teitz, M. B. (1968). Toward a theory of public facility location. Papers in Regional Science, 21, 35–51. Wang, Y., Hao, J. K., Glover, F., & Lü, Z. (2014). A tabu search based memetic search for the maximum diversity problem. Engineering Applications of Artificial Intelligence, 27, 103–114. Wang, Y., Lü, Z., Hao, J. K., & Glover, F. (2012). Path relinking for unconstrained binary quadratic programming. European Journal of Operational Research, 223(3), 595–604. Woodruff, D., & Zemel, E. (1993). Hashing vectors for tabu search. Annals of Operations Research, 41, 123–137. Wu, Q., & Hao, J. K. (2015). A review on algorithms for maximum clique problems. European Journal of Operational Research, 242(3), 693–709. Xu, H., Lü, Z., Yin, A., Shen, L., & Buscher, U. (2014). A study of hybrid evolutionary algorithms for single-machine scheduling with sequence-dependent setup times. Computers & Operations Research, 50, 47–60.