François Larroche
email: francois.larroche@imt-atlantique.fr

Odile Bellenguez

Guillaume Massonnet

Clustering-based solution approach for a capacitated lot-sizing problem on parallel machines with sequence-dependent setups

Keywords: Lot-sizing, sequence-dependent setups, parallel machines, clustering

Clustering-based solution approach for a capacitated lot-sizing problem on parallel machines with sequence-dependent setups

François Larroche, Odile Bellenguez, Guillaume Massonnet

Introduction

In this paper, we focus on production problems that arise in practical cases from the food industry and involve the manufacturing of several items on multiple production lines. The main goal of the industrial application behind this study is to reach "good" feasible solutions in a short computational time, providing the planners with a decision support tool that they can use on a daily basis. Unfortunately, those problems generally combine multiple lot-sizing and scheduling constraints, which makes them difficult to handle using a straightforward Mixed Integer Linear Programming (MILP) modeling. In particular, when the setup times between different types of items heavily depend on the production sequence, such approaches lead to large mathematical formulations that turn out to be intractable by commercial solvers. Thus, fast and reliable dedicated methods are necessary in order to meet practitioners' requirements.

We call this problem the multi-item capacitated lot-sizing problem with lost sales, safety stock, overtimes, and sequence-dependent setups on parallel machines, abbreviated CLSSD-PM in the remainder of this paper. It extends many existing lot-sizing models that have been developed in the production planning literature following the work of [START_REF] Wagner | Dynamic Version of the Economic Lot Size Model[END_REF]. The ever-growing need to propose customized models that allow to precisely describe and solve a wide range of industrial production processes has shifted the interest of manufacturers towards different generalizations of the original problem, see for instance [START_REF] Jans | Modeling industrial lot sizing problems: a review[END_REF] for a comprehensive survey of these applications. The version that we consider in this paper regroups some of these well-known extensions, along with new ones specific to our industrial case. In particular, it includes several aspects related to the Capacitated Lot Sizing Problem (CLSP) and the Asymmetric Travelling Salesman Problem (ATSP). Both problems are proven NP-hard (see [START_REF] Florian | Deterministic Production Planning: Algorithms and Complexity[END_REF] and [START_REF] Karp | Reducibility among Combinatorial Problems[END_REF], respectively), hence it is straightforward that the CLSSD-PM is itself NP-hard.

The capacity constraints that appear in the practical applications we consider are similar to the ones found in the CLSP, see [START_REF] Karimi | The capacitated lot sizing problem: a review of models and algorithms[END_REF] or [START_REF] Quadt | Capacitated lot-sizing with extensions: a review[END_REF] for reviews of extensions and solution approaches. Lost sales is another aspect of the CLSSD-PM that frequently appears in the lot-sizing literature. In particular, it is studied by [START_REF] Absi | The multi-item capacitated lot-sizing problem with setup times and shortage costs[END_REF], who introduce new classes of valid inequalities. In [START_REF] Absi | Heuristics for the multi-item capacitated lot-sizing problem with lost sales[END_REF], a Lagrangian heuristic is developed and provides good results for the CLSP with lost sales. On the other hand, the safety stock is seldom considered in the deterministic production and inventory literature. [START_REF] Loparic | The uncapacitated lot-sizing problem with sales and safety stocks[END_REF] define the safety stock as a lower bound on the number of units that must be held in the inventory in each period. However, in many real-life cases, it is mainly seen as a target inventory level, set to mitigate the risk of shortage rather than an hard constraint. Following this approach, [START_REF] Absi | The multi-item capacitated lot-sizing problem with safety stocks and demand shortage costs[END_REF] choose to penalize missing units from the safety stock in the objective function and develop a Lagrangian heuristic approach to solve the problem.

The combination of production planning and scheduling problems is relevant in practice, since setup times are often sequence-dependent. This observation has motivated the introduction of integrated models generally referred to as lot-sizing and scheduling problems, that consider both levels of decision to properly approximate the capacity consumption in each time period. For instance, [START_REF] Gicquel | A tight MIP formulation for the Discrete Lot-sizing and Scheduling problem with parallel resources[END_REF] present a discrete lot-sizing and scheduling problem (DLSP) based on a small bucket formulation, i.e. the production in each period is limited to only one type of item, with the additional constraint that any production must be at full capacity. In contrast, the so-called big bucket formulations consider that time periods are long enough to produce multiple types of items. [START_REF] Fleischmann | The general lotsizing and scheduling problem[END_REF] are among the first to study a sequence-dependent version of this problem, called the General Lot-sizing and Scheduling Problem (GLSP), and consider two time structures, namely a set of fixed periods (big buckets) that contain sequences of smaller variable production periods (small buckets). For a state of art of these different variants of the lot-sizing and scheduling problem, the reader can refer to the extensive reviews of [START_REF] Guimarães | Modeling lotsizing and scheduling problems with sequence dependent setups[END_REF] and [START_REF] Copil | Simultaneous lotsizing and scheduling problems: a classification and review of models[END_REF].

In most cases, complex lot-sizing problems with capacity restriction, sequencedependent (SD) setups and/or multiple machines are tackled using decomposition and constructive heuristics to obtain feasible solutions. [START_REF] Almada-Lobo | Single machine multi-product capacitated lot sizing with sequence-dependent setups[END_REF] propose an iterative heuristic to solve the capacitated lot-sizing problem with sequencedependent setup times and costs (CLSD) and compare its result with a lower bound from a strengthened formulation. The well-known Relax-and-Fix is also a popular decomposition heuristic that usually offers a good trade-off between solution quality and runtime. It is frequently combined with the local search heuristic Fix-and-Optimize to improve the solutions. We refer the reader to [START_REF] Absi | MIP-based heuristics for multi-item capacitated lot-sizing problem with setup times and shortage costs[END_REF] for a thorough presentation of this class of methods, including an extensive parameter analysis. [START_REF] Toso | Lot sizing and sequencing optimisation at an animal-feed plant[END_REF] apply this type of procedure to the GLSP, while [START_REF] Clark | Production setup-sequencing and lot-sizing at an animal nutrition plant through atsp subtour elimination and patching[END_REF] propose a solution procedure based on subtour elimination and a patching heuristic. [START_REF] Lang | Fix-and-optimize heuristics for capacitated lot-sizing with sequence-dependent setups and substitutions[END_REF] use a Fix-and-Optimize heuristic to address the single machine version of CLSD including overtimes. [START_REF] Fiorotto | Impact analysis of setup carryover and crossover on lot sizing problems[END_REF] apply a Relax-and-Fix and Fix-and-Optimize procedure to obtain feasible solutions for difficult instances of a single machine problem with sequence-independent setup times including crossover (i.e. start a setup and end it in the next period) and carryover (i.e. continue the production of an item from one period to the next without incurring a new setup).

In addition, problems that include multiple resources are divided into two broad categories: identical [START_REF] Beraldi | Rollinghorizon and fix-and-relax heuristics for the parallel machine lot-sizing and scheduling prob-lem with sequence-dependent set-up costs[END_REF]) and non-identical [START_REF] Kaczmarczyk | Modelling Set-up Times Overlapping Two Periods in the Proportional Lot-Sizing Problem with Identical Parallel Machines[END_REF]) parallel machines. The problem CLSSD-PM that we address considers non-identical production machines. It extends the case of sequence-dependent setup times and costs on parallel machines presented in [START_REF] James | Single and parallel machine capacitated lotsizing and scheduling: New iterative MIP-based neighborhood search heuristics[END_REF], for which the authors develop new MIP-based neighborhood search heuristics and compare them to other solution procedures. Several methods have been experimented for similar frameworks. [START_REF] Almada-Lobo | Multiple machine continuous setup lotsizing with sequence-dependent setups[END_REF] consider a version of the problem in which a given machine can only produce one type of item per period and propose to combine the addition of valid inequalities with a Lagrangian heuristic to obtain good solutions. [START_REF] Mateus | Capacitated lot sizing and sequence dependent setup scheduling: an iterative approach for integration[END_REF] solve an integrated version of the problem with an iterative procedure, while [START_REF] Xiao | MIP-based fix-andoptimise algorithms for the parallel machine capacitated lot-sizing and scheduling problem[END_REF] develop an original Fix-and-Optimize algorithm to address a problem from the semiconductor manufacturing industry. Other works such as Fiorotto and de Araujo (2014) and [START_REF] Xiao | A hybrid Lagrangian-simulated annealing-based heuristic for the parallel-machine capacitated lot-sizing and scheduling problem with sequence-dependent setup times[END_REF] rely on Lagrangian decomposition heuristics to deal with parallel machines problems. More recently, Armas and Laguna (2020) study a similar problem, from the pipe-insulation industry, and propose a solution procedure that runs in two phases: In the first one, a simplified version of the problem with sequence-independent setups is solved, before re-scheduling the production using a post-processing heuristic in the second phase.

The remainder of this paper is organized as follows: Section 2 describes the notations and assumptions of the CLSSD-PM problem, including a presentation of the industrial context that explains the introduction of production sequence. In Section 3, we propose two different mathematical formulations for our problem, respectively referred to as the aggregate and the facility location ones. In Section 4, we succinctly presents the so-called Relax-and-Fix and Fix-and-Optimize heuristics adapted to the CLSSD-PM problem in order to obtain feasible solutions in a short computational time on large instances. Then we develop an approach based on a clustering procedure to approximate the sequence-dependent setup times. The created clusters are used to reduce the problem size and obtain good solutions within a limited computational time. In the following, we introduce a two-phases approach: The first one aims at finding a solution to the modified clustered-problem, which is then used in the second one to build optimal production sequences for each machine in each period. This approach ensures that any solution to the modified problem remains feasible in the original framework. Finally, Section 5 presents a sensitivity analysis for parameters tuning and computational results apply on instances inspired by practical data sets from our industrial partners. Numerical results suggest that even for medium size instances, our approximation of the setup times enables us to quickly find better solutions to the CLSSD-PM, regardless to the type of formulation used.

Problem description

Notations and assumptions

The problem CLSSD-PM is an extension of the CLSP with lost sales and safety stock on (unrelated) parallel machines. The goal is to plan the production of N different items over a discrete finite horizon of T periods, each one considered as a big bucket, on M non-identical parallel machines (or other production resources). For conciseness, we denote N , M and T the set of items, machines and periods, respectively (see Table 1). In each period t ∈ T , there is a deterministic demand d i t for each item i ∈ N , which is either lost or satisfied from on-hand inventory. For notation convenience, we also let D i tt = t u=t+1 d i u be the cumulative demand from period t + 1 to t for all t < t . While it is often advised to minimize the physical stock at all time, production managers usually set a safety stock S i t for each item i and period t in order to absorb uncertain perturbations of the production process (e.g. machine failures). For all i ∈ N and m ∈ M, we let τ i m represent the per-unit production time of item i on machine m. Each machine m ∈ M has a time-dependent production capacity C mt , which corresponds to the expected production time necessary in period t, estimated beforehand by the S&OP (Sales and Operations Planning). However, the planner can choose to exceed this planned capacity for machine m up to a hard limit C mt which is the "true" available time capacity in period t. In that case, C mt -C mt corresponds to the maximum overtime allowed. Finally, we also impose that such a production is greater than a minimum production quantity q i min . For a given machine m, the setup time required to switch from one type of item i to another type j = i depends on the pair (i, j) and is denoted λ ij m . We assume that the sequence-dependent setup times satisfy the triangle inequality.

All the above operations induce costs. Specifically, producing one unit of item i in period t generates a per-unit production cost p i mt . In addition, using machine m for production or setup during period t incurs a cost c mt per time-unit of usage. Whenever an overtime is requested, each extra unit of time needed is charged at an additional cost cmt . Recall that the production planner aims at keeping the stock of each final product as close as possible to its target inventory level. We penalize any deviation from this objective with a holding cost h i+ t that applies to each excess unit of item i over the safety stock S i t at the end of period t. On the other hand, failing to meet this safety stock requirement also induces a penalty cost h i- t for each missing unit. The capacity restrictions that apply to the production process may also lead to some unsatisfied demands. Every unmet unit of demand for item i induces a shortage cost l i t (lost-sales) in period t. The objective is to minimize the sum of inventory holding costs, shortage costs, production costs and line usage costs. The different input data are summarized in Table 2.

Production sequence

We focus on a version of the problem in which the periods are considered as independent. Machines are reset between two consecutive periods, hence no carryover is considered and each machine is regarded as idle at the beginning of each period. Moreover, we assume that the first setup is performed outside of the time slot available for production in each period. This assumption is rare in related works on lot-sizing and scheduling. To the best of our knowledge, it only appears in [START_REF] Clark | Production setup-sequencing and lot-sizing at an animal nutrition plant through atsp subtour elimination and patching[END_REF]. In our case, this is in line with many food-factory practices, in which a team is assigned to the setup of machines before the actual production period begins.

We model this idle state of a machine using a virtual item indexed 0, which is always considered as active, i.e. in production, at the beginning of a period. For notational convenience we also let N 0 = N ∪ {0} and artificially define null setup times between item 0 and any other item i ∈ N , i.e. λ 0i m = 0 for all m ∈ M. One of the main challenges of the CLSSD-PM problem lies in the influence of the production sequence on the setup times. The so-called sequence-dependent setups extension of lot-sizing problem reconciles production planning and scheduling problems, in the sense that the two levels of decision are integrated. In the food-industry, as in many other fields of application, the sequence-dependent setups are determined on the basis of some items characteristics. One can think of the transition from non-organic to organic products which requires a long clean-up to meet the sanitary standards, while clean-up time is negligible in the opposite sequence, as an example of such features likely to influence the setup times. In the CLSD case, the main works that deal with sequence-dependent setups are related to the ATSP problem. Indeed once the production mix is determined in a given period, the optimal production sequence corresponds to the optimal tour for an ATSP on the selected subset of items. Figure 1 presents two different sequences for a problem with 4 items on a single machine. The two solutions are equivalent in terms of quantity produced in the period considered, but differ in their respective cost due to a longer production overtime in the upper sequence. As in the ATSP, the difficult part of this problem largely comes from the subtour elimination constraints. For a quick representation of the different subtour configurations, see [START_REF] Almada-Lobo | Single machine multi-product capacitated lot sizing with sequence-dependent setups[END_REF]. Note that in the case of the CLSSD-PM problem, this issue becomes even more significant with the multiplication of production resources, each requiring to solve an ATSP problem for any subset of items affected to it.

Mathematical formulations

In the remainder of this section, we present two different formulations for the CLSSD-PM. For conciseness, we first introduce decision variables and constraints that are shared by both formulations, before giving a more detailed description for each of them. In all of our models, we use the binary variables w ij mt ∈ {0, 1} to indicate if item i is the immediate predecessor of item j in the production sequence on machine m in period t. We also use U mt and O mt ∈ R + to represent the total production plus setup time and the total overtime of machine m in period t, respectively.

Several constraints, relative to capacity restrictions and subtours elimination, are present in all of our formulations, namely :

-Enforce the total production time of machine m in period t to be lower than the capacity allocated to it, plus the possible overtime fixed by the planner:

U mt ≤ C mt + O mt ∀m ∈ M, ∀t ∈ T (1)
-Prevent the overtime to exceed the difference between the physical (hard) capacity and the planned capacity in a given period:

O mt ≤ C mt -C mt ∀m ∈ M, ∀t ∈ T (2)
-Impose that an item can only be setup if there exists a setup from the neutral state to any real item:

i∈N w 0i mt ≥ i∈N0 w ij mt ∀j ∈ N , ∀m ∈ M, ∀t ∈ T (3)
-Allow a setup from an active item only if it belongs to the production sequence of the period:

j∈N0 w ji mt ≥ j∈N w ij mt ∀i ∈ N , ∀m ∈ M, ∀t ∈ T (4)
-Forbid more than one setup per item:

j∈N w ij mt ≤ 1 ∀i ∈ N 0 , ∀m ∈ M, ∀t ∈ T (5)
-Prevent disconnected subtours: The goal is to ensure that in each period t the sequence variables w ij mt defines a single cycle that starts and ends in the neutral item 0.

Multiple options are available from the literature to satisfy this property: One of the most frequently seen relies on the so-called MTZ constraints, see e.g. [START_REF] James | Single and parallel machine capacitated lotsizing and scheduling: New iterative MIP-based neighborhood search heuristics[END_REF]. More recently, Guimarães, Klabjan, and Almada-Lobo (2014) introduces the Single Commodity Flow (SCF) constraints and provides numerical evidences of their superior performances compared to MTZ.

Formulations based on the SCF constraints turned out to be more efficient in our numerical experiments, hence we focus on the latter in the remainder of the paper. For all i ∈ N 0 ,j ∈ N m ∈ M and t ∈ T , we introduce flow variables f ij mt ∈ R + that represent the flow of one commodity from item i to item j on machine m in period t. The following set of constraints prevents the subtours in the sequence :

j∈N f 0j mt = i∈N0 j∈N w ij mt ∀m ∈ M, ∀t ∈ T (6) i∈N0 f ij mt = i∈N0 w ij mt + i∈N f ji mt ∀j ∈ N , ∀m ∈ M, ∀t ∈ T (7) 0 ≤ f ij mt ≤ N w ij mt ∀i ∈ N 0 , ∀j ∈ N , ∀m ∈ M, ∀t ∈ T (8)
Constraints ( 6) ensure that for each period and machine, the total commodity flow from the source corresponds to the number of items produced in the period on each machine. Constraints (7) define the flow balance. Constraints ( 8) impose an upper bound on the amount of flow. -Define the domain of every variable :

w ij mt ∈ {0, 1} ∀i ∈ N 0 , ∀j ∈ N , ∀m ∈ M, ∀t ∈ T (9) U mt , O mt ≥ 0 ∀m ∈ M, ∀t ∈ T (10)
We are now ready to introduce the two mathematical formulations for the CLSSD-PM.

Aggregate formulation MIP-AGG

We first present an aggregate formulation, which is the most common to describe lot-sizing problems. Specific decision variables for this model are summarized below:

x i mt ∈ R + Quantity of item i produced in period t on machine m L i t ∈ R + Quantity of lost sales for item i in period t I i t ∈ R +
Inventory of item i on hand at the end of period t I i+ t ∈ R + Overstock (based on safety stock value) of item i at the end of period t I i- t ∈ R + Safety stock deficit of item i at the end of period t z i t ∈ {0, 1} Binary variable equals to 1 if the stock of item i is null at the end of period t The aggregate formulation, inspired by [START_REF] Clark | Production setup-sequencing and lot-sizing at an animal nutrition plant through atsp subtour elimination and patching[END_REF], is expressed with the following MIP:

min t∈T m∈M (c mt U mt + c mt O mt ) + t∈T i∈N h i+ t I i+ t + h i- t I i- t + l i t L i t + m∈M p i mt x i mt (11) s.t. U mt = i∈N (τ i m x i mt + j∈N λ ji m w ji mt ) ∀m ∈ M, ∀t ∈ T ( 12
)
I i t = S i t + I i+ t -I i- t ∀i ∈ N , ∀t ∈ T ( 13
)
I i t = I i t-1 + m∈M x i mt + L i t -d i t ∀i ∈ N , ∀t ∈ T ( 14
)
I i t ≤ D i tT (1 -z i t ) ∀i ∈ N , ∀t ∈ T (15)
L i t ≤ z i t d i t ∀i ∈ N , ∀t ∈ T (16) x i mt ≤ D i t-1,T j∈N0 w ji mt ∀i ∈ N , ∀m ∈ M, ∀t ∈ T (17) x i mt ≥ q i min j∈N0 w ji mt ∀i ∈ N , ∀m ∈ M, ∀t ∈ T (18)
Constraints ( 1)-( 10)

I i+ t ≤ D i tT -S i t ∀i ∈ N , ∀t ∈ T (19) I i- t ≤ S i t ∀i ∈ N , ∀t ∈ T (20) L i t , I i t , I i+ t , I i- t ≥ 0 ∀i ∈ N , ∀t ∈ T (21) x i mt ≥ 0 ∀i ∈ N , ∀m ∈ M, ∀t ∈ T (22)
Constraints (12) define the total working time of each machine which is equal to the production time and the setup times. Constraints (13) defines the inventory on hand for item i at the end of period t, which is equal to the safety stock value plus the overstock minus the deficit. Constraints ( 14) are the inventory flow conservation equations through the planning time horizon. Constraints ( 15) and ( 16) ensure that lost sales for item i occur only when the corresponding stock is null. Constraints (17) use the cumulative demand over the remainder of the horizon as an upper bound on the quantity of each item produced on each machine in a given period. Constraints (18) force the production to be higher than its minimum value. Finally constraints ( 19) and ( 22) define the domain of specific variables.

Facility location formulation MIP-FL

The second formulation we present is inspired by the facility location problem. This idea was originally developed in [START_REF] Krarup | Plant location, Set Covering and Economic Lot Size: An 0 (mn)-Algorithm for Structured Problems[END_REF] for the uncapacitated lotsizing problem and adapted later to several of its extensions. It requires to define new variables x it ms ∈ R + that represent the quantity of item i produced in period s on machine m to satisfy a demand in period t. In this formulation, linear costs are written using new parameters denoted by H it ms that aggregate all the costs induced by the production in period s till the demand satisfaction in period t, i.e. the production cost and holding costs of the corresponding units. We introduce additional notations to define the value of H it ms . For a given demand i and period t, we let θ i t ≤ t be the first period in which the safety stock contains units of item i that serve demand d i t . Note that since θ i t is uniquely defined for all pair (i, t), we implicitly assume that the safety stock is defined in such a way that it covers the demand of an integer number of future periods. This effectively turns the safety stock into a safety cover, which is the case in practice. Using this new notations, one can express the inventory costs in relation with the safety stock using only the production period s and the consumption period t. Formally, we define parameters H it ms as follow:

H it ms = p i ms -l i t + θ i t -1 u=s h i+ u - t-1 u=max{θ i t ,s} h i- u ( 23
)
Note that the coefficients in front of the lost-sales and safety stock penalty costs in (23) are negative. This is due to the fact that they both represent a deficit to a target level (demand and safety stock, respectively), which is reduced through the production of units of item i. In contrast, production and excess holding costs both increase with the production of new units, which explains their positive coefficients in (23). Therefore, one simply needs to add the constant i∈N t∈T l i t + t-1 u=θ i t h i- u d i t to compare the solutions values to the ones obtained with the MIP-AGG formulation.

We define the MIP-FL formulation as follows :

min s∈T m∈M c ms U ms + c ms O ms + i∈N T t=s H it ms x it ms (24) s.t U ms = i∈N   τ i m T t=s x it ms + j∈N λ ij m w ij ms   ∀m ∈ M, ∀s ∈ T (25) x it ms ≤ d i t j∈N0 w ji ms ∀i ∈ N , ∀m ∈ M, ∀s ∈ T (26) m∈M t s=1 x it ms ≤ d i t ∀i ∈ N , ∀t ∈ T (27) T t=s x it ms ≥ q i min j∈N0 w ji ms ∀i ∈ N , ∀m ∈ M, ∀s ∈ T (28)
Constraints ( 1)-( 10)

x it ms ≥ 0 ∀i ∈ N , ∀m ∈ M, ∀t ∈ T , ∀s = 1, . . . , t (29)
The objective function ( 24) is still the minimization of the total cost incurred. Constraints ( 25) is just a reformulation of ( 12) with the new variables x it ms . Constraints (26) impose that an item can only be produced if the corresponding machine is setup for it. Constraints ( 27) ensure that it is impossible to produce more that the demand in period t from productions up to that period. Constraints (28) set the lower bound on production quantities. Finally, constraints (29) define the specific variables.

Solution approach

The industrial specifications relative to this problem impose to solve the MIP presented in the previous section in a short execution time. Unfortunately, both formulations are too large and complex to quickly obtain good (or even feasible) solutions on any realistic data set with any of the existing commercial solvers. In this section, we propose alternative solution approaches to tackle this problem. A classical approach in such cases is to use heuristics to obtain good, albeit often suboptimal, solutions in a short computational time. One example of heuristic that is frequently encountered in the lot-sizing literature is the so-called Relax-and-Fix (RF) procedure. In a nutshell, it is a constructive heuristic that builds a feasible solution using a decomposition strategy over the time horizon. It is often combined with a second step, called Fix-and-Optimize (FO), that aims at improving the quality of the initial solution using a local search procedure. Both approaches are quickly introduced and adapted to the CLSSD-PM in section 4.1.

It turns out that a simple adaptation of a RFFO heuristic does not yield satisfactory results in the allocated time on large instances. Since sequence-dependent setup times make the problem significantly harder to solve, we then take advantage of the data structure to overcome these limitations. Specifically, we introduce in section 4.2 a clustering approach that allows us to derive a variation of the original problem. Thus, it becomes easier to solve, through a significant reduction of the number of sequence variables w ij mt and subtour elimination constraints. Combining this approach with the RFFO procedure, we are then able to obtain good solutions to this modified problem in a limited amount of time. In addition, the production plans obtained remain provably feasible for the original CLSSD-PM.

4.1.

A Relax-and-Fix and Fix-and-Optimize procedure

We present two well-known decomposition-based heuristics: The Relax-and-Fix (RF) and the Fix-and-Optimize (FO). Both procedures rely on the resolution of multiple "easier" subproblems. The goal of each subproblem is to set the value of a small subset of integer variables, while the remaining ones are either relaxed or fixed.

Relax-and-Fix with period-based decomposition

The RF heuristic constructs a feasible solution iteratively by solving several subproblems with fewer binary variables. Several types of such decomposition exist in the literature, however we focus on the period-based decomposition which has been widely applied and is generally considered the most efficient method (see [START_REF] Helber | A fix-and-optimize approach for the multi-level capacitated lot sizing problem[END_REF] for a comparison with the product-based decomposition). The basic principle of the procedure consists in moving a decision window of δ periods through the planning horizon. The periods before and after the decision window are said to be frozen and approximate, respectively.

Specifically, during the first iteration the algorithm solves the MIP-AGG (or MIP-FL) problem considering binary variables w ij mt in the decision window (t = 1, . . . , δ) and their relaxed counterpart in the approximate periods, i.e. 0 ≤ w ij mt ≤ 1 for all t = δ + 1, . . . , T . In the second iteration, the decision window is shifted to the right by δσ periods, where σ < δ is the number of overlapping periods in the decision window between two consecutive iterations. The setup variables from the previous approximated periods that enter the decision window are reverted to binary variables. Binary variables that leave the decision window, i.e. whose time index belongs to {1, . . . , δ -σ}, become frozen and keep their value from the previous resolution step. The algorithm repeats the procedure of shifting and fixing variables until the decision window reaches the end of the horizon. Figure 2 is inspired from [START_REF] Absi | MIP-based heuristics for multi-item capacitated lot-sizing problem with setup times and shortage costs[END_REF] and represents the different time windows during two consecutive iterations of the heuristic. We define sol k (x) as the value of the decision variable x in the solution obtained at iteration k of the algorithm and denote

a k = 1 + (k -1)(δ -σ) and b k = kδ -(k -1)σ
the first and the last period of the associated decision window, respectively. Let R k be the corresponding relaxed subproblem, which is solved using either the aggregate or the facility location formulation introduced in Section 3. Problem R k is similar to the original one, except for the variables definition constraints (9) that are replaced by the following ones :

w ij mt ∈ {0, 1} ∀i, ∀j, ∀m, ∀t ∈ [a k , b k ] (30) w ij mt ∈ [0, 1] ∀i, ∀j, ∀m, ∀t ∈ [b k + 1, T ] (31)
w ij mt = sol k-1 (w ij mt ) ∀i, ∀j, ∀m, ∀t ∈ [1, a k -1], ∀k ≥ 2 (32)
In order to control the total time necessary to obtain a solution with the RF procedure, we add two parameters that limit the computational time allocated to each subproblem R k . Specifically, we first set a target gap γ for the solution of R k , i.e. in each iteration k the solver returns the best solution found with a gap lower than or equal to γ. To avoid situations in which reaching γ induces a long computational time, we also set a time limit ρ for each iteration, i.e. the algorithm returns the best feasible solution obtained in ρ units of time, regardless of the gap achieved.

Note that the RF heuristic we describe in Algorithm 1 uses a forward period decomposition. A backward version, which starts at the end of the horizon and finishes at the first period, has also been tested in the literature [START_REF] Toso | Lot sizing and sequencing optimisation at an animal-feed plant[END_REF]). Since the two versions displayed equivalent performances during our preliminary tests, we focus on the forward version in our numerical experiments.

Algorithm 1 Relax-and-Fix heuristic

1: procedure RF(δ, σ, γ, ρ) 2: k ← 1, a k ← 1, b k ← δ, sol k ← ∅ 3: while b k < T do 4:
sol k ← Solve R k with maximum gap γ within ρ units of time 5:

a k+1 ← b k -σ + 1, b k+1 ← min{b k + δ -σ, T }, k ← k + 1 6: end while 7: Solve R k 8:
Return sol k 9: end procedure

Fix-and-Optimize

The RF approach is frequently combined with a local search procedure called Fix-and-Optimize to improve the feasible solution found. As for the RF heuristic, the procedure relies on a decomposition of the original problem into simpler subproblems, where the most classic decomposition is again based on periods. Starting from an initial reference solution, the idea is to "release" the subset of its binary variables whose time index belongs to a given decision window and reoptimize over them in order to improve the solution. If a better solution is found, it becomes the new reference and is used as such in the next iterations. During each iteration, binary variables that do not belong to the subset mentioned above are fixed either to their initial value or to a new value set in a reference solution from a previous iteration. At each iteration, only the binary variables (i.e setup variables) outside the decision window are fixed when others have to be reoptimized.

RFFO applied on the CLSSD-PM

The combination of RF and FO heuristics, denoted RFFO, is particularly efficient for complex lot-sizing problems and usually offers a good trade-off between execution time and solution quality. For instance, James and Almada-Lobo (2011) present these two approaches for the capacitated lot-sizing problem with sequence-dependent setups on parallel machines and obtain good results in short computational time. In both case, the size of the decision window greatly influences the computational time of the procedure. The different steps of the proposed RFFO procedure are decribed below :

Initialization Define the parameters δ, σ, γ, ρ for the RF and FO steps.

Step 1 Apply the Relax-and-Fix heuristic RF(δ, σ, γ, ρ) to obtain a feasible solution.

Step 2 Apply the Fix-and-Optimize heuristic on the current solution with the same parameters δ, σ, γ, ρ.

An analysis of the impact of the different parameters (i.e the decision window size, the overlap between two consecutive decision windows, the limiting gap and time during each iteration) is presented in Section 5. In the remainder of the paper, we extend the notation introduced in the RF algorithm and denote RFFO(δ, σ, γ, ρ) the RFFO procedure using the parameters δ, σ, γ, ρ in Step 1 and Step 2.

Clustering approach

As mentioned in the previous section, SD setups are computationally challenging when one wishes to obtain a good solution in a reasonable time, even when they rely on fast heuristics. In fact, the previous results obtained in the simpler framework of a singlemachine problem provide evidences that SD setups are much harder to tackle than the traditional item-dependent setup. In contrast with this setting, we call sequenceindependent (SI) a setup time that only depends on the next item to produce on a given machine, i.e. where λ ij m = λ j m is independent of i for all (i, j, m) ∈ N 2 × M. Our main idea is to find a way to transform an instance I of CLSSD-PM into a smaller instance Ĩ, for which some of the setups are considered sequence-independent. Since this operation may remove crucial information on the true setup times associated with a production sequence, we shall ensure that (i) there exists a feasible production sequence on each machine in each period that contains the same quantity of each item within the capacity available, and (ii) the SI setup times associated with a production plan for Ĩ remain comparable to their SD counterpart in I.

We build our conversion procedure on a simple observation on the setup times characteristics that is verified in many real-life manufacturing problems encountered in food-industries. Indeed, the majority of switching times between items observed in practice are either "long" (40 minutes to 2 hours) or "short" (1 to 10 minutes) and satisfy the triangle inequality. For instance, let's consider the case of a manufacturer that produces three types of cookies. Type 1 cookies contain allergen ingredients while types 2 and 3 are both free of allergens. The cleaning time between type 1 and any of the two others will likely be significantly longer than the one between types 2 and 3, and extend the setup times accordingly. Similar asymmetries may also occur from differences in oven temperatures, organic/non-organic product, etc.

In an attempt to reduce the global setup times and maximize resources utilization, practitioners often adopt an intuitive scheduling policy and try to regroup items with similar ingredients or recipe characteristics in the production sequences, with the objective to limit the number of long setup times in their production sequences. The main driver of the total setup time then becomes the arrangement of these group of items, rather than the sequence of items within a particular group. Our approach takes advantage of this specificity and use the data available on the setup times to build logical families of items. This allows us to define approximate instances for which we only consider sequence-dependent switching times between families, while intra-family item setup times are approximated by virtual, sequence-independent ones. In Robinson and Lawrence ( 2004), a problem including major family setups and minor setups are presented for the coordinated replenishment lot-sizing problem but the authors do not consider the sequence-dependent features. This approach reduces significantly the number of sequence-based decision variables, therefore lowering the size of the problem.

A clustering approach

Our procedure consists in partitioning the items into clusters that gather items together into families, such that the setup time between any two items in the same cluster is negligible compared to the setup time between two items that belong to different clusters. To this end, we rely on a partitioning algorithm called PAM (Partitioning Around Medoids) sometimes referred to as k-medoids, whose only inputs are an integer k and the setup times between all pairs of items. Hence PAM only requires a matrix of distances for each pair of points, a desirable feature compared to alternative procedures such as the well-known k-means that requires to define the coordinates of every point in a normed space. Note that the original PAM algorithm considers symmetric distances to determine which medoid is the closest to a given point. To meet this requirement while ensuring the robustness of our approach, we define a virtual setup time λij m = max{λ ij m , λ ji m } for all i, j ∈ N and m ∈ M. Thus two items with asymmetric changeover times (caused by organic components for instance), are likely to end up in two different clusters if the maximum value of their changeover times is large.

Although PAM is a well-known algorithm in the Machine Learning community, we quickly describe it below to make the paper self-contained. In a nutshell, the procedure builds clusters in an iterative fashion around specific points (items), considered as the medoid (center) of their respective cluster. At any given stage of the algorithm, the quality of the current partition on a given machine m ∈ M is evaluated using a cost function E m that is defined from the virtual setup times as follows:

E m = i∈I j∈Ci λij m (33)
, where I ⊆ N is the set of medoids in the current solution and for all i ∈ I, C i is the set of items that belong to the cluster associated with i. The procedure starts with an initial partition based on k randomly selected medoids as the reference and proceeds by performing successive swaps between one of the current medoid and other (nonmedoid) items. It then assigns every point to its closest medoid and takes this new configuration as the new reference if it improves the quality of the clusters obtained (based on the value of E m ). The procedure repeats the previous steps until no further improvement is observed. Algorithm 2 presents the PAM algorithm.

As mentioned earlier, one of the inputs required by PAM is the number k of clusters to define. It thus remains to estimate a value of k that ensures a relevant partition of items. To this end, we use an indicator based on the so-called silhouette score for Return I 19: end procedure every item, which measures its matching quality with the other items that belong to the same cluster. Formally, let m ∈ M and i ∈ I: The score obtained by an item j ∈ C i relies on a combination of the average distance between j and other items in C i , versus the average distance between j and the items in the cluster C i closest to j for i = i. Let us define a m (j) as the average setup time from item j to other items in cluster C i on machine m:

a m (j) = 1 |C i | -1 j ∈Ci,j =j
m (j) = min i =i 1 |C i | j ∈C i λjj m (35)
The silhouette score of item j, noted s m (j) is then obtained from (34) and ( 35):

s m (j) = b m (j) -a m (j) max{a m (j), b m (j)} (36)
Note that the closest to 1 the silhouette score is, the better the item is positioned in the cluster. One often considers having a convenient clustering for average silhouette higher than 0.75. To obtain more details on the meaning of the silhouette score in a clustering approach, see [START_REF] Laan | A new partitioning around medoids algorithm[END_REF] and [START_REF] Campello | A fuzzy extension of the silhouette width criterion for cluster analysis[END_REF].

In order to select the best number of clusters k * m for each machine m, we loop over the possible values of k and apply the PAM algorithm at each iteration. For a given k and a given machine m ∈ M, we consider the output of the procedure and compute the average silhouette criterion as the mean of the silhouette score over all items:

S m (k) = 1 |N | i∈N s m (i) (37)
The higher this average score is, the better is the clustering quality for the machine we consider. We select the clustering with the greatest average silhouette, which is considered as the best possible on the machine for the given instance. In Section 5, we show that S m (k) is a good a priori indicator to obtain the best trade-off between computational time and solution quality. For a quick description of the procedure, we present the pseudo-code of the clustering pre-processing approach in Algorithm 3. Note that for a CLSSD-PM instance, we run this algorithm |M| times since we consider machine-dependent setups.

Note that we use PAM as a pre-processing procedure before the definition of the MIP model. Its runtime is negligible (less than 15 seconds for the largest instances) compared to the one necessary to the solver to obtain a solution to the CLSSD-PM. Return bestClustering 12: end procedure

Algorithm 3 Clustering

Cluster-based mathematical formulation

We now define a virtual setup time v j m = max k∈Ci {λ kj m } for each item j ∈ C i , which only depends on the end point of the corresponding arc in the directed path representing the production sequence. For each pair of medoids (i, i ) and for each machine we also define a new sequence-dependent setup time f ii m = max j∈Ci,j ∈C i {λ jj m -v j m } incurred when an item of C i follows an item of C i in the production sequence. Note that this setup time excludes the end item individual setup time defined above, since the latter is added to the total setup time if the item is indeed included in the production sequence.

With these new virtual setup times, we are now able to only consider sequencedependent setup times when the production switches from one cluster to another. In particular we no longer take into account the production sub-sequences within a cluster and therefore replace the transition variables w ij mt with traditional setup binary variables y j mt . In our modified model, the switching decision between two clusters is represented by new binary variables z ii mt , equals to 1 if cluster i is setup just after cluster i on machine m in period t and 0 otherwise. On the other hand, whenever an item is produced on machine m in a given time-period, we only account for its item-specific setup time v i m in the total active time of the machine. Specifically, constraints (3) to ( 8) can be directly adapted to take into account the clusters instead of the items in the sequences of production. Constraints ( 12) are replaced by :

Ũmt = j∈N (τ j m x j mt + v j m y j mt ) + (i,i )∈I 2 f ii m z ii mt ≤ C mt + O mt ∀m = 1, . . . , M, ∀t = 1, . . . , T (38)
The constraint (25) for the facility location formulation is also replaced in an equivalent fashion. We add (39) to force a cluster setup if we want to produce an item in the cluster :

i∈I z ii mt ≥ y j mt ∀i ∈ I, ∀j ∈ C i , ∀m = 1, . . . , M, ∀t = 1, . . . , T (39)
We denote MIP-AGG-C and MIP-FL-C the modified versions of problems MIP-AGG and MIP-FL, respectively. It remains to prove that one can easily build a feasible solution to MIP-AGG (resp. MIP-FL) from a feasible solution to MIP-AGG-C (resp. MIP-FL-C) without any overcost. We prove this property for MIP-AGG, but a similar analysis can be conducted to obtain the same results for MIP-FL. Let (x, z, y, O) be a solution to MIP-AGG-C. For all j ∈ N , define µ(j) as the medoid of the cluster of item j, i.e. j ∈ C µ(j) , obtained after applying the PAM clustering algorithm. We propose the following simple procedure to create a solution to MIP-AGG: For all m ∈ M, t ∈ T , we define the following transition variables:

w 0j mt = 1 if j = min{j ∈ N : z 0µ(j ) mt = 1} 0 otherwise ( 40
)
w jj mt =      1 if j = min{j ∈ C µ(j)
, j > j : y j mt = 1} or z µ(j)µ(j ) mt = 1 and j = min C µ(j ) 0 otherwise (41)

The following proposition states that the above transformation yields a feasible solution to MIP-AGG: 40) and ( 41) is a feasible solution to MIP-AGG that incurs at most the total cost of s C .

Proof. Let m ∈ M and t ∈ T and compute the total active time on machine m in period t in solution s

U mt = j∈N τ j m x j mt + i∈N λ ij m w ij mt = j∈N τ j m x j mt + i∈N (λ ij m -v j mt + v j mt )w ij mt ≤ j∈N τ j m x j mt + v j mt i∈N w ij mt + i∈N max i ∈Cµ(i) {λ i j m -v j mt }w ij mt (42) = j∈N   τ j m x j mt + v j mt y j mt + i / ∈Cµ(j) max i ∈Cµ(i) {λ i j m -v j mt }w ij mt   (43) = j∈N τ j m x j mt + v j mt y j mt + j∈N i / ∈Cµ(j) max i ∈Cµ(i) {λ i j m -v j mt }z µ(i)µ(j) mt = j∈N (τ j m x j mt + v j m y j mt ) + (i,i )∈I 2 f ii m z ii mt = Ũmt ≤ C mt + O mt
Inequality (42) comes from the fact that max i ∈Cµ(j) {λ i j m -v j mt } = 0, while ( 43) follows from the definition (41) of w ij mt . To conclude the proof, notice that all variables with nonzero coefficients in the objective function are the same in s and s C , except U mt ≤ Ũmt whose coefficient is positive. Hence the total cost incurred by s is at most equal to the total cost incurred by s C .

In general, the partition of items into clusters reduces significantly the number of decision variables related to the sequence of production. In the cluster version of CLSSD-PM, the latter becomes quadratic in the number of clusters, as opposed to quadratic in the number of items in the original version.

Our definition of the virtual setup times ensures that any production sequence of items intra and inter-cluster will be lower than or equal to Ũmt , as defined in (38), hence it is always possible to construct a feasible solution with at most a production time Ũmt for all m ∈ M and t ∈ T . In particular, any algorithm that reconstructs intra-cluster production sequences (random, closest neighbor, etc.) from a feasible solution to MIP-AGG-C and then concatenates them according to the sequence of clusters corresponding to variables z ii mt is guaranteed to be feasible for MIP-AGG. In our computational experiments, we use a greedy algorithm to reconstruct the optimal production sequences. Finally, note that in addition to the flexibility that the method offers to the scheduler, one may also take advantage of such reductions in production time to decrease production costs or even to serve demands that are considered lost in MIP-AGG-C.

Computational experiments

The data that our industrial partner could gather from its clients was neither sufficiently complete nor precise enough to perform a thorough numerical study. We therefore chose to generate a new set of instances based on the information provided and use them to compare the different methods introduced above. Data sets are inspired by different real cases from food-industry producers and are built to cover a great variety of situations that one can expect to encounter in practice. One common characteristic of these instances lies in the structure of the setup times, which permits to apply the clustering approach described in the previous section. Although this assumption may seem restrictive, it reflects the distinctive feature of production processes in the food-industry.

The algorithms are coded in Java and we use IBM Ilog CPLEX 12.10.0 version as the MIP solver The experiments are performed using four threads on an Intel Xeon X5650 @ 2.57 GHz. For CPLEX and CPLEX-C, we set a computational time limit of 900 seconds. For each method tested, we evaluate the solution quality based on the deviation of its objective value from the lower bound LB obtained after 4 hours with CPLEX using the MIP-FL formulation, known to yield better lower bound. If U B corresponds to the best feasible solution found by the method tested, the gap is then defined as follows:

Gap = 100 • U B -LB LB (44)
We use the following abbreviations to identify the algorithms used in our test :

-CPLEX : branch The RFFO results presented in this section use the best settings obtained in section 5.2, i.e. the ones reaching the best trade-off between solution quality and computational time. After preliminary tests, we also observed that the solution methods that rely on the MIP-AGG formulation obtain better feasible solutions compared to their MIP-FL counterpart. Hence, we focus on the former and only use the latter to compute the best possible lower bound.

Instance generation

Problem size. Recall that our main goal is to obtain good solutions in a short computational time on large instances. The time limit is an important aspect for our industrial partner and is loosely set to a few minutes. We agreed afterwards to set it to 180 seconds, with a tolerance to up to 300 seconds for some instances (namely with more than 100 different products). In order to determine the influence of problem size (based on the number of items, machines and periods) on both the quality of the solution and the computational burden, we generated several instances for different combinations of values for N , M and T . We then partitioned the instances obtained into three classes (small, medium and large), presented in Table 3. In order to test the performances of the clustering approach in different settings, we use two different approaches that generate symmetric and asymmetric sequence dependent setup times, respectively. We create six instances for every combination of N , M and T , giving us a total of 336 instances. Definition of problem parameters. We derive our instances from the information provided by our industrial partner, with the objective to build data sets as close as possible to real-life cases. The parameters are defined as follows:

Time Capacity. The production time capacity on each machine, in each period is either 8, 16 or 24 hours. The maximum overtime is equal to 2% of the total capacity but is not considered in the latter case. Processing time for multiple machines. For instances with 2 machines or more, we perform for each item i an independent Bernouilli trial to decide whether it can be produced on several machines or not, with a success probability drawn uniformly in [0.25, 0.5]. In case of success, we draw uniformly between 2 and M the number M i of machines compatible with i and randomly select their indices as a subset of {1, . . . , M } of size M i . Whenever a machine can produce a given item, the per-unit processing time is drawn uniformly in [0.05, 0.5] minutes. Demand. We first determine how the total cumulative demand is allocated across the different items. Specifically for each item i, we draw a random number from a uniform distribution U (0, 1) and normalize the values obtained so that they sum up to 1. We then define the demands in each period t in an iterative fashion, assigning each unit to a specific item according to the probabilities previously defined. The procedure adds units of demand one after another as long as i min m {τ im }d i t /Util < m C mt , where Util ∈ {0.6, 0.8} represents the capacity utilization of the machines (decided beforehand as an input to the instance generator). Finally, we consider that 25% of the items present a periodic demand pattern, in which periods with null demand are regularly spaced in the planning horizon. Setup times. As mentioned earlier, we define two types of setup times, expressed in minutes. In what follows, "short" setup times are drawn uniformly in the interval [1, 10], while "long" ones are drawn uniformly in the interval [40,120].

To ensure they satisfy the triangle inequality, we initialize their values to zero and instantiate them sequentially: In each step of the procedure, we draw

λ ij m in the interval [max{l ij , |λ kj m -λ ik m |}, min{u ij , λ ik m + λ kj m }]
for any triplet i, j, k ∈ N , where l ij and u ij represent the minimum and the maximum of the corresponding (i.e. short or long) domain. We distinguish two classes of instances depending on how the setups are generated:

-Instances in the first set are strongly related to real applications. In the food industry, an item is often characterized by several features (for instance organic, specific oven temperature or food allergens) that influence the setup times. In practice, it is then the change in features between an item and its successor in the production sequence that effectively defines the length of this time interval. Based on this observation, we define 4 features for each item. Each feature has either 2 (e.g. organic and non organic) or 3 (e.g. oven temperatures) possible states. We then draw the setup times associated to each features randomly in the "short" or the "long" interval, depending on the nature of the feature considered. For the same feature, the switching time between two states may be asymmetric: A typical example arise when we set up a machine for allergen-free items after a production that contains allergen, a longer operation than the one performed for the opposite transition. As the setup operations corresponding to each specification are usually performed in parallel, we define the changeover time between two items as the maximum duration among all the switching times necessary for their respective features states. -We build the instances of the second set with the following procedure: (1)

We (randomly) choose a number of subsets, (2) We randomly assign each item to a subset and (3) We assign short setup times between all pair of items within the same subset and long setup times when they belong to different ones, ensuring in the process that the triangle inequality is always satisfied. This procedure guarantees that we obtain the clustering structure that is encountered in the food industry and for which our clustering approach behaves well. Safety stock. The safety stock of an item i in period t corresponds to the sum of demands for i over future periods t + 1, . . . , t + ∆ it . The number ∆ it of future periods considered is often called safety cover and is chosen according to a discrete uniform distribution U (1, 3). It can evolve over time with a low probability, but always includes the cover of the previous period (i.e ∆ i,t+1 ≥ ∆ it -1 for all i, t). Minimum production quantity. This quantity is generated from a uniform distribution U (1, 3).

Table 4 describes the generation of cost parameters. We consider an overcost equal to 50% of the regular line use cost. The lost sales cost for a given demand is built from the greatest possible cumulated production cost and holding cost incurred to produce the corresponding units. This ensures that serving a demand is always profitable whenever it is possible.

l i t max s≤t,m,j {p i ms + (cmt + cms)(τ i m + γ ji m q i min ) + θ i t -1 u=s h i+ u }

RFFO parameters analysis

The performances of the RFFO procedure heavily depends on the values δ, σ, γ and ρ. Hence it is crucial to determine the combination of parameters that provides the best expected trade-off between the cost of the solutions obtained and the average computational time for each class of instances. After an empirical preliminary analysis, we decided to focus on a limited number of settings. We set the number of overlapping periods σ to 1 since increasing its value does not improve significantly the solutions but severely deteriorates the computational time for medium and large instances. Even for small instances, we did not observe a substantial gain in performances for σ > 1. We tested different combinations of decision window length δ ∈ {2, 3}, gap tolerance (given by CPLEX) γ ∈ {1, 3, 5, 8, 10}% and time limit per iteration ρ ∈ {30, 60, 90, 120} seconds.

In order to decide which stopping criteria γ and ρ perform best, we analyzed their impact on the solutions. Note that this analysis was made on RFFO-C as RFFO often failed to find a feasible solution during the RF part on many instances. The results for medium size instances are shown in Appendix A. As expected, the gap of the final solution increases with γ, but this influence is mild when γ ≤ 5% (Figure A1) On the other hand the influence of γ on the computational time is more pronounced, for instance setting γ = 3% reduces significantly the computational time compared to γ = 1%. The former value thus provides the best compromise between the final gap and computational time. A similar analysis is conducted for different values of ρ. On the other hand, we observe in Figure A3 that the time limit is never exceeded for instances with N = 50. For other instances, the results suggest that ρ has a low impact on both the total computational time and the final gap. Hence, according to the results, ρ = 60 seems to be consistent as the time limit is seldom the limiting parameter when γ is adjusted appropriately.

The tables in Appendix B presents the results using RFFO on each class of instances. For each RFFO settings presented, we let σ = 1 and ρ = 60. We observe that δ = 3 leads to a better gap on small and medium instances for all the γ values we tested.

-For N ≤ 50, γ = 1% appears to produce the best solutions under 180 seconds in general. Hence we decided to apply RFFO-C(3, 1, 1, 60) on instances with N ≤ 50. -For larger instances, γ = 3% turns out to be a good choice to configure RFFO-C, except when N = 100 and T = 30, where γ = 5% gives a better trade-off between performances and computational time. -The effect of a bigger decision window is more visible on the largest instances (N = 125), with a significant deterioration of the gaps obtained on instances with (N = 125, M = 6, T = 30). In fact when N = 125 and T = 15, setting δ = 2 appears to be a more robust choice in general. We thus decided to retain RFFO-C(2, 1, 5, 60) for (N = 125, T = 15) and RFFO-C(3, 1, 8, 60) for (N = 125, T = 30) as the preferred configuration for the largest instances.

A summary of the RFFO configurations for the different types of instance is presented in Table 5.

Clustering approach results

To investigate the benefits of our clustering approach, we performed computational experiments using CPLEX directly and applying the RFFO procedure on both MIP-AGG and MIP-AGG-C. The results are presented in Table 6 and7. As expected, the direct use of a commercial solver on the MIP formulation turns out to be inadequate, even for some instances with 30 items. We allow CPLEX a maximum computational time of 900 seconds, which already exceeds the industrial requirements of 3 to 5 min-utes. While the solver is able to derive an optimal solution on some of the instances, the quality of the production plans obtained with this approach is not satisfying on average. Based on our observation, each acceptable solution from CPLEX was found towards the end of that time limit, which increases the variability of the solutions and makes it incompatible with an industrial approach.

In contrast, the clustering pre-processing reduces the size of the problem and allows to find significantly better solutions in most scenarios, with a faster convergence rate. The idea is that the loss in solution quality due to the approximation is compensated by a faster resolution. The results obtained suggest that it is the case for instances with N ≥ 40, for which CPLEX-C reaches better solutions than its sequence-dependent counterpart. We observe that this tendency becomes more and more pronounced as N increases. Finally, it appears that the length of the planning horizon T affects significantly the solution quality for both approaches.

We also test whether the silhouette score introduced in section 4.2 is a relevant criterion to drive the number of families considered. That is, we evaluate the solutions obtained by MIP-AGG-C for several number of clusters k and compare their quality with the corresponding average silhouette criterion S m (k) as defined in (37). The general conclusion is that the cost of the solution obtained decreases as this number grows to its "best" value k * = arg max{S m (k)}. When the number of clusters grows beyond k * , the solutions obtained either deteriorate significantly or show mild improvement. Indeed as the number of families increases, we introduce more sequence-dependent setup variables in the formulation MIP-AGG-C, which makes the problem closer to the original one but harder to solve. In our case, the silhouette criterion thus appears to be a suitable and easy-to-compute indicator of when this trade-off is favourable. Figure 3 illustrates this tendency on the worst case we observed across the instances generated for M = 2 and N = 100 for different number of clusters. The red and green curves show the average silhouette score S m (k) obtained on machine m ∈ {1, 2} for different number of clusters k. The blue curve provides an indicator of the solution quality, expressed as the ratio with the worst solution obtained across all the possible number of families considered. The quality of the solution clearly improves as the number of clusters increases towards k * . In the case of the RFFO resolution, we see that slightly better solutions are found for a larger number of clusters, but the increase in computational time is significant. On the contrary for a straightforward CPLEX resolution in 900 seconds, we observe that beyond a certain point, the solution quality deteriorates as we consider more clusters. We conclude that this criterion is a viable candidate to choose a good number of clusters a priori in order to obtain a good solution.

The experiments for the RFFO procedures are performed with the values δ, σ, γ and ρ adjusted for each class of instances according to the analysis in the previous section. In Table 6, we clearly observe a reduction of the average gap compared to a straightforward MIP resolution. On small instances, the gap obtained by RFFO combined with the clustering approach remains low compared to the others methods. In general, we observe that combining these two approaches significantly reduces the computational time while providing acceptable to good solutions on all the instances.

Despite performing better than CPLEX, the RFFO procedure directly applied to the CLSSD-PM problem generally leads to insufficient results or time overrun even on instances with N = 40. For larger instances, it often comes from incomplete iterations during the RF part of the heuristic, which cannot find a solution within the allocated time. A solution could be to increase the time limit ρ allowed for each iteration but this comes at the expense of a longer total computational time, which then frequently exceeds the acceptable threshold set by the practitioners. In contrast both the computational time and the average gap obtained remain acceptable when we combine the RFFO procedure with the clustering approach. In Table 7 we compare the results obtained with CPLEX and RFFO on the modified problems obtained after the clustering pre-processing is performed. The clustering pre-processing defines a number of clusters varying between 5 and 20. The results on medium size instances, shown in Table 7, indicate that RFFO-C reaches better average gaps while maintaining the computational time well below 3 minutes on average. However the minimum gap (Gap min ) and maximum gap (Gap max ) obtained show that there is a great variation in the solution qualities even within a single class of instances. Regardless, RFFO-C has the best worst case performance among the two methods, with a maximum gap of 204.19% compared to the one of CPLEX-C which exceed 7000%. In addition, we observed that the average computational time for RFFO-C is sensitive to the different parameters and size of the problem. Finally, for the 125-6-30 problem size instances, we exceed the computational time objective, which calls for additional refinements of the methods to obtain a quicker convergence. times are sequence-dependent. Since the number of setup variables strongly influences the computational effort necessary to solve such problems, we take advantage of the characteristics of those changeover times to apply a clustering algorithm and derive approximate instances of reduced size. The loss due to the approximation is regained thanks to the effectiveness of the method on large industrial cases. When we combine this approach with an adaptation of the classical Relax-and-Fix and Fix-and-Optimize heuristics, we are able to derive acceptable solutions in a very short computational time relative to the problem complexity. Our experiments allow us to test several combinations of parameters and determine which configuration is susceptible to provide the best results for different classes of instances. In addition, we provide evidences that the silhouette score is a good indicator to define an appropriate number of clusters with respect to the final solution. Computational experiments confirm that our procedure performs well on medium to large industrial instances, for which we were able to obtain feasible production plans in less than 3 and 5 minutes, respectively. Our approach remains useful on well structured instances in which setup times are not uniformly distributed. In the general case however, it would likely be irrelevant since we observe that a poor silhouette score results in a poor solution.

The approximate version of the problem obtained after our pre-processing clustering algorithm remains quite rich and includes many industrial constraints. In this study, we mainly focus on the contribution of the clustering on the resolution time and quality, but we could refine our conclusions by examining in details the different influencing factors. In particular, investigating the impact of the setup times or the parameter Util on the solutions obtained may help to better explain the variability observed on the gaps for similar problem sizes. From an academic standpoint, evaluating the impact of each type of constraint such as the inclusion of safety stocks or overtimes may help to identify which subproblems benefit the most from the clustering transformation. One could also simplify even further the problem by defining sequence-independent setup times between clusters of items. However, a direct application of our method may loose too much accuracy since inter-cluster changeover times contribute a lot more than their intra-cluster counterpart to the total machine usage time. The virtual setup times that are currently defined as the maximum over all the possible changeover times towards a given item are likely to be too far from the true value achieved. Hence, one should seek a different upper bound that retains more information on the production sequence to refine the approximation.

Finally, it is clear that all the solution methods presented in this paper are likely to take advantage of a strengthened MIP formulation. First, deriving a better lower bound would enable us to draw sharper conclusions on the quality of the solution obtained, especially on large instances. In addition, tightening the formulation is likely to speed up each iteration of the RFFO procedure while improving the gaps of the intermediate solutions. A natural research direction is therefore to explore the impact of adding valid inequalities on both the computational time and the final gap obtained. This approach may indirectly benefit from the clustering approximation developed in this paper since the structure of the modified problem may be compatible with a wider range of possible cuts compared to its sequence-dependent counterpart.

Figure 2 .

 2 Figure 2. Horizon decomposition in the RF heuristic

  Choose randomly k items ∈ N to be the first medoids and update I 4: Associate each item i ∈ N \ I to the closest medoids Compute E m for the set of medoids I = I ∪ {o} \ {i}

  b m (i) as the average setup time from the closest cluster of j on machine m: b

  1: procedure Clustering(maxClusters, N , m) 2:actualClustering ← ∅, bestClustering ← ∅, actualScore ← 0, bestScore ← 0, 3: for k = 2, . . . , maxClusters do 4: actualClustering ← Run PAM(k, N , m)

Proposition 4. 1 .

 1 Let s C = (x, z, y, O) a solution to MIP-AGG-C. Then s = (x, w, O) with w = w ij mt i,j∈N ,m∈M,t∈T defined from equations (

  -and-cut by CPLEX on MIP-AGG -CPLEX-C : branch-and-cut by CPLEX on MIP-AGG-C -RFFO : Relax-and-Fix & Fix-and-Optimize heuristic on MIP-AGG -RFFO-C : Relax-and-Fix & Fix-and-Optimize heuristic on MIP-AGG-C

Figure 3 .

 3 Figure 3. Effect of the number of cluster on the solution quality and the silhouette score

Table 1 .

 1 Sets and constants N : Number of items N = {1, . . . , N } M : Number of machines M = {1, . . . , M } T : Number of periods T = {1, . . . , T }

Table 2 .

 2 Models

		parameters
	d i t	Demand for item i in period t
	Cmt	Capacity (time) available on machine m in period t
	Cmt	Total capacity available on machine m in period t, including possible overtimes
	τ i m min q i S i t λ ik m h i+ t h i-t l i t p i mt	Time necessary to produce one unit of item i on machine m Minimum production quantity of item i Safety stock for item i in period t Changeover time from item i to item k on machine m Per-unit holding cost for excess inventory of item i in period t Per-unit penalty cost for item i in period t applied to inventory shortfall under S i t Per-unit lost-sales cost for item i in period t Per-unit production cost for item i on machine m in period t
	cmt	Cost of usage time for machine m in period t, per time unit
	cmt	Cost of overtime for machine m in period t, per time unit

Table 3 .

 3 Class

		instance size
	Class	N	M	T	Max. time (s)
	Small Medium	{20,30,40} {1, 2} {50}	{15, 30}	180
	Large	{75,100} {125}	{2,4} {4,6}		300

Table 4 .

 4 Cost

		parameters
	Parameter	Settings
	cmt	U (500, 1000)
	cmt	0.5cmt
	h i+ u	U (1, 3)
	h i-u	U (15, 30)
	p i mt	U (0.2, 0.5)

Table 5 .

 5 RFFO

			configuration
	N	20		30		40		50		75		100	125
	T	15	30	15	30	15	30	15	30	15	30	15	30	15	30
	γ (%)				1						3		5		8
	δ						3							2	3

Table 6 .

 6 Summary of the results obtained with each method on small instances

	Problem Size	Avg. Gap (%)	Sol. Time (s)
	N-M-T	CPLEX CPLEX-C RFFO RFFO-C	CPLEX CPLEX-C RFFO RFFO-C

Table 7 .

 7 Summary of the results obtained with CPLEX-C and RFFO-C on medium and large instances Gap min (%) Gap max (%) Avg. Gap (%) Sol. Time (s) Gap min (%) Gap max (%) Avg. Gap (%) Sol. Time (s)

	Problem Size CPLEX-C

Acknowledgement

The authors thank the two anonymous reviewers for their helpful comments. This work has been partially financed by VIF Software.

Data availability statement

The data that support the findings of this study are available from the corresponding author, Larroche, F, upon reasonable request.

Appendix A. Analysis of the stop criteria in RFFO A.1. Analysis of the gap limit γ


François Larroche

Odile Bellenguez

Guillaume Massonnet

Three-phase method for the capacitated lot-sizing problem with sequence dependent setups

Keywords: lot-sizing, iterative method, scheduling

This paper focuses on an industrial lot-sizing and scheduling problem that arises in the food industry and includes lost sales, overtimes and sequence-dependent setups on parallel machines. We propose a preliminary version of a three-phase iterative approach to optimize separately the affectation, the sequencing and the production of items. Our first numerical results suggest that with some additional improvements, this approach could be use in real-life by planners to reduce their costs.

Introduction

This paper presents a new heuristic to solve a lot-sizing and scheduling problem that arises in practical cases from the food industry. The problem we address is large and complex, which real-life applications often involve multiple production machines and specific constraints, which makes this problem complex to solve and motivates the planners to use optimization tools in order to reduce the costs incurred during their production process. However, modeling such complex systems leads to large mathematical formulations that are (in general) intractable with commercial solvers for Mixed Integer Programs (MIP). As a consequence, we aim to develop new approaches that enable us to obtain solutions to this problem in a quick and efficient fashion. Our work focuses on a multi-item capacitated lot-sizing problem including lost sales, safety stock, overtimes and sequence dependent setups on parallel machines. The discrete lot-sizing is a classical problem in the Operations Research literature since the work of Wagner and Whitin [START_REF] Wagner | Dynamic Version of the Economic Lot Size Model[END_REF] and it has since been enriched with various additional considerations in order to model more accurately practical situations. Limited production capacity is among the most popular extensions of the original problem and has been the topic of multiple surveys, see e.g. [START_REF] Karimi | The capacitated lot sizing problem: a review of models and algorithms[END_REF] and [START_REF] Quadt | Capacitated lot-sizing with extensions: a review[END_REF]. The model we consider further generalizes this setting to include lost sales and shortage costs, as in [START_REF] Absi | The multi-item capacitated lot-sizing problem with setup times and shortage costs[END_REF] where the authors introduce new classes of valid inequalities. Parallel resources production is more and more common in the literature, with a distinction between identical ( [START_REF] Beraldi | Rolling-horizon and fix-and-relax heuristics for the parallel machine lot-sizing and scheduling problem with sequencedependent set-up costs[END_REF]) and non-identical ( [START_REF] Kaczmarczyk | Modelling Set-up Times Overlapping Two Periods in the Proportional Lot-Sizing Problem with Identical Parallel Machines[END_REF]) machines. We focus on the latter in this paper.

The integration of lot-sizing and scheduling constraints greatly increases the complexity of the problems, which often leads researchers to propose heuristic solutions. Notable ones include [START_REF] James | Single and parallel machine capacitated lotsizing and scheduling: New iterative MIP-based neighborhood search heuristics[END_REF], who develop a MIP-based neighborhood search heuristics to address a lot-sizing and scheduling problem on parallel machines. [START_REF] Clark | Production setup-sequencing and lotsizing at an animal nutrition plant through atsp subtour elimination and patching[END_REF] present an industrial problem in which setup times depend on the sequence of production and propose a solution procedure based on subtour elimination and patching. [START_REF] Guimarães | Modeling lotsizing and scheduling problems with sequence dependent setups[END_REF] reviews the main modeling techniques for this class of problems and compare the efficiency of several solution methods. Iterative heuristics has been studied in the literature for the capacitatetd lot-sizing and scheduling problem in [START_REF] Mateus | Capacitated lot sizing and sequence dependent setup scheduling: an iterative approach for integration[END_REF]. In [START_REF] Dauzère-Péres | Integration of lotsizing and scheduling decisions in a job-shop[END_REF], the authors used a two-level method to deal with the lot-sizing and the scheduling phases separately. In the related field of Production Routing problem, [START_REF] Absi | A Two-Phase Iterative Heuristic Approach for the Production Routing Problem[END_REF] use a two-phase method to dissociate the production and the distribution part.

In this paper, we propose a three-phase iterative method for the production problem described above, that has been introduced in [START_REF] Larroche | Clustering-based solution approach for a capacitated lot-sizing problem on parallel machines with sequence-dependent setups[END_REF] and is referred to as CLSSD-PM for multi-item capacitated lot-sizing problem with lost sales, safety stock, overtimes, and sequence-dependent setups on parallel machines. Section 2 describes the notations and assumptions. In Section 3, we present an original three-phase iterative method, with a quick description of each phase. In Section 4, we succinctly present the preliminary results obtained on a set of instances generated from industrial data. Finally, we propose some ideas to improve the procedure as well as perspectives for future research in the Section 5.

Notations and assumptions

The problem studied in this paper is directly related to the problem originally presented in [START_REF] Larroche | Clustering-based solution approach for a capacitated lot-sizing problem on parallel machines with sequence-dependent setups[END_REF]. For an extensive presentation, the reader can refer to this paper. In this problem, we consider N different items that can be produced over a discrete finite horizon of T periods and M non-identical parallel machines. We denote N , M and T the set of items, machines and periods, respectively. For all items i ∈ N , we consider a deterministic demand d i t in each period t ∈ T , which is either satisfied from on-hand inventory or lost. We denote τ i m the production time of one unit of item i on machine m. Sequence-dependent setup times are denoted λ ij m for a given machine m and for each pair (i, j) of items. Note that we allow asymmetric setup times but assume they satisfy the triangle inequality. For each machine m ∈ M, we consider a time-dependent capacity C mt corresponding to the total time available in period t. However, in some circumstances this capacity can be exceeded up to a larger value C mt which is the "true" available time capacity in period t. In that case, C mt -C mt corresponds to the maximum overtime allowed. Finally, we also impose that every production of item i have to be greater than a minimum production quantity denoted q i min . The objective is to minimize the costs incurred in the production problem. We denote p i mt the cost of producing one unit of item i in period t on machine m. In addition, we denote c mt the per-unit of time usage cost of machine m during period t. For each unit of overtime, an additional cost cmt is applied. We define the safety stock as an "optimal" inventory level to reach in each period. Inventory-related cost are then defined with respect to its value, namely a deficit cost h i- t and an overstock cost h i+ t for each unit below or above the safety stock level. Finally, for each unmet unit of demand for item i in period t, a shortage cost l i t is incurred. Throughout the remainder of this paper, we use the following decision variables: In this section we present a three-phase iterative method to solve the problem under study. The central idea is to decompose the original problem into three (smaller) subproblems we can solve sequentially, where the output of a phase serves as an input for the following one. Below is an overview of one iteration of the heuristic:

x i mt ∈ R + Quantity of item i produced in period t on machine m y i mt ∈ {0, 1} binary variables equal to 1 if item i is affected to machine m in period t U mt ∈ R + Time usage of machine m in period t O mt ∈ R + Time in overtime of machine m in period t L i t ∈
-The first phase minimizes the costs induced by production and inventory management, considering virtual sequence-independent setup times per item. Capacity restrictions are also relaxed based on the information provided by the other phases. -The second phase proposes production sequences for the items affected to each machine and period. -The last phase decides the quantity to produce and pushes its output information into the first phase of the next iteration as an additional constraint.

The procedure loops through these different phases until a stopping criteria (either the time limit or a given number of iterations without improvement) is met.

First phase: Assignment. We describe the first phase with a mathematical model based on an aggregate formulation. In this phase, we do not consider the sequence dependent setup times, which are approximated. For that purpose, we introduce ST i mt the setup time for item i for machine m at period t (set to 0 at the first iteration). The aggregate formulation is expressed with the following MIP:

min t∈T m∈M (c mt U mt + cmt O mt ) + t∈T i∈N h i+ t I i+ t + h i- t I i- t + l i t L i t + m∈M p i mt x i mt (1) s.t. U mt = i∈N (τ i m x i mt + ST i mt y i mt ) ∀m ∈ M, ∀t ∈ T (2)
U mt ≤ C mt + O mt ∀m ∈ M, ∀t ∈ T (3)
I i t = S i t + I i+ t -I i- t ∀i ∈ N , ∀t ∈ T (4)
I i t = I i t-1 + m∈M x i mt + L i t -d i t ∀i ∈ N , ∀t ∈ T (5)
I i t ≤ D i tT (1 -z i t ) ∀i ∈ N , ∀t ∈ T (6)
L i t ≤ z i t d i t ∀i ∈ N , ∀t ∈ T (7)
x i mt ≤ D i t-1,T y i mt ∀i ∈ N , ∀m ∈ M, ∀t ∈ T (8)
x i mt ≥ q i min y i mt ∀i ∈ N , ∀m ∈ M, ∀t ∈ T (9)
L i t , I i t , I i+ t , I i- t ≥ 0 ∀i ∈ N , ∀t ∈ T (10)
U mt , O mt ≥ 0 ∀m ∈ M, ∀t ∈ T (11)
y i mt ∈ {0, 1} ∀i ∈ N , ∀m ∈ M, ∀t ∈ T (12)
z i t ∈ {0, 1} ∀i ∈ N , ∀t ∈ T ( 13
)
x i mt ≥ 0 ∀i ∈ N , ∀m ∈ M, ∀t ∈ T (14)
The objective (1) is to minimize the cost of production including line usage costs, lost sales and inventory stocks. Constraints (2) define the total working time of each machine which is equal to the production time and the approximated setup times. Constraints (3) define the time usage and overtime variables. Constraints (4) define the inventory stock from the safety stock, the overstock and the deficit. Constraints (5) are the inventory flow conservation equations through the planning time horizon. Constraints ( 6) and ( 7) ensure that lost sales for item i occur only when the corresponding stock is null. Constraints (8) use the cumulative demand over the remainder of the horizon as an upper bound on the quantity of each item produced on each machine in a given period. Constraints (9) force the production to be higher than its minimum value. Finally constraints [START_REF] Larroche | Clustering-based solution approach for a capacitated lot-sizing problem on parallel machines with sequence-dependent setups[END_REF] and ( 14) define the domain of specific variables. An additional constraint is necessary to link the information from the following phases after the first iteration. We denote y i mt and x i mt as the value of the corresponding variables y i mt and x i mt decided in the previous iteration. In addition, we also define for each machine m and period t the residual capacity Cap res mt as the remaining available machine time after the previous iteration of the algorithm. Before the first iteration, we initialize these values as follows : y i mt = 0, x i mt = 0 and Cap res mt = C mt . The following constraints (15) then indi-cate that for each item, the production time is bounded by the last solution but can be increased or decreased depending on items that are added to/removed from the assignment derived in the previous iteration and the residual capacity. Parameters σ i mt then correspond to the proportion of the production time available that is allocated to item i.

τ i m x i mt + (1 -y i mt )ST i mt y i mt ≤ τ i m x i mt + ( j =i x j mt =0
(1 -y j mt )(τ j m x j mt + ST j mt )

+ Cap res mt -

j =i x j mt =0 (τ j m q j min + ST j mt )y j mt )σ i mt ∀i ∈ N , ∀m ∈ M, ∀t ∈ T (15)
Second phase: Sequencing. From the first phase, we obtain the assignment of items to each pair machine-period (m, t). The goal of this phase is then to solve M × T sequencing problems with a straightforward mathematical formulation adapted from the TSP problem. For conciseness we do not detail it in this paper.

Third phase: Production. The last phase solves the production problem. We denote Λ mt the sum of setup times on machine m in period t obtained from the previous phase. If Λ mt + i∈N q i min y i mt > C mt , i.e. the capacity is not sufficient to carry out the setup times and the minimum production for the current affectation, we set Λ mt =0 and y i mt = 0 for all i ∈ N . Otherwise, the values Λ mt and ỹi mt remain the same. Note that since the assignment decisions ỹi mt are fixed in the first phase, they are considered as parameter in this step and thus the only binary variables considered are z i t to ensure a First-Come First-Serve discipline, which make the problem easy to solve. min t∈T m∈M

(c mt U mt + cmt O mt ) + t∈T i∈N h i+ t I i+ t + h i- t I i- t + l i t L i t + m∈M p i mt x i mt (16) s.t. U mt = i∈N τ i m x i mt + Λ mt ∀m ∈ M, ∀t ∈ T (17)
U mt ≤ C mt + O mt ∀m ∈ M, ∀t ∈ T (18)
O mt ≤ C mt -C mt ∀m ∈ M, ∀t ∈ T (19)
Constraints ( 4)-( 7)

x i mt ≤ D i t-1,T y i mt ∀i ∈ N , ∀m ∈ M, ∀t ∈ T (20) x i mt ≥ q i min y i mt ∀i ∈ N , ∀m ∈ M, ∀t ∈ T (21)
Constraints ( 10),( 11),( 13),( 14)

Update the parameters of the first phase. At the end of the third phase, there is a crucial update step on the input parameters of the first phase. The algorithm updates the individual setup times based on their marginal contribution to the production sequence computed in the second step. That is, for each item i in the sequence of production on machine m in period t, the setup time ST i mt corresponds to the additional time needed to insert this item is the sequence based on its direct predecessor and successor. On the other hand, we use the cheapest insertion time for all items that are not already included in the current sequence. This idea is described more in details in [START_REF] Absi | A Two-Phase Iterative Heuristic Approach for the Production Routing Problem[END_REF]. Current affectation decisions are stored in the variables ỹi mt , while a priori production quantities xi mt are directly extracted from the production decisions x i mt in the third phase. The three-phase algorithm is presented in flowchart 1.

Experimental results

In this section, we present the experimental results obtained with our threephase method. Our instances are designed from actual data from industrial cases. We derived 48 different instances combining the following parameters : N ∈ {20, 30, 40}, T ∈ {15, 30}, M ∈ {1, 2}. Note that these preliminary tests are not representative of the industrial reality, where the number of different items can be significantly larger than this benchmark. We use IBM Ilog CPLEX v12.10.0 to solve each phase and set a time limit of 900 seconds as a stopping criterion. If U B corresponds to the best feasible solution found by the method tested, the gap is defined as follows:

Gap = 100 • U B -LB LB (22)
, where LB is the best lower bound obtained after 4 hours with CPLEX. Table 1 presents a Gap comparison between the 3-phase method and a straightforward MIP resolution in 900 seconds by CPLEX. We set σ i mt = 1/N in constraints (15). Although the results suggest that our method needs refinement before we consider applying it on larger instances, it shows some promises as a decision support tools for practitioners. In particular for larger instances, the 3-phase method seems to be more robust to obtain adequate solution. We observe that the computational time is mainly decided by the first and the second phase while the third one is solved quickly. For instances with N ≥ 30, the second phase take a significant amount of the total computational time.

Perspectives

We study a problem encountered in the food industry. Since the problem is too complex to obtain good solution using a straightforward MIP resolution, we base our approach on an iterative resolution from a decomposition into three smaller subproblems. Our first preliminary results suggest that some enhancement must be implemented to improve the method. To that end, we could consider the addition of diversification mechanisms in the first phase to speed up the convergence.

Another direction under investigation rely on the clustering approach introduced in [START_REF] Larroche | Clustering-based solution approach for a capacitated lot-sizing problem on parallel machines with sequence-dependent setups[END_REF] to simplify the second phase and compute efficiently good production sequences. Another research direction is to use heuristics or MIP-heuristics in the first and second phases to speed up the resolution.

Table 1 .

 1 Average Gaps obtained with the 3-phase method

	Start
	Initialize all parameters :
	x i mt , y i mt ST i mt , Cap res mt