European Journal of Operational Research 132 (2001) 22-38 www.elsevier.com/locate/dsw Theory and Methodology A Tabu search heuristic for the generalized assignment problem Juan A. Diaz !, Elena Fernandez * Dpt. d'Estadistica i Investigaciö Operativa, Universität Politecnica de Catalunya, Pau Gargallo 5, 08028 Barcelona, Spain Received 18 August 1998; accepted 20 March 2000 Abstract This paper considers the generalized assignment problem (GAP). It is a well-known NP-hard combinatorial optimization problem that is interesting in itself and also appears as a subproblem in other problems of practical importance. A Tabu search heuristic for the GAP is proposed. The algorithm uses recent and medium-term memory to dynamically adjust the weight of the penalty incurred for violating feasibility. The most distinctive features of the proposed algorithm are its simplicity and its flexibility. These two characteristics result in an algorithm that, compared to other well-known heuristic procedures, provides good quality solutions in competitive computational times. Computational experiments have been performed in order to evaluate the behavior of the proposed algorithm. © 2001 Elsevier Science B.V. All rights reserved. Keywords: Generalized assignment problem; Metaheuristics; Tabu search; Local search 1. Introduction The generalized assignment problem (GAP) is a well-known NP-hard combinatorial optimization problem. It considers a situation in which n jobs have to be processed by m agents. The agents have a capacity expressed in terms of a resource which is consumed by job processing. The decision maker seeks the minimum cost assignment of jobs to * Corresponding author. Tel.: +34-93-4016948; fax: +34-934015855. E-mail addresses: jadiaz@eio.upc.es (J.A. Diaz), ele-na@eio.upc.es (E. Fernandez). 1 Partially supported by grant 110598/110666 from Consejo Nacional de Ciencia y Tecnologia (CONACYT), Mexico. agents such that each job is assigned to exactly one agent subject to the agents' available capacity. Assignment of jobs to computers in a computer network, storage space allocation, design of communication networks with node capacity constraints, scheduling variable-length television commercials into time slots, etc., are examples of practical applications for the GAP. It also appears as a subproblem in many real-life problems such as vehicle routing, plant location, flexible manufacturing systems and resource scheduling. In this paper, we propose a Tabu search algorithm for the GAP. Two distinctive features of this algorithm have to be remarked: its simplicity and flexibility. These two characteristics result in an algorithm that provides good quality solutions 0377-2217/01/$ - see front matter © 2001 Elsevier Science B.V. All rights reserved. PII: S0377-2217(00)00108-9 J.A. Diaz, E. Fernandez I European Journal of Operational Research 132 (2001) 22-38 23 that are competitive with other well-known heuristic procedures. The flexibility is obtained from two sources: a strategic oscillation scheme that results from allowing the search to explore the infeasible solutions' space and an adaptively modified objective function. During the search, the solutions' space is enlarged to permit considering infeasible solutions. This gives the search a higher level of lexibility and some feasible solutions, that otherwise would not have been considered, can be reached by successively crossing the feasible/infeasible boundary. This strategy has already been used for the GAP by other authors (see [16,25-27]). As explained later on, the difference is that in our approach no solution is qualitatively preferred to others in terms of its structure so a priori all considered solutions are equally desired. This further increases the degree of lexibility of the procedure. Our algorithm considers a modified objective function that includes a penalty term associated with infeasible solutions. The weight given to such penalty is dynamically adapted according to the history of the search. A similar idea has been proposed independently by Yagiura et al. [25], where the search parameter is also adjusted automatically according to the result of the previous iteration. However, in our approach the modified objective function plays a central role since it is the element that really takes care of guiding the search. This is achieved by using a combination of recency-based and medium-term memory for adjusting the penalty weight. The combination of these two memory components results in an objective function that is continuously adapting itself to the stage of the search. This is an issue that we believe is interesting since, traditionally memory has been widely used to explore complex neighborhood structures, whereas it has been used in a very limited fashion to obtain automated mechanisms to modify the objective function that guides the search process. The adaptability of the proposed algorithm permits not having to resort to complex neighborhood structures. This, in its turn, has an important efect on the overall speed ofthe procedure. This is particularly true with the neighborhood structures that are explored (the classical shift and swap neighborhoods) and with the long-term memory component considered which is the same (the assignments' frequency) both for the intensi-ication and diversiication phases. The result is an algorithm where the search strategy is based on (a) exploring large areas of the infeasible solutions' space and (b) dynamically adapting the objective function to the stage of the search. The proposed algorithm has been tested on diferent sets of problems publicly available in Beasley OR-library (http:IImscmga.ms.ic.ukIjebl orlibIgapinfo.html) and in Yagiura library (http:II www-or.amp.i.kyoto-u.ac.jplyagiuralgap). It has also been compared to other well-known existing methods proposed by other authors in [7,16,21,25,27]. The results prove that the performance of the algorithm is remarkable since it produces high-quality solutions in small computational times. This paper is organized as follows: in Section 2 a model for the GAP is given and a relaxation RGAP that will be used throughout the paper is presented. Section 3, reviews solution methods for the GAP found in the literature. The Tabu search heuristic for the GAP is explained in detail in Section 4: neighborhood structures, the neighbor generation mechanism, and the memory functions used are detailed in that section. Section 4 also presents the mechanism for dynamically adjusting the penalty parameter. In Section 5, computational results are presented and our approach is compared to other existing algorithms. Finally, Section 6 includes some conclusions and remarks. 2. The model Let / = {1,m} be the set of agents and J = {1,..., n} the set of jobs. A standard integer programming formulation for the GAP is the following: (GAP) min z =^ J2C,JX,J> s.t. J.A. Diaz, E. Fernandez I European Journal of Operational Research 132 (2001) 22-38 24 5>,7 = 1 V/ e J, (2) ie/ ^a,7x,/ < bi Vi e /, (3) jeJ e{0,1}Vi e / V/ e J, (4) where cj is the assignment cost of job j to agent i, ai/ the resource required for processing job j by agent i,and bi is the available capacity of agent i. Decision variables xi/ are set to 1 if job j is assigned to agent i, 0 otherwise. Constraints (2), together with the integrality conditions on the variables, state that each job is assigned to exactly one agent. Constraints (3) insure that the resource availabilities of agents are not exceeded. In our approach, the assignment costs have been expressed in terms of their best assignment Aij = - < min, where c min/ = min{ciy : i e /} is the cost of the best assignment for job j. AiJ can be interpreted as the relative assignment cost of job j to agent i with respect to its best assignment cost. The use of the /dy-'s allows fair comparisons for deciding which jobs have to be assigned to a given agent, particularly when the assignment costs have different ranges for the various jobs. Then, the objective function can be rewritten as min z = K + min ^ dyX/ >, I ie/ jeJ J where K = J2jejC min/, represents a fixed assignment cost which has to be incurred. Since K is constant we can remove it from (1) and consider the following equivalent objective function: ie/ /eJ Throughout the paper, we consider a relaxed formulation for the GAP in which the capacity constraints (3) are dropped and a penalty term is added to the objective function: (RGAP) min z'' = ]T ^ '-x.- + p, (1'') ie/ /eJ Xi/ = 1 V/ e J, (2) ie/ Xjj e{0,1}VJ e J, Vi e /. (4) The term P in the objective function evaluates the infeasibility of the solutions to RGAP, with respect to constraints (3) of GAP. It takes the form ie/ where st = max^e j «,.,f.. — bi, 0} measures the amount by which the capacity ofagent i is violated by a given solution. The parameter p is a penalty factor. The reasons for considering RGAP are next explained. In the formulation of the GAP constraints (3) are 'difficult' and, usually, they are not easy to fulfill. For this reason these constraints may conine feasible solutions to a fairly narrow region. When exploring the neighbors of a given solution in a heuristic based method, restricting the search to feasible solutions using simple moves tends to generate very few trial solutions, forcing the process to terminate very quickly, often with low quality solutions. On the other hand, allowing some violation of feasibility has two advantages: irst, it permits the execution of moves that are less complex than might otherwise be required. Second, it gives the search a higher level of flexibility and much larger areas of the solutions' space can be explored. The deviation from feasibility can be measured and controlled through a penalty term in the objective function. Violation of feasibility is allowed in several recent works for the GAP like the Tabu search approach considered in [16], the VDS proposed by Yagiura et al. [27], the VDS with branching search considered in [26] and the ejection chain approach proposed in [25]. However, these references, each one in its own context, follow the same a priori criterion that, in some sense, guides the search that is performed. Namely: among feasible solutions to RGAP, feasible solutions to GAP (those that fulill constraints (3)) are always preferred. This is an important diference with respect to our approach. We consider all feasible solutions to RGAP equally attractive in J.A. Diaz, E. Fernandez I European Journal of Operational Research 132 (2001) 22-38 25 terms of their structure. This is further explained in Section 4.2. 3. Related work In this section, we present a brief review of solution methods found in literature for GAP. A comprehensive treatment of these methods can be found in [6] and in [21]. Exact methods for the GAP differ from each other in the way the lower bounds are computed. These bounds are derived from relaxations of the assignment constraints (constraints (2)), the knapsack constraints (constraints (3)) or the integrality conditions on the variables (constraints (4)). See [3,18,20,22-24]. Large-sized instances of the GAP have been tackled by means of approximate methods since exact methods have only proved to be efective for small-sized instances where the capacity constraints are loose. Different Lagrangean relaxations have been proposed in [2,9,13-15,17]. Other types of heuristic methods focused on providing good quality feasible solutions can be found in the literature for the GAP. Martello and Toth [19] propose a two phase heuristic method MTH. In the first phase an initial solution is constructed using a greedy function to measure the 'desirability' of assigning job j to agent i. In the second phase the initial solution is improved using a simple local exchange procedure. Cattrysse [4] considers a heuristic that incorporates a variable reduction procedure in order to reduce the size of the problem. This reduction procedure is based on the addition of valid inequalities (facets) of the knapsack polytopes associated to the capacity constraints. It starts solving the LP relaxation of the GAP. Then, violated valid inequalities for each knapsack constraint are identiied and added to the formulation. This process continues until no further valid inequalities can be generated. Once the reduction step inishes, all variables with value 1 are ixed and the capacity of each agent adjusted. In a second step, Simulated Annealing is applied to solve the reduced problem. Amini and Racer [1] developed a Variable Depth Search Heuristic (VDSH) for the GAP. Two phases are considered. In the irst phase, an initial solution and the LP lower bound are generated. A doubly-nested iterative reinement procedure is performed in the second phase. Feasible shift moves of a job from one agent to another one, and feasible swap moves between jobs assigned to diferent agents are considered during the reinement procedure. Trick [24] considers a heuristic approach based on the linear programming relaxation of the GAP. It is an iterative method which consists of three steps. In the irst step all 'useless' variables are removed (a variable is considered useless if ctjj > b,). In the second step the LP-relaxation is solved. During the third step all variables with value 1 are fixed. When the LP based heuristic inishes, an improvement phase that consists in swapping pairs of jobs is applied. Cattrysse et al. [5] reformulate the GAP as a set partitioning problem. A heuristic procedure based on this reformulation is proposed. Each column represents a feasible assignment pattern of jobs to an agent. For each agent, columns are generated solving a knapsack problem in which the coefficients are obtained from the vector of dual variables of the LP-relaxation. Since this problem is degenerate, a dual ascent procedure is used to obtain the dual variables. Then, a subgradient optimization procedure is applied to improve the lower bound. Heuristic solutions may be found by coincidence during the subgradient optimization or by search among the columns generated. Osman [21] proposes a hybrid heuristic SA/TS which combines Simulated Annealing and Tabu search. It uses a A-generation mechanism which describes how a solution can be altered to generate another neighbor solution. The SA/TS algorithm is called hybrid because it combines the nonmonotonic oscillation strategy of Tabu search with the simulating annealing philosophy. This nonmonotonic cooling scheme is used within the framework of simulating annealing to induce an strategic oscillation behavior in the temperature values. Also a Tabu search algorithm is proposed for the GAP. Both, the hybrid SA/TS and the Tabu search algorithms use a frequency-based 26 J.A. Diaz, E. Fernandez I European Journal of Operational Research 132 (2001) 22-38 memory that records information used for diversification purposes. A genetic algorithm (GA) for the GAP is proposed by Beasley and Chu in [7]. Instead of using a binary representation for the solutions, they use a n-dimensional vector of integer alphabets in the set /. These integer alphabets identify the agent to which each job is assigned. Fitness-unfitness pair evaluations are used to deal with both, feasible and infeasible solutions. Apart from using mutation and crossover operators, a two-phase heuristic improvement operator is also used: in the first phase this operator tries to recover feasibility by reducing the unfitness score. The second phase tries to improve the cost (fitness) of the solution without further violating the capacity constraints. Laguna et al. [16] introduce a Tabu search algorithm for the multilevel assignment problem MGAP based on ejection chains. MGAP differs from GAP in that agents can perform tasks at different efficiency levels. An ejection chain is an embedded neighborhood construction that considers the neighborhoods ofsimple moves to create more complex and powerful moves. Multiple dynamic Tabu lists and a long-term memory component are used within the framework of Tabu search. Also a strategic oscillation scheme that allows searching paths to cross the capacity-feasibility boundary is used. Yagiura et al. [27] have considered a Variable Depth Search algorithm VDS for the GAP. They incorporate an adaptive use of modiied shift and swap neighborhoods where some moves are forbidden in order to avoid cycling. The method also allows the search to visit infeasible solutions modifying the objective function to penalize the violated capacity of the agents. In Yagiura et al. [26], a branching search processes to construct the neighborhoods are incorporated in order to improve the performance of the VDS algorithm [27]. Simultaneously to our work, Yagiura et al. [25] have proposed a Tabu search algorithm where ejection chains are used to create more complex and powerful moves. The algorithm maintains a balance between visits to feasible and infeasible regions using an automatic mechanism for adjusting the parameter that weights the penalty term in the objective function. 4. Tabu search In this section, we present the components of our Tabu search algorithm for the GAP. First, the solution representation is described. Second, the neighborhood structures and their associated moves are presented. Then, a strategic oscillation scheme that permits the search to alternate between feasible and infeasible solutions is explained. Finally, the characteristics of the short-term memory and the frequency-based memory components of the algorithm are detailed. Tabu search, introduced in Glover [11] and [12], exploits a collection of principles for intelligent problem solving. It uses: (a) flexible memory structures; (b) an associated mechanism of control, to be used along with the memory structures, to deine Tabu restrictions and aspiration criteria and (c) records of historical information for diferent time spans that permit the use of strategies for intensiication and diversiication of the search process. Its power is based on the use of adaptive memory which is used for guiding the search process to overcome local optimality and to obtain near-optimal solutions. The use of short-term memory constitutes a form of aggressive search that tries to make good moves subject to the constraints implied by Tabu restrictions. Tabu restrictions are used to prevent cycling and aspiration criteria deine rules to ignore Tabu restrictions under certain conditions. Intensiication and diversii-cation strategies are attained by means of long-term memory. Intensiication strategies exploit features historically found desirable while diver-siication strategies encourage the search to explore unvisited regions. 4.1. Solutions A solution 7i for both GAP and RGAP is represented with a vector of n elements, where the kth component is the agent to which job k is assigned, i.e. 7 =(7i1,..., nn), where 77j G /; 77j = i ■<=>■ xij = 1, J.A. Diaz, E. Fernandez I European Jour 4.2. Neighborhood structures and generation mechanism Given the set S of feasible solutions for an optimization problem, a neighborhood structure defines, for each solution n G S, a set Sn c S of solutions that in some sense are 'close' to n.We have considered two simple classical neighborhood structures: shift and swap. Given a solution 7 , the shift neighborhood, Nshift, is defined to be the set of solutions that can be obtained by reassigning one job from one agent to another (see Fig. 1). AWt(n) = {(n1,..., n'n)\3j* G J s.t. n'f = nf, nj = nj Vj = The swap neighborhood, Nswap, is the set of solutions that can be obtained by interchanging the assignments of two jobs, originally assigned to diferent agents (see Fig. 2). of Operational Research 132 (2001) 22-38 27 Nswap(n) = {(n1, ... , n'n)\3j1, j2 G J, nj1 = nj2, s.t. nj1 = nj2> n^2 = nj1 > nj = nj Vj = j1> j2}. It is important to note here that when only swap moves are used the solutions have a special structure. In particular, solutions have a ixed number of jobs assigned to each agent. On the other hand, shift moves are more restrictive in terms of feasibility. For these reasons we have used N = iVshift U ^swap as neighborhood structure which is a more lexible option. A neighbor generation mechanism is a set of rules for selecting a solution from a given neighborhood. The following rules deine the generation mechanism used in our Tabu search implementation: • Moves in Nshift U Nswap that are not Tabu-active are considered admissible. • For the current assignment n =(rc1;..., nn), jobs are processed by decreasing values of zlnjj. • For a given job j, shift and swap moves are Fig. 2. Swap neighborhood. 28 J.A. Diaz, E. Fernandez I European Journal of Operational Research 132 (2001) 22-38 evaluated and the best admissible move with respect to objective function (1") is associated to the job. • The candidate generation mechanism halts if the best admissible move for the job being processed improves the objective function value (1").At each iteration, the candidate moves are regenerated. • Otherwise, if all the jobs have been processed but none of the associated moves improve objective function (1"), the admissible move with the lowest increment is selected and performed. 4.2.1. Remarks (1) Note that both, the considered neighborhood and the selection rules, are extremely simple. It is interesting to note that since the search is first improving, in most cases only a partial exploration of candidate moves is performed at each iteration. It is highly probable to ind an improving move with respect to objective function (1") before all jobs have been processed. (2) As already mentioned, using the above strategy, the search is guided by the objective function (1") and no priority is given to solutions that are feasible for the GAP. This does not happen in Laguna et al. [16], where the strategy of the search is guided towards feasibility through a qualitative criterion which makes prefer feasible solutions in most cases. Neither does this happen in Yagiura et al. [27], where the search strategy changes each time that feasibility is attained. In Yagiura et al. [25] the ejection chain moves do not consider solutions that increase infeasibility. Thus, their strategy also guides the search towards feasibility. (3) It can also be observed that with the above selection rules, no priority is given to shift moves with respect to swap moves and vice versa. Once more the element that takes care of guiding the search is the modified objective function (1"). This allows the solutions obtained at the various iterations to change their structure and to have different numbers of jobs assigned to each agent. To a great extent this does not happen in [27] which uses a more complex neighborhood structure also based on Nshift U iVswap. In that work, at some given iterations a shift move is forced and then only swap moves are permitted until feasibility is recovered. 4.3. Strategic oscillation As mentioned previously, the objective function used in the relaxed formulation RGAP includes a penalty to evaluate the violation of the agents' capacity constraints for the diferent solutions to RGAP. This penalty term is of the form P = s, , where p is a penalty parameter. Since the value of P is null for feasible solutions to the GAP, objective function (1") in RGAP provides the original evaluation (1') for such solutions and a modified evaluator for infeasible solutions. This modiied evaluator, together with the (previously described) rules of movement guide the search towards the domain of feasibility when dealing with solutions far away from such domain. Therefore, considering objective function (1'') enables the algorithm to have a strategic oscillation behavior that permits alternating between feasible and infeasible solutions. In this context, the weight assigned to the feasibility violation penalty, p, plays a central role to furnish the process with an equilibrium that permits alternating between feasibility and flexibility. In our algorithm this is done dynamically according to the expression p := p • a(infeas/(niter-1)-1); where infeas is the number of infeasible solutions to GAP obtained during the last niter iterations and a ^ 1 is a parameter that varies as explained later in this section. The above expression can be interpreted in terms of incorporating a recency-based memory (infeas) and a medium-term memory (a) to adapt the parameter p to the current stage of the search. The rationale behind the expression that we use is the following: the search strategy of our algorithm relies on the usefulness of a wider exploration of the infeasible solutions' space. For this J.A. Diaz, E. Fernandez I European Journal of Operational Research 132 (2001) 22-38 29 reason, we consider that variations on the value of p, although necessary, should be performed very smoothly. This permits the search to operate according to a uniform criterion for a series of iterations. In other words, it allows the paths of successive solutions to be generated with respect to homogeneous rules. For a given a, the value of p ranges in [a-1 • p, a1/^niter_1^ • p]. Note that, since the values of a will always be greater or equal to one, p will only increase when all the solutions found in the last niter iterations are infeasible. However, even in that case, the upper limit on the interval makes the increase in p to be very small. In the other cases, p is reduced (or maintained) increasing (or not decreasing) the importance assigned to the quality of the solutions versus their closeness to feasibility. If, in the above expression, a were fixed to 1, the value of p would also be fixed during all the procedure and the history of the search would not be taken into account. When a = 2, the expression obtained is very similar to the one used in [10] for Routing Problems with stochastic demands and customers. Observe, however, that when the value of a is ixed, only recency-based memory is used. The value assigned to a influences considerably the progress of the search. The bigger the value of a the larger the range of variation for p. Therefore, when only infeasible solutions are visited in the last iterations, the rate of increase of p is faster for higher values of a. Since most of the time the considered solutions will be infeasible, small values of a tend to give smaller p's and vice versa. Consequently, a medium-term memory component has been incorporated to see if good feasible solutions have been found with a similar rate of increase/ decrease for p. When this is so, the value of a is maintained to enhance the search under the same conditions on a given area. Otherwise, the value of a is increased guiding the search towards the feasible subregion of the area being explored. Again, the variation of a is extremely moderate to avoid sudden changes in the search. In particular, when the best feasible solution known so far has not been improved for the last 100 iterations, the value of a is increased by 0.005 every 10 iterations until a new better solution is found. Each time a new best feasible solution is found, the value of a is reset to the value 2. 4.4. Short-term memory phase Fig. 3 outlines the short-term memory phase of the algorithm. The short-term memory component of Tabu search incorporates lexible memory functions to forbid moves in order to avoid previously visited solutions. This memory is often managed via solution's attributes that are forbidden during some period of time. Short-term Memory Phase short-term jphase(k, -nk) Let 7r* be the best solution found so far, 7rk the current solution and k the total number of trial solutions generated while (stopping-criterion not satisfied ) do begin Select a neighbor solution 7rfc+1 £ N(7rk) (Using the neighborhood structure N and the neighborhood generation mechanism detailed in the previous section) if (cos*(7r*+1) < cost(7v*)) then Update short-term memory (tabu restrictions detailed below) Update long-term memory (frequenr y-based memory detailed below) k := k + 1 end Fig. 3. Short-term memory phase. 30 J.A. Diaz, E. Fernandez I European Journal of Operational Research 132 (2001) 22-38 The attributes we have considered are the assignments of jobs to agents. They are represented by ordered pairs (i, j) meaning that job j is assigned to agent i. In our implementation of Tabu search for the GAP some of these attributes are forbidden for the solutions to be generated in the next iterations. In other words, these attributes will be Tabu-active during a given number of iterations. More specifically, a matrix T is used to record the forbidden attributes for the future solutions. If during a given iteration k, a job j assigned to agent i1 is reassigned, then, any solution having the pair j) as attribute will be forbidden for the next t iterations. That is, 7]/ = k +1, records the last iteration number for which the attribute (z'i, j) will remain Tabu-active. In the case of swap moves only one pair is recorded as Tabu-active in the matrix T. From the two assignments being considered it is the one with the higher Anjj. Recency-based memory is managed using dynamic Tabu tenures. The parameter t is randomly selected within a range [tmin, fmax]. A standard aspiration criterion is used to bypass a Tabu restriction for a move. If, after performing the move, a feasible solution better than the best solution known so far would be obtained, the move is selected and performed. As usual, the short-term memory phase is terminated after a ixed number of iterations without improvement. 4.5. Long-term memory phase Long-term memory functions are used within the framework of Tabu search to furnish the search with a higher level of flexibility in order to achieve regional intensification and/or global di-versiication. The role of long-term memory is to record features of a selected number of trial solutions. Our TS algorithm uses a frequency-based memory scheme, both, for intensiication and diversification purposes. The frequency matrix fr records the number of times that job j has been assigned to agent i in the total number ofsolutions visited up to that point. The intensiication phase recovers the best solution found so far and ixes the current assign- ments of jobs to agents with a frequency of at least 85%. This assignments can be considered 'good' given that: first, they are also present in the best solution found so far. Second, if they appear in most of the solutions generated so far either they have been selected many times or after being selected they have not been changed for a long period of time. Fixing some assignments is equivalent to reducing the size of the problem given that some decision variables are set to the value 1. Once this ixing process inishes, the short-term memory phase restarts for the resulting reduced GAP with the unfixed jobs and adjusted agents' capacities. In the diversiication phase, the frequency-based memory is also used but in a different way. The cost matrix A is replaced by the matrix fr + A. Now the elements of the matrix fr can be viewed as penalties for selecting the most frequent assignments of jobs to agents. This is done to encourage the assignments that have seldom been selected and to discourage those assignments that have been frequently used. The matrix A is added to the frequency matrix for tie breaking purposes. During this phase it is intended to perform a number of high inluence moves in order to produce solutions with diferent structures from the ones generated up to this point. After performing such moves, the original cost matrix A is recovered and the short-term memory phase reinitiated. Our experience has shown that this diversiication scheme is efective since diferent high quality solutions from the ones generated to that point are usually obtained after a reduced number of iterations in the diversiication phase. Fig. 4 depicts the proposed Tabu search algorithm. It starts with one application of the short-term memory phase to improve the initial solution obtained using the MTH heuristic proposed by Martello and Toth [19]. Then, a cycle that alternates between intensiication and diversiication phases is performed during l iterations. 5. Computational results Three diferent sets of problems have been used to test the proposed Tabu search algorithm. The irst two sets of problems are publicly available in J.A. Diaz, E. Fernandez I European Journal of Operational Research 132 (2001) 22-38 31 Tabu Search Algorithm TabuSearch() Let 7r = (7ri,.. . , 7rn) be the best solution found so far and k the total number of trial solutions generated k := 0 Use MTH (Martelio &; Toth heuristic) to obtain an initial solution so Apply short-term.phase(k,7Tk) for / iterations do begin {Intensification Phase} Retrieve tt* and set ixk — 7r* if {.fv^j > 0.85 * k) then fix the value of xn* 1 to 1 3 Apply short-term-phase(k, 7rk) {Diversification Phase} Unfix all fixed assignments Restart short-term memory phase for very few iterations using the matrix cost A + fr Recover the original cost matrix A Apply short-term jphase(k, 7rk) end Fig. 4. Tabu search algorithm for the GAP. Beasley OR-library (http:IImscmga.ms.ic.ac.ukIjebI orliblgapinfo.html). The third set of problems is publicly available in Yagiura library (http:IIwww-or.amp. Lkyoto-u.ac.jpImembersIyagiuraIgapI). The first set, S1, contains 60 'small-sized' maximization problems. These problems were used to test the set partitioning heuristic of Cat-trysse et al. [5], the Simulated Annealing heuristic FSA of Cattrysse [4] as well as the hybrid Simulated Annealing/Tabu search heuristic SA/TS and the Tabu search heuristic TS proposed in Osman [21]. The test problems of this set have the following characteristics: 1. The number of agents m is set to 5, 8 and 10. 2. The ratio r = m/n is set to 3, 4, 5 and 6 to determine the number of jobs n. 3. aij values are integers generated from a uniform distribution U(5,25). 4. cij values are integers generated from a uniform distribution U(15,25). 5. The bi values are set to ((0.8/m) ^2j€j aij)• The test problems are divided into 12 groups (each containing ive problems), gap1 to gap12, according to the size of the test problems as is depicted in Table 1. The second set of problems, S2, contains 24 'large-sized' minimization problems used to test Table 1 Dimensions for S1 problems Problem group mn gap1 5 15 gap2 5 2O gap3 5 25 gap4 5 3O gap5 8 24 gap6 8 32 gap7 8 4O gap8 8 48 gap9 1O 3O gaplO 1O 4O gap11 1O 5O gap12 10 60 the GA proposed in [7], the Variable Depth Search algorithm proposed in [27], the VDS with branching search considered in Yagiura et al. [26], and the ejection chain approach proposed in Yagiura et al. [25]. The third set of problems, S3, contains: 3 Type C problems with n = 400, 3 Type D problems with n = 400 and 9 Type E problems used in [25]. Dimensions of S2 and S3 problems can be seen in Table 5. The problems in S2 and S3 are divided into ive classes according to the way they were generated: 1. Type A. aij are integers generated from an uniform distribution U( 5 , 25) , cij are integers 32 J.A. Diaz, E. Fernandez I European Journal of Operational Research 132 (2001) 22-38 generated from an uniform distribution U(10, 50) and Z>, = 0.6 x(n/m)x 15 + 0.4R, where R ma^ ) an '6/ . 4^ . j6J Jj=i and Ij = mm{i\cij < ckj Vk e I}. 2. Type B. aij and cij generated as in Type A problems and bi is set to 70% of the value given for Type A. 3. Type C. aij and cij generated as in Type A problems and bi = 0.8 EjeJ a*//m. 4. Type D. In this case aij are integers generated using an uniform distribution U (1,100), cij = 111 — aij + e, where e are integers generated using an uniform distribution U(—10,10) and bi = 0.8 E^a^/m. 5. Type E. alS = 1 — 10 x ln(Uniform(0,1]), clj = 1000/ay- — 10 x Uniform^,1], bi = max(0.8/m)x J2"=1, aij, max{aij Vj e J}. Type E problems were generated using the procedure detailed in Laguna et al. [16]. A fairly simplistic assumption with respect to the cost/resource relationship is assumed for Types A, B and C problems. However Types B and C problems are more difficult than Type A since the capacity constraints are tighter. Types D and E problems are known to be more difficult to solve. Results for Type A problems are not given since the MTH heuristic produces the optimal values in all cases. The experiments have been performed on a Sun Sparc Station 10/30 with 4 HyperSparc at 100 MHz using only one processor. The algorithm is coded in C language. For each of the test problems 30 runs were performed. After some tuning the following parameters' values were used: number of intensification and diversification iterations, 1 = 6; the Tabu tenures values range from ;min = 2to tmax = 6. The stopping criterion for the short-term memory phase has been adapted to the size of the test problems. In particular, for problems in S1, this parameter has been set to 350 iterations without improvement. For large-sized problems in ^bestp S2 and S3, the stopping criterion has been fixed to 1500 iterations without improvement. The value of a is initially set to 1. Starting with a smaller value of a has the effect of producing good quality feasible solutions at an earlier stage of the search. Once a feasible solution is found, the value of a is reset to the value 2 and it is permitted to take values between 2 and 3. Its value is reset to 2 every time a better solution is found. Arbitrarly, the value of p is initialized to 1. For each of the 60 problems of S1 the following values were computed: • Average percentage deviation from the optimal value, ffoptp = E3=1(Op - Spr)/30. • Average best-solution time, řbestpr/30. • Average execution time, řtotalp = E3=1 /30. Op denotes the optimal solution for problem p, and Spr is the solution obtained for problem p during the rth execution. řbestpr and itotal denote the times needed to reach the best solution and the total running time for problem p during run r, respectively. Additionally, the following values were calculated for each problem group (gap1 to gap12): • Average percentage deviation from the optimal value, ffgroup = E/L ffoptp/5. • Average best-solution time, řbests = J2p=1 tbestp/5. • Average execution time, = Ep=1 /5. Table 2 summarizes some characteristics of the methods used for comparison in our computational experiments. We have focused on the neighborhoods that are explored, the possibility of violating feasibility and the use of recent, medium-term and long-term memory. Tables 3 and 4 show the results obtained with our Tabu search algorithm (denoted TSf) for problems in S1. Results are compared to the hybrid Simulated Annealing/Tabu search SA/TS and the Tabu search TS1 heuristics proposed in [21] and also to the GA proposed in Beasley and Chu [7]. For SA/TS and TS1, Table 3 shows the average percentage deviation from optimum of the best solution found for each problem, over all the problems in each group (see [21]). For the GA and the TSf approaches, this table shows the average percentage deviations from optimum (ffgroup) taken over all the executions of the problems in each J.A. Diaz, E. Fernandez I European Journal of Operational Research 132 (2001) 22-38 33 Table 2 Algorithms' characteristics Moves Infeasibility Objective function Long term memory Shift move Swap move Ejection chain Allow search to visit infeasible solutions Moves towards feasibility are preferred Automatic penalty adjustment Recency memory Medium term memory Diversification Intensification SA/TS TS1 x x GA x x BVDS TS(MGAP) x x x xx xx TSEC x x x TS x x xx 1 it. 10 it 100 it x x x x x x Table 3 Average percentage deviation from optimal for S1 Problem set SA/TS TS1 GA TSf gap1 0.00 0.00 0.00 0.00 gap2 0.00 0.10 0.00 0.00 gap3 0.00 0.00 0.00 0.00 gap4 0.00 0.03 0.00 0.00 gap5 0.00 0.00 0.00 0.00 gap6 0.05 0.03 0.01 0.01 gap7 0.02 0.00 0.00 0.00 gap8 0.10 0.09 0.05 0.01 gap9 0.08 0.06 0.00 0.00 gaplO 0.14 0.08 0.04 0.03 gap11 0.05 0.02 0.00 0.00 gap12 0.11 0.04 0.01 0.00 Tall 0.210 0.070 0.009 0.004 °best 0.04 0.03 0.00 0.00 No. Optimal 39 45 60 60 solutions group. Also, for each of the approaches, the following averages over all groups of problems are depicted in Table 3: • ffall: Global average percentage deviation from optimum. • ffbest: Average percentage deviation from optimum of the best solutions found. • Number of problems for which the optimal solution was found at least once. Table 4 shows the average best-solution time, tbest, and the average execution time, ttotal, for each compared method and for each group of problems. Times are measured in seconds. It can be seen that our Tabu search algorithm provides the optimal values for all the problems. In this sense it outperforms SA/TS and TS1 and is equivalent to GA. With respect to the deviation from optimal values it never exceeds 0.03% and, therefore, it outperforms all the compared methods. Although the average best-solution times and average execution times are not comparable since the experiments were carried out in diferent computers, it must be noticed that for TSf the values tbest and ttotal are extremely good (recall that we only use one processor). We next show the results obtained for problems in S2 and S3. We compare our results to the best known solution for these problems. They are shown in the last column of Table 5. For ive Type C instances, one Type D instance and six Type E instances this value corresponds to the optimal solution value obtained with the branch-and-bound algorithm proposed by Nauss [20]. These instances are marked with an asterisk. For the other instances, the best known values are not known to be optimal and have been taken from [8] in the case of Type B instances and from [25] for the remaining Types C, D and E instances. Now, for each of the problems the following values were computed: 34 J.A. Diaz, E. Fernandez I European Journal of Operational Research 132 (2001) 22-38 Table 4 CPU Times for Sla Problem group SA/TSb TS1b GA 0 TSf d 4ota] 4ota] ^best 4ota] ^best 4ota] gap1 n.a. 0.22 n.a. 0.73 0.40 72.38 0.04 1.09 gap2 n.a. 0.91 n.a. 1.44 1.68 79.96 0.04 1.48 gap3 n.a. 1.82 n.a. 2.85 2.18 85.06 0.08 2.05 gap4 n.a. 2.32 n.a. 4.86 4.16 92.36 0.13 2.54 gap5 n.a. 2.48 n.a. 3.97 5.52 100.36 0.22 2.38 gap6 n.a. 5.67 n.a. 8.85 16.36 130.02 0.41 3.53 gap7 n.a. 13.97 n.a. 12.72 10.54 130.34 0.35 4.92 gap8 n.a. 19.53 n.a. 22.99 49.90 184.42 2.38 7.64 gap9 n.a. 5.73 n.a. 12.19 13.82 129.32 0.77 3.72 gap10 n.a. 15.55 n.a. 18.62 33.70 167.66 1.07 5.90 gap11 n.a. 36.96 n.a. 34.46 32.60 193.00 1.53 8.02 gap12 n.a. 65.96 n.a. 47.07 39.08 260.34 1.83 10.67 an.a. = Not available. bCPU times in seconds on a VAX-8600 computer. cCPU times in seconds on a Silicon Graphics Indigo R4000 (100 MHz). dCPU times in seconds on a Sun Sparc Station 10/30 with 4 HyperSparc 100 MHz (SPECint95 2.35). • Average percentage deviation from best known solution, V - Bp)/30. • Average best-solution time, ?bestr. • Average execution time, ftotalp. Bp is the best known solution value for problem p. Tables 5 and 6 show the results obtained. Now, TSf is compared to the GA proposed by Beasley and Chu in [7], to the Variable Depth Search with Branching algorithm BVDS proposed by Yagiura et al. [26], to the Tabu search algorithm for the MGAP, denoted TS(MAGAP), proposed by Laguna et al. [16], and to the Ejection Chain algorithm TSEC proposed by Yagiura et al. [25]. Results for BVDS and TS(MGAP) were taken from [25]. For each problem, Table 5 shows the best values (best) and the average percentage deviation from the best known solution (