Ilhem Slama

Oussama Ben-Ammar
email: oussama.ben-ammar@mines-ales.fr

Simon Thevenin

Alexandre Dolgui

Faouzi Masmoudi

Stochastic program for disassembly lot-sizing under uncertain component refurbishing lead times

Keywords:

Stochastic refurbishing lead time Monte-Carlo sampling Rolling horizon strategy a b s t r a c t

Planning disassembly operations for a given demand in components is challenging in practice because the quality of recovered components is very uncertain, and thus the duration of refurbishing operations is unpredictable. In this paper, we address the capacitated disassembly lot-sizing problems under uncertain refurbishing durations. More precisely, we consider a two-level disassembly system with a single type of end-of-life product, a dynamic demand, and stochastic refurbishing lead times for all components. To deal with the static decision frameworks, this problem is modeled as a two-stage stochastic Mixed-Integer Linear Program (MILP), where the objective is to minimize the expected total cost. To alleviate the scalability issues, we propose a reformulation of the inventory constraint that significantly reduces the number of scenarios. In addition, to solve large scale problems, we couple this reformulation with Monte-Carlo sampling. We provide a rolling horizon approach to deal with the static decision framework, where disassembly decisions are updated when new information unfolds. Experimental results show the effectiveness of the proposed models and the convergence of the resulting Sample Average Approximation (SAA) estimator.

Introduction

Governments around the world are imposing more and more environmental regulations on manufacturers and suppliers to minimize waste and recover resources. The directive on waste electrical and electronic equipment (directive 2002/96/EC) is a good example of the implementation of an environmental management policy [START_REF] Wang | A two-stage robust programming approach to demand-driven disassembly planning for a closed-loop supply chain system[END_REF]. Members of the European Union were invited to set up different systems for returning products at the end of their life, reducing waste, and carrying out recovery and recycling processes. All of these actions aim to reduce electrical equipment wastes and their negative environmental impacts. According to [START_REF] Pini | Environmental performance of waste management in an italian region: How LCI modelling framework could influence the results[END_REF] , Europe produces up to 3 billion tons of waste every year. In this context, recovery processes provide financial opportunities, and they allow manufacturers to engage in sustainable practices by extending the life of an end-of-life (EoL) product. In fact, the refurbishing activity is not only attractive from an environmental perspective, it can provide economic benefits by obtaining valuable materials for reuse or energy production. Refurbished components, which typically undergo a thorough cleaning process are equivalent in performance to the new ones, yet they sell for only 40% of the original price, moreover the cycle of resale and refurbishing can be repeated several times [START_REF] John | The profit-making allure of product reconstruction-MIT SMR store[END_REF]. Currently, in the United States, refurbishing is $50 billion industry with direct employment of about 480,0 0 0 in 73,0 0 0 firms' [START_REF] Gutwski | Remanufacturing and energy savings[END_REF].

The recovery process encompasses all activities related to the collection, disassembly, refurbishing, repair, recycling, disposal, of EoL products. Once collected, they are disassembled into components, and these components are then cleaned and tested. Based on their condition, these components undergo one of the two value recovery activities (refurbishing, or recycling) or are disposed of [START_REF] Langella | Heuristics for demand-driven disassembly planning[END_REF] . As the reverse flow grows, many firms seek opportunities to reduce EoL product management costs. Among various tasks in the recovery process, disassembly and refurbishing are the most significant operations in terms of costs. In addition, these operations are difficult to manage due to the uncertainty of various refurbishing parameters such as quantity, quality, and timing.

According to a review proposed by [START_REF] Jaehn | Sustainable operations[END_REF], the EoL product recovery process that considers random input is a promising area of research for sustainable operations in production. After disassembly, the components undergo several renovation processes (such as repair, cleaning, etc.) to meet the required quality level, and this renovation process can take several weeks. In a remanufacturing context, these refurbishing lead times (RLT) are difficult to predict, because the quality of the components is only known once the EoL item is disassembled. Managing the uncertainty of refurbishing lead times plays an important role in disassembly planning and inventory management since bad management can cause stock-outs and unnecessary stocks.

In this paper, we provide a methodology for disassembly planning of EoL products under stochastic refurbishing time for components, known component demand, and when the capacity of the disassembly line is limited. The goal is to decide the size of the lots to disassemble and to plan the amount of overtime required on the disassembly line. According to the classification proposed by [START_REF] Ji | Capacitated disassembly scheduling with parts commonality and start-up cost and its industrial application[END_REF], the studied problem enters into the category of a Capacitated Disassembly Lot Sizing (CDLS) problem that considers the best use of disassembly resources to satisfy the demand on components over a given planning horizon. Note that we assume the disassembly duration is known, and it is shorter than a period, whereas the refurbishing duration may last several weeks and is unknown. In this setting, the disassembly duration corresponds to process duration, and it consumes the production capacity. In the production planning literature, the time difference between the start of processing of an operation and the availability of the resulting parts is called lead time. Therefore, the considered problem is a CDLS under uncertain refurbishing lead times.

As we aim to extend the current knowledge about DLS problem resolution by integrating the uncertainty of refurbishing lead times, we consider a single type of EoL product, a two-level disassembly system, and a multi-period planning. In addition, the customer's demand and the disassembly duration are known, whereas the refurbishing lead times are discrete random variables that follow known probability distributions. This tactical problem aims to determine the optimal quantity of EoL products to disassemble in each period over a given planning horizon which guarantees the best trade-off between costs related to disassembly operations and the sum of expected inventory and backlog costs. The contributions and innovations in this study are the following:

1. We formulate the considered problem as a stochastic mathematical program. The resulting MILP accounts for a set of possible scenarios that represent possible refurbishing lead times for components obtained by disassembly. The proposed model is flexible enough to account for the case where different disassembly lots cross in time due to the uncertain refurbishing lead times. In this situation, if the disassembly of a lot A occurs before the disassembly of a lot B , the components resulting from B may be received before the components of A . This is an important contribution since a large part of the literature on classical planning under lead time uncertainty requires the restrictive assumption that orders do not cross in time. To the best of our knowledge, we are the first to propose a method able to account for order crossover in the DLS problem. 2. Our experiments show that CPLEX cannot solve instances with more than 30 components and 10 periods in less than 3600 seconds of computation, because the number of scenarios becomes too large. To get around this obstacle, we developed a scenario aggregation technique. More precisely, we propose a new formulation of the inventory constraint that significantly reduces the number of scenarios. The resulting model significantly improves the scalability of the approach, and it can solve instances with a large number of components and peri-ods when order crossover is forbidden, or if the refurbishing lead time range (i.e., the difference between the maximum and the minimum refurbishing lead time value for a component) is lower than 6 periods when order crossovers can occur. 3. Even though the scenario aggregation (SA) approach reduces the number of scenarios drastically, this number remains exponential in the refurbishing lead times range. As a consequence, the resulting method cannot solve instances with a large range of refurbishing lead times. To solve these large-scale problems, we couple the scenario aggregation technique with Monte-Carlo (MC) sampling. 4. To deal with the dynamic decision frameworks, where the disassembly plan is updated when new information unfolds, we employ a rolling horizon strategy. 5. Extensive computational experiments assess the effectiveness of the proposed approaches. In particular, the resulting sampling method provides a good approximation, with an average percentage deviation from the optimal solutions of less than 0.5%. In addition, our tests show the convergence of the Sample Average Approximation (SAA) estimator.

This paper is divided into six sections. Section 2 provides a brief overview of the literature review. Section 3 gives the problem description and formulation. The proposed approaches are described in Section 4 . Section 5 gives numerical results. Finally, conclusions and some perspectives for future work are drawn in Section 6 .

Literature review

Since the seminal work of [START_REF] Gupta | Scheduling disassembly[END_REF] , the DLS problems have attracted a lot of attention, and recent reviews are presented in [START_REF] Kim | Disassembly scheduling: Literature review and future research directions[END_REF] and Slama, Ben-Ammar, Masmoudi, & Dolgui (2019a) . [START_REF] Gupta | Disassembly modeling for assembly, maintenance, reuse and recycling[END_REF] highlighted the importance of disassembly lot sizing problems. Compared to the assembly process, the disassembly process is characterized by more complex operational and physical properties [START_REF] Gungor | Issues in environmentally conscious manufacturing and product recovery: A survey[END_REF]Slama, Ben-Ammar, Dolgui, & Masmoudi, 2020a ). The most important difference is the number of demand sources. In an assembly system, the parts tend to converge to a single demand source related to the final product. Under a disassembly system, parts tend to diverge to multiple demand sources corresponding to each disassembled component. For more details on these divergences, the interested reader is referred to [START_REF] Gungor | Issues in environmentally conscious manufacturing and product recovery: A survey[END_REF] and [START_REF] Kim | Disassembly scheduling: Literature review and future research directions[END_REF] .

Different variants of the DLS problem exist, and they differ according to the product structure (single or multi EoL products, two or multi-echelon bill of materials (BOM), with or without parts commonalities) and by the consideration of disassembly process resources (capacitated or uncapacited). Most of the literature on disassembly systems concerns the deterministic uncapacitated problem (e.g., [START_REF] Godichaud | Economic order quantity for multistage disassembly systems[END_REF]2019;[START_REF] Gupta | Disassembly modeling for assembly, maintenance, reuse and recycling[END_REF][START_REF] Kang | Disassembly leveling and lot sizing for multiple product types: A basic model and its extension[END_REF]Kim, Doh, & Lee, 2018a;[START_REF] Kim | Two-phase heuristic for disassembly scheduling with multiple product types and parts commonality[END_REF]Kim, Lee, Xirouchakis, & Kwon, 20 09;Langella, 20 07;[START_REF] Lee | A two-stage heuristic for disassembly scheduling with assembly product structure[END_REF][START_REF] Neuendorf | Disassembly scheduling with parts commonality using petri nets with timestamps[END_REF]Pour-Massahian-Tafti, Godichaud, & Amodeo, 2019a;2019b;[START_REF] Pour-Massahian-Tafti | Single product disassembly lot sizing problem with disposal[END_REF]2020a;2020b;[START_REF] Prakash | Constraint-based simulated annealing (CBSA) approach to solve the disassembly scheduling problem[END_REF][START_REF] Taleb | Disassembly of complex product structures with parts and materials commonality[END_REF], while a small number of studies focus on the deterministic capacitated version (e.g., [START_REF] Ji | Capacitated disassembly scheduling with parts commonality and start-up cost and its industrial application[END_REF]Kim, Lee, & Xirouchakis, 2006a;[START_REF] Kim | Capacitated disassembly scheduling: Minimizing the number of products disassembled[END_REF][START_REF] Lee | Disassembly scheduling with capacity constraints[END_REF]Slama et al., 2020a;[START_REF] Ullerich | Flexible disassembly planning considering product conditions[END_REF]. Table 1 summarizes the literature on the DLS problem, and indicates the resolution approach, the BOM structure, whether capacity is considered, and the type of uncertainties.

In practice, the product recovery process is sensitive to various sources of uncertainty, such as uncertain demand from cus- tomers, recovery rates, disassembly lead times, etc. [START_REF] Guide | An evaluation of order release strategies in a remanufacturing environment[END_REF] . As in supply planning for assembly systems, these uncertainties create disruptions in the disassembly plan, and leads to unmet customer demand. Many researchers have applied stochastic algorithms to the DLS problem (e.g., [START_REF] Tarim | An efficient computational method for a stochastic dynamic lot-sizing problem under service-level constraints[END_REF][START_REF] Tarim | Constraint programming for computing non-stationary (R , S) inventory policies[END_REF][START_REF] Vargas | An optimal solution for the stochastic version of the Wagner-Whitin dynamic lot-size model[END_REF], and they focus on demand, disassembly yield, and disassembly lead time (DLT) uncertainty [START_REF] Suzanne | Towards circular economy in production planning: Challenges and opportunities[END_REF]. Therefore, to position our research in the existing literature, we successively review the workson stochastic DLS problems under demand, yield, and disassembly lead time uncertainty. Demand uncertainty: Barba-Gutiérrez & Adenso-Díaz (2009) propose a fuzzy logic approach for the single EoL item, multi-period, multi-echelon DLS problem under stochastic demand. In this work, every unsatisfied demand in each period is back-ordered and postponed for future periods. [START_REF] Kim | Capacitated disassembly scheduling with random demand[END_REF] consider the capacitated multi-item DLS problem where each product has a two-echelon BOM structure. To solve the problem, the authors develop a stochastic inventory model, and they propose a Lagrangian relaxation-based heuristic that optimizes the sum of expected inventory holding, setup, and penalty costs. [START_REF] Fang | A stochastic production planning problem in hybrid manufacturing and remanufacturing systems with resource capacity planning[END_REF][START_REF] Quezada | A multi-stage stochastic integer programming approach for a multi-echelon lot-sizing problem with returns and lost sales[END_REF] ; [START_REF] Wang | A two-stage robust programming approach to demand-driven disassembly planning for a closed-loop supply chain system[END_REF] consider the multi-level, multi-product, multi-period stochastic lot-sizing problem in remanufacturing systems. [START_REF] Wang | A two-stage robust programming approach to demand-driven disassembly planning for a closed-loop supply chain system[END_REF] propose a two-stage stochastic programming model to find a compromise between the solution robustness and the expected cost. [START_REF] Fang | A stochastic production planning problem in hybrid manufacturing and remanufacturing systems with resource capacity planning[END_REF] formulate a multi-stage stochastic programming approach, and they propose a Lagrangian relaxation algorithm to solve the problem. Recently, [START_REF] Quezada | A multi-stage stochastic integer programming approach for a multi-echelon lot-sizing problem with returns and lost sales[END_REF] focused on a multi-echelon system with disassembly, refurbishing, and reassembly operations. [START_REF] Quezada | A multi-stage stochastic integer programming approach for a multi-echelon lot-sizing problem with returns and lost sales[END_REF] extend the multi-stage stochastic programming approach proposed by [START_REF] Fang | A stochastic production planning problem in hybrid manufacturing and remanufacturing systems with resource capacity planning[END_REF] by considering the uncertainties related to the quantity and quality of returned products, the customer demand, and the production costs. To solve the problem, the authors develop an exact solution approach based on a branch-and-cut algorithm.

The yield uncertainty occurs when the number of components obtained by disassembling a unit of product is uncertain [START_REF] Inderfurth | How yield process misspecification affects the solution of disassemble-to-order problems[END_REF]. To the best of our knowledge, only two papers consider the DLS problem under yield uncertainty. These two papers are restricted to the multi-products case with two-echelon BOMs in a one-period disassembly to-order (DTO) environment. [START_REF] Inderfurth | Heuristics for solving disassemble-to-order problems with stochastic yields[END_REF] give a heuristic to reduce the expected disposal, purchasing, and disassembly operation costs, whereas [START_REF] Inderfurth | How yield process misspecification affects the solution of disassemble-to-order problems[END_REF] propose a mathematical model for the specific case with three end-components and two sub-assemblies. The authors study the impact of process specification on the solution of the DTO problem. Finally, a statistical analysis shows that the yield uncertainty has a large impact in engine remanufacturing.

Uncertainty of demand and yield: to the best of our knowledge, [START_REF] Liu | Capacitated disassembly scheduling under stochastic yield and demand[END_REF] is the only publication investigating simultaneously the uncertainty of component demand and disassembly yield. The yield uncertainty is modeled with a uniform distribution, whereas the demand follows a normal distribution. The authors provide a mixed integer nonlinear programming (MINLP) for the multi-period and two echelon disassembly system, and they solve it with an outer approximation-based solution algorithm.

Uncertainty of disassembly lead time (DLT) : while yield and demand uncertainties impact the disassembly operation planning, remanufacturing facilities also face uncertain disassembly lead times. In fact, Guide (20 0 0) indicates that the average time to re-manufacture and disassemble a recycled EoL product varies from a low level of 5.54 hours to a high level of 300 hours, with coefficients of variation as high as 5.0. The literature on disassembly planning plan under random DLT is scarce. [START_REF] Slama | Scenario-based stochastic linear programming model for multi-period disassembly lot-sizing problems under random lead time[END_REF] propose a stochastic model for the multi-period, single product type, and two-level disassembly system. Recently, Slama, Ben-Ammar, Dolgui, & Masmoudi (2020b) proposed a generalization of the discrete Newsboy formulae to find the optimal release date when the time of disassembling the EoL product is a random variable. In this paper, the authors investigated the effect of uncertainty of DLT on the disassembly planning problem. Their analyses show that disassembly lead time variability has a significant effect on the system performance. The limitation of this work is that it considers a mono-period disassembly planning problem with unlimited capacity. Slama, Ben-Ammar, Dolgui, & Masmoudi (2020c) extends [START_REF] Slama | Scenario-based stochastic linear programming model for multi-period disassembly lot-sizing problems under random lead time[END_REF] by considering the capacitated case. Recently, Slama, Ben-Ammar, Dolgui, & Masmoudi (2021b) proposed a genetic algorithm to solve the multi-period, single EoL product, and two-level BOM. However, the authors assume the disassembly lead time is uncertain, where all the components of an end-item are received in the same period. On the contrary, the present work considers refurbishing lead time uncertainty, where the uncertainty affects each component independently.

As discussed in [START_REF] Kim | Disassembly scheduling: Literature review and future research directions[END_REF], the disassembly process is very sensitive to the uncertainty of refurbishing lead time of each component. Guide (20 0 0) states that this uncertainty can interrupt the delivery of disassembled parts, and it increases the inventory and backlog costs. Nevertheless, to the best of our knowledge, the uncertainty of the refurbishing lead time (i.e., the time required for component recovery) has not yet been investigated. In fact, after each disassembly operation, the part must be refurbished before being shipped to the customer. As the refurbishing lead time (RLT) depends on the quality of the disassembled items, it cannot be estimated accurately.

To fill this gap in the literature, the present paper investigates the DLS problem under component RLT uncertainty for the multiperiod planning, single EoL item, with a two-echelon disassembly BOM. Here, the demand is considered dynamic and we consider order crossovers. As far as we know, order crossovers have not been investigated in the DLS problem.

In the field of assembly systems, several studies have focused on lead time uncertainty with order crossovers. These works can be split into three groups according to the number of levels in the BOM, the number of periods in the planning horizon, and the aspect of the demand: (i) one-level assembly systems, multi-period planning with a one demand (e.g., [START_REF] Hnaien | Robust single-item lot-sizing problems with discrete-scenario lead time[END_REF][START_REF] Jansen | Setting optimal planned leadtimes in configure-to-order assembly systems[END_REF]2018 ), (ii) one-level assembly systems for a multi-period planning and a constant demand (e.g., [START_REF] Ould-Louly | Calculating safety stocks for assembly systems with random component procurement lead times: A branch and bound algorithm[END_REF]2013;[START_REF] Shojaie | A study on MRP with using leads time, order quality and service level over a single inventory[END_REF], and (iii) multi-level assembly systems for a mono-period planning and a one-demand (e.g., [START_REF] Axsäter | Planning order releases for an assembly system with random operation times[END_REF][START_REF] Ben-Ammar | A hybrid genetic algorithm for a multilevel assembly replenishment planning problem with stochastic lead times[END_REF]Ben-Ammar & Dolgui, 2018;[START_REF] Sakiani | Multi-objective supply planning for two-level assembly systems with stochastic lead times[END_REF]. However, most of these studies are based on analytic models that require problem-specific assumptions, and they do not apply to our case. Interested readers can refer to the literature review proposed by Díaz-Madroñero, Mula, & Peidro (2014) ; [START_REF] Golmohammadi | Review of supplier diversification and pricing strategies under random supply and demand[END_REF] .

Finally, as one of the main contributions of the present work is a reformulation of the stochastic program that leads to scenario aggregation, we mention that scenario aggregation techniques exist in stochastic optimization (e.g., Dupa čová, [START_REF] Dupa Čová | Scenario reduction in stochastic programming[END_REF][START_REF] Henrion | Scenario reduction in stochastic programming with respect to discrepancy distances[END_REF][START_REF] Hewitt | Decision-based scenario clustering for decision-making under uncertainty[END_REF][START_REF] Kuhn | Aggregation and discretization in multistage stochastic programming[END_REF]. These approaches seek a good representation of the uncertain parameters with a restricted set of scenarios. More precisely, the problem of optimal scenario reduction determines a scenario subset of prescribed cardinality (and a probability measure based on this set) that is the closest to the initial distribution. For instance, [START_REF] Hewitt | Decision-based scenario clustering for decision-making under uncertainty[END_REF] rely on clustering approaches to group scenarios that have a similar impact on the decision making. The resulting approach leads to highquality upper and lower bounds for stochastic models. [START_REF] Henrion | Scenario reduction in stochastic programming with respect to discrepancy distances[END_REF] rely on a discrepancy measure to build the scenarios set. The idea is to build a probability distribution represented by a few scenarios and that deviates at least as possible from the original distribution. [START_REF] Kuhn | Aggregation and discretization in multistage stochastic programming[END_REF] aggregate the scenarios based on upper and lower bounds of the optimal value of the original problem. Dupa čová et al. ( 2003) have proposed a scenario reduction approach for a convex stochastic programming problem with a discrete initial probability distribution. However, the existing methods are different from the one proposed in the present work. Here, we rely on the properties of the considered problem to reformulate the inventory constraint and aggregate the scenarios.

Problem description and formulation

Problem description

In the present work, we study the capacitated disassembly lotsizing problem in a two-level bill of materials. Fig. 1 shows the disassembly and refurbishing process considered in this work. The EoL products (first level) are disassembled into | N | components (second level), where the number in parentheses represents the disassembly yield (i.e., the number of components obtained by the disassembly operation). Such a BOM structure is common in practice, for instance, in the disassembly of an EoL vehicle ( Pour-Massahian-Tafti, [START_REF] Pour-Massahian-Tafti | New models and efficient methods for single-product disassembly lot-sizing problem with surplus inventory decisions[END_REF]. Once the factory has disassembled the End of Life item, each lot of component is sent to a different refurbishing facility, where they are repaired, cleaned, and refurbished. The refurbishing and repair duration varies significantly from a lot of component to another because it depends on the quality of the components. Each component undergoes different refurbishing operations depending on its type and its quality.

In this context, planning the disassembly operation is complex, and it requires to account for the uncertain refurbish and repair time. More precisely, once disassembled, the components i are available after a stochastic refurbishing and repair lead time L i . That is, if any disassembly operation is performed in period t, the components i are available at period t + L i . For sake of simplicity, we use the term refurbishing lead times to refer to refurbishing and repair lead times in the rest of the paper. As the refurbishing and repair times are specific for each component, the refurbishing lead times are non-identically distributed. In addition, we consider that the supports of the refurbishing lead time distributions are discrete and finite. More precisely, the refurbishing lead times L i are random discrete variables characterized by known probability distributions and varying between L - i and

L + i , ∀ i ∈ N , L - i , L + i ∈ N * ,
Here, N * is the set of non-negative integers except zero. Such an assumption is consistent with the literature on supply planning under uncertain lead times (see Ben-Ammar, Dolgui, & Wu, 2018 , and related works). Note that the disassembly orders can cross in time; i.e., if a lot A is disassembled before lot B , the components resulting from the disassembly of B can arrive before the components resulting from A .

The problem is to decide the disassembly lot sizes Z t in each period t, and the amount of extra capacity O t required in the disassembly facility. In most remanufacturing facilities, the disassembly operation is done manually, and it is a bottleneck. Consequently, we consider that the disassembled quantity is limited by a given capacity in each period C t . The capacity consumption corresponds to the disassembly time F which is known. As the disassembly line is the bottleneck, the disassembly lot-sizing plan must use its capacity correctly. Extra capacity can be added if a cost is paid. A classic solution is to use overtime, but it is also possible to bring in an extra worker for instance by using temporary workers. As in [START_REF] Voß | Introduction to computational optimization models for production planning in a supply chain[END_REF] , we assume 'overtime' which encompasses all these capacity additions with a cost penalty u t . Note that we consider that the disassembly time is known whereas the repair/refurbishing times are uncertain. This assumption is realistic since the refurbishing/repair time depends on the quality of the item. We assume that the time required to disassemble an item is known because the items are disassembled by batch, and the capacity consumption of each individual item averages out over the lot. Our objective is to manage the capacity on the disassembly line, whereas the refurbishing operations are performed in other parts of the shop floor.

The demand D it for each disassembled component i is known in each period t of the planning horizon. The stochastic refurbishing lead times can lead to holding costs h i and backlog costs b i whenever the demand of a disassembled component i is not satisfied on time. The objective is to optimize the disassembly quantities at level 1 to minimize the sum of the expected holding and backlog costs of all components in level 2, as well as the setup and overtime costs over a finite horizon. Also, the release of a disassembly order in period t leads to a setup cost s t in that period.

In this work, we consider the situation where the disassembly plan is decided over a planning horizon of several periods, and it is frozen. This decision process is required because the capacity planning decision (such as the schedule of the employees) cannot be modified in one period for the next one. Frozen horizons are commonly used in the industry to avoid nervousness. The academic literature (e.g., [START_REF] Sridharan | Freezing the master production schedule under demand uncertainty[END_REF][START_REF] Sridharan | Freezing the master production schedule under rolling planning horizons[END_REF][START_REF] Sridharan | Freezing the master production schedule: Implications for fill rate[END_REF][START_REF] Xie | Freezing the master production schedule under single resource constraint and demand uncertainty[END_REF][START_REF] Zhao | Freezing the master production schedule for material requirements planning systems under demand uncertainty[END_REF] shows that freezing the production planning has a notable impact on total cost, the service level of the system, and schedule instability. In addition, freezing the decisions is necessary when they must be made several months in advance. Such decisions include capacity adjustment since workers cannot be hired/fired at the last moment.

As the decisions are frozen, even though the refurbishing lead times are not revealed at once, the problem corresponds to a two-stage decision model. Nevertheless, the approaches proposed in this work can also be used in a dynamic decision framework, where the decisions are updated in each period after observing the quantities refurbished in each period. Solving the dynamic decision framework exactly requires a multi-stage model that would be too complex to solve. In fact, with uncertain refurbishing lead times, the setup decision should be adjustable to make good decisions. With the classical scenario tree representation, such an approach would lead to an exponential number of binary variables, and the resulting model would not be solvable for realistically sized instances. In this context, the two-stage model provides a heuristic that can be used in a rolling horizon framework and can make the problem more tractable from the computational point of view. Such an approach leads to close to optimal solutions for lot-sizing under uncertain demand [START_REF] Thevenin | Material requirements planning under demand uncertainty using stochastic optimization[END_REF].

Before moving on to the core topic of this paper, we briefly investigate the complexity of the problem. The deterministic single item capacitated lot-sizing problem (SCLP) is a special case of CDLS with a single disassembled item i and L i equals to 0 with probability 1. As the SCLP is NP-hard [START_REF] Bitran | Computational complexity of the capacitated lot size problem[END_REF][START_REF] Florian | Deterministic production planning: Algorithms and complexity[END_REF], the CDLS is NP-hard too.

Two-stage stochastic programming formulation

This section gives the two-stage stochastic programming model corresponding to the CDLS problem with stochastic refurbishing lead times. The first-stage decisions correspond to the disassembly plan and the disassembly capacity. These decisions are made before the uncertainty is revealed, and they are frozen for the entire planning horizon. The second-stage decisions (control decisions) correspond to the inventory and backlog level which are taken after the uncertainty is revealed. These second stage decisions are computed once the random refurbishing lead times are revealed.

Before introducing the solution approaches proposed to solve the studied problem, the full list of notations used throughout this paper is given in Table 2 .

In this work, we provide a scenario based stochastic optimization formulation. A scenario represents a possible realization of the refurbishing lead time for each component in each period. Classically, in the literature, the set of scenarios contains all possible combination of the L it 's value (e.g., [START_REF] Borodin | Component replenishment planning for a single-level assembly system under random lead times: A chance constrained programming approach[END_REF], such as | | = i ∈N t∈T (L + i -L - i + 1) . Typically, in an assembly problem, an end-item is assembled when all its components arrive. Consequently, the inventory level of the end-items depends on the lead time of each of its components. Conversely, in a disassembly operation, the refurbishing lead time of a component i only affects the inventory and backlog level of component i . Consequently, the expected inventory level can be computed independently for each end-item i , and thus the scenarios can be generated independently for each component i . Definition 3.1. Let i be the set of all possible scenarios for each component i . Each scenario ω corresponds to the realizations from periods 1 to t of the refurbishing lead times of component i . More precisely, if any disassembly operation starts in period t, the component i is available in period t + L ω it in scenario ω, and the probability of scenario ω is P ω i . If the refurbishing lead times in each period are independent, P ω

i = t∈T P i (L it = L ω it ) .
Definition 3.2. Let A ω it be the set that regroups the indices of quantities that arrive at the least before the end of period t, i.e., A ω

it := τ ∈ T | τ + L ω iτ ≤ t .
Such an approach drastically reduces the total number of sce-

narios to | i | = t∈T (L + i -L - i + 1) , ∀ i ∈ N .
The CDLS problem with stochastic refurbishing lead times can be expressed with the following two-stage stochastic MILP model.

E (T

C) = min t∈T i ∈N ω∈ i P ω i h i H ω it + b i B ω it + s t .Y t + u t .O t (1)
H ω it -B ω it = I i 0 + τ ∈A ω it R i Z τ - t τ =1 D iτ i ∈ N , t ∈ T , ω ∈ i (2) Z t ≤ K.Y t t ∈ T (3) F .Z t ≤ C t + O t t ∈ T (4) Z t , O t ≥ 0 t ∈ T (5) Y t ∈ { 0 , 1 } t ∈ T (6) H ω it , B ω it ≥ 0 i ∈ N , t ∈ T , ω ∈ i ( 7
)
where

K = max i ∈N ∀ i ∈N ∀ t∈T D it R i ( 8
)
The objective function (1) minimizes the expected value of the total cost over the planning horizon. Constraints (2) define the inventory balance for each component i at the end of period t under scenario ω. Constraints (3) guarantee that a setup cost is generated in period t if any disassembly operation needs to be performed in that period, where the big K is equal to the maximum cumulative demand of all components over the planning horizon (see Eq. ( 8) ). Constraints (4) give the disassembling capacity constraints in each period t. Finally, constraints ( 5) -( 7) provide the domain of the decision variables.

Note that the model ( 1) -( 7) is flexible enough to represent the situation where order crossover is allowed or forbidden. In the case without order crossover, no part of a component disassembled in period t can be available before the parts of a component disassembled in period t with t < t.

To better explain the notion of non-crossing orders, we present a small example (with a single component and a planning horizon of 4 time periods). The lead times realisation L 1 t = (3 , 1 , 1 , 1) correspond to a case with order crossover since the components disassembled in period t = 1 are received after the reception period of the component disassembled in period t = 2 . On the contrary, if the value of refurbishing lead times L 1 t are (2,1,1,1), there is no order crossover since the components disassembled in period t = 1 are received at the same time as the component disassembled in period t = 2 . To account for the situation where order crossover does not occur, the set of scenarios must be generated such the reception period ( t + L ω it ) of any order passed in period t is later than the reception period ( t + L ω it ) of any order passed in a period t earlier than t.

Solution approaches

The stochastic model describes the uncertain refurbishing lead time through all possible scenarios. However, the size of the scenario set increases exponentially with the number of periods, and the resulting MILP requires too much memory and time resources to solve. Consequently, we investigate two techniques to reduce the number of scenarios, namely Monte-Carlo (MC) sampling and the scenario aggregation (SA).

Monte-Carlo simulation and mixed integer linear programming (MC-MILP)

Model ( 1) -( 7) is difficult to solve because the size | i | of the scenario set increases exponentially with the number of periods. Therefore, we introduce the well-known Monte-Carlo (MC) sampling to estimate the expected total cost [START_REF] Fishman | Monte Carlo: Concepts, algorithms, and applications[END_REF]. MC relies on discrete samples to approximate a probability distribution that cannot be calculated analytically. is, for any small value : lim

ϑ i → + ∞ P (| T C -E (T C) | ≥ ) = 1 .
Therefore, the model ( 1) -( 7) can be converted to the following MC-MILP model:

T C = min ∀ t∈T ∀ i ∈N ω∈ ϑ i P ω i h i H ω it + b i B ω it + s t Y t + u t O t (9)
subject to (3) -( 6) and:

H ω it -B ω it = I i 0 + τ ∈A ω iτ R i Z τ - t τ =1 D iτ i ∈ N , t ∈ T , ω ∈ ϑ i (10) H ω it , B ω it ≥0 i ∈ N , t ∈ T , ω ∈ ϑ i (11)
The objective in ( 9) is now an estimation of the exact value in (1) . Constraints ( 10) and ( 11) provide an estimation of the inventory level for each component i at the end of period t based on the random samples.

Scenario aggregation (SA) approach

This section describes an aggregation strategy that reduces the number of scenarios as much as possible. The SA approach creates a reduced scenario set that keeps the stochastic properties of the original. This approach provides an efficient scenario reduction technique for stochastic programming problems involving both discrete and continuous uncertain parameters, and it can assist researchers to solve the two-level and/or multi-product disassembly/assembly systems under uncertain lead times. The basic idea is to analyze the characteristics of the scenarios, to identify the essential features that affect the decisions, and to remove redundant scenarios.

After releasing the disassembly operation of EoL product in a period τ (with t > τ ), each component is received after a random refurbishment lead time L iτ . As the support of the refurbishing lead time distribution is finite, the probability of receiving an order passed before period t -L + i is one, and the probability of receiving an order passed after period t -L - i is zero. Consequently, the distribution of the backlog and inventory level depends only on the refurbishing lead times of orders passed within a specific time interval. Fig. 2 represents all possible cases of component reception.

In cases 1 and 3, the disassembly occurs in period τ | τ ≤ t -L + i , and component i is received before the end of period t. In case 4, the disassembly occurs after period t -L - i , and the resulting components i are not available in period t. In the case 2, as the disassembly occurs in the interval [ t -

L + i + 1 , t -L - i ] , the components can arrive in the time interval [ τ + L - i , τ + L + i ] ,
and the inventory and backlog levels depend on the realization of lead-times. Definition 4.1. Let R t be the set of disassembly periods that may be involved in the calculation of the backlog or inventory level in period t. A period τ ∈ R t verifies:

R t := τ ∈ T | t -L + i + 1 ≤ τ ≤ t -L - i , ∀ i ∈ N .
Remark 4.2. The cardinality of R t is constant, and it does not depend on the value of |T | , i.e.:

|R t | = (L + i -L - i ) , ∀ i ∈ N , t ∈ T .
Definition 4.3. Let α iτ t be a stochastic parameter indicating if the quantity of component i resulting from a disassembly in period τ arrives before period t. That is,

α iτ t = 1 if L iτ ≤ t -τ 0 if L iτ > t -τ.
Proposition 4.4. Eq. ( 2) are equivalent to the following inventory equations for all components i , periods t, and scenarios ω:

H ω it -B ω it = I i 0 + t-L + i τ =1 R i Z τ + t-L - i τ = t-L + i +1 α ω iτ t R i Z τ - t τ =1 D iτ . ( 12
)
The parameter α ω iτ t (see Definition 4.3 ) defines the quantities of EoL to be disassembled at period τ and received before the end of period t under scenario ω. Thus, ∀ i ∈ N , ∀ t ∈ T and ∀ ω ∈ i , we have:

A ω it := { τ ∈ T | τ + L ω iτ ≤ t} := { τ ∈ T | α ω iτ t > 0 } .
The inventory equation shown in constraints (2) can be replaced by:

H ω it -B ω it = I i 0 + t τ =1 α ω iτ t R i Z τ - t τ =1 D iτ ∀ i ∈ N , t ∈ T , ω ∈ i . ( 13
)
The quantities ( Z τ ) disassembled in period τ | τ ≤ t -L + i are refurbished in period t with probability 1 (see cases 1 and 3 in Fig. 2 ). Thus, for any scenario ω, component i , and period t:

α ω iτ t = 1 ∀ τ ≤ t -L + i . ( 14
)
Similarly, the quantities ( Z τ ) disassembled in period τ | τ > t -L - i are not available in period t with probability 1 (see the cases 2 and 4 in Fig. 2 ). Thus, for any scenario ω, component i , and period t:

α ω iτ t = 0 ∀ τ > t -L - i . (15)
Finally, the components of an EoL disassembled in period τ | ≥ t -L + i + 1 can be received before or after t (see case 2 in Fig. 2 ).

For these periods, α ω iτ t may take value 1 or 0. More precisely, for any scenario ω, component i , and period t :

α ω iτ t ∈ { 0 , 1 } ∀ τ ≥ t -L + i + 1 . ( 16
)
We can easily deduce expression (12) by setting α ω iτ t in Eq. ( 13) to the values given in ( 14) -( 16) .

We explain below how Proposition 4.4 helps to identify the scenario that leads to the same amount of inventory or backlog for a component in a period. Here, the idea is to propose an aggregation strategy in order to reduce the number of possible scenarios by regrouping similar scenarios to form only one.In fact, the inventory and backlog levels of component i in period t are equal for any two scenarios ω 1 and ω 2 such that α ω 1 iτ t = α ω 2 iτ t , ∀ τ ∈ R t . Therefore, for each component i and period t, such scenarios ω 1 and ω 2 can be aggregated.

The resulting set it of aggregated scenarios contains all possible values for the α iτ t 's with τ ∈ R t . Clearly, P (α iτ t = 1) = P (L iτ ≤ t -τ ) , and P (α iτ t = 0) = P (L iτ > t -τ ) . Therefore, the occurrence probability P φ it of scenario φ ∈ it corresponds to the joint probability of these refurbishing lead times. In particular, if the refurbishing lead times in each period are independent, the occurrence probability P φ it of scenario φ ∈ it is:

P φ it = τ ∈R t α φ iτ t × P (L iτ ≤ t -τ ) + (1 -α φ iτ t ) × P (L iτ > t -τ )
When the probability distribution is more complex to express (e.g., the case where order crossover is forbidden), the distribution can be described with the set , and the joint probability distribution can be computed with P φ it = ω∈ φ P ω , where P ω is the probability of scenario ω, and φ is the set of scenarios in with the correct value for the indicator α iτ t ( α ω iτ t = α φ iτ t ) :

φ := { ω ∈ | L ω iτ ≤ t -τ ∀ τ with α φ iτ t = 1 and L iτ > t -τ ∀ τ with α φ iτ t = 0 }
Definition 4.5. The cardinality of the set of the aggregated scenar-

ios ( it ) is | it | = 2 |R t | , ∀ i ∈ N , t ∈ T .
More precisely, it contains all possible values for the vector

[ α i (t-L + i +1 ) t , . . . , α i (t-L - i ) t ] .
Proposition 4.6. The cardinality of the set of the aggregated scenarios ( it ), without considering the crossover of orders:

| it | = L + i -L - i + 1 ∀ i ∈ N , t ∈ T Proof. If orders cannot cross, α iτ -1 t ≥ α iτ t .
That is, if the order passed in period τ arrived in period t ( α iτ t = 1 ), then the order passed in period τ -1 arrived ( α iτ -1 t = 1 ). As a consequence, the vector [ α i (t-L + i +1 ) t , . . . , α i (t-L - i ) t ] contains only 1 at the beginning followed by only 0:

[ α i (t-L + i +1 ) t , . . . , α i (t-L - i ) t ] = [1 , . . . , 1 , 0 , . . . , 0] . ( 17
)
There is only

L + i -L - i + 1 such vectors.
Based on the scenario aggregation (SA) technique and Proposition (4.4) , a model denoted as SA-MILP is proposed to solve the studied problem with and without orders crossover. Therefore, the model ( 1) -( 7) is equivalent to the following SA-MILP model:

E (T C) = min t∈T i ∈N φ∈ it P φ it h i H φ it + b i B φ it + s t Y t + u t O t (18)
subject to (3) -( 6) and:

H φ it -B φ it = t-L + i τ =1 R i Z τ + τ ∈R t α φ iτ t R i .Z τ + I i 0 - t τ =1 D iτ i ∈ N , t ∈ T , φ ∈ it (19) H φ it -B φ it > 0 i ∈ N , t ∈ T , φ ∈ it (20)
Although the SA approach reduces significantly the number of scenarios, this number remains exponential in the refurbishing lead time range when the orders can cross in time. In the next section, we show that SA-MILP can be coupled with MC sampling to solve large instances. Note that in the case where order crossover is forbidden, the number of aggregated scenarios generated by the SA technique is small enough to directly solve large instances (without the need to sample with MC).

Monte-Carlo simulation and scenario aggregation (MC-SAMILP)

A straightforward approach to couple the scenario aggregation with MC sampling is to sample a large number of refurbishing lead time scenarios (as presented in Section 4.1 ), before aggregating these scenarios (as presented in Section 4.2 ). However, such an approach prevents the control of the number of resulting scenarios, and thus the computation time would be unpredictable. In this section, we propose to generate a scenario set it by sampling directly from the zero-one matrix ( α iτ t ) and to determine the production plan with the minimal estimated total cost T C . The steps of the MC-SAMILP approach are given below.

• Step 1: Generate random values for α γ iτ t by following the distribution of α iτ t . More precisely, generate a random value X γ iτ t following a uniform distribution in [0 ;1] for each component i , each period t and τ , and each scenario γ . Then, set α γ iτ t to 1 or 0 according to the distribution of α iτ t , as follows:

α γ iτ t = 1 if X γ iτ t ≤ P (α iτ t = 1) 0 otherwise
• Step 2: Set the probability of each scenario γ ∈ it as follows:

P γ i = 1 | it | ∀ i ∈ N , ∀ t ∈ T , ∀ τ ∈ R t .
• Step 3: Estimate the expected total cost expressed in (18) using the following MC-SAMILP model:

T C = min t∈T i ∈N γ ∈ it P γ it h i H γ it + b i B γ it + s t Y t + u t O t (21)
subject to (3) -( 6) and:

H γ it -B γ it = t-L + i τ =1 R i Z τ + τ ∈R t α γ iτ t R i .Z τ + I i 0 - t τ =1 D iτ i ∈ N , t ∈ T , γ ∈ it (22) H γ it , B γ it ≥0 i ∈ N , t ∈ T , γ ∈ it (23)

Sample average approximation (SAA) approach

The sample average approximation algorithm, which is based on the MC simulation technique, is a method that gives a robust solution for a stochastic optimization problem by solving the average sample approximated problems under a large random sample of the stochastic parameters ( Kim, Park, & Lee, 2018b;[START_REF] Kleywegt | The sample average approximation method for stochastic discrete optimization[END_REF].

The stochastic problem can be approximated by a deterministic optimization problem expressed by MC-MILP (respectively MC-SAMILP). However, the costs resulting from the optimization over a sample of scenarios in MC-MILP and MC-SAMILP are biased because the decisions are tailored for these specific scenarios. Therefore, we use the SAA framework developed in [START_REF] Lamiri | A stochastic model for operating room planning with elective and emergency demand for surgery[END_REF] to compute an unbiased approximation of the expected total costs (see expressions ( 1) and ( 18) ). The basic steps of the SAA algorithm based on MC simulation are given in Algorithm 1 .

Algorithm 1: The SAA algorithm for MC-MILP (respectively MC-SAMILP).

1 Input: The model MC-MILP (respectively MC-SAMILP) and the related parameter ϑ i (respectively it ).

2 Step 1 : Solve the MC-MILP (respectively MC-SAMILP) model; 3Step 2 : Let X * and T C * be the optimal solution and the 'estimated' optimal cost obtained by each model, respectively;

4 Step 3 : Generate B i large set of independent random samples of L 1 it , . . . , L B i it (respectively B it of α 1 iτ t , . . . , α B it iτ t ); 5
Step 4 : Evaluate the exact cost of the optimal solution T C B for the MC-MILP using (24) (respectively ( 25), for the MC-SAMILP).

6 Output: The approximate cost T C B for each model.

Algorithm 2: Rolling horizon (R-H) sampling method.

1 Generate the sets of μ it optimization scenarios with Monte Carlo: α 1 iτ t , . Observe the arrival of components sent to repair. That is,

the value α o iτ t unfolds ∀ i ∈ N , τ ≤ t 7 Modify the MC-SAMILP model to set Z t = Z o t , Y t = Y o t , O t = O o t , and α ω iτ t = α o iτ t ∀ i ∈ N , τ ≤ t, 1 ≤ ω ≤ μ it 8
Solve MC-SAMILP model Corollary 4.7. The expected total cost expressed in ( 1) can be approximated by the following SAA model:

T C B = t∈T i ∈N ω∈B i P ω i h i H ω it + b i B ω it + s t .Y t + u t O t ( 24
)
subject to ( 4) , ( 5) and ( 10) , ( 11) , where B i is a large set of independent random samples of L 1 it , . . . ,

L B i it such that B i ⊂ i , ∀ i ∈ N .
Corollary 4.8. The expected total cost in Eq. ( 18) can be approximated by the following SAA model:

T C B = ∀ t∈T ∀ i ∈N γ ∈B it P γ it h i .H γ it + b i .B γ it + s t .Y t + u t .O t ( 25
)
subject to ( 4) , ( 5) and ( 19) , ( 20) , where B it is a large set of independent random samples of α 1 iτ t , . . . , α B it iτ t , ∀ i ∈ N .

We note that the models presented in Corollaries 4.7 and 4.8 do not contain the equations related to the binary indicator of disassembly variables formulated in (3) and ( 6) . Here, we consider an evaluation of the disassembly plans provided by the MC-MILP and MC-SAMILP models.

Rolling horizon framework

As mentioned earlier, the proposed MILP, MC-MILP, SA-MILP, and MC-SAMILP models are developed to deal with the static decision framework. In the static decision framework, the disassembly decisions are made in period 1, and they are frozen. In this section, we develop a rolling horizon (R-H) strategy to deal with the dynamic decision framework. In the dynamic decision framework, the disassembly decisions are updated in each period t of the horizon after observing the value of the refurbishing lead time up to period t.

We extend the SA-MILP approach to revise the decisions at each period. More precisely, the simulation model (MC-SAMILP) optimizes the disassembly plan in period 0 by considering the first H period, and the decisions of period 0 are implemented. Then, the value of RLT in period 0 is revealed, and the inventory levels and on-hand inventory can be updated. The model is updated to account for the initial and on-order inventory, and the disassembly plan is re-optimized on the horizon from period 1 until the last period.

Note that the T C considered in this simulation is not the objective function of the optimization model, but the sum over time of the true cost incurred by the resolution of the MC-SAMILP model over a true scenario. In the next section, all the proposed approaches are evaluated in order to verify their effectiveness.

Numerical results

This section presents the results of experiments conducted to study the behavior of the proposed optimization approaches. The proposed models are implemented in C with Concert Technology and solved with IBM CPLEX 12.5 on a PC with processor Intel (R) Core TM i7-5500 CPU @ 2.4 GHz and 8 Go RAM under Windows 10 Professional.

Instances generation

The instance generation methods extends the one used by [START_REF] Kim | Capacitated disassembly scheduling with random demand[END_REF] to account for the specificity of the CDLS with stochastic refurbishing lead times. Table 3 shows the generation method for each parameter, where D ∼ U (a, b) indicates that the parameter is randomly generated following a discrete uniform distribution with support { a, . . . , b} .

The numerical experiments are performed with two test beds. The first contains small and medium size instance to compare the SA approach with the classical stochastic MILP model, and the second contains large instances to show the scalability and the limits of the proposed methods. The small and medium size test bed consists of eight sets of instances. Each set contains one randomlygenerated problem. Table 4 gives the refurbishing lead time range ( L + i -L -), the number of components and the number of periods of each instance set. The large test beds consists of 50 problems for Table 3 Characteristics of data sets [START_REF] Kim | Capacitated disassembly scheduling with random demand[END_REF].

Parameters h i s t I i 0 b i D it Values D ∼ U (5 , 10) D ∼ U (50 0 , 10 0 0) D ∼ U (20 , 100) D ∼ U (10 0 , 20 0) D ∼ U (50 , 200) Parameters R i F C t ut L it Values D ∼ U (1 , 4) D ∼ U (1 , 4) D ∼ U (280 , 480) D ∼ U (20 , 40) D ∼ U (L - i , L + i ) Table 4 Set combinations. Set L it |N | |T | 1 D ∼ U (1, 4) 10 5 2 D ∼ U(1, 4) 20 5 3 D ∼ U(1, 3) 10 7 4 D ∼ U(1, 3) 20 7 5 D ∼ U(6, 10) 10 15 6 D ∼ U(6, 10) 20 15 7 D ∼ U(1, 6) 10 20 8 D ∼ U(1, 6) 20 20
each combination of 3 levels of the number of components (10, 20 and 30) and 3 levels of the number of periods (10,20, and 30). For each level of the number of periods, 4 sizes of refurbishing lead time range are considered (1, 6, 9 and |T | ).

The size of the random samples ϑ i and it takes value in [2,10,30,50,150,200]. For each ϑ i and it , we perform 10 independent runs to provide 10 Monte-Carlo optimum solutions X * and their 'estimated' optimal costs T C * . The SAA evaluation of X * uses 10 4 random independent samples. However, to keep the computational time reasonable, the instances with 30 components and 30 periods are evaluated with 10 3 random samples. The time limit is set to one hour, and we report the following performance measures throughout the numerical study:

• I * : the number of optimal solutions (out of 50 instances) obtained by CPLEX solver in one hour; • CPU (s): the average CPU times in second; • G * : the average integrality gap given by CPLEX;

• G + : the average percentage deviation between the estimated total cost ( T C ) and the expected total cost ( E (T C) ). More precisely, G + is calculated as follows:

G + = T C -E (T C) E (T C)
× 100 (26)

Effectiveness of the scenario aggregation approach

We start our experiments by showing the good computational performance of the proposed SA-MILP model compared to the MILP model. Table 5 reports the performance of SA-MILP and MILP on the small and medium-size testbed, whether order crossover is considered or not. More precisely, for each instance set, Table 5 reports the number of scenarios, and the CPU times. For some instances, the MILP model is too large to be generated, and it is indicated with a "-" in such cases.

As the number of scenarios in the stochastic MILP is exponential in the number of periods, the MILP consumes too much memory for instances with more than 10 periods. On the contrary, the number of scenarios remains small using our SA technique. As a consequence, the SA method requires significantly less computation time. While MILP requires several minutes to solve instances with less than 10 periods, SA-MILP solves these instances in less than a second. Besides, the special case where orders do not cross in time is even easier to solve with SA-MILP since CPLEX requires less than 0.50 seconds for this special case versus 9.50 seconds for the general case. Note that, for this test bed, the optimality gap G * returned by CPLEX is null for all instances.

Tables 6 and7 report the performance of MILP and SA-MILP on the large testbed, when the orders can and cannot cross in time, respectively. As previously stated, the MILP returns an optimal solution only for small-size instances with 10 periods when the refurbishing lead time range does not exceed 1. Table 6 provides further evidence of the effectiveness of the scenario aggregation technique for the case without order crossover since SA-MILP solves all instances in few seconds. For instances with order crossover, Table 7 shows that the performance of SA-MILP depends on the size of the support of the refurbishing lead times probability distribution, the number of periods, and the number of components. Among these parameters, the lead-time range has the most significant impact on the computation time, and once the size of the interval of uncertainty passes a certain value (i.e., L + i -L - i ≥ 9 ), SA-MILP cannot solve any instance with 30 periods and/or components.

Effectiveness of scenario aggregation combined with scenario sampling (case with order crossovers)

The numerical results reported in Section 5.2 show that SA-MILP outperforms MILP. Besides, SA-MILP solves large instances in a few seconds when order crossover is forbidden. However, SA-MILP cannot solve instances with large refurbishing lead time ranges when orders can cross. Therefore, the present section investigates the performance of SA-MILP for the case with crossover, when the total costs are approximated with Monte-Carlo Sampling, as presented in Eqs. (1) and (18) . Figs. 3-5 illustrate the average convergence of the estimated values towards the exact values under different number of samples. On the one hand, Figs. 3-5 show that for ϑ i (respectively it ) equal to 200, the MC-MILP approach (respectively MC-SAMILP) can generate a good approximation of the stochastic problem for all tests on average. In addition, the results show that the proposed MC-SAMILP model results in a good approximation with less scenarios than the MC-MILP. For example, for problems with 10 periods, the curves in Fig. 3 shows that the MC-SAMILP provides a good approximation with 50 scenarios, whereas MC-MILP requires at least 100 scenarios.

The optimization results of the MC-MILP and MC-SAMILP obtained under 200 random samples are detailed in Table 8 . This table indicates that the computation time of the MC-MILP approach decreases with the size of the refurbishing lead time range. Regarding the method combining Monte-Carlo simulation and scenario aggregation approach (MC-SAMILP), the same table clearly shows that this approach can solve all problems within a reasonable time frame.

To investigate further the quality of the approximation obtained with optimization, in Table 9 , we analyze the percentage deviations G + between the optimal cost obtained with SA-MILP, and the solution obtained with the sampling approach MC-MILP (respectively, MC-SAMILP). For instances with a large number of components, a large number of periods, and a large refurbishing lead times range, the model based on scenario aggregation (SA-MILP) becomes too large to be generated. In such a case, SA-MILP cannot provide a solution, the gap cannot be computed, and this case is indicated by a " * ". The most remarkable observation is that the percentage deviation between MC-MILP and SA-MILP solutions does not exceed 1.76% on average for all tests, while that of MC-SAMILP

0 - 0 - 0 - 0 - 0 - 0 - 30 0 - 0 - 0 - 0 - 0 - 0 - Mean 0 - 0 - 0 - 2 5
2.5 25 20.9 25 34.5

does not exceed 1.05% on average. This observation confirms that MC-SAMILP provides better approximations than MC-MILP. Note that, the optimality gap G * returned by CPLEX is null for all instances.

To investigate the computation time and the average optimality deviation obtained by solving the SAA approach for both MC-MILP and MC-SAMILP, the sampling results are presented in Table 10 . Note that the costs for MC-MILP (respectively, MC-SAMILP) are unbiased since they are evaluated through a simulation over 10 4 / 10 3 scenarios different from the scenarios used for optimization. Table 10 shows that the SAA can guarantee the convergence of the exact solution with a G+ and a resolution time that do not exceed 0.5% and 110 seconds on average for the two models, respectively. From all these analyses, we conclude that the model that combines scenario aggregation and Monte-Carlo simulation (MC-SAMILP) is quite efficient and can provide a good approximation with 200 independent random samples. Therefore, we can state that our study provides additional support for solving the stochastic CDLS problem and validates the usefulness of the SA, MC and SAA approaches.

+ i -L - i G + G + G + G + G + G + (
L it |N | |T | | it | E (T C) | it | T C [ 1 ,

Effectiveness of R-H method (case with order crossovers)

We continue our experiments by highlighting the good computational performance of the proposed R-H method compared to SA-MILP model. Table 11 reports the performance of these two approaches on the small and medium-size testbed for the case with order crossovers. More precisely, for each instance set, Table 11 reports the number of scenarios, the expected total costs E (T C) , the estimated total costs T C and the percentage deviation between them. Results from Table 11 suggest that the disassembly plan obtained by the R-H strategy is significantly better than the disassembly plan provided by a two-stage stochastic programming model. We note that cost reduction mainly comes from the decrease in the backlog costs obtained when the R-H approach makes dynamic and corrective decisions.

Impact of the time between order

The time Between Order (TBO) plays an important role on the computational complexity of lot-sizing problems. In Section 5.1 , numerical results are presented based on the benchmarks of [START_REF] Kim | Capacitated disassembly scheduling with random demand[END_REF] that contains instances with a TBO equal to 1 (an order at each period). In this section, these benchmarks are modified to obtain a TBO greater than or equal to 2. We modify the ratio related to setup cost, inventory holding costs and backlog costs. Eq. ( 27) presents the TBO formula adapted from [START_REF] Godichaud | Economic order quantity for multistage disassembly systems[END_REF], 2019) ;[START_REF] Pour-Massahian-Tafti | New models and efficient methods for single-product disassembly lot-sizing problem with surplus inventory decisions[END_REF] . We define N as the set of components and s , R , h , b , D as the mean of the discrete uniform distribution used to randomly generate setup, yield, holding inventory, backlog and demand, respectively. In addition, the instance of [START_REF] Kim | Capacitated disassembly scheduling with random demand[END_REF] consider rather small values for u t (cost of outsourcing) and rather large values for the backlog cost b i . In this section, we investigate the impact of instances with values more in line with the lot-sizing literature.

The parameters given in Table 12 are generated in such a way that TBO ≥ 2. The time needed to obtain optimal solutions presented in Table 13 is significantly reduced. This is explained by the fact that increasing TBO reduces the number of orders whom the number of scenario depends.

Benefit of stochastic modelling

In this subsection, we seek to assess the practical performance of the stochastic programming model (SA-MILP) by comparing it with three deterministic models which completely ignores uncertainty. These deterministic versions consider the minimum, average, and maximum lead time, respectively. More precisely, in the first version, the value of L i corresponds to L - i ( L i = L - i ). In the second version, the value of L i is the expected refurbishing lead time ( L i = E (L i ) ). Finally, in the most conservative model, the value L i is set to the maximum lead time L + i ( L i = L + i ). The parameters given in Table 12 and the MILP model (1-8) are used to evaluate the "exact" costs of the solutions obtained with this deterministic model. Results from Table 14 show that the solution provided by the stochastic optimization approach significantly outperforms the deterministic methods, since the stochastic approach reduces the costs by up to 31.14%

Conclusion

In this paper we have proposed a stochastic model for the capacitated disassembly lot sizing problem (CDLS) that can allow the orders to cross in time. The time necessary to refurbish the components is considered as a stochastic variable under the interval representation of uncertainty. The planning problem identifies how much end-of-life (EoL) product to disassemble during each period in order to minimize the expected total cost. Our contribution can be summarized as follows: first, a formulation of the stochastic CDLS problem by a two-phase stochastic mathematical program is proposed. The first-stage decisions correspond to the quantity of EoL product that has to be disassembled as well as the extra capacity to be added in each period. The secondstage decisions correspond to the inventory and the backlog levels. This model describes the random parameter through the set of all possible scenarios. Second, due to the large number of scenarios, a scenario aggregation (reduction) technique is developed. The effectiveness of the SA approach has been successfully tested and proved, especially for the case when order crossover is forbidden. Third, an approach combining a Monte-Carlo simulation and mixed integer programming (respectively scenario aggregation) is proposed in order to study the convergence of the expected cost. Fourth, an almost optimal solution is obtained from a SAA based on the Monte-Carlo optimization method. Fifth, to deal with the dynamic decision frameworks, we introduce a rolling horizon approach.

Our investigations into this area are still ongoing. First, most studies treat the problem under uncertainty on two level DLS problems. Since in real life, most DLS are multi-levels, a promising future research area could be solving the multi-level CDLS problem. Second, studying the uncertainty of demand and/or yield is also an important consideration. Finally, most existing research focuses only on developing effective and efficient approaches for solving stochastic models. However, it is important to evaluate the

Table 14

Performances of the stochastic approach. effectiveness and the robustness of the solution methods. Performance indicators for evaluation should be introduced in future research.

L i = L - i L i = E (L i ) L i = L + i #

Fig. 1 .

 1 Fig. 1. Two level disassembly system.

Fig. 2 .

 2 Fig. 2. Aggregation technique.

  Calculate the T C over the planning horizon.

Fig. 3 .

 3 Fig. 3. Convergence of problems with 10 periods on average.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Convergence of problems with 20 periods on average.

  2 . sN . R . h . D ≤ T BO ≤ 2 . s N . R . h . D h + b i b i (27)

Table 1

 1 Summary of relevant literature.

	Authors	Optimisation approach	#levels	Cap.	PC	Uncertainty
	Gupta & Taleb (1994)	RMRP	Multi
	Taleb et al. (1997) Lee & Xirouchakis (2004)	RMRP Heuristic	Multi Multi		√	--
	Inderfurth & Langella (2006) Kim et al. (2006b) Kim et al. (2006a)	Heuristic Heuristic LH	Two Multi Multi	√ √		Yield --
	Kim, Lee, & Xirouchakis (2008) Kim et al. (2009) Kim & Xirouchakis (2010)	DP Branch and bound LH	Two Multi Two	√ √		--Demand
	Kim & Lee (2011) Prakash et al. (2012) Ji et al. (2016) Godichaud, Amodeo, & Hrouga (2015)	Heuristic CBSA LH GA	Multi Multi Two Multi	√	√ √	----
	Fang et al. (2017)	LH	Two			Demand
	Kim et al. (2018a)	Heuristic	Multi			-
	Godichaud & Amodeo (2018)	EOQ model	Two			-
	Liu & Zhang (2018)	Outer-approximation	Two			Demand, yield
	Tian & Zhang (2019)	PSO	Two			-
	Slama et al. (2019b) Piewthongngam, Chatavithee, & Apichottanakul (2019)	S-MIP MILP	Two Two		√	DLT -
	Godichaud & Amodeo (2019)	EOQ model	Multi			-
	Pour-Massahian-Tafti et al. (2019a)	AF, DF	Two			-
	Pour-Massahian-Tafti, Godichaud, & Amodeo (2019b)	Heuristics	Two			-
	Pour-Massahian-Tafti, Godichaud, & Amodeo (2019c) Slama et al. (2020a)	Heuristics MIP	Two Multi	√	√	--
	Pour-Massahian-Tafti, Godichaud, & Amodeo (2020a) Pour-Massahian-Tafti et al. (2020b) Slama et al. (2020c) Slama et al. (2021b)	EOQ model Heuristics SAA GA	Two Two Two Two	√ √	√	--DLT DLT
	Slama, Ben-Ammar, & Dolgui (2021a) Current paper	Analytical model, recursive formula SA-MILP	Multi Two	√		DLT RLT
	Cap., Capacity; PC, Parts Commonalities; NP, Number of Disassembled products; RMRP, Reverse Material Requirements Planning; LP, Linear Pro-
	gramming; DP, Dynamic Programming; LS-RMRP, Lot-Sizing-RMRP; CBSA, Constraint-Based Simulated Annealing; LH, Lagrangian Heuristics; EOQ,
	Economic Order Quantity; AF, Aggregate Formulation; DF, Disaggregate Formulation; GA, Genetic Algorithm; SA-MILP, Scenario Aggregation-Mixed
	Integer Linear Programming.

Table 2

 2 Notation. , ∀ i ∈ N ω index for scenario ω for each component i , ω ∈ i probability of scenario ω for component i , with ω∈ i P ω i = 1 , ∀ i ∈ N , ∀ ω ∈ i positive inventory level of component i at the end of period t under scenario ω. H ω it is equal to I ω it if I ω it > 0 and zero otherwise B ω it negative inventory level of component i at the end of period t under scenario ω. B ω

	Indexes
	t	index of period t, ∀ t ∈ T
	i index of component i Parameters
	T	set of time periods of the planning horizon
	N	set of components
	i	set of possible scenarios for each component i
	R i	disassembly yield of component i
	D it	external demand for component i in period t
	I i 0	starting inventory of component i
	L ω it	refurbishing lead time for component i in period t under scenario ω
	h i	per-period inventory holding cost of one unit of component i
	s t	setup cost in period t
	b i	per-period backlog cost of one unit of component i
	F	disassembly operation time
	C t	available capacity in period t
	ut	cost of adding a unit of extra capacity in period t
	K	a large number
	Functions
	E (. )	expected value
	P ω
	First-stage decision variables
	Z t	disassembly quantity ordered in period t
	Y t	binary indicator of disassembly in period t
	O t	disassembly overtime in period t
	Second-stage decision variables
	I ω it H ω	inventory level at the end of period t

i it it is equal to -I ω it if I ω it < 0 and zero otherwise

subject to :

Table 5

 5 SA performance.

					Without order crossovers		With order crossovers
					MILP		SA-MILP	MILP				SA-MILP
	L it	|N |	|T |	| i |		CPU(s)	| it |	CPU(s)	| i |		CPU(s)	| it |	CPU(s)
	[ 1 , 4 ]	10	5	1024	9.17	4	0.04	1,024	30.28	8	0.10
	[ 1 , 4 ]	20	5	1024	22.59	4	0.09	1,024	98.65	8	0.15
	[ 1 , 3 ]	10	7	2,187	699.56	3	0.07	2,187	457.03	4	0.17
	[ 1 , 3 ]	20	7	2187	1,790.62	3	0.10	2,187	2095.56	4	0.24
	[ 6 , 10 ]	10	15	5 15		-	5	0.12	5 15		-		1 6	0.35
	[ 6 , 10 ]	20	15	5 15		-	5	0.18	5 15		-		1 6	5.28
	[ 1 , 6 ]	10	20	6 20		-	6	0.21	6 20		-		3 2	1.42
	[ 1 , 6 ]	20	20	6 20		-	6	0.38	6 20		-		3 2	9.32
	Table 6
	Performances of the SA-MILP: without orders crossover.
			MILP							SA-MILP
	#components	10			20		30		10			20	30
	L + i -L -i		I *	CPU	I *	CPU	I *	CPU	I *	CPU		I *	CPU	I *	CPU
	(a) Problem with 10 periods.
	1		50	59.7	50	210.8	50	552.1	50	0.09	50	0.38	50	0.61
	6		0	-		0	-	0	-	5 0	0.14	50	0.83	50	1.28
	9		0	-		0	-	0	-	5 0	0.17	50	0.31	50	1.5
	10		0	-		0	-	0	-	5 0	0.33	50	1.5	50	1.73
	Mean		12.5	14.9	12.5	52.7	12.5	138	50	0.2		50	0.8	50	1.3
	(b) Problem with 20 periods.
	1		0	-		0	-	0	-	5 0	0.31	50	0.34	50	0.44
	6		0	-		0	-	0	-	5 0	0.83	50	1.08	50	1.98
	9		0	-		0	-	0	-	5 0	1.19	50	1.63	50	2.63
	20		0	-		0	-	0	-	5 0	2.16	50	5.17	50	6.19
	Mean		0	-		0	-	0	-	5 0	1.12	50	2.1	50	2.9
	(c) Problem with 30 periods.
	1		0	-		0	-	0	-	5 0	0.55	50	1.03	50	1.88
	6		0	-		0	-	0	-	5 0	0.75	50	1.59	50	2.20
	9		0	-		0	-	0	-	5 0	1.23	50	2.97	50	3.42
	30		0	-		0	-	0	-	5 0	6.47	50	13.39	50	29.02
	Mean		0	-		0	-	0	-	5 0	2.25	50	4.74	50	9.13
	Table 7
	Performances of the SA-MILP: with orders crossover.
		MILP							SA-MILP
	#components	10			20			30		10			20	30
	L + i -L -i	I *	CPU	I *		CPU	I *	CPU	I *	CPU		I *	CPU	I *	CPU
	(a) Problem with 10 periods.
	1	50	510.8	50		829.9	50	2336.5	50	0.3		50	0.07	50	0.4
	6	0	-	0		-	0	-	5 0	2.7		50	8.9	50	35.1
	9	0	-	0		-	0	-	5 0	4.9		50	99.6	50	341.2
	10	0	-	0		-	0	-	5 0	9.5		50	1620.2	0	-
	Mean	12.5	127.7	12.5	207.5	12.5	584.2	50	4.3		50	432.2	37.5	94.2
	(b) Problem with 20 periods.
	1	0	-	0		-	0	-	5 0	0.4		50	0.5	50	0.7
	6	0	-	0		-	0	-	5 0	8.4		50	71.3	50	89.3
	9	0	-	0		-	0	-	5 0	622.3	0	-	0	-
	20	0	-	0		-	0	-	0	-		0	-	0	-
	Mean	0	-	0		-	0	-	37.5	157.7	25	17.9	25	22.5
	(c) Problem with 30 periods.
	1	0	-	0		-	0	-	5 0	1.2		50	1.9	50	2.6
	6	0	-	0		-	0	-	5 0	7.8		50	81.7	50	135.4
	9

Table 8

 8 Performances of the MC simulation.

				MC-MILP					MC-SAMILP
		#components	10		20		30		10		20		30
		L + i -L -i	I *	CPU	I *	CPU	I *	CPU	I *	CPU	I *	CPU	I *	CPU
	Table 9
	Percentage deviation from MC simulation.
		MC-MILP		MC-SAMILP
	#components	10	20	30	10	20	30
	L

Table 10

 10 Performances of the SAA on average.

		MC-MILP					MC-SAMILP
	#components	10		20		30		10		20		30
	periods	CPU	G +	CPU	G +	CPU	G +	CPU	G +	CPU	G +	CPU	G +
	10	24.1	0.62	48.2	0.17	83.1	0.33	22.3	0.29	53.9	0.12	85.4	0.64
	20	53.9	0.43	117.3	0.51	211.7	0.09	69.6	0.35	81.2	0.05	120.3	0.07
	30	88.4	0.04	157.3	0.52	27.2	0.74	77.7	0.64	136.3	0.26	22.9	0.17
	Mean	55.5	0.36	107.6	0.4	107.3	0.38	56.5	0.42	90.46	0.14	76.3	0.29

Table 12

 12 Characteristics of data sets (TBO ≥2). U (0 . 3 , 0 . 5) D ∼ U (350 0 , 450 0) D ∼ U (20 , 100) D ∼ U (0 . 6 , 1) D ∼ U (0 , 160)

	Parameters	h i	s t	I i 0	b i	D it
	Values D ∼ Parameters R i	F	C t	ut	L it
	Values	D ∼ U (1 , 3)	D ∼ U (1 , 4)	D ∼ U (280 , 480)	D ∼ U (150 , 200)	D ∼ U (L -i , L + i )

Table 13

 13 Performances of the SA-MILP according to TBO.

		TBO = 1		TBO ≥ 2
	#components	10	20	30	10	20	30
	periods	CPU	CPU	CPU	CPU	CPU	CPU
	10	0.15	0.18	0.34	0.10	0.21	0.26
	20	0.54	0.76	1.75	0.46	0.578	1.56
	30	22.0	25.65	42.69	5.17	10.79	23.89
	Mean	7.57	8.87	14.93	1.92	3.86	8.57