A COMPARATIVE EVALUATION OF EXACT AND HEURISTIC METHODS FOR TRANSFER LINE BALANCING PROBLEM

A COMPARATIVE EVALUATION OF EXACT AND HEURISTIC METHODS FOR TRANSFER LINE BALANCING PROBLEM

A COMPARATIVE EVALUATION OF EXACT AND HEURISTIC METHODS FOR TRANSFER LINE BALANCING PROBLEM Olga Guschinskaya and Alexandre Dolgui Division for Indust...

164KB Sizes 0 Downloads 98 Views

A COMPARATIVE EVALUATION OF EXACT AND HEURISTIC METHODS FOR TRANSFER LINE BALANCING PROBLEM Olga Guschinskaya and Alexandre Dolgui Division for Industrial Engineering and Computer Sciences Ecole des Mines de Saint Etienne 158, Cours Fauriel, 42023 Saint Etienne Cedex 2, France [email protected], [email protected]

Abstract: Ten best exact and heuristic methods for Transfer Line Balancing Problem (TLBP) are compared in a computational experiment. The TLBP deals with the optimization of design solutions for serial machining lines. Such lines consist of a sequence of unit head machines. The operations are grouped into blocks at every station. These blocks are executed sequentially and the operations within each block are performed simultaneously by the same multi-spindle head. The objective is to assign the operations to blocks and the blocks to stations minimizing the total number of station and spindle heads. The challenge is to minimize the line cost and time for line design. Experimental results are presented. They help to choose the best optimization method for each situation. Copyright © 2006 IFAC Keywords: Machining transfer lines, Line balancing, Optimization, Graph theory, Mixed integer programming, Heuristics, Random search.

machines are: high productivity, reduced space and cost.

1. INTRODUCTION The considered problem is referred as Transfer Line Balancing Problem (TLBP) taking into account its proximity to the Assembly Line Balancing Problem (ALBP), see e.g. Rekiek et al. (2002). )

)

)

)

)

)

)

)

)

)

)

0

)

The TLBP appears before a serial machining line is put into production (new line is designed or existing line is reconfigured). Finding a good (and if possible the best) solution reduces line costs and space requirements. A correctly configured line offers a competitive advantage to a company (Dashchenko, 2003). The considered serial machining lines consist of a number of linearly ordered machines that perform the operations required to manufacture parts. All the machines are linked by an automated material handling device without buffers between machines. Each machine can be equipped by one or several multi-spindle heads activated sequentially by changing the active spindle head or by moving the part to successive spindle heads. Each spindle head may have several tools that perform the operations simultaneously, so the spindle head processing time is equal to the maximal operation time of operations assigned to it. The advantages of multi-spindle

There is at least one multi-spindle head at each machine. Each spindle head incurs an investment cost and requires a certain working space at a machine, therefore the total number of spindle heads at a machine is limited. All spindle heads of the same machine are activated sequentially, so the machine processing time is equal to the sum of spindle head processing times. The objective is to minimize the total number of spindle heads and machines (i.e. the line cost which is composed of fixed cost per machine and additional costs for spindle heads). This objective can be reached by an optimal assignment the operations to spindle heads and to machines. The solution must provide a desired productivity; as well as satisfy all technological constraints. In this paper, 10 best exact and heuristic methods for TLBP are compared in a computational experiment. Experimental results are presented. They help to choose the best optimization method for each situation.

The paper is organized as follows. Section 2 deals with the input data. Section 3 presents investigated methods. Section 4 is dedicated to experimental results. Concluding remarks are given in Section 5. 2. INPUT DATA As with standard terminology, the following notions are used: “station” for “machine” as well as for a set of operations assigned to machine and “block” for “spindle head” as well as for a set of operations assigned to a spindle head. The following input data assumed to be known (Dolgui et al., 2005b): • N is the set of all operations needed for a part machining; • tj is the processing time of the operation j∈N. • T0 is the maximal admissible line cycle time (desired productivity); • C1 is the relative cost of one station; • C2 is the relative cost of one block; • m0 is the maximal authorized number of stations; • n0 is the maximal number of blocks per station. The precedence constraints between the operations define non-strict partial order relation over operation set N. They can be represented by digraph G = (N, D). The arc (i, j) ∈ N×N belongs to set D if and only if operation j cannot precede operation i. Some groups of operations must be assigned to the same machine, because of a required machining tolerance. Such constraints are represented by ES , the family of subsets from N, where all operations of each subset must be assigned to the same station. The incompatibilities between groups of operations interdict the assignment of such groups of operations to the same spindle head or machine. These constraints are modeled by the following families of subsets of N: • ES , where the operations of the same subset cannot be assigned to the same station; • EB , where the operations of the same subset cannot belong to the same block.

K(j) is the set of indices k (stations) where operation j can be assigned; e is a set of operations which is an element of ES , ES or EB ; j(e) is some fixed operation from the set e. Variables: Xjq is a binary decision variable (1 if operation j is assigned to block q and 0 otherwise); Fq ≥ 0 is an auxiliary variable for determining the time of block q; Yq ∈ {0, 1} is an auxiliary variable that indicates if the block q exists or not; Zk ∈ {0, 1} is an auxiliary variable that indicates if the station k exists or not. The variables Yq and Zk are used to count the number of blocks and stations, respectively.

The precedence constraints may be formulated in four different ways. In this paper, each creates a separate MIP model, so we have four different MIP models to be investigated: MIP0 with constrains (2), MIP1 with constrains (2'), MIP2 with constrains (2''), MIP3 with constrains (2''') and (2''''). Constraints (3)(4) deal with the impossibility of grouping certain operations in one block or execution certain operations at the same station, respectively. Constraints (5) determine the necessity of grouping certain operations in the same station. Constraints (6) reflect the fact that each operation must be assigned to exactly one block. Constraints (7) determine the block processing times. Constraint (8) is the constraint on the cycle time. Constraints (9) ensure that the block q exists in the design decision if and only if Xjq =1 for any j. Constraints (10) ensure that the station k exists in the design decision if and only if Yq =1 for any q∈S(k). Constraints (11) guarantee that block q is created in station k only if block q-1 exists for this station. Constraints (12) ensure that station k can be created only if station k-1 exists at the line. q0

k =1

q =1

subject to: ∑

∑ X iq ' ≥ X jq Pred ( j )

(2)

X is ≥ X jq

(2')

i∈Pr ed ( j ) q '≤ q



s ≤ q, s∈Q (i )

3. METHODS INVESTIGATED



3.1 Mixed Integer Program

q∈Q ( j )

qX jq ≥



In (Dolgui et al., 2000), a mixed integer program (MIP) was suggested.

s ≤ q, s∈Q (i )



s∈Q (i )

(2'')

sX is

X is ≥ X jq , q∈Q(j)

qX jq ≥

(2''')

sX is

(2'''')

∑ X jq ≤ e − 1 , e∈ EB

(3)



Model notations: Pred(j) is the set of direct predecessors of j∈N; q is the index of block, e.g., for the block l of the station k, q=(k-1)n0 +l; q0 is the maximal possible value of q, q0=m0n0; S(k)={(k-1)n0 +1, …, kn0} is the set of block indices for the station k; Q(j) is the set of indices q (blocks) where operation j can be assigned;

m0

Min W(Zk,Yq)= C 1 ∑ Z k + C2 ∑ Yq , (1)

q∈Q ( j )



s∈Q (i )

j∈e





j∈e q∈S ( k )





j∈e \{ j (e)} q∈S ( k )

X

X jq ≤ e − 1 , e∈ ES jq

=

(

e −1 )

e ∈ ES



q∈S ( k )

(4)

X

j (e ) q

,

(5)

q0

∑ X jq = 1

(6)

Fq ≥ tj Xjq

(7)

q =1



q∈S ( k )

Fq ≤ T 0

(8)

Yq ≥Xjq

(9)

Zk ≥Yq, q=(k-1) n0+1

(10)

Yq −1 − Yq ≥ 0 , q=(k-1) n0+1

(11)

Z k −1 − Z k ≥ 0 , k =2, …, m0,

(12)

where j ∈ N, i∈Pred(j), q = 1,2, …, q0 , k =1,…, m0 . To reduce the number of decision variables and constraints, the intervals of possible values of block and station numbers for each operation are obtained from the analysis of all the constraints (for details, see Dolgui et al., 2000 and Dolgui et al., 2005b).

decision P(x) = {{u1\u0, ..., u j1 \u j1 −1 }, { u j1 +1 \u j1 , …, u j2 \u j2 −1 }, …, { u jm ( x ) −1 +1 \u jm ( x ) −1 , …,

u jm \u jm ( x ) −1 }}. Thus, the initial problem can be reduced to the following constrained shortest parameterized path problem:

Min W(x)=C1m(x)+C2 l(x)

(12)

x∈X

(13)

t (ui \ ui −1 ) ≤T0

(14)

j

∑i =r j

r −1 +1

e ⊄ ( u j r \ u j r − 1 ),

e∈

ES

(15)

jr - jr-1 ≤ n0, j0=0

(16)

m(x)≤ m0,

(17)

where r=1,...,m(x). The algorithm developed in (Dolgui et al., 2000) finds a shortest constraint path in the above digraph.

3.2 Shortest path approach In (Dolgui et al., 2000), the following exact method was suggested. A feasible design decision can be represented by the collection P={{N11, …, N1 n1 } , …, {Nm1, …,

Nm nm }}, where m - number of stations; nk - number of blocks for station k; Nkl (k=1, 2, …, m; l=1, 2, …, nk) - assignment of the set N of operations to stations and blocks. Let P be a set of collections P. The set vkl = ∪ kr =−11 ∪ qr=1 N rq ∪ ∪lq=1 N kq is the set of operations n

executed on the part after machining it by l-th block of the k-th station. Let V be the set of all vkl for all P∈P, including v0 = ∅, vN=N. For each vertex v∈V a parameter Γ(v) is introduced. Γ(v) = 1 if for all e∈ES either e ∩ v = ∅ or e ⊆ v, otherwise Γ(v) = 0. The value 1 means that there are no operations in v linked with operations not included in v by a constraint ES. A digraph H = (V, F) is constructed in which the arc (v′, v″) ∈ F if and only if v′ ⊂ v″, and e ⊄ N″ = v″\v′ for all e∈ ESB = EB ∪ ES , i.e. the arc (v′, v″) represents the block v″\v′ of operations. For each arc (v′, v″), we assign a time parameter t(v′, v″) = max{tj | j∈N"}. To each design decision P∈P we can associate a parameterized path x(P) = ((u0 = v0, ..., uj-1, uj, ..., ul(x) = vN), (γ 0 = 1, ..., γ j, ..., γ l(x) = 1)) in the digraph H from the vertex v0 to the vertex vN. The parameter γ j is equal to 1 if the block uj\uj-1 is the last block of the corresponding station of P. Let X be the set of all parameterized paths x in the digraph H from v0 to vN where γ j takes the value 0 if Γ(uj )=0 and 0 or 1 otherwise. In the parameterized path x∈X, there is a sequence j0 = 0, j1, j2, ..., jm(x) = l(x) of indices jr with γ jr = 1 for r = 0, ..., m(x). Then, for each parameterized path x∈X, a design

3.3 FSIC heuristic The random search algorithm named FSIC (First Satisfy Inclusion Constraints) was suggested in (Dolgui et al., 2005a). It is based on the COMSOAL method (Arcus, 1966), i.e. several operation assignments are randomly generated, and the best assignment is kept. In FSIC algorithm, all the operations located by the constraints at one station (and their non-assigned predecessors) are processed first. Let list L1 be a list of operations that can be assigned to the current station. Let Op be the randomly chosen operation from the list L1. The algorithm uses an additional list L2. This list groups all the operations that must be assigned to the same station as operation Op, i.e., the list is made of operation Op and all the operations from each set e∈ES including operation Op and all their non-assigned predecessors. The algorithm tries to assign all the operations from the list L2 to the current station. If impossible, the algorithm then creates a new station and tries to assign all the operations from the list L2 to the new station. The algorithm stores the state of the current station before the first assignment attempt of the list L2. It restores its state before the second try if the list L2 contains more than one operation. When all the operations from L2 have been assigned, L1 is rebuilt. An iteration of this algorithm gives a feasible solution for the transfer line (or the conclusion that the solution of this iteration is not feasible). The algorithm makes TRtot iterations and returns a best assignment in the end. 3.4 Decomposition methods The concept of decomposing the initial problem to several sub-problems on basis of the graph of precedence was suggested in (Dolgui et al., 2005b). An approximate solution can be found by

partitioning the initial set N of operations into several disjoint subsets (N1, N2…, Np) and a step by step solution of the obtained problems with N1, N2…, Nu…, Np. The subset Nu is chosen in a way that: Pred(j)⊆Nu∪Nu-1 for any j∈Nu, either e∩Nu=∅ or e⊆Nu for any e∈ ES , where Nu-1= ∪ur=−11 N r . In the model (1)-(12), the set N is replaced by the set Nu; each set e containing operations from N\Nu is removed from the sets ES, EB , and ES ; a new order (precedence) relation is defined over the set N′(u-1) and introduced into digraph Gr. At the u-th step, the operations of the current subset Nu are assigned to blocks and stations. By varying both the partition of the set N and the solution chosen at each step, various results can be obtained, the best is kept. 3.5 Hybrid method Another method based on the decomposition of initial problem was suggested in (Guschinskaya et al., 2005c). At first, the initial problem is solved by FSIC heuristic to obtain a feasible solution Pheur, i.e. an assignment of all operations to mheur stations. This solution is a sequence Sheur of stations, the hybrid approach is based on the decomposition of the sequence Sheur into several subsequences in such a way that: • Each subsequence r includes a number of consecutive stations; • There are no intersections between subsequences; • The union of all subsequences is the initial sequence Sheur. Each subsequence r of mr stations is used to generate a new sub-problem SPr which is the same type as the initial problem but with a reduced size, where: • The set of operations N is reduced to Nr ⊂ N, where Nr is the set of operations assigned to the mr stations of the subsequence r; • Constraints are transformed by replacing the set N with the subset Nr, and then removing the constraints related to operations from N\ Nr. • m0 is replaced by mr. Each sub-problem SPr is then solved by shortest path approach. Since, the optimal solution of the subproblem Cr_opt is found, the algorithm replaces the old sequence of mr stations by the obtained optimal sequence. The resolution of all sub-problems SPi, i = 1, .., knum obtains a new feasible solution Phybr that consists of mhybr stations. Then the heuristic is used to find a new feasible solution that will be decomposed and so on until the end of available calculation time. 3.6 Aggregate solving Each sub-problem defined by Nu for the deterministic decomposition or by Nr for the hybrid

method can be solved independently (as shown in sections 3.4 and 3.5) or considering the fact that all operations of the subset Nu-1 (Nr-1) have been already assigned. This last resolution is called here aggregate solving, because the results of previously solved subproblems are included as macro-operations in the resolution of a consecutive sub-problem. Let Nu(q) be a set of operations of Nu assigned to the block q after u steps, N′(u) = {Nu(q) ≠ ∅ | 1 ≤ q ≤q0}. Each subset N u-1(q)∈N′(u-1) is treated as a macrooperation, i.e.: - Digraph Gr is modified (all arcs (i, j) are deleted from Gr if i, j∈N u-1(q); all arcs (i, j) are replaced by one arc for operations i∈N u-1(q); multiple arcs are replaced by one arc); - Each operation of subsets ES, EB , and ES is replaced by the macro-operation which contains this operation; - In the obtained family EB , all subsets which contain other subsets from this family are deleted, idem for the family ES . Both independent solving and aggregate solving for deterministic decomposition and hybrid method are included in the evaluation. 4. COMPARATIVE RESULTS 4.1 Experimental Conditions In order to compare the 10 optimization procedures described above 5 data sets were generated, each of which contains 50 transfer line balancing problems. The data sets differ by two parameters: the number of operations |N| and the order strength OS of the precedence graph. They represent wide range of possible transfer line balancing problems. Therefore, the performances of studied methods are investigated on the problems of different difficulty. Note: the order strength is defined as the density of the transitive closure of the precedence graph (Scholl and Klein, 1998). Each unique transfer line balancing problem was generated at random. The operation times are derived from the uniform distribution (a=10, b=20). The precedence graph is generated at random by adding arcs between randomly chosen vertexes until getting desired order strength. All data sets were solved by the following algorithms: 1) SP, shortest path approach; 2) M0, mixed integer approach, model MIP0; 3) M1, mixed integer approach, model MIP1; 4) M2, mixed integer approach, model MIP2; 5) M3, mixed integer approach, model MIP3; 6) FSIC heuristics; 7) Deterministic decomposition algorithm with independent solving of sub-problems; 8) Deterministic decomposition algorithm with aggregate solving of sub-problems;

9) Hybrid method with independent solving of sub-problems; 10) HY, hybrid method with aggregate solving of sub-problems. But even for the simpler data set, the results obtained by deterministic decomposition algorithm both with independent and aggregate solving were not competitive. Thus they were excluded from further experimentation, leaving 8 optimization methods to be evaluated. Furthermore, the hybrid method with independent solving was eliminated because for all data sets it was surpassed by the method with aggregate solving. So, only the results obtained by the remaining 7 methods are pertinent. Experiments were carried out on Pentium IV, 3GHZ. The obtained results are reported in Tables 1-5. For the presentation, the following notations are used: NS – number of solved tests of data set; NO – number of obtained optimums; ∆ max , ∆ av , ∆ min – maximal, average and minimal deviation of the criterion obtained by a method from the best value, respectively; Tmax, Tav, Tmin – maximal, average and minimal solution time, respectively. Note: for heuristic algorithms solution time is not equal to the time of finding the best solution, it is equal to available computational time. Note: if an algorithm does not solve all tests the parameters ∆ max , ∆ av , Tmax , Tav are calculated only for solved instances. 4.2 Small size problems For the two first data sets, the available computational time was limited by 1800 sec for the exact methods and by 90 sec for the heuristics. The first data set contains 25-operations problems with OS= 50%. The results are reported in Table 1. Table 1 Results for 1st data set NS NO ∆ max

SP 50 50 0

M0 50 50 0

M1 50 49 4

M2 50 50 0

M3 50 49 4

FSIC 50 23 10.5

HY 50 39 5.26

∆ av

0

0

0.08

0

0.08

3.14

1

∆ min

0 11.4 1.4 0.03

0 1800 169 1.39

0 1800 164 1.2

0 403 57.6 0.16

0 1800 190 1.5

0 90 90 90

0 90 90 90

Tmax Tav Tmin

For this data set, the shortest path method is the best both in terms of computational time and the quality of provided solutions. The best MIP model is MIP2, because it was able to find all optimums but not as fast as shortest path method. The hybrid method found only 39 optimums but its performance is better than FSIC heuristics. The second data set contains 25-operations problems with OS= 15%. The results are in Table 2.

Table 2 Results for 2nd data set NS NO ∆ max

SP 50 50 0

M0 50 50 0

M1 50 49 4

M2 50 50 0

M3 50 49 4

FSIC 50 27 10.5

HY 50 35 11.1

∆ av

0

0

0,08

0

0.08

2,8

1.7

∆ min

0 1638 292 3.77

0 1044 147 0.72

0 1199 75.5 0.63

0 479 62 0.20

0 1800 286 1.0

0 90 90 90

0 90 90 90

Tmax Tav Tmin

Here, MIP2 method is the best in terms of computational time. Therefore, even for small size problems, the performance of the shortest path method is not as good as MIP for the problems with small value of OS. The hybrid method found only 35 optimums but its performance is better than FSIC heuristics. 4.3 Medium size problems Here, the available computational time was limited by 1800 sec for the exact methods and by 300 sec for the heuristics. The third data set contains 50-operations problems with OS= 90%. The results are reported in Table 3. For the third data set the shortest path method is the best both in terms of computation time and the quality of provided solutions. Between MIP methods, MIP2 method found the maximum of optimal solutions, but it was not able to solve one instance. MIP0 method solved all instances, but it found only 44 optimal solutions. The hybrid method found only 34 optimums but its performance is better than FSIC heuristics. The fourth data set contains 50-operations problems with OS= 45%. The results are reported in Table 4. Here, MIP algorithms found only some solutions (< 10), this is why they were excluded. The shortest path algorithm found only 14 optimum solutions. So, one can conclude that even for medium size problems with small value of OS the two heuristic methods are preferable if the computational time is limited. Between heuristics, the hybrid method found 47 optimums and its performance is much better than FSIC. Table 3 Results for 3d data set NS NO ∆ max

SP 50 50 0

M0 50 44 31.1

M1 44 41 31.25

M2 49 47 7.84

M3 46 42 16

FSIC 50 18 6

HY 50 28 4.2

∆ av

0

1.25

0.88

0.21

1

2.13

1

∆ min

0 0.09 0.04 0.01

0 1800 672 0.01

0 1800 569.5 0.03

0 1800 613.3 0.03

0 1800 982.7 0.01

0 300 300 300

0 300 300 300

Tmax Tav Tmin

4.4 Large size problems The fifth data set contains 100-operations problems with OS= 25%. The results are reported in Table 5. Table 4 Results for 4th data set

Table 5 Results for 5th data set

∆ max

SP 14 14 0

FSIC 50 28 6

HY 50 47 2.7

FSIC 50 31 16.2

∆ av

0

1.38

0.16

1.31

0

∆ min

0 694 1406 13

0 300 300 300

0 300 300 300

0 600 600 600

0 600 600 600

NS NO

Tmax Tav Tmin

HY 50 50 0

The available computational time was limited by 5700 sec for the exact methods and by 600 sec for the heuristics. Exact algorithms found only some solutions (< 10) for these data sets; this is why they were excluded. So, for large size problems heuristics are recommended. The hybrid method with aggregate solving is much better than FSIC heuristics. In about 38% of instances it found a better solution.

heuristic approaches. The experiments show that the heuristics performances depend on the characteristics of the problem (order strength, constraints of compatibility). On the basis of our analysis the following general recommendations are formulated. For small size problems (less than 25 operations), exact methods should be used. If the order strength of precedence graph is relatively large (>40%) then shortest graph method is recommended, otherwise mixed integer programming method providing model M2 should be used. For medium size problems, for which the value of the order strength of precedence graph is relatively large (>60%) shortest graph method can be used, otherwise as for large size problems (more than 50 operations), hybrid method with aggregate solving of sub-problems is advised. Further investigations will concern improving the hybrid and deterministic decomposition methods as well as mixed integer linear programming models. The heuristic algorithms can be also used for finding “good” initial solutions and upper bounds for the value of an objective function. ACKNOWLEDGMENT

5. CONCLUSIONS An important industrial problem was considered: optimal configuration of machining lines. A line configured properly means competitive advantage later. The solution of this problem is also a challenge taking into account the exigency of reducing the time from design to manufacturing. The optimization of line configuration is very complex and it often requires a great computation time. The development of efficient methods and choice of a method adapted to each concrete problem is an important issue. Both exact and heuristic algorithms for the configuration a transfer line with sequentially activated spindle heads at stations are evaluated. The exact algorithms are based on shortest path and MIP approaches. The heuristic algorithms use different techniques: random search; deterministic decomposition and decomposition based on a feasible solution with independent and aggregate solving of sub-problems. The algorithms have been tested on 5 data sets each of which contained 50 randomly generated problems. Four MIP models have been presented for using with MIP solvers in order to find an optimal design decision. The experiments show that MIP2 model is more preferable from the point of view of the quality of the obtained decisions and the running time. So far, it is possible to conclude that hybrid method with aggregate solving of sub-problems is the best heuristic method from the point of view of the quality of the obtained decisions. This conclusion is based on its experimental comparison with other

This work is financially supported by INTAS Project 03-51-5501. REFERENCES Arcus, A.L. (1966) COMSOAL: A computer method of sequencing operations for assembly lines. International Journal of Production Research, 4, 259-277. Dashchenko A. I. (Ed) (2003) Manufacturing Technologies for Machines of the Future 21st Century Technologies, Springer. Dolgui, A., Guschinsky, N. and Levin, G. (2000). Approaches to balancing of transfer line with block of parallel operations, Preprint No. 8, 42 pages, Institute of Engineering Cybernetics/ University of Technology of Troyes, Minsk. Dolgui A., Finel B., Guschinsky N., Levin G., and Vernadat F. (2005a) A heuristic approach for transfer lines balancing. Journal of Intelligent Manufacturing, 16 (2), 159-171. Dolgui A., Finel B., Guschinsky N., Levin G., and Vernadat F. (2005b) MIP approach to balancing transfer lines with blocks of parallel operations. IIE Transactions, (accepted). Guschinskaya O., Dolgui A., Guschinsky N., Levin G., (2005c). Hybrid approach for optimization of a class of machining lines, European Journal of Operational Research, (submitted). Rekiek, B., A. Dolgui, A. Delchambre and A. Bratcu (2002). State of art of assembly lines design optimisation. Annual Reviews in Control, 26(2), 163-174 Scholl, A. and Klein, R. (1998) Balancing assembly lines effectively: a computational comparison. European Journal of Operational Research, 114, 51-60.