C Gicquel
email: celine.gicquel@ecp.fr

N Miègeville

M Minoux

Y Dallery

Optimizing glass coating lines: MIP model and valid inequalities

Keywords: Integer programming, Branch & Bound, Valid inequalities, Production line design, Glass coating

Glass coating is a specific transformation aiming at improving glass performance. The work presented in this paper deals with the determination of the optimal configuration of the production lines used to perform this operation. We propose a first MIP formulation of the problem and then discuss several types of valid inequalities for improving it. The main idea is to exploit explicit or implicit binary exclusion constraints to derive stronger valid inequalities: the maximal clique constraints. Efficient (polynomial time) separation algorithms exploiting special structure of the problem are described, giving rise to a cutting-plane generation procedure for strengthening the initial formulation. The computational study carried out shows that, with the enhanced formulation, good solutions can be obtained within reasonable computation times using currently available integer programming software.

Introduction

This work is motivated by an industrial problem arising in the glass industry in connection with a specific transformation of flat glass called glass coating.

Glass coating consists of depositing in vacuum thin layers of metal on the surface of glass sheets. As a general rule, the process involves several layers of distinct metals. This aims at giving the glass additional properties such as a better thermal insulation: see e.g. [START_REF] Arnaud | Industrial production of coated glass: future trends for expanding needs[END_REF] for an overview on the applications of coated glass. According to the sequence and thickness of the layers, the property obtained is different: hence production managers have to cope with some product diversity.

Glass coating can be done on specific production lines called "soft-coating lines" using a process called "cathodic sputtering" [START_REF] Suzuki | State of the art in large area vacuum coatings on glass[END_REF]). Basically, these lines are made of a number of metallic cathodes, each being used to spray or "sputter" a specific metal on the glass sheets. Each sheet can go only once through the production line: during this single passage, all the metal layers to be deposited on the sheet must be sputtered following the sequence imposed by the product specifications.

The cathodes are ordered along the line: a configuration of the line corresponds to a sequence of cathodes. A cathode contains a finite volume of a single metal. Once the metal of a cathode has been used up, the cathode must be changed. But, due to technical reasons, this requires a line shutdown during several days. Because of these time-consuming changeovers, soft-coating lines are operated according to the following organization. All cathodes on the line are changed together during a line shutdown. After this, production takes place continuously with this configuration during the next production run, the duration of which is typically about one month. When the run is over, all cathodes are changed and a new configuration is set up.

The problem addressed in the present paper concerns the determination of the optimal configuration to be set up between two line shutdowns. This decision can be based on reliable future demand forecasts: the requested products and the anticipated surface to be coated are assumed to be perfectly known. The configuration set up at the beginning of a production run should be able to process all needed products in the quantity requested until the next production shutdown. In this context, determining the configuration to be set up consists of selecting among a set of available cathodes the ones to be placed on the line, ordering them along the production line and deciding how to use them to process the requested products. Because of its limited capacity, a cathode may not be sufficient to sputter the entire volume needed to process a given layer. Thus we have to consider the situation where a layer is sputtered by several cathodes placed at different positions on the line. The objective is to minimize the number of cathodes to be placed on the line. Indeed, the larger the number of cathodes to be placed on the line during a setup, the greater the changeover operations will be and the more time will be lost for useful production. Investigation of models and algorithms for solving the resulting discrete optimization problem is the subject addressed in the present paper.

The problem under study shares some common features with a string processing problem called the Shortest Common Supersequence problem (see e.g. [START_REF] Maier | The complexity of some problems on subsequences and supersequences[END_REF]). It is however significantly different due to various extra constraints which must be taken into account, one of the most significant being the limitations imposed on cathode capacity, which frequently result in the use of a significant number of additional positions.

The problem of optimizing glass coating lines can also be related to the "Assembly Line Design Problem" (ALDP). In the ALDP, a production line is described as a series of workstations, each being responsible for performing a specific set of assembly tasks. The problem consists of selecting a piece of equipment for each workstation and deciding which tasks should be performed by which workstation. Recent overviews on the literature on the ALDP can be found in [START_REF] Becker | A survey on problems and methods in generalized assembly line balancing[END_REF] and [START_REF] Boysen | A classification of assembly line balancing problems[END_REF].

Nevertheless, the glass coating line problem is different from the ALDP studied in most papers (see e.g. [START_REF] Pinnoi | Assembly system design: a branch and cut approach[END_REF] and [START_REF] Bukchin | Design of flexible assembly line to minimize equipment cost[END_REF]). The main reason is that in the ALDP, each assembly task is performed exactly once, i.e is assigned to a single workstation on the line, whereas on a glass-coating line, a layer can be sputtered by several cathodes placed at different positions on the line. The problem of optimizing glass coating lines can thus be seen as an extension of the ALDP to the case where a task can be assigned to more than one workstation ("parallelized" in the terminology of [START_REF] Boysen | A classification of assembly line balancing problems[END_REF]). To the best of our knowledge, the only solution approach already available to deal with this particular extension of the "Assembly Line Design Problem" can be found in [START_REF] Boysen | A versatile algorithm for assembly line balancing[END_REF] who propose a flexible heuristic search procedure that can be modified to solve various extensions of the "Assembly Line Balancing Problem". In their paper, the authors assume that the processing time of a parallelized task is equally allocated to the chosen workstations. On the contrary, on a glass coating line, the volume of a layer sputtered by several cathodes can be unequally divided among the various cathodes so that we have to decide about the allocation of the metal volume to be sputtered among the chosen cathodes. Moreover, their solution approach is purely heuristic whereas ours being based on a mixed integer linear programming (MIP) model is intended to provide exact optimal solutions. The paper is organized as follows. In section 2, we introduce an initial mathematical formulation of this problem as a mixed integer linear program.

In section 3, we consider several ways to strengthen this initial formulation by adding valid inequalities of various types. In section 4, we discuss the results of some computational experiments showing the practical usefulness of the proposed valid inequalities at improving the efficiency of a Branch & Bound type procedure. Conclusions and perspectives for future work are presented in section 5.

Problem formulation

We wish to determine the optimal configuration of a glass coating line to be set up between two production shutdowns. In this section, we introduce an initial formulation for this optimization problem as a mixed integer linear program. To describe the problem precisely we introduce the following notation.

The set of anticipated requirements is supposed to involve M metals and P distinct final products. Each metal type is indexed by m: m = 1, 2, ..., M .

Each product, indexed p = 1, 2, ..., P , is made of a glass sheet on which O p layers are to be sputtered. For a given product p, a layer o = 1, 2, ..., O p is made of a specific metal denoted m po and its thickness is given by e po . The anticipated surface of product p to be processed during the production run, S p , being known, the volume of metal m po needed to sputter the o th layer of product p can be deduced as:

V po = e po * S p .
Possible positions of cathodes on the production line are indexed by i = 1, 2, ..., N . These positions are ordered according to the orientation of production flow. The cathode i + 1 is located immediately after the cathode i along the line.

Available cathodes correspond to C types of cathodes. For each type of cathode c = 1, 2..., C, we assume that we know m c , the corresponding metal, V c , the volume of available metal in each cathode and ν c , the number of cathodes belonging to this type. We agree to use an additional type c = 0 (the empty cathode) to represent free positions on the line. For each metal m, we denote C(m) the subset of cathode types c such that m c = m. The complementary subset is denoted C(m) = {c = 1, 2, ..., C st m c = m}.

First formulation of the problem as a MIP

Here we first provide a mathematical statement of the problem involving the following decision variables:

-z i c = 1 if a cathode of type c is placed in line position i, z i c = 0 otherwise. -y i po = 1 if the cathode placed in position i is used to sputter the o th layer of product p, y i po = 0 otherwise. -x i po gives the proportion of V po sputtered by the cathode placed in the i th position. Thus all the x i po are continuous variables in [0; 1].

Considering the criterion of minimizing the total number of positions used, the formulation proposed is:

min N i=1 C c=1 z i c (1) ∀i, C c=0 z i c = 1 (2) ∀c, N i=1 z i c ≤ ν c (3) ∀i, ∀p, ∀o, y i po + c∈C(mpo) z i c ≤ 1 (4) ∀p, ∀o, N i=1 x i po = 1 (5) ∀i, ∀p, ∀o, x i po ≤ y i po (6) ∀p, ∀(o, o ) st o > o , ∀(i, i ) st i < i , y i po + y i po ≤ 1 (7) ∀i, ∀m, c∈C(m) V c z i c - (p,o) st m po =m x i po V po ≥ 0 (8) ∀i ∈ [1; N -1], z i 0 ≤ z i+1 0 (9)
The objective expressed by (1) is to minimize the total number of cathodes placed on the line. Constraints (2) ensure that at most one cathode is placed in position i. z i 0 = 1 means that the i th position on the line is free. Constraints (3) ensure that no more than the number of available cathodes of type c, ν c , are placed on the line. Constraints (4) guarantee the compatibility between the metal m po for layer o of product p and the metal of the cathode placed in position i: we cannot open the connection y i po = 1 if the cathode in position i contains a metal other than m po . Equalities (5) ensure that the demand is perfectly met: all the volume of each layer should be sputtered. Constraints (6) link the continuous variables x i po with the binary variables y i po : some volume of the o th layer of product p can be sputtered in position i only if

the connection is open. Precedence constraints (7) force compliance with the order according to which the layers of a product should be sputtered: for a given product p, the layer o which is above the layer o should not be processed with a cathode i placed before the cathode i' if the latter is used to sputter o . Inequalities (8) guarantee that the limited capacity of the cathodes is not exceeded: for each type of metal, the volume remaining at the end of the production run in the cathode placed in position i should be nonnegative. Constraints (9) are used to enforce consecutive empty positions at the end of the line in case all positions are not used.

A small illustrative example

Problem P0 is a small instance we use in order to illustrate the problem and its resolution. P0 involves M = 4 metals, P = 3 products made of 3 or 4 layers and N = 12 positions on the line. Table 1 gives the numerical data relative to this example. The optimal configuration in this case is a sequence of Z * = 7 cathodes. Table 2 gives this sequence as well as the optimal use of cathodes to process the 3 products. We may notice that the first layer of product p = 3 is sputtered by two cathodes made of gold (placed at positions 2 and 4). This is due to the fact that the volume of metal needed to sputter this layer exceeds the capacity of a single cathode made of gold. In the sequel, P0 is used to illustrate various features of the proposed resolution method.

Valid inequalities

The formulation introduced in section 2 enables us to solve exactly only small instances: computation times for industrial problems of larger size using one of the best currently available commercial MIP solver are prohibitively long as can be seen from table 4. A possible explanation for this lies in the ob-

i Cathode p = 1 p = 2 p = 3 1 Ag, 1000
(1,1) 210 (2,1) 410 2 Au, 3000

(1,2) 400 (3,1) 1200 3 Ti, 2000

(1,3) 100 (2,2) 800 4 Au, 3000

(2,3) 200 (3,1) 2800 5 Pt, 2000

(3,2) 1000 6 Ti, 2000

(3,3) 1000 7 Ag, 1000

(1,4) 100 servation that the linear relaxation of the problem (1)-( 11) only provides a poor approximation to the exact optimal integer solution values. In order to address this issue, we investigate below several ways of strengthening the initial formulation (i.e. of reducing the integrality gap). The enhancements discussed here focus on various aspects of the problem under study, namely:

-available cathodes have a limited capacity, -only one metal can be assigned to each position on the line, -precedence constraints between layers of a given product must be respected.

In section 4, computational experiments will be reported showing that, thanks to these enhancements, the linear relaxation is tightened and instances of significantly larger size can be solved exactly with standard integer linear programming tools.

Valid inequalities from limited capacity of available cathodes

For each metal, we can compute a lower bound on the number of cathodes containing this metal to be placed on the line. This gives M valid inequalities ( 12) that can be added to the formulation.

∀m,

N i=1 c∈C(m) z i c ≥ P p=1 o=1...Op st mpo=m V po M ax{V c , c ∈ C(m)} (12)
Namely, for each metal m, dividing the global volume of metal needed to process all final products by the volume contained in the maximum capacity cathode containing metal m and rounding up gives the minimal number of cathodes of type c ∈ C(m) to be placed on the production line.

Valid inequalities from metal compatibility constraints

In this subsection, we discuss another family of valid inequalities to further strengthen the formulation. In a first step, we derive a series of binary exclusion constraints. These constraints are logical consequences of the formulation ( 1)-( 11). In a second step, we exploit the special structure of these constraints to derive stronger valid inequalities which correspond to maximal clique constraints in the underlying graph. (See e.g. [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF]). We note here that similar approaches have been used on other optimization problems such as assembly line design [START_REF] Pinnoi | Assembly system design: a branch and cut approach[END_REF]), harvest scheduling [START_REF] Goycoolea | Harvest scheduling subject to maximum area restrictions: exploring exact approaches[END_REF]), cellular telecommuni-cations networks design [START_REF] Kalvenes | Hierarchical cellular network design with channel allocation[END_REF] or air line crew scheduling [START_REF] Zeghal | Modeling and solving a crew assignment problem in air transportation[END_REF]). We observe however that in all the abovementioned references the structures of the underlying constraint graphs were significantly different from those studied in the present paper, leading to clearly distinct separation algorithms. In particular, in [START_REF] Pinnoi | Assembly system design: a branch and cut approach[END_REF], the separation of clique constraints is carried out using either complete enumeration or a greedy heuristic whereas our separation algorithms are exact and polynomial.

We first state various families of binary exclusion constraints. These constraints are implied by the constraints ( 2)-( 11) of the initial formulation but their explicit statement turns out to be useful with respect to strengthening.

They link pairs of binary variables related to the same position i on the production line, but to different products, layers or types of cathodes:

∀i, ∀c, ∀p, ∀o st m c = m po , z i c + y i po ≤ 1 (13) ∀i, ∀p, ∀o, ∀p , ∀o st m po = m p o , y i po + y i p o ≤ 1 (14) ∀i, ∀p, ∀(o, o ) st o = o , y i po + y i po ≤ 1 (15) ∀i, ∀c, ∀c st c = c , z i c + z i c ≤ 1 (16)
Constraints ( 13) state that for a given position, there is an incompatibility between a cathode and a given layer if the corresponding metals are different.

Similarly, constraints (14) state that two layers made of distinct metals cannot be sputtered at the same position. Constraints (15) are a consequence of the precedence constraints: they guarantee that two layers belonging to a given product will not be sputtered at the same position on the production line. Constraints ( 16) ensure that two distinct cathodes will not be placed at the same position on the production line.

We next investigate a strengthened formulation for the constraints ( 13)-( 16) based on the analysis of the associated constraint graph and the use of valid inequalities deduced from maximal cliques.

In the constraint graph G = (V, A), a node v ∈ V represents either a type of cathode placed in a given position (i.e a variable z i c ) or a metal layer sputtered in a given position (i.e a variable y i po ). There is an edge a ∈ A between two pairs of nodes if the corresponding variables are linked by one of the binary exclusion constraints ( 13)-( 16). Constraints ( 13)-( 16) all deal with variables related to a same position i on the production line. Therefore, there is no edge in graph G between two nodes corresponding to different positions on the line. In addition, for a given position, the binary variables involved as well as the exclusion relations linking them are identical. G is thus seen to decompose into N independent subgraphs with identical structure:

G i = (V i , A i ) with V i
the subset of nodes related to position i and A i the subset of edges linking these nodes. In the remainder of this subsection, we will study one of these graphs G i for an arbitrary choice of i.

A set C ⊂ V i is called a clique if each pair of nodes in C is connected
by an edge. A maximal clique is a clique which is not properly contained in another clique. Each maximal clique in G i thus gives rise to a valid inequality called a maximal clique constraint stating that the sum of corresponding binary variables should be less than or equal to 1. In the following, those are referred to as type I valid inequalities. We observe that constraints (2) and (4) of the initial formulation are among the clique constraints we can obtain thanks to the study of one of the subgraphs G i . But not all of them are maximal clique constraints. In the sequel, we show how to exploit the special structure of the graphs G i to strengthen these constraints and find other maximal clique constraints, in particular constraints linking variables related to various distinct products.

The structure of a constraint graph G i is close to that of a complete multipartite graph, i.e a graph with node set partitioned into clusters such that any two nodes belonging to different clusters have an edge connecting them and that there is no connection between nodes within a single cluster.

Here, the set of nodes V i can be divided into M + 1 subsets V i m which will be referred to as clusters. The cluster V i m ⊂ V i is the subset of nodes in G i related to metal m. There is an additional cluster, denoted V i 0 , made of a single node, namely the node related to the variable z i 0 . Thanks to constraints (13), ( 14) and ( 15), there is an edge between any two nodes belonging to different clusters. But, because of the constraints ( 15) and ( 16), there are some additional edges between nodes belonging to the same cluster, namely in the cases where the corresponding variables are related to two layers of the same product or to two distinct types of cathodes.

The clusters V i m can be seen to have a special structure. Namely let π(m) ⊂ {1, ..P } be the set of indices of products p using metal m in at least one layer. The cluster V i m is made of 1 + |π(m)| disjoint cliques (possibly containing a single node):

-K 0 m , the clique containing the nodes related to cathode types c ∈ C(m), -for each p ∈ π(m), K p m , the clique containing the nodes related to the layers of product p made of metal m.

Proposition 1 is a direct consequence of this particular graph structure.

Proposition 1. A maximal clique in G i is the union of M + 1 cliques, each clique belonging to a different cluster V i m of the graph.

Figure 1 illustrates the structure of a graph G i for problem P0 introduced in section 2.2. We show a maximal clique containing 6 nodes and yielding the valid inequality: z i 0 +z i 2 +y i 1,1 +y i 1,4 +y i 2,3 +y i 3,2 ≤ 1. For the sake of simplicity, only the edges linking the nodes of this maximal clique are displayed.

Using proposition 1, we can compute the number of maximal cliques in a graph G i as the number of possible combinations obtained by choosing in

each cluster V i m one clique out of (1 + |π(m)|) cliques:
Proposition 2. The number of maximal cliques in a graph G i is given by:

m=1..M (1 + |π(m)|).
This number can be quite high: e.g for the industrial problem P20 (see appendix), G i has 2880 maximal cliques so that 2880 * 30 = 86400 type I valid inequalities should be added to the formulation. Due to this large number, all maximal cliques constraints cannot be added a priori to the model. They can, however, be generated as needed according to a cutting-plane strategy.

In order to do this, we need to address the so-called separation problem.

The separation problem here can be stated as follows: "given (z * , y * , x * ) the optimal solution of the linear relaxation of the problem, find a violated type I valid inequality or decide that (z * , y * , x * ) satisfies all type I valid inequalities ". To solve this problem, we use the following separation algorithm:

(SEP1) Given (z * , y * , x * ) the optimal solution of (1)-( 10 found. In the sequel, algorithm (SEP1) is used to generate type I violated inequalities in order to strengthen the initial formulation.

Valid inequalities from precedence constraints between layers

We now focus on another subset of constraints in our problem: the precedence constraints between layers of a given product. As in subsection 3.2, we exploit the special structure of these binary exclusion constraints to derive a family of stronger valid inequalities corresponding to maximal clique constraints and to further strengthen the formulation.

We first explain how this family of stronger valid inequalities is derived.

We have two families of binary exclusion constraints related to a single product p: constraints (7) of the original model stated in section 2 and valid inequalities (15) stated in the previous subsection. We define the corresponding constraint graph G p = (V p , A p ). A node v ∈ V p refers to a binary variable y i po and can thus be indexed by (i, o). There is an edge a ∈ A p between two nodes of V p if there is a binary exclusion constraint (7) or (15) linking the corresponding variables. Figure 2 shows the graph G 1 obtained for product p = 1 of problem P0. Only a fraction of edges is presented, the edges drawn as dotted lines connect the nodes belonging to a maximal clique.

Maximal cliques in graphs G p have special features that can be exploited as shown by the following result: Proposition 3. A maximal clique in G p consists of exactly O p nodes, each node related to a different layer of product p.

Proof. No two nodes in V p related to the same layer are connected so that the cardinality of a clique in G p cannot be greater than O p , the number of layers of product p.

In addition, a clique in G p cannot be maximal if it does not include a node related to each layer 1...O p . Namely, suppose K is a clique containing O p -1 nodes. All nodes relate to a different layer so that all layers 1, 2, ...O p except layer o are present. K is thus the subset of nodes

K = {(i 1 , 1), (i 2 , 2), ..., (i o-1 , o -1), (i o+1 , o + 1), ..., (i O p , O p )} with i 1 ≥ i 2 ≥ ... ≥ i o-1 ≥ i o+1 ≥ ... ≥ i Op .
We now show that K cannot be a maximal clique. Consider a node (i o , o)

such that i o-1 ≥ i o ≥ i o+1 .
(i ω , ω). So K {(i o , o)} is a clique containing K : K is not a maximal clique of G p .
Each maximal clique in G p provides a valid inequality for our problem.

These valid inequalities will be referred to as type II valid inequalities.

We can compute the number of maximal cliques in a graph G p by induction, as stated below: Proposition 4. Let G p (N, L) be the graph for a product p made of L layers and a production line with N positions. We denote by µ(N, L) the number of maximal cliques in G p (N, L). We have:

(i) ∀N, µ(N, 1) = N (ii) ∀N, ∀L ≥ 2, µ(N, L) = N i=1 µ(N -i + 1, L -1) (The proof is left to the reader).
The number of maximal cliques in G p grows very fast with the problem size, in particular with the number N of positions. With the recurrence given above, the reader can easily check that e.g. for the product p = 5 in problem P20, there are more than 38 billion type II valid inequalities. Hence it is not possible to include directly all type II valid inequalities in the formulation. This is why we propose two ways of using them to strengthen the formulation.

First, we remove constraints (7) from the formulation and replace them by a much smaller number of type II valid inequalities. This involves finding a minimal subset of type II valid inequalities such that every binary exclusion constraint of type ( 7) is implied by at least one valid inequality belonging to this subset, i.e. to find a minimal subset of maximal cliques in G p such that each edge a ∈ A p is covered by at least one maximal clique belonging to this subset. The following heuristic procedure (REP) was devised in order to achieve this. It is based on the idea that the edges of graph G p can be covered in a systematic way by relying on the angle they make with the horizontal axis. More precisely, for each node (i, O p ) in V p , we generate the As shown by the computational experiments to be presented in section 4, the use of procedure (REP) results in a substantial reduction on the total number of constraints in the model as well as in an enhancement of the formulation.

Second, in order to further strengthen the formulation, we generate additional type II valid inequalities according to a cutting-plane strategy. This involves solving the following separation problem for type II valid inequalities: "given (z * , y * , x * ) the optimal solution of the linear relaxation of the problem, find a type II violated valid inequality or decide that (z * , y * , x * ) satisfies all type II valid inequalities". In order to solve it, we will make use of proposition 5 below. We first build the oriented graph Gp = ( Ṽp , Ãp ) in which nodes in Ṽp correspond to binary variables y i po and are indexed by (i, o). There is an oriented arc from node (i , o ) to node (i, o) if i ≤ i and o = o + 1. We define a path as a sequence of nodes linked by arcs directed from a node to the following one. A maximal path is a path which is not contained in another path. Figure 3 shows the graph G1 obtained for product 1 in problem P0. Only a fraction of all arcs is presented. different layer. P is a subset of nodes in Ṽp :

P = {(i 1 , 1)..., (i o , o)..., (i Op , O p )} with i 1 ≥ ... ≥ i o ≥ ... ≥ i O p .
Thanks to proposition 3, we know that the corresponding subset of nodes in V p is a maximal clique of G p . Thus a maximal path of Gp corresponds to a maximal clique in G p . The converse is straightforward.

Thanks to proposition 5, solving the separation problem for type II valid inequalities reduces to the solution of a number of longest path problems in an acyclic graph, leading to the following separation algorithm: (SEP2)

Given (z * , y * , x * ) the optimal solution of ( 1)-( 10), for p = 1...P : 1. Assign to each node (i, o) in Gp a weight equal to the value of the corresponding variable y i po . 2. Find the maximal weight path P max in the acyclic oriented graph Gp with a standard longest path algorithm.

3. Let W max be the weight of P max .

-if W max > 1, the valid inequality given by the maximal clique K max corresponding to the path P max is violated, -else all type II valid inequalities are satisfied for product p.

If W max ≤ 1 for all products p = 1...P , then all type II valid inequalities are satisfied, otherwise we have found at least one violated inequality. In the ILOG (2002) for more details).

We used an industrial data set available in [START_REF] Miegeville | Supply Chain Optimization in the process industry. Methods and Case Study of the Glass Industry[END_REF] to build 20 test problems. P20 is the industrial problem presented in [START_REF] Miegeville | Supply Chain Optimization in the process industry. Methods and Case Study of the Glass Industry[END_REF] and described in table 7 (see Appendix). P0 is the simple example introduced in section 2, P1 to P19 are simpler versions of P20. These instances were obtained by using one or several of the following simplifications: removal of possible positions along the production line; removal of some products; removal of some layers; removal of a metal; removal of some cathode types.

Table 3 displays the following information for each problem tested: the number N of possible positions along the line; the total number p O p of layers to be sputtered; the number C of cathode types; the number Var of binary variables and the number Const of constraints in the initial formulation. The computational results obtained with the initial and enhanced formulations are displayed in table 4. For both series of results, we provide:

-Const: the number of constraints in the formulation. For the enhanced formulation, this is the value obtained after applying the procedure (REP).

-Gap 0 : the initial gap, i.e. the relative difference between the lower bound provided by the linear relaxation of the problem and the best integer solution found after at most 8 hours of computation. For the enhanced formulation, we use the value obtained after the strengthening procedure has stopped.

-Nodes: the number of nodes of the search tree explored before the optimal solution is found or the computation time limit of 8 hours is reached.

-CP U IP : the time in seconds required to find the optimal integer solution when it has been found.

-Gap: the gap obtained after at most 8 hours of computation between the best integer solution found and the best lower bound found.

For the enhanced formulation, we also provide:

-CutsI and CutsII : the number of type I and type II cuts added to the formulation during the strengthening procedure, -CP U cuts : the time in seconds spent to find the type I and type II violated inequalities.

As can been seen from table 4 (columns 2-6), using the initial formulation, only 7 of the 21 problems can be solved exactly within the computational limits. Despite long computation times (8 hours), non-optimal integer solutions are found for 13 problems and in these cases, the remaining gaps obtained remain quite large (16% on average). In addition, no feasible integer solution can be found for problem P20.

We compare these results with the ones obtained while using the enhanced formulation to solve the problem. The results from table 4 (columns 7-14)

show that computation times for small instances are decreased and that more instances (11 out of 21 problems) are solved exactly. In addition, using the enhanced formulation, a feasible integer solution is found for all test problems and, in case the optimal integer solution could not be found after 8 hours of computation, the remaining gap is smaller (9.6 % on average).

Comparison between the results obtained with the two formulations thus

shows that the enhanced formulation improves the efficiency of the Branch & Bound procedure. The main explanatory factor for this is that the lower bounds provided by the linear relaxation of the enhanced formulation (

Influence of cathode capacity

We carried out some additional numerical tests to evaluate the influence of cathode capacity on the algorithmic performance. We considered problems P1 to P5 described in table 3 and we modified the data relative to the cathodes. More precisely, we considered only one type of cathodes per metal and we built instances with various cathode capacity values:

-infinite capacity, -large capacity: for each metal, the available cathode is the cathode with the largest volume among those described in table 7, -medium capacity: for each metal, the available cathode is the cathode with the second largest volume among those described in table 7, -small capacity: for each metal, the available cathode is the cathode with the third largest volume among those described in table 7.

We used the enhanced formulation to solve these instances. Table 5 displays the computational results. We provide Gap 0 , Gap, CP U IP and Nodes as defined in subsection 4.1 and Opt the number of instances that could be The symbol " #" indicates that a guaranteed optimal solution could not be found within 8 hours of computation. solved to optimality within the computation limit. These results suggest that instances with medium or small capacity cathodes are more difficult to solve than instances with infinite or large capacity cathodes. Namely, all instances using infinite or large capacity could be solved to optimality within 2 hours of computation whereas only 2 out of the 10 instances using medium or small capacity cathodes could be solved to optimality within 8 hours of computation. Moreover no feasible solution could be found for 2 out of the 5 instances using small capacity cathodes.

Influence of product composition

We finally discuss the results of some experiments carried out to evaluate the influence of product composition, i.e. of the sequence of metal layers to be deposited on the glass sheets. We built 15 instances involving M = 5 metals, P = 5 products made of 6 layers, N = 20 positions on the line, C = 5 infinite capacity cathodes. They differ only with respect to the sequence of metal layers:

-In E1 to E5, there is a basic sequence of metal defined by product 1.

Products 2 to 5 are obtained by a simple modification of this sequence (switch between two consecutive layers or modification of the metal for one layer).

-In R1 to R5, the sequences of metal layers are randomly generated from a discrete uniform DU (1, 5) distribution. If two consecutive layers are made of the same metal, we repeat the random generation until a product is obtained without any identical consecutive layers.

-In H1 to H5, the sequence of layers for each product are chosen in order to obtain supposedly difficult instances (products made of reverse sequences of metal, products made of sequences with no common metal...).

In order to compare the generated instances, we introduce a measure aiming at evaluating the difference between the products of a given instance with respect to the sequence of metal layers. This difference denoted d is defined as:

d = P p 1 =1 P p 2 =p 1 +1 d(p 1 , p 2 ) where d(p 1 , p 2 ) = SCS(p 1 , p 2 ) - LCS(p 1 , p 2 ). SCS(p 1 , p 2
) is defined as the minimum number of cathodes needed to sputter products p 1 and p 2 and LCS(p 1 , p 2 ) is the maximum number of cathodes that can be used to sputter layers from both p 1 and p 2 . SCS(p 1 , p 2 ) and LCS(p 1 , p 2 ) can be computed by a dynamic programming algorithm as respectively the Shortest Common Supersequence containing p 1 and p 2 and the Longest Common Subsequence contained in p 1 and p 2 . We used the enhanced formulation to solve these instances. The computational results are displayed in table 6. These results suggest that instances with a large value of d are more difficult to solve than instances with a small value of d. Namely, all instances E1-E5 could be solved to optimality within one hour of computation whereas the mean computation time for the instances R1-R5 and H1-H5 is above 4.5 hours. Moreover no feasible solution could be found for 2 out of the 5 instances H1-H5. It is worth pointing out that for the instance P20 presented in Appendix d = 5.8. This seems to indicate that in a industrial situation, the products to be made on the glass coating line are quite different with respect to the sequence of metal layers to be deposited, leading to an additional difficulty to solve the problem.

Conclusion and perspectives

In this paper, we studied an optimization problem arising in the context of the glass industry in connection with a specific transformation of flat glass called glass coating. In order to improve an initial MIP formulation, three families of valid inequalities have been discussed: valid inequalities from limited capacity constraints; valid inequalities from metal compatibility constraints (type I valid inequalities); valid inequalities from precedence constraints between layers of a given product (type II valid inequalities). The results of our computational experiments confirm the positive impact of the proposed enhancements on the computation times and solution quality.

Among the possible research directions suggested by the present work, it might be worth exploring other optimization criteria such as minimizing the volume of unused metal remaining in the cathodes at the end of the production run. Indeed, partially consumed cathodes at the end of a production run represent a cost, either as a direct loss because of the unused metal or as additional constraints for the forthcoming production run because they will impose the use of a set of initial reduced capacity cathodes. Looking for other families of valid inequalities in order to further improve the formulation might also be an interesting research direction.

Appendix

Figure 1 :

 1 Figure 1: Constraint graph G i for problem P0

Figure 2 :

 2 Figure 2: Constraint graph 1 for product p = 1 of problem P0: the edges drawn as dotted lines connect nodes belonging to a maximal clique.

  maximal cliques made up by the edges forming an angle α with the horizontal axis such that tan(α) = b a where a = 1...N -i and b = 1...O p -1. (REP) 1. For p = 1...P , for i = 1...N , for a = 1...N -i, for b = 1...O p -1, generate the following type II valid inequality: i, k, a) = max(i + a(k -1); N ) and O(O p , k, b) = min(O p -bk + 1; 1). 2. Check wether each binary exclusion constraint of type (7) is covered by at least one generated type II valid inequality. If an uncovered binary exclusion constraint is found corresponding to layers o and o < o and positions i and i > i, generate the following type II valid inequality: each generated type II valid inequalities, check wether all binary exclusion constraints it replaces are covered by more than one type II valid inequalities. If so, eliminate the corresponding type II valid inequality.

Proposition 5 .Figure 3 :

 53 Figure 3: Oriented graph G 1 for product p = 1 of problem P0: the arcs drawn as dotted lines are the arcs connecting nodes belonging to the maximal path corresponding to the maximal clique shown in figure 2.

Table 1 :

 1 Problem P0: data on products and available cathodes

	Product p = 1	o	1	2	3	4
	Metal m po	Ag	Au	Ti	Ag
	Volume V po 210 400 100 100
	Product p = 2	o	1	2	3
	Metal m po	Ag	Ti	Au
	Volume V po 410 800 200
	Product p = 3	o	1	2	3
	Metal m po	Au	Pt	Ti
	Volume V po 4000 1000 1000
	Cathode c	1	2	3	4
	Number ν c	5	5	5	5
	Metal m c	Ag	Ti	Au	Pt
	Volume V c 1000 2000 3000 2000

Table 2 :

 2

Problem P0: optimal solution Line Volume sputtered for the layer (p, o)

Table 3

 3

				: Test problems
	N	p O p C Var Const Remarks
	P0 12	10	4 180	1117 small example in section 2.2
	P1 20	17	10 560	8826 P20 without products 4 and 5
	P2 20

  table4 column 11) appear to be stronger than the ones provided by the linear relaxation of the initial formulation (table 4 column 3). Indeed, the integrality gap (i.e the relative difference between Z LP and Z IP ) is reduced on average from around 22% with the initial formulation to about 7.1% with the enhanced formulation. Moreover, it is worth pointing out here that the results provided in table 4 strongly suggest that the automatic cutting-plane generation procedures embedded in the CPLEX software do not seem able to identify the type I and type II valid inequalities exhibited and discussed in section 3.

Table 4 :

 4 Comparison of the initial and the enhanced formulations: results

		Gap		0	0	0	0	0	0	5.9	5.9	16.6	5.3	9.6	0
		Gap 0 Nodes CP U IP	(s)	14.5 46 1.2	0 933 393	13.3 53131 8334	7.1 42141 10027	6.7 112190 25344	0 3942 1500	5.9 39543 #	5.9 38352 #	16.6 30136 #	5.3 26011 #	10.5 24531 #	0 12355 5604
	Enhanced formulation	cuts Const CutsI CutsII CP U	(s)	0.2 777 12 16	72 3100 129 604	30 3465 114 174	122 3634 116 461	108 3999 157 292	48 3846 131 59	716 6342 187 1110	686 6660 226 700	660 6884 184 632	179 6660 189 74	644 7454 214 264	275 6817 227 299	5700 193 364
		Gap		0	0	0	0	14.5	0	22.2	21.5	21.2	19.4	25.6	13.8	5.88
	Initial formulation	Const Gap 0 Nodes CP U IP	(s)	P0 1117 35.2 351 1.8	P1 8826 23.9 10547 688	P2 11378 30.3 143329 10501	P3 11799 31.6 163415 13706	P4 14351 30.8 240214 #	P5 13170 21.5 133718 17456	P6 20333 25.6 119944 #	P7 22484 22.9 89052 #	P8 24335 23.3 92031 #	P9 22484 22.2 98344 #	P10 27137 26.7 73411 #	P11 20108 13.8 105293 #	P12 13004 11.1 266188 #

Table 5 :

 5 Influence of cathodes capacityCathode capacity Gap 0 Opt Gap CP U IP (s) #Nodes

	infinite	11.4	5	0	2617	8781
	large	5.4	5	0	10930	30983
	medium	11.9	1	8.8	> 28800	105480
	small	2.6	1	2.6	> 28800	39088

Table 6 :

 6 Influence of product compositionInstances d Gap 0 #Opt Gap CP U IP (s) #Nodes

	E1-E5	2.8 22.4	5	0	2184	7405
	R1-R5	5.5 13.1	4	6.5	14642	32325
	H1-H5	6.4 11.1	3	0	19075	29792

Table 7 :

 7 Problem P20: data on products and cathodes

	p = 1	o	1	2	3	4
			m po Ag	Au	Ti	Ag
			V po 2200 2400 2000 2000
	p = 2	o	1	2	3	4	5	6	7
			m po Ag	Au	Ti	St	Au	Ti	Ag
			V po 2100 1300 1000 1000 700 2000 4000
	p = 3	o	1	2	3	4	5	6
			m po Au	Pt	Ti	St	Au	Pt
			V po 2000 1000 1000 2400 1000 2000
	p = 4	o	1	2	3	4	5	6	7
			m po	St	Ti	St	Au	Pt	Ti	St
			V po 1500 1000 2400 1000 750 500 1000
	p = 5	o	1	2	3	4	5	6	7	8
			m po Au	Pt	Ti	St	Ag	St	Ti	St
			V po 1000 750 500 1000 2000 1500 1000 2400
	c	1	2	3	4	5	6	7	8	9	10
	ν c	10	10	10	10	10	10	10	10	10	10
	m c Ag Ag Ag	Ag	Ag	Ag	Ti	Ti	Ti	Ti
	V c 300 500 1000 2000 3000 4000 4000 3000 2500 1000
	c	11	12	13	14	15	16	17	18	19	20
	ν c	10	10	10	10	10	10	10	10	10	10
	m c	Ti	Au Au	Au	Au	Au	St	St	St	St
	V c 400 100 500 1000 2500 3500 500 750 1000 1500
	c	21	22	23	24	25	26
	ν c	10	10	10	10	10	10
	m c	St	Pt	Pt	Pt	Pt	Pt
	V c 2000 500 750 1000 1500 2000

sequel, algorithm (SEP2) is used to generate type II violated inequalities in order to strengthen the initial formulation.

Computational results

In this section, we discuss the results of the computational experiments carried out to evaluate the impact of the formulation enhancements presented in section 3. We also present an empirical study carried out to show the impact of some problem parameters on the algorithmic performance.

Comparison of the initial and enhanced formulations

In order to evaluate the impact of the proposed formulation enhancements, we solved the problem with a standard MIP software (CPLEX 8.1.0) using either the initial formulation described in section 2 or the enhanced formulation. More precisely, the strengthened formulation is obtained thanks to the following procedure:

1. We use procedure (REP) to replace precedence constraints of type ( 7) by a subset of type II valid inequalities and we add the M valid inequalities (12) to the formulation. We solve the linear relaxation of the problem.

3. We use the separation algorithm (SEP1) to add type I violated valid inequalities.

4. When no more type I violated valid inequalities can be found, we look for type II violated valid inequalities using the separation algorithm (SEP2).

5. When no more type II violated valid inequalities can be found, we go back to step 3 and repeat until no more violated valid inequalities (whether of type I or type II) can be generated. All the tests were run on a Pentium 4 (2.8 GHz) with 504 Mo of RAM, running under Windows XP. We used the default settings of CPLEX solver. This means that some cutting planes, among which are clique cuts, cover cuts and Gomory fractional cuts, are added automatically to the model (see


Céline Gicquel
email: celine.gicquel@ecp.fr

Nicolas Miègeville

Michel Minoux

Yves Dallery

N M Ièg

Ev Ille

V Ig Nes

F S Ch Atenay -M Alabry

R Aint-Gobain

Ech Erch

Q Uai Lucien-Lefranc

Discrete lot sizing and scheduling using product decom position into attrib utes

Keywords: P ro d uctio n P la n n in g, L o t S iz in g, M ix ed

HAL is

Introduction

A wide variety of models for production planning and inventory management h as b een investigated. Among th em, capacitated lot sizing models aim at determining a minimal cost production sch edule complying with given capacity restrictions and such th at demand for all products is satisfi ed with out b ack logging. Lot sizing prob lems can b e b roadly classifi ed b y th e numb er of products, capacitated resources and product levels involved.

T h e present work pertains to multi-item single-resource single-level lot-sizing models: multiple items assumed to b e independent to one anoth er (i.e. not link ed b y parentcomponent relationsh ips) are to b e produced on a single constrained resource. Moreover all prob lem parameters are supposed to b e deterministically k nown.

A furth er important distinction among deterministic lot-sizing models can b e made b ased on considerations on th e demand (see e.g. [START_REF] Kuik | Batching decisions: structure and models[END_REF]):

• O n one h and, th ere are models such as th e E conomic Lot Sch eduling Prob lem (E LSP) wh ere th e demand rate is supposed to b e constant and a continuous time scale and infi nite time h orizon are used. E arly contrib utions on th e E LSP can b e found for instance in [START_REF] Bomberger | A dynamic programming approach to a lot size scheduling problem[END_REF] and [START_REF] Elmaghraby | The Economic Lot Scheduling problem (ELSP): review and extensions[END_REF].

• O n th e oth er h and, some models assume a dynamic, time-varying demand and use a discrete time scale and fi nite time h orizon. In such models, time-dependent demand is specifi ed on a period-b y-period b asis. R ecent overviews on th e literature ab out th is line of research can b e found in [START_REF] Drexl | Lot sizing and scheduling -Survey and extensions[END_REF], [START_REF] Degraeve | Meta-heuristics for dynamic lot sizing: a review and comparison of solution approaches[END_REF], [START_REF] Degraeve | Modelling industrial lot sizing problems: a review[END_REF] and [3 2] .

In th e present paper, we focus on a capacitated dynamic lot-sizing prob lem k nown as th e Discrete Lot sizing and Sch eduling Prob lem (DLSP). T h e DLSP is b ased on several k ey assumptions (see e.g. [START_REF] Fleischmann | The discrete lot sizing and scheduling problem[END_REF]):

• Demand for products is deterministically k nown and time-varying.

• T h e production plan is estab lish ed for a fi nite time h orizon sub divided into several discrete periods.

• At most one item can b e produced per period (" small b uck et" model) and th e facility processes eith er one product at full capacity or is completely idle (" all-ornoth ing" assumption).

• Costs to b e minimized are th e inventory h olding costs and th e ch angeover costs.

In lot sizing prob lems, th e ch angeover costs to b e incurred wh en th e production of a new lot b egins can depend eith er on th e next item only (sequence-independent case)

or on b oth th e previous and th e next items (sequence-dependent case). W e consider here the (more difficult) case of sequence-dependent changeover costs. In this paper, we assume zero changeover times. H owever the integration of positive changeover times is an important extension of our model and is currently being worked on.

The reader is referred to [START_REF] Zhu | Scheduling and lot sizing with sequence-dependent setup: A literature review[END_REF] for a recent literature review on scheduling and lot sizing with sequence-dependent setup. In the sequel, we discuss more specifically contributions on the DLSP with sequence-dependent changeover costs. In [START_REF] Fleischmann | The discrete lot sizing and scheduling problem with sequencedependent set-up costs[END_REF] and [START_REF] Salomon | Solving the discrete lotsizing and scheduling problem with sequence dependent set-up costs and set-up times using the Travelling Salesman Problem with time windows[END_REF], the problem is reformulated as a Travelling Salesman Problem with Time Windows and is solved using either Lagrangian relaxation or a dynamic programming-based algorithm. The authors of [START_REF] Jordan | Discrete Lotsizing and Scheduling by Batch Sequencing[END_REF] show the equivalence between the DLSP with sequence-dependent changeover costs and the B atch Sequencing Problem (B SP) and use a specific branch-and-bound type algorithm for solving the B SP to optimality. In these papers, the number of items considered in the computational experiments is relatively small (no more than 10 items) whereas the horizon length can be up to 10 0 periods. In a recent paper, [START_REF] Wolsey | Solving multi-item lot-sizing problems with an MIP solver using classification and reformulation[END_REF] proposes to strengthen a basic MIP (mixed-integer programming) formulation of the DLSP with sequence-dependent changeover costs using both a reformulation of the changeover variables and valid inequalities. Thanks to this strengthened formulation, the lower bounds provided by the linear relaxation of the problem are significantly better, enabling a B ranch & B ound type procedure to solve the problem more efficiently. H owever, as pointed out by [START_REF] Belvaux | Modelling practical lot-sizing problems as mixed-integer programs[END_REF], the large number of variables needed in the reformulation to handle changeovers is an important drawback of this approach.

The purpose of the present paper is to propose a new way of modelling the production system to be planned by properly exploiting a multi attribute product structure frequently encountered in industrial applications. The basic idea of the new model originates from the observation that products can usually be described in terms of a set of physical attributes such as color, dimension, quality level... If this is possible, each item to be produced will be identified, not only by a unique index as it is usually done, but also by a M -tuple, each component of which indicates the value of the corresponding attribute for the given item. When such a multi attribute product structure can be exhibited in the industrial context under study, we propose to exploit it to reduce the size of the mixed integer program to be solved. This is achieved by looking at changeovers at an aggregate level using the relevant physical attributes instead of considering each individual changeover between items. B y doing so, we are able to significantly reduce the number of changeover variables and associated constraints in the MIP formulation, while maintaining the quality of the bounds provided by the linear relaxation of the problem.

We extend the approach used by [START_REF] Wolsey | Solving multi-item lot-sizing problems with an MIP solver using classification and reformulation[END_REF] to derive valid inequalities for the resulting mixedinteger linear program. Computational results show that exact optimal solutions can in general be obtained more efficiently with the new model as compared with previously described approaches.

The paper is organized as follows. In section 2, we introduce the strengthened reformulation proposed by [START_REF] Wolsey | Solving multi-item lot-sizing problems with an MIP solver using classification and reformulation[END_REF] for the DLSP with sequence-dependent changeover costs. In our computational experiments, we use it as a reference for comparison with our model.

In section 3, we describe our proposal to model the production system using product attributes and discuss our main assumptions. When these assumptions hold, we propose in section 4 to exploit them to formulate the DLSP with sequence-dependent changeover costs with a reduced number of variables and constraints. Results of computational experiments carried out on a large number of randomly generated instances to validate our approach are reported in section 5 . Section 6 presents the concluding remarks and discussions for future investigations.

A strong formulation for the DLSP with sequencedep endent changeov er costs

In this section, we recall a strong formulation for the DLSP with sequence-dependent changeover costs. This formulation was presented by [START_REF] Karmarkar | The deterministic dynamic product cycling problem[END_REF] for the variant of the DLSP referred to as CSLP (Continuous Setup Lot sizing Problem), where the all-or-nothing assumption is relaxed. More recently, [START_REF] Belvaux | Modelling practical lot-sizing problems as mixed-integer programs[END_REF] and [START_REF] Wolsey | Solving multi-item lot-sizing problems with an MIP solver using classification and reformulation[END_REF] propose to use it to solve the DLSP with sequence-dependent changeover costs. We next discuss a further strengthening of this formulation obtained by exploiting valid inequalities developed by [START_REF] Van Eijl | On the discrete lot-sizing and scheduling problem with Wagner-Whitin costs[END_REF].

We wish to optimize the production schedule for a set of N items over an horizon featuring T planning periods. A period is indexed by t = 1, ..., T , an item by k = 0, .., N .

We agree to use item k = 0 to represent idle periods.

We use the following notation for the parameters:

• D kt : demand (in units) for item k in period t,

• P kt : production capacity (in units per period) for item k in period t,

• h k : holding costs per unit and period for item k,

• c kl : changeover costs from item k to item l.

Decision variables are defined as follows:

• I kt : inventory level corresponding to item k at the end of period t,

• y kt : setup variables. y kt = 1 if the resource is setup for item k in period t, and 0 otherwise,

• w klt : changeover variables. w klt = 1 if the resource is switched from item k to item l at the beginning of period t, and 0 otherwise.

With this notation, [START_REF] Belvaux | Modelling practical lot-sizing problems as mixed-integer programs[END_REF] propose to formulate the DLSP with sequence-dependent changeover costs as follows:

(DLSP1)

m in N k=1 T t=1 h k I kt + N k=0 N l=0 T t=1 c kl w klt (1)
∀k, ∀t,

I kt = I k,t-1 + P kt y kt -D kt (2)
N k=0

y k0 = 1 (3) ∀k, ∀t, y k,t-1 = N l=0 w klt (4) ∀l, ∀t, y lt = N k=0 w klt (5)
∀k, ∀l, ∀t, w klt ≥ 0 (6) ∀k, ∀t, I kt ≥ 0 [START_REF] Elmaghraby | The Economic Lot Scheduling problem (ELSP): review and extensions[END_REF] ∀k, ∀t, y kt ∈ { 0, 1}

The objective, minimizing the sum of inventory holding costs and changeover costs, is expressed by [START_REF] Belvaux | Modelling practical lot-sizing problems as mixed-integer programs[END_REF]. Changeover costs c kl are incurred between two successive production batches of item k and item l, in the first period of production of item l. Constraints [START_REF] Bomberger | A dynamic programming approach to a lot size scheduling problem[END_REF] express the inventory balance. The "all-or-nothing" assumption is enforced by the term P kt y kt in the equality: if the resource is setup for k in period t, then all the available capacity is used and the production quantity of item k must be equal to P kt . ( 3) is also linked to the "all-or-nothing" assumption: together with constraints (4)-( 5), they ensure that in each period, the resource either produces a single product at full capacity, or is idle (i.e y 0t = 1). Equalities (4) and ( 5) link the setup variables with the changeover variables. (4) guarantee that item k can be produced in period t -1 if and only if a changeover from k to another item l (possibly l = k) takes place at the beginning of period t. Similarly, (5) guarantee that item l can be produced in period t if and only if a changeover from another item k (possibly k = l) to item l takes place at the beginning of period t. (6) state the non-negativity of the changeover variables: observe, as pointed out by [START_REF] Belvaux | Modelling practical lot-sizing problems as mixed-integer programs[END_REF], that thanks to constraints (3)-( 5) and ( 8), there is no need to define variables w klt as binary variables. The set of constraints ( 2) and [START_REF] Elmaghraby | The Economic Lot Scheduling problem (ELSP): review and extensions[END_REF] ensure that demand for each item is fulfilled without backlogging. The binary character of the setup variables is represented by [START_REF] Fleischmann | The discrete lot sizing and scheduling problem[END_REF].

As suggested by [START_REF] Wolsey | Solving multi-item lot-sizing problems with an MIP solver using classification and reformulation[END_REF], the formulation DLSP1 can be further strengthened using a family of strong valid inequalities developed by [START_REF] Van Eijl | On the discrete lot-sizing and scheduling problem with Wagner-Whitin costs[END_REF] for the single-item DLSP with Wagner-Whitin costs, constant capacity and no backlogging. U nder these assumptions, demand can be measured in terms of how many units can be produced during one production period, i.e. the production capacity and demand quantity can be normalized to one unit per period without loss of generality: d kt ∈ {0, 1} and P kt = 1

In order to present the valid inequalities developed by [START_REF] Van Eijl | On the discrete lot-sizing and scheduling problem with Wagner-Whitin costs[END_REF], we first introduce some additional notation:

• Dk,t,τ : cumulated demand for item k in the interval {t, ..., τ }. Thanks to the normalization, demand on item k is binary so that Dk,t,τ is equal to the number of positive demand periods for k in {t, ..., τ }.

• S k,q : q th positive demand period for item k. Note that S k, Dk,1,t +q denotes the q th period with positive demand after period t.

We also introduce the start-up variables z kt defined as:

z kt =       
1 if the production of a new lot of item k begins at period t, i.e. if a start-up for item k takes place at the beginning of period t, 0 otherwise.

These variables are linked to the changeover variables by the equations:

∀k, ∀t, z kt = l:l =k w lkt (9)
With this notation, the following inequations [START_REF] Gao | A comparative study of modeling and solution approaches for the coordinated lot-size problem with dynamic demand[END_REF] are valid inequalities for the DLSP with sequence-dependent changeover costs:

∀t, ∀k, ∀p ∈ {0... Dk,t+1,T }, I kt ≥ p q=1 1 -y k,t+q - S k, Dk,1,t +q τ =t+q+1 z kτ (10)
We briefl y explain the underlying idea. F irst note that y k,t+q + S k, Dk,1,t +q τ =t+q+1 z kτ = 0 if and only if the resource is not setup for item k in period t + q and no startup for k takes place between the period t+q+1 and the period where the q th demand after period t occurs, i.e. if and only if no production of item k is possible in the interval {t + q, ..., S k, Dk,1,t +q }. In this case, the quantity needed to satisfy the q th demand after period t should be in stock at the end of period t. Thus we see that constraints [START_REF] Gao | A comparative study of modeling and solution approaches for the coordinated lot-size problem with dynamic demand[END_REF] force an increase of the stock of item k at the end of period t by one for each index q for which no production occurs in the interval {t + q, ..., S k, Dk,1,t +q }. The reader is referred to [START_REF] Van Eijl | On the discrete lot-sizing and scheduling problem with Wagner-Whitin costs[END_REF] for a detailed proof of the validity of [START_REF] Gao | A comparative study of modeling and solution approaches for the coordinated lot-size problem with dynamic demand[END_REF]. In our computational experiments, we add these valid inequalities to the formulation according to a cutting-plane generation strategy (for details see section 4.2). The resulting strengthened formulation is denoted DLSP1* .

As pointed out by [START_REF] Belvaux | Modelling practical lot-sizing problems as mixed-integer programs[END_REF], an important drawback of the DLSP1 formulation is that the number of variables needed in the formulation to handle changeovers, (N + 1) 2 T , grows very rapidly with problem size. In the sequel, we present one way to avoid this issue in certain situations, namely when products can be described as combinations of a number of physical attributes.

3

The DLSP with products described as combina-

tions of physical attributes

In most papers dealing with lot sizing problems, each individual item to be produced is described with a single index (k in the formulation presented above) and is considered independently of the other items. However, in many industrial situations, the items to be produced diff er in their physical characteristics or attributes (e.g. color, dimensions, shape, mixture composition, quality level...). As explained by [START_REF] Oh | Planning production on a single processor with sequencedependent setup part 1: determination of campaigns[END_REF] for the case of a metals plant, such physical attributes often contribute significantly to the sequence dependency of changeover costs. Moreover it is frequently the case that many items share a common value for some attribute so that we can define a (small) finite number of possible values for each attribute. In the sequel, we propose to describe each item to be produced as well as the resource setup state using the physical attributes of the products and to look at changeovers at an aggregate level using the relevant physical attributes instead of considering each individual changeover between items. By doing so, we are able to simplify the production planning process with respect to:

• the evaluation of changeover costs. Changeovers costs are most often estimated by human production experts. In the presence of a multi attribute product structure, it may be easier for them to evaluate changeover costs attribute by attribute rather than item by item. This was experienced by the authors in the industrial application described in [START_REF] Miègeville | Supply Chain Optimization in the process industry. Methods and Case Study of the Glass Industry[END_REF].

• the resolution of the optimization problem. Thanks to the use of a multi-attribute product structure, we are able to reduce the size of the obtained mixed integer program and thus to make it easier to solve with a standard MIP software.

The idea of exploiting a multi attribute product structure may be used in a variety of modelling environments. We focus here on its application on a particular type of lot sizing problems, namely the DLSP, but the present work could be extended to deal with a closely related dynamic lot sizing problem known as the General Lot sizing and Scheduling Problem or GLSP (see e.g. [START_REF] Meyr | Simultaneous lotsizing and scheduling by combining local search with dual reoptimization[END_REF]).

In section 4, we derive a new formulation for the DLSP which is likely to be solved more efficiently using standard MIP software. In order to do this, we need to make several assumptions on the production system.

Main assumptions

1. We suppose that each item to be produced can be described by a set of physical attributes, each of them takes a finite number of discrete values. We also suppose that each item is uniquely identified thanks to a M -tuple, each component of which gives the value of the corresponding attribute for the given item.

2. We assume that the setup state of the resource can also be described using product attributes. Thus, we will not describe the setup state of the resource by indicating the item that the resource is able to produce, but by indicating, for each attribute, for which value of this attribute the resource is setup. The resource setup state will therefore also be described by a M -tuple, each component of which gives the value of the corresponding attribute for the present state of the resource. To ensure consistency, it should be understood that a given item can be produced on the resource if and only if the resource is setup with the correct value for every attribute.

3. We assume that we are able to evaluate the changeover costs on the resource for each attribute separately . This means that given an attribute and two possible values for this attribute, we are able to evaluate the cost of a changeover from one value to the other and that this cost does not depend on the setup state of the resource with respect to the other attributes.

4. We need to specify how the costs relative to different attributes will combine, i.e. how we will compute changeover costs when changeovers for different attributes happen simultaneously on the resource. In this paper, we consider the case where the global changeover costs is the sum of all individual changeover costs for the different attributes. Another possible assumption is that global changeover costs equal the max imum of the individual costs.

Thanks to these assumptions, we will be able to decide about the production plan on the resource using the product attributes. In this case, the production plan will consist of a set of parallel sequences, one for each attribute. Each of these sequences indicates, for every planning period, for which value of the corresponding attribute the resource is setup. Thus, in each planning period, combining the values for the different attributes, we will be able to deduce the item for which the resource is setup. A detailed mathematical programming formulation is proposed in section 4, but in order to illustrate the usefulness of the new model, we first discuss some industrial situations where it appears to be well suited.

Possible industrial applications

In order to illustrate the practical relevance of the proposed model, we provide examples of industrial situations found in the literature where using physical attributes to describe the products is appropriate.

• In [START_REF] Dobson | The cyclic lot scheduling problem with sequence-dependent setups[END_REF] and [START_REF] Smith-Daniels | A mixed integer programming model for lot sizing and sequencing packaging lines in the process industries[END_REF], a production planning problem for a bottle filling line is considered.

For this type of production line, two physical attributes of the products have to be taken into account: the size or shape of the bottle and the liquid used to fill it.

Hence each item can be described by means of two attributes: the bottle size/ shape and the liquid to be used. Each individual item would be described by a pair (i 1 , i 2 )

where i 1 is the index of the corresponding bottle size/ shape and i 2 the index of the corresponding liquid.

• [START_REF] Dos Santos-Meza | A lot-sizing problem in an automated foundry[END_REF] discuss a lot sizing problem they found in an automated foundry. Each item to be produced can be described by two attributes: the type of metal alloy it is made of and the shape it takes from the used mould. Here we could use as well a pair (i 1 , i 2 ) where i 1 would give the index of the alloy type and i 2 the index of the mould shape.

• [START_REF] Silva | Heuristic lot size scheduling on unrelated parallel machines with applications in the textile industry[END_REF] study a production planning problem arising in the textile industry in a company producing acrylic fibers. The authors report that two physical characteristics of the products have an impact on the scheduling of the plant spinning unit, namely the fiber composition and their diameter. Thus, we could use a pair (i 1 , i 2 ) to describe each item: i 1 would refer to the fiber composition and i 2 to its diameter.

• [START_REF] Miègeville | Supply Chain Optimization in the process industry. Methods and Case Study of the Glass Industry[END_REF] considers the production planning problem for a float glass production line.

Here each item (a glass sheet) can be described using several physical characteristics: glass color and quality, dimensions of the sheet (thickness, width and length).

An item could thus be described using a 5-tuple, with components corresponding to the color, quality and dimensions of the corresponding glass sheet. In fact, the pro-duction system model studied in the present paper was suggested in this particular context of application.

Though far from being exhaustive, the above list is a good indication of the wide applicability of the model proposed here. Moreover the proposed model could also prove useful for the frequently encountered situation where changeover costs have a majorminor setup structure. This type of structure is mainly found when the items to be produced can be grouped into families. In this case, there is a large changeover involved with setting up the resource to produce a given family of items (major setup) and an additional smaller changeover is required between two items belonging to the same family (minor setup). This situation is described for instance in [START_REF] Gao | A comparative study of modeling and solution approaches for the coordinated lot-size problem with dynamic demand[END_REF], [START_REF] Ho | A new MRP/ GT lot sizing heuristic: a simulation study[END_REF], or [START_REF] Ep | Coordinated capacitated lot-sizing problem with dynamic demand: a lagrangian heuristic[END_REF] for the case of a time-varying demand (i.e. for dynamic lot sizing problems such as the DLSP) as well as in [START_REF] Dobson | The cyclic lot scheduling problem with sequence-dependent setups[END_REF], [START_REF] Mohamed | An efficient model for multifamily lot-sizing and scheduling: application to a real life problem[END_REF], [START_REF] Oh | Planning production on a single processor with sequencedependent setup part 1: determination of campaigns[END_REF] or [START_REF] Wagner | A search heuristic for the sequence-dependent economic lot scheduling problem[END_REF] for the case of a constant demand (i.e. for the ELSP). Such a major-minor setup structure could be modelled using a multi attribute product structure.

Namely, each item could be described using two attributes: the first one corresponding to the product family and the second one being used to identify the individual item within the family.

Our model is also capable of handling the situations where there is a change of magnitude within the values of a given physical attribute m, i.e. when the indices i ∈ [0, V m ] correspond to values of m belonging to different orders of magnitude. Namely, in this case, the corresponding changeover costs are likely to be sequence-dependent: there will be a large changeover between two values of m belonging to different orders of magnitude and a smaller changeover between two values belonging to the same order of magnitude.

In our model, we allow sequence-dependent changeover costs between values of a given physical attribute and are thus able to represent such situations without making any additional assumptions.

4

A formulation for the DLSP with product attributes and sequence-dependent changeover costs

B asic formulation

We now present a formulation for the DLSP with product attributes and sequencedependent changeover costs. This formulation can be used to solve the DLSP when a product description using physical attributes is possible and when the assumptions discussed in section 3.1 hold.

We use the same notation as in section 2 for the parameters relative to items:

• D kt : demand for item k in t,

• P kt : production capacity for item k in t,

• h k : holding costs per unit and period for k.

We assume that each item can be described using M physical characteristics or attributes.

Correspondence between products and attributes is given by a matrix A of dimensions M × (N + 1). A mk represents the value of the attribute m for product k and the k th column of A gives the M -tuple describing product k in terms of product attributes. For each attribute m, we have:

• a set of possible values: i ∈ [0, V m ] • a changeover cost matrix: C m . C m ij is the cost of a transition from the value i ∈ [0, V m ] to the value j ∈ [0, V m ] of attribute m.
We agree to use the M -tuple (0,0,....,0) to describe the item k = 0: i.e ∀m, A m0 = 0.

We use the following decision variables:

• I kt : inventory level corresponding to item k at the end of period t,

• y kt : setup variables at the item level. y kt = 1 if the resource is setup for item k in period t, and 0 otherwise,

• w m ijt : changeover variables at the attribute level. w m ijt = 1 if a switch from the value i to the value j of attribute m takes place at the beginning of period t, and 0 otherwise.

Under the assumption that changeover costs related to different attributes are added whenever two transitions occur simultaneously (see assumption 4 in section 3.1), the DLSP can be formulated as follows:

(DLSP2) min N k=1 T t=1 h k I kt + M m=1 V m i=0 V m j=0 T t=1 C m ij w m ijt (11) ∀k, ∀t, I kt = I k,t-1 + P kt y kt -D kt (12)
N k=0

y k0 = 1 (13) ∀m, ∀i ∈ [0, V m ], ∀t, k st A mk =i y k,t-1 = V m j=0 w m ijt (14) ∀m, ∀j ∈ [0, V m ], ∀t, k st A mk =j y kt = V m i=0 w m ijt (15) ∀m, ∀(i, j) ∈ [0, V m ] × [0, V m ], ∀t, w m ijt ≥ 0 (16) ∀k, ∀t, I kt ≥ 0 (17) ∀k, ∀t, y kt ∈ {0, 1} (18)
The objective, minimizing the sum of changeover costs and inventory holding costs, is expressed by [START_REF] Ho | A new MRP/ GT lot sizing heuristic: a simulation study[END_REF]. Note that inventory holding costs are computed item by item whereas changeover costs are computed attribute by attribute. Constraints ( 12) express the inventory balance. Combined with the non negativity constraints [START_REF] Meyr | Simultaneous lotsizing and scheduling by combining local search with dual reoptimization[END_REF], they prevent any backlogging. ( 13), together with constraints ( 14)-( 15), guarantee that in each period the resource either produces a single item or is idle.

Equalities ( 14) and ( 15) link the setup variables with the changeover variables. First note that the term k st A mk =i y kt equals 1 if and only if an item k requiring the resource to be setup for the value i of the attribute m is produced in period t, i.e. if and only if the resource is setup for the value i of attribute m in period t. Thus ( 14) guarantee that the resource is setup for the value i of attribute m in period t-1 if and only if a changeover from value i to another possible value j of attribute m (possibly j = i) takes place at the beginning of period t. Similarly, [START_REF] Karmarkar | The deterministic dynamic product cycling problem[END_REF] guarantee that the resource is setup for the value j of attribute m in period t if and only if a changeover from another possible value i of attribute m (possibly i = j) to value j takes place at the beginning of period t. The non negativity of the changeover variables is stated by ( 16) and the binary character of the setup variables is expressed by [START_REF] Meyr | Simultaneous lotsizing and scheduling on parallel machines[END_REF].

The formulation DLSP2 can be easily modified to consider the other possible assumption about the combination of costs relative to different attributes, i.e. the assumption that global changeover costs equal the maximum of the individual costs. This can be done by defining additional continuous variables C t to evaluate the changeover costs to be incurred at the beginning of each period t. In this case, the DLSP can be formulated as follows:

(DLSP2 MAX ) min N k=1 T t=1 h k I kt + T t=1 C t s.t. ∀m, ∀t, C t ≥ V m i=0 V m j=0 C m ij w m ijt
and ( 12) -(18)

In the sequel, we assume that global changeover costs equal the sum of the individual costs and thus use the formulation DLSP2. However similar results could be obtained with the other assumption.

Let us now compare the number of changeover variables in the formulations DLSP1 and DLSP2. For the sake of simplicity, we do not take into account the item k = 0 in the comparison. As shown in section 2, in this case, DLSP1 includes N 2 T changeover variables, one for each possible pair of items and for each period. Note that when the product description using attributes is possible, we can compute the number of products as the number of possible combinations obtained by choosing for each attribute m one value out of V m . Thus we have:

N 2 = M m=1 V m 2
. Now, as can be seen above, in th e D L S P 2 form u lation, th ere are M m=1 V m2 T ch ang eover variables, one for each p air of p ossible valu es of each attribu te and for each p eriod . In m ost cases wh ere th e p rod u ct d escrip tion u sing attribu tes will be im p lem ented , we will h ave:

M m=1 V m2 M m=1 V m 2
, thus leading to a significant reduction in the number of changeov er v ariables needed in the formulation.

In word s, in th e p rop osed m od el, we d o not consid er each ind ivid u al ch ang eover between item s, bu t rath er look at ch ang eovers at a m ore ag g reg ate level u sing p rod u ct attribu tes. B y d oing so, we are able to sig nifi cantly red u ce th e siz e of th e m ix ed integ er linear p rog ram m to be solved (e.g . u sing a B ranch & B ou nd p roced u re).

Strengthening the formulation with valid inequalities

A s for th e D L S P 1 form u lation, th e D L S P 2 form u lation can be fu rth er streng th ened u nd er th e assu m p tion of W ag ner-W h itin costs, constant cap acity and no back log g ing . T h is can be ach ieved by ex tend ing th e ineq u alities (1 0 ) to th e D L S P 2 form u lation. In ord er to d o th is, we d efi ne two new sets of variables:

Y m it = k st A mk =i y kt =       
1 if th e resou rce is setu p for th e valu e i of attribu t m in p eriod t, 0 oth erwise.

Z m it = j∈[0 ,V m ] st j =i w m jit =       
1 if a startu p for th e valu e i of attribu t m tak es p lace at th e beg inning of p eriod t, 0 oth erwise.

W ith th is notation, we h ave:

Proposition 1. All feasible solutions of DLSP2 satisfy:

∀t, ∀k , ∀p ∈ {0... Dk,t+1,T }, ∀m = 1...M ,

I kt ≥ p q=1 1 -Y m A mk ,t+q - S k, Dk,1,t +q τ =t+q+1 Z m A mk ,τ (19 )
Proof. Before the proof, which extends the one given in [START_REF] Van Eijl | On the discrete lot-sizing and scheduling problem with Wagner-Whitin costs[END_REF] for the DLSP1 formulation, we briefl y explain the idea underlying [START_REF] Miègeville | Supply Chain Optimization in the process industry. Methods and Case Study of the Glass Industry[END_REF]. Y m A mk ,t+q + S k, Dk,1,t +q τ =t+q+1 Z A mk = 0 if and only if the resource is not setup in period t+ q for the value A mk of attribute m needed to produce item k and no startup for this value occurs between the period t + q + 1 and the period where the q th demand after period t occurs, i.e. if and only if no production of item k is possible in the interval {t + q, ..., S k, Dk,1,t +q }. In that case, the quantity needed to satisfy the q th demand on item k after period t should be in stock at the end of period t.

Now consider an arbitrary integral feasible solution of DLSP2, say (y, w, Y, Z). We arbitrarily choose an item k, a period t, a demand occurrence p ∈ {0... Dk,t+1,T } and an attribute m and we show that the feasible solution chosen satisfies the corresponding valid inequality. In the sequel, for the sake of simplicity, we drop the item index k and we denote i = A mk the value of attribute m for item k.

We denote R q the q th production period for this item in the feasible solution considered. By definition, we have R 1 < R 2 < ... < R q < ... < R D1,T . M oreover, because backlogging is not allowed, the q th production period must occur before the q th demand period: ∀q, R q ≤ S q . Let q 0 be the highest index such that R D1,t +q < t + q. Then we have:

• ∀q ≤ q 0 , R D1,t +q < t + q. The q th demand after period t is produced before period t+ q .

• ∀q > q 0 , t + q ≤ R D1,t +q ≤ S D1,t +q . The q th demand after period t is produced between t+ q and the period S D1,t +q where it occurs. In this case, the resource must be setup for the value i of attribute m at least once in the interval {t+q, ..., S D1,t +q }.

Thus we have:

∀q > q 0 , Y m i,t+q + S D1,t +q τ =t+q+1 Z m i,τ ≥ 1 (20)
Hence,

t t=1 y t + p q=1 Y m i,t+q + S D1,t +q τ =t+q+1 Z m i,τ ≥ t t=1 y t + q 0 q=1 Y m i,t+q + p q=q 0 Y m i,t+q + S D1,t +q τ =t+q+1 Z m i,τ (21)
≥ t t=1 y t + q 0 q=1 y t+q + p -q 0 (22)
≥ D1,t + q 0 + pq 0 (23)

≥ D1,t + p (24)
(21) comes from the fact that q 0 q=1 S D1,t +q τ =t+q+1 Z m i,τ ≥ 0. To obtain [START_REF] Özdamar | Hybrid heuristics for the capacitated lot sizing and loading problem with setup times and overtime decisions[END_REF], we use the fact that y t+q ≤ Y m i,t+q as well as the inequalities [START_REF] Mohamed | An efficient model for multifamily lot-sizing and scheduling: application to a real life problem[END_REF]. F inally, ( 23) is true because, by definition of q 0 , the cumulated demand D1,t + q 0 is satisfied by the cumulated production before t + q 0 , t+q 0 t=1 y t , so that t+q 0 t=1 y t ≥ D1,t + q 0 .

As t t=1 y t -D1,t is the inventory level of item k at the end of period t, this establishes the validity of [START_REF] Miègeville | Supply Chain Optimization in the process industry. Methods and Case Study of the Glass Industry[END_REF].

The number of valid inequalities [START_REF] Gao | A comparative study of modeling and solution approaches for the coordinated lot-size problem with dynamic demand[END_REF] and [START_REF] Miègeville | Supply Chain Optimization in the process industry. Methods and Case Study of the Glass Industry[END_REF] grows quite fast with the problem size and the production capacity utilization: e.g. for the instances involving 30 products and 100 periods with a capacity utilization of 90% (see section 5 ), there are more than 7 000 valid inequalities [START_REF] Gao | A comparative study of modeling and solution approaches for the coordinated lot-size problem with dynamic demand[END_REF] for the DLSP1 formulation and more than 21000 valid inequalities [START_REF] Miègeville | Supply Chain Optimization in the process industry. Methods and Case Study of the Glass Industry[END_REF] for the DLSP2 formulation. Hence it is not possible to include directly all valid inequalities in the formulations DLSP1 and DLSP2. In the computational experiments to be presented in section 5 , the following cutting-plane generation strategy has been implemented to strengthen both formulations:

1. We solve the linear relaxation of the problem using the formulation DLSP1 (respectively DLSP2).

2. We look at each valid inequality of type [START_REF] Gao | A comparative study of modeling and solution approaches for the coordinated lot-size problem with dynamic demand[END_REF] (respectively [START_REF] Miègeville | Supply Chain Optimization in the process industry. Methods and Case Study of the Glass Industry[END_REF]). If it is violated by the current continuous solution, we add it to the formulation.

3. If at least one violated inequality is found in step 2, we go back to step 1 and repeat until no more violated valid inequalities can be generated.

The resulting strengthened formulation is denoted DLSP1* (respectively DLSP2* ).

A small illustrative example

We use a small example to illustrate the proposed model and to show an application of the DLSP2 formulation. We consider a bottle filling line where 4 items can be produced. An item is described by the corresponding bottle size (attribute 1 with two possible values) and the composition of the liquid to be used (attribute 2 with two possible values). Table 1 shows how each of the 4 items can be described using the two attributes. We agree to use the item k = 0 described by the 2-tuple (0,0) to represent the idle period. Table 2 gives the changeover costs for each attribute and table 3 provides the demand for each product.

Figure 1 shows the optimal production plan obtained while using the DLSP2 formulation. The first two lines give the sequence of setup states for each attribute. In each planning period, we can deduce from these sequences the item for which the resource is setup. The changeover costs to be incurred between each lot are also shown in the bottom line of figure 1. C onsistently with our assumption, changeover costs relative to diff erent attributes are added whenever transitions for diff erent attributes occur simultaneously. This is the case here at the beginning of periods 1, 4, 9 and 10 where both the bottle size and the liquid composition are changed.

C omputational results

In this section, we discuss the results of some computational experiments carried out to compare the two formulations presented in sections 2 and 4. We created 5 sets of randomly generated instances. The instances diff er with respect to the following characteristics:

• Problem dimension: The problem dimension is represented by the number of products N and the number of periods T. We use three diff erent combinations:

(N , T ) ∈ {(10, 6 0), (25, 50), (30, 100)}.

• Multi-attribute p roduct structure: The product structure is described by the number of attributes M and the number of possible values V m for each attribute m. We use five diff erent combinations, leading to 5 sets of instances. Table 4 gives the characteristics of the generated instances for each set. For set C and E instances, we have M m=1 V m > N . Therefore we used the following procedure to generate matrix A:

1. We generated a matrix A with M m=1 V m columns. A describes all possible combinations of the attribute values.

2.

For each column i of A , we randomly generated a weight w i from a discrete uniform DU(1, M m=1 V m ) distribution. 3. The matrix A is generated by selecting the N columns of A corresponding to the N smallest weight w i .

• Inventory holding costs: For each item, inventory holding costs have been generated randomly from a discrete uniform DU(5, 10) distribution.

• Production capacity utilization: Production capacity utilization ρ is defined as the ratio of the total cumulated demand on the total cumulated available capacity. We experimented different values for ρ: 0.5, 0.7 and 0.9.

• Demand pattern: Binary demands for each item have been randomly generated according to the procedure described in [START_REF] Salomon | Solving the discrete lotsizing and scheduling problem with sequence dependent set-up costs and set-up times using the Travelling Salesman Problem with time windows[END_REF].

• C hangeover costs: For each attribute m, changeover costs C m lk have been randomly generated from a discrete uniform DU(C m min , C m ma x ) distribution. We tested several possibilities: the changeover costs for all attributes can either be taken from the same interval or the changeover costs for the first attribute are greater than for the other(s). In our study, we define the ratio r as: r = For each possible combination of multi-attribute product structure, production capacity utilization and changeover costs ratio, 5 problems were generated, resulting in 5 × 3 × 5 × 5 = 375 instances. All tests were run on a Pentium 4 (2.8 G hz) with 505 Mb of R AM, running under Windows X P. We used a standard MIP software (CPLEX 8 .1.0) with the solver default settings, using either formulation DLSP1 or formulation DLSP2.

Tables 56789show the computational results obtained with formulations DLSP1* and DLSP2*, for each set of instances. As the value of the ratio r appears to have an impact on the results quality, we grouped the instances with respect to the value of r so that each line corresponds to the average value for 15 randomly generated instances (5 instances for each value of production capacity utilization). For both series of results, we provide:

• # O pt: for set A, B and C instances, the number of instances out of the corresponding 15 instances that could be solved to optimality within 30 minutes of computation.

• # F eas: for set D and E instances, the number of instances out of the corresponding 15 instances for which a feasible solution could be found within 30 minutes of computation.

• Gap: for the instances that could not be solved to optimality, the average relative gap value obtained after 30 minutes of computation between the best integer solution (if one could be found) and the best lower bound found.

We now compare the results obtained with formulations DLSP1* and DLSP2*. The results from tables 5-9 show that:

• for high values of ratio r (r ≥ 5), i.e. when one attribute has corresponding changeover costs clearly higher than the other(s) attribute(s), the results obtained with formulation DLSP2* are better. This can be seen as:

a feasible solution could be obtained for all instances, -more instances could be solved to optimality within 30 minutes of computation, -when a guaranteed optimal solution could not be found within 30 minutes of computation, the residual gap is smaller.

• for small values of the ratio r (r ≤ 2), formulation DLSP1* provides better results for medium-sized instances (sets A, B and C). However, this is not the case for the larger instances in sets D and E. Namely, for these instances, -a feasible solution could not always be found with formulation DLSP1* whereas at least one feasible solution could be found for each instance with formulation DLSP2*.

the residual gap is significantly smaller on some instances with formulation DLSP2*.

Comparison between the results obtained with the two formulations thus shows that using formulation DLSP2*, we are able to improve the effi ciency of the Branch & Bound procedure, especially for the high values of ratio r and for the largest instances. This can be explained by two main factors:

• U sing formulation DLSP2*, the problem size (i.e. the number of variables and constraints) is significantly reduced. As a consequence, the time spent at each node of the Branch & Bound tree to solve the linear relaxation is shorter and more nodes can be explored within 30 minutes of computation.

• The formulation enhancement obtained thanks to the valid inequalities adapted for formulation DLSP2 gives better results when ratio r has a high value. More precisely, for high values of r, the lower bounds provided by formulation DLSP2* are higher than the ones provided by formulation DLSP1*. O n the contrary, for small values of r, the lower bounds provided by the formulation DLSP1* are higher than the ones provided by formulation DLSP2*.

Thus the combined advantages of a reduced problem size and of tighter lower bounds enable the formulation DLSP2* to outperform the formulation DLSP1* on many instances.

Conclusion

We presented a new formulation for the DLSP with sequence-dependent setup costs.

The main idea is to use a possible product description as combinations of a number of physical attributes. When such a structure is present in the industrial context under study, we have shown how to exploit it to reduce the size of the mixed integer linear program to be solved and improve the efficiency of the solution process. Computational experiments show that the proposed DLSP2 formulation performs better than the tight DLSP1 formulation we chose as a reference for comparison, especially in cases where one of the physical attributes has corresponding changeover costs higher than the other(s) attribute(s).

Finally, future research on this topic could be aimed at extending the proposed MIP formulation DLSP2* to take into consideration additional relevant industrial concerns such as the presence of multiple resources that need to be planned simultaneously (see e.g. [START_REF] Dastidar | Scheduling injection molding operations with multiple resource constraints and sequence dependent setup times and costs[END_REF], [START_REF] Meyr | Simultaneous lotsizing and scheduling on parallel machines[END_REF], [START_REF] Özdamar | Hybrid heuristics for the capacitated lot sizing and loading problem with setup times and overtime decisions[END_REF] or [START_REF] Segerstedt | A capacity-constrained multi-level inventory and production control problem[END_REF]) or the integration of positive changeover times (see e.g. [START_REF] Jordan | Discrete Lotsizing and Scheduling by Batch Sequencing[END_REF], [START_REF] Salomon | Solving the discrete lotsizing and scheduling problem with sequence dependent set-up costs and set-up times using the Travelling Salesman Problem with time windows[END_REF], [START_REF] Trigeiro | Capacitated lot sizing with setup times[END_REF]or [START_REF] Vanderbeck | Lot-sizing with start-up times[END_REF]). It would also be worth investigating whether a multi attribute product structure could be used to solve lot sizing problems such as the ELSP where different modelling assumptions are used.

  denotes the mean of interval [C m min , C m ma x ]. We tested several values for r : 1, 2, 5, 10 and 30. In all instances, the resulting changeover costs between two items belong to the interval [0,200].