Bernardo Almada-Lobo
email: almada.lobo@fe.up.pt

Diego Klabjan

Maria Antónia Carravilla

Jose Fernando

José F Oliveira

Single Machine Multi-product Capacitated Lotsizing with Sequence-dependent Setups

Keywords: Sequence-dependent setup, capacitated lotsizing and scheduling capacitated lotsizing and scheduling, sequence-dependent setup, integer programming, heuristics

come   L'archive ouverte pluridisciplinaire

Production planning arising in the glass container industry is a complex process. It is a semicontinuous manufacturing process, which includes the glass production (the continuous part) and the container manufacturing (the discrete part). Furnaces are operated continuously and machine lines operate on a twenty four hour seven days a week basis. During a product changeover on a molding machine, the furnace keeps feeding the machine with molten glass, however, the gobs are discarded and melted down again in the furnace (wasting a huge amount of energy, the main production cost). Fast job changes with fewer job change parts lower the energy costs. Therefore, there is a sequence dependent setup time (and proportional cost) in a product changeover. Since machine setups consume the furnace capacity, it is essential to carry the setup from one period to the next. This production planning problem results in a capacitated multi-item deterministic lotsizing problem.

The lotsizing problem is to find production orders or lots in order to satisfy customer demand. The objective is to minimize holding and setup costs. This problem has been intensively studied and, therefore, many models were proposed to describe it, involving different features and assumptions. The majority of the models consider a finite planning horizon divided into discrete time periods. The single stage, multi-item lotsizing problem with capacity constraints is referred to as the capacitated lotsizing problem (CLSP). CLSP is considered to be a large-bucket problem, because several products/setups may be produced/performed per period. On the other hand, in small-bucket problems, at most one setup may be executed during a period. Different models are distinguished with respect to the number of products that can be produced within a period and to the period capacity utilization (all-or-nothing assumption or partial). All these models determine simultaneously the size and sequence of lots. Despite allowing more precise plans, sub-dividing the (macro-) periods of CLSP into several (micro-) periods increases the complexity of these mixed integer linear programming models. [START_REF] Drexl | Lot sizing and scheduling -survey and extensions[END_REF] present and explain the main differences of these formal models.

We study single machine CLSP with sequence dependent setup times and costs. In each time period several products can be produced. Whenever a switch of a product is performed, a certain amount of time is consumed and the setup cost is incurred. This time and cost depend on the two products. In each time period at most a given amount of time can be allotted to switchover times and to the actual production time. This leads to sequence dependent setup times and costs. A further complicating matter is in the fact that a product can be set up at the end of a time period and the actual production starts in the next time period. We call this property setup carryover. We present two novel linear mixed integer programming formulations for single machine CLSP with sequence-dependent setup times and costs, incorporating all the necessary features of setup carryovers. The compact formulation has a polynomial number of constraints. The efficiency of this model is benchmarked against the one proposed by [START_REF] Gupta | The capacitated lot-sizing and scheduling problem with sequencedependent setup costs and setup times[END_REF], which is actually not a completely accurate formulation [START_REF] Almada-Lobo | A note on "The capacitated lot-sizing and scheduling problem with sequence-dependent setup costs and setup times[END_REF]). The second formulation uses an exponential number of constraints, which makes it impractical to be solved directly. We generate constraints dynamically by presenting a polynomial time separation algorithm. We compare the strength of the linear programming relaxations of both formulations. The two formulations are improved by adding valid inequalities that lead to good lower bounds. Since this problem is NP-hard, we also present a five-step heuristic to solve large problem instances. The flexibility of the heuristic's first two steps allows us to find feasible solutions even under tightly capacitated scenarios. The last three steps focus on solution quality. We compare the heuristic solution against the lower bound generated by the strongest model and against other known heuristics. The main contributions of our work are as follows. To the best of our knowledge, the two presented formulations are the first exact formulations for the problem under consideration. We also formally prove a relationship between the two linear programming relaxations. Another major contribution of our work is the underlying heuristic.

In Section 2 we first present the compact formulation of the problem, which is also benchmarked against the model presented by [START_REF] Gupta | The capacitated lot-sizing and scheduling problem with sequencedependent setup costs and setup times[END_REF]. We then propose the stronger formulation and prove that it provides a better linear programming relaxation. In Section 3 we explain the heuristic for finding good solutions. Numerical experiments are given in Section 4. We conclude this introduction with a literature review.

Literature Review

Standard CLSP does not sequence or schedule products within a period. In addition, it assumes that setup costs occur for each lot in a period, even if the last product to be produced in a period is the first one in the period that follows. For manufacturing environments with considerable setup times with respect to the length of a period and with tight capacity, which is the case in the glass container industry, this model may not provide any feasible production plans. Therefore, setup carryovers need to be incorporated in CLSP. The setup carryover occurs when a product is the last one to be produced in a period and the first one in the following period. In this case no setup is required in the latter period. Moreover, setups are preserved over idle periods and an idle time at the end of a period may be used to perform a setup for the first product to be produced at the beginning of the next period.

Despite CLSP being intensively studied, modeling setup carryover in CLSP has not received much attention due to model complexity and computational difficulty, [START_REF] Sox | The capacitated lot sizing problem with setup carry-over[END_REF]. The efficiency of solving CLSP depends on the structure of setups. Different studies have demonstrated that proper accounting for setup times and carryovers decreases the number of setups and also frees a significant amount of production capacity, [START_REF] Porkka | Multiperiod production planning carrying over set-up time[END_REF]. [START_REF] Gopalakrishnan | A framework for modeling setup carryover in the capacitated lot-sizing problem[END_REF] developed a modeling framework for formulating CLSP with sequence and product independent setup times and setup carryovers. In their model, a fixed charge is incurred whenever a product is produced in a period. The authors conclude that the model complexity could be reduced by looking for alternative ways to model setup carryover. In [START_REF] Gopalakrishnan | A modified framework for modelling set-up carryover in the capacitated lotsizing problem[END_REF], a modified framework for modeling setup carryover in CLSP is presented, incorporating sequence independent and product dependent setup times and costs. Nevertheless, setup carryover is modeled in a similar way as in [START_REF] Gopalakrishnan | A framework for modeling setup carryover in the capacitated lot-sizing problem[END_REF] and decision variables have the same interpretation. [START_REF] Sox | The capacitated lot sizing problem with setup carry-over[END_REF] present a more efficient mixed integer linear programming model to CLSP with sequence independent setup costs and no setup times. [START_REF] Kang | Lotsizing and scheduling on parallel machines with sequencedependent setup costs[END_REF] propose a different approach, consisting of breaking the entire schedule into smaller segments, called split-sequences. A drawback of this model is that the number of split-sequences L t in period t is a parameter. Therefore, multiple runs with different values for L t are needed for optimization. Moreover, setup times are not considered in this model. [START_REF] Porkka | Multiperiod production planning carrying over set-up time[END_REF] modify the mixed integer linear programming based setup carryover model of [START_REF] Sox | The capacitated lot sizing problem with setup carry-over[END_REF] by using sequence independent setup times. To avoid the complexity of separate setup costs and times, explicit setup costs are excluded from the model. [START_REF] Suerie | The capacitated lot-sizing problem with linked lot sizes[END_REF] present a new model formulation for CLSP with sequence independent setup costs and times. The authors start with the standard CLSP model, use the standard facility location reformulation, and introduce new sets of variables and constraints to model the setup carryover. [START_REF] Gopalakrishnan | A tabu-search heuristic for the capacitated lot-sizing problem with set-up carryover[END_REF] also address CLSP with setup carryover and present two effective tabu-search heuristics.

All the aforementioned manuscripts address CLSP with constant sequence independent setup times and/or setup costs, with a setup carryover. In fact, these models do not consider lot sequencing within a period. They focus only on determining the products produced last and first in two consecutive periods, and also the configuration of the machine at the end of the period. Due to the "correct" calculation of the setup costs by linking lots of adjacent periods, these models are usually referred to as CLSP with linked lot-sizes [START_REF] Haase | Lot sizing and scheduling for Production Planning[END_REF]).

Notwithstanding, there are several examples of process industries in which sequence-dependent setup times and costs are considerable, such as glass container and some chemical industries. [START_REF] Clark | Rolling-horizon lot-sizing when set-up times are sequence-dependent[END_REF] model CLSP with sequence-dependent setup times using a new mixed integer programming formulation. They assume that exactly a given number of setups occur in a time period between any two given products, independently of their demand patterns. [START_REF] Haase | Lot sizing and scheduling with sequence-dependent setup costs and times and efficient rescheduling opportunities[END_REF] propose a model for CLSP with sequence dependent setup times and costs in which the efficient product sequences are pre-determined. Therefore, it is to decide which sequences will be 3 used in each period. The authors also assume that the inventory of an item must be null at the beginning of a period to allow its production in the period. [START_REF] Gupta | The capacitated lot-sizing and scheduling problem with sequencedependent setup costs and setup times[END_REF] extended the framework proposed by [START_REF] Gopalakrishnan | A modified framework for modelling set-up carryover in the capacitated lotsizing problem[END_REF] to incorporate sequence dependent setup times and setup costs. The setup carryover is modeled in the same way as in [START_REF] Gopalakrishnan | A framework for modeling setup carryover in the capacitated lot-sizing problem[END_REF] and [START_REF] Gopalakrishnan | A modified framework for modelling set-up carryover in the capacitated lotsizing problem[END_REF].

2 New Models for CLSP with Sequence-dependent Setup Costs and Times

The First Formulation

Throughout the exposition, t denotes time periods, which range from 1 to T , and i and j index products, which are labeled from 1 to N . We denote by [M ] the set {1, 2, . . . , M }. We consider a general single-stage model involving multiple items to be scheduled on a single machine, whose data consist of: h i cost of carrying one unit of stock of product i from one period to the next, p i processing time of one unit of product i, d it demand for product i at the end of period t, C t capacity of the machine in period t (measured in time units), s ij time needed to set up the machine from product i to product j, c ij cost incurred to set up the machine from product i to product j.

In addition, let M it = min Ct p i , T u=t d iu be an upper bound on the quantity of product i to be produced in period t.

The total cost consists of the holding cost and the machine setup cost. We do not allow stockouts, which is common in the deterministic demand setting. We assume that the triangle inequality with respect to the setup cost and time holds, i.e., c ik ≤ c ij + c jk and s ik ≤ s ij + s jk for all products i, j, and k. This assumption holds in many practical settings, and definitely in our glass container case.

In order to capture the lotsizes and the resulting inventory we need the following decision variables: X it quantity of product i produced in period t, I it stock of product i at the end of period t.

In order to capture the setup cost and time, we introduce

T ijt =     
1 if a setup occurs on the machine configuration state from product i to product j in period t, 0 otherwise.

We assume that T iit = 0 for every product i. Recall that we allow the machine to be set up with product i at the end of a time period and then start the next time period with the same product but not occurring any extra setup cost. Since such a setting involves two consecutive time periods, a new decision variable needs to be introduced. Let We also need the auxiliary variables V it that assign product i in period t.

Our first formulation (F 1 ) for CLSP with sequence-dependent setup costs and times, and setup carryover is as follows:

min i j t c ij • T ijt + i t h i • I it (1)
I it =I i(t-1) + X it -d it i ∈ [N ], t ∈ [T ] (2) i p i • X it + i j t s ij • T ijt ≤C t t ∈ [T ]
(3)

X it ≤M it •   j T jit + α it   i ∈ [N ], t ∈ [T ] (4) i α it =1 t ∈ [T ]
(5)

α it + j T jit =α i(t+1) + j T ijt i ∈ [N ], t ∈ [T ] (6) V it + N • T ijt -(N -1) -N • α it ≤V jt i ∈ [N ], j ∈ [N ] \ {i}, t ∈ [T ] (7) (X it , I it , α it , V it ) ≥ 0, T ijt ∈ {0, 1}. (8)
The objective function (1) minimizes the sum of sequence-dependent setup costs and the holding cost. Constraints (2) represent the inventory balances and (3) ensure that the total production and setup time in each period does not exceed the available capacity. Requirement (4) guarantees that a product is produced only if the machine has been set up for it. Constraints ( 5)-(7) determine the sequence of products on the machine in each period and keep track of the machine configuration state by recording the product that a machine is ready to process (setup carryover information is thereby tracked). We next discuss these constraints further. Hereafter we refer to T jit as an input setup for product i in period t and to T ijt as an output setup for product i in period t.

The model does not explicitly have binary variables representing whether a product is produced or not in a certain period. Production of product i can occur in period t if the machine is set up for i at the beginning of t or if at least one input setup is performed for product i. These conditions are guaranteed by constraints (4).

Constraints ( 5) ensure that α it is one for exactly one product in a given period, i.e., that the machine is set up for exactly one product at the beginning of each period.

Constraints (6) ensure a balanced network flow of the machine configuration states and carry the setup configuration state of the machine into the next period. If no setup is performed in period t, which happens when T ijt = 0 for every i and j, configuration state of the machine is carried from period t -1 to period t + 1. The idle time of a period t may also be used to perform an empty setup for the first product to be produced at the beginning of the next period. If there is an input setup and no output one for product i in period t, it means this setup was the last one to be performed on the machine in period t and, consequently, the machine is configured for product i at the beginning of the next period, i.e., α i(t+1) = 1. On the other hand, if there is an output setup and no input one, it means that the machine is configured for product i at the beginning of period t, i.e., α it = 1.

In general, constraints (6) impose that flow in equals flow out. We note that constraints (2)-( 6) and ( 8) allow solutions with cycles that are not disjoint. However, such a scenario will never occur in an optimal solution since it is easy to see that the triangle inequality assumption implies that the in flow of every node is at most 1.

Let us define a digraph

G = ([N ], A = [N ] × [N ]
), where nodes correspond to products and an arc a = (i, j) corresponds to the setup from product i to product j. It follows that the set {(i, j)|T ijt = 1} corresponds to a collection of disjoint paths and cycles. Constraints ( 2)-( 6) and ( 8) admit solutions that have at most one path, but potentially many cycles, called also subtours, see Figure 1. The left configuration represents a solution that entails one path from product i 1 to product i 4 and 3 disconnected subtours. The right configuration illustrates a solution with a cycle that starts and ends in product i 1 and 2 disconnected subtours. In order to obtain a feasible solution, we must impose that we have only a single connected component linked to α's. It can be shown that constraints (7) achieve this by using auxiliary variables that value the machine state through any sequence. It is also not difficult to see that no feasible solution is excluded by (7). A similar modeling trick is known in the context of the traveling salesman problem, see, e.g., [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF].

We conclude that the underlaying subgraph associated with an optimal T is either a cycle or a path. Thus it looks like the one depicted in Figure 1 except that there are no disconnected subtours.

Note that the integrality condition of α it is not necessary. Let us assume that α i1 = 1. From (6) it follows that α j2 is an integer for every j. Together with (5) we obtain that α j2 is binary. This reasoning rolls over to the subsequent periods. If the machine is not set up for any product in the first time period, then the integrality requirements of α j1 's need to be added. As an alternative to avoid it, the machine can be set up for a virtual product at the beginning of the planning period. No setup times and costs would be incurred if a setup occurs between this product and any other one.

The presented model assumes the single-machine case. The extension of this model to multiple machines is straightforward, also considering other usual features, such as shortages, backlogging 6 Example 1. Consider the following small data set provided in Table 1 for a three-product, threeperiod problem.

Table 1: Data for the three product, three-period problem Setup times are s ij = 5 time units for i = j and zero otherwise. Moreover, C t = 100 for every t and p i = 1 for every i.

d it c ij t = 1 t = 2 t = 3 j = 1 j = 2 j = 3 h i i =
We solved the model by using a commercial solver. Figures 2 and3 show the optimal solution for this instance, when the machine is set up for product 3 at the beginning of the planning period, i.e., α 31 = 1. Note that there is idle time at the end of period 3. at the beginning of period 2. In period 2 only product 3 is produced and no setup is performed. Therefore, the setup state of the machine for product 3 is carried from period 1 to period 3. At the end of the planning horizon the machine is ready to produce product 2.

Benchmarking the First Formulation

In this section we compare our model F 1 against Gupta and Magnusson's model. The experiments were conducted on an IBM machine with a 3.2 GHz processor and 1GB of random access memory and CPLEX 9.0 as a mixed integer programming solver.

In Table 2 we conduct a benchmark on the complexity (in terms of the number of binary variables, continuous variables and constraints) of F 1 against the models presented in [START_REF] Clark | Rolling-horizon lot-sizing when set-up times are sequence-dependent[END_REF] and [START_REF] Gupta | The capacitated lot-sizing and scheduling problem with sequencedependent setup costs and setup times[END_REF]. Here S is the maximum allowable number of setups per period used in Clark and Clark's model.

F 1 N 2 T 4N T N T (2 + N ) + 2T Clark and Clark [2000] N 2 (ST -1) + N SN T N T (SN + 2S + 2) +1 -N Gupta and Magnusson [2005] N 2 T + 4N T 2N T N T (5 + 3N ) + 8T
As we can observe from the number of binary variables, F 1 requires less variables than the other two. In Example 1, F 1 has 27 binary variables, whilst the Clark and Clark's model has 75 (S has to be set to 3 to reach the same optimal solution) and the Gupta and Magnusson's model has 63 binary variables.

In order to computationally evaluate the models, we generated random instances based on the parameters used by [START_REF] Haase | Lot sizing and scheduling with sequence-dependent setup costs and times and efficient rescheduling opportunities[END_REF]. The machine is always set up for product 1 at the beginning of the planning horizon. The number of products N ranges from 3 to 10, the number of periods were 3, 4, 5, 7, 10, and 15, and the capacity utilization (Cut) defined as i d it /C t takes the value 0.6. Processing times (p i ) are unitary for all of the products. Demand (d it ) is selected based on the uniform distribution U [40, 60], holding costs (h i ) based on U [2, 10], and setup times (s ij ) based on U [5, 10] (the triangle inequalities hold). The setup costs c ij are proportional to the setup times, i.e., c ij = θ • s ij with θ = 50. Note that θ implicity defines a relationship between the setup and holding costs. In the case of the company producing glass containers that inspired this work, θ falls in the interval [50,75].

For each triplet N , T , and Cut, three different instances were generated and the average running times (in seconds) were calculated. Table 3 gives the computational times to obtain an optimal solution by model F 1 . An empty field means that at least one of the instances could not be solved optimally within a 3600 seconds limit.

A comparison between the efficiency of F 1 and the Gupta and Magnusson's model (denoted by GM ) for Cut = 0.6 is shown in Figure 4. We did not compare F 1 to the Clark and Clark's model since multiple runs with different values for S would be needed for reaching the same optimal solution. The vertical axis denoted by F 1 /GM represents the ratio of the respective computational times. It is clear that F 1 model is always more efficient than the GM model. Moreover, as the number of periods increases, the ratio F 1 /GM tends to 0, i.e., the performance of F 1 is significantly better with increasing T . It turns out that the number of products does not influence significantly the relative efficiency.

The Compact Formulation

In this section we propose an alternative equivalent formulation (F 2 ) to CLSP with sequence dependent setup times and setup costs.

Since the triangle inequality holds for both setup costs and setup times, constraints are always satisfied in an optimal solution. The following inequalities are an alternative formulation to constraints (7):

j T ijt ≤1 i ∈ [N ], t ∈ [T ] (9) j T jit ≤1 i ∈ [N ], t ∈ [T ] (10
ϑ it = α i(t+1) + j T ijt i ∈ [N ], t ∈ [T ] ( 11
) i∈S j∈S T ijt ≤ i∈S ϑ it -ϑ kt + α kt t ∈ [T ], k ∈ S, S ⊆ [N ]. ( 12
)
By substituting ( 11) in ( 12), we obtain the equivalent requirement

i∈S j ∈S T ijt + i∈S α i(t+1) ≥ j T jkt t ∈ [T ], k ∈ S, S ⊆ [N ]. ( 13
)
If there is a subtour on a subset S of nodes, the left-hand side of ( 12) equals to |S|. If α it = 0 for every i ∈ S, this constraint is violated since its right-hand side equals to |S| -1 for every product k in the cycle. Note that ϑ it = 1 for every i ∈ S.

Now we argue that a feasible solution is not excluded by (13). To this end, let (X, I, α, T ) be feasible to F 1 . It suffices to consider the case in which the right-hand side of (13) equals one. For a fixed S and k ∈ S, let h ∈ [N ] be such that T hkt = 1. We distinguish two cases: h ∈ S or h ∈ S. h ∈ S: If subset S contains the node set of the entire cycle or path corresponding to the solution, then i∈S α i(t+1) = 1 and ( 13) is satisfied. On the other hand, if this is not the case, then we have to consider three different scenarios.

(a) If the solution entails a cycle, then clearly i∈S j ∈S T ijt ≥ 1 since the cycle must "enter" and "leave" S.

(b) If the solution consists of a path that ends at node k, then α k(t+1) = 1 and, therefore, i∈S α i(t+1) = 1. (c) If the solution consists of a path that does not end at node k, then either the path ends at another node in S (and i∈S α i(t+1) = 1) or it ends at a node not in S and, in this case, i∈S j ∈S T ijt ≥ 1.

Each of these scenarios yields constraint (13).

h ∈ S: This case is similar to the above case.

The F 2 formulation reads min 9), ( 10),(13) (X it , I it , α it ) ≥ 0, T ijt ∈ {0, 1}.

i j t c ij • T ijt + i t h i • I it subject to (2) -(6), (
Let P F 1 and P F 2 denote the feasible sets of the linear relaxations of the formulations F 1 and F 2 , respectively. In the following theorem we show that F 2 is at least as strong as F 1 .

Theorem 1. We have P F 2 ⊆ proj(P F 1 ), where proj denotes the projection of P F 1 on the space of (X it , I it , T ijt , α it ). Proof. We show the statement by exhibiting the appropriate V s. They are defined as shortest path distances with respect to predefined weights. Since shortest path distances are well defined on networks with no negative weight cycles, we rely on (12) and appropriately selected weights to show this.

Let (X * it , I * it , T * ijt , α * it ) ∈ P F 2 be arbitrary. We find V * it such that (7) hold, i.e.,

(X * it , I * it , T * ijt , α * it , V * it ) ∈ P F 1 .
Consider time period t and let w ij = (N -1) + N • α * it -N • T * ijt for every i, j. Consider network G defined in Section 2.1, where each arc has weight w ij . Let us fix an arbitrary source node s and let dist i denote the length of the shortest path from the source node to node i. We show later that these distances are well defined. They satisfy the following optimality conditions:

dist i + w ij ≥ dist j (i, j) ∈ A(G).
Then, by setting V it = -dist i , we obtain

V it -(N -1) -N • α * it + N • T * ijt ≤ V jt . (14)
To ensure that V it are non-negative, we define

V * it = V it + max j,t
|V jt |. Clearly, from ( 14) it follows that (V * ijt , T * ijt ) satisfy ( 7). We conclude that (X * it , I * it , T * ijt , α * it , V * it ) ∈ P F 1 . In order to complete the proof, we need to show that dist i are well defined. This is equivalent to showing that there are no negative cost cycles with respect to w (see, e.g., [START_REF] Ahuja | Network Flows: Theory, Algorithms, and Applications[END_REF]). We need to show that (i,j)∈C w ij ≥ 0 for any cycle C. Let S be the node set of cycle C and k ∈ S. From (12) and nonnegativity of T * , we have

(i,j)∈C T * ijt ≤ i,j∈S T * ijt ≤ i∈S ϑ it -ϑ kt + α * kt . (15)
Since

(i,j)∈C w ij = (N -1) • |S| + N • i∈S α * it -N • (i,j)∈C T * ijt (16)
and by using (15), we obtain

(i,j)∈C w ij ≥ (N -1) • |S| + N • i∈S α * it -N • i∈S ϑ it + N • ϑ ik -N • α * ik . (17)
Replacing ( 11) into (17), yields

(i,j)∈C w ij ≥ (N -1) • |S| -N • i∈S\{k} j T * jit .
Taking into account (10), leads to Note that F 2 has an exponential number of constraints. To introduce the most violated (t, S, k) inequalities ( 13) we need to solve the separation problem. We next show how to solve this problem.

(i,j)∈C w ij ≥ (N -1) • |S| -N • (|S| -1) ≥ 0, since |S| ≤ N .
Let (X * it , I * it , T * ijt , α * it ) satisfy ( 2) -( 6), ( 9) and ( 10). Consider network G with capacity T * ijt on arc (i, j) in A. Next we define a new network G 0 by adding sink node z to G and arcs (i, z) for every product i with capacity α * i(t+1) . Finally, we solve the min (s, z) cut problem for every product s. The algorithm hinges on the fact that the value of an s -z cut S in G 0 equals to i∈S j ∈S T * ijt + i∈S α * i(t+1) . If there exists a k ∈ S such that j T * jkt > Z * with Z * being the value of the minimum s -z cut, then (t, S, k) inequality is violated.

Since finding a minimum cut in a directed graph can be performed in polynomial time, this separation problem is polynomial.

Valid Inequalities

In this section we present two classes of valid inequalities (cuts) that strengthen F 2 .

The first class of valid inequalities is not a family in the original space, but we need to consider a higher dimensional space. Let W t denote a binary variable that equals 0 if at least one setup is performed in period t and 1 otherwise. Then, the following sets of valid inequalities, denoted as (W t ) inequalities,

α it ≤ j T ijt + W t i ∈ [N ], t ∈ [T ] α it ≤ j T ji(t-1) + W t-1 i ∈ [N ], t ∈ [T ] 1 -W t ≤ kj T kjt t ∈ [T ] W t ≤1 - kj T kjt N t ∈ [T ]
are valid.

To argue that they are valid, it suffices to consider the case α it = 1. We distinguish the following four cases. α i(t-1) = 0, α i(t+1) = 0: In this case j T ji(t-1) ≥ 1, j T ijt ≥ 1, and W t = W t-1 = 0 and therefore the inequalities are satisfied.

α i(t-1) = 1, α i(t+1) = 0: In this case j T ijt ≥ 1 and W t = 0. If kj T kj(t-1)
≥ 1, then also j T ji(t-1) ≥ 1 (and W t-1 = 0) and they are valid. On the other hand, if kj T kj(t-1) = 0, then W t-1 = 1 and again they are valid. α i(t-1) = 0, α i(t+1) = 1: This case is similar to the above. α i(t-1) = 1, α i(t+1) = 1: It is easy to adapt the second case to this case.

The introduction of W t does not increase significantly the computational time since there are only T of them, and, furthermore, their integrality is implied. It is also helpful to use the well known (l, S) inequalities for uncapacitated single-item lotsizing [START_REF] Barany | Strong formulations for multi-item capacitated lot sizing[END_REF]). The (l, S) inequalities adjusted for the multi-item case of our problem are expressed as

d i1l ≤ t∈[l]\S d itl • j T jit + α it + t∈S X it l ≤ T, S ⊆ [l], (18)
where d ipg = g t=p d it . Since there are an exponential number of these constraints, we introduce them as cutting planes (instead of adding all of them a priori). To find the most violated (l, S) inequalities, we use the separation algorithm introduced by [START_REF] Barany | Strong formulations for multi-item capacitated lot sizing[END_REF], which requires O(N T 2 ) time steps.

A Heuristic for CLSP

In this section we propose a heuristic for finding good solutions to CLSP, i.e., upper bounds on the optimal value. The heuristic consists of five basic steps. The first two attempt to find an initial feasible solution, whilst the last three are geared toward improving the quality of the solution.

In the first step the machine is loaded in each period with production gross requirements. Such a solution typically violates the capacity requirements. In the second step the production is sequenced from period T to 1. Whenever the capacity constraint is violated, we perform a procedure to eliminate overtime, which is the amount of the capacity violation. The third step is geared toward eliminating setups by anticipating productions, and the forth step tries to reduce holding costs by postponing productions. Finally, the last step attempts to find a better coordination between periods, eliminating, whenever possible, empty setups.

The heuristic starts with a lot-for-lot forward pass, allocating to each period the demand for that period (X it = d it , for every i, t), without considering the capacity constraints.

The second step works backwards and performs in each period both sequencing and amending procedures. In a given time period, the sequencing procedure is basically a minmax algorithm; in each pass a new arc (u, v) is selected due to either the variability of the input setup times of its initial node u (∆I u ) or the variability of the output setup times of its terminal node v (∆O u ). The variability is measured as the difference between the two smallest setup time values. The node with the maximum variability among the head and tail nodes is chosen. Then, the other node of the arc is the argument that minimizes the setup time of the selected arc. Since feasibility is placed over optimality in this step, the sequencing procedure tries to minimize the sequence dependent setup times. This procedure is detailed in Appendix A, together with an example. The option of considering the difference between the smallest and the second smallest setup time is preferred over the standard greedy approach of basing the selection only on the smallest setup time to "hedge" against the possibility of making a bad decision and later having very little room for correction. In addition, computational experiments have revealed this to be a superior criterion. For t < T , the last product to be scheduled is the first product in period t + 1 (the sink node is specified before generating the path in Algorithm 1).

Let

O t = max 0, i p i • X it + ij s ij • T ijt -C t denote
the amount of overtime in period t. After the sequencing procedure, there can be time periods with positive overtime, i.e., with violated capacity. In the amending pass, we push overtime towards the previous period. While working on period t, we might create a new overtime at period t -1. In case O t > 0, in order to determine the product or set of products whose (partial or entire) lots will be moved backwards, the following 13

p i • X it ≥ O t and X i(t-1) > 0} and k = arg min i∈S h i . Move O t /p k of X kt to period t -1. (b) Let S = {i ∈ [N ] : 0 < p i • X it < O t
and X i(t-1) > 0}, and let f s it denote the minimum savings in setup time by not producing product i in period t. Formally,

f s it =        j T ijt • s ij if α it = 1 and α i(t+1) = 0 and X it > 0, min X kt >0 ,X jt >0 k =i, j =i [s ki + s ij -s kj ] if α it = 0 and α i(t+1)
= 0 and X it > 0, 0 otherwise. Note that α i(t+1) = 1 is not considered because it would modify the production sequence in downstream periods.

Let k = arg min i∈S (p i • X it + f s it : p i • X it + f s it ≥ O t ). If k is not defined, then k = arg max i∈S (p i • X it + f s it ). Move all X kt to period t -1 and return to (a) if O t > 0. (c) Let S = {i ∈ [N ] : p i • X it ≥ O t } and k = arg min i∈S h i . Move O t /p k of X kt to period t -1. (d) Let S = {0 < p i • X it < O t } and k = arg min i∈S (p i • X it + f s it : p i • X it + f s it ≥ O t ). If k is not defined, then k = arg max i∈S (p i • X it + f s it ).
Move all X kt to period t -1, and return to (a) if O t > 0.

We note that steps (c) and (d) are similar to (a) and (b), respectively, however, they lead to a higher capacity consumption in period t -1 since a new setup has to be performed. This is the reason for placing step (c) after step (b) as opposed to execute (a) and (c) simultaneously. Proceeding from the last period to the first period, the solution is either feasible or there is overtime in period one (and no feasible solution is generated). We point out that finding a feasible solution is NP-complete, Trigeiro et al. [1989]. If the first two steps do not find a feasible solution, then the entire heuristic fails.

The following steps try to improve the quality of the initial solution. Let φ = {(i, t) : X it > 0}, let f c it denote the minimum savings in setup costs by not producing product i in period t (its expression is identical to f s it with setup costs instead of setup times as one of the arguments), and let Ca t = max 0, C ti p i • X itij s ij • T ijt denote the available (idleness) capacity in period t. The third step of the heuristic is a backwards pass in time that seeks to avoid the cost and capacity consumption of a setup, by moving an entire lot of product i * from period t to the closest previous period t * f in which product i * is produced and there is enough idle capacity (Ca t * f > p i • X i * t ) to absorb the new lot, where

(i * , t * f ) = arg max (i,t f )∈φ,(i,t)∈φ t f <t,Cat f >p i •X it [f c it -h i • X it • (t -t f )].
Note that these operations may destroy feasibility which needs to be recovered at the end. The fourth step is a forward pass that seeks to reduce inventory holding costs, by shifting forward a fraction or an entire lot of product i from a source period t to a target period t f in which product i is also produced (X it f > 0). This move is similar to the one presented by [START_REF] Hindi | An effective heuristic for the CLSP with set-up times[END_REF], but with a different purpose. Our objective is optimality, while their is feasibility. The candidate products have currently a positive inventory level (I it > 0). In each iteration, the production quantity Q t i * t * f of product i * is carried over from period t into period t * f , where

(i * , t * f ) = arg max i∈[N ],t f ∈{1,...,T -1} h i • Q t it f • (t f -t) : X it f > 0 , Q t it f = min{X it , I it , I i(t+1) , . . . , I i(t f -1) , C t f }.
Finally, the last step looks for improvements in the links between adjacent periods, in a forward pass. The sequencing procedure of the second step may generate a schedule in which the machine is set up to produce product i at the beginning of period t + 1 (α i(t+1) = 1), without producing i in that period (X i(t+1) = 0). If there is a product j such that X jt > 0, X j(t+1) > 0, α j(t+1) = 0, j = i, then by assigning α j(t+1) = 1, α i(t+1) = 0 a setup is eliminated in period t + 1. Let S = {j ∈ [N ] : X jt > 0, X j(t+1) > 0, α jt = 0, α j(t+2) = 0} be the set of candidate products. The last two conditions in the definition of S restrict the current analysis to periods t and t + 1 (otherwise other periods would have to be rescheduled). To select product i * from S we compute for every j ∈ S variables ∆I j and ∆O j , in the same way as in the sequencing procedure, except that all quantities are computed with respect to c ij (see Appendix A). To compute ∆O j and ∆I j we only take into account products that are produced in periods t and t + 1, respectively. The selected product has the minimum aggregated variability i * = arg min j∈S (∆O j + ∆I j ). Then, both periods t and t + 1 are rescheduled with the input α i * (t+1) = 1.

From all the steps described before, only the third and the fifth steps may generate unfeasible schedules due to the creation of overtime in a period. Whenever this occurs, we try to recover feasibility by applying the four options of step 2 from the last period with overtime backwards.

Computational Results

All computational experiments in this section were performed on an IBM personal computer with a 3.2 GHz CPU and 1GB of random access memory. CPLEX version 9.0 from ILOG was used as the mixed integer programming solver. We used the instance generator described in Section 2.1.1 to generate random instances. For each type of instance, characterized by the quadruple N , T , Cut and θ, we present the average of the values obtained in 10 randomly generated instances.

We first compare the LP relaxations of models F 1 and F 2 , both strengthened with the (l, S) and (W t ) valid inequalities. Let F 1LP and F 2LP denote the values of the linear relaxations of F 1 and F 2 , respectively, and let F 1LP and F 2LP denote the values of the linear relaxations of F 1 and F 2 including all the (l, S) and W t cuts, respectively. The separation algorithms for inequalities (13) and ( 18) were coded in OPL version 3.7 from ILOG. We solved the minimum cut problem of the disconnected subtour separation algorithm with the Ford-Fulkerson's algorithm. The running times of the linear relaxation of F 2 strengthened with the (l, S) and (W t ) valid inequalities were always less than 150 seconds for all instances. As an example, with T = 10 and N = 25, model F 1 has 6760 columns and 4135 rows. 2LP be the value of the linear relaxation of F 2LP including the (W t ) cuts. The computational results are displayed in Table 4. Both (l, S) and (W t ) inequalities improve the lower bound, with an effect that increases (decreases) with the number of periods (products). Clearly, the improvement provided by the (l, S) cuts is especially noteworthy.

Table 5 reports the average deviation of the optimal solution from the lower bound F 2LP , i.e., its integrality gap. As in Section 2.1.1, the empty buckets mean that at least one of the instances could not be solved optimally within a 3600 seconds time limit. It is clear that the quality of the lower bound is significantly improved as the number of products (N ) increases, and the quality seems to be independent of the number of periods (T ). The results show that, on average, the integrality gap of the instances with θ = 50 is lower than the integrality gap of the instances with θ = 100. Nevertheless, this difference becomes smaller as the number of products in a problem instance increases. In both cases, from 25 onwards the gap is below 1%.

To generate an upper bound, we implemented the heuristic from Section 3 in C++ by using Visual Studio 6.0 as development environment. Tables 6 and7 2.7 3.1 2.5 5.6 6.7 5.9 7 1.5 2.1 2.3 4.3 4.5 4.6 10 1.5 0.9 0.9 2.5 3.1 15 1.0 0.9 1.5 25 0.8 1.0 0.7 maximum gap of the heuristic solution from the lower bound, for θ = 50 and θ = 100, respectively. The running times of the heuristic were always less than 0.4 seconds for all of the instances. 2.3/6.0/9.9 2.7/7.0/12.7 2.6/6.6/10.3 10 5.6/9.5/18.0 4.7/8.9/13.4 3.3/7.7/14.9 15 4.7/9.7/14.6 4.9/10.0/18.9 6.3/10.1/12.4 25 6.5/9.9/11.9 6.0/11.1/14.9 8.6/12.0/16.3 minimum / average / maximum gap (%) The heuristic found a feasible solution for all problem instances. The flexibility of the four options of step 2 of the heuristic enables us to find feasible solutions, even for tightly capacitated instances. The results indicate that for both θ = 50 and θ = 100 the heuristic performance deteriorates as the number of products increases. This situation is pronounced when the setup costs are considerable higher than the inventory costs (θ = 100). The performance of the heuristic is not influenced by the number of periods. The amplitude of the gap between the lower and upper bounds tends to decrease as the number of product increases. Recall that in our production setting from the glass industry, typically θ = 50 and thus Table 6 is more relevant. Based on Table 5, we also conclude that the most significant portion of the gap comes from the heuristic (and not from 17 To the best of our knowledge, [START_REF] Gupta | The capacitated lot-sizing and scheduling problem with sequencedependent setup costs and setup times[END_REF] is the only paper in the literature to publish results (gaps between lower and upper bounds) for this problem. However, since inventories and production quantities are expressed in units of capacity (normalized to one), they adapted the parameters of the instance generator used by [START_REF] Haase | Lot sizing and scheduling with sequence-dependent setup costs and times and efficient rescheduling opportunities[END_REF]. By comparing our parameters with the parameters of their generator, their instances are similar to the case with θ = 50. The authors present the performance of their heuristic from N = 2 to N = 15 for T = 4, and from T = 2 to T = 15 for N = 4. For problem instances with 15 products and 4 periods, [START_REF] Gupta | The capacitated lot-sizing and scheduling problem with sequencedependent setup costs and setup times[END_REF] report an average gap between lower and upper bounds of approximately 20%, and for problem instances with 4 products and 15 periods an average gap of approximately 30%. From Table 6 it clearly follows that our models and the heuristic are better.

Concluding Remarks

We have presented two novel exact formulations for modeling setup carryover in the challenging CLSP problem with sequence dependent setup times and costs. The models are simpler than other available in the literature. Note that all the models published so far only perform a setup whenever there is sufficient time available to do it entirely. In the case of very large setup times (as the ones observed in the glass container industry that may reach 120 hours), setup operations may start at the end of one period and finish at the beginning of the following period. This new feature has not yet been incorporated in CLSP models.

We have formally shown that one of the formulations is stronger than the other and we have implemented valid inequalities that enable us to reduce the integrality gap of the LP relaxation with a desirable property: the integrality gap reduces as the number of products increases (it is independent on the number of periods). It is clearly important to find strategies to combine even stronger valid inequalities based on the polyhedral structure of this problem with tighter reformulations. These inequalities should take into account the capacity constraints and the potential cycles and paths of products sequence in each time period.

We have described a five-step heuristic that is effective both in finding a feasible solution and in producing good solutions to these problems (measured by the gap between the lower and upper bounds). The simplicity of this heuristic enables its implementation as a building block of more complex heuristics and metaheuristics. It is worth pointing out that it may effectively generate good initial solutions and repair infeasibilities due to capacity requirements of solutions generated by any local search procedure.

Figure 1 :

 1 Figure 1: Possible configurations

  and upper bounds on lot sizes.

Figure 2 :

 2 Figure 2: Gantt chart of the optimal solution

Figure 3 :

 3 Figure 3: Network representation of the optimal solution

Figure 4 :

 4 Figure 4: Comparison of run-time performance between F 1 and GM models for Cut = 0.6

  10

  taken into account in the amending procedure (an option is only considered if the preceding one has not eliminated O t ).(a) Let S = {i ∈ [N ] :

FigureFigure 5 :

 5 Figure 5: Comparison of F 1LP and F 2LP for Cut = 0.6 and θ = 100

  display the minimum, average,

  Figure 6 gives an example of the above procedure. Instead of setting the time of already considered arcs to Θ, we delete such arcs.

Table 2 :

 2 Benchmark of models' complexity

	Model	Binary	Continuous Constraints
		variables	variables

Table 3 :

 3 Average running times for Cut = 0.6 and θ = 50 (model F 1 )

					T
	N	3	4	5	7	10	15
	3	0.1 0.1	0.1	0.2	0.2	14.7
	4	0.1 0.1	0.3	2.8	1.3	34.7
	5	0.2 0.3	0.7	7.3	7.9	300.1
	6	0.2 0.3	0.7	2.3	24.5	593.6
	7	0.5 1.2	2.0	5.1	310.7 1628.0
	8	0.7 1.1	4.7	37.7	162.9
	9	0.7 2.1	9.4	55.5	967.8
	10	1.4 1.2 39.7 313.1 2875.7

Table 4 :

 4 Comparison (%) of F

				(l,S) 2LP and F	(Wt) 2LP with F 2LP for Cut = 0.6 and θ = 100
		(F	(l,S) 2LP -F 2LP )/F 2LP	(F	(Wt) 2LP -F 2LP )/F 2LP
				T				T
	N	5	7	10		5	7	10
	5	102.1 177.1 276.0	16.1 44.1	86.1
	7	96.1 148.8 256.3	3.8 16.4	41.2
	10	87.1 142.7 231.9	0.9	3.3	14.3
	15	81.5 133.0 204.0	0.4	0.6	2.4
	25	79.5 124.5 191.3	0.2	0.2	0.3

Table 5 :

 5 Integrality gap (%) for Cut = 0.6

			θ = 50			θ = 100
			T			T
	N	5	7	10	5	7	10
	5

Table 6 :

 6 Gap (%) between the lower and upper bounds for Cut = 0.6 and θ = 50

			T
	N	5	7	10
	5	2.8/6.4/12.7	2.5/8.3/19.2	2.0/6.4/15.6
	7

Table 7 :

 7 Gap (%) between the lower and upper bounds for Cut = 0.6 and θ = 100

	T

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.ukInternational Journal of Production Research

Acknowledgments

The authors are grateful to Professor Stephen C. Graves from the Sloan School of Management for his helpful comments.

The first author is also grateful to the Portuguese Foundation for Science and Technology for awarding him the PhD grant SFRH/BD/23987/2005.

A The Sequencing Heuristic

In this section, we present a heuristic to sequence N products in order to minimize the sequence dependent setup times s ij . Since cycles are not taken into account and each product is produced, 18


J M Framinan

R Leisten
email: rainer.leisten@uni-due.de.

Keywords: SCHEDULING, FLOW SHOP SCHEDULING, META-HEURISTICS Scheduling, flow shop, due dates, greedy algorithm

Total tardiness minimisation in permutation flow shops: a simple approach based on a variable greedy algorithm

Introduction

A flowshop is a common manufacturing layout where the operation sequence for all jobs is the same. Additionally, it is usually assumed that jobs cannot overtake one another from machine to machine and therefore, once a sequence is fixed for all jobs on the first machine, this sequence is maintained for all machines (permutation flowshop). 3 Several objectives may be considered when scheduling in a flowshop. Without any doubt, the most widely studied is the makespan objective (see e.g. Reza [START_REF] Hejazi | Flowshop-scheduling problems with makespan criterion: a review[END_REF][START_REF] Ruiz | A comprehensive review and evaluation of flowshop heuristics[END_REF] for recent reviews on this problem). In contrast, significantly less attention has been paid until recent years to other objectives, particularly those involving due dates. Among due-date related objectives, minimisation of total tardiness aims to find schedules that satisfy the due dates committed to customers and it is therefore of significance and importance to real life considerations [START_REF] Hasija | Scheduling in flowshops to minimize total tardiness of jobs[END_REF].

In this paper, we propose an algorithm with the objective of minimizing the total tardiness of jobs. The algorithm is based on a variable greedy algorithm concept. A salient feature of the algorithm is that it employs a tailored design for the problem, which is based on the slack-concept, i.e. jobs with minimum slacks are removed from a given solution in order to find for them a better position in the schedule. Besides, the algorithm is parameter-free (with the exception of the running time), so it does not require any parameter setting and parameter tuning. This simplicity and the lack of control parameters is a desired characteristic of any heuristic approach as makes it more user-friendly [START_REF] Hansen | Variable neighbourhood search: Principles and applications[END_REF].

The remainder of the paper is organised as follows: in the next section we review in more detail the problem under consideration and the contributions up-to now. Next, we present the proposed algorithm for the problem. We then discuss and set up an extensive computational experiment in order to validate the suggested approach and to compare with the best-so-far ones. Finally, we present the main conclusions and some future research lines. The problem under consideration refers to the scheduling of n jobs in a permutation flow shop consisting of m machines. Each job i has a processing time on each machine j, which is denoted by p ij . Additionally, a due date d i for each job i is given. In a permutation flowshop, the completion time of the i-th job in a given schedule σ on a machine j, C [i]j (σ), is given by the following recursive equations:

C [i]j (σ) = max{ C [i-1]j (σ) ; C [i]j-1 (σ) } + p ij C [0]j (σ) :=0; C [i]0 (σ):=0;
With these equations, C i (σ), the completion time of a job i (i.e. the earliest time when this job can be can released to the customer), can be computed. T i , the tardiness of a job i is then defined as the maximum of the difference between its completion time and its due date, and zero, i.e.: T i = max{ C i (σ)d i ;0}. The total tardiness TARD is the sum of the tardiness across all jobs: TARD = Σ T i .

Clearly, tardiness minimisation is related to maintaining high service levels, as it aims to respect the committed due dates. Research on the topic has been carried out by [START_REF] Gelders | Four simple heuristics for scheduling a flowshop[END_REF], [START_REF] Kim | Heuristics for flow shop scheduling problems minimizing mean tardiness[END_REF], Parthasarathy and Rajendran (1998), [START_REF] Armentano | Tabu search for total tardiness minimization in flowshop scheduling problems[END_REF], and [START_REF] Hasija | Scheduling in flowshops to minimize total tardiness of jobs[END_REF].

Gelders and Sambandam develop constructive heuristics for the problem, while

Kim suggests an algorithm based on the NEH heuristic [START_REF] Nawaz | A heuristic algorithm for the m-machine, njob flow-shop sequencing problem[END_REF] for makespan minimisation. His proposal is not compared to these from Gelders and Sambandam. Parthasarathy and Rajendran present two versions of Simulated Annealing (SA) for the problem and compare them with those of Gelders and Sambandam, and that of Kim. They found their proposal to outperform all these heuristics. Armentano Ronconi using an adaptation of the benchmark problems by [START_REF] Taillard | Benchmark for basic scheduling problems[END_REF]. The results of this computational experience show that (at least for the test-bed employed), the heuristic by Hasija and Rajendran clearly outperforms the existing ones and henceforth it should be considered the best heuristic for the problem up to date. In their computational experience, the heuristic by Parthasarathy and Rajendran also outperforms the one by Armentano and Ronconi.

As a summary, it can be concluded that the most efficient heuristics for the problem under consideration are PR and HR. Both heuristics use the same heuristic for generating the initial sequence, i.e. the so-called Earliest Apportioned Due Date Method (EADD). In EADD, jobs are ordered for each machine in ascending order of a value proportional to the processing times of the job and its due date. In this manner, one order of jobs is obtained for each machine. Among these, the ordering with the least value of total tardiness is chosen. PR employs a Random Insertion Perturbation Scheme (RIPS) for generating new solutions from a given sequence in the following manner: for each job in the seed sequence, two random positions are generated (one on its right and other on its left) to insert it, thus obtaining two new sequences. Obviously, the first (last) job cannot be inserted on any position on its left (right), hence only one new sequence can be obtained from these extreme positions. In total, 2•(n -1) new sequences are obtained. Among them, the one with the least total tardiness is selected. Step Swap (PSS) perturbation scheme. Using JIBIS, each job (starting from job 1 to job n) in a seed sequence is selected for insertion in all possible positions of the sequence. If the total tardiness of the so-obtained sequence is lower than the one for the seed, it replaces the seed. Note that JIBIS obtains a solution which is not worse than the initial seed. In the JSB scheme, a job in the seed sequence is chosen based on a probabilistic function and then inserted either on the right or on the left of its original position (the same probability is assigned to the insertion on the left and on the right).

The job in the seed sequence is chosen according to the relative displacement of the job's position in the current sequence from the job's position in the best-so-far sequence. Finally, PSS generates sequences from a seed sequence by probabilistically swapping jobs close to each other, and far from each other. Among them, the one with the least total tardiness is chosen. In Figures 1 and2 we present the detailed pseudocode of both heuristics, keeping the same steps and the programme structure as in the original descriptions.

[INSERT FIGURE 1 ABOUT HERE] [INSERT FIGURE 2 ABOUT HERE]

Proposed heuristic

In recent years, a number of general-purpose heuristic approaches (also known as meta-heuristics) have been introduced to tackle combinatorial optimisation problems. π is better than π ). In addition, a local search procedure can be implemented before this step in order to improve ' π .

IG is closely related to another technique, Iterated Local Search (ILS), in which the algorithm iterates over a local search neighbourhood (see e.g. [START_REF] Lourenço | Iterated Local Search[END_REF].

To the best of our knowledge, IG has been only applied to few problems, including crew scheduling [START_REF] Marchiori | An evolutionary algorithm for large set covering problems with applications to airline crew scheduling[END_REF] and flowshop scheduling with makespan objective [START_REF] Ruiz | A simple and effective iterated greedy algorithm for the permutation flowshop scheduling problem[END_REF]. Salient features of IG are its simplicity and the possibility to use a tailored heuristic to obtain good solutions for the problem under consideration.

Clearly, an obvious drawback of the method is the stagnation due to the 'greedy' nature of the procedure. In order to overcome this problem, [START_REF] Ruiz | A simple and effective iterated greedy algorithm for the permutation flowshop scheduling problem[END_REF] adopt a simulated annealing-based acceptance criterion, i.e.: accepting worse solutions with some probability depending on the number of iterations of the algorithm. A similar approach is also adopted in [START_REF] Marchiori | An evolutionary algorithm for large set covering problems with applications to airline crew scheduling[END_REF]. In this paper, we propose a different procedure for overcoming local optima, which is based on the VNS concept [START_REF] Mladenović | Variable Neighborhood Search[END_REF]. VNS uses a finite set of pre-selected neighbourhood structures that we denote by N k (k= 1, …, k max ). Let us denote by N k (x) the set of solutions in the neighbourhood of x. The pseudo-code of the basic structure of a VNS algorithm can be described as follows:

(1) Generate an initial solution x;

(2) While termination conditions are not met, repeat:

(2.1) Set k ← 1; (2.2) Until k = k max , repeat: (2.2.1) Find the best neighbour x' of x ( x' ∈ N k (x) ); (2.2.2) If x' is better than x, then set x ← x' and set k ← 1; otherwise set k ← k + 1;
Sophisticated features can be added to this most basic procedure. For these, we refer the reader to [START_REF] Hansen | Variable neighbourhood search: Principles and applications[END_REF]. VNS has shown to be efficient for a number of NP-hard problems, including some scheduling problems (see e.g. [START_REF] Gagné | A new hybrid Tabu-VNS metaheuristic for solving multiple objective scheduling problems[END_REF][START_REF] Schuster | Approximative procedures for no-wait job shop scheduling[END_REF][START_REF] Tasgetiren | Particle swarm optimization algorithm for single machine total weighted tardiness problem[END_REF]). In our proposal, we use the VNS concept of varying the neighbourhood, but we apply it to the destruction and construction phases of the IG algorithm. To do so, a greedy mechanism specifically tailored for the problem is described in the next section.

Slack-based greedy algorithm

For a schedule σ := ( σ 1 , ...), let us define ( ) i s σ the slack of job σ i in the sequence σ in the following manner: , ( ) : Based on this idea, it is possible to define a greedy algorithm to obtain a sequence σ' from another sequence σ in the following manner:

i i i m s d C σ σ σ = - . It is clear that ( ) i s σ is sequence- F o r P e e r
(1) Remove from σ the k jobs with lowest ( )

i s σ . Store them in k R .
Let µ be the remaining subsequence from σ after removing the k jobs.

(2) For each job r k R σ ∈ , repeat:

(2.1) Insert r σ in all possible slots of µ (starting from the first position). Let ν be the best subsequence obtained, and b the position of r σ in ν .

(2.2) Perform an adjacent pairwise exchange among all jobs in positions b+1 to (nk + 1). Let π be the best sequence obtained by this approach.

(2.3) If π improves ν , then set µ π ← ; otherwise set µ ν ← .

(3) Set ' σ µ ← .

According to this approach, it is possible to define k different types of movements, depending on the number of jobs to be removed. Aside, this construction mechanism is related to the heuristic by Kim, which in turn is based on the NEH heuristic.

Nevertheless, we added step (2.2) as we try to allow some small perturbation of the solution for those jobs after the one for which the minimum tardiness has been found. In addition, it is possible to introduce a speed up in order to reduce the computational effort. When executing step (2.1), instead of calculating the completion times from scratch for every slot j (j = 1, …, nk + 1) in which the job r σ is inserted, the completion time in all machines for jobs preceding job (j -1) do not need to be

Insertion-based local search

In order to improve the solution provided by the slack-based greedy algorithm, a local search in the neighbourhood of the solution is accomplished. Usual neighbourhood definitions for flowshop scheduling problems include adjacent pairwise interchange (exchanging the positions of two adjacent jobs), general pairwise interchange (exchanging the positions of each two jobs), and insertion (removing one job from its position and inserting it into a different one). Among these possibilities, the experiences conducted produced the best result for the last one. Consequently, an insertion-based neighbourhood is adopted.

The local search algorithm then obtains n positions at random (without repetition), and for each position, the job occupying this position is removed and inserted in all possible slots. If the best solution out of the n insertions is better than the current solution, then it replaces the latter as the best incumbent solution. The pseudocode of this procedure is as follows:

(1) Let σ be the current solution.

(2) For i = 1 to n do:

(2.1) Obtain k = random (1, n) (without repetition).

(2.2) Remove k σ , the job in position k, in σ and insert it in all possible slots.

Let ' σ be the best solution found by this way.

(2.3) If ' σ improves σ, then set ' σ σ ← .

(3) Return σ as solution.

Pseudo-code of the algorithm

The algorithm operates as follows:

(1) While computation time does not exceed a given time limit, repeat:

(1.1) Generate a solution x at random.

(1.2) Set k ← 1;

(1.3) Until k = (n -1), repeat:

(1.3.1) Obtain x' from x according to the procedure explained in Section 3.1. i.e.

using the slack-based greedy algorithm.

(1.3.2) Obtain x'' from x' according to the procedure explained in Section 3.2. , i.e. using the insertion-based local search.

(

1.3.3) If x'' is better than x, then set x ← x'' and set k ← 1; otherwise set k ← k + 1;
The algorithm starts with a random solution. Then, this solution is improved by a stepwise approach that systematically enhances the capabilities of the greedy algorithm until it shuffles all jobs. If a better solution is found, then the process starts again.

Otherwise, if the greedy algorithm has been expanded without reaching a better solution (i.e., k = n -1), then a new random solution is tried. This allows a diversification of the exploration once all greedy structures have been tried. Note that this algorithm is parameter-free, with the exception of the CPU time.

With respect to the stopping criterion, it has been chosen as the decision interval for scheduling problems is usually very short and hence the total running time of a heuristic is a critical issue. However, replacing the original criterion by a different one 12 based e.g. on the number of solutions evaluated, the number of iterations without improvement, etc. is straightforward.

Computational experiments

In this section we present the results of the computational experiments carried out in order to validate our proposed procedure. We first discuss the experimental design and particularly the setting of due dates in the problem instances. We then compare our proposal with the best-so-far algorithms for the problem under consideration. Finally, we analyse how the different design decisions in the algorithm (namely the design of the destruction phase, and the concept of 'variable destruction') affect its performance.

Experimental design

It is clear that the design of the experimental test-bed may have a huge influence on the results. Particularly, being a due-date problem, several types of problems can be distinguished, each one with a different degree of difficulty. Henceforth, we first review the different settings for due dates that have been adopted in the literature.

Rather different procedures for establishing due dates in permutation flowshops can be found in the papers by [START_REF] Gelders | Four simple heuristics for scheduling a flowshop[END_REF], [START_REF] Armentano | Tabu search for total tardiness minimization in flowshop scheduling problems[END_REF], [START_REF] Demirkol | Benchmarks for shop scheduling problems[END_REF], and Hasija and Rajendan (2004). In all these, the processing times of the problems are obtained from a uniform distribution. In fact, the where T and R are the tardiness factor of jobs and the dispersion range of due dates, respectively, and P is a lower bound of the makespan, defined as follows:

1 This due-date setting would result in approximately 50% of the jobs being late, and it has been used also by [START_REF] Chakravarthy | A heuristic for scheduling in a flowshop with the bicriteria of makespan and maximum tardiness minimization[END_REF], and [START_REF] Framinan | A heuristic for scheduling a permutation flow shop with makespan objective subject to maximum tardiness[END_REF] due to its ability to produce tight due dates.

1 1 1 1 1 max max min min , max j n m m ij il il ij i i j m i i l l j j P p p p p - ≤ ≤ = = = + =       = + +           ∑ ∑ ∑ ∑ With the

Comparison with existing heuristics

As mentioned in the introduction, the procedure proposed by Hasija and Rajendran (HR in the following) is the best algorithm so far for the problem under consideration, Moreover, it is not possible to establish a fair comparison based on the number of solutions evaluated by each algorithm, as our proposal constructs (complete) solutions by exhaustively selecting and exploring partial sequences, therefore the number of final solutions evaluated is much lower than in PR and HR (less than 1% in the DMU testbed, as we will show later). Consequently, the fairest comparison with HR is running this algorithm exactly with the same parameter setting described by the authors, and then compare its results to those obtained by the new heuristic proposal, using as termination criterion the same CPU time taken by the HR algorithm for this problem.

Each instance is solved by 24 runs of each algorithm. The average and the best value of the total tardiness for each algorithm and instance are obtained. Then, the quality of the solutions obtained for each instance and each heuristic i is measured in terms of RPD (Relative Percentage Deviation) with respect to the average and the best solution, in the following manner:

* * ( ) 100 i i T T PRD T - = × * * * * ( ) 100 i i T T PRD T - = ×
where i T is the average tardiness obtained over the 24 runs of algorithm i, * i T is the best tardiness among the 24 runs, and * T is the best known solution for this instance.

For all instances of a given problem size, the average and the standard deviation of both

[INSERT TABLE 1 ABOUT HERE] [INSERT TABLE 2 ABOUT HERE] [INSERT TABLE 3 ABOUT HERE]

As can be seen from Tables 1,2, and 3, the algorithm proposed exhibits a very good performance for the three test-beds. Both the best solution (out of 24 runs) and the average obtain better results than those of HR and PR.

In order to test the performance of our proposal for different CPU times, we change the original design of PR and HR in order to introduce CPU time as stopping criterion. Then, we run all three algorithms with fractions of the CPU time invested for the original design. More specifically, let TIME be the CPU time required (average over 24 runs) for a given instance in the experiments presented in Tables 1, 2, and3 4 in terms of the best solution and the average over 24 runs of each algorithm.

[INSERT TABLE 4 ABOUT HERE]

As can be seen from Table 4, the algorithm outperforms HR and PR for the different CPU times. It is interesting to note that, for all three algorithms, the results for different CPU times are similar with respect to the best solution. Regarding the average solution, there is a steady improvement for our proposal, but not for HR and PR, which may indicate that both algorithms could stagnate for longer CPU times.

Finally, in Table 5 we give the average number of solutions explored by HR and our proposal in the DMU test bed for each problem size (the CPU allowed is that of the original HR description, see Table 1). It has to be noted that this value is not an indicator of the computational effort of each algorithm, as in the case of our proposal, most of the effort is devoted to insert the jobs removed in the destruction phase in the different slots of the remaining schedules. Therefore, the number of (complete) solutions explored by our proposal is very low as compared to HR, which use most of the CPU time to explore and evaluate new (complete) solutions.

[INSERT TABLE 5 ABOUT HERE]

Effectiveness of the variable greedy design

The proposed algorithm is based on a variable neighbourhood depending on k, which determines the number of elements in the sequence to be destroyed and later reconstructed in a range from 1 to (n -1). Therefore, it is of interest to assess the Are there alternative upper values of k (other than n -1) for which the algorithm yields a higher performance?

In order to answer the first question, we analyse the contribution of each value of k to improve the current solution of the algorithm. More specifically, for a given problem instance, we record the value of k for which a new best-so-far solution is obtained. We then calculate the frequency of each value of k and plot the results. In Figure 3, we show the plot for different problems in the T-HR test-bed. Each line shows the frequency of the improvements of k for 20 runs of 10 instances of a given (n,m).

[INSERT FIGURE 3 ABOUT HERE]

In Figure 3, it is shown that all values of k in the range contribute to the improvement of the solution. The frequency decreases with k, which is not surprising, since in the algorithm, every time a new solution is found, the value of k is set to 1 and starts increasing once the solutions found with this destruction level are worse than the best-so-far sequence. The only exception is k = 1, indicating the relatively low ability of this destruction level to find better solutions as compared to k = 2. Nevertheless, the second highest frequencies are obtained for k = 1. Finally, it has to be noted that the pattern of the frequencies is the same for different number of machines. In a different experiment also aimed to the first question, we compare the results obtained by this algorithm in the DMU test-bed with those obtained using different (fixed) numbers of components in the destruction phase. For all runs, the stopping criterion is the same as in the previous sub-section, i.e. the CPU time required for each problem by the HR heuristic. Consequently, the CPU time is again the same for all algorithms. The results are shown in Table 6 in terms of the APRD with respect to the best solution found by the different settings.

[INSERT TABLE 6 ABOUT HERE]

In view of the results in Table 6, the best overall results are obtained by the variable destruction design. While some results with a fixed k are very good, these are strongly dependent on the particular problem size.

Once the convenience of using a variable k as compared to employing fixed values has been established, we check whether there exist other upper values of k yielding better results than the original proposal (n -1). To do so, we run different versions of our heuristic with several upper levels of k for the DMU test-bed. More specifically, we test n/10, n/5 and n/2. All versions are allowed to use the same CPU time (that employed for the comparison of the different heuristics in Section 4.2., see Table 1).

The results are shown in Table 7 in terms of the RPD.

[INSERT TABLE 7 ABOUT HERE]

As can be seen from Table 7, employing other values as upper limit of the level of destruction is less efficient than the original proposal. Particularly, incorporating large confirming the results in Figure 3 showing that also large k contribute to improve current best solutions.

Effectiveness of the slack-based greedy algorithm

In this subsection we investigate whether the mechanism adopted for selecting the k components to be destroyed (i.e. those k with least slack) is more efficient than, for instance, picking these k at random, or picking them in reverse order (descending order of their slacks). The results are shown in Table 8 in terms of the ARPD with respect to the best results obtained for the three approaches. In this table, 'NEW' denotes the approach proposed, while 'random' and 'inverse' indicates picking the components in the destruction phase at random, or with descending order of the slacks, respectively.

[INSERT TABLE 8 ABOUT HERE]

As can be seen from Table 8, the slack-based definition of the greedy algorithm is, on the overall results, more efficient than picking the solutions at random or in descending order of their slack. It is to note that, for the problem size n = 50 and m = 15, the inverse ordering obtains slightly better results. However, this small difference is not supported for the rest of the problem sizes.

Conclusions

We have presented a new algorithm for the problem of minimising total tardiness in permutation flowshops. A salient feature of the proposed method is that it does not require any parameter, with the exception of the computation time. Step 1: Obtain one initial sequence using heuristic EADD and assign it to S and B.

Step 2: Initialize SA parameters, i.e. T:=475; ACCEPT := 0; TOTAL:=0; FR_CNT:=0.

Step 3: If (FR_CNT = 5) or (T ≤ 20) then go to Step 14.

Step 4: Invoke perturbation routine RIPS to obtain S'. Step 5: S:= S'; ACCEPT := ACCEPT + 1

Step 6: If ( ') ( ) 100 ( )

TARDINESS S TARDINESS B TARDINESS B

-⋅ > 0 then go to Step 10.

Step 7: B := S'; FR_CNT := 0. Go to Step 10.

Step 8: If exp( ) T ∆ -< Random[0;1], then go to Step 10.

Step 9: S:= S', and ACCEPT := ACCEPT + 1.

Step 10: TOTAL:= TOTAL + 1.

Step 11: If (TOTAL > 2 n) or (ACCEPT > n/2) go to Step 12; else go to Step 4.

Step 12: If (ACCEPT• 100/TOTAL ≤ 15) then FR_CNT := FR_CNT + 1

Step 13: Set T := T • 0.9; ACCEPT := 0; TOTAL := 0. Go to Step 3.

Step 14: The algorithm is frozen. B contains the final solution. Step 1: Obtain one initial sequence using EADD. Assign it to S and B.

Step 2: Apply JIBIS to S and B.

Step 3: Z = TARDINESS(S), Z B := TARDINESS(B)

Step 4: Initialize SA parameters, i.e. T :=540, COUNT := 0; REPLACE :=0; FREEZE :=0.

Step 5: Generate and evaluate nine random sequences and store them (together with B) in an archive of the best ten sequences.

Step 6: If (FREEZE = 5) or (T ≤ 20) then go to Step 21.

Step 7: Set Z'' a large number.

Step 8: Do the following steps twice: Do the following steps n times:

Invoke JSB to obtain S'. Z':= TARDINESS(S').

If Z' is less than the maximum value of the objective function of a sequence amongst all sequences in the archive, then replace that sequence and its objective function by S' and Z'. If Z' < Z'' then S'':=S' and Z'':=Z'.

Step 9: Do the following steps twice:

Invoke perturbation scheme PSS to obtain S'; Z':=TARDINESS(S')

If Z' is less than the maximum value of the objective function of a sequence amongst all sequences in the archive, then replace that sequence and its objective function by S' and Z'.

If Z' < Z'' then S'':=S' and Z'':=Z'.

Step 10: ( ') ( ) 100 ( )

TARDINESS S TARDINESS S TARDINESS S -∆ = ⋅

Step 11: If ∆ > 0 then go to Step 15.

Step 12: S:= S''; Z:=Z''; REPLACE = REPLACE + 1.

Step 13: If ( '') ( ) 100 ( )

TARDINESS S TARDINESS B TARDINESS B

-⋅ > 0, then go to Step 17.

Step 14: B:=S''; Z B := Z''; FREEZE = 0. Go to Step 17.

Step 15: If exp( ) T ∆ -< Random[0;1], then go to Step 17.

Step 16: S:= S''; Z:=Z''; REPLACE := REPLACE + 1.

Step 17: COUNT := COUNT + 1.

Step 18: If (COUNT ≤ 2 n) go to Step 7.

Step 19: If REPLACE • 100/COUNT ≤ 15, then FREEZE := FREEZE + 1.

Step 20: T = T • 0.9; COUNT := 0; REPLACE := 0. Go to Step 6.

Step 21: Implement JIBIS on all sequences stored in the archive. The best amongst them is returned as solution.

  schemes. They compare their algorithm with the most efficient version of the heuristic

  mechanisms for the generation of new solutions, namely Job Index Based Insertion Scheme (JIBIS), Job Shift-Based (JSB) scheme, and Probabilistic

  this procedure, the computation times are reduced roughly to the half of the time.

  above setting, they define four scenarios depending on the values of T and R. This way of generating due dates is nearly identical to the one employed by Demirkol et al. who build a test-bed for several due date problems, including 160 problem instances of flowshops with due dates. They employ different T and R factors that produce much tighter due dates as compared with the paper by Armentano and Ronconi. Additionally, from a bigger subset of problems, the 160 hardest instances were collected. Gelders and Sambandam propose a different method for generating due dates, consisting of employing the following uniform distribution: 5 × mean processing time of a job on machine m].

Finally, we describe

  the due date setting in the test-bed developed byHasija and Rajendran, which is based on Taillard's test-bed[START_REF] Taillard | Benchmark for basic scheduling problems[END_REF].Taillard's test-bed consists of a collection of 120 flowshop problem instances chosen so as 'good' makespan values are difficult to obtain. In this test-bed, (n,m) ∈ {(20,5), (20,10), 20)}. For each problem size, 10 instances are available. In order to provide due dates d i for the problem under consideration, Hasija and Rajendran propose the following method to generate them: where u is a real random number uniformly distributed in the interval [0,1]. In view of the different contributions, we select three different test-beds for the comparison of the heuristics. For these, we use the following acronyms: DMU: The test-bed by Demirkol et al. As mentioned before, this test-bed consists of 160 instances, where (n,m) ∈ {(20,15), (20,20), (30,15), (30,20), (40,15), (40,20), (50,15), (50,20)}. For each problem size, there are 20 instances. T-GS: The test-bed obtained using the processing times of the test-bed of Taillard, and due dates generated according to the method proposed by Gelders and Sambandam. In view of the prohibitive computation times required to solve the 10 biggest instances of Taillard's test-bed (instances ta110-ta120, n = 500) by using the existing algorithms for the problem, we exclude them from T-GS. Therefore, the test-bed under consideration consists of 110 instances. T-HR: The test-bed obtained using the processing times of the test-bed of Taillard, and due dates generated according to the method proposed by Hasija and Rajendran. Analogously to the previous test-bed, the 10 biggest instances by Taillard are removed.

  procedure of Parthasarathy and Rajendran (PR in the following). In order to compare them with our algorithm, we employ the three test-beds described in the previous section. Note that neither HR nor PR have a simple termination criterion (such as the computation time, or the number of iterations without improvement), but use a rather sophisticated combination of parameters depending on the temperature (number of iterations, number of accepted moves, and number of sequences generated).

  are obtained, and shown in Tables 1, 2, and 3 for the DMU, T-GS, and T-HR test-beds, respectively. The tables show the quality of the solutions obtained in terms of the i ARPD (Average Relative Percentage Deviation) of heuristic i for each problem size. tardiness obtained by the application of the heuristic i to the problem instance, and min RPD is the minimum tardiness found by any heuristic for this instance.

  . Then we solve the DMU test-bed with all three algorithms for a CPU time of TIME/2, and The quality of the solutions obtained is shown in Table

  design to the performance of the algorithm. More specifically, we are interested in answering the following questions: a) Does the variable design contribute to improve the solution with respect e.g. fixed value of k? b)

  suitable to obtain better results for a given CPU time,

Figure

  Figure 1. Pseudocode of PR

FigureFigure 3 .

 3 Figure 2. Pseudocode of HR

  it is possible to express the total tardiness in terms of the slacks, i.e.:

	dependent and that ( ) max[ ( );0] i TARD s σ σ = -
	i
	R
	e
	v i e
	w
	O n l
	y
	9

∑

. Consequently, if a job i σ has a negative value for ( ) i s σ , then this job is contributing to increase the total tardiness.

  [1, 99] distribution is employed for all test-beds except for the one byDemirkol et al., where processing times are drawn from a[1, 200] uniform distribution. Therefore, in the following we focus our discussion on how due dates are established in these papers.

	[ (1 Uniform P	T --	2 R	), (1 P	T -+	2 R	)]
	F o r
		P
		e
			e r
					R
						e
							v i e
							w
							O n l
							y
	Armentano and Ronconi suggest using the following distribution to generate due
	dates:
							13

Table 1 .

 1 need to adjust or tune it for different problems. The computational experience carried out shows the proposal to outperform the best-so-far available algorithms. The algorithm is shown to be efficient for different computational efforts, and its design is checked to be more suitable as compared to different options regarding a variable/fixed destruction level, different upper limits for k, and other mechanisms for selecting the components to be destroyed. Further sophisticated features may be added to the suggested algorithm, such employing a more complex criterion for accepting new solutions. However, these enhancements must be contrasted with the simplicity and the speed of the current design. ARPD obtained in DMU test-bed (all heuristics use the same CPU time)

	This means that

Table 4 .

 4 ARPD obtained in for different CPU times (average on 24 runs)

	n	m	HR	NEW
	20	15	109152,5	1947,6
	20	20	108766,9	1796,2
	30	15	243373,2	2223,3
	30	20	243452,5	2014,5
	40	15	431356,5	2449,7
	40	20	432307,4	2235,5
	50	15	673924,0	2771,9
	50	20	673929,0	2474,2

Table 5 .

 5 Number of solutions explored by HR and NEW in the DMU-testbed

	n	m k=1 k=2 k=10 k=19 k=20 k=29 k=30 k=39 k=40 k=49 NEW
	20 15 2.14 1.79 0.39 0.42	0.40
	20 20 1.68 1.42 0.45 0.17	0.37
	30 15 3.68 3.88 1.88 1.10 1.00 0.82	0.84
	30 20 3.16 2.81 0.84 0.61 0.62 0.64	0.58
	40 15 4.56 4.79 1.94 1.64 1.69 1.92 1.14 1.68	1.51
	40 20 3.83 3.88 1.66 1.23 1.12 1.11 1.33 1.14	1.36
	50 15 6.07 6.78 3.19 1.87 2.21 1.96 1.12 2.00 1.06 1.66 2.28
	50 20 4.29 4.30 1.67 1.69 1.73 1.17 1.49 1.51 1.44 1.58 1.15
	Average
	:	3.68 3.71 1.50 1.09 1.40 1.27 1.27 1.58 1.25 1.62 1.06

Table 5 .

 5 ARPD for different greedy algorithms for the instances in the DMU test-bed.

				Upper limit of k
	n	m	n/10	n/5	n/2 n-1 (NEW)
	20	15	8.165	3.833	0.291	0.225
	20	20	6.956	3.473	0.500	0.203
	30	15	11.503	4.137	0.465	0.544
	30	20	7.938	3.488	0.594	0.235
	40	15	13.850	5.017	0.917	0.515
	40	20	9.831	3.478	0.519	0.511
	50 50 F r o	15 20	14.022 9.773 10.255	5.428 4.439 4.161	0.602 0.769 0.582	0.535 0.601 0.421
	P
	e
		e r
			R
			e
			v i e
				w
				O n l
					y

Table 6 .

 6 ARPD for several upper limits of k for the instances in the DMU test-bed (n-1 corresponds to the original proposal).

	n	m	NEW	inverse	random
	20	15	0.390	0.469	1.605
	20	20	0.147	0.230	0.481
	30	15	0.341	0.622	1.807
	30	20	0.258	0.597	1.803
	40	15	0.977	1.023	3.225
	40	20	0.407	0.603	1.490
	50	15	1.110	0.949	3.573
	50	20	0.530	0.704	2.497
	Average:	0.520	0.650	2.060
	Table 7. ARPD for different mechanism for picking the components to be removed
			(DMU test-bed).

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.ukInternational Journal of Production Research

Acknowledgements

The authors wish to thank the referees for the insightful comments on earlier versions of the paper. This research was accomplished while the first author was with the University of Duisburg-Essen (Germany) as recipient of an Alexander Von Humboldt Fellowship.