Edith Archambault

de l'Enseignement

Relevance

The Ligue de l'Enseignement (Education League) is one of the oldest nonprofit organizations in France; it is also the main champion of the French conception of secularity or laïcité, which is very different from the Anglo-Saxon one. According to this concept, in the public sphere, especially inside public schools and peri-school activity, no reference to faith should be made and no religious signs should be visible. Ligue de l'Enseignement is very close to the public education network and it is the main receiver of public funding.

Introduction

Education League is an umbrella of 30,000 French associations of popular education providing services in the fields of culture, sports, recreation, and continuing education. It advocates for secularity, ideological and religious neutrality, and tolerance inside public schools and more generally inside the public sphere.

Brief History

Founded in 1866 by Jean Macé, a republican journalist, the League fought for a public, lay, free of charge, and compulsory primary education that was finally passed into law in 1882. After this first success in the long-lasting conflict of influence on education between the state and the Catholic Church, the League reoriented toward peri-school activities and civic and adult education. During World War II, the League was suppressed by the authoritarian government of Vichy and refounded after the war. It remained the champion of the lay/secular camp and the public school during the twentieth century when in the 1950s the government began to subsidize the private schools, mainly Catholic, by paying the teachers. After the cooling down of the school war since the 1980s, the League fought racism, discrimination, and social inequality and promoted solidarity, ecology, and sustainable development.

The first mission of the League is to complement the education given by the public school during the postschool hours and in the holidays and vacation periods. Its first target group, of course, is youth inside and outside the school but it addresses the adult population as well, especially those who could not complete their school curricula. It advocates also for secularity/laity, fights discrimination against immigrants and other minorities, and promotes civic behavior and European citizenship.

The main focus areas of the League are: affordable sports and recreation activities with an explicit reference to an access to culture for everybody; education toward solidarity, civics, and tolerance; mainstreaming of immigrant youth and their families; education toward the respect of environment; and sustainable development.

Activities

To fulfill its mission, the League provides a wide range of services for youth and prepares publications intended either for the teachers or for the public at large.

The League gives youth the opportunity to perform arts and sports during their idle hours, inside the school or in community centers, free of charge or at a low price. The League also runs many large facilities used as holiday camps during the vacation periods or welcoming discovery classes during school time. It also organizes linguistic exchanges, the twinning of French and foreign schools or classes. It delivers 10,000 diploma of professional and vocational training for sociocultural animators (BAFA). Since 1999, the League has developed a volunteer association called "read and make read" to prevent the desaffection of children of the digital era for reading books. Lastly, the League supports development projects, mainly in Africa.

Through its website, an annual forum on education, and many debates and publications, the League gives information and references on education, citizenship, secularism, social inclusion, and many other current issues.

Structure and Governance

The Ligue de l'Enseignement is an umbrella of 103 departmental federations grouping 30,000 local associations. It includes also specialized subsidiaries for sports or social tourism. It enjoys strong ties with the main teachers' labor union and the main parent-teacher association. The indirect membership is 1,800,000 families. The League employs 5,000 employees but local associations rely mainly on 500,000 volunteers, either parents or teachers, and 18,000 civic service volunteers. Its board is composed of 36 persons, and its executive bureau includes 17 persons.

The League is affiliated with European platforms on education and with networks or INGOs, such as Solidar and the European Civic Forum.

Funding

The two main resources of the League are membership dues and users' fees on the one hand, and public funding, either subsidies or contracts, on the other hand. Public funding includes the payment of teachers seconded by government to the League. Private donations play a more symbolic role in the League's resources.

Major Accomplishments

The impact of the Ligue de l'Enseignement on the French civil society is deep-rooted as the League played a major role in defining the French concept of secularism, which was difficult to explain to foreigners: the fight against wearing an Islamic scarf inside the school is one example. More generally, the League was, and still is, a major actor in popular education, democratization of culture and human rights defense, and one of the leaders of the nonprofit world in France.

Two Nobel Peace Prizes were awarded to the second and the third presidents of the League, in 1920 and 1927, respectively: Léon Bourgeois, who was also the first president of the League of Nations, and Ferdinand Buisson, who founded the Human Rights League.

Cross-References


Edith Archambault

Introduction

Secours Catholique is a service of the Catholic Church in France. It refers to the social teaching of the Church. Its annual report on poverty in France is a main source of information on this issue.

Brief History

Secours Catholique was founded in 1946 as a declared association by a priest, Jean Rhodain, who was struck by the great poverty in France after World War II. Secours Catholique was dedicated to the alleviation of poverty and operated mainly by volunteers. In 1950, it joined the newly created confederation, Caritas Internationalis, and in 1962, it was acknowledged as a public utility association. Its international action expanded in the 1970s and 1980s.

In 1990, Secours Catholique created a subsidiary, Cités du Secours Catholique, whose mission is to welcome and accommodate, house, and accompany nearly 11,000 people a year, in a situation of social exclusion and/or disability, on the path of integration and autonomy.

In 2009, it was at the origin of Fondation Caritas, a public utility foundation with the same mission in France and in the world as Secours Catholique. In addition, Fondation Caritas allows families to create their own foundations. In 2019, about 100 family foundations are sheltered by the foundation that affords legal advice and financial management.

Mission

According to the bylaws, the mission of Secours Catholique is to expand Christian charity in France and in the world without any national or sectarian particularism and to afford relief, material, and moral help to the poor. Its key values are justice, solidarity, fraternity, and respect for human dignity.

Activities

The activities are shared between France (2/3) and the Third World (1/3).

The main program is to help the poor in France by providing advice and in-kind and financial assistance to children and families. Secours Catholique runs social facilities and services to welcome the homeless and immigrants night and day. It also helps the unemployed to create their own enterprises in collaboration with Caritas Foundation. It affords emergency relief during times of natural disaster or social crisis, either alone or in partnership with public agencies and other nonprofit organizations. On refugees and immigrant issues, it works in cooperation with CIMADE, the main Protestant organization.

Every year about 1,400,000 persons are benefited by Secours Catholique in France. In 2015, Secours Catholique created through Fondation Caritas an impact investment fund intended towards innovative charitable projects in France or Third World.

The international activities of Secours Catholique are mainly located in Africa first and then Asia, besides East Europe and Latin America. Emergency relief and assistance to refugee camps are operated with an aim of achieving sustainable development. The development programs are run through local organizations; they give priority to education and women's empowerment.

The major publication is the 600,000 edited copies of the monthly Messages, sent mainly to donors.

Structure and Governance

Secours Catholique is a confederation of 104 departmental organizations grouping 3,500 local teams. Volunteers are its main human resource: 68,000 volunteers work inside the local teams and their number is growing; however, 950 salaried employees are also included in the paid staff. The governance at the upper level is the responsibility of a board of 16 persons, chaired since 2014 for the first time by a woman, Véronique Fayet, a former politician.

Funding

The annual budget was €135 million in 2017 (plus the money value of volunteering: €203 million). The main financial resource through bequests and gifts comes to about 78% of the annual budget; public subsidies are less than 8%. The net assets are €155 million.

Major Contributions

Secours Catholique is among the most popular charitable organizations and a masterpiece of the Catholic Church in France. Every year, it reports on extreme poverty in France and its evolution. This report gives interesting information to the media, the public at large, and even to the government. Through Fondation Caritas, Caritas Europa, Caritas Internationalis, and other networks, its contribution to the alleviation and the knowledge of poverty and social inequality is decisive.

Cross-References

▶ Caritas Internationalis ▶ Federations, Nonprofit ▶ Foundations, Family ▶ Volunteers


Jing Chang

Wenrui Wang

Damien Chablat
email: damien.chablat@cnrs.fr

Fouad Bennis
email: fouad.bennis@ec-nantes.fr

Evaluating the Risk of Muscle Injury in Football-kicking Training with OpenSim

Keywords: Football, football injury, football kicking, OpenSim, muscle fatigue, fatigue model

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

As one of the most popular and physically demanding games around the world, football is reported to bring about 12 to 35 injury incidences per 1000 hours outdoor games [START_REF] Junge | Soccer injuries[END_REF]. Most of the injury is caused by impact trauma [START_REF] Junge | Soccer injuries[END_REF], followed by overuse injury, which takes 34% of all [START_REF] Nielsen | Epidemiology and traumatology of injuries in soccer[END_REF]. Researches on retired football players have shown that they are faced with twice the incidence of coxarthrosis compared with normal people [START_REF] Klünder | Osteoarthritis of the hip and knee joint in retired football players[END_REF][START_REF] Lindberg | Prevalence of coxarthrosis in former soccer players: 286 players compared with matched controls[END_REF]. In fact, the lifelong occurrence of musculoskeletal problems in former football athletes is reported to be even higher than that of former long-distance runners or shooters [START_REF] Räty | Lifetime musculoskeletal symptoms and injuries among former elite male athletes[END_REF]. Kicking football is a basic training specific for football players. A normal football, weighing about 0.44kg could reach up to 200km/h in a shot, which may bring heavy load to athletes. Body parts that form the kicking-related kinetic chain such as the back, hip, and knee are reported to be the most fragile parts in terms of musculoskeletal injuries [START_REF] Nielsen | Epidemiology and traumatology of injuries in soccer[END_REF]. To prevent related impact and overuse injuries, it is necessary to have an evaluation on the muscle loads as well as a prediction on the fatigue process in football-kicking training.

OpenSim is one of the main biomechanical human modeling platforms [START_REF] Delp | OpenSim: open-source software to create and analyze dynamic simulations of movement[END_REF]. Various kinematical and dynamical calculations are available to simulate human motion as well as to estimate body loads. Furthermore, algorithms such as Static optimization (SO), Compute muscle control (CMC) were developed to estimate the activities of individual muscles [START_REF] Delp | OpenSim: open-source software to create and analyze dynamic simulations of movement[END_REF]. In previous researches, the OpenSim platform has been used to study muscle contributions to sports such as running [START_REF] Hamner | Muscle contributions to propulsion and support during running[END_REF], to determine the risky muscle groups in physical works such as overhead drilling [START_REF] Chang | A full-chain OpenSim model and its application on posture analysis of an overhead drilling task[END_REF], and to study the ways of human-machine interaction from the aspects of body load distribution, such as to compare the unilateral and bilateral crutch gaits [START_REF] Chang | Evaluating the Effect of Crutch-using on Trunk Muscle Loads[C]//International Conference on Human-Computer Interaction[END_REF]. OpenSim will be effective to study the muscle contribution in football-kicking.

Muscle fatigue is an important contributing factor to overuse muscle injuries as well as musculoskeletal disorders. To date, very few studies have investigated the muscle fatigue induced by football-kicking. The muscle fatigue process has been described by the decline of muscle maximal voluntary contraction force (MVC) [START_REF] Ma | A new simple dynamic muscle fatigue model and its validation[END_REF]. As illustrated by the Dynamic muscle fatigue model, the rate of force decline is dependent on the current muscle loads as well as the subject-specific trait. By the time, the muscle MVC declines to a level close to its current load, the risk of overuse injury increases.

In this study, we seek to understand how different muscle groups contribute to football-kicking training as well as to predict the fatigue process in this training. Null hypotheses include: 1) Muscle groups near the acting point of external force (i.e., the calf and knee) contribute more than that of the farther (i.e., the hip and trunk muscles); 2) Muscles of left and right sides differ in contribution. The dominant side that applies football-kicking contributes more than the other.

3) The lower leg and knee muscles contribute especially to accelerating the football.

Results of this study will expand our knowledge about the kinetics of footballkicking and help prevent both impact and overuse injuries in football.

Methods

This is a biomechanical simulation study. A football-kicking motion is simulated and muscle activations with and without the action force of football are calculated with OpenSim. Muscle fatigue process is indicated by muscle MVC declines in constant repetitive football-kicking training.

OpenSim model

In this study, a generic OpenSim model, gait2354 is used for simulation. It is one of the most widely-used models with satisfactory effectiveness (for details see https://simtkconfluence.stanford.edu/display/OpenSim/Gait+2392+and+2354+Models). Muscles are categorized by four groups. The trunk group includes three pairs of muscle (the internal and external obliques, and the erector spine); the hip group includes 18 pairs of muscle that connects the pelvis and leg; the knee group include seven pairs of muscle that cross the knee joint; the calf group includes 4 pairs of muscles.

Data and simulation

For kicking football, athlete lifts the dominant foot and swings forward to reach the ball, during which ground reaction force acts on the supporting foot and footballaccelerating force on the dominant foot. In this study, we compare the case with or without the accelerating force. Motion and force data come from OpenSim document (https://simtk.org/frs/index.php?group_id=679, copyright open from Stanford University), which represents a subject kicking football with the right foot. The kicking motion lasts for 1.5 second.

Figure 1. The simulated football-kicking motion at t = 0.45s (left) and t = 0.9s (right).

Muscle contribution is indicated by computed muscle activation. Inverse kinematics and inverse dynamics are computed to estimate joint positions and joint moments; then Residual reduce algorithm and CMC are applied to estimate the muscle activations.

Muscle fatigue analysis

The Dynamic muscle fatigue model is presented as Eq.1. In this study, we approximate the maximal muscle ability (F max (t)) as the MVC. With the computed muscle activation from OpenSim simulation, the profile of maximal muscle force is predictable. As illustrated in Introduction, the muscle would be at risk for overuse injury by the time its maximal force declines near the required muscle load. In this study, we analyze the muscle fatigue process in repetitive football-kicking motions to check athlete's maximal endurance time. The subject-specific muscle fatigability is set to 1 min -1 , indicating average fatigue resistance [START_REF] Chang | A full-chain OpenSim model and its application on posture analysis of an overhead drilling task[END_REF].

-

With  F cem (t) N Current existing maximal muscle ability

 F load (t) N Muscle load  F max (t) N Maximal muscle ability  k min -1 Subject-specific fatigability  F(t)
1 Muscle activation 3 Results

Muscle contribution in football-kicking

The estimated muscle activations are summarized in Table 1. In general, muscle group that contribute the most are the knee and the hip, instead of the calf; the first hypothesis is not supported. Muscles on the right side (dominant side in this simulation) are more activated, which supports the second hypothesis. However, significant difference only presents in the knee and calf. Finally, contrary to the third hypothesis, acting force from football only makes significant difference in the trunk muscle group. Both the left and right trunk are much more activated in simulation with football (2.18% vs. 24%), suggesting that football acceleration force mainly depends on the trunk muscles.

The muscle force profiles

Muscles that activated the most in each muscle group are chosen to study the fatigue process. They are the iliacus on the hip, the biceps femoris-long head (bifemih) on the knee, the tibialis anterior (tib_ant) of the calf, and the internal oblique (intobl) of the trunk. Mean activation of these muscles are shown in Table 2, and their force profiles are presented in Figure 2. The fatigue process in 200 repetitive football-kicking motions of the four typical muscles are shown in Figure 3456. Red line indicates the muscle force ability while green lines indicates the muscle force requirement at each moment. The intersection of the two indicates the time when current existing muscle force ability declines to the requirement of the motion, so-called the maximal endurance time (MET). After that, the muscle would enter into a risky situation for overuse injuries. It can be found from the figures that the right iliacus in kicking-with-football motion has the shortest MET (12.39s, about 8 cycles), followed by the right introbl in kicking-with-football motion (22.85s, about 15 cycles).

The fatigue process and maximal endurance time

Discussion

Contribution of different groups of muscle in the football-kicking motion

The first surprising result about muscle loads in kicking motion is that the hip muscle group is highly activated (more than 30% on average). This result may explain the high incidence of hip and groin injury in football players [START_REF] Werner | UEFA injury study: a prospective study of hip and groin injuries in professional football over seven consecutive seasons[END_REF][START_REF] Langhout | Hip and groin injury is the most common non-timeloss injury in female amateur football[END_REF][START_REF] Werner | UEFA injury study: a prospective study of hip and groin injuries in professional football over seven consecutive seasons[END_REF][START_REF] Feeley | Hip injuries and labral tears in the national football league[END_REF]. In fact, the hip and groin injury is reported to be the most common non-time-loss injury in football players [START_REF] Langhout | Hip and groin injury is the most common non-timeloss injury in female amateur football[END_REF]. Results of this study also show that the loads of hip muscles present little difference between the swing side and supporting side, or between motions with and without football, which indicates that the hip mainly contributes to posing and stabilizing the legs instead of accelerating the football. Besides, the calf muscles, although near the acting point of football reaction force, are not highly activated. The result emphasizes the importance of a holistic view in examining workloads in ergonomics [START_REF] Chang | The risk assessment of work-related musculoskeletal disorders based on OpenSim[END_REF].

The comparison between simulated muscle loads of the kicking motion with and without football emphasizes the importance of core muscles for football players. The trunk muscle activation is the only that aggravates significantly with the presence of football (about 2% vs. 24%), which indicates that it is the main muscle group contributing to football acceleration. Therefore, stronger trunk muscles may generate better kicking performance. This result expands our previous knowledge about the relationship between strength and football velocity [START_REF] Mcguigan | Strength training for athletes: does it really help sports performance?[END_REF][START_REF] Young | Enhancing foot velocity in football kicking: the role of strength training[END_REF], from whom a new strategy for football kicking training may come out.

Fatigue process for football-kicking training

The fatigue process analysis of the typical muscles has outlined the fragile muscles, in which the iliacus of the hip proves to be at the highest risk of overuse injury with only height continuous football-kicking exercises sustainable. This muscle is responsible for stabilizing the pelvis in contralateral hip extension and posing the legs [START_REF] Andersson | The role of the psoas and iliacus muscles for stability and movement of the lumbar spine, pelvis and hip[END_REF]. Another muscle at high risk for overuse injury is the internal obliques of the trunk muscle group, with about 15 football-kicking exercises sustainable. Previous research has shown that a larger core muscle reduces the number of missed matches due to injury [START_REF] Hrysomallis | Injury incidence, risk factors and prevention in Australian rules football[END_REF]. In summary, results of this study suggest that the strength training for the hip and the core should be ahead of specific football-kicking training to prevent injuries.

Conclusions

In this study, we examine the muscle behavior in kicking motion with and without football using OpenSim, then evaluate the fatigue process and MET of repetitive football-kicking training by applying the Dynamic Muscle Fatigue model. It is concluded that muscles of the hip and the knee joints are highly activated in kicking motion, and the football acceleration force mainly depends on the trunk muscles. Muscle at the highest risk for fatigue and overuse injury is the iliacus of the hip joint, which, according to the settings of this study, sustains only height continuous footballkicking cycles. The results may offer a reference for the coaches, athletes, and physiotherapists about the football-kicking training and related injuries.

Figure 2 .

 2 Figure 2. The estimated muscle exertions in a football-kicking motion. (a) iliacus of the hip; (b) biceps femoris-long head (bifemlh) of the knee; (c) tibialis anterior (tib_ant) of the calf; (d) internal oblique (interobl) of the trunk.

Figure 3 .

 3 Figure 3. The fatigue process of right and left iliacus on the hip.

Figure 4 .

 4 Figure 4. The fatigue process of right and left bifemih of the knee.

Figure 5 .

 5 Figure 5. The fatigue process of right and left tib_ant on the calf.

Figure 6 .

 6 Figure 6. The fatigue process of right and left intobl on the trunk.

Table 1 .

 1 The calculated muscle activations in football-kicking training. Muscles are compared between left and right side.

	Muscle group	Number of muscles	Muscle activation (%) R-with R-without football football	L-with football	L-without football
	Hip	18	35.64	36.77	33.23	32.11
	Knee	7	40.54	41.53	28.65	27.20
	Calf	4	18.22	17.49	6.85	5.48
	Trunk	3	2.18	24.21	2.18	23.97

Table 2 .

 2 The typical muscle activation of each muscle group in football-kicking training.

	Muscle group	Name of muscle	R-with football	Muscle activation (%) R-without L-with football football	L-without football
	Hip	Iliacus	67.1	74.9	65.0	65.5
	Knee	Bifemih	35.8	33.1	25.9	26.9
	Calf	Tib_ant	31.8	32.3	15.0	8.9
	Trunk	Intobl	2.1	29.8	2.1	29.1


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

Philippe Castagliola
email: philippe.castagliola@univ-nantes.fr

Alexandre Dolgui
email: alexandre.dolgui@imt-atlantique.fr

Faicel Hnaien
email: faicel.hnaien@utt.fr

A hybrid genetic algorithm for a multilevel assembly replenishment planning problem with stochastic lead times

Keywords: Assembly systems, Assemble-to-order, Replenishment planning, Stochastic lead-times, Genetic algorithms, Discrete optimization

) by using a recursive procedure and a Branch and Bound algorithm. Here, in order to decrease the computation time, a novel generalized probabilistic model based on an iterative approach is developed. The proposed model calculates the expected total cost, which is composed of the inventory holding cost for components and the backlogging and inventory holding costs for the finished product.

An iterative approach and a hybrid genetic algorithm are introduced to determine the planned order release dates for components at the last level of the bill of materials that minimizes the expected total cost. Experimental results show that the proposed optimization algorithm efficiently finds good-quality approximate solutions regardless of the type of assembly system, the number of components at the last level and the variability of the finished product-related costs.

Introduction

In an assemble-to-order (ATO) environment, the manufacturer needs several components to assemble a finished customized product. The demand is not known in advance and there is no stock of finished products to anticipate it.

The customer asks for a product composed of a given set of standard or personalized components. The producer uses information on the lead-times of the components (time to order them from suppliers or to make and assemble them) and the assembly process to set the delivery time to the customer (d'Avino et al., 2013). This delivery time is often used as the due date [START_REF] Hammami | Supply chain design to guarantee quoted lead time and inventory replenishment: model and insights[END_REF]. The producer's Manufacturing Resource Planning (MRP) system uses the input data (e.g., demand and its due date, assembly times, delivery times) to calculate the assembly order release dates for components. If the actual lead times are stochastic, the calculation of planned lead times, and thus, release dates for the components, is a complex problem.

ATO strategies are widely used in industry. However, ATO systems are very sensitive to uncertainty, which can disrupt the ATO process and its MRP system [START_REF] Milne | Optimizing planned lead times for enhancing performance of MRP systems[END_REF]. The assembly process can be interrupted due to supply variations [START_REF] Flynn | On theory in supply chain uncertainty and its implications for supply chain integration[END_REF][START_REF] Simangunsong | Supply-chain uncertainty: a review and theoretical foundation for future research[END_REF][START_REF] Wazed | Uncertainty factors in real manufacturing environment[END_REF] or disruptions [START_REF] Snyder | OR/MS models for supply chain disruptions: a review[END_REF][START_REF] Speier | Global supply chain design considerations: Mitigating product safety and security risks[END_REF][START_REF] Kleindorfer | Managing disruption risks in supply chains[END_REF]. In practice, machine breakdowns can interrupt the assembly process. Component customization and replenishment lead-times may also be uncertain and significantly longer or shorter than planned. This leads to stochastic component replenishment lead times. If only one component is delayed, the entire assembly process is stopped. The holding cost for other components, already delivered, increases, and delivery of the finished product may be delayed.

MRP logic with deterministic assumptions often proves to be too limited in ATO environments. The majority of the existing academic literature over the past 30 years has argued that MRP logic with deterministic assumptions about lead time is too restrictive. In fact, the random variability of lead-times drastically decreases the system's performance and this leads to a need to estimate planned lead times far more precisely [START_REF] Bandaly | Impact of lead time variability in supply chain risk management[END_REF]. A possible solution to cope with this difficulty consists of introducing safety lead-times: planned lead-times are assumed to be equal to the contractual plus additional safety lead times [START_REF] Dolgui | Supply planning under uncertainties in MRP environments: A state of the art[END_REF]. Nevertheless, in several survey papers [START_REF] Dolgui | A state of the art on supply planning and inventory control under lead time uncertainty[END_REF][START_REF] Dolgui | Supply planning under uncertainties in MRP environments: A state of the art[END_REF][START_REF] Damand | Parameterisation of the MRP method: automatic identification and extraction of properties[END_REF][START_REF] Koh | Uncertainty under MRP-planned manufacture: Review and categorization[END_REF] which have investigated how MRP systems address lead-time uncertainties, it was concluded that a vast amount of literature highlighting the use of safety stocks to handle them, and that the safety lead-times have long been neglected [START_REF] Van Kampen | Safety stock or safety lead time: coping with unreliability in demand and supply[END_REF]. Furthermore, there have been very few studies focusing on how to minimize the additional costs caused by the lead-time uncertainty [START_REF] Dolgui | Supply planning under uncertainties in MRP environments: A state of the art[END_REF] or how to anticipate them [START_REF] Jansen | Lead time anticipation in supply chain operations planning[END_REF]. For more details on replenishment planning models under uncertainty, interested readers can refer to [START_REF] Aloulou | A bibliography of non-deterministic lot-sizing models[END_REF], [START_REF] Díaz-Madroñero | A review of discrete-time optimization models for tactical production planning[END_REF], [START_REF] Dolgui | A state of the art on supply planning and inventory control under lead time uncertainty[END_REF], [START_REF] Dolgui | Supply Chain Engineering : Useful Methods and Techniques[END_REF], [START_REF] Ko | A review of soft computing applications in supply chain management[END_REF], [START_REF] Peidro | Quantitative models for supply chain planning under uncertainty: A review[END_REF], Mula et al. (2006b) and [START_REF] Koh | Uncertainty under MRP-planned manufacture: Review and categorization[END_REF].

Despite the importance of this issue, as far as we know, only a small number of researchers have studied multi-level assembly systems in an ATO environment with non-deterministic lead-time behavior [START_REF] Atan | Assemble-to-order systems: A review[END_REF]. In this paper, this gap is addressed by investigating an ATO supply network to assemble a given tailored finished product under lead time uncertainty. The customer's request defines this finished product and the set of customizable components needed to customize it. In the considered case, there are no product or component stocks available to anticipate this demand and so it is necessary to set a due date for client delivery. To do this, the planners use information about the tailored product, customizable components, supplier availabilities and the assembly process in order to design the whole supply network, define the list of partners (local assembly units or external suppliers) and estimate the planned lead-times. As highlighted by [START_REF] Golini | Moderating the impact of global sourcing on inventories through supply chain management[END_REF], [START_REF] Chandra | Inventory management with variable lead-time dependent procurement cost[END_REF], [START_REF] Berlec | Predicting order lead times[END_REF] and [START_REF] Arda | Inventory control in a multi-supplier system[END_REF], this information is essential at the negotiation step to design the supply network, estimate lead-times and release the orders.

There is usually little information available on the capacity and state of resources and production management at supply chain partner enterprises (suppliers). The only information available to help assess the lead-time of each partner is its probability distribution based on past statistics. These include not only processing times but also additional times related to workload, capacity constraints, local planning decisions, variations, etc. Due to the complexity and, often, the impossibility (e.g. confidentiality) of considering all the details on the state of supply chain partners, this approach has escaped criticism and has been highly recommended in the literature on contracting and planning under uncertainty [START_REF] Chen | Outsourcing contracts and ordering decisions of a supply chain under multidimensional uncertainties[END_REF][START_REF] Jabbarzadeh | Closed-loop supply chain network design under disruption risks: A robust approach with real world application[END_REF][START_REF] Giri | Improving performance by coordinating a supply chain with third party logistics outsourcing under production disruption[END_REF][START_REF] Carvalho | Supply chain redesign for resilience using simulation[END_REF][START_REF] Berlec | Predicting order lead times[END_REF][START_REF] Fiala | Information sharing in supply chains[END_REF][START_REF] Song | Contract assembly: Dealing with combined supply lead time and demand quantity uncertainty[END_REF].

In this paper, a multi-level assembly system with stochastic component lead-times is studied. The aim is to calculate planned lead-times for all components of a given product order at a given due date, in the case of random variability of real component lead-times. The study is composed of the following steps:

• a new approach and efficient model for performance evaluation (expected cost calculation) is suggested;

• an overall discrete optimisation model is proposed;

• an efficient approximate algorithm to solve the proposed model is developed;

• numerical tests and analysis of model behaviour are carried out;

• conclusions and managerial insights are given.

The remainder of the paper is organized as follows: Section 2 reviews the relevant literature; Section 3 describes the problem; Section 4 presents the generalized probabilistic model; Section 5 describes the optimization approach; Section 6 gives performance tests; and Section 7 rounds off with a conclusion; managerial insights; and perspectives for future research.

Related publications

The literature features various probabilistic approaches for dealing with lead-time uncertainty in the case of an MRP environment (readers can refer to [START_REF] Díaz-Madroñero | A review of discrete-time optimization models for tactical production planning[END_REF]; [START_REF] Dolgui | A state of the art on supply planning and inventory control under lead time uncertainty[END_REF]; [START_REF] Dolgui | Supply planning under uncertainties in MRP environments: A state of the art[END_REF]Mula et al. (2006b)). Analytical probabilistic approaches often model the simplest structures (e.g. for one-level systems) and give way to simulation methodology for more complex structures (multi-level systems, multi-period planning, etc.).

Several studies have modeled one-or two-level assembly systems under a probabilistic distribution of lead-times (see Table 1). The proposed models consider continuous or discrete random lead-times, but remain limited:

• To one-level in the bill of materials (BOM), with constant demand for all periods;

• Or to two-levels in the BOM, with mono-period planning (including one customer demand).

For the case of multi-level BOMs, [START_REF] Ben-Ammar | Optimization approaches for multi-level assembly systems under stochastic lead times[END_REF], 2018) were the first to propose an analytical approach to model a one customer demand. There are also studies based on fuzzy logic. In the area of fuzzy reasoning, statistical data are assumed to be either too unreliable or insufficient to obtain reliable probability distributions. Mula et al. (2006b) presented an excellent review of supply and production planning under uncertainty, in which they noted that most studies have focused on handling the uncertainty of demand [START_REF] Guillaume | Management of the risk of backorders in a MTO-ATO/MTS context under imperfect requirements[END_REF][START_REF] Peidro | Quantitative models for supply chain planning under uncertainty: A review[END_REF][START_REF] Grabot | Integration of uncertain and imprecise orders in the MRP method[END_REF]Mula et al., 2006a). This under-focus on the supply side is explained by the difficulty of modeling leadtime variability. Most studies have considered lead-times as nil or fixed. Without attempting an exhaustive review of the literature, we specifically focus here on probabilistic approaches. [START_REF] Dolgui | A state of the art on supply planning and inventory control under lead time uncertainty[END_REF] presented a literature review classifying the techniques and approaches in the field of supply planning and inventory control systems under stochastic lead-times. As mentioned earlier and described in Table 1, most models have assumed (i) one-level assembly systems and mono-or multi-period planning or (ii) two-level assembly systems and mono-period planning. Yano (1987a,b) were among the first to note the effects of the variability of lead-times. They proposed the case of serial production systems and mono-period planning, which was later generalized by [START_REF] Elhafsi | Optimal leadtimes planning in serial production systems with earliness and tardiness costs[END_REF]. In a third study [START_REF] Yano | Stochastic leadtimes in two-level assembly systems[END_REF] studied a two-level assembly system and one-period planning. One component is assembled from two components and lead-times follow stochastic Poisson and negative binomial distributions. The objective function is the expected total cost (ETC), which is composed of tardiness and inventory holding costs. An algorithm based on an exploration of the properties of the objective function determines the optimal planned lead-times. [START_REF] Kumar | Component inventory costs in an assembly problem with uncertain supplier lead-times[END_REF] studied the case of one-period planning and one-level assembly systems. Like Yano (1987a), he minimized the total expected cost and proposed optimal order release dates. An exact analysis was proposed with assumptions of special types of distribution (normal, uniform and exponential). A few years later, [START_REF] Chu | Supply management in assembly systems[END_REF] looked at the case of a one-level assembly system: they proved the convexity of the mathematical expectation of total cost, and proposed an iterative algorithm to minimize it. To solve this issue for the case of a two-stage model production model, [START_REF] Hegedus | Setting procurement safety lead-times for assembly systems[END_REF] proposed an approach based on a Newsboy-like analytical formulation. The main weakness in these studies is that they are limited to one-or two-level assembly systems.

Next, [START_REF] Axsäter | Planning order releases for an assembly system with random operation times[END_REF] considered a three-level assembly network with independent stochastic operation times. The objective was to choose starting times for different operations in order to minimize the ETC, composed of the inventory holding costs of components and the tardiness cost of end items. An approximate decomposition technique, based on the repeat application of the solution of a single-stage problem, was suggested. In the numerical result for the first example of a two-stage problem, the error was only 1% error, i.e. the relative cost increase due to approximation.

However, for the second three-stage example, the error was about 10%, which cannot be disregarded.

Earlier, [START_REF] Dolgui | Planification de systèmes d'assemblage avec approvisionnements aléatoires en composants[END_REF] and Dolgui & Ould-Louly (2002) had developed an approach to multi-period planning based on coupling simulation with integer linear programming to model one-level assembly systems. Several types of finished product were produced and assembled from different types of components. An inventory holding cost for each component and a backlogging cost for the finished products were introduced. For each period, the authors optimized the quantity of finished products to be assembled and the quantity of components to be ordered. [START_REF] Proth | Supply management for cost minimization in assembly systems with random component yield times[END_REF] studied the same problem and added a heuristic to select which products were to be assembled and which components were to be ordered.

For the case of multi-period models and one-level assembly systems, Ould-Louly & Dolgui (2002a,b, 2004, 2009, 2011, 2013), Ould-Louly et al. (2008b,a) and [START_REF] Shojaie | A study on MRP with using leads time, order quality and service level over a single inventory[END_REF] presented mathematical models with random leadtimes, a known finished product demand (the same for all periods) and unlimited assembly capacity. [START_REF] Hnaien | Robust single-item lot-sizing problems with discrete-scenario lead time[END_REF] recently considered the min-max robust lot-sizing problem with discrete lead-time scenarios, and provided a complexity analysis proving that robust lot-sizing problems are NP-hard even when there are two scenarios. They showed that several optimality conditions for the deterministic cases provided in [START_REF] Wagner | Dynamic version of the economic lot size model[END_REF], as well as a classic facility location-based model, were no longer valid. To analyze the effect of lead-time uncertainty, they measured different indicators (backlogging costs, inventory-holding costs, service level and set-up cost), modeled some policies (periodic order quantity (POQ), lot for lot (L4L), economic order quantity (EOQ)) and used several optimization approaches to find optimal release dates for components.

In the cases of one-period planning and two-level assembly systems, [START_REF] Tang | The detailed coordination problem in a two-level assembly system with stochastic lead times[END_REF] modeled both process times and lead-times for components. They considered a fixed demand, unlimited capacity, and a known due date. A Laplace procedure was used to calculate the optimal safety lead-times to minimize backlogging and inventory holding costs. Later, [START_REF] Hnaien | Genetic algorithm for supply planning in two-level assembly systems with random lead times[END_REF] and Fallah-Jamshidi et al. ( 2011) treated the same problem using genetic algorithms (GAs) to minimize the expected value of the same total cost. [START_REF] Hnaien | Genetic algorithm for supply planning in two-level assembly systems with random lead times[END_REF] supposed that components at level 1 of the BOM were stored and that the finished product was only assembled after the given due date. Fallah-Jamshidi et al. ( 2011) explored the same problem but in a multi-objective context, and reinforced the GA by a reliable technique called the electromagnetism-like mechanism.

Hnaien et al. ( 2016) studied the case of a one-period inventory model for a one-level assembly system under stochastic demand and lead-times. An analytical model and a Branch and Bound (B&B) approach were found to optimize the component quantities and planned lead-times. Although the authors limited their model to one period, it can be extended to multi-period planning to account for possible trade-offs between stocks from different periods.

Building on this work, [START_REF] Borodin | Component replenishment planning for a single-level assembly system under random lead times: A chance constrained programming approach[END_REF] proposed a new approach based on chance-constrained programming and a linear model to solve the same problem. However, they limited themselves to one-level assembly systems because of the dependency between levels. This obstacle can easily be overcome, as in this paper, using separate lead-time uncertainty models at each level. [START_REF] Atan | Setting planned leadtimes in customer-orderdriven assembly systems[END_REF] recently considered a final assembly process fed by a multi-stage parallel process. Each stage had a stochastic throughput time and the system was controlled by planned lead-times at each stage. The authors developed an iterative heuristic to optimize the planned lead-times for different stages and to minimize the expected cost of a customer order. The problem studied contained a one-assembly process which, if modeled differently by adding planned assembly dates at intermediate levels, could be considered a special case of our study in this paper.

Therefore, it could be modeled in this way and resolved using our optimization approach. More recent papers [START_REF] Jansen | Setting optimal planned leadtimes in configure-to-order assembly systems[END_REF][START_REF] Jansen | Newsvendor equations for production networks[END_REF] feature an interesting investigation into the combination of planned lead-time and commitment lead-time decisions. In Ben-Ammar et al. ( 2018), the authors provided a generalized probabilistic model and a B&B algorithm to optimize multi-level assembly systems for the case of one customer demand. The objective function is composed of inventory holding costs for components at all levels, and holding and backlogging costs for the finished product. This total cost is expressed as a recursive procedure taking into account inter-level dependency. However, the limitation of this earlier study was that each lead-time varies between 1 and a known upper limit, whereas we observe that the recursive function requires a recursively enumerable domain and depends on the lead-time distribution functions. The major drawback of this recursive procedure is its influence on computation time.

The purposes of this paper, based on the modeling approach presented in Ben-Ammar et al. (2018), are (i) to develop a more general approach than is found in the models introduced by Yano (1987a), [START_REF] Chu | Supply management in assembly systems[END_REF], Tang

C 1,m C i,m C N m ,m C 1,2 C i,2 C N 2 ,2 C 1,1 C i,1 C N 1 ,1
Order release dates

Problem assumptions

To be closer to real-world industrial planning methods, we consider a discrete temporal environment and integer decision variables. Figure 1 shows that the finished product is produced from components that are themselves obtained from the following level's components and so on.

To satisfy the customer's request for a tailored finished product, the planner configures the whole supply network according to the desired delivery date and so has to know when the overall processes should be released. Note that our approach is focused on the case of contracting with our customer in an ATO environment. We suppose that all supply chain partners are independent enterprises and that the supply chain will be managed in a decentralized manner.

As in Ben-Ammar et al. ( 2018), we only determine order release dates for components at the last level of the BOM and we assume that there is no decision variable for internal levels (no possibility for accounting for future local decisions by supply chain partners, and we only know the lead-time probability distributions of our partners). The production cycle thus extends from order release dates of components at the last level to finished-product delivery date. The assembly system is composed of m levels. At each level, we consider ∀l = 1, . . . , m, N l components. In this way, ∑ m l=1 N l components are needed to assemble the finished product. Let us assume that: (i) only components c i,m (i = 1, . . . , N m ) at the last level m are ordered from external suppliers and their order release dates X i,m (i = 1, . . . , N m ) have to be defined; (ii) in local assembly units, components c i,l (i = 1, . . . , N l ) at intermediate levels (l = 1, . . . , m-1) are assembled without any decision possible on their start dates; (iii) suppliers and local assembly units are independent.

As in Ben-Ammar et al. ( 2018), without loss of generality, we introduce the following assumptions: the finished product demand D is known and equal to 1, and to assemble it, only one unit of each component is required. The unit backlogging cost b and the unit inventory holding cost r for the finished product, and the unit inventory holding cost h i,l for the component c i,l , are known.

The following notations are used in this paper:

Parameters T Due date for the finished product, T > 0 D Demand (known) for the finished product at the date T , without loss of generality, let

D = 1 l Level in a bill of material (BOM), l = 1, . . . , m c i,l Component i of level l of the BOM N l Number of components of level l S i,l Set of components needed to assemble component c i,l L i,l Random lead-time for component c i,l t i,l Minimum value of c i,l u i,l Maximum value of L i,l ; each L i,l varies in [t i,l , u i,l ] U i,m
The longest time between the release date for component c i,m and T . Equal to the maximum value of

∑ m v=1 L i v ,v ; ∀θ ∈ [1, m -1] and ∀c i θ +1 ,l ∈ S i θ ,l-1 T i,m
The shortest time between the release date for component c i,m and T . Equal to the minimum value of

∑ m v=1 L i v ,v ; ∀θ ∈ [1, m -1] and ∀c i θ +1 ,l ∈ S i θ ,l-1 h i,l
Unit holding cost for component c i,l per period We use the following notations to simplify several expressions:

• Assembly date for c i,m-1 : M i,m-1 = max c k,m ∈S i,m-1 (X k,m + L k,m ) • Assembly date for c i,l-1 , ∀ 2 ≤ l ≤ m -1, ∀i ∈ [[1, N l ]]: M i,l-1 = max c k,l ∈S i,l-1 (M k,l + L k,l )
• Assembly date for the finished product:

M FP = max i∈[[1,N 1 ]] (M i,1 + L i,1 )
• Maximum between M FP and due date T :

M + FP = max(M FP , T )
• Minimum between M FP and due date T :

M - FP = min(M FP , T ) • ∑ N l i=1 H i = ∑ N l i=1 h i,l -∑ c k,l+1 ∈S i,l h k,l+1 • H = ∑ N 1 i=1 h i,1 + b • R = r -∑ N 1 i=1 h i,1

Expected cost calculation

The main issue for this problem is how calculate the expected cost for different values of decision variables. The total cost C(X, L) is the sum of the inventory holding or backlogging cost for the finished product and the inventory holding cost for components. Figure 2 gives an example of when the finished product is assembled after the due date.

We note that p i, j is the planned lead time for component c i, j and PM i, j is its planned assembly date.

Because of the lead-time uncertainties:

• There is a backlog for the finished product if at least one type of component at level 1 is delivered after the due date T . Then the corresponding backlogging cost is equal to b(M + FP -T ).

• If all components c i,1 ,∀i = 1, . . . , N 1 , at level one are available before T , the finished product may be assembled and stored. The corresponding inventory holding cost is equal to r(T -M - FP ).

• The components c k,m at the last level m are ordered at date X k,m and delivered at date L k,m + X k,m . Assembly of component c i,m-1 begins when all the necessary components c k,m ∈ S i,m-1 are available, i.e. at the date M i,m-1 .

The holding cost of components c k,m at the last level m is equal to:

N m-1 ∑ i=1 ∑ c k,m ∈S i,m-1 h k,m M i,m-1 -(X k,m + L k,m )
• There are inventories for components c i,l (l = 2, . . . , m -1 and c k,l ∈ S i,l-1 ) during the time period between their arrival at M k,l + L k,l and M i,l-1 , which is the assembly date for component c i,l-1 . The corresponding holding cost is equal to: • There are inventories for components c i,1 , ∀i = 1, . . . , N 1 during the time period between their arrival at M i,1 +L i,1 and M FP , which is the assembly date for the finished product. The corresponding inventory holding cost is equal to:

m-1 ∑ l=2 N l-1 ∑ i=1 ∑ c k,l ∈S i,l-1 h k,l M i,l-1 -(M k,l + L k,l ) Time Components c 1,3 X 1,3 p 1,3 L 1,3 • PM 1,2 M 1,2 PM 1,1 X 2,3 • X 3,3 • • • M 1,1 T M FP X 5,3 c 5,3 X 4,3 M 2,1 PM 2,
N 1 ∑ i=1 h i,1 M FP -(M i,1 + L i,1 )
The total cost C(X, L) is a random discrete variable (because the lead-times L i, j , ∀i = 1, . . . , N j and ∀ j = 1, . . . , m, and assembly dates M i, j are random variables). Its explicit form is demonstrated in Ben-Ammar et al. ( 2018) and it reads as follows:

C(X, L) = N 1 ∑ i=1 h i,1 M FP + m-1 ∑ l=1 N l ∑ i=1 H i M i,l - m ∑ l=1 N l ∑ i=1 h i,l L i,l - N m ∑ i=1 h i,m X i,m + b(M + FP -T ) -r(T -M - FP ) (1)
with L = (L 1,1 , . . . , L i,1 , . . . , L N 1 ,1 , . . . , L 1,m , . . . , L i,m , . . . , L N m ,m ) and X = (X 1,m , . . . , X i,m , . . . , X N m ,m ).

203

An explicit form of the mathematical expectation of total cost E(C(X, L)) was proposed in Ben-Ammar et al.

(2018). Here we develop a new explicit form that avoids using a recursive function and allows us to model actual lead-times L i,l whose variations are between any two limits t i,l and u i,l . Note that integrating these two limits, which are time intervals in which assembly dates M i,l vary, is calculated to substantially reduce computation times. Let Γ be a positive random discrete variable with a finite number of possible values and F Γ (.) be its cumulative distribution function. Its expected value is equal to:

E(Γ) = ∑ s≥0 (1 -P(Γ ≥ s)) = ∑ s≥0 (1 -F Γ (s)) (2)
In [START_REF] Hnaien | Genetic algorithm for supply planning in two-level assembly systems with random lead times[END_REF][START_REF] Hnaien | Multi-objective optimization for inventory control in two-level assembly systems under uncertainty of lead times[END_REF]

E(M i,l ) = ∑ s≥0 (1 -Q + (L i,l , s, l)) (3) E(M - FP ) = ∑ 0≤s≤T -1 (1 -Q + (L i,1 , s, 1)) (4) E(M + FP ) = T + ∑ s≥T (1 -Q + (L i,1 , s, 1)) (5)
Nevertheless, as mentioned in the Appendix in Ben-Ammar et al. ( 2018), calculating the expected value for an assembly date of a given item (component or finished product) requires using the probability distributions of all the components of which it is composed. Here we propose an iterative approach to calculate the cumulative distribution function for each assembly date.

Proposition 4.1.

F i,m-j-1 (s) =              ∏ c k,m ∈S i,m-1 F k,m (-X k,m + s) i f j = 0 ∏ c k,m-j ∈S i,m-j-1 ∑ o j +w j =s o j +w j ∈N P(L k,m-j = o j )F k,m-j (w j ) i f j ∈ [[1, . . . , m -1]] ( 6
)
where F 0,0 (s) is the cumulative distribution function of the finished-product assembly date.

Proof. see Appendix.

In [START_REF] Hnaien | Genetic algorithm for supply planning in two-level assembly systems with random lead times[END_REF][START_REF] Hnaien | Multi-objective optimization for inventory control in two-level assembly systems under uncertainty of lead times[END_REF]

and ∀i ∈ [[1, N 2 ]],
varied between 1 and a fixed upper limit. These studies considered assembly dates as random discrete variables s varying between 0 and T +U -2 with U = max (U i,2 ) and U i,2 = u k,2 + u i,1 the maximum value of L k,2 + L i,1

(c k,2 ∈ S i,1 , L k,2 + L i,1 ∈ [[2,U k,2 ]]
). In Ben-Ammar et al. ( 2018),the authors generalized these limits and defined s as a natural number such as 0 ≤ s ≤ T + m (u -1). Here we suppose that each lead-time

L i,l , ∀l ∈ [[1, m]] and ∀i ∈ [[1, N m ]]
varies between t i,l and u i,l . The mathematical expectation of the total cost, introduced in Ben-Ammar et al. ( 2018), can be simplified using expression (6).

Proposition 4.2. The mathematical expectation of the total cost, noted E(C(X, L)), can be simplified and written as follows:

E(C(X, L)) = R • T -1 ∑ s=α 0,0 F 0,0 (s) - m-1 ∑ l=1 N l ∑ i=1 H i • β i,l - β i,l -1 ∑ s=α i,l F i,l (s) + H • β 0,0 - β 0,0 -1 ∑ s=T F 0,0 (s) -b • T - m ∑ l=1 N l ∑ i=1 h i,l E(L i,l ) - N m ∑ i=1 h i,m E(X i,m ) (7)
where:

α i,l = max c k,l+1 ∈S i,l α k,l+1 + t k,l+1 with α i,m = X i,m and α 0,0 = max i∈[[1,N 1 ]]
(α i,1 + t i,1 ).

β i,l = max c k,l+1 ∈S i,l β k,l+1 + u k,l+1 with β i,m = X i,m .
Proof. see Appendix.

The intervals T -U i,m ≤ X i,m ≤ Tm are the initial search space, which depends on maximum and minimum lead-times and on the number of levels. To reduce their upper limits, the multi-level assembly system is decomposed to N m (the number of components at level m) multi-level linear supply chains. A finished product is delivered by each

linear chain i, i ∈ [[1, N m ]
] on a specified delivery date ψ i . Two costs are taken into account: (i) if a given finished product is delivered after the due date T , a backlogging cost is considered, and (ii) if it arrives before T , an inventory holding cost is assumed. The optimal order release date, noted X * * i , for one linear chain is used to reduce the initial search space for the corresponding component release date at the last level in the BOM.

T -U i,m ≤ X i,m ≤ X * * i ∀i ∈ [[1, N m ]]
where the optimal order release date X * * i satisfies the optimality condition for the discrete Newsboy model:

F (T -X * * i -1) ≤ b b + r ≤ F (T -X * * i ) ∀i ∈ [[1, N m ]] (8)
and where F (.) is the cumulative distribution function of the total lead-time L associated with linear chain i.

In order to solve this non-linear probabilistic problem, in Ben-Ammar et al. ( 2018) the authors developed bounds and proposed a B&B procedure. Its efficiency depends on the number of levels and the ratio between b (the unit backlogging cost of the finished product per period) and r (the unit inventory holding cost for the finished product per period). In this work, we observed that optimal solutions are computed exactly for small and medium-sized problems with at most 40 components and 3 levels in the BOM. Nevertheless, this branch-and-bound-based approach may not be practical for assembly systems with more than 3 levels in the BOM, as it carries the key limitation that it depends on both (i) backlogging and inventory holding costs and (ii) the number of levels in the structure.

In Ben-Ammar et al. ( 2013), the authors tried to develop metaheuristics to determine good quality approximate solutions for larger problems. To achieve a more effective search, a first optimization approach was based on the integration of the classic genetic algorithm (CGA). However, we observed major drawbacks with this approach: (i)

CPU time depends heavily on the number of levels in the assembly system, and (ii) a CGA is not guaranteed to converge to a global minimum -even though the stopping criterion is fixed to 1000 iterations and the number of components at the last level m (size of each individual). In Ben-Ammar et al. ( 2016), several techniques were proposed to reinforce the CGA, such as a reduction in the initial research space (RSR), a local search (LS), and an integration of perturbation (P). These techniques only seem to be efficient enough to optimize two-level assembly systems with less than 100 components at the last level of the BOM, even though the authors did not determine the set of parameter values that yields the best performance of the hybrid genetic algorithm (HGA).

Even though these works improved the efficiency of the HGA, most of the improvements were scattered across several studies that are not realistically comparable in terms of methodology and results. Nonetheless, it is possible to improve the efficiency of the approach further, and group the results. With this goal in mind, the next section describes the techniques used to reinforce the CGA, and explains the experiments carried out in order to study the behavior of the HGA.

Optimization model and approach

The goal is to find the optimal order release dates for the components at level m to minimize the total expected cost E (C (X, L)). The optimisation problem is as follows:

min(E (C (X, L)))

Subject to constraints:

E (M FP ) = β 0,0 - β 0,0 -1 ∑ s=T F 0,0 (s) - T -1 ∑ s=α 0,0 F 0,0 (s) (9)
E M + FP = β 0,0 - β 0,0 -1 ∑ s=T F 0,0 (s) (10) E M - FP = T - T -1 ∑ s=α 0,0 F 0,0 (s) (11) E M i,l = β i,l - β i,l -1 ∑ s=α i,l F i,l (s) ∀i ∈ [[1, N l ]], ∀l ∈ [[1, m -1]] (12) T -U i,m ≤ X i,m ≤ X * * i ∀i ∈ [[1, N m ]] (13)
X i,m ∈ N ∀i ∈ [[1, N m ]] (14)
where:

• α i,l = max c k,l+1 ∈S i,l α k,l+1 + t k,l+1 with α i,m = X i,m and α 0,0 = max i∈[[1,N 1 ]] (α i,1 + t i,1 ), • β i,l = max c k,l+1 ∈S i,l β k,l+1 + u k,l+1 with β i,m = X i,m .
The objective function which expresses the expected cost is nonlinear with integer variables, see ( 7). The assembly date for the finished product is equal to max

i∈[[1,N 1 ]] (M i,1 + L i,1
). Constraint (9) gives its mathematical expectation.

Constraints ( 10) and ( 11), respectively, gives the mathematical expectation of the maximum/minimum, between the assembly date M FP and the due date T . Constraints (12) express the mathematical expectation of assembly dates for components c i,l . The search spaces for order release dates are defined in constraints (13). Finally, constraints ( 14) are the integrity constraints.

This study made partial use of the optimization approach introduced by Ben-Ammar et al. ( 2016). In our new approach, the CGA presented by [START_REF] Hnaien | Genetic algorithm for supply planning in two-level assembly systems with random lead times[END_REF] is reinforced by several techniques such as a reduced space of research (RSR), a local search (LS), and the integration of perturbation (P).

Let us use the following set of parameters:

• x 1 Maximum iteration number,

• x 2 Population size,

• x 3 Crossing-over probability,

• x 4 Mutation probability,

• x 5 Number of individuals from the initial population that can be a local minimum,

• x 6 Number of iterations that have passed without the best solution being improved (to apply the perturbation approach),

• x 7 Number of iterations that have passed without the best solution being improved (to apply the local search).

The various operations of the proposed HGA are described in Algorithm 1. The stop criterion is that the maximum iteration number, noted x 1 , is reached. The initial population is formed by x 2 individuals. Note that the chromosome representation, the mutation, the crossover and the neighborhood search are done in the same manner as in [START_REF] Hnaien | Genetic algorithm for supply planning in two-level assembly systems with random lead times[END_REF]. However, the only differences concerning these operations are that (i) the length of a given chromosome is N m , (ii) a given chromosome contains a sequence of order release dates X i,m (integer decision variables), (iii) the crossover and mutation probabilities, noted x 3 and x 4 , will be recalibrated.

As shown in Algorithm 1, the proposed HGA includes several steps. The following subsections discuss some of these steps in detail.

14 S Pop ←-Initial Population(x 2 , x 5 ); Choosing the right population, of the right size, is crucial because it can affect the performance of the algorithm.

foreach j ∈ [[1, x 1 ]] do /*

294

We use the RSR (see Definition 4.2) exploring the property related to the distribution function of lead-times, and the 295 initial population is generated randomly with uniform distribution as follows:

296 X i,m = Rnd (T -U i,m , X * * i ) ∀i ∈ [[1, N m ]]
where Rnd (T -U i,m ; X * * i ) is a discrete random value between T -U i,m and X * * i .

In this suggested algorithm, we assume that x 5 individuals can be a local minimum. The procedure, introduced in [START_REF] Hnaien | Genetic algorithm for supply planning in two-level assembly systems with random lead times[END_REF] and applied one time to obtain the best neighbor solution (Best In Neighborhood), is considered.

Starting out from one individual, all neighbouring solutions are explored and the best one is selected to undergo the same operation again. This procedure is applied several times until a local minimum is met. Note that this procedure (Best In Neighborhood) is applied to the best solution at a given iteration if, and only if, (i) the best solution is not a local minimum, and (ii) x 7 iterations have passed without the best solution being improved.

Perturbation approach

The perturbation consists of replacing the solutions which have the same cost by solutions undergoing a special mutation (using block mutation). Each duplication X = (X 1,m , . . . , X i,m , . . . , X N m ,m ) undergoes a modification. To do so, an item (finished product or component) is selected randomly. The order release dates for components c v,m , . . . , c w,m at the last level m needed to assemble this item undergo a block mutation. Note that this approach is applied to duplications at a given iteration if, and only if, x 6 iterations have passed without the best solution being improved.

Computational experiments

In this section, we first present the data generation process, then we illustrate the main contribution of this paper and go on to describe the experiments carried out in order to study the effect of multiple factors on HGA performance.

We finish by analyzing the performance of the HGA.

Random test instances

We compared our HGA against a heuristic on 140 randomly generated test instances, defined in the following way.

First, we consider 14 different BOM, and randomly generate 10 instances for each BOM. For different BOMs, the total number of levels is equal to 8, and the total number of components (Comp total ) at all levels is equal to 39. We distinguish two types of BOM: B m and C m . In the intermediate levels 3 gives the BOM of each assembly system.

(l ∈ [[1, m -1]]) of B m , the number of components at each level N l is equal to 2, for C m , N l ≥ 2. Table
Each component lead-time L i, j , ∀ j ∈ [[1, m]], ∀i ∈ [[1, N j ]
] varies between two limits t i, j and u i, j with a discrete uniform probability distribution. These limits are generated randomly between 1 and 10 with t i, j < u i, j . The unit inventory

holding costs h i, j , ∀ j ∈ [[1, m]], ∀i ∈ [[1, N j ]]
for components c i, j and the unit inventory holding and backlogging costs for the finished product are generated as follows: (i) for components at the last level m, h i,m are random between 1 and 40; (ii) for components at level l ∈ [[1, m -1]], h i,l is random between α = ∑ c k,l ∈S i,l h k,l and 1.2α; (iii) for the finished product, the unit inventory holding cost r is random between α = ∑ c k,l ∈S FP h k,l and 1.2α; and (iv) the backlogging cost b is made to vary in {0.1r, 0.5r, r, 5r, 10r}. Compt total 39 39 39 39 39 39 39 39 39 39 39 39 39 39

B 2 C 2 B 3 C 3 B 4 C 4 B 5 C 5 B 6 C 6 B 7 C 7 B 8 C 8 N 1 2 4 2 2 2 4 2 2 2
Table 3: The BOM of each assembly system

Comparison between fitness functions

Before we go on to explain the theoretical basis of our HGA parameter setting, we first present the main contribution of this paper. In order to evaluate the efficiency of the proposed model, we carried out several experiments based on the test instances described above.

In this paper, the fitness function, which is an important part of the GA, is the expected total cost E(C(X, L)). A first explicit form of this expected value, noted V (.), was proposed in Ben-Ammar et al. ( 2018). Here we propose a new explicit form, noted V * (.) (see Propositions 4.1 and 4.2), which avoids using a recursive function, and can model real-world lead-times L i,l whose variations are logically between any two limits t i,l and u i,l .

We design a CGA to compute each instance 5 times using V (.) or V * (.) as the fitness function. Without preliminary testing or parameter calibration, the following parameters are fixed: maximum iteration number; noted x 1 ; is equal to 1000; population size x 2 is equal to 60 chromosomes; crossing-over probability x 3 is fixed to 0.95; and mutation probability x 4 is equal to 0.05. The CGA and the mathematical models, which calculate the fitness functions V (.) and V * (.), are coded in C++. The experiments were carried out on a computer with a 1.66 GHz Intel Core 2 Duo CPU and 4 GB of RAM. Figure 4 shows that the CGA using our new explicit form of E(C(X, L)) as fitness function is less sensitive to the increase in number of levels m than the CGA which uses the earlier form given in our last work. Thus, V * (.) will be used as the fitness function in the HGA.

Calibration of parameters and performance measures.

In order to find the set of HGA parameters that guarantees optimal performance for our algorithm, we decided to use a 'design of experiments' approach based on a standard central composite design (CDD). Based on a reduced number of experiments (see below for details), this kind of 3 level experimental design allows us to fit efficiently a quadratic multivariate regression model whose optimum can be searched for. The set of parameters being investigated is x 1 , x 2 , . . . , x 7 , which are set out in Section 5. The min, middle and max levels for these studied parameters are given in Table 4:

V (B) V (C) V * (B) V * (C)
Parameters Levels Min Middle Max
x 1 50 100 150

x 2 40 60 80

x 3 0.85 0.90 0.95

x 4 0.05 0.10 0.15

x 5 0 1 2

x 6 0 10 20

x 7 0 10 20 Note that the backlogging cost b is equal to r and a single set of instances was considered. There are 42 instances in each set (3B 2 , . . . , 3B 8 , 3C 2 , . . . , 3C 8 ). Using each combination of parameter values, the algorithm was applied in the following way. The set of instances was executed 5 times, and, for each execution, we determined the same single performance indicator i.e. the average expected total cost E(TC) of the best solution found. In order to obtain the p = 7 optimal parameters x * = (x * 1 , . . . , x * 7 ) of HGA for E(TC), we used a standard CCD with n = 79 experiments and k = 5 replicates, i.e. 395 observations. The CCD structure thus generated allows us to estimate the following quadratic model with 36 coefficients (a 0 , a 1 , . . . , a 35 ):

y = a 0 + a 1 x 1 + . . . + a 7 x 7 + a 8 x 1 x 2 + a 9 x 1 x 3 + . . . + a 28 x 6 x 7 + a 29 x 2 1 + . . . + a 35 x 2 7 + ε
where ε is an error term. Once this model was obtained, we searched for the solution x * minimizing this quadratic model in the hypercube defined by the min and max levels listed in Table 4. For E(TC), Table 5 (column headed 'coeff') gives the coefficients (a 0 , a 1 , . . . , a 35 ).

All the rows in Table 5 marked with 'a = 0' correspond to the regression coefficients ('coeff') a 0 , a 1 , . . . , a 35 for which the p-value (p-val) is larger than the standard threshold value α = 0.05(5%). As these coefficients are considered statistically equal to 0, their corresponding parameters or parameter combinations should therefore have no influence on response y. For example, concerning E(TC), parameters x 3 , x 4 and x 6 do not seem to be influential. The estimated model appears to fit the experiments well, as the coefficient of determination R 2 = 0.98216 is very close to 1. Thus, the vector of optimal parameters is x * = (150, 80, 0.9057, 0.141, 1.20, 11.87, 10.93). Knowing that x 5 , x 6 and x 7 must be integer values, they are then rounded to the nearest whole the number. Finally, x * = (150, 80, 0.9057, 0.141, 1, 12, 11).

Performance analysis of HGA

For each instance, 100 tests were carried out. The B&B method proposed in Ben-Ammar et al. ( 2018) requires more than 10 hours to find exact solutions. The CPU times depend not only on the number of components in the last level m, the number of levels m and the different costs, but are also significantly dependant on t i,l and u i,l , the lower and upper limits of the actual lead-times L i,l .

Note that the HGA is reinforced by several techniques, such as a reduced space of research (RSR), a local search (LS) and the integration of perturbation (P). We analyzed the influence of several parameters: (i) type of BOM; (ii) number of levels; and (iii) ratio of backlogging-to-inventory costs (b/r) for the finished product. To analyze the performance of the optimization approach, several notions were introduced:

• Average number of iterations in which the best solution is found; .100;

• Average CPU time of the HGA when the best solution is found.

Table 6 shows that the number of iterations needed to find the best known solution (BKS) is unaffected by the type of BOM and the number of levels. The efficiency of RSR translates through the fact that a big backlogging-to-inventory costs ratio (b/r) slightly reduces the average number of iterations needed to find the BKS. Note too that regardless of the variation of parameters, fewer than 60 iterations are required.

Table 7 shows the evolution of the average gap as a function of all parameters. This evolution becomes substantial in the case of either a small ratio b/r and/or small number of levels in the BOM. This is explained by the fact that the search space is biggest for a small b/r, and the number of components at the last level is big for the case of a small m.

Table 8 shows that, for assembly systems composed of 40 items, the HGA seems to be efficient and is not parameterdependent. Nevertheless, its performance needs to be further evaluated for more complex systems involving more than 40 items, and for a small backlogging cost for the finished product. Algorithm 2: Proposed heuristic

1 Classi f y(X i m ,m , ω i m ); 2 V 1 ←-/ 0; V 2 ←-/ 0; 3 foreach k ∈ [[1, N m ]] do 4 while E(C(Φ, L)) is decreasing & φ k ≤ X * * k m do 5 V 1 ←-E(C(Φ, L)); 6 φ k ←-φ k + 1; 7 end while 8 while E(C(Ψ, L)) is decreasing & ψ k ≥ T -U k m do 9 V 2 ←-E(C(Ψ, L)); 10 ψ k ←-ψ k -1; 11 end while 12 end foreach 13 UB -→ min(V 1 ,V 2 );
To analyze the performance of the optimization approaches, we introduce the following notions: .100;

• ACPU times (s), the average execution time of all CPU times (s) with the HGA when the best solution is found;

• ACPUH times (s), the average execution time of the proposed heuristic.

In Table 10, looking at the performances of HGA and the proposed heuristic, the b/r ratio seems to influence the quality of the best solutions and the CPU times. Nevertheless, the solutions proposed by the HGA are the best, they are still less than 1% from the best-known solutions found among the 100 tests, and less than 6% for the heuristic. For execution times, both approaches require less than a second (on average) to propose the best solution.

Conclusion and future research

In this work, a generalized probabilistic model, and algorithms to optimize the planned lead-times for multi-level ATO systems working in an MRP environment under uncertainty of actual lead-times are proposed. Planned lead-times were determined for the case of multi-level assembly systems with a one-period approach, i.e. for a given customer demand. An infinite assembly capacity at all levels was assumed, a known and fixed demand for the finished product, The same problem was treated in Ben-Ammar et al. ( 2018), where a recursive procedure was introduced to account for the dependence among levels, and a B&B algorithm was introduced to determine optimum solutions. However, the limitation of that study was that each lead-time varied between one and a known upper limit, and it was observed that the recursive function requires a recursively enumerable domain and depends on the distribution functions. The major drawback of this recursive procedure is its large influence on computation time, which becomes far too time-intensive.

In this paper, an iterative approach ensures a significant reduction in the time required to calculate the cumulative distribution function of each assembly date. To optimize the parameters, various additional techniques introduced in Ben-Ammar et al. ( 2016) were used. They are based on hybrid GAs (HGAs) to find good planned order release dates.

The results show that the HGA obtained can very efficiently find good-quality approximate solutions, regardless of the type of assembly system, the number of components at the last levels, or the variability of finished product-related costs related.

Our approach, for a given demand and due date, determines the optimal order release dates for the components, and is therefore readily adaptable to similar ATO environments. Furthermore, our analytical approach can be employed to generalize the models proposed in Ould-Louly et al. (2008b,a); Ould-Louly & Dolgui (2002a, 2013); [START_REF] Dolgui | Supply planning under uncertainties in MRP environments: A state of the art[END_REF]; [START_REF] Shojaie | A study on MRP with using leads time, order quality and service level over a single inventory[END_REF], and can also be applied in disassembly systems under uncertainties [START_REF] Bentaha | A sample average approximation method for disassembly line balancing problem under uncertainty[END_REF].

The current study was not designed to model multi-period planning for multi-level assembly systems. It also carries several limitations. First, it only considers an ATO environment with one demand. Second, we only consider release dates of the entire supply chain from the final assembly point of view, as we assume it is impossible to interact with supply chain partners to optimize all release decisions globally at all intermediate levels. Third, we do not explicitly treat the choice of suppliers. Fourth, we do not integrate supplier-related costs such as purchasing cost and ordering cost. Research is underway to overcome these limitations. We intend to focus on developing mathematical formulations for dynamic supplier selection strategies in multi-period supply planning for assembly systems under stochastic leadtimes. Furthermore, our future work will extend this model and the various techniques, in particular treating some new case studies calculating the planned lead-times when a company has to deal with production and replenishment time uncertainties.

This paper offers techniques for replenishment planning in dynamic ATO environments with stochastic lead-times.

discrete variables with a finite number of possible values. Knowing that:

P(M i,m-2 ≤ s) = P max c k,m-1 ∈S i,m-2 (X k,m-1 + L k,m-1 ≤ s)
and ∀k = 1, . . . , N m-1 , the random variables (M k,m-1 + L k,m-1 ), are independent, then: Using expression (15), we obtain the cumulative distribution function of M i,m-2 :

F i,m-2 (s) = ∏ c k,m-1 ∈S i,m-2 ∑ o 1 +w 1 =s o 1 +w 1 ∈N P(L k,m-1 = o 1 ) ∏ c k,m ∈S i,m-1 F k,m (-X k,m + w 1 ) (16)
At the third step, the iterative process starts by calculating F i,m-3 (s), the cumulative distribution function of M i,m-3 the assembly dates for components at level m -3. These assembly dates M i,m-3 are positive random discrete variables with a finite number of possible values. In the same way and using expression ( 16), we can deduce this cumulative distribution function: Proof of Proposition 4.2.

In Ben-Ammar et al. (2018), the general expression of the total expected cost, which is noted E(C(X, L)) was given and written as follows:

E(C(X, L)) = N 1 ∑ i=1 h i,1 E(M FP ) - m-1 ∑ l=1 N l ∑ i=1 H i • E(M i,l ) - m ∑ l=1 N l ∑ i=1 h i,l E(L i,l ) - N m ∑ i=1 h i,m E(X i,m ) + b • E(M + FP ) -T + r • T -E(M - FP ) (17)
Based on expression (2), the authors gave the expressions of E(M i,l ), E(M FP ), E(M + FP ) and E(M - FP ) (see Definition 4.1). They were calculated independently using the recursive function Q + (L i, j , s, j). Here, the cumulative distribution functions F . (.) of assembly dates (found by an iterative process; see Proposition 4.1) and the limits t i,l and u i,l of lead-times L i,l are used to simplify the calculations. β k,l+1 + u k,l+1 with β i,m = X i,m . Thus, M i,l varies between α i,l and β i,l (see Fig. 4), and expression (3) can be simplified as follows:

E M i,l = β i,l - β i,l -1 ∑ s=α i,l
Note that E M - FP is calculated in the same way and is equal to:

E M - FP = T - T -1 ∑ s=α 0,0 F 0,0 (s) (20)
where α 0,0 = max (α i,1 + t i,1 ).

The expression of E (M FP ) is calculated using expressions (19-20) and is equal to:

E (M FP ) = E M + FP + E M - FP -T = β 0,0 - β 0,0 -1 ∑ s=T F 0,0 (s) - T -1 ∑ s=α 0,0 F 0,0 (s) (21)
Then, using expressions (18, 19, 20 and 21), the mathematical expectation of the total cost can be found directly.

Figure 1 :

 1 Figure1: A multi-level assemble-to-order system.

Figure 2 :

 2 Figure 2: Composition of total cost (in the case of a backlog).

  ; Fallah-Jamshidi et al. (2011); Sakiani et al. (2012); Hnaien et al. (2016); Borodin et al. (2016); Guiras et al. (2019); Ben-Ammar & Dolgui (2018), this expression (2) is used to calculate E(M i,l ), E(M FP ), E(M + FP ) and E(M - FP ) for two-level assembly systems. In Ben-Ammar et al. (2018), the authors were the first to propose a recursive function that expressed the dependence among levels and enabled us to study assembly systems with more than two levels. Definition 4.1. (Ben-Ammar et al., 2018) Let Q + be the recursive function that serves to express the dependence among levels. All assembly dates M i,l , M - FP and M + FP are positive random discrete variables with a finite number of possible values. Their expected values are as follows:

  ; Fallah-Jamshidi et al. (2011); Sakiani et al. (2012); Hnaien et al. (2016); Borodin et al. (2016); Guiras et al. (2019); Ben-Ammar & Dolgui (2018), the authors assumed that all lead-times L i,l , ∀l ∈ [[1, 2]]

  Definition 4.2.(Ben-Ammar et al., 2018)

Algorithm 1 :

 1 Proposed genetic algorithm Function Best Subset(A, n): return S ⊆ A, S = n and s ∈ S, s ∈ A \ S, Fitness(s ) < Fitness(s); End Function card(S LO) ←x 5 ;

Fig. 3

 3 Fig. 3 illustrates the average CGA execution time after 1000 iterations. V (B) and V * (B) use the first type of BOM (B m ) and V (C), and V * (C) use the second type (C m ).Figure4shows that the CGA using our new explicit form of

Figure 3 :

 3 Figure 3: Comparison between fitness functions.

  and independent and identically distributed discrete component procurement times were considered. Therefore, the proposed analytical model calculates the mathematical expectation of the total cost. Our model is a generalization of those proposed inYano (1987a);[START_REF] Chu | Supply management in assembly systems[END_REF];[START_REF] Tang | The detailed coordination problem in a two-level assembly system with stochastic lead times[END_REF];[START_REF] Hnaien | Genetic algorithm for supply planning in two-level assembly systems with random lead times[END_REF][START_REF] Fallah-Jamshidi | A hybrid multi-objective genetic algorithm for planning order release date in two-level assembly system with random lead times[END_REF].

F

  k,m-1 = o 1 )P(M k,m-1 ≤ w 1 ) = ∏ c k,m-1 ∈S i,m-2 ∑ o 1 +w 1 =s o 1 +w 1 ∈N P(L k,m-1 = o 1 )F k,m-1 (w 1 )

F

  k,m-2 = o 2 )F i,m-2 (w 2 )Using the same iterative process, the cumulative distribution function of M i,l , for l ∈ [[0, m -4]] are calculated in the same way.

For

  l ∈ [[1, m]], each L k,l varies in [[t k,l , u k,l ]]. Let α i,l = max c k,l+1 ∈S i,l α k,l+1 + t k,l+1 with α i,m = X i,m and β i,l = max c k,l+1 ∈S i,l

F

  l ∈ [[α i,l , β i,l ]]

Figure 4 :

 4 Figure 4: Time interval in which the assembly date M i,l varies.By using the cumulative distribution function defined in expression (6), and knowing that M + FP varies between T andβ 0,0 = max k∈[[1,N m ]] X k,m +U k,m , E M +FP , given in expression (5), can be written as follows:

Table 2 :

 2 Model notation and definition

	b	Unit backlogging cost of the finished product per period
	r	Unit inventory holding cost for the finished product per period
	Variables
	X i,m	Decision variable: release date for component c i,m (this type of variable is only defined for components
		at level m), X i,m ∈ [T -U i,m ; T -T i,m ]
	Functions
	P(.)	Probability value
	E(.)	Expected value
	F(.)	Cumulative distribution function

Table 4 :

 4 Parameters and levels tested.

•

  Average gap between the best solution in the initial population and the best solution best sol 150 found by the algorithm after 150 iterations: gap = Average gap * between best sol 150 and the best-known solution best BKS found among the 100 tests: gap * =

	best sol 0 -best sol 150 best sol 150	.100;
	• best sol 1000 -best BKS
	best BKS

Table 5 :

 5 Standard CCD results for E(TC).

	m

Table 6 :

 6 Average number of iterations in which the best solution is found.

					m
	b/r 2	3	4	5	6	7	8
	0.1 5.77	4.96 5.36 3.44 2.96 3.39 0.45
	0.5 12.76 1.50 2.19 1.33 1.97 0.03 0.30
	1	0.00	0.16 0.51 0.09 1.61 0.09 0.08
	5	2.12	1.99 0.21 1.23 0.52 0.92 0.00
	10	1.40	0.70 0.05 0.00 0.00 0.01 0.11
	B
	25	3.52	0.41 0.19 0.00 0.02 0.00 0.00
	50	0.64	0.12 0.43 0.87 0.46 0.16 0.08
		0.00	0.17 0.12 0.00 0.33 0.11 0.17
		0.00	0.00 0.00 0.00 0.00 0.00 0.00
		0.00	0.00 0.19 0.36 0.09 0.00 0.00
	0.1 3.61	4.10 0.78 0.00 0.14 0.20 0.42
	0.5 9.25	3.81 1.78 2.48 6.83 0.09 0.01
	1	3.15	0.57 0.94 0.16 0.19 0.74 1.16
	5	0.95	0.00 0.17 0.00 0.00 0.00 0.00
	10	2.94	0.86 0.64 0.18 1.54 0.00 0.21
	C
	25	1.04	0.00 0.01 0.00 0.00 0.29 0.37
	50	0.00	0.00 0.02 0.03 0.00 0.02 0.06
		0.00	0.25 0.45 0.26 0.37 0.00 0.01
		0.00	0.00 0.08 0.00 0.00 0.16 0.41
		0.00	0.00 0.02 0.00 0.00 0.17 0.21

Table 7 :

 7 Average gap (%).

					m
	b/r	2	3	4	5	6	7	8
	0.07 0.14 0.29 0.34 0.34 0.66 0.42 0.45
	0.07 0.07 0.12 0.24 0.21 0.17 0.24 0.30
	0.01 0.02 0.08 0.15 0.19 0.28 0.20 0.08
	0.03 0.02 0.10 0.15 0.32 0.26 0.23 0.00
	0.02 0.03 0.04 0.04 0.08 0.11 0.35 0.11
	B
	0.02 0.03 0.09 0.06 0.11 0.08 0.16 0.00
	0.01 0.02 0.05 0.11 0.34 0.22 0.27 0.08
	0.02 0.04 0.06 0.04 0.08 0.15 0.16 0.17
	0.02 0.03 0.02 0.09 0.06 0.07 0.10 0.00
	0.03 0.03 0.02 0.05 0.17 0.06 0.22 0.00
	0.08 0.13 0.29 0.13 0.44 0.43 0.35 0.42
	0.08 0.19 0.41 0.43 0.66 0.32 0.39 0.01
	0.04 0.10 0.18 0.17 0.24 0.24 0.44 1.16
	0.02 0.02 0.11 0.11 0.09 0.15 0.18 0.00
	0.01 0.03 0.10 0.16 0.25 0.11 0.19 0.21
	C
	0.01 0.03 0.08 0.07 0.11 0.22 0.23 0.37
	0.02 0.02 0.05 0.08 0.08 0.11 0.23 0.06
	0.07 0.04 0.12 0.11 0.13 0.10 0.14 0.01
	0.07 0.01 0.05 0.06 0.08 0.14 0.38 0.41
	0.07 0.04 0.07 0.10 0.08 0.12 0.19 0.21

Table 9 :

 9 CPU time (s).

Table 10 :

 10 Performances of HGA and the proposed heuristic.

	b/r	0.1	0.5	1	5	10	25	50	100 200 500
	AGap * (%)	0.74 0.37 0.14 0.04 0.01 0.00 0.01 0.03 0.00 0.00
	AGapH * (%)	5.19 3.73 2.95 2.87 2.27 1.45 1.17 0.83 0.84 0.64
	ACPU times (s)	0.29 0.26 0.17 0.13 0.11 0.09 0.12 0.09 0.08 0.09
	ACPUH times (s) 0.03 0.03 0.03 0.03 0.03 0.02 0.02 0.02 0.02 0.02

1.21 0.00 0.14 1.73 0.86 0.89 0.16 0.45 0.66 1.26 0.12 0.00 0.33 0.00 0.00 0.30 0.00 0.16 0.00 0.01 1.61 0.07 0.03 0.08 0.00 0.00 0.00 0.09 0.45 0.00 0.00 0.00 0.00 0.05 0.02 0.00 0.00 0.00 0.00 0.11 0.00 0.00 0.00 0.00 0.02 0.00 0.00 0.00 0.00 0.00 0.15 0.00 0.00 0.00 0.00 0.08 0.00 0.00 0.02 0.00 0.15 0.00 0.06 0.17 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 C 3.60 1.26 0.14 0.00 0.02 0.00 0.33 0.42 2.51 0.00 0.24 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.07 0.00 0.00 0.00 0.00 1.16 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.05 0.00 0.00 0.00 0.00 0.21 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.37 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.06 0.00 0.00 0.18 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.03 0.00 0.00 0.00 0.02 0.41 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.21 Table 9 presents the execution times. Note that the best solutions are always found in less than a second and that the convergence to a good solution does not depend on number of levels, nor the type of the BOM, nor the backloggingto-inventory holding costs ratio.

For this problem, we used a heuristic and compared its performance against the HGA. It is the upper bound introduced in Ben-Ammar et al. (2018) and detailed in Algorithm 2. The solution proposed by this heuristic is equal to the minimum between two variables V 1 and V 2 . As mentioned in Definition 4.2, we decompose the multi-level assembly system to N m (the number of components at level m) multi-level linear supply chains. The different X i m ,m are ranked (Classify(X i m ,m , ω i m )) in descending order according to the costs of the linear chains. So, the first X 1 m ,m has the largest cost

We start by delaying the order release date φ 1 (by advancing ψ 1 ), and the same operation is executed until the E(C(Φ, L)) no longer decreases further.

We then repeat the same operations for order release date φ 2 of the next component.

To our knowledge, in only one previous study (Ben-Ammar et al., 2018) the authors develop a generalized probabilistic model to study the case of one-period planning for multi-level assembly systems. The newly-developed model and optimization algorithm can, in just a few seconds, compute optimal replenishment release dates for ATO supply networks with more than seven levels. All previous results have been limited to three levels. The proposed model can be used to reduce the impact of lead-time uncertainty in enterprises, especially by selecting appropriate planned lead-time parameters in their MRP systems. This is especially critical in manual assembly systems with several levels and components (typically automotive sub-modules or electrical appliances, for example). For small and medium sized problems, the models furnish optimal solutions. For more complex systems, the approach proposed in this paper can generate good-quality solutions within reasonable computation times.

From the practitioner's standpoint, the proposed approach can be successfully used for a number of applications.

For example, we have been working with ZF Friedrichshafen AG, a German car parts maker, at its assembly plant for gear boxes located in Saint-Étienne, France. Based on statistical data, the company allocated, for all suppliers, safety coefficients greater than 1 and used them in calculation of planned lead-times. A coefficient is calculated for each supplier according to its previous delivery performances. To set its planned lead-time in the MRP system, its contractual lead-time is multiplied by the corresponding coefficient. This strategy is applied to anticipate delays and to assess supplier reliability better. The closer the coefficient is to 1, the more reliable the supplier is. The main limitation of this empirical strategy is its inability to provide good parameters since suppliers are considered independently, the coefficients are calculated empirically and inventory synchronization aspects and costs are neglected. The proposed model and algorithms offer better estimations of these coefficients by taking into account inventory and backlogging costs, the independence (synchronization) of suppliers via the assembly operations, and probability distributions of supplier lead-times.

Appendix

Proof of Proposition 4.1.

Components c k,m at the last level m are ordered from several suppliers and used to assemble components c i,m-1 at level m -1. Then, at the first step, the iterative process starts by calculating F i,m-1 (s) the cumulative distribution functions of M i,m-1 , i.e. the assembly dates for components at level m -1. These assembly dates M i,m-1 are positive random discrete variables with a finite number of possible values. Knowing that L k,m + X k,m , for c k,m ∈ S i,m-1 , are independent, and that:

Then:

Then, at the second step, the iterative process starts by calculating F i,m-2 (s), the cumulative distribution function of M i,m-2 , i.e. the assembly dates for components at level m -2. These assembly dates M i,m-2 are positive random