Two-stage solution-based tabu search for the multidemand multidimensional knapsack problem

The multidemand multidimensional knapsack problem (MDMKP) is a signiﬁcant generalization of the popular multidimensional knapsack problem with relevant applications. In this work we investigate for the ﬁrst time how solution-based tabu search can be used to solve this computationally challenging problem. For this purpose, we propose a two-stage search algorithm, where the ﬁrst stage aims to locate a promising hyperplane within the whole search space and the second stage tries to ﬁnd improved solutions by exploring the reduced subspace deﬁned by the hyper-plane. Computational experiments on 156 benchmark instances commonly used in the literature show that the proposed algorithm competes favorably with the state-of-the-art results. We analyze several key components of the algorithm to highlight their impacts on the performance of the algorithm.


Introduction
Given a set V = {1, 2, . . ., n} of n items, a set R = {r 1 , r 2 , . . ., r m } of m resources with a capacity upper limit b i for resource r i (1 ≤ i ≤ m), where each item j of V is associated with a profit c j and consumes a given quantity a ij for each resource r i (i ∈ {1, 2, . . ., m}), the popular NP-hard 0-1 multidimensional knapsack problem (MKP) involves selecting a subset of items from V such that the resource consummation of the selected items does not exceed the given capacity upper limit for each resource in R (knapsack constraints), while maximizing the total profit of the selected items.Formally, the MKP can be written as follows: x j ∈ {0, 1}, ∀j ∈ {1, 2, . . ., n} where c j ≥ 0, a ij > 0, b i > 0, ∀i ∈ {1, 2, . . ., m}, ∀j ∈ {1, 2, . . ., n} and Eq.
(3) indicates that the binary decision variable x j (1 ≤ j ≤ n) takes the value of 1 if the item j is selected, 0 otherwise.
Clearly, the classic MKP is a special case of the MDMKP when q equals 0 and the profit c j of item j takes a nonnegative value (i.e., c j ≥ 0, ∀j ∈ V ).
Like the MKP, the MDMKP has a number of practical applications [5] like obnoxious and semiobnoxious facility location [4,22,24], capital-budgeting, and portfolio-selection [3], among others.On the other hand, the MDMKP is computationally challenging, given that it generalizes the NP-hard multidimensional knapsack problem.Consequently, there is no polynomial-time algorithm for the MDMKP, unless P = NP.
Unlike the MKP that has been subject of intensive studies in the past decades (see e.g., [9,11,13,15,20,21,23,[25][26][27][28]), the MDMKP receives much less attention until now.Still, there exist several exact and heuristic approaches in the literature.For example, general mixed integer programming solvers like CPLEX can be used to solve instances with n ≤ 100 to optimality within an acceptable time.However, it is usually difficult for the existing exact approaches to find an optimal solution for larger instances.As a result, several heuristic algorithms have been proposed to solve large instances approximately.
Specifically, in 2005, Cappanera and Trubian presented a nested-tabu-search heuristic [5], which combines a standard attribute-based tabu search with an oscillation method presented in [13].In 2006, Arntzen et al. proposed an adaptive memory search method called Almha in [1], which uses a dynamical tabu search mechanism and a weighting scheme to handle infeasible solutions.Their computational results show the Almha algorithm outperforms the previous best MDMKP methods and can be viewed as one of the best performing MDMKP algorithms in the literature.In 2007, Hvattum and Løkketangen investigated the behavior of scatter search on the MDMKP [16].In 2009, Gortázar et al. introduced a black box scatter search method for general classes of binary optimization problems, and assessed their method on the MDMKP and some other binary problems [14].In particular, their method uses a static penalty approach proposed in [32] to handle the constraints of the MDMKP.In 2010, Hvattum et al. proposed an alternating control tree (ACT) search framework for the MDMKP [17], which can lead to an exact algorithm or heuristic algorithm by choosing the routine of solving subproblems.Their computational results show that the associated ACT algorithms have a high performance compared to a previous tabu search algorithm and scatter search algorithm.At the same year, Balachandar proposed a dominance principle based heuristic for the MDMKP [2].
In addition to these studies, there exist some theoretical investigations dedicated to the MDMKP in the literature.For example, Delissa investigated the existence and usefulness of equality cuts for the MDMKP [10], while Wishon and Villalobos studied robust efficiency measures for the MDMKP [30].
To enrich the solution arsenal for the MDMKP, we present in this work the first study of using solution-based tabu search [6,7,31] to effectively solve the MDMKP.Actually, unlike the popular attribute-based tabu search approach [12], solution-based tabu search began to attract attention only very recently.Interestingly, this approach already showed excellent performances on several binary optimization problems as reported in [19,20,29].This work aims thus to investigate the interest of the solution-based tabu search approach for the MDMKP.Compared to attribute-based tabu search, solution-based tabu search has at least two appealing features.First, this approach ensures a stronger intensification ability, which is crucial for locating good local optima.Second, this approach makes the notion of tabu tenure irrelevant, thus simplifying the design of the algorithm and reducing the number of required parameters.
We summarize the contributions of this work as follows.First, based on the solution-based tabu search approach, we introduce an effective two-stage search algorithm for the MDMKP.The first search stage aims to identify a promising hyperplane within the whole search space while the second search stage tries to find improved solutions by examining both feasible and infeasible solutions on the identified hyperplane.For both stages, the solution-based tabu search strategy is employed, which relies on a one-flip and swap neighborhoods and a hash-based mechanism to efficiently determine the tabu status of neighbor solutions.Second, we assess the performance of the proposed algorithm based on 96 benchmark instances commonly used in the literature (n = 100, 250, m = 5, 10 and q = 2, 5, 10).The computational results show that the algorithm improves and matches the best known solutions for 17 and 71 instances respectively.Moreover, we report detailed computational results of the proposed algorithm for 60 additional instances with a large number of constraints (with n = 100, 500, m = q = 30).Third, given that the ideas of the two-stage search framework and solution-based tabu search developed in this work are quite general, they could be applied to solve other related binary optimization problems.
The remaining parts of the paper are structured as follows.In the next section, the proposed two-stage tabu search algorithm is described.In Section 3, we present the computational assessment of the proposed algorithm and report experimental results on the well-known benchmark instances.In Section 4, several essential ingredients of the algorithm are investigated to show how they affect the performance of the algorithm.In the last section, we summarize the present work and provide research perspectives.
Our two-stage solution-based tabu search (TSTS) algorithm combines two search procedures working on two different search spaces.The first stage of the algorithm performs an exploratory search within the whole search space to find a feasible solution as good as possible.Starting from this solution, the second stage carries out a focused exploitation within the reduced space composed of candidate solutions with exactly k selected items (k being identified by the final solution of the first search stage).To explore both search spaces, TSTS relies on two solution-based tabu search procedures guided by a penaltybased evaluation function.One notices that the two-stage search strategy has been used with success to solve other knapsack problems like the Quadratic Knapsack Problem [8] and the classic MKP [26].Specifically, the TSTS algorithm first generates randomly an initial solution by its initialization procedure (Section 2.2).Then the algorithm enters the first search stage where the initial solution is improved by the solution-based tabu search procedure presented in Sections 2.3 (line 4).During this search stage, the algorithm explores both feasible and infeasible solutions within the whole search space to find a high-quality feasible solution.At the end of its search stage, the best (feasible) solution found and the consumed time t are returned.At this point, the second search stage is triggered, which starts from the solution returned by the first stage and uses another solution-based tabu search procedure (Section 2.4) to seek improved solutions (line 5).During this stage, the search is limited to the hyperplane Ω [k] (k being the number of the selected items in the returned solution of the first stage, see Section 2.4.1).Finally, the whole algorithm terminates when a given time limit t max is met, and the best solution found during the search process is returned as the final result of the algorithm (line 6).

Initial Solution
Algorithm 2: Procedure of Generating Initial Solution Function InitialSolution() Input: An instance I, the size of instance (n) Output: An initial solution s = (x 1 , x 2 , . . ., x n ) begin for i ← 1 to n do The initial solution of the TSTS algorithm is generated by a randomized procedure whose pseudo-code is given in Algorithm 2. Specifically, given an instance with n items, the initialization procedure assigns randomly a value from the set {0, 1} to each component x i (i = 1, 2, . . ., n) to obtain an initial solution s = (x 1 , x 2 , . . ., x n ).This random initialization has the advantage of being simple and fast.However, an initial solution generated in this way can be infeasible.If this is the case, its feasibility will be established during the first search stage described below.

Tabu Search Method of the First Search Stage
The first search stage is ensured by a solution-based tabu search algorithm (denoted by T abuSearch 1 ()), whose pseudo-code is given in Algorithm 3.After initiating its tabu lists (lines 3-5), T abuSearch 1 () performs a number of iterations to improve the current solution until 1) a feasible solution is found and no improvement can be observed during α consecutive iterations, where α is a parameter called the depth of tabu search, or 2) the allowed maximum time limit t max is reached (lines 8-17 The ingredients of this tabu search algorithm, including the search space, the neighborhood structures and the tabu strategy, are respectively described in the following subsections.

Search Space and Neighborhood
The search space Ω explored by the T abuSearch 1 () procedure is composed of all feasible and infeasible solutions of the given problem instance, i.e., The neighborhood used by T abuSearch 1 () is a combined neighborhood composed of two basic neighborhoods, namely the one-flip neighborhood N 1 and the swap neighborhood N 2 .The one-flip neighborhood N 1 is defined by the one-flip operator (denoted by F lip(•)).Given a solution s = (x 1 , x 2 , . . ., x n ), a one-flip move F lip(q) consists of changing the value of a variable x q to its complementary value 1 − x q .As such, the neighborhood N 1 (s) of solution s includes all possible solutions that can be obtained by applying the one-flip operator to s. Formally, the N 1 (s) can be written as follows: The neighborhood N 2 is defined by the swap operator (denoted by Swap(•, •)).Given a solution s = (x 1 , x 2 , . . ., x n ), let I 1 = {q : x q = 1 in s} and I 0 = {q : x q = 0 in s}, the swap neighborhood N 2 (s) can be written as follows: The first tabu search algorithm explores the union of these two neighborhoods, i.e., N (s) = N 1 (s)∪N 2 (s), whose size equals to n+|I 1 |×|I 0 |.At each iteration of the algorithm, a best non-tabu solution from N (s) according to the extended evaluation function defined by Eq. ( 17) in Section 2.3.3 and the tabu strategy explained in Section 2.3.2 is selected to replace the current solution s.

Tabu Strategy
In the present tabu search method, we adopt the solution-based tabu strategy to determine the tabu status of neighbor solutions.Specifically, the tabu lists are based on three hash vectors H 1 , H 2 , and H 3 of length of L, where each position of them represents a binary variable, and each hash vector H t is associated with a hash function h t .In particular, the effect of hash functions is to map a candidate solution of the search space Ω to an index of H t , i.e., h t : Ω → {0, 1, 2, . . ., L − 1}.
Based on these hash vectors and the associated hash functions, we determine the tabu status of candidate solutions by the following rule.Given a candidate solution s, the hash vectors H t (t = 1, 2, 3) and the associated hash functions h t , s is identified as a tabu solution if Otherwise, s is determined as a non-tabu solution.
Let s = (x 1 , x 2 , . . ., x n ) be a candidate solution, our hash functions h t (t = 1, 2, 3) are given by: where γ t is a parameter that is used to define each hash function and L is the length of hash vectors that is empirically set to 10 8 in this work.
Given a solution s and its hash value h(s), the hash value of its neighbor solutions can be determined in O(1) according to Eq. ( 15).Thus, the time complexity of determining the tabu status of a neighbor solution is O(1).

Extended Evaluation Function
Since the search space explored by the first tabu search algorithm contains both feasible and infeasible solutions, we devise an extended evaluation function F which uses a penalty function P to assess constraint violations.
Let s be a candidate solution in Ω, the penalty value P (s) is defined as the summation of all constraint violations, i.e., Thus, a small (large) function value P (s) means a weak (strong) constraint violation in s.In particular, P (s) = 0 means that s is a feasible solution.
Given this penalty function, the extended evaluation function F (s) is defined as a linear combination of the objective function f (s) in Eq.( 4) and P (s): where λ is a weighting factor that is empirically set to 10 2 in this work.
For any two solutions s ′ and s ′′ in As shown in our experimental results (Section 3, Section 4.2 and the Appendix), the first search stage equipped with the extended evaluation function always ends up with a feasible solution for all tested instances, including those with 30 demand constraints and 30 knapsack constraints.In other words, the first solution-based tabu search algorithm is typically able to locate a promising valid hyperplane that is further explored by the second search stage.

Tabu Search Method of the Second Search Stage
The second optimization stage of the TSTS algorithm uses another tabu search algorithm (denoted by T abuSearch 2 (), Algorithm 4) to examine candidate solutions of a given hyperplane Ω [k] (see below).Unlike the first tabu search algorithm, T abuSearch 2 () explores only solutions that contain exactly k selected items.T abuSearch 2 () first initializes the hash vectors (lines 3-5), and then performs a number of iterations to improve the current solution (lines 7-14).At each iteration, the algorithm replaces the current solution by a best non-tabu neighbor solution s ′ in terms of the evaluation function in Eq. ( 17).
During the search, the best feasible solution encountered s * is updated each time a better feasible solution is found, and the hash vectors are accordingly updated by the new solution (line 13).Finally, the algorithm terminates if the time limit t max is reached, and then returns the best feasible solution s * found during the search process.

Search Space, Tabu Strategy and Evaluation Function
The search space Ω [k] explored by T abuSearch 2 () is composed of both feasible and infeasible solutions with a fixed number of k selected items.In other words, Additionally, like the first tabu search algorithm, T abuSearch 2 () uses the solution-based tabu strategy described in Section 2.3.2 to determine the tabu status of neighbor solutions, and employs the extended evaluation function in Eq. (17)

Neighborhood Structure
To search effectively the hyperplane Ω [k] , T abuSearch 2 () uses a constrained swap neighborhood N 3 (s).Formally, given a solution s = (x 1 , x 2 , . . ., x n ), the neighborhood N 3 (s) is given by: where f (•) is the objective function of the MDMKP in Eq. ( 4), s * is the best feasible solution found so far in the current tabu search run, I 1 and I 0 denote respectively the sets of indices having the value of 1 (selected items) and 0 (non selected items) in s.Clearly, the size of this neighborhood is bounded by It is worth noting that we constraint the neighborhood by f (s ⊕ Swap(i, j)) > f (s * ) to eliminate non-promising neighbor solutions.Similar idea was previously investigated for the related MKP in [26].

Space and Time Complexities of the Algorithm
At each stage of our TSTS algorithm, in addition to three hash vectors (i.e., H 1 , H 2 and H 3 ) with a length of L, we maintain three solutions (i.e., s, s ′ , and s * ) to follow the search process, where each solution is stored by two vectors (i.e., I1 and I 0 ) with a maximum length of n and a vector W = (w 1 , w 2 , . . ., w m+q ) where w i = n j=1 a ij x j (j ∈ V ) holds.Thus, the space complexity of our TSTS algorithm is bounded by O(n + m + q + L).
In addition, for each neighbor solution in the search space, the time complexity of evaluating its quality is bounded by O(m + q), since there are m + q constraints needed to be checked.Hence, for each iteration, the time complexities of the first and second tabu search stages are respectively bounded by according to the size of neighborhoods explored by the algorithm (see Sections 2.3.1 and 2.4.2).

Experimental Results and Comparisons
We evaluate the proposed TSTS algorithm by conducting extensive computational experiments based on four sets of benchmark instances commonly used in the literature and by making a comparison between our results and the state-of-the-art results in the literature.

Benchmark Instances
In this study, we employed four sets of benchmark instances to assess the performance of our TSTS algorithm, where the first two sets of benchmark instances are available at http://www.optsicom.es/binaryss,and the third and fourth sets of benchmark instances are available in OR-Library 1 .The first set contains 48 small instances with n = 100, and the second set contains 48 larger instances with n = 250.In addition, for the instances in the first two sets, the number of knapsack constrains m varies from 5 to 10, and the number of demand constraints q belongs to {2, 5, 10}.The third set includes 30 instances with n = 100, where both the number of knapsack constrains m and the number of demand constraints q equal to 30.The fourth set is composed of 30 large instances with n = 500, m = 30 and q = 30.Our TSTS algorithm employs four parameters, including γ 1 , γ 2 and γ 3 that are used to define the hash functions, and the depth α of tabu search used in the first search stage.The parameters γ 1 , γ 2 and γ 3 are set as in Table 1 according to the analysis shown in Section 4.3 while α is empirically set to 5 × 103 .

Parameter Settings and Experimental Protocol
In addition, our algorithm was implemented in C++ and compiled by g++ compiler with -O3 option2 .All computational experiments were carried out on a computer with an Intel E5-2670 processor (2.5 GHz and 2G RAM), running the Linux operating system.Moreover, when running the DIMACS machine benchmark procedure dmclique 3 , our processor requires 0.19, 1.17, and 4.54 seconds to solve the graphs r300.5, r400.5, and r500.5, respectively.Finally, due to the stochastic nature of the algorithm, we independently ran the algorithm 30 times to solve each instance, where the time limit t max for each run was set to 60 seconds for instances with n ≤ 250 according to [14].
For the large instances with n = 500, the time limit t max was set to n seconds, where n is the size of instances.

Computational Results and Comparison
In this section, we report the computational results of our TSTS algorithm on the first two sets of benchmark instances.We provide in the Appendix the computational results of TSTS on the third and fourth sets of benchmark instances for which no detailed results are available for existing algorithms in the literature.
The computational results of our TSTS algorithm on the first set of benchmark instances with n = 100 are summarized in Table 2, together with the results of the Almha algorithm implemented in [14].Column 1 gives the names of instances.Columns 2 and 3 indicate respectively the best known objective  values (BKV ) and the results of the Almha algorithm, which were reported in [14] and available at http://www.optsicom.es/binaryss.The results of our TSTS algorithm are reported in columns 4-7, including the best objective value obtained over 30 runs (f best ), the average objective value (f avg ), the worst objective value (f worst ), and the standard deviation (σ f ) of objective values.The row Avg.shows the average result over all instances tested for each column.The rows #Better, #Equal, and #W orse respectively show the number of instances for which our best result f best is better than, equal to, worse than the BKV.Moreover, in terms of f best , our improved results (new lower bounds) are indicated in bold and our worse results are indicated in italic compared to the BKV .Finally, to verify whether there exists a significant difference between our best results (f best ) and the BKV, we provide in the last row the p-values (∈ [0, 1]) from the non-parametric Friedman test, where a p-value less than 0.05 means a significant difference between the compared results.This test was also performed to compare the results of the Almha algorithm and our average results (f avg ) 4 .
Table 2 shows that our TSTS algorithm matches the best known results for 47 out of 48 instances, and improves the best known result for the remaining one instance, leading to an improved Avg. value compared to the averaged BKV (18108.96vs. 18108.10).When comparing the average objective values f avg of our TSTS algorithm with the results of the Almba algorithm, one can find that both algorithms have a very similar performance, which is confirmed by a large p-value of 0.87.In addition, regarding the average results over all instances (Avg.), the value of f avg of our TSTS algorithm is 18088.27,which is slightly better than that of Almba (i.e., 18083.17).These outcomes show that our TSTS algorithm performs similarly on these small instances n = 100 compared to the state-of-the-art Almba algorithm.
The second experiment aims to evaluate our TSTS algorithm on the set of 48 larger instances with n = 250, and the computational results are summarized in Table 3, where we report the same statistics as in Table 2.We observe from Table 3 that our algorithm matches the best known results for 18 out of 48 instances, improves the best known results for 12 instances, and misses the best known results for the remaining 18 instances.In terms of Avg., the value of f best of TSTS is slightly worse than the value of BKV (i.e., 49821.10 vs. 49836.48),and slightly better than the result of Almha (i.e., 49821.10 vs. 49717.81),which is however slightly better than the value of f avg of TSTS (i.e., 49687.60 vs. 49717.81).This experiment indicates that under the short time limit t max = 60 seconds, TSTS performs globally well on these larger instance with n = 250 especially by finding 12 improved best solutions.Additionally, we observe from Table 3 that TSTS performs particularly well for instances with a large number of constraints and achieves a better result than Almha for most of the 12 instances with 10 demand constraints and 10 knapsack constraints.Finally, we mention that the results of TSTS can be further improved by increasing the time limit (see the detailed results in the Appendix), implying that the current time limit (t max = 60) is too short for the TSTS algorithm on these instances.

Analysis and Discussions
To shed light on the functioning of the proposed algorithm, we now analyze and discuss several essential components of the TSTS algorithm.

Effectiveness and Robustness Analysis for Two-Stage Strategy
To study the effectiveness of the underlying two-stage search strategy, we summarize in ).
The results of small instances with n = 100 are reported in the first 4 columns of Table 4, including the name of instances, the k value of the best solution obtained over 30 runs (k best ), the average k value of solutions obtained, and the standard deviation of k values obtained (σ k ).The results of larger instances with n = 250 are reported in columns 5-8, with the same statistics as in columns 1-4.In addition, the row Avg.shows the average results of standard deviations σ k of k values over all tested instances of each test set.
Table 4 shows that the value of k avg is very close to that of k best for most tested instances, which means that the first search stage of the TSTS algorithm is able to find a hyperplane that is very close to the best hyperplane containing the current best known solution.On the other hand, we observe that the standard deviations σ k of k values obtained are very small for most instances.
In particular, the average standard deviations of k values are respectively 0.22 and 0.66 for the set of small instances with n = 100 and the set of larger instances with n = 250.Hence, this experiment confirms to some extend the effectiveness and robustness of the two-stage search strategy employed by the TSTS algorithm.

Effects of Two Stages on the Performance of Algorithm
Table 5 Comparison between the two search stages on the instances with n = 250.To investigate the respective role of the two stages of our algorithm, we carried out an experiment based on the instances with n = 250.We ran our TSTS algorithm 30 times to solve each instance according to the experimental protocol of Section 3.2.The average results from the first stage and the second stage over 30 independent runs are summarized in Table 5.The first column gives the name of instances, columns 2 and 3 show the objective value (f 1 ) obtained by the first stage and computation time (t 1 ) in seconds needed to reach f 1 .Columns 4 and 5 show the objective value (f 2 ) obtained by the second stage and the computation time (t 2 ) needed to reach f 2 from f 1 .The last two columns indicate the gap between f 2 and f 1 and the improvement ratio (ρ) of f 2 over f 1 , which is calculated as ρ = 100 × (f 2 − f 1 )/f 1 .
We observe from Table 5 that the first search stage of the TSTS algorithm is able to obtain a high-quality feasible solution for each instance and the solutions obtained in the first stage can be further improved during the second search stage.Furthermore, the improvements of f 2 over f 1 are significant with an average improvement ratio ρ close to 1%.On the other hand, regarding the computation times needed by the two search stages, we observe that most computational efforts are required by the second search stage and that t 2 is about five times larger than t 1 .Of course, this proportion depends also on the setting of the parameters α and t max .These outcomes indicate that both search stages of the TSTS algorithm are indispensable for the high performance of the algorithm.The first search stage is able to generate high-quality feasible solutions while the second search stage is able to further improve the solutions by performing an intensified search in the reached hyperplane.

Sensitivity Analysis of Hash functions
Now, we investigate the impacts of hash functions on the performance of the algorithm and discuss the sensitivity of the associated parameters.For this purpose, we carried out an additional experiment based on 30 representative instances in terms of the numbers of knapsack and demand constraints.We ran our TSTS algorithm 30 times for each instance and each parameter combination of (γ 1 , γ 2 , γ 3 ), where γ i (i = 1, 2, 3) are the parameters used to define the hash functions h i (see Section  number of instances for which the associated setting of parameters yielded the best results and the average results over all instances tested. We observe from Table 6 that the algorithm is sensitive to the setting of the parameters γ 1 ,γ 2 , and γ 3 .For the parameter combinations containing two small parameter values, the TSTS algorithm performed badly.For example, the algorithm with the combination (1.1, 1.5, 1.8) yielded the worst results in terms of Avg.However, when all parameters in (γ 1 , γ 2 , γ 3 ) take a relatively large value, the algorithm obtained a much better performance.Taking (1.9, 2.1, 2.3) as an example, the algorithm achieved the best results on 9 instances.In summary, the parameters γ 1 ,γ 2 , and γ 3 have an important impact on the performance of the algorithm, and parameter combinations containing at least two large parameter values lead usually to a good performance of the algorithm.

Effectiveness of Solution-based Tabu Strategy
The solution-based tabu strategy is an essential ingredient of our TSTS algorithm.To show its effectiveness with respect to the popular attribute-based tabu strategy, we created a variant A-TSTS of the TSTS algorithm by replacing the solution-based tabu strategy with the popular attribute-based tabu strategy, while keeping the other TSTS components unchanged.In A-TSTS, the adopted tabu strategy can be simply described as follows.Given an incumbent solution s = (x 1 , x 2 , . . ., x n ), once a 0-1 variable x i (1 ≤ i ≤ n) is flipped as x i ← 1 − x i , x i is forbidden to be flipped again for the next tt iterations (tt is the tabu tenure) and the associated neighbor solutions are excluded for consideration during the period identified by the tabu tenure.We empirically set tt = C + rand[0, 2] where C is a parameter that takes the value of 20 and rand[0, 2] is a random integer in [0, 2].Finally, the tabu status of a variable is disabled if flipping the variable leads to a solution better than all previously visited solutions (this is the so-called aspiration criterion).
To compare TSTS and A-TSTS, we carried out an experiment based on the set of 48 large instances with n = 250, where both algorithms were run 30 times according to the experimental protocol given in Section 3.2.The computational results are summarized in Table 7 where we show for each algorithm the best results (f best ) obtained over 30 runs, the average results (f avg ), and the worst results obtained (f worst ).In addition, the row 'Avg.' shows the average results for each performance indicator.Finally, to check whether there exists a significant difference between the two algorithms in terms of f best , f avg and f worst , we provide the p-values from the non-parametric Friedman test in the last row, where a p-value smaller than 0.05 implies a significant difference between the compared results.We observe from Table 7 that the solution-based algorithm TSTS dominates the attribute-based algorithm A-TSTS in terms of all indicators, by reporting better results in terms of f best , f avg and f worst on all the instances.Furthermore, the small p-values indicate that the performance differences between the compared results are statistically significant.Therefore, this experiment confirms that under the two-stage framework of this work, the solution-based tabu strategy is much more suitable than the attribute-based tabu strategy for solving the MDMKP.

Discussion about the Solution-based and Attribute-based Tabu Search Approaches
Attribute-based tabu search is a popular approach, whose key idea is to prevent one attribute or a combination of several attributes of a solution from being changed during a number of iterations since their last changes.For such a method, a tabu attribute or a combination of several attributes is usually associated with a number of tabu solutions.However, unlike attribute-based tabu search, solution-based tabu search tries to record all the visited solutions and prevent them from being revisited during the following search process.Thus, solution-based tabu search ensures a stronger intensification search ability.
Recent studies on several binary optimization cases including two dispersion problems (minimum difference dispersion [29] and maximum min-sum dispersion [19]) and the classic multidimensional knapsack problem [20] demonstrate that solution-based tabu search is more suitable than attribute-based tabu search.On the other hand, our experience on another dispersion problem (max-mean dispersion [18]) does not confirm the advantage of the solutionbased tabu search approach over the attribute-based tabu search approach.So one interesting question concerning the solution-based and attribute-based tabu search approaches is under what circumstances one approach will be more suitable than the other.Given that solution-based tabu search has been investigated only very recently, little knowledge is currently available, which makes it difficult to provide a meaningful guidance on the choice between these two approaches.To fully characterize these approaches and understand the relationships between these approaches and the optimization problem under consideration, more studies are clearly needed, which constitutes an interesting research perspective.
Finally, the ideas of the present solution-based tabu search algorithm being quite general, they could conveniently be tested on other binary optimization problems, by adjusting the γ parameters of the hash functions (Eq.( 15), Section 2.3.2) or by increasing the number of the hash vectors and the associated hash functions.
In this work, we investigated the NP-hard multidemand multidimensional knapsack problem, by proposing a two-stage tabu search (TSTS) algorithm that combines two solution-based tabu search procedures and a penalty-based evaluation function to explore different search spaces.Computational results on 156 benchmark instances showed that the proposed algorithm is competitive compared to the state-of-art results in the literature, especially for those instances with a large number of knapsack and demand constraints.
The experimental analysis showed the usefulness of the two-stage search strategy and the respective impacts of two stages on the performance of the algorithm.The first search stage is able to reach a promising hyperplane containing high-quality solutions, and the second search stage is able to find elite solutions by an intensified examination of the given hyperplane.We also showed that the hash functions used by the tabu search algorithms are a key component that significantly influences the performance of the algorithm.
This work enriches the existing tools for effectively solving the MDMKP and invites more research and attention on solution-based tabu search for solving binary optimization problems in the future.Specifically, given that the ideas of the two-stage strategy and the solution-based tabu search procedures developed in this work are quite general, it would be interesting to investigate their effectiveness on other problems, especially those related to subset selection problems for which the search space can be divided into a series of hyperplanes.It is also interesting to study search strategies mixing solution-based and attribute-based approaches.As a more fundamental research perspective, studies on a characterization of both solution-based and attribute-based search approaches are needed, which could ease the choice of one or the other approach to solve additional binary optimization problems.

A Appendix
This Appendix presents the detailed computational results of the proposed TSTS algorithm for the third and fourth sets of benchmark instances for which no detailed results are available in the literature.These instances have 100 or 500 items, 30 knapsack constraints and 30 demand constraints, making them more difficult to solve.For each of these instances, the TSTS algorithm was run 30 times under the condition presented in Section 3.2 and the computational results are summarized in Table A.1, where the symbols have the same meanings as in the previous tables.In addition, we report in Table A.2 the computational results of our TSTS algorithm on instances with n = 250 under a long time limit of t max = 300 (instead of t max = 60 used in Section 3).The results reported in this Appendix can serve as references for future comparative studies of new MDMKP algorithms.
We observe from Table A.1 that TSTS is able to reach the best result (f best ) with a success rate of 100% for 15 out of 30 small instances with n = 100, which means a good robustness of our TSTS algorithm on these instances.Note that for some small instances with n = 100, m = 30, and q = 30, it is difficult for some state-of-the-art algorithms in the literature [17] to obtain a feasible solution.However, it is very easy for our TSTS algorithm to obtain a feasible solution for all these instances.For other instances, the standard deviation (σ f ) of objective values obtained by the TSTS algorithm is relatively small, which indicates a good robustness of the algorithm.Regarding the value of k which represents the number of items selected, it can be seen that the gap between k best and k avg and the standard deviation (σ k ) of k values obtained are very small, implying that the two-stage strategy of the algorithm is very robust and effective.
The proposed TSTS algorithm is thus composed of two optimization stages, where the first stage identifies a suitable hyperplane Ω [k] (see Section 2.4.2 for the definition of hyperplane) that is exploited intensively during the second optimization stage to locate improved solutions (see Algorithm 1).
ProcedureAlgorithm 1: General procedure of two-stage tabu search algorithm for the MDMKP Function TSTS() Input: Instance I, time limit t max Output: The best solution s * Initial solution s, extended evaluation function F , hash vectors H 1 , H 2 , H 3 of length L, hash functions h 1 , h 2 , h 3 , depth of tabu search α, time limit t max Output: The best solution s * found, the time t consumed by the search The tabu search method used in the second search stage Function T abuSearch 2 () Input: Initial solution s, extended evaluation function F , penalty function P , hash vectors H 1 , H 2 , H 3 of length L, hash functions h 1 , h 2 , h 3 , time limit t max , and time (t) consumed in the first search stage.
to evaluate the solutions in the search space Ω [k] .Algorithm 4: Find in N 3 (s) a best non-tabu solution s ′ in terms of the extended evaluation function F in Eq. (17)/* N 3 (s) is defined in Eq. (18), and the tabu rule is given in Section 2.4.1.

Table 1
Settings of parameters

Table 2
Computational results and comparison on the 48 instances with n = 100.In terms of f best , the improved results are indicated in bold compared to the best known objective values (BKV ).

Table 3
Computational results and comparison on the 48 instances with n = 250.In terms of f best , the improved results are indicated in bold and the worse results are indicated in italic compared to the best known objective value (BKV ).

Table 4
the computational results about the k values returned by the TSTS algorithm, where the results are based on the experiments in Section 3.3 and the value of k represents the number of selected items in the solution found.It is worth noting that the purpose of the first search stage is just to discover a promising hyperplane Ω [k] that contains high quality solutions, and the second search stage aims to locate improved solutions in the given hyperplane.Hence, the two-stage search strategy can be considered to be relevant and robust if the first search stage is able to reach very stably the identical or close hyperplane to the best known solution (i.e., Ω [k best]

Table 4
Statistical results over 30 runs in terms of the number k of items in the obtained solutions.
The experimental results are summarized in Table6, where the first column gives the name of instances, the second row shows the settings of parameters, and the average objective values (f avg ) obtained over 30 runs are reported in columns 2-11 for each parameter combination and each instance, respectively.In addition, the rows #Best and Avg. of the table indicate respectively the

Table 6 .
Influence of the hash functions on the performance of algorithm.Each instance was independently solved 30 times for each parameter combination, and the average objective values (f avg ) over 30 runs are reported.

Table 7
Comparison between the solution-based and attribute-based tabu strategies on the instances with n = 250.The better results between the two algorithms are indicated in bold in terms of f best , f avg and f worst .

Table A .
2 shows that our TSTS algorithm improves the best known results for 16 out of 48 instances, matches the best known results for 24 instances, and misses the best known results for only 8 instances.These results indicate Table A.1 Computational results of the TSTS algorithm on the instances with a large number of constraints under the condition presented in Section 3.2.For each of these instances, there are 30 knapsack constraints and 30 demand constraints.