Hybridization of tabu search with feasible and infeasible local searches for the quadratic multiple knapsack problem

Computers & Operations Research 66 (2016) 199–214

Contents lists available at ScienceDirect

Computers & Operations Research journal homepage: www.elsevier.com/locate/caor

Hybridization of tabu search with feasible and infeasible local searches for the quadratic multiple knapsack problem Jin Qin a, Xianhao Xu a, Qinghua Wu a,n, T.C.E. Cheng b a b

School of Management, Huazhong University of Science and Technology, No. 1037, Luoyu Road, Wuhan, China Department of Logistics and Maritime Studies, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong

art ic l e i nf o

a b s t r a c t

Available online 14 August 2015

The quadratic multiple knapsack problem (QMKP) concerns assigning a set of objects, which interact among themselves through paired profit values, to a set of capacity-constrained knapsacks such that the overall profit of the objects included in the knapsacks is maximized and the total weight of the objects in each knapsack does not exceed the capacity of the knapsack. In this paper we present a highly effective tabu search (TS) approach for QMKP that incorporates a hybridization scheme combining both feasible and infeasible local searches. The feasible local search focuses its search on the most relevant feasible solutions, while the infeasible local search explores a large search space composed of both feasible and infeasible solutions, and employs several tailored move selection rules to keep the infeasible solutions close to the feasible regions located in promising areas. Extensive computational results on a set of 60 benchmark instances in the literature illustrate that the new TS approach compares very favorably with the current state-of-the-art solution methods for QMKP. In particular, the TS approach finds improved best solutions for ten instances. We also analyze the hybridization scheme in the TS approach to ascertain its effect on the performance of the solution method. & 2015 Published by Elsevier Ltd.

Keywords: Tabu search The quadratic multiple knapsack problem Infeasible local search

1. Introduction The quadratic multiple knapsack problem (QMKP) [12] is an extension of two well-known combinatorial optimization problems, namely the multiple knapsack problem (MKP) and the quadratic knapsack problem (QKP). In QMKP, we are given a set of n objects N ¼ f1; 2; …; ng and a set of K knapsacks K ¼ f1; 2; …; Kg. Each object i A f1; 2; …; ng has a positive linear profit pi and a positive weight wi, each pair of objects (i and j) has a positive quadratic profit pij, and each knapsack k A f1; 2; …; Kg has a capacity Ck. The profit achieved by assigning objects to knapsacks is the sum of the linear profit of the selected objects and the quadratic profit obtained by pairs of objects belonging to the same knapsack. QMKP is concerned with allocating each object to at most one knapsack in such a way that the total weight of the objects placed in each knapsack does not exceed its capacity and the overall profit of the objects contained in the knapsacks is maximized. Formally, let xik be the decision variable such that xik ¼ 1 if object i is assigned to knapsack k and

n

Corresponding author. E-mail addresses: [email protected] (J. Qin), [email protected] (X. Xu), [email protected] (Q. Wu), [email protected] (T.C.E. Cheng). http://dx.doi.org/10.1016/j.cor.2015.08.002 0305-0548/& 2015 Published by Elsevier Ltd.

xik ¼ 0 otherwise. Then QMKP can be formulated as the following integer program [12]: Maximize

n X K X

xik pi þ

i¼1k¼1

subject to

n X

n 1 X

n K X X

xik xjk pij

ð1Þ

i ¼ 1 j ¼ iþ1 k ¼ 1

xik wi r C k ;

8 k A f1; 2; …; Kg

ð2Þ

i¼1

and

K X

xik r 1;

8 iA f1; 2; …; ng:

ð3Þ

k¼1

In the above formulation, the objective function (1) is to maximize the total profit of the objects contained in the knapsacks, the constraint (2) stipulates that the total weight of the objects in each knapsack k does not exceed its capacity Ck, and the constraint (3) ensures that no object is allocated to more than one knapsack. Evidently QMKP is NP-hard as it reduces to the NP-hard 0/1 Knapsack Problem [19] when all the quadratic values pij are set to zero and the number of knapsacks sets to one. In addition to its theoretical challenge as an intractable combinatorial problem, QMKP has many important practical applications in different contexts, including the assignment of workmen to different tasks

200

J. Qin et al. / Computers & Operations Research 66 (2016) 199–214

when their ability to cooperate may affect the results [20], the determination of the locations of the earth stations for communication satellites with a budget constraint, and the determination of the sites of railway stations, freight handling terminals, and airports [19]. Hiley and Julstrom [12] first considered QMKP. Since then, considerable efforts have been made by researchers to address the problem. While effective exact methods have been devised to solve the related knapsack problem, MKP, and QKP [3,4,18], to the best of our knowledge, we are unaware of any exact method proposed to solve QMKP in the literature. Instead, heuristics are often employed to find approximate solutions for QMKP in reasonable computing time. The representative heuristics include greedy constructive heuristics, local search methods, and population-based evolutionary algorithms. For instance, Hiley and Julstrom [12] proposed three methods for QMKP, namely a greedy heuristic, a stochastic hill-climber, and a genetic algorithm, in which the crossover operator preserves the assignments of objects to knapsacks that are common to both selected parents. Singh and Baghel [20] presented a steady-state grouping genetic algorithm to address QMKP. The crossover employed by the genetic algorithm iteratively selects one of the two parents, transmits the knapsack with the largest profit value from the selected parent to the child offspring, and updates the remaining knapsacks associated with both the parents by removing all the objects belonging to the selected knapsack from both the parents. In [19], the authors addressed QMKP with another genetic algorithm using a specialized crossover that interchanges the object assignments between two randomly selected parents. In [21], an artificial bee colony algorithm combined with local search (SS-ABC) was introduced for QMKP. More recently, GarcíaMartínez et al. [7] applied the tabu-enhanced iterated greedy method (TIG-QMKP) to tackle QMKP. This method alternates between constructive and tabu-enhanced destructive phases linked by an improvement process based on exchanging objects from different knapsacks. Later, García-Martínez et al. [5] successfully applied the strategic oscillation method (SO-QMKP) to solve QMKP. This approach iteratively applies an oscillation method that permits the critical boundary to be crossed and uses two specific steps (a divergent step and a convergent step) to generate a new promising solution, a descent-based improvement method to refine the new solution, and an acceptance criterion that decides which configuration becomes the new current solution. Most recently, Chen and Hao [2] presented the iterated responsive threshold search (IRTS), a very effective heuristic approach, for QMKP. This approach alternates between a threshold-based exploration phase and a descent-based improvement phase, and triggers a dedicated perturbation when the search is judged to be stagnating. For QMKP, the knapsack's capacity is imposed as a strict constraint. As for strictly constrained problems, a natural issue is how to effectively explore the search space so as to enhance the performance of the algorithm [8,16]. Indeed, accounting for problem constraints in the definition of the search space often restricts the search process too much and can lead to mediocre solutions [8]. In this paper we propose a highly effective tabu search (TS) method that incorporates a hybridization scheme combining both feasible and infeasible local searches to effectively explore the search space of QMKP. The proposed TS method embraces several features that contribute to its effectiveness: (1) a multi-neighborhood local search procedure that intensively explores the feasible search space through a joint use of four complementary neighborhoods. (2) An oscillation-based local search procedure that permits searches through the infeasible search space to bring more freedom to the search. (3) An effective hybridization scheme that follows the general idea presented in

[22], and combines both feasible and infeasible local search methods in the TS algorithm for solving QMKP. In order to evaluate the performance of our proposed TS algorithm, we conduct experimental tests on a set of wellknown QMKP instances. The computational results indicate that our TS algorithm yields highly impressive outcomes for these instances by finding improved best known solutions for ten instances and matching the best known results for 49 instances. Only for one instance does the TS algorithm fail to reach the previous best known result. Furthermore, we analyze the effect of the hybridization of feasible and infeasible local search methods in the TS paradigm and establish its key role in boosting the performance of the proposed TS algorithm. The remainder of this paper is organized as follows: in Section 2 we present in detail our proposed TS algorithm for solving QMKP. In Section 3 we discuss the computational results and comparisons with state-of-the-art approaches on a set of benchmark instances, followed by analyses of the important components of the TS algorithm in Section 4. We conclude the paper and suggest topics for future research in Section 5.

2. A tabu search algorithm for QMKP 2.1. Main framework Algorithm 1. The framework of the tabu search algorithm for QMKP. Require: A QMKP instance Ensure: The best solution Sbest found and its objective value fbest 1. Begin n 2. S’initial_solutionðÞ / Construct a feasible initial solution, n Section 2.3 / 3. Sbest ’S /n Sbest records the best solution found so far n/ 4. f best ’f ðSÞ /n fbest records the best objective value reached so far n/ 5. While stop condition is not met do 6. /nthe feasible local search phase n/ 7. ðS1 ; Slocal  best Þ’feasible  local  searchðS; Ωf Þ /n Section 2.5 n/ 8. if f ðSlocal  best Þ 4 f best then 9. Sbest ’Slocal  best /n Update the recorded best solution n/ 10. f best ’f ðSlocal  best Þ 11. end if 12. /nthe infeasible local search phase n/ ðS; Slocal  best Þ’infeasible  local  searchðS1 ; ΩÞ /nSection 2.6

13. n

/

if f ðSlocal  best Þ 4 f best then Sbest ’Slocal  best /n Update the recorded best solution n/ f best ’f ðSlocal  best Þ end if if S is an infeasible solution then S’repair  solutionðSÞ /n Solution repair procedure, Section 2.7 n/ 20. end if 21. end while 22. return Sbest and fbest 23. End

14. 15. 16. 17. 18. 19.

Algorithm 1 presents the main framework of our proposed TS approach for QMKP that combines both feasible and infeasible local searches. Essentially, it alternates between a feasible local search (FLS) phase and an infeasible local search (ILS) phase. FLS and ILS play different roles in our proposed method. FLS is used to

J. Qin et al. / Computers & Operations Research 66 (2016) 199–214

intensify our TS algorithm by focusing its search on the most relevant feasible solutions, while ILS is used to diversify the search of the TS method by introducing more freedom during the search process [9,15]. Iteratively applying these two complementary search phases, the algorithm is able to explore various zones of the search space without being easily trapped in local optima. Computational experiments show that the combination of these two local search methods constitutes a very effective refinement method for obtaining high quality solutions. Starting with an initial solution S generated by the greedy procedure given in Section 2.3, the algorithm first enters the feasible local search phase. In order to discover the most promising feasible solutions in the feasible search space, FLS uses a multineighborhood local search procedure to make an intensified examination in the feasible search space through a joint use of four complementary neighborhoods. Specifically, in each iteration of FLS, the algorithm jointly examines four neighborhoods induced by four move operators and selects the most favorable neighboring solution that is subject to the tabu conditions to generate the next solution. However, restricting the search process only to feasible solutions can easily make the algorithm blocked, since the capacity constraints are sometimes too tight to allow moving objects beneficially between knapsacks. In such cases, constraint relaxation is an attractive strategy [5,10,11,14,16,26,28], since by visiting some intermediary infeasible solutions, despite violating the capacity constraints, the algorithm can easily switch from one high quality feasible solution to another (see Fig. 1 for example) [16]. Thus, when FLS is judged to be trapped in a deep local optimum, the algorithm switches to the infeasible local search phase, which is based on the oscillation-based local search procedure [9]. The basic idea of the oscillation-based local search is to allow the algorithm to temporarily violate the capacity constraints and let the algorithm explore in a larger search space consisting of both feasible and infeasible solutions. Furthermore, in order to keep the search process close to the boundaries of the feasible region, it uses a specialized move operator that always moves objects from heavy weight knapsacks to light weight knapsacks, and a normalized move gain criterion that takes into

consideration both the move gain value and the weight of the selected object. Finally, if the final solution at the end of the ILS phase is an infeasible one, a simple repair procedure (see Section 2.7) is invoked to transform this solution into a feasible one before the algorithm switches back to the FLS phase. The algorithm then iterates the above two search phases until the timeout limit is reached. 2.2. Preliminary definitions In this section we provide some preliminary definitions that are useful for a precise description of the proposed TS algorithm. Definition 1. For a given candidate solution S, regardless of its feasibility, its quality is evaluated by its objective function, i.e., P P P P f ðSÞ ¼ Kk ¼ 1 ð ni¼ 1 xik pi þ ni ¼ 11 nj¼ i þ 1 xik xjk pij Þ. Definition 2. Given a solution S A Ω, the value contribution ðΔði; kÞÞ of object i relative to a knapsack k A K is defined as the sum of the profit value of object i and its joint values with the objects in P knapsack k, i.e., Δði; kÞ ¼ pi þ j A k pij . Definition 3. Given a solution S A Ω, the value density ðDði; kÞÞ of object i relative to a knapsack k A K is defined as the value contribution Δði; kÞ of i relative to k divided by the weight of i, i.e., Dði; kÞ ¼ Δði; kÞ=wi . Definition 4. . For a given infeasible solution S ¼ fI 1 ; …; I K g, the degree of infeasibility of the solution S (denoted by DI(S)) is measured by the total overloaded parts of all the knapsacks associated with the solution. 2.3. Generation of initial solutions Our algorithm starts with an initial feasible solution S A Ωf and then uses the tabu search procedure (Sections 2.5 and 2.6) to improve S by maximizing its overall profit. Specifically, we use a greedy construction method similar to that used in [2,5,7,12,19– 21] to generate the initial solution of our tabu search procedure.

2

C1

C2

C3

2

2

2

S

25

1

2

4

2

12

5

4

12

11

5 7

12

13

k1

k2

k3

feasible

k1

2

2

2

S

2

2

2

12

4 12

k1

5

11 7

k2

k1

feasible

13

k3

2

S

2

2

11 7

12 2

k3

12

5

2

k3

12

11

4 12

k1

7

k2

infeasible

13

k2

5

feasible

13

k2

4

S'

11

12

7

S1

201

feasible

13

k3

Fig. 1. Enlarging the solution space through capacity violation. The number inside each object provides the weight of the corresponding object. (a) Infeasible solutions are forbidden. Solution S2 (a more promising solution) cannot be accessed from the current solution S1. (b) Infeasible solutions are allowed. Solution S2 can be easily accessed from S1.

202

J. Qin et al. / Computers & Operations Research 66 (2016) 199–214

Initially, all the knapsacks are empty and the greedy construction method fills the K knapsacks one by one. To fill a knapsack k A f1; 2; …; Kg, an unassigned object is first randomly selected to put into the knapsack. Then in each construction step, we identify an unassigned object that fits the current knapsack (i.e., P wi þ j A Ik wj r C k ) and has the highest relative value density D(i, k) (see Definition 3 in Section 2.2) with respect to the objects already in the knapsack, and assign this identified object to the knapsack. This process is repeated until it becomes impossible to add any more objects to the knapsack. After knapsack k is constructed, the remaining knapsacks are filled using the same method to complete the construction of an initial solution. Now we consider the time complexity of the initial solution procedure. In each construction step of the initial solution procedure, one unassigned object with the highest D(i,k) value is allocated and then the density value D(i,k) is updated. Given that D(i,k) can be updated in O(n) time [2] and the initial solution procedure can be executed for at most n steps, the initial solution procedure requires Oðn2 Þ time in the worst case [2].

2.4. Basic move operators and neighborhoods In a local search, a neighborhood is typically defined by a move operator that transforms a given solution S into a neighboring solution S0 . Our TS algorithm jointly employs four basic move operators to generate its neighborhoods during the two search phases of the algorithm. Other heuristic approaches have used these basic move operators [2,5,7], which we briefly describe in the following.

2.4.1. Extraction move operator Given a solution S, regardless of its feasibility, the extraction move operator removes an already assigned object i from its knapsack k. Let extr k-0 denote such a move (we use a dummy i knapsack 0 to include all the unassigned objects) and S0 be the neighboring solution obtained by applying extr k-0 to S. Then the i change in the objective function value between the two solutions δextrk-0 ¼ f ðS0 Þ  f ðSÞ, also known as the move gain value, can be i efficiently calculated as follows: δextrk-0 ¼  Δði; kÞ:

ð4Þ

i

In our implementation, we employ an n  K matrix Λ with elements Λik ¼ Δði; kÞ recording the value contributions. After an extraction move operation is performed, we only need to update a subset of the value contributions affected by the move as follows: Δðj; kÞ ¼ Δðj; kÞ  pij ;

8 j ai:

ð5Þ

Evidently, the update of the value contributions after an object extraction operation can be implemented in O(n) time.

2.4.2. Insertion move operator Given a solution S A Ω, the insertion move operator (denoted by ) inserts an unassigned object i into a knapsack k. The inse0-k i move gain value of an insertion move inse0-k can be efficiently i computed by δinsek-0 ¼ Δði; kÞ:

ð6Þ

i

0

δrealk-k0 ¼ Δði; k Þ  Δði; kÞ:

ð8Þ

i

k-k0 reali

After a re-allocation move is performed, the value contributions can be efficiently updated using the following formula: ( 0 0 Δðj; k Þ ¼ Δðj; k Þ þ pij ; 8 ja i ð9Þ Δðj; kÞ ¼ Δðj; kÞ  pij ;

8 ja i:

ð10Þ

2.4.4. Exchange move operator This move operator swaps two assigned objects i and j in two different knapsacks ki and kj (case 1) or exchanges an object i in knapsack ki with an unassigned object j (case 2). Let exchi2j denote such a move. Then the move gain value of δexchi2j can be easily computed as follows: ( Δði; kj Þ  Δði; ki Þ þ Δðj; ki Þ  Δðj; kj Þ  2pij caseð1Þ ð11Þ δexchi2j ¼ caseð2Þ: ð12Þ Δðj; ki Þ  Δði; ki Þ  pij It is noted that an exchange move exchi2j can be decomposed into k -k a re-allocation move reali i j followed by another re-allocation k -k move realj j i (case 1), or an extraction move extr iki -0 followed by i an insertion move inse0-k (case 2). So, in order to update the j contribution matrix after an exchange move exchi2j , we could first k -k use (11) and (12) to update the change induced by reali i j , then k -k apply (11) and (12) to update the change induced by realj j i (case 1), or first use (7) to update the change induced by extr iki -0 , then i apply (9) to update the change induced by inse0-k (case 2). It is j clear that updating the contribution matrix after an exchange move can also be performed in O(n) time. A neighboring solution S can be generated by applying one of the above four move operators. The size of the neighborhood induced by the insertion, extraction, or re-allocation move operator is bounded by O(n), while that of the neighborhood induced by the exchange move operator is bounded by Oðn2 Þ. The calculation of the objective function value for each neighboring solution can be implemented in Oð1Þ time with the contribution matrix. The updating of the contribution matrix after application of each of the four move operators takes O(n) time. 2.5. The feasible local search phase Our feasible local search procedure is based on the general tabu search framework and combines several neighborhoods during its search process, so it can be viewed as a multi-neighborhood tabu search for QMKP [24]. In our tabu search, FLS is used to intensively examine the feasible search regions so as to find the most promising solutions in these regions. FLS jointly explores four neighborhoods N1, N2, N3, and N4, and restricts its search only within the feasible search space Ωf. The four neighborhoods explored by FLS are generated by applying the insertion, extraction, re-allocation, and exchange move operators introduced in the last section, respectively, to the incumbent solution under the capacity constraints. Algorithm 2. The feasible local search phase.

The value contributions after an insertion move insek-0 can be i efficiently updated using the following formula and can also be implemented in O(n) time: Δðj; kÞ ¼ Δðj; kÞ þ pij ;

2.4.3. Re-allocation move operator 0 k-k The re-allocation move (denoted by reali ) re-allocates an 0 object i from its current knapsack k to another knapsack k . Similarly, the move gain value of re-allocating i from knapsack k 0 to k can be easily computed by the following formula:

8 j ai:

ð7Þ

Require: An initial solution S, integer Ncons (the search depth) Ensure: The final resulting solution S, the best solution found during this phase Slocal  best 1: Begin

J. Qin et al. / Computers & Operations Research 66 (2016) 199–214

2: n’0 /n n is the number of consecutive iterations during which f(S) is not improvedn/ 3: Slocal  best ’S /nSlocal  best records the best solution found during the phase*/ 4: While ðn o N cons Þ do 5: Initialize tabu  list 6: Construct neighborhoods N1; N2; N3; N4 A Ωf from S 7:

Choose an overall best allowed neighbor S1 A N1 [ N2 [ N3 [ N4 according to the maximum gain criterion

8: S’S1 =n Perform the best move n/ 9: Update tabu  list 10: if f ðSÞ 4 f ðSlocal  best Þ then 11: Slocal  best ’S 12: n’0 /n Reset counter for the number of consecutive non-improving iterations n/ 13: else 14: n’n þ1 15: end if 16: end while 17: return S and Slocal  best 18: End

When several neighborhoods are available, the issue of finding an effective way to exploit the neighborhoods arises. In the literature, the above four neighborhoods are typically combined in the token-ring way [2]. To strengthen the search ability of our FLS in the feasible regions, in this work we propose to explore the four neighborhoods in a combined manner by examining the union of the four neighborhoods N ¼ N1 [ N2 [ N 3 [ N 4 . Specifically, in each iteration of the multi-neighborhood exploration, the algorithm examines all the insertion, extraction, re-allocation, and exchange moves, and selects among them the best admissible (non-tabu or globally improving) solution that yields the largest move gain value to generate the next solution (see Lines 6–8, Algorithm 2). Obviously, such a combined neighborhood enables our algorithm to make a more thorough and intensified examination of the feasible search space, so increasing the chances of finding feasible solutions of better quality. To prevent the search from entering short-term cycling, each time an object i is removed from a knapsack k (by an extraction, re-allocation or exchange move), it is forbidden to assign i back to its original knapsack k for the next tt iterations (called tabu tenure). The tabu tenure tt is tuned according to the number of objects n as follows: tt ¼ n  α;

ð13Þ

203

constraints and this confines the search process to a localized portion of the search space [10,11]. In order to introduce more diversity into the search process, thereby enhancing the capacity of the algorithm to search unexplored areas in the search space, the algorithm switches to the infeasible local search phase when the feasible local search is judged stagnating. The basic idea of ILS is to relax the capacity constraints and allow the algorithm to explore in a larger search space containing both feasible and infeasible solutions. By relaxing the capacity constraints, more beneficial object moves between knapsacks become available, which enables a much easier transition from one promising feasible region to another. The proposed ILS is also based on the tabu search framework and can be viewed as the oscillation-based tabu search procedure. Essentially, it uses two types of move operators to explore the search space Ω consisting of both feasible and infeasible solutions. For the first type of move operators, the capacity constraints are strictly enforced and these move operators are applied only to feasible solutions. For the second type of move operators, we relax the capacity constraints and use a simple re-allocation move operator to explore both feasible and infeasible solutions. This move operator not only helps optimize the objective function value (overall profit of the objects), but also takes care of the infeasibility degree of the solution. Typically, the re-allocation move operator always moves objects from heavy weight knapsacks to light weight knapsacks (we assume that all the knapsacks have the same capacity). In this way, the weight of each knapsack can be maintained close to its capacity, and more solutions will be generated inside the feasible region and also near its boundaries, so yielding better chances of finding high quality feasible solutions. Given the incumbent solution S ¼ fI 1 ; I 2 ; …; I k g, these two types of move operators are used in the oscillation-based search under specific conditions as explained below.

 If S is a feasible solution: ILS first examines all the insertion, re-



allocation, and exchange moves under the capacity constraints (the first type of move operators), and identifies among them the best move with the largest move gain value. If this move produces a neighboring solution S0 better than S in terms of solution quality, it accepts this move and replaces the incumbent solution S by S0 . Otherwise, it applies the second type of move operators to S. If S is an infeasible solution or a local optimal feasible solution: ILS applies the second type of move operators to S, which reallocates objects to knapsacks with smaller weights without considering the capacity constraints. This move operator can be described more precisely as follows: first, choose randomly a knapsack k, if W k o C k . Then select an object i with the largest normalized move gain value δ realk0 -k whose current knapsack is i 0 0 k such that W k0 Z W k . Move object i from k to k. Otherwise, an object i ðiA I k Þ with the largest normalized move gain value δ k-k″ and satisfying W 00k o W k is selected and moved from reali ″ knapsack k to another knapsack k .

where α is a parameter that takes a value in the range ½0; 1. To accompany this tabu rule, we apply a simple aspiration criterion: a tabu move leading to a solution better than the best solution found so far is always accepted. The multi-neighborhood exploration stops if no improved solution is found for Ncons consecutive iterations (Ncons is called the search depth). In this case, the search is judged to be trapped in a deep local optimum and more diversification is needed to bring the search to a new search region. To escape from such a situation, the algorithm switches to the infeasible local search phase to explore in a larger search space.

of the second type, we use the For a given move reali normalized move gain value δ realk-k0 to evaluate its quality, which i is defined as follows:

2.6. The infeasible local search phase

δ realk-k0 ¼

The feasible local search phase ensures an aggressive exploration of the most relevant feasible search space. However, it may be trapped in deep local optima. Indeed, when restricting the search process only to feasible solutions, too many beneficial object moves between knapsacks are prohibited by the capacity

where β is a parameter taking a value in ð0; 1Þ. The purpose of using the normalized move gain value to evaluate a move of the second type is to strike a balance between the quality and the degree of infeasibility of the solution by taking into account both the move gain value and the weight of the moved object.

k-k0

0

i

Δði; k Þ  Δði; kÞ ; wβi

ð14Þ

204

J. Qin et al. / Computers & Operations Research 66 (2016) 199–214

Algorithm 3. The infeasible local search phase. Require: A feasible solution S obtained in the feasible local search phase, integer M (the allowed maximum number of iterations M for the infeasible local search phase) Ensure: The final resulting solution S, the best feasible solution found during this phase Slocal  best 1: Begin 2: m’0 3: Slocal  best ’S /nSlocal  best records the best feasible solution found during this phase*/ 4: Initialize tabu  list 5: While m oM do 6: if S is feasible and is not a local optima induced by the first type of move operators then 7: /n apply the first type of move operators n/ 8: Construct neighborhoods N2; N3; N4 A Ωf from S using the first type of move operators 9:

Choose an overall best allowed neighbor S1 A N2 [ N3 [ N4 according to the maximum gain criterion S’S1 /n Perform the best-improving move n/ Update tabu  list if f ðSÞ 4 f ðSlocal  best Þ then Slocal  best ’S end if else /n apply the second type of move operatorsn/ Randomly choose a knapsack k of solution S if W k o C k then 0 Select an object i A I k0 , k A fkb A Kj W kb Z W k g such

10: 11: 12: 13: 14: 15: 16: 17: 18: 19:

that: maxi A Ik0 fδ realk0 -k g i

20: Move object i to knapsack k to generate a neighboring solution S1 S’S1 else Select an object i A I k and a knapsack

21: 22: 23: ″

k A fks A Kj W ks oW k g such that: maxi A Ik fδ

k-k″

reali

g

″

24: Move object i to knapsack k to generate a neighboring solution S1 25: S’S1 26: end if 27: Update tabu  list 28: if S is a feasible solution and f ðSÞ 4 f ðSlocal  best Þ then 29: Slocal  best ’S 30: end if 31: end if 32: m’m þ1 33: end while 34: return S and Slocal  best 35: End

Algorithm 3 summarizes the general scheme of ILS. The procedure starts with a feasible solution generated in the FLS phase. Then in each iteration of ILS, it determines and applies a move operator according to the current search context. If the incumbent solution S is a feasible solution and S is not a local optimum, it applies the first type of move operators to S (see Lines 8–15 of Algorithm 3); otherwise, it applies the second type of move operators to S (Lines 18–36). During the search process, if a feasible solution better than the current best feasible solution Sbest in terms of solution quality is generated, Sbest is updated and replaced by S accordingly. Furthermore, to avoid short-term

cycling, each time an object i is moved from its original knapsack k (by the first type or the second type of move operators) to another knapsack, it is forbidden to bring i back to knapsack k for the next tt iterations. The tabu tenure tt is tuned using the same method as described in Section 2.5. The procedure stops when an allowed maximum number of iterations M are reached. 2.7. Solution repair procedure We use the solution repair procedure to transform an infeasible solution into a feasible one so as to bring the search process back to the feasible region. Specifically, in each step of the solution repair procedure, it examines all the extraction, re-allocation, and exchange moves without considering the capacity constraints, then among all the moves improving the infeasibility measure DI (S), it selects a best move with the largest move gain value to generate the next solution. This process is repeated until the incumbent solution becomes feasible.

3. Computational experiments We conduct all the experiments on a computer with an Intel Xeon E5440 processor (2.83 GHz CPU and 2GB RAM). We code the proposed TS algorithm in C þ þ and compiled the source code by the gcc compiler without using any optimization option, although using the compiler's optimization options, e.g., -O3 flag, slightly improves the speed of the algorithm. 3.1. Benchmark instances To evaluate the performance of the proposed TS algorithm, we carry out extensive experiments on a set of 60 well-known QMKP benchmark instances that are frequently used to evaluate algorithms for QMKP. These instances were originally proposed by Billionnet and Soutif for the purpose of testing solution algorithms for QKP and are available at: http://cedric.cnam.fr/soutif/QKP/QKP. html [1]. These instances are characterized by three significant features: the number of available objects n A f100; 200g, the number of knapsacks K A f3; 5; 10g, and the density d A f0:25; 0:75g (i.e., the proportion of non-zero pij). For each instance, we set the capacities of the knapsacks to 80% of the sum of the object weights divided by the number of knapsacks. 3.2. Parameter settings In this section we examine the effects of the parameters on the performance of the algorithm in order to obtain a fine-tuned tabu search algorithm. Table 1 shows the settings of the parameters of our TS algorithm. We determine the parameter values by performing a preliminary experiment on the whole set of the 60 QMKP instances. Since these parameters are independent of one another, we study them separately. For each parameter, we tested different Table 1 Settings of the parameters. Parameter Section Ncons M α β

Section 2.5 Section 2.6 Section 2.5 Section 2.6

Description

Value

The search depth of each FLS phase

1000

Maximum number of iterations of each ILS phase Tabu tenure management factor Weight factor

600 0.1 0.7

J. Qin et al. / Computers & Operations Research 66 (2016) 199–214

Table 2 Results of the Friedman test on parameter settings. Parameters

Ncons

M

α

β

p Value Diff ?

0.223 No

0.287 No

3.73e 11 Yes

8.6 e  05 Yes

Table 3 The post hoc test for solution sets obtained by varying the value of α. Values less than 0.05 are presented in bold. α¼

0.05

0.10

0.15

α ¼ 0.10 α ¼ 0.15 α ¼ 0.20 α ¼ 0.25 α ¼ 0.30

0.0125 0.2526 0.8183 0.1741 0.0095

0.9518 0.1702 0.0163 1.77e 04

0.2906 0.001 2.36e  05

0.20

0.0010 6.71e  06

0.25

205

with small values of α A ½0:05; 0:2 are statistically better than those obtained with larger values of α A ½0:25; 0:30. In view of these observations, we choose 0.1 as the default value of α, which yields the best average objective value among all the tested α values. (2) For parameter β, the results of the post hoc test in Table 5 disclose that when β varies in the range [0.1, 0.8], the performance differences among the algorithm with different settings of β are not significant, implying that the algorithm is not particularly sensitive to β with values from 0.1 to 0.8. Furthermore, the statistical results in Table 6 reveal that the algorithm with β ¼ 0:7 performs statistically better than with other β values. Following this observation, we choose 0.7 as the default value of β. In summary, we adopt Ncons ¼ 1000, M ¼600, α ¼ 0:1, and β ¼ 0:7 as the default settings of our TS algorithm.

9.08e  08

Table 4 Statistical results obtained when varying α in the given range. We mark the best result in each row in bold. α¼

0.05

0.10

0.15

0.20

0.25

0.30

Tot. #Best

5,022,045 57

5,022,618 59

5,022,498 57

5,022,294 55

5,021,623 54

5,020,883 52

values for the parameter in a reasonable range while fixing the rest of the parameters to their default values given in Table 1. For each instance and each parameter setting, we ran our TS procedure 40 times, with each run being limited to 15 s for small instances (n¼ 100) and 90 s for large instances (n¼ 200). Specifically, we tested values for Ncons in the set f100; 500; 1000; 1500; 2000g, M in the set f300; 450; 600; 750; 900g, α in the set f0:05; 0:10; 0:15; 0:20; 0:25; 0:30g, and β in the set f0:1; 0:2; 0:3; 0:4; 0:5; 0:6; 0:7; 0:8; 0:9g. Then we perform the Friedman test to see whether there is any difference in the performance of our TS algorithm in terms of its average solution value when varying the value of a single parameter as mentioned above. Table 2 summarizes the results of the Friedman test for the four parameters with 0.05 as the significance level. From Table 2, we observe that there are no statistical differences in performance (with p value 4 0.05) when varying the values of Ncons and M, implying that our algorithm is not particular sensitive to these two parameters. On the other hand, the Friedman test indicates that the TS algorithm is sensitive to α and β since statistical differences are detected by the test with p values of 3.73e  11 and 8.6e  5, respectively. Hence, we perform the post hoc test to examine the statistical difference between each pair of settings of the parameters α and β. We show the results of these analyses in Tables 3 and 5, where each table entry indicates the resulting p value for two sets of the average results obtained with two different values of the corresponding parameter. In addition, we show in Tables 4 and 6 the sum of the average objective values and the number of best solutions obtained with different values of α and β to support our choices of the adopted parameter settings. (1) For parameter α, we observe from Table 4 that there are significant differences between sets of the average solutions obtained with small values of α A ½0:05; 0:2 and those obtained with larger values of α A ½0:25; 0:30. Furthermore, the statistical results in Table 4 show that the set of results generated

3.3. Comparison of results with a long time limit We report the computational results of our TS algorithm under two different stopping criteria. One is the relaxed time condition, which was proposed by Chen and Hao [2], to evaluate and compare their IRTS algorithm; while the other is a truncating time limit, which has been widely adopted in the literature to assess different algorithms for QMKP. Our first experiment is to assess our algorithm under the relaxed time condition, where for each run, 15 s and 90 s are allowed for instances with 100 and 200 objects, respectively. The main purpose of this experiment is to understand better the search potential of our TS algorithm by permitting a longer computing time. To further evaluate the relative effectiveness and efficiency of our proposed algorithm, we also compare our TS algorithm with three currently best-performing heuristic algorithms for QMKP [2] in the literature, as well as the state-of-the-art integer programming solver CPLEX 12.4 based on the 0-1 quadratic model of QMKP formulated in (1)–(3). Since the results of the SS-ABC algorithm (artificial bee colony algorithm) [21] under long time limit are unavailable in the literature, we did not introduce it in our comparison when long run times are allowed. The three heuristic approaches for comparison are as follows:

 IRTS: Iterated responsive threshold search [2].  TIG-QMKP: Tabu-enhanced iterated greedy algorithm for QMKP [7].

 SO-QMKP: Strategic oscillation for QMKP [5]. Moreover, these three reference algorithms were recently tested and compared in [2] under the same relaxed time limit on an Intel Xeon E5440 processor (2.83 GHz CPU and 2GB RAM), which is exactly the same type of machine used to test our TS algorithm. Tables 7 and 8 summarize the computational results of our TS and the other four reference approaches on the 60 QMKP instances with d ¼0.25 and d¼ 0.75, respectively. Given the stochastic nature of the four heuristic approaches (TIG-QMKP, SOQMKP, IRTS, and our TS), each instance is independently solved 40 times with different random seeds by these heuristic approaches, with each trial being limited to the same relaxed time stopping condition as mentioned above. In Tables 7 and 8, the first four columns 1–4 give the characteristics of the instances, including the number of objects (n), density (d), number of knapsacks (k), and instance identity number. Column 5 (fbk) presents the previous best known results reported in the literature, most of which are established in [2] under exactly the same experiment protocol as the one used here. Columns 6–8 present the detailed results obtained by our TS algorithm, including the best objective value

206

J. Qin et al. / Computers & Operations Research 66 (2016) 199–214

Table 5 The post hoc test for solution sets obtained by varying the value of β. Values less than 0.05 are presented in bold. β¼

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

β ¼ 0.2 β ¼ 0.3 β ¼ 0.4 β ¼ 0.5 β ¼ 0.6 β ¼ 0.7 β ¼ 0.8 β ¼ 0.9

0.113 0.9544 0.5197 0.4101 0.6307 0.8820 0.3932 0.0616

0.4418 0.6476 0.7582 0.6540 0.3822 0.1449 0.0111

0.4599 0.2762 0.6247 0.6880 0.2527 0.0054

0.6822 0.9970 0.5370 0.2268 0.0331

0.1912 0.5009 0.0987 0.0107

0.1075 0.0162 0.0008

0.0034 4.35e 06

7.68e  06

Table 6 Statistical results obtained when varying β in the given range. We mark the best result of each row in bold. β¼

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Tot. #Best

5,022,530 52

5,022,014 51

5,022,515 52

5,022,574 54

5,022,022 54

5,022,316 54

5,022,630 57

5,021,355 57

5,020,053 55

over the 40 independent trials (Best), the average objective value over the 40 trials (Avg), and the standard deviation of the 40 objective values. In the following six columns, we present the results obtained by the three heuristic algorithms; for each algorithm, we report its overall best values and the average results over the 40 trials. The last three columns give the results obtained by the CPLEX 12.4 solver under a time limit of one hour, including the lower bound (LB), the upper bound (UB), and the gap (GAP). Note that the results of the three heuristic algorithms and the CPLEX solver are directly extracted from [2]. From Tables 7 and 8, we observe that for all the 60 instances, except the last one in Table 8, our TS algorithm under the relaxed time condition attains the previous best known results. Remarkably, for 10 out of the 60 instances, our TS algorithm is able to further improve the previous best known results. When compared with the three reference heuristic algorithms, we see that the results obtained by our TS algorithm remain competitive. Indeed, in terms of the best results achieved by these algorithms, TS attains a better result than IRTS, SO-QMKP, and TIG-QMKP for 10, 43, and 41 instances, while the converse is true only for 1, 0, and 0 instance, respectively. Moreover, the average results attained by TS are better than those obtained by IRTS, SO-QMKP, and TIG-QMKP for 37, 57, and 54 instances, whereas IRTS, SO-QMKP, and TIG-QMKP outperform TS for only 11, 1, and 1 instances, further confirming the effectiveness of our proposed algorithm. From Tables 7 and 8, we also notice that the lower bounds established by the CPLEX 12.4 solver for the 60 QMKP instances are much worse than the best results obtained by our TS algorithm. In addition, we perform non-parametric statistical tests on the experimental results to see whether there exist significant performance differences between our TS algorithm and the reference methods with respect to the best and average values. The Friedman test shows that there are significant differences among the compared methods with p values of 5.41e  24 and 2.05e  27 for the best and average results, respectively. Furthermore, the Wilcoxon test outcomes shown in Table 9 indicate that the proposed TS algorithm with long time limit yields results better than the SOQMKP and TIG-QMKP algorithms (the p values are less than 0.05) not only in terms of the best solution values, but also in terms of the average solution values. However, the Wilcoxon test outcomes disclose that in terms of the best results achieved, there is no significant performance difference (p value ¼0.0718) between the proposed TS algorithm and IRTS, though TS is able to attain a better result than IRTS for ten instances. To understand why our TS algorithm shows a better

performance for certain problem instances, we highlight two main components of the proposed algorithm that boost its performance. (1) The hybridization scheme combining both feasible and infeasible local searches plays a key role in affecting the performance of our algorithm and is particularly effective for instances with d¼ 0.75 (refer to Section 4.1 for a more detailed analysis), while IRTS restricts its search only within the feasible search space. (2) In the context of tabu search, our neighborhood union method for exploring the feasible search space is more effective than the token-ring method employed by the IRTS algorithm, especially for instances with d ¼0.25, though both TS and IRTS rely on the same four move operators (refer to Section 4.2). Despite the statement of the no free lunch (NFL) theorem [23], it has been shown recently that the real-world does not satisfy the hypothesis of this theorem [6,13], and this supports the possibility that our TS algorithm might show in general a better performance than other algorithms. 3.4. Comparison of results with a truncating time limit To further assess the performance of the proposed TS algorithm, we compare it with TIG-QMKP, SO-QMKP, and IRTS with a truncating time limit. In addition to these three comparative algorithms, we also include another recent algorithm SS-ABC (artificial bee colony algorithm) [21] in our comparison. To perform a fair comparison, we adopted exactly the same experiment protocol as that used in [2,5,7,21] and ran our algorithm 40 times on each of the 60 instances, with each run being limited to an indicated time, which was first reported in [21] and adopted in [2,5,7] for comparisons. Moreover, as indicated in [2], the computers for testing TIG-QMKP, SO-QMKP, and SS-ABC are almost as fast as the computer used to run IRTS, which is exactly the same type of machine used to test our TS algorithm. Tables 10 and 11 show the best and average results of our TS algorithm compared with the reference algorithms for the QMKP instances with d ¼0.25 and d¼ 0.75, respectively. The results of the four reference algorithms are directly extracted from the original papers [2,5,7,21]. In Tables 10 and 11, the first five columns present the characteristics of the tested instances, including the truncating time for each instance (column Time). Columns 6–16 summarize the results obtained by the five comparative algorithms. From Tables 10 and 11, we observe that in terms of the best solution found, our TS algorithm outperforms SS-ABC, TIG-QMKP, and SO-QMKP since it obtains absolutely no worse results than SSABC and gets better results than TIG-QMKP and SO-QMKP for 47

Table 7 Comparison of results of TS with respect to CPLEX, TIG-QMKP, SO-QMKP, and IRTS under a long running time limit for the QMKP instances with d ¼0.25. The best results of the compared algorithms are presented in bold. Instance

fbk d

K

I

100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 200 200 200 200 200 200 200 200 200 200 200 200 200 200 200

25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25

3 3 3 3 3 5 5 5 5 5 10 10 10 10 10 3 3 3 3 3 5 5 5 5 5 10 10 10 10 10

1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

29,286 28,491 27,179 28,593 27,892 22,581 21,704 21,239 22,181 21,669 16,221 15,700 14,927 16,181 15,326 101,471 107,958 104,589 100,098 102,311 75,623 80,033 78,043 74,140 76,610 52,293 54,830 53,661 51,297 53,621

TIG-QMKP

SO-QMKP

IRTS

Cplex

Best

Avg

SD

Best

Avg

Best

Avg

Best

Avg

LB

UB

GAP(%)

29,286 28,491 27,179 28,593 27,892 22,581 21,704 21,239 22,181 21,669 16,221 15,700 14,927 16,181 15,326 101,471 107,958 104,589 100,136 102,311 75,623 80,033 78,043 74,140 76,610 52,293 54,830 53,678 51,302 53,621

29,286.00 28,491.00 27,179.00 28,593.00 27,892.00 22,579.20 21,687.90 21,239.00 22,181.00 21,669.00 16,201.80 15,700.00 14,892.00 16,181.00 15,326.00 101,396.00 107,958.00 104,579.00 100,103.00 102,311.00 75,479.50 80,033.00 78,043.00 74,024.20 76,607.70 52,114.80 54,539.70 53,594.80 51,097.60 53,605.20

0.00 0.00 0.00 0.00 0.00 28.57 8.80 0.00 0.00 0.00 8.07 0.00 12.68 0.00 0.00 72.13 0.00 2.19 12.49 0.00 124.70 0.00 0.00 66.78 26.77 72.16 70.05 28.32 46.42 30.78

29,286 28,491 27,095 28,593 27,892 22,413 21,678 21,181 22,181 21,669 16,157 15,700 14,832 16,181 15,289 100,372 107,927 104,532 99,000 101,999 74,682 79,604 77,795 73,189 76,137 51,592 54,290 52,985 50,577 53,337

29,027.90 28,470.70 27,015.90 28,593.00 27,885.33 22,273.98 21,648.00 21,099.30 22,180.42 21,663.85 16,057.60 15,557.68 14,736.23 16,168.50 15,189.45 100,207.00 107,814.00 104,445.00 98,836.80 101,877.00 74,361.52 79,459.33 77,720.58 72,984.40 75,905.14 51,298.40 54,077.30 52,791.49 50,282.50 52,856.38

29,286 28,491 27,179 28,593 27,892 22,509 21,678 21,188 22,181 21,669 16,065 15,617 14760 16,159 15,196 101,218 107,958 104,538 99,559 102,084 74,665 79,473 77,695 73,405 76,037 51,389 54,102 52,841 50,371 52,596

29,201.70 28,488.32 27,175.20 28,580.75 27,821.98 22,403.50 21,622.43 21,153.00 22,164.32 21,567.00 15,996.83 15,446.40 14,648.43 16,082.68 15,094.89 100,776.00 107,663.00 104,365.00 99,170.50 101,792.00 74,389.82 79,244.40 77,570.82 73,005.00 75,829.90 51,043.93 53,831.20 52,483.48 50,002.82 52,395.10

29,286 28,491 27,179 28,593 27,892 22,581 21,704 21,239 22,181 21,669 16,221 15,700 14,927 16,181 15,326 101,471 107958 104,589 100,098 102,311 75,623 80,033 78,043 74,140 76,610 52,293 54,830 53,661 51,297 53,621

29,286.00 28,491.00 27,179.00 28,593.00 27,892.00 22,530.68 21,667.00 21,235.95 22,180.90 21,656.42 16,200.53 15,665.65 14,852.00 16,181.00 15,293.00 101,441.00 107,958.00 104,559.00 100,098.00 102,310.00 75,554.10 80,023.40 78,028.95 74,061.29 76,597.62 52,158.50 54,666.25 53,588.28 51,078.20 53,532.24

28,077 28,169 26,492 27,793 27,058 21,194 20,725 19,674 20,644 20,054 14,804 14,191 13,560 14,630 14,142 94,992 105,978 98,214 93,693 94,818 67,464 71,876 70,259 65,940 69,523 43,054 45,774 44,187 41,608 43,847

53,963.23 52,942.22 51,010.72 53,382.60 53,083.61 55,784.57 54,647.80 52,614.27 55,098.88 54,887.58 57,710.44 56,294.00 54,274.15 56,871.67 56,609.11 222,051.26 222,937.93 219,794.08 219,675.07 217,794.90 226,198.30 226,791.30 223,882.05 223,859.76 222,056.15 230,387.32 231,059.61 227,911.74 228,108.72 226,045.81

92.2 87.94 92.55 92.07 96.18 163.21 163.68 167.43 166.9 173.7 289.83 296.69 300.25 288.73 300.29 133.76 110.36 123.79 134.46 129.7 235.29 215.53 218.65 239.49 219.4 435.11 404.78 415.79 448.23 415.53

J. Qin et al. / Computers & Operations Research 66 (2016) 199–214

n

TS

207

Instance

fbk d

K

I

100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 200 200 200 200 200 200 200 200 200 200 200 200 200 200 200

75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75

3 3 3 3 3 5 5 5 5 5 10 10 10 10 10 3 3 3 3 3 5 5 5 5 5 10 10 10 10 10

1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

69,977 69,504 68,832 70,028 69,692 49,421 49,365 48,495 50,246 48,753 30,296 31,207 29,908 31,762 30,507 270,718 257,288 270,069 246,993 279,598 185,493 174,836 186,774 166,990 193,310 113,139 105,807 114,596 99,106 117,309

TIG-QMKP

SO-QMKP

IRTS

Cplex

Best

Avg

SD

Best

Avg

Best

Avg

Best

Avg

LB

UB

GAP(%)

69,977 69,504 68,832 70,028 69,692 49,421 49,400 48,495 50,246 48,753 30,296 31,207 29,908 31,762 30,507 270,718 257,288 270,069 246,993 279,598 185,493 174,836 186,782 167,142 193,310 113,324 105,959 114,860 99,309 117,110

69,977.00 69,504.00 68,832.00 70,028.00 69,692.00 49,421.00 49,386.00 48,495.00 50,246.00 48,753.00 30,232.40 31,056.70 29,896.30 31,710.90 30,450.90 270,681.00 257,288.00 270,069.00 246,881.00 279,598.00 185,443.00 174,836.00 186,732.00 166,918.00 193,263.00 113,061.00 105,800.00 114,599.00 98,850.30 116,940.00

0.00 0.00 0.00 0.00 0.00 0.00 12.50 0.00 0.00 0.00 57.07 34.63 9.84 19.80 30.50 75.78 0.00 0.00 64.16 0.00 34.35 0.00 31.58 79.24 46.27 72.72 86.24 89.75 112.76 94.19

69,935 69,504 68,832 70,028 69,692 49,421 49,360 48,495 50,246 48,752 30138 31092 29812 31672 30188 270,718 257,288 270,069 246,882 279598 184,984 174,776 186,674 166,832 193,255 112,591 105,297 114,237 98,556 116,725

69,935.00 69,504.00 68,816.20 70,028.00 69,681.30 49,295.60 49,266.80 48,474.20 49,966.60 48,735.20 29,900.84 30,969.15 29,662.00 31,491.82 30,046.10 270,718.00 257,099.00 270,069.00 246,684.00 279,598.00 184,774.00 174,642.00 186,507.00 166,487.00 193,002.00 112,330.00 105,064.00 113,930.00 98,219.93 116,266.00

69,935 69,504 68,832 70,028 69,692 49,363 49,320 48,495 50,246 48,752 30,018 30,973 29,765 31,634 30,348 270,718 257277 270,069 246,882 279,598 184,882 174,682 186,619 166,584 193,138 112,457 105,260 114,007 98,285 116,298

69,935.00 69,497.40 68,813.00 70,028.00 69,652.22 49,238.83 49,226.60 48,360.85 50,124.20 48,718.38 29,897.80 30,914.00 29,638.80 31,481.30 30,055.42 270,697.00 256,931.00 270,028.00 246,555.00 279,598.00 184,641.00 174,445.00 186,352.00 166,246.00 192,836.00 112,258.00 104,947.00 113,717.00 97,885.95 116,031.00

69,977 69,504 68,832 70,028 69,692 49,421 49,365 48,495 50,246 48,753 30,296 31,207 29,908 31,762 30,507 270,718 257288 270,069 246,993 279,598 185,493 174,836 186,774 166,990 193,310 113,139 105,807 114,596 99,106 117,309

69,977.00 69499.6 68,832.00 70,028.00 69692 49,365.98 49,350.60 48,495.00 50,141.50 48,749.10 30,240.20 31,095.80 29,894.75 31,706.50 30,458.50 270,685.00 257,273.00 269,926.00 246,877.00 279,570.00 184,904.00 174,688.00 186,674.00 166,747.00 193,217.00 112,809.00 105,437.00 114,367.00 98,851.55 116,947.00

69,010 68,157 67,681 69,717 68,638 48,270 48,643 44,474 48,756 47,286 28,124 29,436 27,340 27,916 27,003 258,740 248,139 262,806 227,193 270,979 164,519 160,577 175,943 158,807 173,698 95,514 88,469 94,768 88,580 98,288

157,543.51 158,390.17 157,395.07 154,667.08 160,393.28 161,306.19 162,654.44 161,403.85 159,342.88 164,701.44 165,782.37 167,227.44 165,637.23 163,759.49 168,992.47 644,499.87 645,348.35 649,880.92 633,617.30 655,098.56 654,784.78 656,046.05 661,330.18 643,762.33 665,834.00 665,339.06 667,546.03 672,239.84 653,941.84 676,442.05

128.29 132.39 132.55 121.85 133.68 234.17 234.38 262.92 226.82 248.31 489.47 468.11 505.84 486.62 525.83 149.09 160.08 147.29 178.89 141.75 298 308.56 275.88 305.37 283.33 596.59 654.55 609.35 638.25 588.22

J. Qin et al. / Computers & Operations Research 66 (2016) 199–214

n

TS

208

Table 8 Comparison of results of TS with respect to CPLEX, TIG-QMKP, SO-QMKP, and IRTS under a long running time limit for the QMKP instances with d ¼0.75. The best results of the compared algorithms are presented in bold.

J. Qin et al. / Computers & Operations Research 66 (2016) 199–214

209

Table 9 Wilcoxon test of the best and average results obtained with a long time limit. Algorithm pair

TS TS TS TS

versus versus versus versus

Best results

CPLEX TIG-QMKP SO-QMKP IRTS

Average results

Rþ

R

p - value

Diff?

Rþ

R

p - value

Diff?

1830.00 1735.00 1753.50 1159.50

0.00 95.00 76.50 670.50

0.0000 0.0000 0.0000 0.0718

Yes Yes Yes No

– 1807.50 1821.50 1274.50

– 22.50 8.50 555.50

– 0.0000 0.0000 0.0142

– Yes Yes Yes

Table 10 Comparison of results of TS with respect to SS-ABC, TIG-QMKP, SO-QMKP, and IRTS under a truncating time limit for the QMKP instances with d ¼ 0.25. The best results of the compared algorithms are presented in bold. Instance

SS-ABC

TIG-QMKP

SO-QMKP

IRTS

TS

n

d

K

I

Time

Best

Avg

Best

Avg

Best

Avg

Best

Avg

Best

Avg

SD

100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 200 200 200 200 200 200 200 200 200 200 200 200 200 200 200

25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25

3 3 3 3 3 5 5 5 5 5 10 10 10 10 10 3 3 3 3 3 5 5 5 5 5 10 10 10 10 10

1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

3.8 3.69 2.87 3.74 3.57 2.7 2.93 2.4 3.26 3.18 1.42 1.84 1.58 2.1 1.82 23.99 18.61 29.85 38.93 29.22 19.88 16.75 22.86 28.07 20.74 10.37 8.48 11.15 12.83 10.99

29,139 28,443 26,901 28,568 27,849 22,390 21,584 21,093 22,178 21,301 15,953 15,487 14,339 15,807 14,719 100,662 107,958 104,521 98,791 102,049 74,922 79,506 77,607 73,181 76,022 49,883 53,298 52,281 49,210 51,921

28,753.00 28,004.00 26,585.33 28,109.03 27,073.67 22,117.12 21,224.03 20,771.12 21,767.50 20,875.47 15,573.65 14,896.35 14,027.83 15,397.00 14,376.80 100,103.02 107,545.20 104,006.98 98,344.32 101,406.48 74,132.95 79,073.32 77,069.52 72,607.25 75,455.98 49,079.47 51,831.55 51,324.28 48,190.60 51,437.97

29,138 28,473 27,013 28,593 27,892 22,264 21,580 21,100 22,180 21,669 16,118 15,525 14,773 16,181 15,150 100,218 107,787 104,479 98,896 101,973 74,239 79,480 77,700 73,173 75,884 51,413 54,116 52,735 50,221 52,651

28,894.60 28,224.35 26,905.40 28,573.60 27,721.68 22,126.00 21,430.38 21,015.03 22,043.00 21,397.58 15,863.00 15,398.43 14,554.00 16,089.95 15,023.45 100,056.23 107,644.98 104,251.50 98,557.40 101,635.43 73,977.78 79,234.28 77,420.50 72,477.65 75,693.18 50,845.78 53,608.45 52,456.28 49,656.40 52,328.38

29,234 28,491 27,179 28,593 27,892 22,509 21,678 21,188 22,181 21,669 16,065 15,510 14,663 16,159 15,130 101,100 107,805 104,538 99,559 102,041 74,559 79,400 77,632 73,327 75,996 51,323 53,975 52,841 50,190 52,470

29,159.25 28,467.00 27,155.35 28,570.40 27,811.53 22,337.43 21,540.35 21,104.75 22,136.00 21,462.88 15,886.58 15,359.95 14,568.38 16,013.00 15,021.33 100,653.50 107,607.15 104,271.68 99,003.63 101,667.73 74,237.40 79,153.55 77,452.25 72,884.03 75,751.38 50,862.70 53,649.03 52,337.73 49,802.43 52,211.58

29,286 28,491 27,179 28,593 27,892 22,581 21,678 21,239 22,181 21,669 16,213 15,700 14,860 16,181 15,326 101,445 107,958 104,575 100,098 102,311 75,596 80,033 78,043 74,111 76,610 52,115 54,716 53,646 51,176 53,568

29,276.80 28,491.00 27,159.15 28,591.72 27,880.60 22,463.43 21,614.00 21,170.90 22,168.95 21,606.80 16,144.13 15,544.30 14,757.24 16,160.00 15,187.25 101,410.00 107,957.00 104,544.00 100,098.00 102,307.00 75,491.48 79,914.90 77,992.93 73,999.22 76,532.10 51,862.20 54,454.75 53,401.80 50,832.48 53,406.60

29,286 28,491 27,179 28,593 27,892 22,581 21,704 21,239 22,181 21,669 16,221 15,700 14,927 16,181 15,326 101,465 107,958 104,589 100,098 102,311 75,623 80,033 78,043 74,140 76,610 52,293 54,688 53,661 51,302 53,616

29,255.40 28,491.00 27,167.50 28,593.00 27,881.30 22,522.10 21,679.70 21,209.30 22,174.00 21,658.80 16,191.10 15,649.20 14,862.40 16,175.80 15,260.50 101,306.00 107,957.00 104,559.00 99,821.90 102,254.00 75,261.20 79,945.00 77,927.80 73,861.80 76,461.80 51,880.80 54,221.50 53,343.10 50,984.90 53,428.40

53.65 0.00 25.58 0.00 28.26 37.39 22.94 47.46 18.46 30.54 16.89 50.39 31.31 10.98 22.36 206.26 8.90 16.86 433.60 281.88 195.31 228.12 60.32 186.87 135.42 172.08 288.25 256.38 134.93 201.43

and 45 instances, respectively. Furthermore, our algorithm also compares very favorably with IRTS, which seems to be currently the most effective method for QMKP, since our TS algorithm obtains better results than IRTS for 20 instances, while worse results for only 8 instances. Considering the average results, our TS algorithm remains competitive with SS-ABC, TIG-QMKP, and SOQMKP. Indeed, our TS algorithm is able to attain better results than SS-ABC, TIG-QMKP, and SO-QMKP for 60, 49, and 48 instances, respectively, and yields worse results than SS-ABC, TIG-QMKP, and SO-QMKP for only 0, 11, and 12 instances, respectively. However, TS is outperformed by IRTS in terms of the average results since it obtains better results than IRTS for 26 instances while worse results than IRTS for 32 instances. To further validate our conclusion, we apply the Wilcoxon test to check the statistical differences between the TS algorithm and the compared approaches and report the results in Table 12. Table 12 indicates that the differences between TS and the four reference algorithms are statistically significant. Furthermore, the Wilcoxon test conducted on the average results obtained by the five algorithms also clearly indicates that there are statistically

significant differences among the average results obtained by the five methods under comparison.

4. Analysis of the main components of TS In this section we investigate two important components of the proposed TS algorithm that contribute to its effectiveness, namely the hybridization scheme and the neighborhood combination method. 4.1. Effect of the hybridization scheme combining feasible and infeasible local searches Most of the heuristic approaches for QMKP in the literature like IRTS [2], TIG [7], and SS-ABC [21] explore the search space of QMKP only in the feasible search space. Compared with these heuristic approaches, one salient feature of our proposed TS algorithm for solving QMKP is hybridization of the feasible and infeasible local search methods. By integrating both the feasible and infeasible local

210

J. Qin et al. / Computers & Operations Research 66 (2016) 199–214

Table 11 Comparison of results of TS with respect to SS-ABC, TIG-QMKP, SO-QMKP, and IRTS under a truncating time limit for the QMKP instances with d ¼ 0.75. The best results of the compared algorithms are presented in bold. Instance

SS-ABC

TIG-QMKP

SO-QMKP

IRTS

TS

n

d

K

I

Time

Best

Avg

Best

Avg

Best

Avg

Best

Avg

Best

Avg

SD

100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 200 200 200 200 200 200 200 200 200 200 200 200 200 200 200

75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75

3 3 3 3 3 5 5 5 5 5 10 10 10 10 10 3 3 3 3 3 5 5 5 5 5 10 10 10 10 10

1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

2.07 1.86 1.86 1.88 2.12 2.07 1.96 1.71 1.83 2.01 1.22 1.45 1.3 1.4 1.42 14.11 16.27 11.87 30.64 10.46 12.34 12.34 12.1 27.03 14.06 7.64 9.96 8.04 14.81 8.21

69,721 69,462 68,635 69,986 69,679 4,922 49,313 48,472 50,199 48,710 29,875 30,939 29,465 31,663 30,219 269,736 256,195 268,890 246,205 279,490 184,448 173,575 185,107 165,273 192,764 111,000 103,540 112,509 96,859 115,125

69,373.00 69,041.00 67,960.05 69,687.68 69,136.40 48,937.47 48,908.05 47,874.50 50,017.93 48,409.75 29,429.20 30,697.80 28,983.78 31,218.85 29,736.47 267,117.92 253,916.75 267,079.03 244,618.40 276,605.00 183,046.65 171,738.85 185,059.52 164,042.20 190,474.27 109,624.73 102,603.18 111,388.20 95,681.70 113,909.60

69,977 69,504 68,811 70,028 69,692 49,345 49,316 48,495 49,866 48,752 29,877 30,980 29,695 31,550 30,096 270,718 257,090 270,069 246,704 279,598 184,909 174,682 186,443 166,358 193,084 112,262 105,092 113,868 98,252 116,513

69,936.05 69,442.78 68,811.00 70,019.88 69,638.48 49,218.10 49,081.53 48,327.48 49,866.00 48,619.60 29,548.68 30,832.15 29,439.95 31,333.45 29,895.40 270,525.90 256,764.98 269,974.03 246,356.53 279,572.30 184,500.80 174,239.48 186,170.68 166,159.55 192,712.25 111,889.75 104,669.83 113,510.55 97,807.73 115,856.30

69,935 69,504 68,832 70,028 69,653 49,363 49,305 48,495 50,246 48,752 29,931 30,973 29,730 31,587 30,229 270,718 257,026 270,069 246,882 279,598 184,822 174,682 186,526 166,584 193,053 112,354 105,151 113,869 98,028 116,298

69,925.73 69,442.48 68,812.58 70,028.00 69,640.53 49,197.53 49,137.38 48,287.88 50,025.03 48,653.18 29,788.48 30,829.05 29,519.48 31,392.48 29,918.70 270,617.48 256,852.30 269,955.03 246,473.13 279,562.43 184,529.00 174,267.00 186,216.75 166,165.38 192,702.25 112,043.68 104,781.50 113,563.08 97,747.55 115,807.53

69,977 69,504 68,832 70,028 69,692 49,421 49,365 48,495 50,246 48,753 30,290 31,079 29,908 31,682 30,465 270,718 257,288 270,069 246,993 279,598 185,221 174,836 186,678 166,990 193,253 112,834 105,807 114,567 99,006 117,092

69,970.13 69,465.50 68,823.00 70,028.00 69,674.60 49,310.80 49,311.80 48,479.60 49,935.17 48,689.20 30,023.00 30,909.93 29,794.85 31,576.68 30,211.50 270,518.00 257,117.00 269,796.00 246,658.00 279,548.00 184,801.00 174,499.00 186,510.00 166,602.00 193,133.00 112,443.00 105,241.00 114,084.00 98,664.60 116,649.00

69,977 69,504 68,832 70,028 69,692 49,421 49,400 48,495 50,246 48,753 30,290 31,207 29,897 31,762 30,507 270,718 257,288 269,744 246,992 279,598 184,810 174,836 186,774 167,142 193,240 113,107 105,790 114,657 99,161 116,953

69,768.80 69,473.00 68,832.00 69,987.10 69,633.10 49,252.40 49,366.80 48,491.60 50,246.00 48,688.60 30,028.60 31,036.70 29,636.30 31,539.60 30,120.70 270,105.00 257,190.00 268,788.00 246,470.00 279,060.00 183,937.00 173,895.00 185,959.00 166,577.00 192,510.00 112,296.00 105,089.00 113,487.00 98,499.80 115,916.00

131.84 25.68 0.00 0.00 62.19 9.97 14.42 21.39 0.00 66.16 128.65 81.86 174.44 225.98 251.14 387.14 149.63 652.16 267.85 515.83 440.76 894.17 265.06 303.49 373.34 475.10 419.93 446.17 260.04 363.04

Table 12 Wilcoxon test for the results obtained with truncated time. Algorithm pair

TS TS TS TS

versus versus versus versus

SS-ABC TIG SO IRTS

Best results

Average results

Rþ

R

p - value

Diff?

Rþ

R

p - value

Diff?

1829.50 1742.00 1728.50 1207.00

0.50 88.00 101.50 623.00

0.0000 0.0000 0.0000 0.0316

Yes Yes Yes Yes

1830.00 1544.00 1492.00 561.50

0.00 286.00 338.00 1268.50

0.0000 0.0000 0.0000 0.0094

Yes Yes Yes Yes

searches into a same search procedure, our TS algorithm is able to explore the search space of QMKP more effectively. In this section we report additional experiments to examine the effect of the hybridization scheme on the performance of our TS algorithm. For this purpose, we compare the performance of the TS algorithm with its two underlying local search components, i.e., the feasible local search procedure and the infeasible local search procedure. To make a fair comparison, we ran each of the three algorithms 40 times on the set of 60 QMKP instances with exactly the same timeout limit, which was set to be 15 s for instances with 100 objects and 90 s for instances with 200 objects. Table 13 reports the results for all 60 tested instances. From Table 13, we observe that our TS algorithm shows an overall better performance than its two underlying FLS and ILS local search procedures in terms of both the best and average solution values for the 60 tested QMKP instances. Especially for most of the instances with d¼0.75, our TS is able to yield better best and average objective values than FLS and ILS. When applying the Wilcoxon test to check the statistical differences between TS and its underlying algorithms in terms of both the best and average solution values, we obtain all p values o 0:05 (Table 14),

indicating that the observed differences between TS and its two underlying local search methods are statistically significant. In view of these observations, we confirm that the hybridization scheme plays a key role in boosting the overall performance of our TS. To augment these observations, we illustrate in Fig. 2 the absolute gaps between the average results obtained by TS, FLS, and ILS, respectively, and the best known results reported in the literature for each of the 60 tested instances. Fig. 2 also clearly discloses that our TS algorithm yields significantly better results than its two underlying local search methods. Furthermore, we also analyze the distribution of objects in knapsacks for some high quality solutions generated during the search process. Based on this analysis, we explain to some extent why our hybridization scheme is effective in tackling some QMKP instances with d ¼0.75. We find that during the search process of the feasible local search, some light weight objects tend to be grouped together to form several (one or two) large knapsacks with many objects (see Table 15 for the distribution of some high quality solutions generated during the search process), since the more objects grouped in a same knapsack, the

J. Qin et al. / Computers & Operations Research 66 (2016) 199–214

211

Table 13 Comparison of results among TS, FLS and ILS. Instance

TS

FLS

ILS

n

d

K

I

Best

Avg

Best

Avg

Best

Avg

100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 200 200 200 200 200 200 200 200 200 200 200 200 200 200 200 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 200 200 200 200 200 200 200 200 200 200 200 200 200 200 200

25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75

3 3 3 3 3 5 5 5 5 5 10 10 10 10 10 3 3 3 3 3 5 5 5 5 5 10 10 10 10 10 3 3 3 3 3 5 5 5 5 5 10 10 10 10 10 3 3 3 3 3 5 5 5 5 5 10 10 10 10 10

1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

29,286 28,491 27,179 28,593 27,892 22,581 21,704 21,239 22,181 21,669 16,221 15,700 14,927 16,181 15,326 101,471 107,958 104,589 100,136 102,311 75,623 80,033 78,043 74,140 76,610 52,293 54,830 53,678 51,302 53,621 69,977 69,504 68,832 70,028 69,692 49,421 49,400 48,495 50,246 48,753 30,296 31,207 29,908 31,762 30,507 270,718 257,288 270,069 246,993 279,598 185,493 174,836 186,782 167,142 193,310 113,324 105,959 114,860 99,309 117,110

29,286.00 28,491.00 27,179.00 28,593.00 27,892.00 22,579.20 21,687.90 21,239.00 22,181.00 21,669.00 16,201.80 15,700.00 14,892.00 16,181.00 15,298.50 101,396.00 107,958.00 104,579.00 100,042.00 102,311.00 75,479.50 80,033.00 77,999.10 74,024.20 76,607.70 52,114.80 54,539.70 53,594.80 51,097.60 53,605.20 69,977.00 69,504.00 68,832.00 70,028.00 69,692.00 49,421.00 49,386.00 48,495.00 50,246.00 48,753.00 30,232.40 31,056.70 29,896.30 31,710.90 30,450.90 270,681.00 257,288.00 270,069.00 246,881.00 279,598.00 185,443.00 174,836.00 186,732.00 166,918.00 193,263.00 113,061.00 105,800.00 114,599.00 98,850.30 116,940.00

29,286 28,491 27,179 28,593 27,892 22,581 21,704 21,239 22,181 21,669 16,221 15,700 14,927 16,181 15,326 101,471 107,958 104,589 100,136 102,311 75,623 80,033 78,043 74,140 76,610 52,293 54,830 53,678 51,302 53,621 69,977 69,504 68,832 69,986 69,692 49,421 49,400 48,495 50,246 48,753 30,290 31,191 29,908 31,634 30,507 269,866 257,013 269,671 246,757 279,331 184,977 174,836 186,337 167,142 192,802 113,304 105,794 114,019 99,072 116,744

29,286.00 28,491.00 27,178.00 28,593.00 27,892.00 22,554.40 21,698.40 21,234.70 22,179.70 21,669.00 16,207.00 15,687.80 14,900.10 16,181.00 15,311.20 101,329.00 107,958.00 104,548.00 99,883.90 102,311.00 75,372.00 79,954.00 77,982.20 73,939.60 76,522.40 52,151.50 54,554.80 53,539.10 51,115.20 53,578.60 69,781.70 69,482.30 68,832.00 69,860.00 69,684.20 48,937.70 49,365.00 48,495.00 50,246.00 48,718.40 29,981.10 31,032.10 29,545.00 31,117.20 30,207.80 269,009.00 256,915.00 268,398.00 246,181.00 276,404.00 184,191.00 173,943.00 185,008.00 166,524.00 191,989.00 112,451.00 105,206.00 112,924.00 98,450.40 115,686.00

29,286 28,491 27,179 28,593 27,892 22,581 21,704 21,239 22,181 21,669 16,213 15,700 14,893 16,181 15,326 101,430 107,958 104,556 100,098 102,311 75,552 80,033 77,967 74,082 76,610 52,164 54,672 53,618 51,219 53,621 69,977 69,504 68,832 70,028 69,692 49,421 49,400 48,495 50,246 48,753 30,290 31,207 29,902 31,762 30,507 270,703 257,288 269,761 246,992 279,598 185,276 174,836 186,782 167,142 193,240 113,324 105,829 114,671 99,047 116,885

29,214.60 28,491.00 27,138.20 28,593.00 27,887.00 22,534.60 21,677.50 21,198.20 22,177.20 21,662.90 16,191.50 15,649.80 14,849.10 16,178.40 15,268.10 101,298.00 107,957.00 104,540.00 99,387.10 102,249.00 75,272.90 80,027.60 77,898.40 73,801.80 76,475.50 51,973.20 54,551.00 53,471.40 51,013.00 53,524.30 69,977.00 69,504.00 68,832.00 70,021.70 69,692.00 49,317.30 49,345.60 48,495.00 50,246.00 48,752.10 30,160.00 31,106.30 29,871.20 31,619.30 30,416.80 270,655.00 257,288.00 269,669.00 246,992.00 279,594.00 184,724.00 174,836.00 186,604.00 166,944.00 192,960.00 113,084.00 105,659.00 114,387.00 98,759.80 116,646.00

Table 14 Wilcoxon test for TS and its two underlying local searches. Algorithm pair

TS versus FLS TS versus ILS

Best results

Average results

Rþ

R

p - value

Diff?

Rþ

R

p - value

Diff?

1357.00 1478.50

473.00 351.50

0.0014 0.0000

Yes Yes

1652.00 1650.00

178.00 180.00

0.0000 0.0000

Yes Yes

212

J. Qin et al. / Computers & Operations Research 66 (2016) 199–214

more non-zero pairwise profits pij will be generated by this knapsack (i.e., the uneven distribution of the objects in knapsacks will generate more effective non-zero pairwise profits, generally leading to a larger total profit, especially for the instances with d ¼0.75). Moreover, the quality of the solution depends critically on these large knapsacks because the profits of these knapsacks contribute significantly to the overall profit of the solution. However, during the search process of the feasible local search, it is very difficult to change the objects (by the four move operators) in these large knapsacks because the removal of any object in them (by the extraction or re-allocation moves) will drastically deteriorate the solution (the removal of an object in a large knapsack will lead to a fast loss in many joint pairwise profits, especially for instances with d¼0.75). This will also prohibit the application of some other highly profitable moves (such as the allocation and exchange moves) to these large knapsacks because, due to the capacity constraints, the application of these moves requires the first removal of objects in the knapsack to make room for new objects. Indeed, we observe in our experiments that for the feasible local search method, few moves are applied to these critical large knapsacks. On the other hand, our hybridization scheme can help our algorithm overcome this shortcoming and apply more effective moves to these critical knapsacks. Finally, we also collect statistics about the proportion of infeasible solutions among all the solutions generated during the search process of our algorithm. Our experiment results show that

3500

TS FLS ILS

Absolute gap to the best result

3000 2500 2000 1500 1000 500 0

0

10

20

30

40

50

60

Instance

Fig. 2. Absolute gaps between the average and the best solution results.

the proportion of infeasible solutions among all the visited solutions is roughly 18.46%, 21.54%, and 22.69% for instances with 3, 5, and 10 knapsacks, respectively. 4.2. Influence of neighborhood combination In order to effectively explore the search space, our proposed FLS algorithm explores four neighborhoods induced by four types of moves, which have been frequently used in previous studies [2,5]. In the context of multiple neighborhood search, there are several ways to combine different neighborhoods. For instance, the four neighborhoods are examined iteratively via the token-ring method in [2], while for our FLS, the four neighborhoods are explored as a union [27]. In this section we investigate the influence of different ways for exploring these neighborhoods on the performance of the FLS algorithm. For this purpose, we compare our neighborhood union method with the token-ring search. In order to isolate the effect of the neighborhood and highlight the difference between the two neighborhood exploring methods, we use only the feasible local search component of our TS algorithm in this experiment. For the token-ring search, we use the same method as that used in [2], where the four neighborhoods are examined iteratively in a sequential manner, i.e., ND (corresponding to N1 in our TS) -CN R (corresponding to N 2 [ N 3 in our TS) -CN E (corresponding to N4 in our TS) -N D -CN R -CN E …. For each neighborhood under consideration, we examine all the neighboring solutions and select the best admissible (non-tabu or globally improving) solution that yields the largest move gain. Furthermore, since the frequent application of the extraction move can rapidly deteriorate the quality of the current solution, the solutions in ND are accepted only with a low probability of 5%. Table 16 shows the statistical results on the solution sets generated with two versions of our feasible local search algorithms, which integrate the neighborhood union method (denoted by FLS_UN) and the token-ring search (denoted by FLS_TR), respectively. To make a fair comparison, we ran each version of the algorithm 40 times on the set of 60 QMKP instances with exactly the same timeout limit, which was set to be 15 s for instances with 100 objects and 90 s for those with 200 objects. From Table 16, we see that on the 60 tested benchmarks, FLS_UN shows an overall better performance than FLS_TR. In terms of the best results achieved by the two algorithms, FLS_UN reaches better results than FLS_TR for 16 instances, while the reverse is true only for four instances. Furthermore, we also notice that our FLS_UN algorithm is especially effective for instances with d¼ 0.25. Indeed, for each of the 30 instances with d ¼0.25, FLS_UN is able to produce equal or even better results in terms of the best results

Table 15 Distribution of objects in knapsacks for some high quality solutions. Columns 5–14 present the number of objects in the knapsack concerned and the corresponding profit gathered. Instance

knapsack

n

d

K

I

knapsack1

knapsack2

knapsack3

knapsack4

knapsack5

knapsack6

knapsack7

knapsack8

knapsack9

knapsack10

100 100 100 200 200 200 100 100 100 200 200 200

25 25 25 25 25 25 75 75 75 75 75 75

3 5 10 3 5 10 3 5 10 3 5 10

1 1 1 1 1 1 1 1 1 1 1 1

39(14,259) 24(6002) 16(3588) 83(56,165) 61(56,165) 37(14,009) 51(51,356) 39(30,498) 26(14,301) 104(207,539) 79(120,552) 54(58,367)

27(8734) 20(5784) 12(2539) 51(26,236) 34(14,595) 19(5794) 21(10,798) 15(6688) 12(3954) 38(33,467) 28(20,346) 22(13,041)

22(6293) 17(4551) 9(1738) 43(19,017) 28(10,478) 17(4784) 17(7823) 13(4947) 9(2610) 36(29,712) 24(15,357) 17(8301)

– 14(3352) 9(1688) – 26(9126) 16(4238) – 11(3528) 8(2214) – 24(14,866) 13(5225)

– 13(2820) 8(1444) – 26(8584) 15(4432) – 11(3706) 7(1747) – 23(14,372) 12(5155)

– – 8(1169) – – 15(4328) – – 6(1435) – – 12(4719)

– – 7(1182) – – 15(4158) – – 5(1090) – – 12(4702)

– – 7(978) – – 15(3772) – – 5(940) – – 12(4607)

– – 6(942) – – 14(3534) – – 5(932) – – 12(4606)

– – 5(953) – – 12(2762) – – 5(905) – – 12(4583)

J. Qin et al. / Computers & Operations Research 66 (2016) 199–214

Table 16 Comparison of results of FLS_UN and FLS_TR for various tested QMKP instances. Instance

FLS_UN

FLS_TR

n

d

K

I

Best

Avg

Best

Avg

100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 200 200 200 200 200 200 200 200 200 200 200 200 200 200 200 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 200 200 200 200 200 200 200 200 200 200 200 200 200 200 200

25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75

3 3 3 3 3 5 5 5 5 5 10 10 10 10 10 3 3 3 3 3 5 5 5 5 5 10 10 10 10 10 3 3 3 3 3 5 5 5 5 5 10 10 10 10 10 3 3 3 3 3 5 5 5 5 5 10 10 10 10 10

1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

29,286 28,491 27,179 28,593 27,892 22,581 21,704 21,239 22,181 21,669 16,221 15,700 14,927 16,181 15,326 101,471 107,958 104,589 100,136 102,311 75,623 80,033 78,043 74,140 76,610 52,293 54,830 53,678 51,302 53,621 69,977 69,504 68,832 69,986 69,692 49,421 49,400 48,495 50,246 48,753 30,290 31,191 29,908 31,634 30,507 269,866 257,013 269,671 246,757 279,331 184,977 174,836 186,337 167,142 192,802 113,304 105,794 114,019 99,072 116,744

29,286.00 28,491.00 27,178.00 28,593.00 27,892.00 22,554.40 21,698.40 21,234.70 22,179.70 21,669.00 16,207.00 15,687.80 14,900.10 16,181.00 15,311.20 101,329.00 107,958.00 104,548.00 99,883.90 102,311.00 75,372.00 79,954.00 77,982.20 73,939.60 76,522.40 52,151.50 54,554.80 53,539.10 51,115.20 53,578.60 69,781.70 69,482.30 68,832.00 69,860.00 69,684.20 48,937.70 49,365.00 48,495.00 50,246.00 48,718.40 29,981.10 31,032.10 29,545.00 31,117.20 30,207.80 269,009.00 256,915.00 268,398.00 246,181.00 276,404.00 184,191.00 173,943.00 185,008.00 166,524.00 191,989.00 112,451.00 105,206.00 112,924.00 98,450.40 115,686.00

29,286 28,491 27,179 28,593 27,892 22,581 21,704 21,239 22,181 21,669 16,221 15,700 14,927 16,181 15,326 101,285 107,958 104,589 99,635 102,311 75,623 80,033 78,043 74,050 76,610 52,293 54,830 53,646 51,215 53,621 69,778 69,504 68,832 69,986 69,692 49,421 49,400 48,495 50,246 48,753 30,290 31,207 29,908 31,634 30,507 268,987 257,288 268,399 246,757 278,370 184,810 174,617 186,337 166,843 192,741 112,874 105,729 114,289 99,193 116,453

29,245.90 28,491.00 27,179.00 28,593.00 27,892.00 22,542.00 21,691.20 21,239.00 22,180.60 21,669.00 16,201.00 15,700.00 14,897.60 16,180.00 15,291.10 100,968.00 107,948.00 104,512.00 99,347.10 102,053.00 75,096.50 79,990.20 77,961.50 73,781.30 76,488.90 52,139.50 54,509.10 53,527.20 51,049.20 53,528.90 69,778.00 69,476.90 68,832.00 69,789.10 69,679.60 49,185.70 49,364.80 48,467.20 50,246.00 48,739.30 29,989.50 31,038.50 29,724.60 31,418.80 30,349.90 267,051.00 255,681.00 267,856.00 243,111.00 276,174.00 183,349.00 174,400.00 185,756.00 166,540.00 192,162.00 112,444.00 105,252.00 113,160.00 98,769.20 115,941.00

achieved than the current best performing heuristic approach IRTS. In addition, the results of the Wilcoxon test also indicate that FLS_UN performs significantly better than FLS_TR in terms of the best values at a p value of 0.0214. The above observation confirms that in the context of our feasible local search, our method for exploring the search space in a neighborhood union manner constitutes an important feature to enhance the performance of the algorithm.

213

5. Conclusion One key issue in designing heuristics to solve strictly constrained combinatorial optimization problems is how to effectively explore the search space so as to strengthen the search ability of the algorithm. In this paper we present a tabu search algorithm that integrates both feasible and infeasible local searches for QMKP. By taking advantage of both the global character of infeasible local search and the more intensive focus typically provided by feasible local search, the proposed tabu search offers the possibility of striking a balance between intensification and diversification within the search procedure. We extensively evaluate the performance of the proposed tabu search algorithm with both short and long run times in solving a set of 60 well-known QMKP benchmark instances. The computational results generated in a short computing time demonstrate that our approach compares very favorably with the current bestperforming algorithms for QMKP in the literature. Moreover, when the running time is prolonged to 15 and 90 seconds for instances with 100 and 200 objects, respectively, our approach succeeds even in improving the previous best known results for ten out of the 60 instances. We also provide an analysis to show the merit of the hybridization of the feasible and infeasible local search schemes within a same search procedure by comparing it with its underlying feasible and infeasible local search components. The outcomes suggest that the hybridization scheme plays a key role in contributing to the effectiveness of the proposed TS algorithm. For future research, it is worth considering to apply our hybrid method of combining both feasible and infeasible local searches to deal with other strictly constrained combinatorial optimization problems such as the max-bisection problem [25].

Acknowledgment We thank the anonymous referees for their helpful comments on an earlier version of our paper. This work was supported in part by the National Natural Science Foundation Program of China under grant numbers 71131004 and 71401059. References [1] Billionnet A, Soutif E. An exact method for the 0-1 quadratic knapsack problem based on Lagrangian decomposition. Eur J Oper Res 2004;157(3):565–75. [2] Chen Y, Hao JK. Iterated responsive threshold search for the quadratic multiple knapsack problem. Ann Oper Res 2014:1–31. [3] Caprara A, Pisinger D, Toth P. Exact solution of the quadratic knapsack problem. INFORMS J Comput 1999;11(2):125–37. [4] Dudziński K, Walukiewicz S. Exact methods for the knapsack problem and its generalizations. Eur J Oper Res 1987;28(1):3–21. [5] García-Martínez C, Glover F, Rodríguez FJ, Lozano M, Martí R. Strategic oscillation for the quadratic multiple knapsack problem. Comput Optim Appl 2014;58(1):161–85. [6] García-Martínez C, Rodríguez FJ, Lozano M. Arbitrary function optimisation with metaheuristics. Soft Comput 2012;16(12):2115–33. [7] García-Martínez C, Rodríguez FJ, Lozano M. Tabu-enhanced iterated greedy algorithm: a case study in the quadratic multiple knapsack problem. Eur J Oper Res 2014;232(3):454–63. [8] Glover F, Kochenberger GA, editors. Handbook of metaheuristics. Berlin, Heidelberg: Springer; 2003. [9] Glover F, Hao JK. The case for strategic oscillation. Ann Oper Res 2011;183 (1):163–73. [10] Hanafi S, Freville A. An efficient tabu search approach for the 0-1 multidimensional knapsack problem. Eur J Oper Res 1998;106(2):659–75. [11] Higgins AJ. A dynamic tabu search for large-scale generalised assignment problems. Comput Oper Res 2001;28(10):1039–48. [12] Hiley A, Julstrom B. The quadratic multiple knapsack problem and three heuristic approaches to it. In: Proceedings of the genetic and evolutionary computation conference (GECCO); 2006. p. 547–552. [13] Igel C, Toussaint M. On classes of functions for which no free lunch results hold. arXiv:http://arXiv.org/abs/cs/0108011; 2001.

214

J. Qin et al. / Computers & Operations Research 66 (2016) 199–214

[14] Lapierre SD, Ruiz A, Soriano P. Balancing assembly lines with tabu search. Eur J Oper Res 2006;168(3):826–37. [15] Liu R, Xie X, Garaix T. Hybridization of tabu search with feasible and infeasible local searches for periodic home health care logistics. Omega 2014;47:17–32. [16] Mezura-Montes E, Coello CAC. Useful infeasible solutions in engineering optimization with evolutionary algorithms. In: MICAI 2005: advances in artificial intelligence; Berlin, Heidelberg: Springer; 2005. p. 652–62. [18] Pisinger D. An exact algorithm for large multiple knapsack problems. Eur J Oper Res 1999;114(3):528–41. [19] Saraç T, Sipahioglu A. A genetic algorithm for the quadratic multiple knapsack problem. Lect Notes Comput Sci 2007;4729:490–8. [20] Singh A, Baghel A. A new grouping genetic algorithm for the quadratic multiple knapsack problem. Lect Notes Comput Sci 2007;4446:210–8. [21] Sundar S, Singh A. A swarm intelligence approach to the quadratic multiple knapsack problem. Lect Notes Comput Sci 2010;6443:626–33.

[22] Talbi EG. A taxonomy of hybrid metaheuristics. J Heuristics 2002;8(5):541–64. [23] Wolpert DH, Macready WG. No free lunch theorems for optimization. IEEE Trans Evol Comput 1997;1(1):67–82. [24] Wu Q, Hao JK, Glover F. Multi-neighborhood tabu search for the maximum weight clique problem. Ann Oper Res 2012;196(1):611–34. [25] Wu Q, Hao JK. Memetic search for the max-bisection problem. Comput Oper Res 2013;40(1):166–79. [26] Lim A, Qin H, Xu Z. The freight allocation problem with lane cost balancing constraint. Eur J Oper Res 2012;217(1):26–35. [27] Qin H, Ming W, Zhang Z, Xie Y, Lim A. A tabu search algorithm for the multiperiod inspector scheduling problem. Comput Oper Res 2015;59:78–93. [28] Zhu W, Qin H, Lim A, Wang L. A two-stage tabu search algorithm with enhanced packing heuristics for the 3L-CVRP and M3L-CVRP. Comput Oper Res 2012;39(9):2178–95.