Mohand Lounes Bentaha
email: mohand.bentaha@univ-lyon2.fr

Alexandre Dolgui

Olga Battaïa

Robert J Riggs

Jack Hu

Orcid Mohand

Lounes Bentaha

Profit-oriented partial disassembly line design: dealing with hazardous parts and task processing times uncertainty

Keywords: product recovery, disassembly, line design and balancing, uncertainty, Monte Carlo sampling, sample average approximation

This paper addresses the problem of profit-oriented disassembly line design and balancing considering partial disassembly, presence of hazardous parts and uncertainty of task processing times. Few papers have studied the stochastic disassembly line balancing problem and existing approaches have focused on heuristic and metaheuristic methods. Most existing work has concentrated on complete disassembly where task times are assumed to be normal random variables and where AND/OR graphs are not considered. The objective of this paper is the design of a serial line that obtains the maximum revenue and then balances the workload under uncertainty. The processing time of a disassembly task is assumed to be a random variable with any known probability distribution. An AND/OR graph is used to model the precedence relationships among tasks. Stochastic programming models and exact-based solution approaches combining the L-shaped algorithm and Monte Carlo sampling techniques are proposed. The relevance and applicability of the proposed models and solution methods are shown by solving efficiently a set of disassembly problem instances from the literature.

Introduction

Product recovery is becoming more important due to its social, environmental and economic benefits. Product recovery preserves resources by reducing the consumption of virgin raw materials, water and energy. In addition, this process plays a key role in minimising the amount of waste sent to landfills and diminishing air and water pollution [START_REF] Ashby | Materials and the Environment, Eco-informed Material Choice[END_REF]. Through recycling, remanufacturing and reuse, product recovery aims to retrieve valuable parts and materials from discarded or End of Life (EOL) products. A mandatory step and most challenging part of product recovery process is disassembly.

Disassembly is a revalorising process, a methodical and organised separation of parts and materials of outdated products for recycling, remanufacturing and reuse. A disassembly line is more suitable than a single workstation or a disassembly cell to carry out disassembly operations with a higher productivity rate and automated disassembly. Figure 1 illustrates forward and reverse logistics flows and positions disassembly in the closed loop logistics [START_REF] Ma | Disassembly Process Planning Algorithms for End of Life Product Recovery and Environmentally Conscious Disposal[END_REF]. This study particularly focuses on design of disassembly systems and more specifically on disassembly line design.

Although disassembly seems to be the reverse of assembly, the disassembly process has unique characteristics which make it more complex than assembly. [START_REF] Gupta | Product Recovery Using a Disassembly Line: Challenges and Solution[END_REF] provide a comparison of operational and technical considerations of assembly and disassembly lines. Indeed, the most obvious difference is the flow process which is convergent for assembly and is divergent for disassembly. In a disassembly environment, a product is broken down into many parts and subassemblies whose qualities, quantities and reliabilities cannot be controlled as in an assembly environment. The assembly process has to be complete while the disassembly process, due to technical and economic restrictions, does not have to be carried out completely and hence disassembly is usually a partial process [START_REF] Lambert | Determining Optimum Disassembly Sequences in Electronic Equipment[END_REF]. For example, irreversible connections of components of a product can be considered as a technical restriction with the disassembly cost being greater than the revenue obtained from retrieved parts as an economic restriction. The structure and quality of EOL products are very uncertain and even the number of components in such products cannot be predicted. Consumers may remove certain components or parts before disposing of the product or these components may be unfastened during their period of use. Moreover, an EOL product may contain certain hazardous parts necessitating special handling at a workstation of a disassembly line. These unique complex characteristics make disassembly processes and especially disassembly lines more challenging. Therefore, Figure 1. Forward and reverse logistics [START_REF] Ma | Disassembly Process Planning Algorithms for End of Life Product Recovery and Environmentally Conscious Disposal[END_REF].

efficient decision-making tools are needed in order to optimise their performance and cost effectiveness. Such tools must take into account the high degree of uncertainty in the structure and the quality of the products to be disassembled.

This paper considers in particular the design and balancing of disassembly lines under uncertainty. Such a line consists of an ordered sequence of workstations connected by a material handling system, which allows the transportation of work-pieces (parts and subassemblies) from one workstation to another. Certain parts or subassemblies may be hazardous and require a particular treatment incurring a supplementary cost. As mentioned above, the disassembly process is typically partial. The optimisation problem considered aims to assign a given set of disassembly operations (tasks) to an ordered sequence of workstations while respecting precedence and cycle time constraints under uncertainty. The case of partial disassembly is studied and the objective is to maximise the profit produced by the line. The line profit is calculated as the difference between the positive revenue generated by recovered parts, and the line operation cost taken as negative revenue. The latter includes the workstation operation costs, additional costs of workstations handling hazardous parts of the EOL product, as well as penalty costs caused by the cycle time constraint violations. This objective allows the minimisation of the number of workstations of the line, which in turn minimises the total idle time. Then, for this minimum number of workstations, an optimal balance is determined under uncertainty for the resulting disassembly line. The obtained balance ensures that idle times are as similar as possible for all the line workstations. This scheme of lexicographic optimisation is motivated for economic reasons. In fact, the revenue generated by EOL products recovery constitutes one of the most important attractions for investment in disassembly lines. Thus, the criterion of line profit is considered as dominant for disassembly lines.

In order to help designers to design economically viable disassembly lines and therefore to facilitate the implementation of disassembly end of life scenarios in practice, this paper develops an efficient procedure to design disassembly lines under uncertainty brought by the EOL products. The added value of this work is the ability to assess the disassembly profit or cost for a product at the end of life. Such an assessment may help to take the decisions not only for EOL options but even at the product design stage. The proposed model is rather complete and is capable of treating at the same time partial disassembly, presence of hazardous material and uncertainty of task processing times.

The remainder of the paper is organised as follows: Section 2 provides an overview of the relevant literature on disassembly line design and balancing. Section 3 follows with a formal description of the considered optimisation problem. Section 4 presents the developed solution algorithm that combines the L-shaped method with Monte Carlo sampling techniques.

Section 5 is dedicated to the analysis of numerical experiments and Section 6 concludes the paper and suggests further research directions.

Previous work

The disassembly line design problem was introduced by [START_REF] Güngör | Disassembly Line Balancing[END_REF]. This problem called DLBP consists of the assignment of disassembly tasks to workstations of a disassembly line with the objective of optimising some measures of effectiveness (time, labour or money).

To deal with the deterministic DLBP, heuristic and metaheuristic approaches were developed. An iterative heuristic using branch and bound technique was developed in [START_REF] Lambert | A Heuristic Solution for the Disassembly Line Balancing Problem Incorporating Sequence Dependent Costs[END_REF] to deal with the line balancing problem subjected to sequence dependent costs. A multi-objective heuristic for U-shaped DLBP was developed in [START_REF] Avikal | A Heuristic for U-shaped Disassembly Line Balancing Problems[END_REF]. The authors considered several performance criteria in a lexicographic order: minimise the workstation idle times, maximise the priority of removing hazardous components and maximise the priority of removing high-demand components. Two multiobjective metaheuristics, a distributed agent ant system and an uninformed deterministic search, for the design and balancing of disassembly lines were developed and compared in [START_REF] Mcgovern | Uninformed and Probabilistic Distributed Agent Combinatorial Searches for the Unary NPcomplete Disassembly Line Balancing Problem[END_REF]. Other multi-objective formulations of the DLBP were presented in [START_REF] Mcgovern | Ant Colony Optimization for Disassembly Sequencing with Multiple Objectives[END_REF], [START_REF] Ding | A New Multi-objective Ant Colony Algorithm for Solving the Disassembly Line Balancing Problem[END_REF] where the objectives are also ordered lexicographically. Ant colony-based approaches were developed for this problem. Several metaheuristic approaches were developed to solve problems of disassembly process planning, disassembly scheduling and sequence dependent disassembly line balancing, respectively, a Petri Net-based approach [START_REF] Tiwari | A Petri Net-based Approach to Determine the Disassembly Strategy of a Product[END_REF], constraint-based simulated annealing approach [START_REF] Prakash | Constraint-based Simulated Annealing (CBSA) Approach to Solve the Disassembly Scheduling Problem[END_REF] and particle swarm optimisation, a tabu search, ant colony optimisation [START_REF] Kalayci | A Tabu Search Algorithm for Balancing a Sequence-dependent Disassembly Line[END_REF]. A beam search-based approach for DLBP was proposed in [START_REF] Mete | A Solution Approach Based on Beam Search Algorithm for Disassembly Line Balancing Problem[END_REF] to minimise the number of workstations. [START_REF] Tang | ASystematicApproach to Design and Operation of Disassembly Lines[END_REF] developed a Petri net approach where a heuristic was employed to maximise the line productivity. Qualitative and quantitative comparisons of different heuristics and metaheuristics for DLBP, including genetic algorithm, ant colony optimisation method, greedy algorithm, greedy/hill-climbing, greedy/2-opt hybrid heuristics and hunter-killer heuristics, were undertaken in [START_REF] Mcgovern | Combinatorial Optimization Analysis of the Unary NP-complete Disassembly Line Balancing Problem[END_REF].

Mathematical programming formulations and exact solution approaches were also proposed for the DLBP. [START_REF] Altekin | Profit-oriented Disassembly-line Balancing[END_REF] developed an integer programming formulation for profit maximisation for the case of partial disassembly. [START_REF] Koc | Two Exact Formulations for Disassembly Line Balancing Problems with Task Precedence Diagram Construction Using an AND/OR graph[END_REF] considered DLBP with the objective of minimising the number of workstations. Two exact approaches based on mixed integer and dynamic programmes were developed. To balance a mixed-model disassembly line, a linear physical programming approach was proposed in [START_REF] Ilgin | Disassembly Line Balancing Using Linear Physical Programming[END_REF]. A comprehensive survey on production line design and balancing particularly on disassembly can be found in the literature [START_REF] Ilgin | Environmentally Conscious Manufacturing and Product Recovery (ECMPRO): A Review of the State of the Art[END_REF][START_REF] Battaïa | A Taxonomy of Line Balancing Problems and Their Solution Approaches[END_REF]Bentaha, Battaïa, and Dolgui 2015a).

Few studies in the literature, however, have studied the DLBP under uncertainty (Bentaha, Battaïa, and Dolgui 2015a). In [START_REF] Riggs | Disassembly Line Balancing under High Variety of End of Life States Using a Joint Precedence Graph Approach[END_REF], the joint precedence graph approach, known from assembly line balancing, was adapted for disassembly line balancing in order to deal with EOL product states. A fuzzy coloured Petri net model with a heuristic solution method was proposed by [START_REF] Turowski | Disassembly Line Design with Uncertainty[END_REF] to deal with the human factors (uncertainty of disassembly task times) and product condition. A collaborative ant colony algorithm for stochastic mixed-model U-shaped disassembly line balancing was developed in [START_REF] Agrawal | A Collaborative Ant Colony Algorithm to Stochastic Mixed-model U-shaped Disassembly Line Balancing and Sequencing Problem[END_REF]. Task times were assumed stochastic with known normal probability distributions. A self-guided ants metaheuristic was proposed in [START_REF] Tripathi | Real World Disassembly Modeling and Sequencing Problem: Optimization by Algorithm of Self-Guided Ants (ASGA)[END_REF] for the disassembly line sequencing problem, where fuzzy optimisation model was developed with the objective of maximising the net revenue of the disassembly process. [START_REF] Tuncel | Solving Large Scale Disassembly Line Balancing Problem with Uncertainty Using Reinforcement Learning[END_REF] used a Monte Carlo-based reinforcement learning technique to solve the multi-objective DLBP under demand variations of the EOL products. [START_REF] Güngör | A Solution Approach to the Disassembly Line Balancing Problem in the Presence of Task Failures[END_REF] proposed a heuristic to deal with task failures caused by defective parts or joints in the EOL product. A MIP-based predictive-reactive approach to deal with task failures was also developed in [START_REF] Altekin | Task-failure-driven Rebalancing of Disassembly Lines[END_REF]. In these two studies, a task failure can be seen as a disruption and the authors proposed remedial actions to deal with such short-term disturbances. A first study to deal dynamically with real time production line disturbances is proposed in [START_REF] Antoine | Iterated Local Search for Dynamic Assembly Line Rebalancing Problem[END_REF]. A binary bi-objective non-linear programme was developed in Aydemir-Karadag and Turkbey (2013) for disassembly line design and balancing under uncertainty of the task times. Disassembly task times were assumed to be independent random variables with known normal probability distributions. Complete disassembly was considered and a genetic algorithm was designed to solve the problem.

This literature review shows that there are few papers that have dealt with stochastic DLBP and these are either restricted to the study of demand fluctuations, task failures, condition of the EOL products or stochastic task times (as normal random variables) with only heuristic/metaheuristic solution methods being proposed without information about the solution quality. Many of these works have not considered the case of partial disassembly and have not used an AND/OR graph to model the precedence relationships among tasks. It is shown in the literature [START_REF] Koc | Two Exact Formulations for Disassembly Line Balancing Problems with Task Precedence Diagram Construction Using an AND/OR graph[END_REF] that the integration of the AND/OR graphs in the DLBP formulation allowed better solutions to be obtained compared to the use of AND precedence diagrams. Most of the AND/OR graphs used in the literature are part-based diagrams, and such graphs do not exploit the precedence relations among tasks and parts of the EOL product. The existing solution methods are unadapted to simultaneously take into account other sources of uncertainty, and only consider task times that are normally distributed.

To bridge the gap, we developed in Bentaha et al. (2014), Bentaha, Battaïa, and Dolgui (2014a, 2014b, 2014c, 2014d), Bentaha, Battaïa, and Dolgui (2015b), Bentaha et al. (2015) mathematical models and exact solution approaches to design disassembly (or assembly) lines under uncertainty of the task processing times. Partial disassembly has been taken into account and complex and/or precedence relationships among tasks have been integrated.

In Bentaha et al. (2014), Bentaha, Battaïa, and Dolgui (2014a, 2014b, 2014c), uncertainty was modelled using workstation expectation times. This problem was introduced in Bentaha et al. (2014). In Bentaha, Battaïa, and Dolgui (2014a), the joint problem of disassembly line design and tasks sequencing was studied. In Bentaha, Battaïa, and Dolgui (2014c), a Lagrangian relaxation was proposed to maximise the disassembly line profit. In Bentaha, Battaïa and Dolgui (2014b), the line balancing problem has been undertaken under the assumption of the fixed number of workstations. In Bentaha, Battaïa, and Dolgui (2015b), Bentaha et al. (2015), the goal was to guarantee a certain operational level defined by the designer. Uncertainty was modelled using joint probabilistic cycle time constraints. In Bentaha, Battaïa, and Dolgui (2014d), uncertainty was modelled using the notion of recourse cost. In other words, this study aimed in minimising the line stoppage costs caused by the task processing time uncertainties. In Bentaha, Battaïa, and Dolgui (2014d), the objective was the line cost minimisation. The line cost includes the operation costs for workstations and penalty costs generated by the cycle time constraint violations. Thus, revenue generated by recovered parts was not considered. However, as mentioned in Section 1, the revenue generated by EOL products recovery constitutes one of the most important attractions for investment in disassembly lines. In addition, handling hazardous parts and balancing the workload of the resulted line in Bentaha, Battaïa, and Dolgui (2014d) were not considered.

The present paper deals with the profit-oriented partial DLBP under uncertainty of the disassembly task times (SP-DLBP) which are assumed to be random variables with known probability distributions. Task times can be modelled using any probability distribution. To deal with this uncertainty, a stochastic programme is developed. An AND/OR graph is used to model the precedence relations among tasks. The objective is to maximise the net revenue produced by the line. As mentioned in Section 1, the line profit is calculated as the difference between the positive revenue generated by recovered parts, and the line operation cost, considered as negative revenue. The cost includes the workstation operation costs, additional costs of workstations handling hazardous parts of the EOL product, as well as penalty costs incurred by the cycle time constraint violations. This work addresses the aforementioned shortcomings in Bentaha, Battaïa, and Dolgui (2014d).

The optimisation procedure proposed in this paper consists of two phases: (1) a minimum number of workstations for the line is found by the maximisation of the line revenue and (2) for this minimum number of workstations, an optimal balance is determined under uncertainty. An exact-based solution method is developed. The method integrates efficient Monte Carlo sampling techniques with the L-shaped algorithm. The developed solution method computes lower and upper bounds for the problem solutions and provides the corresponding confidence intervals as well as optimality gaps. In addition, the proposed solution approach deals with any probability distribution for task times and can be extended to deal with additional stochastic parameters without additional complexities. A detailed problem formulation is given in the next section.

Problem definition and formulation

In this section, the SP-DLBP and the graph used to model disassembly alternatives of an EOL product and relationships among tasks and subassemblies are defined (Section 3.1). The first phase of the SP-DLBP is formulated in Section 3.2 while the second phase is formulated in Section 3.3.

AND/OR graph and disassembly alternatives

The SP-DLBP consists of the assignment of the disassembly tasks, a set I , of an EOL product to an ordered sequence of workstations, a set J , while satisfying precedence and cycle time constraints under uncertainty of the task times. The cycle time (C T ) is the amount of time allocated to each station to complete its assigned tasks. It is the ratio of the planning period length to the number of products that need to be disassembled in order to meet the demand.

The solution method developed in this paper considers the following assumptions: a single type discarded product is to be partially (or completely) disassembled on a straight paced line; the EOL products are sufficiently available and their demand quantities are deterministic; all received EOL products contain all their parts with no addition or removing of components. Disassembly tasks of every product have to be assigned to the workstations in the same manner. All products are assumed to be of the same type. Certain components of a product are hazardous. Task times are assumed to be random variables with known probability distributions. A disassembly task can be performed by all workstations but only by one at a time.

Each component or subassembly of a product has a certain recovery (resale) value which is considered the same for all components or subassemblies of the same type. This is particularly the case if the resale value of a given component or subassembly is only based on the contained materials. In this paper, the study is focused on the variability in disassembly task processing times. A fixed cost is allocated per operating time unit of a workstation with an additional fixed cost per operating time unit of a workstation handling hazardous parts.

All possible disassembly alternatives of an EOL product are modelled using an AND/OR graph (Homem de Mello and Sanderson 1990; [START_REF] Koc | Two Exact Formulations for Disassembly Line Balancing Problems with Task Precedence Diagram Construction Using an AND/OR graph[END_REF]Bentaha, Battaïa, and Dolgui 2015b). Such a graph represents explicitly the precedence relationships existing among tasks and parts or subassemblies (Bentaha, Battaïa, and Dolgui 2014d;Bentaha et al. 2015). The use of an AND/OR graph is illustrated with a rigid caster example (see Figure 2 and Table 1). The product to be disassembled is composed of nine components as shown in Figure 2. The AND/OR graph in Figure 3 corresponds to the rigid caster example and is constructed as follows: each subassembly is represented by a node labelled A k , k ∈ K . For example, node A 0 represents the rigid caster which can be noted as '1:9', A 1 represents the subassembly '1:5,7:9', etc. (see Table 1). Subassemblies with one component are not represented; set K contains the indices of all possible subassemblies that can be generated by the tasks from I . Each node labelled B i , i ∈ I , represents a disassembly task. For instance, node B 1 represents disassembly task '1', B 2 represents disassembly task '2', etc. Two types of arcs define the precedence relations between the subassemblies and disassembly tasks: AND and OR. The first type imposes a mandatory precedence relation and the second type is employed for optional precedence dependencies where only one option is chosen for a final solution. For example, if a disassembly task generates two subassemblies or more, then it is related to these subassemblies by AND-type arcs (in bold in Figure 3). If several concurrent disassembly tasks may be performed on a subassembly, this subassembly is related to these disassembly tasks by OR-type arcs. A sink node s is introduced and linked with dummy arcs to all disassembly tasks. Thus, if the dummy task s is assigned to a workstation, the disassembly process is finished (partial or complete disassembly). Table 1 summarises all the possible disassembly tasks that can be performed on the rigid caster in Figure 2. For each task, the corresponding generated components and/or subassemblies are given.

First phase: line profit maximisation

As mentioned earlier, disassembly task times t i , i ∈ I, are assumed to be stochastic variables with known probability distributions. They are represented by a random vector ξ = (t 1 , t2 , . . . , tN ) from a set ⊂ R N + . Let ti = t i ( ξ ), i ∈ I . A disassembly task i ∈ H ⊂ I is called hazardous if its execution generates a hazardous subassembly or component. Disassembly tasks have to be assigned to the workstations of the line under precedence constraints given by the AND/OR graph, i.e. exactly one alternative sequence of the disassembly process is chosen, and at the same time the line profit is maximised. In this paper, the line profit is the difference between the positive revenue generated by recovered parts and the line operation cost taken as negative revenue. The line cost includes the workstation operation costs, additional costs of workstations handling hazardous parts and penalty costs incurred by workstation overloads. Each time the duration of tasks assigned to one workstation exceeds the given line cycle time C T , a corrective action is needed, and its cost is included in the objective function.

To formulate the SP-DLBP, the following notation is used. L i : set of retrieved parts by the execution of disassembly task B i , i ∈ I ; A k : a subassembly, k ∈ K ; B i : a disassembly task, i ∈ I ; s: the AND/OR graph's sink node; r l : revenue generated by part l, l ∈ L; F c : fixed cost per operating time unit of a workstation, F c > 0; q: time unit fixed recourse cost; C T : cycle time, C T > 0; C h : additional cost per time unit for stations handling hazardous parts, C h > 0;

P k : set of indices for predecessors of A k , k ∈ K , i.e. P k = {i| B i precedes A k }; S k : set of indices for successors of A k , k ∈ K , S k = {i| A k precedes B i }.

Decision variables

x i j =      1, if disassembly task B i is assigned to workstation j; 0, otherwise. h j =     
1, if a hazardous task is assigned to workstation j; 0, otherwise.

x s j =     
1, if the dummy task s is assigned to workstation j; 0, otherwise. y j ( ξ ) 0, j ∈ J , measures the amount of time exceeding C T if there is any.

For the problem described, the following mathematical model has been developed.

Stochastic mixed integer programme

max    i∈I j∈J l∈L i r l • x i j -C T F c • j∈J j • x s j + C h • j∈J h j + E ξ q • j∈J y j ( ξ )    (I) s.t. i∈S 0 j∈J x i j = 1 (1) j∈J x i j 1, ∀i ∈ I (2) i∈S k j∈J x i j i∈P k j∈J x i j , ∀k ∈ K \{0} (3) i∈S k x iv i∈P k v j=1 x i j , ∀k ∈ K \{0}, ∀v ∈ J (4) j∈J x s j = 1 (5) j∈J j • x i j j∈J j • x s j , ∀i ∈ I (6) h j x i j , ∀ j ∈ J, ∀i ∈ H (7) i∈I t i ( ξ ) • x i j -y j ( ξ ) C T , ∀ j ∈ J (8) (x s j , x i j , h j ) ∈ X ⊆ {0, 1} |J |•(|I |+2) (9) y j ( ξ ) 0, ∀ j ∈ J (10)
The terms of the objective function represent, respectively, the earned profit of retrieved parts, the cost of operating workstations, the additional cost for handling hazardous parts and the expected recourse cost E ξ with respect to the distribution of ξ .

If the dummy task s is assigned to workstation j, then j defines the number of processed stations. Constraint (1) imposes the selection of only one disassembly task, i.e. OR-successor, to begin the disassembly process. Without loss of generality, we assume that at least one task is required for the disassembly process. Naturally, this task has to be selected from the alternative disassembly tasks realisable for the initial EOL product A 0 . Only one alternative has to be selected. The following disassembly options will be determined by the corresponding A-node and B-nodes connected to it. Constraint set (2) indicates that a task is to be assigned to at most one workstation. Constraints (3) ensure that only one OR-successor is selected. Constraint set (4) defines the precedence relations among tasks. Constraint (5) imposes the assignment of the dummy task s to one station. Constraints (6) ensure that all disassembly tasks are assigned to lower or equal-indexed workstations than the one to which s is assigned. Constraints (7) ensure the value of h j to be 1 if at least one hazardous task is assigned to a workstation j. Constraints (8) represent the cycle time limitation. The variable y j ( ξ ), j ∈ J measures the amount of time exceeding C T if there is any at each opened workstation. Finally, sets ( 9) and (10) represent constraints on the possible values of the decision variables.

Let Z be the set of all possible solutions of programme (I). An assignment of the disassembly tasks to workstations of the line is always possible and the number of the possible assignments is finite, i.e. Z = ∅ and finite. It is straightforward that 0 E ξ q • j∈J y j ( ξ ) < ∞.

Since Z is not empty and the objective function defines finite values, then problem (I) possesses an optimal solution. Let x be a vector of decision variables x i j , x s j , h j , ∀i ∈ I, ∀ j ∈ J and X = {x| constraints (1)-( 7) and ( 9) are satisfied}. If ξ has a finite-discrete distribution {(ξ ℓ , p ℓ ), ℓ ∈ D, p ℓ > 0, ∀ℓ ∈ D}, where D = {1, 2, . . . , D}, D ∈ N * and p ℓ is the realisation probability of ξ ℓ of ξ , then programme (I) is an ordinary linear programme with a dual decomposition structure. Programme (II) given below represents a deterministic equivalent version of programme (I) and shows its particular block structure; each programme (11) can be seen as a block. The programme (II) is a two-stage stochastic linear mixed integer programme with fixed recourse [START_REF] Birge | Introduction to Stochastic Programming[END_REF]. This same programme, i.e. (II), will be used in the solution method. max

   i∈I j∈J l∈L i r l • x i j -C T F c • j∈J j • x s j + C h • j∈J h j + Q D (x)    (II) s.t. x ∈ X whereQ D (x) = D ℓ=1 p ℓ • Q(x, ξ ℓ ) and Q(x, ξ ℓ ) = min q • j∈J y j (ξ ℓ )| i∈I t i (ξ ℓ ) • x i j -y j (ξ ℓ ) C T , y j (ξ ℓ ) 0, ∀ j ∈ J , ℓ ∈ D (11)
Depending on the number of realisations of ξ , i.e. D, programme (II) may become very large in scale, but its block structure can be exploited efficiently by specially designed algorithms such as the L-shaped algorithm [START_REF] Ahmed | The Sample Average Approximation Method for Stochastic Programs with Integer Recourse[END_REF]. The L-shaped method, proposed by [START_REF] Van Slyke | L-shaped Linear Programs with Applications to Optimal Control and Stochastic Programming[END_REF], is a variant of Benders' decomposition. The initial problem (I) is decomposed into a master problem and a subproblem. In this paper, the 0-1 binary decision variables, x ∈ X, in the master problem are called master variables and the continuous variables, y j (ξ ), ξ ∈ , ∀ j ∈ J , of the subproblem are called subproblem variables.

Second phase: idle time leveling

In the first phase, a number of workstations m * = j * M (where x s j * = 1) and a subset of tasks I * ⊂ I (a disassembly alternative) are determined under uncertainty of the task times. An example of such a disassembly alternative and a selected subset I * of disassembly tasks is illustrated in Figure 4. The selected alternative (in bold) is represented by an AND-graph and I * = {B 2 , B 9 , B 23 }. There are no more OR-relations since a decision of disassembling partially a product is made. For this reason, only direct precedence relationships among tasks are considered, i.e. generated subassemblies {A 0 , A 2 , A 8 } are simply deleted. The simple graph in Figure 5 is then created.

In the second phase, a balance measure is optimised. This measure seeks to assign the disassembly tasks I * to the m * stations ensuring (if possible) similar idle time at each workstation, i.e. total idle time balancing or levelling. Let J * = {1, 2, . . . , m * }, ST j ( ξ ) = i∈I * t i ( ξ ) • x i j , ∀ j ∈ J * and Pred(i) = {i ′ ∈ I * | i ′ precedes i}, i ∈ I * . The maximum of the expectations of the differences among all workstations workloads is minimised in the non-linear programme (III) as follows.

min max ∀ j, j ′ ∈J * , j = j ′ E ξ ST j ( ξ ) -ST j ′ ( ξ ) (III) s.t. j∈J * x i j = 1, ∀i ∈ I * (12) j∈J * j • x i ′ j j∈J * j • x i j , ∀i ∈ I * , ∀i ′ ∈ Pred(i) (13) x i j ∈ {0, 1}, ∀i ∈ I * , ∀ j ∈ J * (14)
Constraints ( 12) ensures that a task is assigned to only one workstation. Constraint set (13) models the precedence relationships among tasks.

Let x be a vector of decision variables x i j , ∀i ∈ I * , ∀ j ∈ J * and X = {x| constraints ( 12), ( 13) and ( 14) are satisfied}.

Let B = {1, 2, . . . , B}, B ∈ N * and S l ( ξ ) = ST j ( ξ ) -ST j ′ ( ξ ), j, j ′ ∈ J * , j = j ′ , l ∈ B, where B = |J * |
2 . The programme (IV) given below represents an equivalent version of programme (III). min Y (IV)

s.t. x ∈ X Y E ξ S l ( ξ ) , ∀l ∈ B
Note that the recourse cost E ξ q • y j ( ξ ) , j ∈ J * , may change for some workstations, caused by the possible reassignment of disassembly tasks I * in the second phase. Thus, the total recourse cost E ξ q • j∈J y j ( ξ ) may increase or decrease. No reassignment of hazardous tasks will be considered. Let x, x ′ be the optimal solutions of problems (I) and (III), and Q D (x) and Q D (x ′ ) the corresponding recourse costs, respectively. The solution x ′ is retained if

η = Q D (x ′ ) -Q D (x) Q D (x) × 100 ̺
Otherwise, x ′ is not considered and x is retained (̺ is a percentage fixed by the decision-maker). In other words, x ′ is retained if the percentage of the recourse cost's increase in phase 2 does not exceed a certain percentage, i.e. ̺, fixed by the decision-maker.

The solution approach proposed is detailed in the next section.

Solution method

Computing the expected value E ξ q • j∈J y j ( ξ ) for a given task assignment to workstations is quite difficult, since it requires numerical integration of implicitly defined probability density functions of variables y j (ξ ), ξ ∈ , ∀ j ∈ J . Even for the discrete distribution of ξ , the number of linear programmes of type (11) to solve may increase exponentially. Although the exact evaluation of the expectation term in (I) is possible, its optimisation presents serious difficulties [START_REF] Santoso | A Stochastic Programming Approach for Supply Chain Network Design under Uncertainty[END_REF]. Indeed, E ξ q • j∈J y j ( ξ ) is implicitly defined. In order to deal with these difficulties, the so-called Sample Average Approximation (SAA) method is used [START_REF] Kleywegt | The Sample Average Approximation Method for Stochastic Discrete Optimization[END_REF]. It combines Monte Carlo sampling techniques detailed in Section 4.1 with the L-shaped algorithm presented in Section 4.2.1.

Monte Carlo sampling

As introduced earlier, disassembly task times are modelled by a random vector ξ varying over a set ⊂ R N + (of a given probability space ( , F, P) introduced by ξ ). The integral

E ξ q • j∈J y j ( ξ ) = E ξ Q(x, ξ ) = q • j∈J y j ( ξ ) d P represents the expected value E ξ Q(x, ξ ) of the function Q(x, ξ ) = min    q • j∈J y j (ξ )| i∈I t i (ξ ) • x i j -y j (ξ ) C T , y j (ξ ) 0, ∀ j ∈ J    , ξ ∈ A Monte Carlo estimation Q λ (x) of the expected value E ξ Q(x, ξ
) is obtained with a random generation of a λ-sample ξ 1 , ξ 2 , . . . , ξ λ of the random vector ξ , using computer generated pseudo-random numbers. We have then:

Q λ (x) = 1 λ λ ℓ=1 Q(x, ξ ℓ ). Random variable Q λ (x, ξ ) defined by Q λ (x, ξ ) = 1 λ λ ℓ=1 Q(x, ξℓ ) represents the Monte Carlo estimator of E ξ Q(x, ξ ) . It is an unbiased estimator for E ξ Q(x, ξ ) .

Solution method of the first phase

L-shaped algorithm

The main idea of the L-shaped method is to approximate the non-linear term in the objective function of the stochastic programme with fixed recourse (I). The fixed recourse is called simple recourse [START_REF] Kall | Stochastic Programming[END_REF] if the matrix of recourse variable coefficients in programme (11) has the following form: (I, -I), where I represents the identity matrix.

The SP-DLBP has a simple recourse. In fact, we have:

A = {a|a = W y = I y ′ -I y, y, y ′ 0} = R |J | , i.e.
the second stage programme is always feasible for any feasible solution of the first phase and for any realisation scenario ξ of ξ ; Q(x, ξ ) = min q T y | I y ′ -I y = (ξ ) -T (ξ )x, y, y ′ 0 , ∀x ∈ X, ∀ξ ∈ , defines the matrix formulation of the second-stage programme, where q = (q, . . . , q) T , (ξ ) = (C T , . . . , C T ) T and T (ξ ) = (T, . . . , T ) T , T = (t 1 , . . . , t N ), ∀ξ ∈ , represent, respectively, the recourse vector, the cycle time vector and the technological matrix.

The result above (i.e. SP-DLBP has a simple recourse) can be proved as follows: ∀x ∈ X, ∀ξ ∈ , we have

i∈I t i (ξ ) • x i j -y j (ξ ) C T ⇐⇒ i∈I t i (ξ ) • x i j -y j (ξ ) + y ′ j (ξ ) = C T ⇐⇒ y ′ j (ξ ) -y j (ξ ) = C T - i∈I t i (ξ ) • x i j y j (ξ ), y ′ j (ξ ) 0, ∀ j ∈ J
Since y j (ξ ), y ′ j (ξ ) 0, ∀ j ∈ J, ∀ξ ∈ , we have then, ∀ j ∈ J, ∀ξ ∈ , y ′ j (ξ )y j (ξ ) ∈ R. It follows from the above that A = {a|a = W y = I y ′ -I y, y, y ′ 0} = R |J | , i.e. ∀ξ ∈ , ∀x ∈ X, the second stage programme

Q(x, ξ ) = min q • j∈J y j (ξ )| i∈I t i (ξ ) • x i j -y j (ξ ) C T , y j (ξ ) 0, ∀ j ∈ J is feasible.
Since the problem addressed in the L-shaped algorithm is a minimisation one and since max f (x) =minf (x) where f (x) is a function of decision variables x, in order to solve problem (I) its equivalent minimisation version min

   C T F c • j∈J j • x s j + C h • j∈J h j - i∈I j∈J l∈L i r l • x i j + E ξ q • j∈J y j ( ξ )    ( 15
)
is used. The optimal value of problem (I) is then the symmetric value of ( 15). Consider the minimisation version of programme (I) and the matrix formulation of problem (II) given below, where c

T x = C T F c • j∈J j • x s j + C h • j∈J h j -i∈I j∈J l∈L i r l • x i j and ℓ = (C T , . . . , C T ) T , ∀ℓ ∈ D. min c T x + Q D (x), x ∈ X where Q D (x) = D ℓ=1 p ℓ • Q(x, ξ ℓ ) and Q(x, ξ ℓ ) = min q T y | W y = ℓ -T (ξ ℓ )x, y 0
In the L-shaped algorithm given below, counters h and υ are used, respectively, for the optimality cuts and for the algorithm iterations; E ν and e ν are defined in the algorithm. In the objective function of the master problem, a new variable ϕ is introduced. It verifies the inequality ϕ Q D (x). Since Q D (x) is implicitly defined by a relatively large number of optimisation problems, the master programme is not directly solved with this inequality. The L-shaped algorithm processes as follows.

L-shaped algorithm

Step 0. Set h = υ = 0

Step 1. Set υ = υ + 1. Solve the master programme:

min c T x + ϕ s.t. x ∈ X E ν • x + ϕ e ν , ν = 1, 2, . . . , h (16)
x binary, ϕ 0 Let (x υ , ϕ υ ) be an optimal solution.

Step 2. For ℓ = 1, 2, . . . , D, solve the subprogramme:

min W = q T • y s.t. W y = ℓ -T (ξ ℓ ) • x υ y 0
Let ω υ ℓ be the simplex multipliers associated with an optimal solution of problem ℓ above; define

E h+1 = ℓ∈D p ℓ • (ω υ ℓ ) T • T ℓ and e h+1 = ℓ∈D p ℓ • (ω υ ℓ ) T • ℓ Let θ υ = e h+1 -E h+1 • x υ . If ϕ υ θ υ , stop:
x υ is an optimal solution. Else, generate an optimality cut of type ( 16), set h = h + 1, add constraint of type ( 16) and return to Step 1.

Sample average approximation procedure

The SAA strategy deals with an approximate computation of the expected recourse cost. It is used to compute lower and upper bounds of the objective function for problem (I). Therefore, the optimality gap and statistical confidence intervals on the quality of the approximate solutions are evaluated.

The SAA strategy is implemented as follows: for a random sample of ξ , generated using MCS and having size , the term E ξ Q(x, ξ ) is approximated by the sample average function Q (x) = 1 • ℓ=1 Q(x, ξ ℓ ). Thus, problem (I) is approximated by problem (II) where p ℓ = 1 , ℓ = 1, 2, . . . , . Then, problem (II) is solved with the L-shaped algorithm.

Let γ , x and γ * , x * be the optimal objective values and optimal solutions of problems (II) and (I), respectively. We have then,

γ * = min x∈X f (x) = c T x + E ξ Q(x, ξ ) and γ = min x∈X f (x) = c T x + 1 • ℓ=1 Q(x, ξ ℓ )
Note that

E f (x) = c T x + 1 • E ℓ=1 Q(x, ξ ℓ ) = f (x)
The strong law of large numbers implies that γ -→ →∞ γ * with probability 1. Moreover, in [START_REF] Kleywegt | The Sample Average Approximation Method for Stochastic Discrete Optimization[END_REF] it is shown that, under certain weak conditions, as increases, x converges to x * with probability approaching 1 exponentially fast. In order to get a good approximate solution to problem (I), the sample size is determined as

3σ 2 max ε 2 log |X| α ( 17
)
where ε > 0 and α ∈]0, 1[, σ 2 max represents the maximum variance of function differences on Q (x) (Kleywegt, Shapiro, and Homem-De-Mello 2001); X, as introduced in constraint (9) of programme (I), defines the set of the possible values of decision variables x.

This sample size is sufficient to get an ε-solution x of problem (I) with a probability at least equal to (1α), i.e. a solution with an absolute optimality gap ε. Even if it may be too conservative in practice, estimate (17) shows that the required sample size is linear in the number of disassembly tasks and workstations: |X|

2 |J |•(|I |+2) = 2 M×(N+2) and |X| 2 M×(N+2) ⇒ log |X| (log 2) • M × (N + 2).
Practically, the SAA strategy consists of generating independent random samples (ξ n 1 , ξ n 2 , . . . , ξ n ), n = 1, 2, . . . , , of a modest sample size . Then, estimates of lower and upper bounds and the optimality gap depending on and are computed. In this study, and are chosen so as to obtain approximate solutions of a reasonable quality.

SAA procedure Step 1. Lower Bound Estimation

Generate independent random samples (ξ n 1 , ξ n 2 , . . . , ξ n ), n = 1, 2, . . . , of size each. Solve the corresponding problem (II) with the L-shaped algorithm, compute an optimal solution x n and the corresponding objective value γ n , for each value of n. Compute a lower bound: L B = 1 • n=1 γ n . Let X * be a set of the optimal solutions of problem (I). Then,

γ min x∈X * f (x) and E γ E min x∈X * f (x) min x∈X * E f (x) = γ * Since E L B = 1 • E n=1 γ n = E(γ ), we have then, L B γ * . The variance σ 2 L B
of L B is estimated as follows:

σ 2 L B = 1 ( -1) • n=1 (γ n -L B ) 2
By the application of CLT, we have:

L B ❀ N E(γ ), σ L B √ , σ L B = √ Var(γ ). A confidence interval of the level
(1α) for E(γ ) is then given by:

L B - z α/2 • σ L B √ , L B + z α/2 • σ L B √
where z α is the quantile of the standard normal distribution N (0, 1) verifying P N (0, 1) z α = 1α.

Step 2. Upper Bound Estimation Generate a sample (ξ 1 , ξ 2 , . . . , ξ λ ) of size λ independent from the generated samples in Step 1. Let x n be a feasible solution obtained in Step 1; x n should be the one for which f λ (x) is minimum. Since x n is a feasible solution of problem (I), then f (x n ) γ * . Compute the following upper bound:

U B λ = f λ (x n ). U B λ is an unbiased estimator of f (x n ): E f λ (x n ) = f (x n ),
we have then U B λ γ * . The variance σ 2 U B λ of U B λ can be estimated with:

σ 2 U B λ = 1 λ(λ -1) • λ ℓ=1 (c T x n + Q(x n , ξ ℓ ) -U B λ ) 2
By the application of the CLT, we have:

U B λ ❀ N f (x n ), σ U B √ λ , σ U B = Var f λ (x n ) .
A confidence interval of the level (1α) of f (x n ) is then given by:

U B λ - z α/2 • σ U B λ √ λ , U B λ + z α/2 • σ U B λ √ λ , P N (0, 1) z α = 1 -α
Step 3. Optimality Gap Estimation Here, the optimality gap is the difference between upper and lower bound values:

OG λ = U B λ -L B . Note that E OG λ = f (x n ) -E(γ ) f (x n ) -γ * .
It follows that OG λ is a biased estimator of the optimality gap, it overestimates f (x n )γ * . Its bias γ * -E(γ ) is monotonically decreasing in [START_REF] Mak | Monte Carlo Bounding Techniques for Determining Solution Quality in Stochastic Programs[END_REF][START_REF] Norkin | A Branch and Bound Method for Stochastic Global Optimization[END_REF]. The variance σ 2 OG λ of OG λ can be estimated with:

σ 2 OG λ = σ 2 L B + σ 2 U B λ
The flowchart in Figure 6 summarises the steps of the SAA method. The random samples are generated using MCS.

Solution method of the second phase

In this subsection, a lower and upper bounding schemes are defined in order to approximate problem (IV).

Upper bound of programme (IV)

In the same way as in Section 4.1, let S l ( ξ ) = |S l ( ξ )| and consider the random variable

S λ l ( ξ ) = 1 λ λ ℓ=1 S l ( ξℓ ), l ∈ B. Then, S λ l ( ξ ), l ∈ B, represents an unbiased estimator of E ξ S l ( ξ ) , l ∈ B: E ξ S λ l ( ξ ) = E 1 λ λ ℓ=1 S l ( ξℓ ) = 1 λ λ ℓ=1 E S l ( ξℓ ) = E ξ S l ( ξ ) Let z ℓl S l (ξ ℓ ) = |S l (ξ ℓ )|, ℓ = 1, 2, . . . , λ, ∀l ∈ B.
Using a λ-sample (ξ 1 , ξ 2 , . . . , ξ λ ) of the random vector ξ , the optimal value of linear programme (V) below is an approximation of the optimal value of programme (IV).

min Y (V) s.t. x ∈ X Y 1 λ • ℓ=λ ℓ=1 z ℓl , ∀l ∈ B -z ℓl S l (ξ ℓ ) z ℓl , ∀l ∈ B, ℓ = 1, 2, . . . , λ (18)
Y 0, z ℓl 0, ℓ = 1, 2, . . . , λ, ∀l ∈ B Recall that, in phase 2, X = {x| constraints ( 12)-( 14) are satisfied}. Constraint set (18) can be replaced with:

u ℓl + S l (ξ ℓ ) = z ℓl , ∀l ∈ B, ℓ = 1, 2, . . . , λ 0 u ℓl 2z ℓl , ∀l ∈ B, ℓ = 1, 2, . . . , λ
Proposition 1 The optimal value of problem (V) is an upper bound of the optimal value of problem (IV).

Proof. Since z ℓl S l (ξ ℓ ), ℓ = 1, 2, . . . , λ, ∀l ∈ B, then

ℓ=λ ℓ=1 z ℓl ℓ=λ ℓ=1 S l (ξ ℓ ), ∀l ∈ B and Y 1 λ • ℓ=λ ℓ=1 z ℓl 1 λ • ℓ=λ ℓ=1 S l (ξ ℓ ), ∀l ∈ B 4.3.2 Lower bound of programme (IV)
The optimal value of programme (VI) given below defines an approximation of the optimal value of programme (IV).

min max l∈B E ξ S l ( ξ ) (VI) s.t. x ∈ X Proposition 2
The optimal value of problem (VI) is a lower bound of the optimal value of problem (IV).

Proof. Obvious since

E ξ S l ( ξ ) E ξ S l ( ξ ) , l ∈ B
The programme (VI') below represents an equivalent version of programme (VI).

min Y (VI')

s.t. x ∈ X -Y E ξ S l ( ξ ) Y, ∀l ∈ B Y 0
Using a λ-sample (ξ 1 , ξ 2 , . . . , ξ λ ) of the random vector ξ , the expectation E ξ S l ( ξ ) can be approximated with 1 λ λ ℓ=1

S l (ξ ℓ ), l ∈ B.

Numerical results

The first phase, i.e. profit maximisation with the SAA method, and the second one, i.e. idle time levelling with programmes (V) and (VI'), were implemented with C++ (under Linux) and ILOG CPLEX 12.4 was used to solve the different models on a PC with a Pentium(R) Dual-Core CPU T4500, 2.30 GHz and 3GB RAM. These two sequential phases have been applied to the problem instance illustrated in Figure 3 and to six instances available in the literature. These six instances are used as benchmarking problems and contain process alternatives for disassembly. The names of the problem instances were, respectively, composed of the first letters of authors' names and the year of publication, i.e. BBD13 (Bentaha, Battaïa, and Dolgui 2014a), BBD13b (Bentaha, Battaïa, and Dolgui 2014c), KSE09 [START_REF] Koc | Two Exact Formulations for Disassembly Line Balancing Problems with Task Precedence Diagram Construction Using an AND/OR graph[END_REF], L99a and L99b from [START_REF] Lambert | Linear Programming in Disassembly/Clustering Sequence Generation[END_REF]) and MJKL11 from [START_REF] Ma | Disassembly Process Planning Algorithms for End of Life Product Recovery and Environmentally Conscious Disposal[END_REF]. BBD13a corresponds to the rigid caster product instance (see Figures 2 and3). The input data for each problem instance are given in Table 2. The columns 'AND-relations' report the number of disassembly tasks with no successor in subcolumn 0, with one AND-type arc in subcolumn 1 and with two AND-type arcs in subcolumn 2. The column 'arcs' gives the total number of AND-type and OR-type arcs. Table 3 reports the optimisation results of phase 1 obtained for the processed instances using Monte Carlo sampling for different values of and λ. The number of generated samples was 10, the recourse cost q was fixed at 5, operating workstation cost F c at 3, additional cost for stations handling hazardous parts C h = 2, the confidence intervals level was 95%, and 25% of tasks were taken as hazardous. The remaining parameters r l , l ∈ L and ti , i ∈ I were randomly generated:

• r l ∈ [0, r max ], where r l is the revenue of part l, l ∈ L and r max is its maximum value; r l , l ∈ L is generated uniformly in the interval [0, r max ]. The value of r max depends on the instance used. As an example, the r max value of BBD13 is 12.

• ti is assumed to be a random variable with a normal distribution of known mean µ i and standard deviation σ i , where σ i = µ i /5; and ti is the processing time of the disassembly task i, i ∈ I .

Although using different probability distributions for task times is possible, they were assumed, without loss of generality, to be normally distributed. Columns 'L B' and 'U B' report the lower and upper bound values and the corresponding confidence intervals, respectively. Column 'OG' reports the optimality gap OG λ and the associated standard deviation σ OG λ . Finally, columns 'o-tasks', 's-tasks', 'n-stations', 'h-stations', 'CPU time' report the original number of tasks of the selected disassembly alternative, the number of selected tasks of the selected alternative, the optimal number of workstations, the number of hazardous workstations with the corresponding rank for each station, and the resolution time in seconds, respectively.

The numerical results in Table 3 show that a higher chosen value of resulted in a higher quality solution. In particular, solutions of instances MJKL11 and L99b represent good quality solutions since their optimality gap values and the corresponding standard deviations are significantly reduced. The resolution time, however, illustrated in Figure 7 increases considerably as the value of increases. This depends on the complexity of the instance (e.g. the number of tasks with 2 AND-type arcs) but is due mainly to the escalating number of subproblems to solve for each master problem. Figure 8 illustrates the convergence of L B and U B λ values for the instance L99b.

In order to show the ability of the proposed approach to handle different probability distributions for the task times, instances BBD12 and L99a are solved for three different cases as illustrated in Table 4. In the first case, each task time ti , i ∈ I is assumed to follow a triangular distribution T (a i , b i , c i ), i ∈ I where b i = µ i , ∀i ∈ I ; a i , c i , i ∈ I are chosen randomly. In the second one, each task time is assumed to follow a continuous uniform distribution U(a i , c i ), i ∈ I where a i , c i have the same values as for triangular distributions. In the third case, a mix of normal, triangular and uniform distributions is considered. Approximatively, 40% of task times are assumed to follow normal distributions, 40% triangular distributions and 20% uniform distributions.

All the returned solutions were the same as the solutions obtained in Table 3 for instances BBD12 and L99a; CPU times are of the same order of magnitude. The only difference is the objective values of the returned solutions (which are slightly greater) and their associated confidence intervals. These values depend on the used probability distributions of the task processing times.

Table 5 reports the optimisation results obtained for phase 2 using Monte Carlo sampling for different values of λ. The values of the remaining parameters were the same as in phase 1. Columns 'Upper bound' and 'Lower bound' report the optimisation results of programmes (V) and (VI'), respectively. Sub-column 'Reassignment' indicates the value 0 if the • 100. The results in Table 5 show that the proposed lower and upper bounding schemes of the line balancing phase returned solutions of good quality, except for instance BBD13b where the solution is of an average quality. For each value of λ, the returned solution was the same for both the upper and lower bounds. CPU time of the lower bound is better due to the reduced number of constraints and variables of programme (V') compared to programme (V). Except for instance L99b, the solution of phase 1 stays unchanged in phase 2. This can be explained by the fact that each unchanged solution is represented by a simple path (see Figure 5) for its precedence graph. In such a case, no reassignment is possible. All instances were solved in a few seconds due to the elimination of hard cycle time constraints. Cycle time constraints are not needed in phase 2 since the number of workstations is fixed. The selected alternative and the corresponding selected tasks for instance L99b are illustrated in Figure 9. Note that this selected alternative is represented by a tree. Table 6 establishes a comparison of the solution quality (in terms of line profit) for different values of λ between the solution of instance L99b in phase 1, i.e. x, and its solution in phase 2, i.e. x ′ . The column 'Estimated recourse cost' reports the recourse cost of workstations 1, 2, 3 in subcolumns 'R W1 ', 'R W2 ', 'R W3 ', respectively and 'R total ' represents the total recourse cost. Variable η and parameter ̺, introduced in Section 3.3, indicate the percentage recourse cost increases in phase 2 and the percentage the recourse cost is allowed to increase, respectively.

The comparison between x and x ′ is illustrated in Figure 10, whose objective is to prove the importance and the added value of phase 2. As shown in Table 6, solution x ′ of phase 2 is notably better than solution x of phase 1. In fact, the three workstations are better balanced with x ′ than with x, and the recourse cost in the second phase is decreased by about 85% when compared to the recourse cost from the first phase. Equivalently, the recourse cost of phase 2 represents only 15% of the recourse cost of phase 1. As a consequence, η ≪ ̺ and x ′ is to be retained as the best and final solution. CPU time increases insignificantly compared to the increasing value of λ, which is due to the continuity of the linear programmes solved for each value of λ. Indeed, if λ = 100, then 100 continuous linear programmes of the same type as in step 2 of the L-shaped algorithm are solved.

Conclusion and future research directions

In this paper, profit-oriented disassembly line design and balancing problems were studied under uncertainty. The case of partial disassembly and the presence of hazardous parts are considered and disassembly task times were assumed to be random variables with known probability distributions. To solve the defined problems with an assessment of solution quality, two phases were addressed. In the first phase, a two-stage stochastic linear mixed integer programme with fixed recourse and the SAA strategy were proposed for line design. The SAA method combined the L-shaped algorithm with Monte Carlo sampling techniques. This procedure provided lower and upper bounds for the optimal value of the optimisation problem solved and associated 95% confidence intervals. The corresponding optimality gaps and deviations were also provided. In the second phase, a min-max formulation and an upper bound mixed integer programme were developed in order to balance, under uncertainty, the line derived from the first phase.

The SAA procedure from the line design phase and the lower and upper bound MIPs from the line balancing phase were evaluated by a set of instances from the literature. The obtained results of phase 1 have shown that high-quality solutions require large Monte Carlo samples and a resolution time which increases considerably with the increase of sample sizes. All instances were solved in less than four hours, except one which was solved in just under seven hours with our previously specified machine. Since the installation of a disassembly line is usually based on a long-term decision, a resolution time of several hours for the problem of disassembly line design and balancing is acceptable. In the line balancing phase, all problem instances were solved in very little time. The numerical results have shown that the solution from phase 1 remained optimal in phase 2 for the majority of problem instances excluding one instance for which the line balance could be improved, resulting in a decrease of the recourse cost. This proved the added value of the line balancing phase.

The obtained results show the applicability of the developed optimisation model in real disassembly context. Indeed, the computational time is acceptable enough for giving to the designer the opportunity to generate different design alternatives depending on the profit expected from the retrieved parts. This model helps to make a decision on the disassembly alternative to be realised in the line and on the assignment of these tasks to workstations. Therefore, the choice between complete or partial disassembly can be made on the basis of the economic arguments. The modelling process presented can be easily adapted for real life cases like End of Life Vehicles or Waste Electrical and Electronic Equipment (WEEE). Undertaking such case studies is one of our next research objectives.

The presented study opens numerous further research directions for problem modelling and algorithmic enhancements. We believe that the computing time at the first solution stage can be reduced by generating new optimality cuts, using a computational grid in order to parallelise the optimisation process or by implementing penalty function techniques. The

Figure 2 .

 2 Figure 2. A rigid caster example.

Figure 3 .

 3 Figure 3. A rigid caster AND/OR graph.

3. 2 . 1

 21 Sets and parameters I : set of disassembly task indices: I = {1, 2, . . . , N}, N ∈ N * ; J : set of workstation indices: J = {1, 2, . . . , M}, M ∈ N * ; H : hazardous disassembly task index set; L: set of part indices: L = {1, 2, . . . , L}, L ∈ N * ; K : set of indices for the generated subassemblies: K = {0, 1, . . . , K}, K ∈ N;

Figure 4 .

 4 Figure 4. A selected disassembly alternative.

Figure 5 .

 5 Figure 5. Obtained simple AND-graph.

Figure 6 .

 6 Figure 6. Flowchart of the SAA procedure.

Figure 7 .

 7 Figure 7. CPU time of phase 1 (in seconds).

Figure 8 .

 8 Figure 8. U B and L B convergence for instance L99b.

Figure 9 .

 9 Figure 9. Selected alternative and disassembly tasks for instance L99b.

Figure 10 .

 10 Figure 10. Idle time levelling and recourse cost decreasing for instance L99b: (a) tasks assignment and profit maximisation in phase 1; (b) tasks reassignment and idle time levelling in phase 2.

Table 1 .

 1 The rigid caster associated disassembly tasks and the corresponding subassemblies and/or components.

	Task	Subassembly	Components	Task	Subassembly	Components
	1	1:5,7:9	6	17	1,3,9	8
	2	1:6,8,9	7	18	1,3,8	9
	3	1:7,9	8	19	1,3,9	7
	4	1:8	9	20	1,2,7	3;4;5;9
	5	1,3,8,9	2;4;5;7	21	1,3,8	2;4;5;7
	6	1:5,7,9	8	22	1,2,7	3;4;5;8
	7	1:5,7,8	9	23	1,3,9	2;4;5;6
	8	1,3,8,9	2;4;5;6	24	1,2,6	3;4;5;9
	9	1:6,9	8	25	1,3,8	2;4;5;6
	10	1:6,8	9	26	1,2,6	3;4;5;8
	11	1:5,7,9	6	27	1,2,7	6
	12	1:6,9	7	28	1,6,6	7
	13	1,2,6,7	3;4;5;9	29	-	1;3;9
	14	1:5,7,8	6	30	-	1;3;8
	15	1:6,8	7	31	-	1;2;7
	16	1,2,6,7	3;4;5;8	32	-	1;2;6

Table 2 .

 2 Instances from the literature used as benchmark problems in the numerical experiments.

	AND-relations

Table 3 .

 3 Obtained results: profit maximisation using normal probability distributions for the task times.

		, λ	L B	U B	OG	o-tasks s-tasks n-stations h-stations CPU time
	MJKL11	100, 200 300, 600 600, 900 800, 1200 1000, 1500 34.027 ± 0.034 34.214 ± 0.192 33.753 ± 0.088 33.942 ± 0.096 33.921 ± 0.093	34.032 ± 0.384 33.865 ± 0.127 33.826 ± 0.085 33.782 ± 0.064 34.063 ± 0.051	-0.182, 2.785 0.112, 1.591 -0.116, 1.307 -0.139, 1.139 0.036, 1.006	7 7 7 7 7	5 5 5 5 5	2 2 2 2 2	(1 : 1) (1 : 1) (1 : 1) (1 : 1) (1 : 1)	101 738 2836 4506 7122
	BBD13a	100, 200 300, 600 600, 900 800, 1200 180.434 ± 0.003 180.423 ± 0.296 -0.011, 5.227 180.400 ± 0.010 180.420 ± 1.778 0.020, 12.831 180.425 ± 0.008 180.446 ± 0.592 0.021, 7.396 180.428 ± 0.003 180.424 ± 0.394 -0.004, 6.037 1000, 1500 180.425 ± 0.002 180.439 ± 0.236 0.014, 4.675	4 4 4 4 4	2 2 2 2 2	1 1 1 1 1	(1 : 1) (1 : 1) (1 : 1) (1 : 1) (1 : 1)	44 333 1103 1973 3057
	L99a	100, 200 300, 600 600, 900 800, 1200 500.176 ± 0.146 499.462 ± 0.863 -0.714, 15.262 0.446, 37.461 499.168 ± 0.402 499.614 ± 5.191 1.542, 21.591 499.628 ± 0.184 501.170 ± 1.727 0.280, 17.624 500.096 ± 0.148 500.376 ± 1.151 1000, 1500 500.136 ± 0.141 500.403 ± 0.691 0.267, 13.649	9 9 9 9 9	7 7 7 7 7	3 3 3 3 3	-----	133 1104 3852 7791 11309
	BBD13b	100, 200 300, 600 600, 900 800, 1200 207.859 ± 0.463 204.477 ± 0.462 -3.382, 8.209 208.532 ± 1.056 209.080 ± 2.779 0.548, 20.125 2.100, 11.599 204.841 ± 0.694 206.941 ± 0.924 207.779 ± 0.476 205.562 ± 0.618 -2.217, 9.491 1000, 1500 207.089 ± 0.268 207.443 ± 0.369 0.354, 7.317	4 4 4 4 4	3 3 3 3 3	2 2 2 2 2	-----	65 549 1863 3069 4826
	KSE09	100, 200 300, 600 600, 900 800, 1200 895.806 ± 0.050 895.910 ± 1.471 896.117 ± 0.142 895.766 ± 8.842 -0.351, 63.802 895.988 ± 0.068 895.582 ± 2.942 -0.406, 36.774 895.692 ± 0.053 895.901 ± 1.961 0.209, 30.018 0.104, 25.992 1000, 1500 895.680 ± 0.053 895.823 ± 1.176 0.143, 23.246	6 6 6 6 6	4 4 4 4 4	2 2 2 2 2	(1 : 2) (1 : 2) (1 : 2) (1 : 2) (1 : 2)	53 457 1393 2415 3845
	L99b	100, 200 300, 600 600, 900 800, 1200 1000, 1500 75.216 ± 0.007 75.261 ± 0.006 75.204 ± 0.008 75.206 ± 0.003 75.232 ± 0.005	75.055 ± 0.742 75.197 ± 0.247 75.246 ± 0.165 75.225 ± 0.123 75.223 ± 0.099	-0.206, 5.353 -0.007, 3.085 0.040, 2.518 -0.007, 2.180 0.007, 1.950	8 8 8 8 8	6 6 6 6 6	3 3 3 3 3	-----	235 2422 8731 14913 24439
	BBD12	100, 200 300, 600 600, 900 800, 1200 1000, 1500 92.912 ± 0.002 92.917 ± 0.006 92.914 ± 0.004 92.910 ± 0.004 92.913 ± 0.003	92.895 ± 0.923 92.894 ± 0.307 92.916 ± 0.205 92.906 ± 0.153 92.932 ± 0.123	-0.022, 6.658 -0.020, 3.837 0.006, 3.132 -0.007, 2.712 0.020, 2.426	3 3 3 3 3	3 3 3 3 3	2 2 2 2 2	(1 : 2) (1 : 2) (1 : 2) (1 : 2) (1 : 2)	40 328 1103 1981 3205
	solution of phase 2 remains the same as in phase 1 and indicates the value 1 if this solution is different. The last column gives
	the optimality gap defined by Gap = U B-L B L B

Table 4 .

 4 Obtained results: profit maximisation using different probability distributions for the task times.

		, λ	LB	UB	OG	o-tasks s-tasks n-stations h-stations CPU time
	L99a	9 9 600,900 506.797 ± 0.108 506.482 ± 1.151 -0.315,17.614 9 100,200 507.069 ± 0.136 508.520 ± 5.188 1.451,37.436 300,600 507.086 ± 0.113 507.597 ± 1.726 0.511,21.578 9 800,1200 506.924 ± 0.089 507.072 ± 0.863 0.148,15.252 1000,1500 506.912 ± 0.066 506.772 ± 0.690 -0.140,13.640 9	7 7 7 7 7	3 3 3 3 3	-----	134 1088 3874 7704 11370	Triangular distribution
	BBD12	100,200 93.251 ± 0.005 93.278 ± 0.923 0.027,6.658 300,600 93.276 ± 0.003 93.274 ± 0.307 -0.002,3.837 600,900 93.269 ± 0.002 93.278 ± 0.205 0.009,3.132 800,1200 93.271 ± 0.002 93.261 ± 0.153 -0.010,2.712 1000,1500 93.264 ± 0.001 93.266 ± 0.123 0.002,2.426	3 3 3 3 3	3 3 3 3 3	2 2 2 2 2	(1:2) (1:2) (1:2) (1:2) (1:2)	42 354 1141 1853 3287
	L99a	9 9 9 800,1200 505.230 ± 0.114 505.063 ± 0.863 -0.167,15.255 9 100,200 504.650 ± 0.370 506.233 ± 5.189 1.583,37.448 300,600 504.987 ± 0.156 505.917 ± 1.727 0.930,21.583 600,900 505.207 ± 0.120 505.489 ± 1.151 0.282,17.617 1000,1500 504.737 ± 0.104 505.224 ± 0.690 0.487,13.644 9	7 7 7 7 7	3 3 3 3 3	-----	132 1070 3911 7650 11251	Uniform distribution
	BBD12	100,200 93.431 ± 0.008 93.447 ± 0.923 0.016,6.658 300,600 93.443 ± 0.005 93.443 ± 0.307 0.000,3.837 600,900 93.434 ± 0.002 93.446 ± 0.205 0.012,3.132 800,1200 93.431 ± 0.002 93.435 ± 0.153 0.004,2.712 1000,1500 93.438 ± 0.002 93.445 ± 0.123 0.007,2.426	3 3 3 3 3	3 3 3 3 3	2 2 2 2 2	(1:2) (1:2) (1:2) (1:2) (1:2)	39 343 1186 1897 3258
	L99a	100,200 503.173 ± 0.462 503.806 ± 5.190 0.633,37.456 300,600 503.434 ± 0.220 503.936 ± 1.727 0.502,21.588 600,900 502.822 ± 0.169 503.877 ± 1.151 1.055,17.623 800,1200 502.990 ± 0.174 503.139 ± 0.863 0.149,15.260 1000,1500 503.034 ± 0.163 503.412 ± 0.691 0.378,13.648	9 9 9 9 9	7 7 7 7 7	3 3 3 3 3	-----	132 1125 3888 7662 11298	Mix of distributions
	BBD12	100,200 92.948 ± 0.008 92.902 ± 0.923 -0.046,6.658 300,600 92.956 ± 0.008 92.977 ± 0.307 0.021,3.837 600,900 92.960 ± 0.003 92.961 ± 0.205 0.001,3.132 800,1200 92.949 ± 0.004 92.954 ± 0.153 0.005,2.712 1000,1500 92.949 ± 0.002 92.964 ± 0.123 0.015,2.426	3 3 3 3 3	3 3 3 3 3	2 2 2 2 2	(1:2) (1:2) (1:2) (1:2) (1:2)	45 356 1118 1938 3160

Table 5 .

 5 Obtained results: idle time leveling.

				Upper bound			Lower bound		Gap(%)
			U B	Reassignment	CPU time λ	L B	Reassignment	CPU time
		500	16.03	0	0.016	15.86	0	0.001	1.07
	MJKL11	1000	16.18	0	0.047	15.74	0	0.020	2.79
		1500	15.82	0	0.031	15.82	0	0.020	0
		500	0	0	0.000	0	0	0.001	0
	BBD13a	1000	0	0	0.000	0	0	0.001	0
		1500	0	0	0.000	0	0	0.001	0
		500	7.679	0	26.22	7.092	0	0.500	8.27
	L99a	1000	7.523	0	90.57	7.011	0	0.08	7.30
		1500	7.644	0	269.19	6.709	0	0.06	13.94
		500	22.95	0	0.078	17.95	0	0.001	27.85
	BBD13b	1000	23.00	0	0.172	18.88	0	0.020	21.82
		1500	23.72	0	0.296	18.72	0	0.001	26.71
		500	6.533	0	0.016	6.532	0	0.020	0.01
	KSE09	1000	6.539	0	0.062	6.451	0	0.001	1.36
		1500	6.526	0	0.034	6.448	0	0.001	1.21
		500	1.904	1	44.73	1.803	1	0.100	5.60
	L99b	1000	1.791	1	198.64	1.789	1	0.06	0.11
		1500	1.827	1	543.09	1.776	1	0.08	2.87
		500	0.498	0	0.030	0.493	0	0.001	1.01
	BBD12	1000	0.499	0	0.020	0.497	0	0.001	0.40
		1500	0.499	0	0.030	0.497	0	0.020	0.40

Table 6 .

 6 Recourse cost analysis of instance L99b.

				Estimated recourse cost
		λ	R W 1	R W 2	R W 3	R total	η(%)	̺(%)	CPU time
		5000 7500	0 0	0.00006 0.00008	0.25006 0.23578	0.25012 0.23586	-85.75	5	24 54
	x	10000 12500	0 0	0.00005 0.00010	0.24539 0.24623	0.24544 0.24633	-85.07	5	94 144
		15000	0	0.00007	0.24874	0.24881	-85.80	5	204
		5000	0.00948	0	0.02617	0.03565			25
	x ′	7500 10000	0.00809 0.00883	0 0	0.02712 0.02601	0.03521 0.03484	-85.24	5	53 94
		12500 15000	0.00873 0.00847	0 0	0.02763 0.02727	0.03636 0.03574	-85.64	5	145 204

developed model should also be extended to more sophisticated line configurations, for example, with parallel workstations and consider multi-product types with different end-of-life states.

Disclosure statement

No potential conflict of interest was reported by the authors.