D Kumar
email: dkumar@laas.fr

J Philippe
email: jphilipp@laas.fr

S Charlot
email: scharlot@laas.fr

S Assié
email: sassie@laas.fr

A Lecestre
email: alecestre@laas.fr

A Ferrand
email: andre.ferrand@insa-toulouse.fr

H Aubert
email: aubert@laas.fr

P Pons
email: ppons@laas.fr

Polyethylene Outgassing Study for MEMS Nuclear Radiation Sensor

Keywords: polyethylene, nuclear radiation sensor, MEMS dosimeters, hydrogen outgassing

In this paper, the out-gassing from PolyEthylene (PE) film is studied under nuclear irradiation for the design and fabrication of passive MEMS dosimeters. The fabrication of sensor needs high temperature process that leads to specific constraints on the PE film. Radiation chemical yield of hydrogen production (G H2 ) from PE under gamma irradiation is verified by mass spectroscopy after temperature annealing in vacuum up to 400°C. Prototypes are fabricated to validate the fabrication of the sensor and then irradiated with high energy gamma radiation (with dose of 20 kGy). Measurements of membrane deflections after irradiation validate the G H2 factor, showing low pollution level of PE during sensor fabrication.

Introduction

Wireless passive (battery-less) and chip-less (without electronic circuits) sensors present interesting solution for harsh environment applications. In the case of nuclear radiation measurement, such solutions allow real time dose monitoring by remotely interrogating nuclear radiation sensors from radar reader [START_REF] Arenas | Wireless and Passive Nuclear Radiation Sensors[END_REF].

The operating principle of the sensor is the following: when PolyEthylene (PE) placed inside a hermetic sealed cavity is subjected to nuclear irradiation, an outgassing (mainly of H 2 molecules) occurs and generates a pressure on a highly resistive silicon membrane. This lead to the membrane deflection, which modify the resonant frequency of a planar microwave resonator placed just below the membrane. From the measurement of the frequency shift, the nuclear radiation dose is derived.

In our previous design (see Figure 1), a double Si/Glass stack was used for the sensor fabrication and leads to a complex technological process and non-sealed cavity for the resonator. The PE was obtained from peeling small ball and the Si/Glass2 bonding was performed under N 2 using a noncollective process Moreover, the sensor sensitivity varied significantly from one sensor to another and was about 20 times higher than the simulation results [START_REF] Augustyniak | Micromechanical high-doses radiation sensor with bossed membrane and interferometry optical detection[END_REF]. The possible reason could be the pollution of PE during sensor fabrication leading to the undesirable out-gassing after irradiation. Consequently, in order to control these additional out-gassing of PE, we report here a simple, more reliable and collective technological fabrication process for passive MEMS dosimeters with a new design (see Figure 2). The process uses a simple Si/Glass stack where the calibrated 200µm thick High Density PolyEthylene (HDPE) film and the resonator are sealed inside the same cavity under vacuum. Moreover, the millimeter-wave excitation of the resonator is performed by using a small radiating aperture [START_REF] Arenas | Wireless and Passive Nuclear Radiation Sensors[END_REF].

Materials and Methods

Prototypes were designed to measure radiation chemical yield of hydrogen production (G H2 ) from sensing devices fabricated with technological process similar to that used for the final sensors (see Figure 3). The HDPE cavity (5mm diameter and 290µm depth) is etched inside a low resistivity 400µm thick SOI wafer. The H 2 out-gassing will provide the deflection of the small silicon membrane (550µm diameter and 5µm thick) that will be monitored with optical profiler.

For G H2 = 4.10 -7 mol/J, the pressure P g generated by the PE is given by equation 1 [START_REF] Augustyniak | Micromechanical high-doses radiation sensor with bossed membrane and interferometry optical detection[END_REF], where V PE and V g denote respectively the PE volume and the volume of gas expansion:

𝑷𝑷 𝒈𝒈 [𝑩𝑩𝑩𝑩𝑩𝑩] = 𝟗𝟗. 𝟓𝟓 𝟏𝟏𝟏𝟏 -𝟑𝟑 𝑽𝑽 𝑷𝑷𝑷𝑷 �𝒎𝒎𝒎𝒎 𝟑𝟑 � 𝑽𝑽 𝒈𝒈 [𝒎𝒎𝒎𝒎 𝟑𝟑 ] ≅ 𝟏𝟏𝟏𝟏 -𝟐𝟐 𝑽𝑽 𝑷𝑷𝑷𝑷 �𝒎𝒎𝒎𝒎 𝟑𝟑 � 𝑽𝑽 𝒈𝒈 �𝒎𝒎𝒎𝒎 𝟑𝟑 � /𝒌𝒌𝒌𝒌𝒌𝒌 (1)
With dimensions given in the Figure 3, P g ≅ 0.045 bar/10kGy which leads a membrane deflection sensitivity of 0.23µm/10kGy.

Results and discussions

PE properties

One of the major issue during the technological fabrication of sensors is the use of high temperature silicon/glass anodic bonding that may degrade the H 2 out-gassing capabilities of HDPE.

The thermo-gravimetric analysis performed under Helium (He) atmosphere shows that the degradation of HDPE occurs only above 400°C (see Figure 4). To provide a security margin, the silicon/glass temperature bonding is then performed at 350°C. HDPE samples are annealed in vacuum at 350°C and then placed in a sealed glass tube under He atmosphere. After irradiation to gamma rays ( 137 Cs 661.7 keV, dose rate 0.3kGy/h, dose 24kGy), the gas composition analysis is finally performed using a quantitative gas mass spectrometer. It

shows that H 2 production is not affected by high temperature process (G H2 ≅ 4 10 -7 mol / J) [START_REF] Seguchi | Hydrogen yield from polyethylene by radiolysis and the application to dosimeter[END_REF].

Prototypes fabrication and characterization before irradiation

For silicon wafer fabrication, the two main steps were: (1) the 5µm thick membrane release and (2) the HDPE cavity fabrication with deep reactive ionic etching using 40µm thick photoresist mask. For the glass wafer, the 10µm depth cavity were obtained from etching in hydrofluorydric bath using Chromium/Gold mask. After 2 minutes annealing at 200°C under vacuum, 200µm HDPE film was cut using a stamping pliers and the 3.5mm diameter samples were placed inside the silicon cavity. Silicon and glass wafer were aligned in AML bonder, heated at 350°C under vacuum during 30min and then bonded at 600V. Prototypes were fabricated using this technological process (see Figure 5). Before silicon/glass bonding, the membrane deflection is small (≅ 0.2µm) and after bonding, the deflection is found to range from 5.1µm to 5.4µm for atmospheric pressure applied on the top of the membrane (see Figure 6). This experimental result is close to the expected value with vacuum inside the cavity (simulated sensitivity was of 5µm/bar). After two weeks, the membrane deflection remains unchanged. As displayed in Figure 6, the membrane sensitivity to pressure variation was also calibrated under pressure between 0.2 bar and 1 bar. For pressure between 0.8 bar and 1 bar, the sensitivity S is around 4 µm/bar.

Prototypes characterization after irradiation

The prototypes are irradiated to gamma rays ( 137 Cs 661.7 keV, dose rate 0.3kGy/h, dose 20kGy). The membrane sensitivity to pressure was also calibrated 3 weeks after irradiation and we checked that no significant change occurred after irradiation (see Figure 6). For atmospheric pressure, the membrane deflection is reduced by 0.4µm. The pressure P g generated inside the HDPE cavity after this irradiation can be then derived from equation (2) :

𝑷𝑷 𝒈𝒈 = ∆𝑷𝑷 𝒐𝒐 -(∆𝑫𝑫 𝑺𝑺 ⁄ ) (2)
where ∆P o denote the atmospheric pressure variation between two measurements of membrane deflection, ∆D is the variation of the membrane deflection related to irradiation dose and S designates the membrane sensitivity around 1 bar derived from Figure 6. The pressure P g is then found to be of 0.05 bar/10kGy and is close to the expected value (simulation results predict a pressure of 0.045 bar/10kGy). After two weeks, the membrane deflection remains unchanged and consequently, we concluded that a good H 2 tightness was achieved.

Conclusions

Prototypes based on silicon membrane were designed and fabricated to validate the fabrication process of a passive MEMS dosimeter. After irradiation to gamma rays (20kGy), membrane deflection analysis shows that the radiation chemical yield of hydrogen production (G H2 ) from HDPE is close to simulation results. So few undesirable out-gassing from polyethylene occurs during prototype irradiation, validating the fabrication process with low level of HDPE pollution.

Figure 1 .Figure 2 .

 12 Figure 1. Previous dosimeter design Figure 2. New dosimeter design

Figure 3 .

 3 Figure 3. Cross section of the prototype for the measurement of G H2

Figure 4 .

 4 Figure 4. Evolution of PE mass versus temperature

Figure 5 .Figure 6 .

 56 Figure 5. Bottom views of fabricated prototypes

Acknowledgments: This work was partially funded by CARANUC project (Occitanie Region, France and Europe) and EDF (Electricité de France). This work was partly supported by LAAS-CNRS micro and nanotechnologies platform, members of the French RENATECH network.

Conflicts of Interest:

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.


J Veyrunes
email: veyrunes@etud.insa-toulouse.fr

J Riondet
email: jriondet@laas.fr

A Ferrand
email: andre.ferrand@insa-toulouse.fr

M Lavayssière
email: maylis.lavayssiere@cea.fr

A Lefrançois

J Luc

H Aubert
email: aubert@laas.fr

P Pons
email: ppons@laas.fr

Dynamic Mechanical Simulation of Miniature Silicon Membrane during Air Blast for Pressure Measurement †

Keywords: air blast, pressure sensor, dynamic behavior, silicon membrane

The development of new ultra-fast sensors for pressure air blast monitoring requires taking into account the very short rise time of pressure occurring during explosion. Simulations show here that the dynamic mechanical behavior of membrane-based sensors depends significantly on this rise time when the fundamental mechanical resonant frequency of the membrane is higher than 10 MHz.

Introduction

The real time and dynamic measurements of pressure during air blasts is very challenging due to the abrupt variation of pressure from the atmospheric pressure to the so-called overpressure peak Pmax (between few bars and several ten of bars) with very short rise time tm (< 100ns) (Figure 1). For accurate measurement of Pmax, sensors with high fundamental mechanical resonant frequency Fo are then required [START_REF] Riondet | Design of air blast pressure sensors based on miniature silicon membrane and piezoresistive gauges[END_REF]. When simulating a pressure sensor for which Fo << 1/tm, the rise-time is assumed to be close to zero. However, we show here that this assumption is no more valid when performing the mechanical simulation of a sensor for which Fo > 10 MHz.

Sensor description and simulation conditions

As shown in Figure 2, the proposed sensor is based on a miniature rectangular silicon membrane (5µm*30µm*90µm) with four piezoresistive gauges located at its center [START_REF] Riondet | Design of air blast pressure sensors based on miniature silicon membrane and piezoresistive gauges[END_REF]. Simulations were performed using Abaqus software [2] for real membrane clamping conditions. Considering the gauges dimensions, the stresses in the gauge areas are close to ones calculated at the membrane center (error lower than 7% on the sensor response). Dynamic behavior of the sensor is then modelled by the differential stress ∆σ given by equation 1 at the center of the membrane, and normalized by the static pressure :

Δ𝜎𝜎 = 𝜎𝜎 𝑙𝑙 -𝜎𝜎 𝑡𝑡 (1)
where 𝜎𝜎 𝑙𝑙 (resp., 𝜎𝜎 𝑡𝑡 ) is the stress applied to the gauge parallel (resp., to the gauge perpendicular) to the current in the gauge. The fundamental mechanical resonant frequency Fo, obtained from harmonic module of Abaqus software, is of 32.7 MHz and consequently 1/Fo is close to 30ns. Due to the short reaction time tr of the sensor (< few µs, see Figure 5), we assume here that the pressure profile (shown on Figure 1) can be modelled by the Heaviside step function (the decrease of pressure is lower than 2% after 1µs).

Moreover, assuming that the acoustic damping is predominant [START_REF] Prak | Q-factor and frequency shift of resonating silicon diaphragms in air[END_REF], the quality factor Q of the membrane is inversely related to the equivalent pressure Pe applied both sides of the membrane as shown in equation 2.

𝑄𝑄 ≅ 95 P 𝑒𝑒 (𝑏𝑏𝑏𝑏𝑏𝑏)

(2) Abaqus simulations indicates that the pressure Pe is the average pressure applied between the two membrane sides. Table 1 reports the Q factor for typical applied pressure on one side of the membrane while the other side is in vacuum. We observe that the Q factor decreases rapidly and is lower than 20 for the applied pressure greater than 10 bars. For a shock wavefront normally incident upon the membrane surface, the dynamic mechanical response to a linear variation of pressure was obtained from the Second Order Transfer Function (SOTF) where the rise time differs from zero (Equation 3). This model is consistent with Abaqus simulations results, as it can be noticed from Figure 3. We can observe that the difference is lower than 4% after the rise time. Consequently, the effects on sensor reaction time can be neglected.

N-type silicon membrane P-type gauges SOI substrate SiO2

(100

) <110> R 1 Ground V a = 10V y x σ t = σ xx R 2 R 3 R 4 V out <1-10> σ l = σ yy σ t = σ yy σ l = σ xx (a) (b) 30µm 90µm 𝑅𝑅 𝑛𝑛 (𝑡𝑡) = 1 𝑡𝑡 𝑚𝑚 • � 𝑡𝑡 -2𝐵𝐵 + 𝑒𝑒 -𝐴𝐴𝑡𝑡 • �2𝐵𝐵 • 𝑐𝑐𝑐𝑐𝑐𝑐(Ω • 𝑡𝑡) + 2𝜉𝜉 2 -1 Ω • 𝑐𝑐𝑠𝑠𝑠𝑠(Ω • 𝑡𝑡)� -𝑈𝑈(𝑡𝑡 -𝑡𝑡 𝑚𝑚 ) �𝑡𝑡 -𝑡𝑡 𝑚𝑚 -2𝐵𝐵 • 𝑒𝑒 -𝐴𝐴(𝑡𝑡-𝑡𝑡 𝑚𝑚 ) • �2𝐵𝐵 • 𝑐𝑐𝑐𝑐𝑐𝑐�Ω • (𝑡𝑡 -𝑡𝑡 𝑚𝑚 )� + 2𝜉𝜉 2 -1 Ω • 𝑐𝑐𝑠𝑠𝑠𝑠�Ω • (𝑡𝑡 -𝑡𝑡 𝑚𝑚 )��� � (3)
where 𝜉𝜉 is the damping factor, 𝐴𝐴 = 𝜉𝜉𝜔𝜔 0 , 𝐵𝐵 = 𝜉𝜉 𝜔𝜔 0

, Ω = 𝜔𝜔 0 �1 -𝜉𝜉 2 , 𝜔𝜔 0 = 1 𝐹𝐹 𝑜𝑜 ⁄ and 𝑈𝑈(𝑡𝑡) = 0 if 𝑡𝑡 < 0, 𝑈𝑈(𝑡𝑡) = 𝑡𝑡 𝑡𝑡 𝑚𝑚 if 0 < 𝑡𝑡 < 𝑡𝑡 𝑚𝑚 , 𝑈𝑈(𝑡𝑡) = 1 if 𝑡𝑡 > 𝑡𝑡 𝑚𝑚

Results

Figure 4 displays an example of the dynamic mechanical response of the membrane for rise time tm up to 100ns and for a Q factor of 19. The response is normalized by the static value. The reaction time tr of the sensor (±5% of the static response) is extracted from these figures for Q factor between 5 and 100 (Figure 5). The following observations can be made:

i. for tm < 1/ Fo, we retrieve the classical damped oscillation for which tr is mainly driven by the Q factor (tr ≅ Q / Fo); ii. as tm increases, the reaction time tr decreases and, due to the decreasing of the fundamental mode amplitude (Figure 6), we obtain tr ≅ tm when tm = 1/ Fo; iii. when tm = n/ Fo where n is a natural number, no oscillation occurs and tr is close to tm; iv. when tm > 1/ Fo, small oscillations appear during the rise time.

The reduction of the reaction time related to the increase of the rise time is given in Figure 7 (reaction time for tm=0 is taken as reference). We can observe that for pressure rise time tm greater than 20ns, the sensor reaction time may be underestimated almost by 20% compared with the zero pressure risetime (tm=0) assumption. This error may reach 100% when tm is close to n/ Fo.

Conclusions

The dynamic mechanical behavior of membrane-based piezoresitive sensors was modelled by a second order transfer function with non-zero pressure rise time. For miniaturized silicon membrane (5µm x 30µm x 90µm) the fundamental resonant frequency Fo is of 33 MHz. The results show that pressure rise time tm larger than 1/(3 Fo) plays a crucial role for the accurate estimation of the sensor reaction time. Consequently, when designing ultra-fast sensors for pressure air blast monitoring, the increase of frequency Fo for a given rise time is expected to provide a significant reduction of the sensor reaction time.

Figure 1 .

 1 Figure 1. Typical dynamic pressure variation during an air blast experiment : (a) Illustration of the different phases and (b) example of pressure measurement at 1m from 1kgTNT by using a Tourmaline piezoelectric sensor.

Figure 2 .

 2 Figure 2. a) Cross sectional view of the sensor; (b) Wheatstone bridge on the rectangular membrane.

Figure 3 .

 3 Figure 3. Comparison of SOTF model and Abaqus simulation results for the dynamic response

Figure 4 .Figure 5 .Figure 6 .Figure 7 .

 4567 Figure 4. Dynamic sensor response for different rise times

Table 1 .

 1 Q factor versus absolute pressure from equation 2.

	Pressure (bar)	2	5	10	20	30	40
	Q	95	38	19	9.5	6.3	4.8

Proceedings 2018, 2, x

Acknowledgments: This work was partially funded by Occitanie Region (France) through the COCNANO project.

Conflicts of Interest:

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.


Mohand Lounes Bentaha
email: mdlounes.bentaha@gmail.com

Olga Battaïa
email: battaia@emse.fr

Alexandre Dolgui
email: dolgui@emse.fr

S Jack Hu

Innovative Applications of O.R

Keywords:

A problem of profit oriented disassembly line design and balancing with possible partial disassembly and presence of hazardous parts is studied. The objective is to design a production line providing a maximal revenue with balanced workload. Task times are assumed to be random variables with known normal probability distributions. The cycle time constraints are to be jointly satisfied with at least a predetermined probability level. An AND/OR graph is used to model the precedence relationships among tasks. Several lower and upper-bounding schemes are developed using second order cone programming and convex piecewise linear approximation. To show the relevance and applicability of the proposed approach, a set of instances from the literature are solved to optimality.

Introduction

Disassembly process plays a crucial role in material and product recovery. It is a required condition for an efficient treatment of endof-life (EOL) products [START_REF] Ilgin | Environmentally conscious manufacturing and product recovery (ECMPRO): A review of the state of the art[END_REF], 2012). The objective of disassembly is to separate EOL products subassemblies and components for recycling, remanufacturing and reuse. To carry out disassembly operations with higher productivity rate, disassembly lines are used [START_REF] Güngör | Disassembly line in product recovery[END_REF].

From practical point of view, disassembly process is more complex than assembly. In fact, in a disassembly environment, a product is broken down into several components and subassemblies whose quality, quantity and reliability cannot be controlled as in an assembly environment. The structure and quality of EOL products are strongly uncertain and even the number of components in such products can not be predicted. Moreover, an EOL product may contain certain hazardous material which necessitates special handling at a workstation of a disassembly line. Due to technical or economic restrictions such as irreversible connections of components of a product and low revenue obtained from retrieved parts, disassembly is usually a partial process [START_REF] Lambert | Determining optimum disassembly sequences in electronic equipment[END_REF].

Because of the peculiarities given above, the design and balancing of disassembly lines (known as DLBP: disassembly line balancing problem), is a hard optimization problem and needs adapted solution methods. A disassembly line consists of an ordered sequence of workstations connected by a material handling system which is used to transport work-pieces from one workstation to another. As aforementioned, certain parts or subassemblies may be hazardous and require a particular treatment incurring a supplementary cost.

The studied optimization problem consists in assigning a given set of disassembly tasks (of an EOL product) to an ordered sequence of workstations, while respecting precedence and cycle time constraints. Cycle time constraints are to be jointly satisfied with at least a certain probability level (1 -α) fixed by the decision maker. Task times are assumed to be independent random variables with known normal probability distributions. The main objective is to maximize the profit produced by the line by optimizing the number of needed workstations of the line and the depth of the disassembly process. Subsequently, the idle times at workstations should be as smooth as possible.

Although the main purpose of this paper is to study stochastic DLBP, it is also shown that the obtained results remain valid for stochastic assembly line design and balancing problem (ALBP).

The paper is organized as follows: Section 2 provides an overview of the relevant literature on disassembly and assembly line design and balancing under uncertainty. A formal description of the studied problem is given in Section 3. Section 4 presents the developed solution approach. Numerical experiments and optimization results are presented in Section 5. Section 6 concludes the paper with future research directions.

Literature review

In this section, papers dealing with line design and balancing under uncertainty of the task processing times for both disassembly and assembly are discussed. In addition, problems that have studied the case of disassembly/assembly processing alternatives are reviewed.

Disassembly line design and balancing

Only limited studies in the literature have taken into account the task processing times variability that characterizes the disassembly context in DLBP. A fuzzy colored Petri net model with a heuristic solution method was proposed in [START_REF] Turowski | Disassembly line design with uncertainty[END_REF] to study the human factors that cause uncertainty of task times. A collaborative ant colony algorithm for stochastic mixed-model U-shaped DLBP 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 to be stochastic with known normal probability distributions. A binary bi-objective non linear program was developed in Aydemir- [START_REF] Aydemir-Karadag | Multi-objective optimization of stochastic disassembly line balancing with station paralleling[END_REF] for DLBP under uncertainty of the task times. Task times were assumed to be independent random variables with known normal probability distributions.

Several mathematical models have also been developed for DLBP under uncertainty of task processing times. In Bentaha, Battaïa, and Dolgui (2014a), uncertainty was modeled using the notion of recourse cost and a sample average approximation method was developed to solve the studied optimization problem. In [START_REF] Bentaha | Dealing with uncertainty in disassembly line design[END_REF], uncertainty was modeled using workstation expectation times. In [START_REF] Bentaha | Disassembly line balancing and sequencing under uncertainty[END_REF], a stochastic program was developed for the joint problem of disassembly line balancing and sequencing under uncertainty. In Bentaha, Battaïa, and Dolgui (2014c), a Lagrangian relaxation was proposed to maximize the disassembly line profit under task times variability where workstation expectation times are considered.

To model the possible disassembly process alternatives and precedence relationships among tasks, some of the existing papers have used directed graphs called AND/OR graphs. There are two types of such graphs: AND/OR graphs constituted of tasks and AND/OR graphs constituted of tasks and subassemblies. The first type is considered in [START_REF] Altekin | Task-failure-driven rebalancing of disassembly lines[END_REF][START_REF] Altekin | Profit-oriented disassembly-line balancing[END_REF][START_REF] Güngör | A solution approach to the disassembly line balancing problem in the presence of task failures[END_REF][START_REF] Lambert | Determining optimum disassembly sequences in electronic equipment[END_REF], the second in (Bentaha et al., 2014a[START_REF] Bentaha | Disassembly line balancing and sequencing under uncertainty[END_REF](Bentaha et al., , 2014c[START_REF] Bentaha | Dealing with uncertainty in disassembly line design[END_REF][START_REF] Koc | Two exact formulations for disassembly line balancing problems with task precedence diagram construction using an AND/OR graph[END_REF][START_REF] Lambert | Linear programming in disassembly/clustering sequence generation[END_REF]. The latter which includes an explicit representation of subassemblies as well as tasks, is used in this paper. It is explained in detail in Section 3.

Assembly line design and balancing

Even if uncertainty level is lower in assembly, however, different sources from the assembly environment may cause the task time variations, as for example, non qualified operators, machine failures, complex assembly tasks, etc. To deal with this uncertainty, the following models were proposed in the literature. Task times were assumed to be random variables with either known continuous probability distributions [START_REF] Zhao | Modeling and analysis of a mixed-model assembly line with stochastic operation times[END_REF], or known or unknown symmetric probability distributions [START_REF] Betts | Identifying multiple solutions for assembly line balancing having stochastic task times[END_REF][START_REF] Raouf | A new method for assembly line balancing having stochastic work elements[END_REF], or known independent normal probability distributions. This third case has received quite some attention: earlier papers have focused on optimizing straight assembly lines where heuristic [START_REF] Carter | A cost-effective approach to stochastic line balancing with off-line repairs[END_REF][START_REF] Chakravarty | A cost minimization procedure for mixed model production lines with normally distributed task times[END_REF][START_REF] Fazlollahtabar | A heuristic methodology for assembly line balancing considering stochastic time and validity testing[END_REF][START_REF] Kao | Computational experience with a stochastic assembly line balancing algorithm[END_REF][START_REF] Lyu | A single-run optimization algorithm for stochastic assembly line balancing problems[END_REF][START_REF] Shin | An efficient heuristic for solving stochastic assembly line balancing problems[END_REF][START_REF] Silverman | A cost-based methodology for stochastic line balancing with intermittent line stoppages[END_REF], metaheuristic [START_REF] Cakir | Multi-objective optimization of a stochastic assembly line balancing: A hybrid simulated annealing algorithm[END_REF][START_REF] Erel | Stochastic assembly line balancing using beam search[END_REF] and exact solution methods [START_REF] Henig | Extensions of the dynamic programming method in the deterministic and stochastic assembly-line balancing problems[END_REF][START_REF] Kao | A preference order dynamic program for stochastic assembly line balancing[END_REF][START_REF] Sarin | A methodology for solving single-model , stochastic assembly line balancing problem[END_REF] were proposed. The case of ALBP with station paralleling was studied in [START_REF] Mcmullen | A heuristic for solving mixed-model line balancing problems with stochastic task durations and parallel stations[END_REF]. Optimization of U-lines was investigated in [START_REF] Bagher | Balancing of stochastic U-type assembly lines: an imperialist competitive algorithm[END_REF][START_REF] Baykasoglu | Stochastic U-line balancing using genetic algorithms[END_REF]Chiang & Urban, 2006;[START_REF] Guerriero | The stochastic U-line balancing problem[END_REF][START_REF] Özcan | A genetic algorithm for the stochastic mixedmodel U-line balancing and sequencing problem[END_REF]. Two heuristic approaches to the assembly line re-balancing problem were developed in [START_REF] Gamberini | A multiple single-pass heuristic algorithm solving the stochastic assembly line rebalancing problem[END_REF][START_REF] Gamberini | A new multi-objective heuristic algorithm for solving the stochastic assembly line re-balancing problem[END_REF]. In [START_REF] Liu | A bidirectional heuristic for stochastic assembly line balancing type II problem[END_REF], the authors studied the problem of minimizing the cycle time of the line to be designed.

Robust balancing of assembly lines with interval task times and stability analysis of optimal solutions for ALBP have been proposed, respectively, in (Gurevsky, Hazır, Battaïa, & Dolgui, 2013b;[START_REF] Hazır | A decomposition based solution algorithm for U-type assembly line balancing with interval data[END_REF] and (Gurevsky, Battaïa, & Dolgui, 2013a;[START_REF] Sotskov | Stability analysis of an optimal balance for an assembly line with fixed cycle time[END_REF]. Robust balancing of an assembly line with uncertain demand has been presented in [START_REF] Chica | A robustness information and visualization model for time and space assembly line balancing under uncertain demand[END_REF]. For cycle time minimization, two robust models and exact solution method for ALBP with interval uncertainty for task times have been proposed in (Hazır & Dolgui, 2013).

Particularly, for the case of task times following known normal probability distributions, exact and heuristic approaches were designed to solve integer linear programs with disjoint probabilistic constraints, for straight and U-lines [START_REF] Agpak | A chance-constrained approach to stochastic line balancing problem[END_REF]Urban & Chiang, 2006) and two-sided lines (Özcan, 2010).

Modeling of process alternatives and precedence relationships among tasks for assembly line balancing is undertaken in [START_REF] Capacho | The ASALB problem with processing alternatives involving different tasks: Definition, formalization and resolution[END_REF], 2008). The authors introduced and defined a new graph using the notion of Alternative Subgraphs. To solve this problem, an exact approach has been proposed in [START_REF] Scholl | Optimally solving the alternative subgraphs assembly line balancing problem[END_REF] and heuristic approaches in [START_REF] Capacho | An evaluation of constructive heuristic methods for solving the alternative subgraphs assembly line balancing problem[END_REF]. It should be noted that Alternative Subgraphs graph is exclusively constituted of tasks and does not represent the possible subassemblies as does AND/OR graph used in this study.

As it can be seen, joint satisfaction of cycle time constraints with a certain probability level has not been considered neither for DLBP nor for ALBP. The next section presents the developed formulation for the former problem that with some reduction can be also applied for latter problem.

Problem statement

The aim is to assign disassembly tasks from set I to an ordered sequence of workstations from set J, while satisfying precedence and cycle time constraints under uncertainty of the task processing times. The value of |J| represents the worst case for the number of workstations of the line. For a given problem instance, |J| corresponds to the number of tasks of the longest disassembly process alternative (longest in terms of number of tasks). The goal is to design a line providing the maximal profit and resulting in a number of stations m * |J|. Cycle time (C t ) constraints for all workstations have to be jointly respected with at least a probability level (1 -α) fixed by the decision maker; 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 following assumptions are used. A single type discarded product has to be partially (or completely) disassembled on a straight paced line. All received EOL items contain all initial parts with no addition or removing of components. Certain components of the EOL products are hazardous. A task can be performed at any workstation but cannot be split between two workstations. Task processing times are independent from the order in which the tasks are performed. Each component or subassembly has a certain non negative resale value but can be 0. A fixed cost per operating a time unit of an opened workstation and an additional fixed cost per operating a time unit for treating a hazardous part are given. (Bentaha et al., 2013a).

Task times are assumed to be mutually independent random variables with known normal probability distributions, i.e. t i ( ξ ) N (μ i , σ i ), i ∈ I. These random variables represent a random vector ξ = (t 1 , t2 , . . . , t|I| ) varying over a set ⊂ R |I| + given a probability space ( , F, P) introduced by ξ . The normal-distribution assumption is motivated by its several theoretical characteristics mainly, here, the central limit theorem. The results obtained for the lower and upper bounds can only be used for normal-distribution assumption of the workstation times. The time of a given workstation is the sum of processing times of the tasks assigned to it. This is particularly true for normal-distribution assumption of each task processing time. In the case of task times assumed to be random variables following a different probability distribution, the resulted time of each workstation is not guaranteed to be normally distributed for any number of tasks assigned to it. A disassembly task i ∈ H ⊂ I is called hazardous if its execution generates a hazardous subassembly or component. The additional time which could be generated by handling hazardous material is assumed to be taken into account while estimating processing times of the corresponding hazardous tasks.

All possible alternatives for the disassembly process and precedence relationships among tasks and subassemblies, for an EOL product, will be modeled by an AND/OR graph [START_REF] Bentaha | A decomposition method for stochastic partial disassembly line balancing with profit maximization[END_REF][START_REF] Koc | Two exact formulations for disassembly line balancing problems with task precedence diagram construction using an AND/OR graph[END_REF]. An example of such a graph for the disassembly of a piston and connecting rod (composed of 16 components) is given in Fig. 1. Each subassembly of an EOL product to be disassembled is represented by a node A k , k ∈ K in the AND/OR graph and each disassembly task by a node B i , i ∈ I. Set K contains the indices of all possible subassemblies that can be generated by the tasks from I.

Two types of arcs define the precedence relations between the subassemblies and tasks: AND-type and OR-type arcs. These arcs are represented in the graph with solid lines. For instance, if a task generates two subassemblies, or more, then it is related to these subassemblies by AND-type arcs. If, for a given subassembly, one or more tasks can be performed, but only one must be selected, this subassembly is related to these tasks by OR-type arcs. Table 1 summarizes all possible disassembly tasks that can be performed on the piston and connecting rod. For each task, the corresponding generated components and/or subassemblies are given. In order to handle the case of partial disassembly, where the product is not necessarily disassembled till 4:9,10,14,15 1;2;3 obtaining single components, all tasks are connected (with dashed lines) to the dummy task s introduced into the precedence graph in Fig. 1 as a sink node. Only a subset I * of set I is selected by the optimization procedure, therefore only the tasks of I * have to be assigned to the workstations of the line. The level of disassembly depends on the profit generated by the corresponding line. The recycling or reuse of certain parts or subassemblies generates profit while the cost of operating opened workstations and additional cost incurred by the treatment of hazardous parts are considered as negative revenue.

The optimization procedure consists of two phases: (1) A number of workstations of the line to be designed is obtained while maximizing line profit. Cycle time constraints are jointly satisfied with at least a probability (1 -α). In this first phase, using piecewise linear approximation techniques, several second-order cone programming formulations are developed providing lower and upper bounds on the optimal profit of the line; (2) for the number of workstations obtained in phase 1, an optimal balance is determined. The goal of this second step is to study the impact of balancing the workstations workload on the probability of cycle time constraint satisfaction resulted from the first step.

Phase 1: maximization of the line profit

Phase 1 consists of maximizing the profit of the line to be designed. The following Chance Constrained 0-1 Binary Program (CCBP) has been developed to formulate this problem.

K = {0, 1, . . . , k}, k ∈ N; A k : a subassembly, k ∈ K; B i : a disassembly task, i ∈ I;
s: the AND/OR graph's sink node; H: hazardous disassembly task index set;

L: parts index set: L = {1, 2, . . . , l}, l ∈ N * ;
L i : set of retrieved parts after execution of a disassembly task B i , i ∈ I; r : revenue generated by part , ∈ L; t i ( ξ ): the processing time of task i, i ∈ I, where t i ( ξ ) N (μ i , σ i ), i ∈ I; C t : cycle time, C t > 0; F c : fixed cost per time unit for workstations; H c : additional cost per time unit for stations handling hazardous parts;

P k : predecessors index set of A k , k ∈ K, i.e. P k = {i| B i precedes A k }; S k : successors index set 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.

x s j = 1, if dummy task s is assigned to workstation j; 0, otherwise.

h j =
1, if a hazardous task is assigned to workstation j; 0, otherwise.

Stochastic binary program

max i∈I j∈J ∈L i r • x i j -C t F c • j∈J j • x s j + H c • j∈J h j (CCBP) 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) P i∈I t i ( ξ ) • x i j C t , ∀ j ∈ J 1 -α (8) x s j , x i j , h j ∈ {0, 1}, ∀i ∈ I, ∀ j ∈ J (9)
The terms of the objective function represent, respectively, the earned profit of retrieved parts, the cost of operating workstations and the additional cost for handling hazardous parts. If the dummy task s is assigned to a workstation j, then j defines the number of stations. Constraint (1) imposes the selection of only one disassembly task (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 realizable for the initial EOL product. 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 ORsuccessor is selected. Constraint set (4) defines the precedence relationships among tasks and subassemblies. Constraint (5) imposes the assignment of the dummy task s to one station. Constraints ( 6) ensure that all selected 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) enforce the station operating time to remain within the cycle time, for all opened workstations, jointly with at least a predefined probability (1 -α). Finally, set (9) represents the possible values of the decision variables.

In program (CCBP), let r = 0, ∀ ∈ L i , ∀i ∈ I, H c = 0 and assume that the EOL product is completely disassembled and only one disassembly alternative exists. Then, the objective function of (CCBP)

becomes max{-C t • F c • j∈J j • x s j } = -min{C t • F c • j∈J j • x s j } and
there is no OR-precedence relations among tasks and subassemblies; the term j ∈ J j • x sj represents the number of opened workstations.

Theoretically, the resulted program defines a simple ALBP of type 1 with stochastic task times and Joint Probabilistic Constraints (JPC-ALBP). Indeed, the simple ALBP of type 1 is defined as the minimization of the workstations number under AND-precedence relationships among tasks and cycle time constraints. Along this paper, several lower and upper-bounding schemes will be developed for program (CCBP). From the results above, it follows that these same schemes will be applicable for the stochastic ALBP with joint probabilistic constraints.

Phase 2: line balancing

In phase 1, a number of workstations m * |J| (where x sm * = 1) and a subset of tasks I * ⊂ I are determined. Note that only the set or a subset of tasks of a disassembly alternative is assigned to the stations of the line. An example of such a disassembly alternative and a selected subset I * of disassembly tasks is illustrated in Fig. 2. This selected alternative (in bold) is represented by an AND-graph and

I * = {B 3 , B 16 , B 22 , B 24 }.
There are no more OR-relations since a decision of disassembling a product (completely in the case of Fig. 2) is made. For this reason, only direct precedence relationships among tasks are considered, i.e. generated subassemblies {A 5 , A 8 , A 10 } are simply deleted. The simple precedence graph in Fig. 3 is then created.

In phase 2, a balance measure is optimized. This measure seeks to assign the disassembly tasks I * to the m * stations enforcing similar idle time at each workstation. 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 differences among Expectations of Workloads of all workstations is Minimized with the following non linear program (MMEW): min max

∀ j, j ∈J * , j = j |E ξ (ST j ( ξ )) -E ξ (ST j ( ξ ))| (MMEW) s.t. j∈J * x i j = 1, ∀i ∈ I * (10)
* j • x i j j∈J * j • x i j , ∀i ∈ I * , ∀i ∈ Pred(i) (11) x i j ∈ {0, 1}, ∀i ∈ I * , ∀ j ∈ J * (12)
To linearize this model, 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 ; i.e. for each value of l is associated a couple (j, j ), j, j ∈ J * , j = j . Then,

ω l = E ξ (S l ( ξ )) = i∈I * μ i • x i j - i∈I * μ i • x i j , j, j ∈ J * , j = j , l ∈ B (13)
Thus, the program (MMEW') given below represents a linearized version of program (MMEW). min

Y MMEW s.t. j∈J * x i j = 1, ∀i ∈ I * (10) j∈J * j • x i j j∈J * j • x i j , ∀i ∈ I * , ∀i ∈ Pred(i) (11) -Y ω l Y, ∀l ∈ B (14) Y 0, x i j ∈ {0, 1}, ∀i ∈ I * , ∀ j ∈ J * (15)
Constraint set ( 14) can be replaced with:

e l + ω l = Y, ∀l ∈ B (16) 0 e l 2Y, ∀l ∈ B (17)
In this phase, the probability of jointly satisfying cycle time constraints may decrease or increase as a consequence of the possible reassignment of disassembly tasks I * . No reassignment of hazardous tasks will be considered. Hence, although the probability mentioned above may decrease or increase, the optimal value of (CCBP) remains unchanged unless this possible decrease or increase of probability will be seen, respectively, as additional cost or profit.

Let x 1 and x 2 be, respectively, optimal solutions of problems (CCBP) and (MMEW), and P 1 , P 2 the corresponding probabilities of jointly satisfying cycle time constraints. The solution x 2 is retained if

= P 1 -P 2 P 1 × 100 (18)
Otherwise, x 2 is rejected and x 1 is maintained; ϖ is a percentage fixed by the decision maker. In other words, x 2 is retained if the percentage of the decrease of the probability in phase 2 does not exceed a certain percentage ϖ.

The developed solution approach is detailed in the next section.

Solution method

In this section, the computational complexity of (CCBP) is first established in subsection 4.1. Then, in order to solve program (CCBP), a lower and several upper-bounding schemes are proposed in subsection 4.2. The development of these schemes is based on convex piecewise linear approximation [START_REF] Boyd | Convex Optimization[END_REF][START_REF] Magnani | Convex piecewise-linear fitting[END_REF] and second order cone programming [START_REF] Alizadeh | Second-order cone programming[END_REF][START_REF] Lobo | Applications of secondorder cone programming[END_REF][START_REF] Nemirovski | Convex approximations of chance constrained programs[END_REF][START_REF] Prékopa | Uniform quasi-concavity in probabilistic constrained stochastic programming[END_REF].

Computational complexity study

Consider the special case of problem (CCBP) defined as follows: task times are deterministic, there is no hazardous task, profits of recovered parts are set to 0, the EOL product is completely disassembled and only one disassembly alternative exists, i.e. precedence relationships among tasks are restricted to AND ones. An instance of such a special case is described with a finite set I of tasks, task times t i ∈ R + * , i ∈ I , a partial order ≺ on I , a number |J | ∈ Z + of workstations and a workstation capacity C t ∈ R + * . This defines the general instance of SALBP-1 [START_REF] Baybars | A survey of exact algorithms for the simple assembly line balancing problem[END_REF]. The problem SALBP-1 is known to be N P-hard [START_REF] Baybars | A survey of exact algorithms for the simple assembly line balancing problem[END_REF][START_REF] Becker | A survey on problems and methods in generalized assembly line balancing[END_REF]. As a consequence, (CCBP) is a N P-hard problem.

Approximation of (CCBP)

Let (1 -α) = ᾱ. Since disassembly task times are assumed mutually independent random variables with known normal probability distributions, then:

P i∈I t i ( ξ ) • x i j C t , ∀ j ∈ J ᾱ ⇐⇒ P i∈I t i ( ξ ) • x i j C t ᾱq j , ∀ j ∈ J, j∈J q j = 1 (19) and ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∀ j ∈ J : P i∈I t i ( ξ ) • x i j C t ᾱq j ⇐⇒ P i∈I t i ( ξ ) • x i j -i∈I μ i • x i j i∈I σ 2 i • x 2 i j C t -i∈I μ i • x i j i∈I σ 2 i • x 2 i j ᾱq j ⇐⇒ P Z j C t -i∈I μ i • x i j i∈I σ 2 i • x 2 i j ᾱq j , Z j N (0, 1) ⇐⇒ i∈I μ i • x i j + -1 ( ᾱq j ) • i∈I σ 2 i • x 2 i j C t (20)
Equivalence ( 19) is obtained in Cheng and Lisser (2012); -1 (t) is the inverse of the standard normal cumulative distribution function

(t) = 1 √ 2π t -∞ e -r 2 2 dr, t ∈ R.
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), ( 9) are satisfied}. Using the result above, problem (CCBP) is equivalent to the following problem (CCBP'):

max i∈I j∈J ∈L i r • x i j -C t F c • j∈J j • x s j + H c • j∈J h j (CCBP ) s.t. i∈I μ i • x i j + -1 ( ᾱq j ) • i∈I σ 2 i • x 2 i j C t , ∀ j ∈ J (21) j∈J q j = 1 (22)
x ∈ X, q j 0, ∀ j ∈ J (23) Inequalities ( 21) are second order convex cone constraints. Indeed,

let (v, w) ∈ R × R l-1 , then Q l = w v v w ( 24
)
defines a unit second-order convex cone of dimension l; • refers to the standard Euclidean norm. Since α < 50 percent (α is a risk and in general α 10 percent), then -1 ( ᾱq j ) > 0, thus (21) are second-order cone constraints of dimension l = |I| + 1 since:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ i∈I μ i • x i j + -1 ( ᾱq j ) • i∈I σ 2 i • x 2 i j C t , ∀ j ∈ J ⇐⇒ 1 2 • x j 1 -1 ( ᾱq j ) • (C t -μ T • x j ), ∀ j ∈ J ⇐⇒ 1 2 -μ T -1 ( ᾱq j ) x j + 0 Ct -1 ( ᾱq j ) ∈ Q |I|+1 , ∀ j ∈ J ( 25
)
where μ = (μ 1 , . . . , μ |I| ), x j = (x 1 j , . . . , x |I| j ) T , ∀ j ∈ J;

1 2 = diag(σ 1 , . . . , σ |I| ) is a diagonal matrix.
The Second Order Cone Mixed Integer Program (SOCMIP) given below represents an equivalent version of problem (CCBP'). max

i∈I j∈J ∈L i r • x i j -C t F c • j∈J j • x s j + H c • j∈J h j (SOCMIP) s.t. v j 1 -1 ( ᾱq j ) • (C t -μ T • x j ), ∀ j ∈ J (26) w i j σ i • x i j , ∀i ∈ I, ∀ j ∈ J (27) v j w j , ∀ j ∈ J (28) j∈J q j = 1 (22) x ∈ X, v j 0, w i j 0, ∀i ∈ I, ∀ j ∈ J (29)
where v j , w ij , ∀i ∈ I, ∀j ∈ J are intermediate variables, w j = (w 1 j , . . . , w |I| j ) T , ∀j ∈ J.

For the sake of compatibility with the literature [START_REF] Cheng | A second-order cone programming approach for linear programs with joint probabilistic constraints[END_REF][START_REF] Jagannathan | Chance-constrained programming with joint constraints[END_REF], to make it easier for readers and since max ϕ(x) =min(-ϕ(x)), where ϕ(x) is a function of decision variables x, in order to approximate (CCBP), minimization version of its objective:

min C t F c • j∈J j • x s j + H c • j∈J h j - i∈I j∈J ∈L i r • x i j (30)
will be used. The optimal value of the objective function of (CCBP) is then equal to the symmetric value of (30). Let this minimization version be called (M-CCBP).

Lower bounding scheme for (M-CCBP)

Problems of chance constrained programming are widely studied in the literature [START_REF] Blackmore | Chance-constrained optimal path planning with obstacles[END_REF][START_REF] Branda | Chance constrained problems: Penalty reformulation and performance of sample approximation technique[END_REF][START_REF] Depaolo | A heuristic algorithm for a chance constrained stochastic program[END_REF][START_REF] Liu | A dual-uncertainty-based chance-constrained model for municipal solid waste management[END_REF][START_REF] Mesfin | A chance constrained approach for a gas processing plant with uncertain feed conditions[END_REF][START_REF] Nemirovski | On safe tractable approximations of chance constraints[END_REF][START_REF] Poojari | Genetic algorithm based technique for solving chance constrained problems[END_REF][START_REF] Reich | A linear programming approach for linear programs with probabilistic constraints[END_REF][START_REF] Shen | Expectation and chance-constrained models and algorithms for insuring critical paths[END_REF][START_REF] Wang | A chance-constrained two-stage stochastic program for unit commitment with uncertain wind power output[END_REF][START_REF] Watanabe | A joint chance-constrained programming model with row dependence[END_REF][START_REF] Zhang | Chance constrained programming for optimal power flow under uncertainty[END_REF][START_REF] Zorgati | Optimizing financial and physical assets with chance-constrained programming in the electrical industry[END_REF]. A special case of linear programs with joint probabilistic constraints has been studied in [START_REF] Cheng | A second-order cone programming approach for linear programs with joint probabilistic constraints[END_REF]. Coefficients of the constraint matrix were assumed to be normally distributed and the vector rows independent. In our case, the coefficients are the task times and each row vector is composed of the processing times of the tasks assigned to a workstation. It is clear that task times of stations are independent of each other since all task times are already assumed to be mutually independent. Therefore, the main results concerning lower and upper bounds of the continuous program studied in [START_REF] Cheng | A second-order cone programming approach for linear programs with joint probabilistic constraints[END_REF] remain valid for the problem (M-CCBP) with 0-1 binary variables.

Piecewise tangent approximation of -1 ( ᾱq )

The piecewise tangent approximation of -1 ( ᾱq ) is needed in order to define a lower bound of program (CCBP).

The function -1 ( ᾱq ), q ∈]0, 1] is approximated using the firstorder Taylor series and input data (q j , -1 ( ᾱq j )), j = 1, . . . , m (q j is a tangent point). Assume that q 1 < q 2 < < q m . A piecewise tangent linear approximation of -1 ( ᾱq ) is given as follows:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ g(q) = max j=1,...,m {a j + b j q}, q ∈]0, 1] b j = ( -1 ) (1) ( ᾱq j ) • ᾱq j ln( ᾱ), j = 1, . . . , m a j = -1 ( ᾱq j ) -b j • q j , j = 1, . . . , m ( -1 ) (1) ( ᾱq j ) = 1 f ( -1 ( ᾱq j )) , j = 1, . . . , m (31)
where f is the standard normal probability density function.

Using approximation (31), program (SOCLB) below is an approximation of (M-CCBP).

min C t F c • j∈J j • x s j + H c • j∈J h j - i∈I j∈J ∈L i r • x i j (SOCLB) s.t. v j C t -μ T • x j , ∀ j ∈ J (32) w i j σ i • z i j , ∀i ∈ I, ∀ j ∈ J (33) v j w j , ∀ j ∈ J (28) z i j a k • x i j + b k • y i j , ∀i ∈ I, ∀ j ∈ J, k = 1, . . . , m (34) j∈J y i j = j∈J o i j , ∀i ∈ I (35) o i j x i j , ∀i ∈ I, ∀ j ∈ J ( 36
)
o i j q j , ∀i ∈ I, ∀ j ∈ J (37)

q j + x i j 1 + o i j , ∀i ∈ I, ∀ j ∈ J (38) j∈J q j = 1 (22)
x ∈ X, v j , q j , y i j , w i j , o i j , z i j 0, ∀i ∈ I, ∀ j ∈ J (39)

Moreover, the optimal value of (SOCLB) is a lower bound of (M-CCBP); it is, then, an upper bound of (CCBP). This lower bound is based on the lower bound defined in [START_REF] Cheng | A second-order cone programming approach for linear programs with joint probabilistic constraints[END_REF] for continuous decision variables.

Upper bounding schemes for (M-CCBP)

In this subsection, four approximations of (M-CCBP) will be defined, the value of each one represents an upper bound (a lower bound for (CCBP)). They are based on Bonferroni's inequality [START_REF] Galambos | Bonferroni inequalities[END_REF]), Jagannathan's approximation [START_REF] Jagannathan | Chance-constrained programming with joint constraints[END_REF], another approximation proposed in [START_REF] Cheng | A second-order cone programming approach for linear programs with joint probabilistic constraints[END_REF] and a variant of Bonferroni's approximation. A piecewise linear approximation of -1 ( ᾱq ) is given first.

Piecewise linear approximation of -1 ( ᾱq ) Since -1 ( ᾱq ), q ∈]0, 1] is convex, for the input data (q j , -1 ( ᾱq j )), j = 1, . . . , m (q j is an interpolation point), it follows a piecewise linear approximation g of -1 ( ᾱq ):

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ g(q) = max j=1,...,m-1 {a j + b j q}, q ∈ ]0, 1] a j = q j+1 -1 ( ᾱq j )-q j -1 ( ᾱq j+1 ) q j+1 -q j , j = 1, . . . , m -1 b j = -1 ( ᾱq j+1 )--1 ( ᾱq j ) q j+1 -q j , j = 1, . . . , m -1 q 1 < q 2 < • • • < q m , q j ∈ ]0, 1], j = 1, . . . , m (40)
The first upper bound approximation (SOCUB1) of (M-CCBP) is defined by replacing the (a k , b k ) values in (SOCLB) by their values defined in (40). This approximation, based on the one given in [START_REF] Cheng | A second-order cone programming approach for linear programs with joint probabilistic constraints[END_REF], defines an upper bound of (M-CCBP) if

γ = j∈J C t -μ T • x j 1 2 • x j 1 -α (41)
The second upper bound approximation (SOCUB2) is most common, it is based on Bonferroni's inequality as follows:

min C t F c • j∈J j • x s j + H c • j∈J h j - i∈I j∈J ∈L i r • x i j (SOCUB2) s.t. v j 1 -1 (1 -α j ) • (C t -μ T • x j ), ∀ j ∈ J (42) w i j σ i • x i j , ∀i ∈ I, ∀ j ∈ J (27) v j w j , ∀ j ∈ J (28)
x ∈ X, v j , w i j 0, ∀i ∈ I, ∀ j ∈ J (29)

where α j , ∀j ∈ J are parameters verifying j∈J α j = α.

The third upper bound (SOCUB3) of (M-CCBP) is defined by replacing (1 -α j ) values in (SOCUB2) by ( ᾱq j ) respectively; q j , ∀j ∈ J are parameters satisfying j∈J q j = 1. The fourth upper bound (SOCUB4) is based on Jagannathan's approximation, it is given below.

min C t F c • j∈J j • x s j + H c • j∈J h j - i∈I j∈J ∈L i r • x i j (SOCUB4) s.t. μ T • x j + o j • β j C t , ∀ j ∈ J (43) x T j • • x j 1 2 β j , ∀ j ∈ J (44) j∈J ln( (o j )) ln(1 -α) (45) x ∈ X, o j 0, ∀ j ∈ J (46)
where (t) = P(Z t ), Z N (0, 1) and β j , ∀j ∈ J are parameters which play a key role in the result of this approximation [START_REF] Jagannathan | Chance-constrained programming with joint constraints[END_REF];

= ( 1/2 ) 2 .
The function ln ( (o j )), j ∈ J has to be approximated using especially convex piecewise linear approximation in order to solve (SOCUB4). The function ln ( (o j )), j ∈ J is concave, hence, ln( 1 45) is equivalently written as:

(o j ) ) is convex. Constraint (
j∈J ln 1 (o j ) ln 1 1 -α (47)
Lemma. For the input data (q j , ln( 1 (q j ) )), j = 1, . . . , m, the function g defined below is a convex piecewise linear approximation of ln( 1 (q) ), q 0:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ g(q) = max j=1,...,m-1 {a j + b j q}, q 0 a j = q j+1 ln 1 (q j ) -q j ln 1 (q j+1 ) q j+1 -q j , j = 1, . . . , m -1 b j = ln 1 (q j+1 ) -ln 1 (q j ) q j+1 -q j , j = 1, . . . , m -1 q 1 < q 2 < • • • < q m , q j 0, j = 1, . . . , m (48)
Proof. Obvious since ln( 1 (q) ), q 0, is convex.

Constraint ( 45) is then replaced by:

j∈J w j ln 1 1 -α (49) w j a k + b k • o j , ∀ j ∈ J, k = 1, . . . , m -1 (50)
w j 0, ∀ j ∈ J (51)

Numerical experiments

Problems (SOCLB), (SOCUB1), (SOCUB2), (SOCUB3), (SOCUB4) and (MMEW') were implemented with MS Visual C++ 2008. ILOG CPLEX 12.4 was used to solve them on a PC with Pentium(R) Dual-Core CPU T4500, 2.30 Gigahertz and 3 Gigabyte RAM. These models have been applied to the problem instance illustrated in Fig. 1 and to five other instances available in the literature. The input data for each instance is given in Table 2. These instances are used as benchmarking problems and contain process alternatives for disassembly. The names of the problem instances are respectively composed of the first letters of authors' names and the year of publication as follows: BBD13 represents a compass (Bentaha et al., 2013c), BBD13a is a piston and connecting rod (Bentaha et al., 2013a), KSE09 is just a sample product created by the authors in [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 are respectively a radio set and a ball-point pen [START_REF] Lambert | Linear programming in disassembly/clustering sequence generation[END_REF] and MJKL11 is an automatic pencil [START_REF] Ma | Disassembly process planning algorithms for end of life product recovery and environmentally conscious disposal[END_REF]. 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 results of phase 1 obtained for the studied instances using the proposed upper bound and the first lower bound of (CCBP), i.e. the symmetric optimal values of (SOCLB) and (SOCUB1), respectively. These two bounds will be called main bounds. The number of points for convex piecewise linear approximation was 15, α = 5 percent, 25 percent of the disassembly tasks were assumed to be hazardous and the first point of input data for piecewise approximation was 0.0001; all sampled points were equidistant. The remaining parameters were randomly generated.

Columns 'UB' and 'LB' report the upper and lower bound values, respectively. Column 'Gap' reports the optimality gap UB-LB LB . Columns 'm * ', '|I * |', 'h-stations' and 'CPU time', in Table 3, report respectively the number of stations, the number of selected tasks, the number of hazardous workstations and the resolution time in seconds; γ = j∈J

( C t -μ T •x j 1 2 •x j
). In all tables, a dash '-' means that the corresponding value does not exist.

The obtained results show that, for each problem instance, the line to be designed is composed of at least two workstations. However, a single workstation can be preferred in some cases as shown with instance KSE09. All instances were solved to optimality in short time.

Table 4 presents lower bounds LB Bon , LB vBon and LB Jag of (CCBP) which represent, respectively, the symmetric optimal values of (SOCUB2), (SOCUB3) and (SOCUB4). The number of points in the convex piecewise linear approximation was fixed at 5. The values LB Bon and LB vBon are computed for α j = α |J| , ∀ j ∈ J and q j = 1 |J| , ∀ j ∈ J, respectively; LB Jag is computed for

β j = (x T j,Bon • • x j,Bon ) 1 2 , ∀ j ∈ J,
where x j, Bon , j ∈ J is an optimal solution of (SOCUB2). This choice of values of parameters β j is made in order to avoid infeasibility in solving program (SOCUB4). Table 4 shows that the returned objective values for all instances were optimal for program (CCBP) except the objective value for SOCUB3) provided optimal values for all processed instances. In order to analyze the impact of the number of points (or segments) of the piecewise linear functions used to approximate the non linear functions of problems (SOCLB) and (SOCUB1), the benchmark problems were solved for 5 points, 10 points and 20 points. The results have shown that even with 5 interpolation points, (SOCLB) and (SOCUB1) returned optimal objective values for (CCBP).

Table 5 reports the optimization results of the line balancing or idle time leveling phase. In this case, the problem (MMEW') was solved. A reassignment of tasks (of line profit maximization phase) in line balancing phase is represented with value 1 in column 'Reassignment', it is 0 otherwise. This value is 1 + if the reassignment leads to

P 2 > P 1 , where P 1 = j∈J ( C t -μ T •x 1 j 1 2 •x 1 j ) and P 2 = j∈J ( C t -μ T •x 2 j 1 2 •x 2 j ).
In such a case, x 2 , i.e. the solution of line profit maximization phase, is retained with no hesitation. The results have shown that, for in-stance BBD13a, the reassignment of tasks has lead to P 2 > P 1 . Therefore, the solution x 2 of line balancing phase should be retained as a final solution.

Optimization results for simple ALBP-1

As for DLBP, two phases were defined for the problem of simple assembly line design with joint probabilistic constraints (JPC-ALBP) introduced in Section 3. The first phase consists in minimizing the number of workstations, the second one in balancing the workload of the line derived from phase 1. The 25 used instances are available for users on the web site <http://alb.mansci.de>. The main parameters of the input data for each instance are given in Table 6; columns '|I|', '| J |' and 'C t ' refer, respectively, to the number of assembly tasks, the lower bound on the number of workstations and the cycle time.

The tests were realized on the same machine as for DLBP where the resolution time was limited by 1 hour.

Task times for all instances were assumed to be independent nor-

( C t -μ T •x j 1 2 •x j
). As shown in Table 7, the optimal lower and upper bounds were obtained for all instances, for most of them in less than 1 hour. However, for some instances like Arcus1 and Arcus2, an optimal solution could not be found in 1 hour. The problem JPC-ALBP was solved to optimality for most instances. In such a case, the lower bound | J * | and upper bound | J * | values were equal.

The optimization results of the idle time leveling phase are grouped in Table 8. Only instances for which an optimal solution of the problem JPC-ALBP is found in 1 hour were considered. For all processed instances, except for instance Mertens, probabilities of jointly satisfying cycle time constraints were increased. In such cases, the solution of the second phase is preferred to the solution of the first phase since it defines the same objective value and a better balance of the stations workload. The obtained results have shown that although P 1 > P 2 for instance Mertens, the relative decrease value ϱ of the probability of cycle time constraints satisfaction is smaller than the fixed level ϖ. Thus, the solution of the second phase of instance Mertens is retained as well.

The results above demonstrate the efficiency of the proposed mathematical models and solution procedures both for disassembly and assembly lines. Limitations of the developed models concern mainly the normal distribution assumption for the disassembly task durations and their mutual independence. The further work should relax these assumptions.

Conclusion

In this paper, a profit oriented disassembly line design and balancing problem was studied. The cases of partial disassembly and the presence of hazardous material were integrated and disassembly task times were assumed to be random variables with known normal probability distributions. To solve the addressed problem with an assessment of solution quality, a two phases approach was considered. In the first phase, a mixed integer program with joint probabilistic constraints, along with one lower bound and four upper bounds, were proposed for the line design.

The developed lower and upper bounding schemes were based on second order cone programming and convex piecewise linear approximation. Cycle time constraints were to be jointly respected with at least a probability level fixed by the decision maker. At this phase, the goal was to maximize the line profit. In the second phase, a 0-1 binary linear programming formulation was developed in order to balance the workload along workstations for the line derived from the first phase.

It is shown, in this study, that under some assumptions the studied problem can be reduced to the assembly line balancing problem. For this reason, the developed optimization procedure can be also applied for the case of assembly.

This paper is the first work dealing with such chance constrained disassembly and assembly problems proposing a lower and several upper bounding schemes and efficient exact models. The developed models were evaluated using a set of instances from the literature. For the disassembly line design phase, all the problem instances were solved to optimality in a short amount of time. As an example, the radio set instance consisting of 30 tasks is solved to optimality in less than 1 second. The results of the line balancing phase showed that it is possible to balance the workload of the designed disassembly line and in the same time improve (increase) the probability of jointly satisfying the cycle time constraints. This fact is occurred with the piston and connecting rod instance, where the probability has passed from 98.0 percent in the line design phase to 98.5 percent in the line balancing phase.

The optimization results in the case of assembly confirmed those of disassembly, mainly, in the line balancing phase. In fact, it is shown that the probability of jointly satisfying cycle time constraints can be increased. In particular, for the 83 tasks instance, this probability has passed from 96.3 percent to 98.2 percent. In the line design phase, the minimal number of workstations for most of the 25 processed instances were determined within reasonable computation time. For example, the 148 tasks instance with cycle time value of 7 is solved in less than 10 minutes.

The optimization results showed that the proposed lower and upper bounding schemes help to solve optimally the defined problems for both assembly and disassembly problem instances.

For future research, a cutting-plane algorithm will be investigated and compared to the default interior point algorithm of the solver Cplex. Other industrial and real life instances for the case of disassembly will also be considered.

Fig. 1 .

 1 Fig. 1. AND/OR graph of the piston and connecting rod (Bentaha et al., 2013a).

ParametersI:

  set of disassembly task indices: I = {1, 2, . . . , n}, n ∈ N * ; J: set of workstation indices: J = {1, 2, . . . , m}, m ∈ N * ; K: set of indices for the generated subassemblies:

Fig. 2 .

 2 Fig. 2. A selected disassembly alternative.

Fig. 3 .

 3 Fig. 3. Precedence AND-graph obtained from the selected disassembly alternative.

  j∈J

Table 1

 1 The piston and connecting rod associated disassembly tasks, components and/or subassemblies.

	Task	Subassembly	Components	Task	Subassembly	Components
			1;2;3	18	4:9	1;2;3
	6	1:4;9,11,14,16	5;6;7;8	19	1:4	5;6;7;8;9
	7	1:9	11;14;16	20	4:9	11;14;15
	8	4:9,11,14,16	10;12;13;15	21	9,10,14,15	4;5;6;7;8
	9	9:16	4;5;6;7;8	22	-	9;11;14;16
	10	4:9,10,14,15	11;12;13;16	23	-	4;5;6;7;8;9
	11	1:9	10;14;15	24	-	1;2;3;4
	12	1:4;9,10,14,15	5;6;7;8	25	-	9;10;14;15
	13

1 1:9,11,14,16 10;12;13;15 14 9,11,14,16 4;5;6;7;8 2 4:16 1;2;3 15 4:9 11;14;16 3 1:4;9:16 5;6;7;8 16 9,11,14,16 10;12;13;15 4 1:9,10,14,15 11;12;13;16 17 9,10,14,15 11;12;13;16 5 4:9,11,14,16

Table 2

 2 Problem instances.

		|I|	| K|	L	Arcs	AND-relations	|J|	Ct
						0	1	2
	MJKL11	37	22	33	76	4	27	6	10	40
	L99a	30	18	28	60	2	26	2	9	50
	BBD13a	25	11	27	49	4	18	3	4	120
	KSE09	23	13	20	47	4	14	5	6	20
	L99b	20	13	23	41	5	9	6	9	10
	B B D 1 3	1 0	5	1 2	1 8	3	6	1	3	0 . 6 1

Table 3

 3 Obtained results: main lower and upper bounds.

		m *	UB	|I * |	h-stations	time (second)	LB	|I * |	h-stations	Time (second)	γ percent	Gap percent
	MJKL11	3	199	6	1	2.277	199	6	1	1.654	99.5	0
	L99a	3	48	7	-	0.218	48	7	-	0.312	98.2	0
	BBD13a	2	13	3	-	0.062	13	3	-	0.094	98.7	0
	KSE09	1	590	4	-	0.001	590	4	-	0.016	97.3	0
	L99b	2	25	6	-	0.062	25	6	-	0.078	99.4	0
	BBD13	2	84	3	1	0.078	84	3	1	0.202	99.6	0

Table 4

 4 Obtained results: additional upper bounds.

		m *	LB Jag	Time (second)	LB Bon	Time (second)	LB vBon	Time (second)
	MJKL11	3	199	0.484	199	0.828	199	0.577
	L99a	3	48	0.047	48	0.079	48	0.077
	BBD13a	2	13	0.016	13	0.015	13	0.031
	KSE09	1	360	0.016	590	0.016	590	0.016
	L99b	2	25	0.031	25	0.015	25	0.015
	BBD13	2	84	0.062	84	0.062	84	0.062

Table 5

 5 Obtained results: idle time leveling.

		m *	Objective	Reassignment	Time (second)
	MJKL11	3	16.00	0	0.016
	L99a	3	7.00	0	0.031
	BBD13a	2	26.00	1 +	0.001
	L99b	2	0.70	0	0.016
	BBD13	2	0.07	0	0.001
	instance KSE09. Programs (SOCUB2) and (

  mal random variables with known means μ i , i ∈ I and standard deviations σ i , i ∈ I. The used means are the same as defined on the web site, while standard deviation values are equal to μ i /10, i ∈ I. The lower bound | J | is computed asi∈I μ i /C t .The first phase of JPC-ALBP is solved for |J| = |J|, if the defined optimization problem with | J | stations is not feasible, then, it is solved for |J| + 1 stations, etc, till obtaining a feasible problem. Note that with this procedure, if the problem is feasible for a number of stations |J| + a, then this number represents the minimal (optimal) number of workstations.

Table 6

 6 Problem instances for JPC-ALBP.

		|I|	|J |	Ct		|I|	|J |	Ct		|I|	|J |	Ct
	1. Mertens	7	2	18	9. Buxey	29	6	54	17. Wee-mag	75	32	47
	2. Bowman	8	4	20	10. Sawyer	30	5	75	18. Arcus1	83	7	10816
	3. Jaeschke	9	3	18	11. Lutz1	32	5	2828	19. Lutz2	89	24	21
	4. Jackson	11	3	21	12. Gunther	35	6	81	20. Lutz3	89	11	150
	5. Mansoor	11	2	94	13. Kilbridge	45	3	184	21. Mukherje	94	12	351
	6. Mitchell	21	3	39	14. Hahn	53	3	4676	22. Arcus2	111	9	17067
	7. Roszieg	25	4	32	15. Warnecke	58	14	111	23. Barthol2	148	25	170
	8. Heskiaoff	28	3	342	16. Tonge	70	7	527	24. Barthold	148	7	805
	25. Scholl	297	25	2787

Table 7

 7 Obtained results for JPC-ALBP: main lower and upper bounds.

			Time (second)	LB: | J * |	Time (second)	γ percent	Gap percent
	1. Mertens	2	0.92	2	0.17	99.97	0
	2. Bowman	6	0.42	5	0.41	96.12	20
	3. Jaeschke	3	0.09	3	0.20	98.14	0
	4. Jackson	3	0.20	3	0.22	99.97	0
	5. Mansoor	3	0.20	3	0.23	96.60	0
	6. Mitchell	3	0.17	3	0.36	96.67	0
	7. Roszieg	5	1.25	5	0.90	99.81	0
	8. Heskiaoff	4	1.78	4	2.82	99.20	0
	9. Buxey	8	3008.46	7	2.17	99.74	14.3
	10. Sawyer	5	1.37	5	1.03	97.79	0
	11. Lutz1	6	10.31	6	6.29	98.51	0
	12. Gunther	8	384.71	7	56.61	98.74	14.3
	13. Kilbridge	4	3.82	4	2.54	99.99	0
	14. Hahn	4	1.33	4	1.28	96.48	0
	15. Warnecke	15	>3600	15	>3600	96.13	0
	16. Tonge	8	>3600	8	2921.43	98.65	0
	17. Wee-mag	33	>3600	33	>3600	97.32	0
	18. Arcus1	8	>3600	8	454.21	96.27	0
	19. Lutz2	25	>3600	25	>3600	98.53	0
	20. Lutz3	13	3072.5	12	>3600	98.91	8.33
	21. Mukherje	13	>3600	13	>3600	96.76	0
	22. Arcus2	10	>3600	10	>3600	98.01	0
	23. Barthol2	26	>3600	26	>3600	97.87	0
	24. Barthold	8	294.51	8	265.61	96.57	0
	25. Scholl	26	>3600	26	>3600	96.28	0

UB : | J * |

Table 8

 8 Obtained results: JPC-ALBP idle time leveling.

		Objective	Reassignment	Time (second)
	1. Mertens	1	1 -	0.03
	3. Jaeschke	1	1 +	0.06
	4. Jackson	1	1 +	0.09
	5. Mansoor	1	1 +	0.05
	6. Mitchell	0	1 +	0.06
	7. Roszieg	4	1 +	0.49
	8. Heskiaoff	0	1 +	0.16
	10. Sawyer	1	1 +	4.05
	11. Lutz1	148	1 +	1.57
	13. Kilbridge	0	1 +	0.31
	14. Hahn	665	1 +	0.95
	16. Tonge	1	1 +	152.45
	18. Arcus1	290	1 +	>3600
	24. Barthold	-	-	>3600

  Table 7 reports the optimization results of the first phase for the main lower and upper bounds. Columns '| J * |', '| J * |', 'Gap' and 'Time (second)' refer, respectively, to the number of workstations of the upper bound, the number of stations of the lower bound, the

	optimality gap	| J * |-|J * | |J * |	and the resolution time in seconds; 'γ ' is
	j∈J