Lagrangian relaxation guided problem space search heuristics for generalized assignment problems

European Journal of Operational Research 182 (2007) 1039–1056 www.elsevier.com/locate/ejor

Discrete Optimization

Lagrangian relaxation guided problem space search heuristics for generalized assignment problems V. Jeet, E. Kutanoglu

*

Operations Research and Industrial Engineering Program, The University of Texas at Austin, Austin, TX 78712-0292, USA Received 23 September 2005; accepted 18 September 2006 Available online 8 December 2006

Abstract We develop and test a heuristic based on Lagrangian relaxation and problem space search to solve the generalized assignment problem (GAP). The heuristic combines the iterative search capability of subgradient optimization used to solve the Lagrangian relaxation of the GAP formulation and the perturbation scheme of problem space search to obtain highquality solutions to the GAP. We test the heuristic using different upper bound generation routines developed within the overall mechanism. Using the existing problem data sets of various levels of difficulty and sizes, including the challenging largest instances, we observe that the heuristic with a specific version of the upper bound routine works well on most of the benchmark instances known and provides high-quality solutions quickly. An advantage of the approach is its generic nature, simplicity, and implementation flexibility.  2006 Elsevier B.V. All rights reserved. Keywords: Assignment; Generalized assignment problem; Heuristics; Problem space search; Lagrangian relaxation

1. Introduction The generalized assignment problem (GAP) has many applications and can be usually found as a subproblem within more complex problems. A non-exhaustive list of such problems includes computer network problems, vehicle routing and scheduling problems, and resource allocation and scheduling problems in manufacturing and service. We use machine scheduling terminology to introduce and present the development of the algorithm tested in the paper; this terminology and the associ-

*

Corresponding author. E-mail address: [email protected] (E. Kutanoglu).

ated approach are directly applicable to other domains with little change. In this machine scheduling environment, the GAP attempts to find the minimum-cost assignment of a set of jobs with certain resource requirements to a given set of capacitated machines, without violating the machines’ capacity limitations. There are many exact algorithms and heuristics developed to solve the GAP. Exact algorithms include branch and bound-based algorithms, branchand-cut and branch-and-price algorithms (e.g., Ross and Soland, 1975; Nauss, 2003; Savelsbergh, 1997). Some heuristics use the linear programming relaxation (Trick, 1992). Others use search techniques such as genetic algorithms (e.g., Chu and Beasley, 1997), and tabu search algorithms (e.g., Higgins, 2001; Dı´az

0377-2217/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2006.09.060

1040

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

and Ferna´ndez, 2001; Yagiura et al., 2004). Lagrangian relaxation (LR) is arguably one of the heavily used methods to solve general integer programming problems (Fisher, 1981). Most papers on GAP also have applied Lagrangian relaxation to the GAP and relax one of the two constraint sets. Some of these papers have proposed combining LR with surrogate relaxation (e.g., Narciso and Lorena, 1999). The Lagrangian decomposition (LD) method has been also applied to GAP (e.g., Haddadi, 1999). More recently, LR and subgradient methods have been integrated in branch and bound schemes (e.g., Haddadi and Ouzia, 2001; Haddadi and Ouzia, 2004). In this paper, we present another algorithm that combines LR and its associated subgradient optimization procedure with a search technique based on problem space search (PSS). PSS uses the problem input data perturbations (as opposed to traditional solution space search, where the search is directly over solutions). The idea of perturbation is not new and has been used in the literature before. Storer et al. (1992) present the first example of PSS within the context of scheduling. Charon and Hudry (1993) use a very similar noising method applied to clique partitioning problem. Their noising methods are based on elementary (or local) transformations applied to solution(s). Codenotti et al. (1996) use the perturbations applied to intermediate solutions of large-scale Euclidean TSP benchmark instances. Here the main difference is that PSS aims the ‘‘noising’’ or the ‘‘perturbations’’ applied to the problem data rather than intermediate solution(s). The main idea is to combine the strengths of LR and PSS to find improved solutions to the GAP. Our goal is to present the general framework of the LR-PSS-based algorithm, especially in the context of the subgradient optimization used to solve the Lagrangian relaxation problem, and to show the effectiveness of the combined approach by applying it to an extensive set of challenging benchmark instances from the literature. To give an overall idea, we first construct a Lagrangian relaxation of the GAP by relaxing the machine capacity constraints in the integer programming formulation of the GAP. The relaxed problem’s solution gives a lower bound on the optimal solution of the GAP, which can be improved using the subgradient optimization procedure. As the lower bound solutions are rarely feasible, the solution of the relaxed problem in each iteration

of the procedure is also used to find heuristic solutions to the original GAP and produce upper bounds. As this ‘‘feasibility restoration’’ in fact does not guarantee capacity-feasible solutions due to its structure, a PSS-based enhancement is introduced: The problem data in each iteration is temporarily perturbed a given number of times to increase the chances of obtaining capacity-feasible solutions with potentially improved objective values. At algorithm termination, the best feasible solution and the best lower bound are reported. The two primary factors affecting the PSS enhancement are (1) the number of perturbations in each iteration, and (2) the amount of perturbation of the problem data. We test the effectiveness of the algorithm using the problem instances from the literature. 2. Literature review Ross and Soland (1975) define the GAP and Fisher and Jaikumar (1981) proves that it is NPhard. However, it is important to solve this problem efficiently as its applications are found in many real life problems addressing process assignments to computer networks (Balachandran, 1976), facility location (Ross and Soland, 1977), vehicle routing and its extensions such as constraints on total travel and delivery times for the vehicles (Fisher and Jaikumar, 1981; Baker and Sheasby, 1999). The GAP is under investigation for more than three decades now and many techniques have been developed to solve this problem. To name a few, Ross and Soland (1975) propose a branch and bound algorithm, Fisher et al. (1986) suggest a multiplier adjustment method, and Trick (1992) introduces a heuristic based on LP relaxation. Cattrysse and Wassenhove (1992) present an extensive survey of algorithms for the GAP published until 1992. Amini and Racer (1994) use statistical analysis based on their experimental design and compare four leading algorithms against their new variabledepth search algorithm. Cattrysse et al. (1994) propose a set partitioning algorithm. Among metaheuristics-based approaches, Osman (1995) introduces simulated annealing methods and Laguna et al. (1995) use Tabu search. Chu and Beasley (1997) test a genetic algorithm for the GAP. Wilson (1997) suggest a simple dual based algorithm, while Savelsbergh (1997) introduces branch and price approach. Yagiura et al. (1998) propose algorithms based on variable-depth search with

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

branching. Narciso and Lorena (1999) suggest to combine Lagrangian and surrogate relaxation for the GAP. Romeijn and Romero-Morales (2000) present a class of greedy algorithms based on several measures for the GAP. Recent work on GAP includes Higgins (2001), a dynamic tabu search algorithm for large-scale instances for the first time, Farias and Nemhauser (2001) show a family of inequalities for the GAP polytope, and Cario et al. (2002) investigate the relationship between problem characteristics and algorithm performance. Dı´az and Ferna´ndez (2001) test a tabu search heuristic with two main advantages, its simplicity and flexibility. More recently, Nauss (2003) combines several ideas in his algorithm including the cuts suggested by Gottlieb and Rao (1990) to improve the Lagrangian lower bound and show the effectiveness of his algorithm for GAP instances of up to 3000 binary variables (e.g., 30 · 100 or 15 · 200, where m · n represents a GAP instance with m machines and n jobs). To our knowledge, the most recent work on the GAP is Yagiura et al. (2006) – they test an algorithm based on path relinking combined with ejection chain neighborhood. Yagiura et al. (2004) also propose a tabu search-based algorithm with ejection chain neighborhood. Yagiura and Ibaraki (in press) provides an overall review of the stateof-the-art metaheuristics and approximation algorithms for the GAP. Nauss (2006) discusses the latest integer programming based algorithms for the GAP. In this paper, we investigate the performance of a problem space search-based enhancement on top of subgradient optimization used to solve the Lagrangian-relaxed GAP. Storer et al. (1992, 1995) introduce and test PSS for sequencing and job shop scheduling problems. Later, Storer et al. (1996) and Naphade et al. (1997) extend the approach to other problem types, including number partitioning and project scheduling. In this paper, we introduce another extension of PSS, for the first time used in conjunction with Lagrangian relaxation and subgradient optimization, to find provably high-quality solutions to the GAP. We instantiate two different implementations of the overall algorithm and a hybrid between the two to calculate upper bounds using the knowledge extracted from problem data perturbations used as part of PSS within the subgradient optimization framework. Hence, the overall algorithm is quite simple, flexible, and generic in its nature. We present our computational results

1041

for all the GAP instances listed by Chu and Beasley (1997) and by Yagiura et al. (2004). Our main comparisons will use the extremely competitive results by Yagiura et al. (2004) and equally challenging results by Dı´az and Ferna´ndez (2001). 3. Generalized assignment problem We can define the GAP as assigning n jobs to a set of m capacitated resources/machines with minimum cost. Each machine i has a limited amount bi of capacity and each job j must be assigned to one and only one machine. Assigning job j to machine i (subject to the availability of machine i) requires aij units of machine i’s capacity, which costs cij when assigned. Let xij (defined for all i and j) be a binary decision variable indicating whether (xij = 1) or not (xij = 0) we assign job j to machine i. We model the problem as the following integer program: Minimize

m X n X i¼1

subject to

n X

ð1Þ

cij xij

j¼1

aij xij 6 bi

8i ¼ 1; . . . ; m;

ð2Þ

j¼1 m X

xij ¼ 1

8j ¼ 1; . . . ; n;

ð3Þ

i¼1

xij 2 f0; 1g 8i ¼ 1; . . . ; m; j ¼ 1; . . . ; n: Objective (1) is the total cost of the assignment. Constraints (2) are the machine capacity constraints, and constraints (3) ensure that each job is assigned to only one of the machines. Similar to the relaxation scheme used by Dı´az and Ferna´ndez (2001), we relax the machine capacity constraints (2) using Lagrange multipliers ki (defined for all i), and obtain the following relaxed model (denoted by LR(k) for given vector k of m ki’s): Minimize

m X n m X X ðcij þ ki aij Þxij  ki bi i¼1

subject to

m X

j¼1

xij ¼ 1

ð4Þ

i¼1

8j ¼ 1; . . . ; n;

ð5Þ

i¼1

xij 2 f0; 1g 8i ¼ 1; . . . ; m; j ¼ 1; . . . ; n:

For given k, we can solve the relaxed problem (LR(k)) efficiently, by inspection: We first define so-called ‘‘weights’’ wij(k) for each (i, j) pair: wij ðkÞ ¼ cij þ ki aij :

ð6Þ

1042

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

We then assign each job j to the machine that offers the minimum weight for the job: ij 2 arg min wij ðkÞ:

ð7Þ

i

We finally set associated xij ;j ¼ 1, and set xij = 0 for all i 6¼ ij . Then, the optimal value of the relaxed problem for given ki’s is vðLRðkÞÞ ¼

n X j¼1

wij ;j ðkÞ 

X

ki bi ;

ð8Þ

i

where v(P) denotes the optimal value of problem P. The relaxed problem’s optimal solution value is a lower bound on the optimal value of the original problem, GAP: v(LR(k)) 6 v(GAP). A practical and efficient method to obtain improving lower bounds is the subgradient optimization procedure, which iteratively solves the Lagrangian-relaxed problem and updates the multipliers in a systematic, but simple way. As the GAP is NP-hard and suffers from strictly positive duality gap, the solutions of the relaxed problem are rarely capacity-feasible, let alone optimal. Hence, we also use the subgradient optimization procedure as a basis for the heuristic development, which tries to restore feasibility in the capacity-infeasible solutions of the relaxed problem in every iteration. Here, we try to obtain feasible, high-quality solutions to GAP, which serve as upper bounds on the value of GAP. A simple and generic idea for feasibility restoration, for which we test two specific versions, is to assign the jobs in a certain order one by one to a machine with enough capacity, making sure that the overall assignment is always feasible. To design a greedy constructive heuristic we need to fix the order of assignment of the jobs and the criterion used to assign one. The key questions to answer are (1) which job do we pick next to assign (the order of jobs that we use to assign them to the machines), and (2) which machine do we assign the picked job. As we assign the jobs one at a time, the assigned jobs determine the remaining capacity for the unassigned jobs. A job that is assigned later in the process may not find a enough leftover capacity in any of the machines (leading to infeasibility) or may have to be assigned to a machine that is particularly expensive (leading to an inferior solution) due to insufficient capacities at other, more favorable machines. To see the effects of specific choices on these issues, we instantiate two different versions of the heuristic. Furthermore, we address both infeasibility and solution quality issues by enhanc-

ing both versions using the problem space search (PSS) mechanism. We also introduce a hybrid approach that aims to be more efficient computationally with a marginal sacrifice in solution quality. We first discuss the versions of the feasibility restoration heuristic and then explain enhancements provided by PSS. 4. Feasibility restoration heuristic The two ideas (denoted by LBR and RGT) for finding feasible assignments and calculating upper bounds are based on two different strategies. LBR is based on finding a good ‘machine’ for a randomly selected ‘job’ where the ‘goodness’ is derived from the lower bound solution, and RGT is based on finding a good ‘ordering’ of ‘job-machine’ pairs where the ordering is derived from a regret (or desirability) calculation. The hybrid version (HBR) combines the two ideas into a single scheme. All versions of the heuristic make use of the weights (wij values) defined in the previous section. Whenever we have good multipliers, the weights carry useful information both from the original assignment costs and from potential ‘‘resource usage costs’’ expressed using the Lagrangian multipliers as the Lagrangian multiplier value times the resource requirement of the assignment. That is, weight wij(k) for a potential assignment between machine i and job j not only represents the original assignment cost cij, but also carries resource load information based on the machine’s Lagrange multiplier ki and the job’s resource requirement aij. In some sense, a weight of an assignment represents how costly an assignment will be and how much additional load it will bring to the involved machine. 4.1. Lower bound replication (LBR) With the lower bound replication (LBR) version of the heuristic, we try to preserve (or replicate) the lower bound allocation (solution x = (xij) of LR(k)) as much as possible without violating the capacity constraints of machines when assigning jobs to machines in a random order. If it is not possible to keep the lower bound assignment due to used capacities of machines, we choose the machine with lowest weight that has sufficient leftover capacity. To use the weights for the heuristics, for given job j, we look for a machine i with the lowest weight wij(k) for the job, among the machines that have sufficient capacity to process the job. There is no trivial

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

way to order the jobs here, due to the difficulty to predict the effects of earlier assignment decisions on both the total cost and the infeasibility. We use a randomly generated sequence for such ordering. The overall idea behind LBR is to modify the lower bound assignments only when necessary, expecting little changes will not increase the lower bound too much over the optimal solution value, hopefully leading to tight upper bounds. Algorithm LBR Input: wij :¼ wij(k) for all (i, j) pairs. Let x be the vector of optimal values of the relaxed (lower bound) problem (LR(k)) for given k, i.e., x = (xij)i=1,. . .,m, j=1,. . .,n. Start: Set J :¼ {1, . . . , n}, b0i :¼ bi for all i, and UB :¼ 0. Step 1: Pick a job (denoted by ^j) randomly from J. Step 2: Let i 0 be the machine that job ^j is assigned in the lower bound solution, i.e., xi0 ;^j ¼ 1. If ai0 ;^j 6 b0i then let ^i :¼ i0 . Else, construct set I 0 ¼ fi 2 I : ai;^j 6 b0i g. If I 0 = ;, then STOP and report ‘‘NO FEASIBLE SOLUTION’’. Else, let ^i 2 arg mini2I 0 fwi;^j g. Step 3: Assign job ^j to machine ^i and update parameters: • Set x^0i;^j :¼ 1 and x0i;^j :¼ 0 for all i 6¼ ^i • Set UB :¼ UB þ c^i;^j • Set b^0i :¼ b^0i  a^i;^j . Step 4: Set J :¼ J n ^j. If J = ; then, STOP and report the solution x0 ¼ ðx0ij Þ with value UB. Else, go to Step 1. 4.2. Greedy regret (RGT) The main idea of this version of the feasibility restoration heuristic is to minimize the ‘‘regret’’ of not allocating a job to its most ‘‘favorable’’ machine (the machine with minimum weight for that job). It orders the jobs on the basis of ‘‘desirability’’ defined using the problem data and the weights. For each job, the difference between the weights of the second most favorite and the most favorite machine is computed, and the jobs are assigned in non-increasing order of this difference. The goal is to minimize the potential ‘‘regret’’ (or maximize the ‘‘desirability’’) by assigning the jobs with large regrets first (otherwise they would have to be assigned to the second most favorite (or worse, to other machines) and increase the total cost unnecessarily).

1043

For given weights wij’s, the desirability of assigning job j to its best machine is given by fwi0 ;j  wij g: qj ¼ max min 0 i

i 6¼i

ð9Þ

Due to capacity constraints on the machines, the desirability for each job is calculated over feasible machines for that job, i.e., those machines that have sufficient capacity to process it. The early versions of algorithms that use regret/ desirability for solving the GAP (which seem to be inspired by Vogel’s method for the transportation problem (Reinfeld and Vogel, 1958)) are originally suggested by Martello and Toth (1981) who also introduce several desirability measures. Romeijn and Romero-Morales (2000) discuss a few extra alternative expressions to be considered as ‘‘desirability,’’ and show the relationship between linear programming relaxation of GAP and these desirability definitions. We use their version as the basis for our specific implementation outlined below: Algorithm RGT Input: wij :¼ wij(k) for all (i, j) pairs. Start: Set J :¼ {1, . . . , n}, b0i :¼ bi for all i, and UB :¼ 0. Step 1: Let F j :¼ fi : aij 6 b0i g for all j 2 J. If Fj = ; for some j, STOP, and report ‘‘NO FEASIBLE SOLUTION.’’ Else, let • ij 2 arg min{wij : i 2 Fj}, for all j 2 J • qj :¼ mini2F j ;i6¼ij fwij  wij ;j g, for all j 2 J. Step 2: Let ^j 2 arg maxj2J fqj g. Let ^i :¼ i^j . Step 3: Assign job ^j to machine ^i and update parameters: • Set x^0i;^j :¼ 1 and x0i;^j :¼ 0 for all i 6¼ ^i • Set UB :¼ UB þ c^i;^j • Set b^0i :¼ b^0i  a^i;^j . Step 4: Set J :¼ J n ^j. If J = ; then, STOP and report the solution x0 ¼ ðx0ij Þ with value UB. Else, go to Step 1. 4.3. Problem space search As we outline above, there is a possibility that the feasibility restoration heuristic stops without finding a feasible solution. This could be viewed as a disadvantage of both versions of the feasibility restoration heuristic. To increase the chances of obtaining a feasible outcome from the heuristic and potentially obtain improved upper bounds in each iteration of the subgradient optimization procedure, we utilize a scheme called Problem Space Search (PSS).

1044

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

PSS is a simple metaheuristic technique that generates a whole new search space within the ‘‘problem space’’ by temporarily perturbing problem instance data instead of traditional strategy of searching within the ‘‘solution space.’’ The idea is to use the problem space search so that, in every iteration of subgradient optimization, we evaluate multiple solutions obtained from multiple executions of the specific algorithms outlined above, each done using temporarily perturbed weights as input. The idea behind PSS is to generate artificial ‘‘neighbor’’ problems by temporarily perturbing the problem input data and to find solutions to the perturbed problems using a fast heuristic. The main premise here is that as we find heuristic solutions to ‘‘similar’’ problems to the original one, we have a chance that we find solutions close to the optimal solution to the original problem. We extend the idea here to the Lagrangianrelaxed GAP and use it in each iteration of subgradient optimization. In our implementation, we find heuristic solutions to a given number of problems generated using the set of perturbed weights. As the weights are the main input to the algorithms, we temporarily modify the weights wij(k)’s in a controlled manner and find the corresponding heuristic solutions to them. We denote perturbed weights by w0ij ðkÞ’s and generate them as follows:  w0ij ðkÞ ¼ wij ðkÞ þ b  U ð1; 1Þ  w

8i; j;

ð10Þ

where b is a parameter that controls the amount of  is the average of all original weights, perturbation, w and U(1, 1) is a random number between 1 and 1 drawn from a continuous uniform distribution. We generate a given number of perturbations of weights (called neighborhood size within the PSS framework) and find algorithms’ solutions for each, keeping the best one (evaluated using the original assignment costs) within the ‘‘neighborhood.’’ In the algorithmic framework, we can view this PSS enhancement as an outer loop outside the iterations of the individual algorithms, feeding each set of the perturbed weights to the specific version (LBR or RGT) of the heuristic. The PSS enhancement to be implemented at each subgradient iteration after the LR problem is solved is as follows: • Set UBp at a very large number. • For a given number of perturbations (neighborhood size) do – Perturb original weights using (10) to obtain w0ij ðkÞ for all (i, j) pairs.

– Input wij :¼ w0ij ðkÞ for all (i, j) pairs, and execute the specific version of the feasibility restoration heuristic (LBR or RGT) and obtain UB 0 with solution x 0 . – Update UBp if the reported UB 0 is an improvement i.e., if UB 0 < UBp then set UBp :¼ UB 0 and xp :¼ x 0 . • Report UBp and xp. Note that xp and UBp are the best solution and its value found for all the perturbed weights for given k in a subgradient iteration. As we update k over the iterations of subgradient optimization, we have the same opportunity to execute this PSS enhancement over the standard runs of the algorithms in each iteration. The previous versions of the greedy regret algorithm produce promising results. Our results also confirm this, but due to additional computations needed for all jobs and all machines every time to choose an assignment at each iteration of the algorithm, the computational burden here is more substantial when compared to that of LBR. We comment on this tradeoff between the improvement in the solution quality obtained RGT over LBR and its extra computational burden later when presenting our results. The hybrid version of the heuristic tries to address this by applying RGT once in every subgradient iteration using the original (unperturbed) weights and applying LBR for the perturbed weights, for a given number of times in every subgradient iteration. The weight perturbation done through problem space search is explained next, followed by a brief discussion of the hybrid approach (HBR). We next outline the overall algorithm within this context.

5. Overall algorithm We have discussed how we compute lower bounds and upper bounds, and PSS-based neighborhood search for given k. As the main input to all the aforementioned components of the algorithm, the quality of Lagrange multipliers has significant effect on the performance of the overall algorithm. A good set of multipliers should produce a good set of weights, which in turn should lead to good lower bound and upper bound solutions. Hence, the overall idea behind Lagrangian relaxation is to update Lagrange multipliers to improve lower and upper bounds and reduce infeasibilities

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

that exist in the lower bound solutions in an iterative fashion. In the algorithmic framework, this introduces yet another outer loop that keeps track of Lagrange multipliers and updates them using subgradient optimization. We calculate a rather loose, initial upper bound using a naı¨ve method (Step 0). If the method in Step 0 cannot compute a feasible upper bound then we can safely declare that the instance is infeasible as the method focuses on feasibility with no regard to its cost. We then execute the main body of the algorithm, which takes in the Lagrange multipliers that are initialized at zero. Input: Problem data: aij, cij for all i and j, and bi for all i. Pick an algorithm for feasibility restoration: LBR, or RGT. Set iteration counter r = 1, set n = number of jobs, and set the number of perturbations (p). Step 0: Calculate an initial UB with the method outlined below: (a) Set j :¼ 1 and UB :¼ 0. Initialize b0i ¼ bi for all i. (b) Find ^i 2 arg mini2I:aij 6b0i faij =b0i g. (c) Set UB :¼ UB þ c^i;j and update ^i’s capacity ðb0i :¼ b0i  a^i;j Þ. (d) Set j :¼ j + 1. If j = n + 1, STOP and go to Step 1. Else, go to Step 0.(b). Step 1: UB* :¼ UB, LB* :¼ 0, a :¼ 2.0, ki :¼ 0 for all i. Step 2: Find solution x of LR(k), and compute LB. (a) Set LB :¼ 0. Set wij(k) :¼ cij + kiaij for all i and j. Set j :¼ 1. (b) Let ^i 2 arg mini2I fwij ðkÞg. (c) Set LB :¼ LB þ c^i;j . (d) Set j :¼ j + 1. If j = n + 1, go to Step 2.e. Else, go to Step 2.b. (e) Keep the best lower bound: If LB > LB* then, set LB* :¼ LB. Go to Step 3. Step 3: Execute the specified feasibility restoration algorithm, using weights wij :¼ wij(k). If a feasible solution x 0 (with value UB) better than UB* is reported, then update UB* and x* : UB* :¼ UB and x* :¼ x 0 . Step 4: Execute PSS for the given number of perturbations p using perturbed weights w0ij ’s as input to the specified algorithm: wij :¼ w0ij ðkÞ. Denote its best upper bound over its current neighborhood as UBp and the corresponding solution as xp. If UBp improves UB*, then update UB* and x*: UB* :¼ UBp and x* :¼ xp.

1045

Step 5: Update the subgradient and PSS-parameters, Lagrangian multipliers and update r :¼ r + 1. Step 6: If a stopping criterion is satisfied, then STOP and report LB*, UB* and x*. Otherwise go to Step 2. To further specify the algorithm’s details, we need to detail (1) how we update step size a used in our subgradient optimization implementation (2) how we stop the algorithm (Step 6), and (3) how we calculate the perturbed weights using perturbation factor, b (Step 4). We give further information on these in the next section. 6. Algorithmic settings The development and successful implementation of new heuristics often require tuning the parameters of the algorithm, as their effects on the algorithmic performance are not exactly known. We too have several parameters in the overall generic feasibility restoration heuristic and in the subgradient optimization procedure. Their proper selection is crucial for the quality of the solutions, the lower bounds, and the overall computational burden and convergence of the algorithm. Once these are set, however, they remain fixed throughout the execution of the overall algorithm. The parameters which affect the quality and speed of the overall algorithm are categorized in three classes below: Subgradient optimization parameters a-halving: A critical strategy on the convergence of the subgradient optimization procedure is to properly decrease the value of scale parameter a used in the step size calculation over the iterations. In an effort not to do this in a premature way, it is usually suggested to halve a within a number of iterations without any LB improvement. Deciding after preliminary experiments, we halve the value of a whenever there has not been any update (improvement) in the LR lower bound LB* for 10 consecutive iterations in our implementation. a-reset: Whenever there is an update (improvement) in the best upper bound UB*, we reset the scale parameter a back to 2. Stopping criterion: We use a three-condition stopping criterion: (1) a specified maximum number

1046

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

of iterations K have been completed, (2) the scale parameter a has become too small (less than a limit, denoted by a 0 ) to cause any changes in Lagrange multipliers, and (3) there has been no update (improvement) in the upper bound for a fixed number of consecutive subgradient iterations (denoted by k). If one of these conditions is satisfied, then we stop the algorithm. We set K = 1500, a 0 = 0.0025, and k = 200 in our implementation. The fourth stopping criterion (based on a pre-specified time limit) is used when we compare the heuristic results and CPU times with those in the literature.

Problem space search parameters b-update: Recall the value of b controls the amount of perturbation applied to the weights. We tested several versions, including a constant factor, and several decreasing beta’s over subgradient iterations. Due to its slightly better performance, we use b = 1/(1 + r) where r is the current subgradient iteration number. The idea behind this is to modify the weights in a rather limited fashion as the Lagrange multipliers become more robust, stable, and close to their ultimate values. Neighborhood size, p: This is the number of perturbations to input to the specific version of the feasibility restoration heuristic in each subgradient iteration. The preliminary experiments show definite improvement due to perturbations, but also indicate that small number of perturbations over more subgradient iterations (outer loop) produce better solutions than a large number of perturbations over a small number of subgradient iterations. Hence, we set p = 10 for all experiments (except for the limited time-runs conducted to compare the performances of multiple heuristics, where we set p = 5). We discuss these in more detail when presenting the computational results.

producing good solutions as compared to LBR. Therefore, we introduce a hybrid version (HBR), which runs RGT once with the original weights (Step 3 in the overall algorithm), and LBR for the rest of the PSS-perturbation iterations (Step 4). 7. Test problems There are five types of benchmark GAP instances (Chu and Beasley, 1997; Laguna et al., 1995) called types A, B, C, D, and E, generally in the increasing order of difficulty. We have two sets of data, differentiated according to their source: (1) Well-known Beasley’s OR-library data set that consists of smaller, once-challenging instances, and (2) the more recent data set from Yagiura et al. (2004) that consists of more challenging and larger instances (downloaded from http://www.al.cm.is.nagoya-u. ac.jp/~yagiura/gap/). The first set has 6 instances (up to 20 machines and 200 jobs) for each instance (type A through E – total of 30 instances). The second set has 9 instances (up to 80 machines and 1600 jobs) of types C, D and E (27 instances overall). A and B-type instances are considered relatively easy instances. Type C, D, and E instances (Yagiura et al.’s E instances are generated slightly differently as compared to the original ones reported by Laguna et al. (1995)) are more challenging instances due to their careful design in their data generation. The following shows how these instances types differ in their designs: Type C: aij’s are random integers from U[5, 25], cij areP random integers from U[10, 50], and bi’s are 0:8 j aij =m. Type D: aij’s are random integers from U[1, 100], cij = 111  aij + e1, where e1 is a random P integer from U[  10, 10], and bi’s are set 0:8 j aij =m. Type E: aij = 1  10 lne2, where e2’s are random numbers from U[0, 1], cij = 1000/aij  10e3, where e3’s arePrandom integers from U[0, 1], and bi’s are set 0:8 j aij =m.

Feasibility restoration algorithm 8. Computational results Heuristic Version: As we have already discussed, we have two versions, namely LBR and RGT, for finding feasible solutions and calculating upper bounds. In preliminary tests, we observe an increase in the computation time of RGT as we increase the problem size while consistently

The LR-PSS-based algorithm is programmed in C++ and run under the Suse Linux (version 2.6.8–24.14-smp) operating system. We perform the computational runs on a PC with an AMD Athlon 1.8 GHz processor. We present the LR-PSS

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

computational results in 3 tables and compare them with the results from Dı´az and Ferna´ndez (2001) and even more recent results obtained by Yagiura et al. (2004). Their efficient Tabu-search-based heuristics outperform almost all previously known heuristics. Hence, we believe that a direct comparison between the results of the LR-PSS algorithm and theirs eliminates the need to compare the results to earlier heuristics individually. In any computational comparison, one challenge is to be able to compare the times obtained on different computers. Table 1 shows the speed comparison of different processors that Yagiura et al. and we use, obtained using benchmark programs at www.specbench.org. SpecInt2000 values are reported to be approximately 10 to 11 times of those of SpecInt95. Using this, we estimate that our processor is approximately 6 times faster than Yagiura et al.’s UltraSprac II 300 MHz processor. This is also consistent with the processor clock speeds as the ratio between the two provides a conservative estimate of their speeds. Table 2 presents the results of all benchmark instances. Each version of the LR-PSS algorithm is run 5 times for each instance, using 5 random number sequences for perturbation generations (called replications). The table lists the best upper bound (UBb) found in 5 replications, the best upper bound (UBn) without PSS component in the search, average duality gap (Gapa) percentage difference between the lower bound and the upper bounds (not shown), and finally (Gapbn) the percentage gap between the upper bounds obtained with and without PSS, all reported for each version of the algorithm. These results are obtained by running each version of the algorithm for a maximum 1500 subgradient iterations (subject to early termination due to small step parameter or no improvement on upper bound for a fixed number of iterations) with 10 perturbations in PSS within each subgradient iteration. (The tables focus on the results of the C, D, and E-type instances. We also solve A and Table 1 Comparisona of different processors Processor Ultrasparc II 300 MHz AMD Athlon 1.8 GHz a

SpecInt95 12.3 N/A

SpecInt2000 N/A 765

b

Estimate 1 6.23c

http://www.specbench.org. SpecInt95 values are not available for new processors Currently available values are SpecInt2000. c Estimates that AMD Athlon 1.8 GHz is 6.23 times faster. b

1047

B-type instances, where the average duality gaps are very small in general and there is little difference between the results of the LR-PSS heuristics and those of the Tabu search heuristics that are compared. These results are not shown in the tables for brevity. We also solve 60, mostly small, maximization instances listed by Chu and Beasley (1997), and the LR-PSS algorithm (using the version RGT) finds the optimal solutions for all 60 instances in very small computation times usually within a second. Those are also not presented here for brevity). It is clear from the results that PSS can enhance the performance of LR significantly. In case of larger instances the percentage gap between the UBs obtained using with and without PSS could be as high as 67%. On average, this effect is magnified in larger instances. We also see in case of LBR, inclusion of PSS degraded the performance of LR a bit for a few instances, but overall PSS truly enhances the search and improves the UB values, regardless of the feasibility restoration heuristic, instance size, or type. All versions of the algorithm (for a fixed set of parameters) are rather robust over the replications. Although there is no magnified difference between the different versions of the algorithm as illustrated by the similar average duality gaps, RGT performs overall slightly better than the other two, and HBR is a close second. Two exceptions to this general behavior are that (1) LBR performs slightly better for small D instances (up to 200 jobs), and (2) HBR is the best for small C instances. Potentially the most interesting result from this table is that all versions of the algorithm produce higher-quality solutions for larger instances (400 or more jobs), as evidenced by extremely small average duality gaps, mostly within 0.5% of the lower bound. Note that these gaps are calculated as average of gaps between respective upper and lower bounds (both over 5 replications); hence, the duality gaps between the best UB and best LB are even smaller. This is especially true for C and D-type large instances. The performance degrades a bit for E-type instances. Recall that heuristics typically perform better for smaller instances (in fact the tabu search heuristics by Dı´az and Ferna´ndez (2001) and Yagiura et al. (2004) are great examples as they consistently produce very high-quality results for small instances). The LRPSS algorithm, especially with RGT and HBR versions, seems to excel in larger instances producing near-optimal results for larger instances. Moreover,

Problem

Algorithm-LBR b

UB

1048

Table 2 Results of all instances using the three versions of the LR heuristic with and without PSS Algorithm-RGT n

UB

a

Gap

bn

Gap

b

UB

Algorithm-HBR n

UB

a

Gap

bn

Gap

UBb

UBn

Gapa

Gapbn

100 100 100 200 200 200

5 10 20 5 10 20

1940 1408 1268 3465 2825 2408

1942 1414 1287 3465 2827 2417

0.98 2.10 4.43 0.48 1.11 1.64

0.10 0.43 1.50 0.00 0.07 0.37

1956 1418 1260 3484 2819 2405

1968 1425 1278 3500 2823 2420

1.80 2.52 3.99 1.02 0.94 1.34

0.61 0.49 1.43 0.46 0.14 0.62

1941 1410 1267 3466 2823 2406

1968 1425 1278 3500 2823 2420

1.07 2.18 4.48 0.48 1.04 1.42

1.39 1.06 0.87 0.98 0.00 0.58

Avg. d d d d d d

100 100 100 200 200 200

5 10 20 5 10 20

6369 6400 6315 12,769 12,491 12,403

6381 6417 6349 12,766 12,491 12,437

1.79 0.49 1.50 3.14 0.30 0.65 1.62

0.41 0.19 0.27 0.54 0.02 0.00 0.27

6376 6396 6337 12,779 12,493 12,394

6427 6436 6366 12,808 12,543 12,466

1.93 0.59 1.32 3.46 0.40 0.74 1.55

0.63 0.80 0.63 0.46 0.23 0.40 0.58

6378 6403 6343 12,782 12,502 12,410

6427 6436 6366 12,808 12,543 12,466

1.78 0.58 1.59 3.36 0.37 0.74 1.65

0.81 0.77 0.52 0.36 0.20 0.33 0.45

Avg. e e e e e e

100 100 100 200 200 200

5 10 20 5 10 20

12,761 11,720 8857 24,978 23,523 22,813

12,788 11,767 8865 24,967 23,505 22,928

1.28 1.10 1.77 6.23 0.26 1.12 2.31

0.21 0.21 0.40 0.09 0.04 0.08 0.50

12,783 11,736 8756 24,969 23,512 22,764

12,959 11,835 8863 25,124 23,617 23,137

1.34 1.42 1.81 5.25 0.20 1.08 2.01

0.52 1.38 0.84 1.22 0.62 0.45 1.64

12,723 11,712 8828 24,978 23,490 22,796

12,959 11,835 8863 25,124 23,617 23,137

1.38 1.00 1.77 5.99 0.32 1.06 2.17

0.44 1.85 1.05 0.40 0.58 0.54 1.50

Avg. c c c c c c c c c

400 400 400 900 900 900 1600 1600 1600

10 20 40 15 30 60 20 40 80

5616 4818 4306 11,372 10,046 9412 18,846 17,206 16,359

5609 4831 4369 11,406 10,817 11,959 19,297 23,084 21,396

2.13 0.53 1.06 1.95 0.38 0.83 1.20 0.27 0.53 0.57

0.18 0.12 0.27 1.46 0.30 7.67 27.06 2.39 34.16 30.79

5610 4801 4266 11,366 10,007 9363 18,820 17,176 16,329

5631 4907 4668 14,515 12,919 12,058 29,199 26,096 21,791

1.96 0.42 0.64 0.97 0.28 0.36 0.58 0.15 0.25 0.33

1.02 0.37 2.21 9.42 27.71 29.10 28.78 55.15 51.93 33.45

5611 4824 4275 11,364 10,024 9380 18,844 17,185 16,344

5631 4907 4668 14,515 12,919 12,058 29,199 26,096 21,791

2.05 0.48 1.07 1.17 0.36 0.49 0.70 0.29 0.27 0.37

0.99 0.36 1.72 9.19 27.73 28.88 28.55 54.95 51.85 33.33

Avg. d d d d d d d

400 400 400 900 900 900 1600

10 20 40 15 30 60 20

25,038 24,743 24,689 55,540 55,123 54,944 98,009

25025 24,779 26,524 55,665 61,955 75,719 100,816

0.81 0.38 0.85 1.50 0.27 0.58 0.81 0.20

11.55 0.05 0.15 7.43 0.23 12.39 37.81 2.86

25,025 24,713 24,590 55,496 55,015 54,857 97,972

25,082 25,323 28,541 59,937 69,581 80,283 156,603

0.44 0.34 0.73 1.05 0.20 0.36 0.59 0.18

26.46 0.23 2.47 16.07 8.00 26.48 46.35 59.84

25,041 24,753 24,677 55,521 55,022 54,871 97,996

25,082 25,323 28,541 59,937 69,581 80,283 156,603

0.58 0.38 0.86 1.42 0.24 0.41 0.66 0.19

26.28 0.16 2.30 15.66 7.95 26.46 46.31 59.81

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

c c c c c c

38.55 1.73

30.81 1.89 22.09 23.43 33.59 42.47 39.73 67.44 55.55 60.75 0.56 0.65 1.35 2.92 0.42 1.23 2.74 0.54 2.73 2.97 46,874 55,486 56,421 137,357 144,505 143,165 303,802 282,582 289,792

1.47

38.82

46,006 45,447 45,712 102,817 101,429 102,458 181,437 181,663 180,277 30.91 1.96 22.43 23.77 33.64 42.63 39.81 67.63 56.19 61.35 0.47 0.59 1.19 2.47 0.38 1.01 2.73 0.39 1.71 2.73 46,874 55,486 56,421 137,357 144,505 143,165 303,802 282,582 289,792

n

bn

b

a

21.06 1.87 Avg.

Average of percentage duality gaps between LB and UB over 5 replications. Best UB over 5 replications with PSS component. UB without PSS component. Percentage gap UBb and UBn.

45,974 45,321 45,586 102,781 10,1316 102,401 181,236 180,917 179,602 13.61 0.07 4.34 15.67 7.39 15.09 41.08 22.61 44.00 39.29 0.62 0.68 1.43 2.97 0.43 1.41 3.23 0.78 2.54 3.35 46,038 47,380 52,795 110,388 116,936 144,422 222,643 260,935 252,879 46,007 45,410 45,644 102,791 101,603 102,371 181,592 181,211 181,550 400 400 400 900 900 900 1600 1600 1600 Avg. e e e e e e e e e

10 20 40 15 30 60 20 40 80

1600 1600 d d

40 80

97,490 97,568

112,596 142,602

0.45 0.56

15.49 46.16

97,387 97,391

161,003 149,455

0.37 0.39

65.32 53.46

97,416 97,453

161,003 149,455

0.38 0.51

65.27 53.36

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

1049

near-optimality is proven as the algorithm produces tight lower bounds as natural part of the process, which is a desirable feature in a heuristic algorithm. Another interesting behavior observed from the results in Table 2 is that the performances of all three versions deteriorate as the number of machines increases for a given number of jobs for all three types of instances. Increasing duality gaps as the number of machines is raised are a sign of this behavior. We think that this observation and the previous one (improved performance for instances with more jobs) are linked: The relative importance of the machines (and the machine capacity constraints) is magnified with instances with a small number of jobs. As we relax the machine capacity constraints as part of the LR scheme used in the algorithm, this is somewhat expected. In other words, the performance of the LR-PSS algorithm will improve as more jobs are added for a given number of machines. Table 3 shows the stability/robustness of the LR-PSS algorithm using three statistics of the duality gaps over 5 replications results of which are shown in Table 2. Data in the table suggest that duality gap does not drastically improves or deteriorates significantly when a different seed is used for generate a sequence of random number stream used in the algorithm. The maximum range with which duality gap could vary is within 1%. Table 4 compares the results of the LR-PSS heuristic with those from Dı´az and Ferna´ndez (2001) (denoted by DF), and from Yagiura et al. (2004) (denoted by YIG). For LR-PSS, we show the RGT and HBR versions, to emphasize the time difference between them. All comparisons are made for the best reported upper bounds (DF best values represent the best of 30 replications, YIG 5 replications). The percentage differences are computed between the specific version of the LR-PSS heuristic and the DF heuristic (called DF% in the table) and the YIG heuristic (called YIG%). The results show that both versions of the LR-PSS heuristic are quite competitive, especially for C and D-type instances, producing solutions usually within 1% of those of DF and YIG for these instances. The main challenge for LR-PSS heuristic is E-type instances. The RGT version seems to improve HBR significantly for E-type instances, especially for large instances. In fact, RGT produces 0.22% improved solutions for large type D instances as compared to the results of DF. The improvement becomes practically insignificant when HBR is compared

1050

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

Table 3 Robustness over replications using different random seed Problem

Algorithm-LBR

Algorithm-RGT

Algorithm-HBR

Mina

Maxb

Rangec

Min

Max

Range

Min

Max

Range

c c c c c c

5 10 20 5 10 20

100 100 100 200 200 200

0.83 1.51 4.02 0.41 1.06 1.31

1.09 2.38 4.84 0.62 1.17 2.02

0.26 0.87 0.82 0.20 0.11 0.71

1.66 2.23 3.37 0.96 0.84 1.18

1.87 2.81 4.43 1.25 0.99 1.56

0.21 0.58 1.07 0.29 0.14 0.38

0.88 1.66 3.94 0.44 0.99 1.22

1.30 2.67 5.01 0.50 1.13 1.65

0.42 1.01 1.07 0.06 0.14 0.42

Avg. d d d d d d

5 10 20 5 10 20

100 100 100 200 200 200

0.37 1.21 2.81 0.26 0.58 1.52

0.56 1.69 3.28 0.34 0.79 1.73

0.49 0.19 0.47 0.47 0.08 0.20 0.21

0.48 1.15 3.17 0.34 0.60 1.44

0.69 1.50 3.62 0.51 0.87 1.66

0.44 0.20 0.35 0.46 0.17 0.27 0.21

0.51 1.26 3.26 0.36 0.67 1.57

0.69 1.83 3.49 0.40 0.79 1.75

0.52 0.17 0.57 0.23 0.04 0.12 0.17

Avg. e e e e e e

5 10 20 5 10 20

100 100 100 200 200 200

0.95 1.53 5.99 0.22 0.98 2.06

1.30 2.02 6.57 0.34 1.23 2.51

0.27 0.36 0.49 0.58 0.12 0.24 0.46

1.12 1.67 4.77 0.19 0.94 1.83

1.50 2.08 5.65 0.26 1.19 2.42

0.28 0.38 0.41 0.89 0.08 0.26 0.59

0.65 1.46 5.62 0.22 0.84 1.97

1.30 2.08 6.29 0.39 1.40 2.44

0.22 0.66 0.62 0.67 0.17 0.55 0.47

Avg. c c c c c c c c c

10 20 40 15 30 60 20 40 80

400 400 400 900 900 900 1600 1600 1600

0.45 0.92 1.75 0.31 0.72 1.04 0.25 0.39 0.47

0.61 1.21 2.10 0.49 0.90 1.26 0.30 0.62 0.61

0.37 0.16 0.29 0.35 0.18 0.18 0.21 0.05 0.23 0.15

0.34 0.56 0.80 0.26 0.32 0.51 0.11 0.21 0.28

0.52 0.69 1.13 0.30 0.41 0.63 0.19 0.33 0.36

0.43 0.18 0.13 0.33 0.04 0.09 0.12 0.08 0.12 0.08

0.36 1.04 1.02 0.24 0.49 0.70 0.24 0.27 0.37

0.53 1.17 1.30 0.44 0.49 0.70 0.32 0.27 0.37

0.52 0.18 0.13 0.28 0.20 0.00 0.00 0.07 0.00 0.00

Avg. d d d d d d d d d

10 20 40 15 30 60 20 40 80

400 400 400 900 900 900 1600 1600 1600

0.33 0.78 1.40 0.25 0.54 0.72 0.19 0.40 0.55

0.44 0.91 1.61 0.29 0.60 0.86 0.21 0.48 0.57

0.20 0.11 0.13 0.21 0.04 0.07 0.14 0.02 0.08 0.02

0.28 0.65 1.00 0.17 0.34 0.56 0.15 0.29 0.37

0.39 0.84 1.11 0.23 0.37 0.63 0.19 0.43 0.41

0.13 0.11 0.18 0.11 0.06 0.03 0.07 0.04 0.14 0.05

0.34 0.82 1.35 0.22 0.35 0.59 0.18 0.32 0.43

0.42 0.88 1.51 0.26 0.45 0.72 0.21 0.43 0.56

0.10 0.08 0.07 0.16 0.04 0.10 0.13 0.03 0.11 0.12

Avg. e e e e e e e e e

10 20 40 15 30 60 20 40 80

400 400 400 900 900 900 1600 1600 1600

0.59 1.24 2.54 0.37 1.19 2.30 0.53 1.67 2.77

0.78 1.53 3.32 0.49 1.61 4.07 0.94 3.44 4.15

0.09 0.19 0.29 0.78 0.13 0.42 1.76 0.41 1.77 1.38

0.51 1.03 2.40 0.36 0.90 2.34 0.33 1.51 1.63

0.62 1.31 2.65 0.43 1.23 3.10 0.49 2.14 3.59

0.09 0.11 0.29 0.25 0.08 0.33 0.77 0.16 0.64 1.96

0.58 1.32 2.67 0.39 1.02 2.39 0.44 1.95 2.02

0.76 1.39 3.10 0.44 1.45 3.01 0.61 3.16 3.61

0.09 0.18 0.07 0.42 0.05 0.42 0.62 0.16 1.21 1.60

Avg. a b c

Minimum percentage duality GAP over 5 replications. Maximum percentage duality GAP over 5 replications. Range of percentage duality GAP over 5 replications.

0.79

0.51

0.53

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

1051

Table 4 Comparison with DF and YIG heuristics in terms of solution quality Problem

Algorithm-RGT UBb

Algorithm-HBR

DF%a

YIG%c

UBb

Time to best (TTB)

DF%a

YIG%c

RGT

HBR

HBR/RGTd

c c c c c c

100 100 100 200 200 200

5 10 20 5 10 20

1956 1418 1260 3484 2819 2405

1.29 1.14 1.37 0.81 0.46 0.59

1.29 1.14 1.37 0.78 0.43 0.59

1941 1410 1267 3466 2823 2406

0.52 0.57 1.93 0.29 0.61 0.63

0.52 0.57 1.93 0.26 0.57 0.63

1.2 1.7 3.2 4.1 6.5 16.8

1.1 2.0 2.5 3.4 3.0 9.3

0.91 1.21 0.79 0.82 0.46 0.56

Avg. d d d d d d

100 100 100 200 200 200

5 10 20 5 10 20

6376 6396 6337 12,779 12,493 12,394

0.94 0.36 0.70 2.04 0.28 0.41 0.94

0.93 0.30 0.65 1.88 0.25 0.29 0.35

6378 6403 6343 12,782 12,502 12,410

0.76 0.39 0.81 2.13 0.30 0.49 1.07

0.75 0.33 0.76 1.98 0.27 0.36 0.48

1.1 3.4 7.1 3.3 9.4 18.2

1.2 1.6 5.9 1.9 6.4 11.1

0.79 1.01 0.48 0.82 0.57 0.68 0.61

Avg. e e e e e e

100 100 100 200 200 200

5 10 20 5 10 20

12,783 11,736 8756 24,969 23,512 22,764

0.79 0.80 1.37 3.76 0.16 0.88 1.72

0.62 0.80 1.34 3.50 0.15 0.83 1.53

12723 11712 8828 24978 23490 22796

0.87 0.33 1.17 4.62 0.19 0.79 1.86

0.70 0.33 1.13 4.35 0.19 0.74 1.67

1.0 3.0 8.4 4.8 15.3 33.1

1.1 2.6 9.3 3.1 7.9 22.2

0.69 1.10 0.88 1.11 0.63 0.52 0.67

Avg. C C C C C C C C C

10 20 40 15 30 60 20 40 80

400 400 400 900 900 900 1600 1600 1600

5610 4801 4266 11,366 10,007 9363 18,820 17,176 16,329

1.45 0.21 0.31 0.42 N/A N/A N/A N/A N/A N/A

1.36 0.23 0.39 0.50 0.23 0.22 0.36 0.09 0.17 0.22

5611 4824 4275 11,364 10,024 9380 18,844 17,185 16,344

1.49 0.23 0.79 0.64 N/A N/A N/A N/A N/A N/A

1.40 0.25 0.87 0.72 0.21 0.39 0.55 0.22 0.22 0.32

28.5 57.1 137.2 293.8 649.4 1220.8 1352.3 2646.0 4084.5

13.0 14.7 59.7 61.1 129.4 380.9 237.8 564.0 812.1

0.82 0.46 0.26 0.44 0.21 0.20 0.31 0.18 0.21 0.20

Avg. D D D D D D D D D

10 20 40 15 30 60 20 40 80

400 400 400 900 900 900 1600 1600 1600

25,025 24,713 24,590 55,496 55,015 54,857 97,972 97,387 97,391

0.32 0.06 0.14 0.47 N/A N/A N/A N/A N/A N/A

0.27 0.19 0.42 0.53 0.11 0.19 0.35 0.10 0.23 0.30

25,041 24,753 24,677 55,521 55,022 54,871 97,996 97,416 97,453

0.55 0.01 0.02 0.12 N/A N/A N/A N/A N/A N/A

0.42 0.26 0.59 0.88 0.16 0.21 0.37 0.13 0.26 0.36

41.6 86.0 299.5 420.6 868.0 1486.9 1632.8 2600.3 5562.8

17.1 32.8 103.0 68.7 241.0 442.4 259.6 613.0 1836.7

0.27 0.41 0.38 0.34 0.16 0.28 0.30 0.16 0.24 0.33

Avg. E E E E E E E E E

10 20 40 15 30 60 20 40 80

400 400 400 900 900 900 1600 1600 1600

45,974 45,321 45,586 102,781 101,316 102,401 181,236 180,917 179,602

0.22 0.42 0.70 1.48 N/A N/A N/A N/A N/A N/A

0.27 0.50 0.97 2.25 0.35 0.87 2.22 0.33 1.46 1.55

46,006 45,447 45,712 102,817 101,429 102,458 181,437 181,663 180,277

0.03 0.49 0.98 1.76 N/A N/A N/A N/A N/A N/A

0.36 0.57 1.26 2.53 0.38 0.98 2.27 0.44 1.88 1.93

60.5 154.2 426.7 640.2 1765.1 3346.8 3470.8 6451.7 15192.9

19.2 61.6 209.4 149.5 434.5 1301.5 600.4 1650.1 6133.1

0.29 0.32 0.40 0.49 0.23 0.25 0.39 0.17 0.26 0.40

0.87

1.17

1.08

1.36

Avg. a b c d

Percentage difference between LR-PSS and TS heuristic of Dı´az and Ferna´ndez (2001). Best UB over 5 runs. Percentage difference between LR-PSS and TS heuristic of Yagiura et al. (2004) (Table 9, best UB). The ratio of the average time to best solution (TTB) of HBR and the average TTB of RGT (over 5 replications).

0.32

1052

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

to DF. Note that DF results are not available for instances larger than 400 jobs (denoted as N/A in the table). Again, we observe the improved performance of the LR-PSS heuristics from small instances (up to 200 jobs) to larger instances (400 jobs or more), as depicted by very small deviations from DF and YIG. Considering that YIG consistently produces the best solution values for all instances considered, and DF is not available for larger instances, this is a rather remarkable result, providing evidence that the LR-PSS could be an alternative with its simplicity and flexibility. The table also shows how much the HBR version improves RGT in terms of computation time. On average, HBR achieves reasonable results within a fraction of time needed to reach the best solution with RGT. More specifically, HBR finds its best about 70–80% of the time needed using RGT for small instances, and about 30% of the time needed with RGT for larger instances. If time is the main concern, and ultimate improvement in the solution value is not essential, one can use HBR without much significant cost. The main downside of using HBR is experienced with type E instances. For example, HBR produces a solution that is over 4.35% away from YIG solution (20 machines and 100 jobs), while RGT is able obtain a solution within 3.5%. Even for type E instances, and all other types, both versions of the LR-PSS heuristic are very competitive, improving relative performance for instances with small number of machines and large number of jobs. In comparison, DF seems to show opposite behavior (see Dı´az and Ferna´ndez, 2001), with no evidence for truly large instances (400 or more jobs), and YIG seems to be the most robust across different types and sizes of instances among the three. Tables 2 and 4 compare the best solution values obtained with the normal completions of the algorithms. These results compare the best case scenarios for all heuristics. To make a more fair comparison among the competing algorithms in terms of their convergence to a good solution, we run the HBR version of the LR-PSS algorithm with time limits that are equivalent to the time limits used in the YIG study. (Here we set the number of perturbations p to 5 instead of 10, so that the limited time is fairly distributed between LR and PSS search.) Table 5 summarizes the results of this time-limited HBR experiment along with those of Yagiura et al. (2004). (The YIG values are average values obtained from Table 8 of Yagiura et al.

(2004), which lists the average of the best upper bounds found up to the corresponding time limit. The time limits shown in the table are already adjusted for the difference between the computers according to the SpecInt discussion in the previous section.) Hence, we run the LR-PSS HBR algorithm for exactly the same amount of time (denoted by ‘‘Time Limit’’ (TL) in the table, adjusted by a factor of 6 for using faster computers). We observe that YIG is very robust, finding highquality solutions consistently. The PSS-LR algorithm with HBR find extremely close solutions to those found by YIG in almost all instances, except large-scale E-type instances. More interestingly, HBR seems to be capable of finding the corresponding best upper bound solutions very quickly, within the TL. The additional improvements from a full run of HBR (instead of a time-limited one) are found only for small C and E-type instances, for which the additional improvement seems to be miniscule (0.02% on average for both instance classes). In some cases, results from limited-time run are even better than those reported in Table 2 since we changed the p value. Hence, if the time is a concern HBR provides an efficient way of obtaining high-quality solutions very quickly for large C and D-type instances. For large instances, RGT seems to be effective (as compared to HBR), especially for E-type instances, but its computation time is not as competitive as that of HBR. Once again, YIG is extremely efficient in finding very high-quality solutions. In fact, most of their solution values are very close to lower bounds; although not proven, they are likely to be optimal. However, the LR-PSS algorithm provides reasonable results, close to the YIG values, within the same computational limits. In addition to being competitive in terms of solution quality and time, the LR-PSS algorithm presents certain other advantages: (1) It produces a Lagrangian relaxation lower bound as part of the algorithm, which can be used to assess the near-optimality of the heuristic solutions. (2) The LR-PSS algorithm combines simple ideas. Its overall procedure is even simpler than the Dı´az and Ferna´ndez (2001) TS algorithm, its main advantage is pronounced to be its simplicity and flexibility. Yagiura et al. (2004) combine sophisticated TS techniques and complex neighborhood structures. In fact, Cordeau et al. (2002) cite simplicity and flexibility (along

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

1053

Table 5 Comparison with YIG results (limited time runs) Problem

Upper bounds HBR

% UB difference YIG

Time

HBR-TL%a

RGT-all%b

HBR-all%c

TLd

C C C C C C

5 10 20 5 10 20

100 100 100 200 200 200

1940 1408 1260 3463 2820 2409

1931.0 1402.0 1243.0 3456.0 2806.2 2391.6

0.47 0.43 1.37 0.20 0.49 0.73

1.29 1.14 1.37 0.78 0.43 0.59

0.52 0.57 1.93 0.26 0.57 0.63

25 25 25 50 50 50

Avg. D D D D D D

5 10 20 5 10 20

100 100 100 200 200 200

6369 6397 6322 12,757 12,487 12,368

6355.6 6359.6 6220.0 12745.6 12445.4 12284.4

0.61 0.21 0.59 1.64 0.09 0.33 0.68

0.93 0.30 0.65 1.88 0.25 0.29 0.35

0.75 0.33 0.76 1.98 0.27 0.36 0.48

25 25 25 50 50 50

Avg. E E E E E E

5 10 20 5 10 20

100 100 100 200 200 200

12,721 11,701 8839 24,967 23,484 22,764

12681.4 11577.0 8440.6 24930.0 23308.0 22384.0

0.59 0.31 1.07 4.72 0.15 0.76 1.70

0.62 0.80 1.34 3.50 0.15 0.83 1.53

0.70 0.33 1.13 4.35 0.19 0.74 1.67

25 25 25 50 50 50

Avg. C C C C C C C C C

10 20 40 15 30 60 20 40 80

400 400 400 900 900 900 1600 1600 1600

5615 4824 4275 11,372 10,014 9380 18,833 17,185 16,344

5597.0 4783.0 4245.0 11341.4 9986.2 9330.4 18804.0 17148.6 16299.4

1.45 0.32 0.86 0.71 0.27 0.28 0.53 0.15 0.21 0.27

1.36 0.23 0.39 0.50 0.23 0.22 0.36 0.09 0.17 0.22

1.40 0.25 0.87 0.72 0.21 0.39 0.55 0.22 0.22 0.32

50 50 50 166 166 166 833 833 833

Avg. D D D D D D D D D

10 20 40 15 30 60 20 40 80

400 400 400 900 900 900 1600 1600 1600

25,027 24,756 24,696 55,504 55,067 55,967 97,990 97,435 97,598

24978.2 24609.2 24466.6 55436.2 54908.8 54677.6 97875.4 97166.0 97105.2

0.40 0.20 0.60 0.94 0.12 0.29 2.36 0.12 0.28 0.51

0.27 0.19 0.42 0.53 0.11 0.19 0.35 0.10 0.23 0.30

0.42 0.26 0.59 0.88 0.16 0.21 0.37 0.13 0.26 0.36

50 50 50 166 166 166 833 833 833

Avg. E E E E E E E E E

10 20 40 15 30 60 20 40 80

400 400 400 900 900 900 1600 1600 1600

46,002 45,415 46,755 102,807 102,935 109,353 181,452 183,980 190,525

45747.2 44889.6 44608.0 102425.0 100445.4 100207.6 180649.0 178319.4 176889.8

0.60 0.56 1.17 4.81 0.37 2.48 9.13 0.44 3.17 7.71

0.27 0.50 0.97 2.25 0.35 0.87 2.22 0.33 1.46 1.55

0.36 0.57 1.26 2.53 0.38 0.98 2.27 0.44 1.88 1.93

50 50 50 166 166 166 833 833 833

3.32

1.17

1.36

Avg. a b c d

Percentage difference between UB of the limited-time run of HBR and Yagiura et al. (2004) (Table 8, Ave. UB). Percentage difference between UB of the full-run of HBR and that of YIG (from Table 4). Percentage difference between UB of the full-run of RGT and that of YIG (from Table 4). Time limit. UBs represent the best reported upper bound found up to the time limit (seconds).

1054

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

with solution quality and speed) as dimensions to be considered when comparing vehicle routing heuristics. More generally, Silver (2004) argues simplicity as a necessary means for ‘‘facilitation of implementation’’ of heuristics. Similar arguments apply here for the GAP as complex algorithms are likely to be rarely implemented, while simple, intuitive algorithms have a better chance to be coded and used. Considering that these algorithms will be used as part of larger schemes to solve more complex problems that contain GAP as subproblems, the simplicity gains even a larger role. Moreover, the simplicity of the LR-PSS algorithm is further magnified as it does not need an optimization package or solver (i.e., one can program the algorithm in any programming language without any LP or IP solver). Finally, the simplicity of the algorithm is strengthened as it has few parameters to adjust. Along with standard subgradient optimization parameters, the algorithm has two main parameters: (i) the amount of perturbation (b), and (ii) the number of perturbations (p). While the literature has significant amount of work suggesting different ways to set and adjust the subgradient optimization parameters, the preliminary experiments are usually enough to set the LR-PSS perturbation parameters. Once we set these parameters, there is no need to adjust them for different problem sizes or types. However, we may need to adjust them when we have only limited time to solve the problem, so that the limited time is distributed in a favorable manner. (3) The LR-PSS algorithm is flexible and generic. It is flexible, in the sense that it can handle additional constraints (to be relaxed along with the machine capacity constraints to keep the simple structure of the relaxed problem). Potentially more importantly, the overall mechanism can be applied to other problems, especially as a PSS enhancement to regular Lagrangian relaxation. Potential application areas are scheduling and logistics problems that naturally come with hard constraints. 9. Conclusions and future research directions We develop and test a heuristic that combines ideas from Lagrangian relaxation, subgradient optimization, and problem space search for the general-

ized assignment problem. We instantiate versions of this generic heuristic using two specific methods to restore feasibility and find high-quality solutions over the iterations of subgradient optimization. Our extensive tests conducted using the most challenging GAP instances in the literature show that the PSS-enhanced Lagrangian heuristic within subgradient optimization is effective finding highquality solutions to small, and especially large instances. Comparing the two methods of feasibility restoration, the greedy regret minimization approach works the best, although this has increased computational burden if applied to all perturbed instances generated in every iteration for larger instances. Hence, we introduce a hybrid heuristic that combines the speed of lower bound replication version for the perturbed weights and the greedy regret version applied for the original weights. This seems to produce a good compromise between solution quality and computation time. Finally, our extensive comparison between the PSS-LR algorithm and the latest, and arguably the best heuristics in the literature (Dı´az and Ferna´ndez, 2001; Yagiura et al., 2004) shows that the heuristic is in general comparable both in terms of solution quality and computation time. We believe that the overall benefit of the LR-PSS heuristic is its flexibility, generality, and simplicity. It is simple, as it can be implemented without a linear programming or integer programming solver, and it makes use of the rather simple iterative subgradient optimization framework. It is flexible, as one can plug in potentially any weight-based GAP heuristic in the same general framework. It can also handle additional constraints. It is also generic, mainly because one can apply the ideas to other hard optimization problem. We outline the following future research directions: • The heuristic can be further improved through a few strategies, such as initializing Lagrange multipliers more intelligently and using enhancements on the traditional subgradient optimization procedure. The former should quickly improve the lower bounds and then hopefully upper bounds, and the latter, an example of which is conditional subgradient optimization, can use the Lagrange multiplier values in last few iterations instead of just the last one, (refer to Camerini et al., 1975). Moreover, other feasibility restoration algorithms can be used instead of the tested ones, LBR, and

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

RGT. During the investigation of these modifications, one needs to make sure that they improve the current heuristic across most, if not all, of the problems. Our rather limited tests for some of these modifications show that they work for some problems but not all, which warrants additional research. • One can investigate ways to make the approach generic enough to be used in conjunction with the already generic Lagrangian relaxation technique, that can be applied to other combinatorial optimization problems. As PSS is successfully applied to many other problems, including hard scheduling problems and others, this can lead to a robust, flexible and generic heuristic approach to solve other hard problems. • We can also improve the PSS module by appending with a more intelligent search strategy such as genetic algorithm or simulated annealing than a constant-size neighborhood search implemented here. This way, one can intensify and diversify search (and its perturbations) depending on the solution quality and overall search strategy. • Finally, one can test the overall heuristic in other problem domains, possibly more challenging domains than the GAP, such as scheduling and network design. One should expect to improve the traditional subgradient optimization. Acknowledgements We acknowledge the detailed and insightful comments of three anonymous referees on earlier versions of this paper. References Amini, M.M., Racer, M., 1994. A rigorous comparison of alternative solution methods for the generalized assignment problem. Management Science 40, 868–890. Baker, B.M., Sheasby, J.E., 1999. Extensions to the generalized assignment heuristic for vehicle routing. European Journal of Operational Research 119, 147–157. Balachandran, V., 1976. An integer generalized transportation model for optimal job assignment in computer networks. Operations Research 24 (4), 742–759. Camerini, P.M., Fratta, L., Maffioli, F., 1975. On improving relaxation methods by modified gradient techniques. Mathematical Programming 3, 26–34. Cario, M.C., Clifford, J.J., Hill, R.R., Yang, J., Yang, K., Reilly, C.H., 2002. An investigation of the relationship between problem characteristics and algorithm performance: A case study of the generalized assignment problems. IIE Transactions on Operations Engineering 34 (3), 297–312.

1055

Cattrysse, D.G., Wassenhove, L.N.V., 1992. A survey of algorithms for the generalized assignment problem. European Journal of Operational Research 60, 260–272. Cattrysse, D.G., Salomon, M., Wassenhove, L.N.V., 1994. A set partitioning heuristic for the generalized assignment problem. European Journal of Operational Research 72, 167–174. Charon, I., Hudry, O., 1993. The noising method: A new method for combinatorial optimization. Operations Research Letters 14, 133–137. Chu, P.C., Beasley, J.E., 1997. A genetic algorithm for the generalized assignment problem. Computers and Operations Research 24, 17–23. Available from: . Codenotti, B., Manzini, G., Margara, L., Resta, G., 1996. Perturbation: An efficient technique for the solution of very large instances of the Euclidean TSP. INFORMS Journal on Computing 8, 125–133. Cordeau, J.-F., Gendreau, M., Laporte, G., Potvin, Y., Semet, F., 2002. A guide to vehicle routing heuristics. Journal of the Operational Research Society 53, 512–522. Dı´az, J.A., Ferna´ndez, E., 2001. A tabu search heuristic for the generalized assignment problem. European Journal of Operational Research 132, 22–38. Farias, I.R. De, Nemhauser, G.L., 2001. A family of inequalities for the generalized assignment polytope. Operations Research Letters 29, 49–51. Fisher, M.L., 1981. The Lagrangian relaxation method for solving integer programming problems. Management Science 27, 1–18. Fisher, M.L., Jaikumar, R., 1981. A generalized assignment heuristic for vehicle routing. Networks 11, 109–124. Fisher, M.L., Jaikumar, R., Wassenhove, L.N.V., 1986. A multiplier adjustment method for the generalized assignment problem. Management Science 32, 1095–1103. Gottlieb, E.S., Rao, M.R., 1990. The generalized assignment problem: Valid inequalities and facets. Mathematical Programming 46, 31–52. Haddadi, S., 1999. Lagrangian decomposition based heuristic for the generalized assignment problem. INFOR 37, 392–402. Haddadi, S., Ouzia, H., 2001. An effective Lagrangian heuristic for the generalized assignment problem. INFOR 39, 351–356. Haddadi, S., Ouzia, H., 2004. Effective algorithm and heuristic for the generalized assignment problem. European Journal of Operational Research 153, 184–190. Higgins, A.J., 2001. A dynamic tabu search for large-scale generalized assignment problems. Computers and Operations Research 28 (10), 1039–1048. Laguna, M., Kelly, J.P., Gonza´lez Velarde, J.L., Glover, F., 1995. Tabu search for the multilevel generalized assignment problem. European Journal of Operational Research 82, 176–189. Martello, S., Toth, P., 1981. An algorithm for the generalized assignment problem. In: Brans, J.P. (Ed.), Operational Research, IFORS. North-Holland, Amsterdam, pp. 589–603. Naphade, K.S., Wu, D., Storer, R.H., 1997. Problem space search algorithms for resource constrained project scheduling. The Annals of Operations Research 70, 307–326. Narciso, M.G., Lorena, L.A.N., 1999. Lagrangian/surrogate relaxation for generalized assignment problems. European Journal of Operational Research 114 (1), 165–177. Nauss, R.M., 2003. Solving the generalized assignment problem: An optimizing and heuristic approach. INFORMS Journal of Computing 15 (3), 249–266.

1056

V. Jeet, E. Kutanoglu / European Journal of Operational Research 182 (2007) 1039–1056

Nauss, R.M., 2006. The generalized assignment problem. In: Karlof, J.K. (Ed.), Integer Programming: Theory and Practice. CRC Press, Boca Raton, FL, pp. 39–55. Osman, I.H., 1995. Heuristics for the generalized assignment problem: Simulated annealing and tabu search approaches. OR Spektrum 17, 211–225. Reinfeld, N.V., Vogel, W.R., 1958. Mathematical Programming. Prentice-Hall, Englewood Cliffs, NJ. Romeijn, H.E., Romero-Morales, D.M., 2000. A class of greedy algorithms for the generalized assignment problem. Discrete Applied Mathematics 103, 209–235. Ross, G.T., Soland, R.M., 1975. A branch and bound algorithm for the generalized assignment problem. Mathematical Programming 8, 91–103. Ross, G.T., Soland, R.M., 1977. Modeling facility location problems as generalized assignment problems. Management Science 24 (3), 345–357. Savelsbergh, M., 1997. A branch-and-price algorithm for the generalized assignment problem. Operations Research 45 (6), 831–841. Silver, E.A., 2004. An overview of heuristic solution methods. Journal of the Operational Research Society 55 (9), 936– 956. Storer, R.H., Wu, S.D., Vaccari, R., 1992. New search spaces for sequencing problems with application to job shop scheduling. Management Science 38, 1495–1509.

Storer, R.H., Wu, S.D., Vaccari, R., 1995. Problem and heuristic space search strategies for job shop scheduling. ORSA Journal on Computing 7, 458–487. Storer, R.H., Flanders, S.W., Wu, S.D., 1996. Problem space local search for number partitioning. The Annals of Operations Research 63, 465–487. Trick, M.A., 1992. A linear relaxation heuristic for the generalized assignment problem. Naval Research Logistics 39, 137– 152. Wilson, J.M., 1997. A simple dual algorithm for the generalized assignment problem. Journal of Heuristics 2 (4), 303–311. Yagiura, M., Ibaraki, T., in press. Generalized assignment problem. In: Gonzalez, T.F. (Ed.), Handbook of Approximation Algorithms and Metaheuristics, Chapman and Hall/ CRC in the Computer and Information Science Series. Yagiura, M., Yamaguchi, T., Ibaraki, T., 1998. A variable depth search algorithm with branching search for the generalized assignment problem. Optimization Methods and Software 10, 419–441. Yagiura, M., Ibaraki, T., Glover, F., 2004. An ejection chain approach for the generalized assignment problem. INFORMS Journal of Computing 16 (2), 133–151. Yagiura, M., Ibaraki, T., Glover, F., 2006. A path relinking approach with ejection chains for the generalized assignment problem. European Journal of Operational Research 169, 548–569.