Loïc Andolfatto
email: loic.andolfatto@eads.net

François Thiébaut

Claire Lartigue

Marc Douilly

Quality-and cost-driven assembly technique selection and geometrical tolerance allocation for mechanical structure assembly

Keywords: assembly process planning, assembly technique selection, geometrical tolerance allocation, multi-objective optimization, genetic algorithm

come

Introduction

The design of the entire manufacturing process of a product -and more particularly the assembly process -has a key impact on the performance of an industrial company. When it comes to mechanical structures composed of a high number of components, like aeronautical structures, choices made for the assembly process account for a large share in the total delivery cost and geometrical quality of the assembled products.

A complete assembly process plan is supposed to describe entirely how the product is assembled out of the given components. As pictured in Fig. 1, it includes the assembly sequence, the selection of the assembly techniques, the geometrical tolerance allocation on the component and the design and organisation of the assembly system. The assembly process planners have to reach several objectives generally expressed through various performance indicators.

Many efforts have been made to develop assisting solution for assembly process planning. Several methods have been proposed in the literature to evaluate (direct problem) or to optimize (inverse problem) assembly process plans. But these methods are generally limited to a single aspect of the general problem -such as sequence planning or tolerancing -and a single performance indicator.

This paper aims at proposing an original method to select assembly techniques and to allocate geometrical tolerances on components by solving a multi-objective optimization problem to minimize cost and maximize quality.

Tolerances on component geometrical deviations

Assembly techniques (and associated capabilities)

Forward problems

Inverse problem

Fig. 1: Evaluation scheme of an assembly process plan, from plan definition to performance indicators (with the proposed approach bordered).

Assembly sequence planning and organisation

Extensive work has been conducted to assist the generation of assembly sequences [START_REF] Wang | Assembly process planning and its future in collaborative manufacturing: a review[END_REF]. Bourjault [START_REF] Bourjault | Contribution a une approche methodologique de l'assemblage automatisé: Elaboration automatique des sequences operatoires[END_REF], De Fazio and Whitney [START_REF] De Fazio | Simplified generation of all mechanical assembly sequences[END_REF] and Dini et al. [START_REF] Dini | Generation of optimized assembly sequences using genetic algorithms[END_REF], among others [START_REF] Latombe | Assembly sequencing with toleranced parts[END_REF][START_REF] Marian | Assembly sequence planning and optimisation using genetic algorithms: Part i. automatic generation of feasible assembly sequences[END_REF][START_REF] Gottipolu | A simplified and efficient representation for evaluation and selection of assembly sequences[END_REF], presented the assembly sequence as the ordered list of components introduced into the assembly. Homem de Mello and Sanderson also considered the concept of attachment [START_REF] Homem De Mello | A correct and complete algorithm for the generation of mechanical assembly sequences[END_REF], what led to the design of assembly sequences defined as the order in which the product's links and joints are made [START_REF] Mantripragada | The datum flow chain: A systematic approach to assembly design and modeling[END_REF].

This joint-based approach reflects from some point of view the assembly task decomposition proposed by Cao and Sanderson [START_REF] Cao | Task decomposition and analysis of robotic assembly task plans using petri nets, Industrial Electronics[END_REF] which is itself close to the issues considered for the organisation of the assembly system commonly treated in operational research [START_REF] Becker | A survey on problems and methods in generalized assembly line balancing[END_REF].

Between those two problems, i.e. the generation of assembly sequences and the organisation of the assembly system, the issues for assembly technique 1 selection [12] are seldom addressed in the literature even though it may have the greatest impact on production cost, according to Abdullah et al. [START_REF] Abdullah | A review of the support tools for the process of assembly method selection and assembly planning[END_REF].

Almost all of the studies on assembly process planning aforementioned aim at minimizing lead time and/or cost, evaluated through various indicators such as, for example, tooling needs, reorientations of sub-assemblies, technological similarities in consecutive operations and so on [START_REF] Wang | Assembly process planning and its future in collaborative manufacturing: a review[END_REF].

Geometrical quality and tolerancing

The concern of geometrical quality of the assembled product is not commonly addressed by the assembly process planning community, even if the description of the geometrical variations of the components to be assembled may be part of the information contained in a comprehensive assembly process plan.

Nevertheless, the combined impact of the component geometrical variations, of the assembly technique capabilities and of the assembly sequence on the product geometrical quality was highlighted in several works related to tolerancing issues [START_REF] Mantripragada | The datum flow chain: A systematic approach to assembly design and modeling[END_REF][START_REF] Camelio | Modeling variation propagation of multi-station assembly systems with compliant parts[END_REF][START_REF] Mounaud | Assembly sequence influence on geometric deviations propagation of compliant parts[END_REF].

Chase et al. dealt with component manufacturing process selection to satisfy geometrical requirements on the assembled product [START_REF] Chase | Least cost tolerance allocation for mechanical assemblies with automated process selection[END_REF]. Adragna et al. proposed a tolerance allocation method that maximises the assembly process capability index [START_REF] Adragna | Inertial tolerancing and capability indices in an assembly production[END_REF]. Ding et al. and Huang et al. presented process-oriented method for tolerance synthesis in multi-station manufacturing environment [START_REF] Ding | Process-oriented tolerancing for multi-station assembly systems[END_REF][START_REF] Huang | Process capability sur-rogate model-based tolerance synthesis for multi-station manufacturing systems[END_REF].

These studies aim at allocating geometrical tolerances to satisfy objectives on indicators of the geometrical quality of the assembled product.

Proposed approach

In the field of aeronautical structure assembly, making tradesoff between automated and manual assembly plans proved to be very complex. The need for decreasing manufacturing cost and increasing delivery rate is usually solved by massive automation in other manufacturing domains, such as automotive industry. But this solution is not always compatible with the high level of geometrical requirements on aeronautical structures.

This paper focuses on a method to select assembly techniques and allocate component geometrical tolerances in order to minimize a cost indicator and to maximize a quality indicator associated with the assembly plan. Considering the assembly technique selection together with the geometrical tolerance allocation allows exploring a wide range of potential solutions.

This multi-objective optimization approach prevents from making decision a priori while it let the decision-making team select the most appropriate assembly plan among several optimal ones. 1 Also called assembly process by some authors.

The assembly sequence (attachment-based defined) is considered to be predetermined -according to the method proposed in [START_REF] Marguet | Contribution à l'analyse des variations géométriques dans les ensembles structuraux en aéronautique: Démarche et outils[END_REF] for example -. The detailed design and organisation of the assembly system take place once the assembly techniques are chosen and are thus not considered in this study. The resulting boundary of the study is illustrated in Fig. 1.

The second section details the data structure used to define a parametric assembly plan, what is required to tackle its optimization. The multi-objective optimization set to solve the problem addressed in this paper is described in the third section.

The section four details the proposed evaluation of the assembly plan. The quality indicator is a conformity rate evaluated thanks to a probabilistic study based on the geometrical variation propagation relation associated with the assembly process plan. The evaluation of the assembly cost indicator combines operations and tolerances cost. The former relies on a simple activity-based analytical model applied on each operation. The latter is obtained thanks to a cost versus tolerance relation adapted from [START_REF] Chase | Least cost tolerance allocation for mechanical assemblies with automated process selection[END_REF]. Cost and quality indicators' evaluation are both intentionally simple to let the reader focus on the global approach proposed for technique selection and tolerance allocation. These indicators can easily be replaced by more realistic ones for actual industrial applications.

The optimization method used is briefly described in section five.

The entire method is illustrated through a simple case study presented in section six where results are exposed and discussed. It is followed by a conclusion.

Parametric assembly plan

Introduction

The inverse problem displayed in Fig. 1 is solved thanks to a multi-objective optimization. Therefore, the performance indicators of the assembly plan -the cost and the conformity rate in this study -must be expressed as mathematical functions (described in section 3). The input parameters of those two functions must be parameters that describe the assembly plan to evaluate. This section explains how this set of input parameters is extracted from the description of the product and from the technical know-how of a company.

Structuro-functional model of a product

An assembled product is a set of components connected to each other thanks to assembly joints. Each joint involves two surfaces, each one belonging to one of the components. The structure of the product can be represented by an Elementary Contact Graph [START_REF] Ballu | Choice of functional specifications using graphs within the framework of education[END_REF] (ECG). The joints can be classified into two categories: mates that pass dimensional constraints from part to part, and contacts that provide support and lead to hyperstatic assembly [START_REF] Mantripragada | The datum flow chain: A systematic approach to assembly design and modeling[END_REF]. This classification helps building Datum Flow Chain (DFC) to identify components and joints involved in the variation of the product's key characteristics (KC) [START_REF] Mantripragada | The datum flow chain: A systematic approach to assembly design and modeling[END_REF]. A key characteristic is a property of a product required to satisfy a function. In this paper, we only consider KC expressed as geometrical conditions between two surfaces, even if the definition applies for any measurable characteristic. Marguet proposed to represent components, joints (mates and contacts) and KC thanks to an Oriented Contact Graph that mixes Elementary Contact Graph and Datum Flow Chain approaches [START_REF] Marguet | Contribution à l'analyse des variations géométriques dans les ensembles structuraux en aéronautique: Démarche et outils[END_REF]. This graph represents the structuro-functional model of a product. It was implemented in a software called GAIA R [START_REF] Falgarone | An innovative design method and tool for structural and functional analysis[END_REF].

The Oriented Contact Graph can also represent temporary components used during the assembly, like tools and jigs. Joints between temporary components and components of the product are called temporary joints. They are all released at the end of the assembly.

An example is given in Fig. 2. Two components are positioned by a tool (tj1 and tj2). The joint J1 is made by drilling the b surfaces located with respect to the tool. The key characteristic KC1 is not impacted by the b surfaces location as illustrated by the DFC.

Finally the structuro-functional model of a product expressed through the Oriented Contact Graph allows one to list the J joints of the product. An assembly technique has to be chosen for each one. It also allows one to list the L links of the datum flow chains (i.e. the component dimensions and the joints involved in the variations of the product's Key Characteristics). A tolerance has to be allocated to each component dimension to ensure that the product's KC are kept within the authorized boundaries.

Assembly sequence representation

As mentioned in subsection 1.1, the assembly sequence can be described either as the order for component insertion or as the order for joint manufacturing. The latter method includes somehow the former, because a joint can only be made if the components involved are already inserted in the assembly. In order to provide a clear representation of the sequence, both approaches are combined in an Assembly Sequence Graph. An example is displayed in Fig. 3. The Assembly Sequence Graph is dual to the Oriented Contact Graph. It shows the order in which components are inserted (and removed when temporary) in the assembly and when the joints are made (and released when temporary).

The assembly sequence is considered as an input for the work presented in this paper. The assembly sequence can be defined by applying a method based on minimizing the length of the datum flow chains as proposed by Marguet [20] for instance.

Assembly technique definition

This paper aims at proposing a method to select an assembly technique to make each joint of the product. For instance, the same riveted joint can be made with manual drilling and rivet fastening or with robots fitted with drilling and fastening endeffectors. These are two possible assembly techniques to make the same joint. Depending on the product, it can be hard to find which one is the most suitable.

An assembly technique is defined as a collection of assembly operations at the end of which a joint is made. For instance, the manual drilling and rivet fastening technique previously introduced require manual drilling, cleaning, sealing, rivet inserting and rivet fastening.

Operation characteristics used in this study are fixed cost, duration, required resources and resource quantities. If an operation is prone to geometrical deviations impacting the product's KC, the probability distribution characterizing its capability is also required in the operation definition. Table 1 gives an example of an assembly technique to make a temporary joint between a tool and a component. The manufacturing company know-how can be formalized by building an assembly technique library. The assembly process planner must define the list of techniques he wants to be explored for each joint of the product instead of manually building and evaluating assembly plans according to the company knowhow not necessarily formalised. It is among this list that an optimal solution will be searched for. Then the search is likely to be more efficient with respect to the high combinatorial that the assembly process planner have to face.

Optimization problem

Decision variables and constraints

An assembly technique has to be chosen for each joint of the product (including the temporary joints represented on the Oriented Contact Graph). Let N t be the number of techniques in the assembly technique library. The list of authorized assembly techniques for the joint j ∈ {1, . . . , J} is denoted T j and is a subset of {1, . . . , N t }. The assembly technique decision variables are gathered in a vector x t defined by eq. ( 1):

x t = x 1 . . . x J T (1)
The assembly technique decision space is then limited to Ω t , defined in eq. ( 2):

Ω t = J j=1 T j (2)
The DFC identified from the sructuro-functional model of the product allow ones to list the L links that impact the product KC. These L links geometrical variations must be allocated tolerances. L t links of the DFCs are already characterized by the assembly technique selected for their related joint, as for tj1 and tj2 in Fig. 2. Geometrical tolerances must be allocated to the remaining L g = L -L t links involved in the variation of the product KC.

Let i ∈ {1, . . . , L g } be the index of one of these links. Its probabilistic distribution is assumed to be fixed. The distribution has n i parameters written {x i1 , . . . , x ini }. For example, a normal distribution has two parameters -the mean value and the standard deviation -and n i equals two. Allocating tolerance comes to setting these two (or n i ) parameters, called geometrical variation decision variable. The geometrical variation decision variables are gathered in a vector x g defined in eq. ( 3):

x g = . . . x i1 . . . x ini . . . T (3)
The dimension of the geometrical variation decision space is denoted N g . Its value is calculated according to eq. ( 4):

N g = Lg i=1 n i (4)
Some constraints can exist in each subset of the geometrical variation decision variables {x i1 , . . . , x ini }. For instance, if the distribution i is uniform, then the lower bound x i1 must be lower than the upper bound x i2 . Those constraints are intrinsic to the distributions. The geometrical variation decision space noted Ω g is the restriction of R Ng satisfying to all the distribution intrinsic constraints.

Finally, the assembly technique selection and component geometrical deviation tolerances can be expressed by n = J + N g decision variables gathered in a vector x defined in eq. ( 5):

x = x t x g T ∈ Ω (5)
with:

Ω = Ω t × Ω g (6)

Multi-objective optimization

As illustrated in Fig. 1, the definition of an assembly process plan is driven by several performance indicators. This work deals with the search for assembly plans that minimize assembly cost indicator and maximize quality indicator (i.e. minimize non-conformity rate defined in subsection 4.1). These objectives are likely to conflict with each other.

One can decide to express the problem into a single objective optimization problem, combining both performance indicators into a single objective function or optimizing a single indicator handling the other one as a constraint. Although this approach seems simpler from a mathematical point of view, its main drawback is that it forces the user to introduce high-level information into the model used for optimization. The tradeoff between several performance indicators can not always be encapsulated into a unique equation.

Solving a multi-objective optimization problem seems to be more adapted within the industrial context. The optimization result is a set of non-dominated solutions (see Fig. 4 for illustration) obtained without modelling anything but the performance indicators. A decision variable vector x * dominates a decision variable vector x according to the N objective functions {f 1 , . . . , f N } if the assumptions of eq. ( 7) are verified [START_REF] Edgeworth | Mathematical physics; an essay on the application of mathematics to the moral sciences[END_REF][START_REF] Pareto | [END_REF]:

f 1 f 2 f 1 (x * ) f 2 (x * )
∀p ∈ {1, . . . , N }, f p (x * ) ≤ f p (x) ∃q ∈ {1, . . . , N }, f q (x * ) < f q (x) (7)
When the set of non-dominated solutions, commonly called the Pareto front, is found, it can be used as a decision aid. Then the final decision on assembly techniques and geometrical tolerance allocation is taken by a decision-making team among the non-dominated solutions. This approach fits common industrial practices better than a blind selection according to a black-box optimization result.

The multi-objective optimization problem dedicated to assembly technique selection and geometrical tolerance allocation is finally formulated in eq. ( 8):

Minimize x∈Ω [NCR(x), C(x)] T (8)
The non-conformity rate function NCR(x) : Ω → R is detailed in subsection 4.1 and the assembly cost function C(x) : Ω → R is detailed in subsection 4.2.

Assembly plan evaluation

Evaluation of the non-conformity rate

The quality indicator used in this paper is called non-conformity rate. The non-conformity rate associated with an assembly process plan figures the probability for a product to have one of its KC out of the boundaries derived from the functional requirements.

The structuro-functional model of the product presented in subsection 2.2 allows one to identify the L links of the DFC. The geometrical variation propagation relation associated with the assembly process plan is assumed to be known. Each KC of the product is expressed as a function of the DFC's links denoted δl i in eq. ( 9):

δKC = f {δl i } i∈{1,...,L} (9)
The geometrical variation propagation relation generically expressed by eq. ( 9) can either be a simple geometrical variation stack-up, or a linear relation integrating component flexibility as in the Method of Influence Coefficient [START_REF] Liu | A parametric study of joint performance in sheet metal assembly[END_REF]. The relation f can even be expressed by a meta-model, such as a trained neural network, that accounts for the non-linearity due to the assembly sequence and the contact interaction between components [START_REF] Andolfatto | On neural network's ability to approximate the geometrical variation propagation in assembly[END_REF].

The probability density function (PDF) of each δKC can be evaluated according to the geometrical variation propagation relation given in eq. ( 9) fed with the probabilistic distribution of the δl i . The distributions of the δl i are either given by the assembly technique definitions related to x t or by the geometrical variation parameters defined in x g .

The probability for the KC deviation indexed r to be smaller than its functional lower bound lb r or greater than its functional upper bound ub r can then be extracted from its PDF. The non-conformity rate of the key characteristic r is expressed in eq. ( 10):

NCR r (x) = P[δKC r < lb r ](x) + P[δKC r > ub r ](x) (10)
The probability for the KC to exceed its bounds is evaluated applying a First Order Reliability Method (FORM) [START_REF]TURNS -Reference Guide -version 1.0[END_REF]. Various input distributions can be considered with this method, making it more generic than approaches such as Root Sum Square. This method can also be more computationally efficient than Monte Carlo based methods when it comes to predict low probabilities.

The non-conformity rate associated with the assembly process plan is defined as the maximum of the non-conformity rate of one of the N KC of the product, as displayed in eq. [START_REF] Becker | A survey on problems and methods in generalized assembly line balancing[END_REF]:

NCR(x) = max r∈{1,...,N KC } NCR r (x) (11)
For example, considering the assembly of a product including two KCs with respective probability 0.1 and 0.2 to be out their functional domain, the associated non-conformity rate will be 0.2. This approach can provide a non-conformity rate lower than the total rate of non-conform products but it allows a separate analysis for each KC.

Evaluation of the cost 4.2.1. Cost sources

The cost associated with an assembly plan can be split into two sources:

• the activity cost, denoted C a , can be evaluated by applying an activity based cost model for all the operations of the process plan;

• the geometrical tolerance cost, written C g , and depending on the geometrical tolerance allocation.

The choice of the assembly techniques directly impacts the activity cost of the assembly process plan. It also indirectly impacts the geometrical tolerance cost as the tolerances allocated may be adjusted for each set of selected assembly techniques. The two cost sources are treated separately in the following but are considered together to evaluate the total cost associated with an assembly process plan:

C(x) = C a (x t ) + C g (x g ) (12)

Activity based cost model

The evaluation of the assembly operations's cost is performed using a generic model applied to each operation. Even if this model is simple, it can easily be replaced by a more realistic one if required.

Let Op α be the operation α of the process plan. Let t α be its duration and F c α its fixed cost. The fixed cost represents, for instance, consumables or tool wear. This operation also consumes the quantities m αβ of the resources {r β } β∈{1,...,R} . Every resource r β has a cost per time unit C tu,β . Then the elementary cost C Op,α associated with the operation Op α is the sum of the fixed cost and the cost of the resources consumed, as expressed in eq.( 13):

C Op,α = F c α + t α × R β=1 m αβ × C tu,β (13)
The assembly recurring cost C req (i.e. the cost required to assemble each product) is assumed to be the sum of the elementary cost of each of the N op operations {Op α } α∈{1,...,Nop} of the assembly plan weighted by the number of occurrences q α of each operation in the plan. This is expressed in eq.( 14) with a matrix notation.

C req = q T • C op (14)
where:

q = q 1 • • • q Nop T (15)
and:

C op =    F c 1 . . . F c Nop    +    t 1 • • • 0 . . . . . . . . . 0 • • • t Nop    •    m 11 • • • m 1R . . . . . . . . . m Nop1 • • • m NopR    •    C tu,1 . . . C tu,R    (16)
With this formulation, the recurring cost is the product of q T , that depends on the choice of assembly techniques x t only, and C op which is a constant vector depending on the assembly technique library.

When it comes to assembly technique selection, not only the recurring cost has to be taken into account, the non-recurring cost must also be evaluated. Let C r,β be the investment to acquire the resource r β . It is assumed that this investment cost is not included in the resource cost per time unit. With q r,β the quantity of the resource r β required by the assembly process plan, then the non-recurring cost C nr associated with the assembly process plan is expressed thanks to eq.( 17):

C nr = q r,1 • • • q r,R •    C r,1 . . . C r,R    = q T r • C r (17)
The non-recurring cost is then the product of q T r , that depends on the choice of assembly techniques only, with C r a constant vector. When a product is manufactured at an industrial scale, there is often a potential market identified for the product and the assembly process can be planned considering the total amount P of product that may be sold. Based on this hypothesis, the non-recurring cost associated with the assembly process plan can be shared among those P products.

Finally, the cost C a associated with the assembly operations of a process plan is expressed in eq.( 18):

C a (x t ) = q(x t ) T • C op + 1 P • q r (x t ) T • C r (18)
The evaluation of the C a associated with an assembly process plan is reduced to the evaluation of q and q r . Resource dimensioning has been investigated in operational research [START_REF] Tan | Integration of process planning and scheduling a review[END_REF] and many possibilities exist to evaluate q r . A simple method is applied: the quantity of a resource is equal to the maximum quantity consumed by an operation of the plan, as expressed in eq.( 19):

∀β ∈ {1, . . . , R}, q r,β = max α∈{1,...,Nop},qα =0 m αβ (19)
For monotonous assembly processes, q can directly be extracted according to the choice of assembly technique for each joint of the product. If some operations of a technique must be repeated several times due to constraints coming from the manufacturing of another related joint -what is a common case in aeronautical structures assembly -, finding q according to the assembly technique choice may require a more elaborated method, like finite-state-machine or Petri net models of the assembly process. In this paper, only monotonous assembly processes are considered.

Cost associated with the geometrical tolerance

The cost associated with the geometrical tolerance allocation of the process plan is assumed to be the sum of the cost required to keep each of the L g dimensions in its allocated tolerance zone.

The relation between the size of the tolerance allocation and the associated cost is commonly investigated in the literature. Chase et al. gathered several simple cost vs. tolerance models in [START_REF] Chase | Least cost tolerance allocation for mechanical assemblies with automated process selection[END_REF]. Most of those cost models can be generalized by the eq. ( 20) in which a, b, m, k and T lim are fixed parameters of the relation.

c = a + b • e -m(T -T lim ) •(T -T lim ) -k (20)
Therefore, for each i ∈ {1, . . . , L g }, a value for T i must be calculated according to the {x i1 , . . . , x ini } parameters of the probabilistic distribution of the dimension i:

• if its probability density function (PDF) is not equal to zero on a finite interval of R (like uniform or truncated normal distribution), T i is the size of this interval;

• else, the value used for T i is chosen equal to the size of the 99.73% coverage interval (6σ for a normal distribution).

The cost associated with the geometrical tolerance allocation is finally expressed in eq. ( 21):

C g (x g ) = Lg i=1 a i + b i • e -mi(Ti-T lim,i ) • (T i -T lim,i ) -ki (21)
Some more elaborated models can easily be included in the framework presented in this paper.

Optimization solution

Choice of an optimization method

Most of the existing methods to solve multi-objective optimization problems rely on population-based evolutionary algorithm. Coello Coello and Reyes Sierra investigated pros and cons several optimization solutions [START_REF] Coello Coello | Recent trends in evolutionary multiobjective optimization[END_REF][START_REF] Reyes-Sierra | Multi-objective particle swarm optimizers: A survey of the state-of-the-art[END_REF]. As it provides a good trade-off between performance and robustness while requiring few parameters to tuned, the NSGA-II algorithm (Non-Sorted Genetic Algorithm -II) proposed by Deb et al. [START_REF] Deb | A Fast Elitist Nondominated Sorting Genetic Algorithm for Multi-objective Optimization: NSGA-II[END_REF] is chosen to solve the optimization problem formulated in eq. ( 8).

Implementation

Fig. 5 presents an overall flowchart of the search for nondominated points in the non-conformity rate versus cost plane.

An initial population of N p assembly plans modelled by N p x decision variable vectors is created randomly. The fitness of each element of the population is evaluated to create N p pairs (N CR, C) according to eqs. 11 and 12. The nondominated points are stored in an archive. Then, based on the selection described in the NSGA-II algorithm, the population evolves thanks to two separate mechanisms.

Crossovers are performed between two parents to create an offspring. The crossover operator used in this paper is called blend crossover. It consists in creating an offspring by computing the arithmetic mean of two parents x p1 and x p2 weighted respectively by α and 1 -α, with α randomly chosen between 0 and 1. This strategy is chosen due to its ability to cope with

i = 1 i < N p ? i NCR = NCR( i x) (11) i C a = C a ( i x t ) ( 18
) i C = i C a + i C g (12) i C g = C g ( i x g ) (21) i = i+1 g < N gen ? g = g+1
Search nondominated points (7)

Store in archive

Load archive

Generate initial population { i x} with N p individuals

Search nondominated points (7)

Return non-dominated points

Create new generation: selection, crossover, mutation continuous variables such as the geometrical part x g of the decision variables vector x. Some candidates x are also subject to Gaussian mutation. It consists in adding a random value drawn with Gaussian distribution to each component of the candidate x. That provides the algorithm with enhanced exploration capability, while crossover alone generally leads to a narrow part of the actual Pareto front -a small set of assembly techniques explored in this study for instance -. A mutation rate r m is defined, and the number of candidates subject to mutation in each generation is r m Ṅp .

When an arbitrarily chosen number N gen of generations has been created and evaluated, the evolution stops and non- The method presented is implemented in Python language. The optimization is performed with the NSGA-II algorithm implemented in the Inspyred Python library [START_REF]Inspyred 1.0: Bio-inspired Algorithms in Python[END_REF].

Use case and results

Structuro-functional model of the product

The methodology described in sections 2 to 5 is illustrated through a simple case study. The product is a box composed of four components: two flanks (right and left), one splice and one cover. Fig. 6 shows a schematic view of this product. The splice is bonded on the two flanks and the cover is riveted on them.

The box has six key characteristics: the bonding areas of the splice (KC1 and KC2), the outcrop of the cover (KC3 and KC4) and the location of the rivets (KC5 and KC6), as pictured in Fig. 6. Fig. 7 presents the assembly sequence and the structurofunctional model of the product after the ninth assembly step. Fig. 8 shows the product state after the sixteenth step of the assembly sequence. Details about the names of the surfaces are given in Fig. 9.

Assembly techniques libraries

The assembly technique library gathering all the assembly techniques available in this study is detailed in Appendix B.

As explained in subsection 3.1, each joint j has its own authorized assembly technique list T j . The Table 2 summarizes the authorized assembly techniques defined for each joint of the product. The temporary joints tj 1 and tj 2 are made by back-toback positioning on a reference plane. So are J 1 -m, J 3 and J 5 , with a few modifications detailed in Table B.1. The temporary joints tj 3 to tj 7 can be made either with a tool, an adjustable tool, a robot used as a tool or a specific automated station. Details are given in Table B.2. Two bondings techniques are available for joints J 1 -c and J 2 . Differences between the two are

Joint j

Authorized assembly techniques list T j tj 1 and tj 2 1-Back-to-back positioning (temporary) tj 3 to tj 7 2-Tool, 3-Adjustable tool, 4-Robot as tool, 5-Specific automated station J 1 -m 6-Back-to-back positioning J 1 -c and J 2 7-Traditional bonding, 8-Rapid bonding J 3 and J 5 9-Back-to-back positioning (with splitting) J 4 and J 6 10-Drilling with grid, 11-Drilling with numerical grid, 12-Robot Drilling, 13-Specific drilling station given in Table B.1. Finally, the joints J 4 to J 6 can be made with either a traditional drilling grid, a numerical grid, a robot fitted with a drilling end-effector or a specific drilling station, as exposed in Table B.3. The assembly technique decision variable vector x t (see eq. ( 1)) has 14 components representing the index of the technique chosen for each of the 14 joints of the product. Some assembly technique decision variables are constrained to be equal because of technical constraints. It is the case for the technique selected for the sets of joints {tj 3 , tj 4 , tj 5 }, {tj 6 , tj 7 } and {J 4 , J 6 }.

Even with these constraints set, the domain to explore is still wide, from manual to highly automated assembly techniques. The aim of the study is to search optimal solution without making any more decision.

Datum flow chains and conformity rate

The datum flow chain of each key characteristic can be built according to the structuro-functional model of the product and its assembly sequence defined in Fig. 7. The datum flow chain of KC3 is displayed in Fig. 8. The geometrical variation propagation problem is kept unidimensional in this use case. The value of KC3 depends on the location of the contact in joints tj6 and tj7. Then, the deviation δKC3 is calculated according to eq. ( 22):

δKC 3 = -δtj 6 + δtj 7 (22)
where δtj 6 and δtj 7 stand for the deviation of the contact locations with respect to their theoretical locations along the x axis.

In the following, δl iab denotes the geometrical deviation of the surface b with respect to the surface a of the component i. The geometrical deviation due to the joint j is written δj. Then, the eq. ( 9) applied to each datum flow chain of the product leads to the system of eqs. [START_REF] Edgeworth | Mathematical physics; an essay on the application of mathematics to the moral sciences[END_REF]: The geometrical deviations of the components impacting the product key characteristics are listed according to eq. ( 23): it consists in δl 1ab , δl 2ab , δl 3ab and δl 4ab . Tolerances are allocated for those deviations in the following. They are assumed to have uniform distributions. Then, the lower and upper bounds of the allowed deviations need to be set. The geometrical variation decision variable vector x g (see eq. ( 3)) has 8 components representing the lower bound and upper bound of the tolerance allocated for each of the four component characteristics impacting the product KC and therefore the conformity rate. Table 3 summarizes their cost vs. tolerance model parameters from eq. ( 21) and Fig. 10 gives a graphical representation of the cost vs. tolerance relation. The geometrical variations of l 1ab and l 2ab are allocated the same tolerance due to the symmetry of the components.

               δKC 1 = δtj 3 -δtj 5 + δl 1ab δKC 2 = -δtj 4 + δtj 5 -δl 2ab + δl 3ab δKC 3 = -δtj 6 + δtj 7 δKC 4 = -δtj 3 + δtj 4 + δtj 6 -δtj 7 -δl 4ab δKC 5 = -δtj 6 + δmJ 4 δKC 6 = -δtj 3 + δtj 4 + δtj 6 -δmJ 6 (23)
The joints subject to deviations impacting the product key characteristics are tj 3 to tj 7 . The surfaces e of the flanks and e and f of the cover are manufactured when the joints J4 and J6 are made. The deviations δmJ 4 and δmJ 6 , that depend on the technique selected for J4 and J6, also impact the key characteristics and the conformity rate.

Implementation

The method presented in this paper is implemented in Python language. The optimization is performed with a population composed of 200 individuals evolving during 20 generations. The size of the population is selected to provide a good trade-off between spread and speed. The number of generations has been set according to empirical results showing that 3. the Pareto front does not significantly moves during further evolutions. Mutation rate is set to set at 0.1 and blend crossover is applied to the entire population. Those settings are providing a good spread along the Pareto front while keeping the computation time below one hour2 . 6.5. Results

Main results

The optimization problem presented in the previous sections is solved considering the number P of products to be assembled equal to 100.

The Pareto front obtained at the end of the evolution of 200 decision variable vector along 20 generations is presented in Fig. 11. Each point of the graph represents the non-conformity rate and cost associated with a set of techniques chosen for the product's joints and a set of tolerances allocated to the components' geometrical deviations. The values in the Table 4 are the U(-0.17, 0.16) decision variables corresponding to the point represented by a blue square in the Fig. 11.

The benefit of the multi-objective optimization can easily be illustrated in Fig. 11. If the problem was solved considering the cost as a single objective and the non-conformity rate constrained to be zero, the result would have been the point in the top left hand corner (cost equal to 4412 cost units). Instead of this, the non-conformity rate could have been constrained to be lower than an arbitrary value decided before running the optimization. Let this value be 5% for the example. It leads to a cost of around 2620 cost units. But running a single objective optimization do not allow the user to know that he can decrease the non-conformity below 2% for around 2700 cost unit per product. The main interest of the multi-objective optimization is that it helps the user decide which solution fits its real objectives the best instead of needing for him to translate his real objectives into blind mathematical constraints.

The Pareto front alone does not provide all the information obtained at the end of the optimization. Four zones are circled in Fig. 11. In each zone, the assembly technique decision variables of the points are identical. In fact, if the assembly tech-niques were all selected before and the optimization was only performed to allocate tolerances, the Pareto front would have a smoothly decreasing shape. Here, each zone corresponds to a smooth front with its own set of assembly techniques selected. In other words, the multi-objective optimization approach shows its ability to let the assembly process planner identify which set of assembly techniques seems to be the most relevant.

It is emphasized in

For all zones, traditional bonding and drilling with grid proved to be the most relevant techniques.

In zone 1, the joints tj 3 to tj 7 are made with the Adjustable tool technique. As expected, the low geometrical variations introduced by assembly techniques helps reaching a non-conformity near zero but the measurement and setting operations make the cost grow noticeably.

In zone 2, the joints tj 3 to tj 5 are made with the Adjustable tool technique and the joints tj 6 and tj 7 are made with the Robot as tool technique. Compared to zone 1, the cost can be significantly reduced if a small loss on the non-conformity is permitted.

In zone 3, only the Robot as tool technique is used. Compared to zone 2, the cost can also be reduced but its reduction will be associated with more significant loss on the nonconformity. This part of the Pareto front provides a lot of information for potential trades-off.

In zone 4, the joints are made with the Tool technique. According to the level of non-conformity reached, this choice of techniques is probably not relevant.

With these results, the assembly process planner can focus on a small amount of solutions. For example, the search for optimal solution can be refined in the zone that best fits its objectives. The actual cost of repairing and/or adjusting products that are not meeting the functional requirements can also been investigated.

The combination of activity-based cost and geometrical tolerance cost into a single indicator can be a challenging work in practical industrial conditions. The identification of several zones in the Pareto front can help overcome this hurdle by damping the effect of an indicator not well balanced between the two cost sources.

The combination of the two cost sources can also be avoided by considering them separately and solving a three-objective optimization problem.

The results presented in Fig. 11 also illustrate the capability of the method in terms of search. The Pareto front obtained is wide enough to show several assembly technique sets with various scenarios of geometrical tolerance allocation. This shows that the optimization algorithm has explored the entire domain.

Illustration with P = 10000

The optimization problem was also solved with a number of products P = 10000 to present the results in a case of manufacturing on a larger scale. The Pareto front obtained is presented in Fig. 12. The zone 4 (in which the joints are realized with the Tool technique) is no longer present in the Pareto front because the investment required to use the Robot as tool technique is compensated by a lower assembly cost.

The separation between zones is less visible than in Fig. 11 but the Pareto front globally has the same shape and the conclusions drawn in the previous subsection still apply.

Conclusion

This paper aims at presenting an original method for computeraided assembly process planning. It addresses the problem of assembly technique selection and of component geometrical tolerance allocation. This problem is solved to satisfy both objectives of minimizing a cost indicator and maximizing a quality indicator. This is achieved by writing a multi-objective optimization problem.

The optimization result, presented as a Pareto front, proved to be an efficient tool for making trades-off once a set of nondominated assembly plans have been found. The identification of zones on the Pareto front helps the assembly process planner to identify relevant sets of assembly techniques. Therefore, the proposed approach is well adapted to an industrial context where decision-aid capabilities are essential.

The method can be used in different ways to satisfy industrial needs. On top of the results illustrated in this paper, it can help evaluate and assess the potential benefits of using innovative assembly techniques or help find the minimum capability required for a technique under development to be competitive, for example.

The use case and the quality and cost indicators presented in this paper are intentionally kept simple to let the reader focus on the global method. More complex models for non-conformity rate or cost evaluation can easily be applied using the same framework. The method is currently under test on industrialscaled use cases.

B. Assembly technique library

The assembly technique library used in the case study presented in section 6 is detailed in Table B

Fig. 2 :

 2 Fig. 2: Example of an Oriented Contact Graph as implemented in GAIA R . The Datum Flow Chain of the key characteristic 1 is represented.

Fig. 3 :

 3 Fig. 3: Assembly Sequence Graph dual to the Oriented Contact Graph in Fig. 2.

Fig. 4 :

 4 Fig.4: Non-dominated elements (red circles) among a population (blue squares) for a two-objective problem.

Fig. 5 :

 5 Fig. 5: Flowchart of the optimization process, with reference to the associated equations.

Fig. 6 :

 6 Fig. 6: Schematic view of the use case.

Fig. 7 :

 7 Fig.7: Structuro-functionnal model of the box (after nine assembly steps) as represented in GAIA R[START_REF] Falgarone | An innovative design method and tool for structural and functional analysis[END_REF].

Fig. 8 :

 8 Fig. 8: Structuro-functionnal model of the box (after sixteen assembly steps) and Datum Flow Chains of KC3.

Fig. 9 :

 9 Fig. 9: Details of the component surfaces.

Fig. 10 :

 10 Fig. 10: Cost vs. Tolerance function with the set of parameters given in Table3.

Fig. 11 :

 11 Fig. 11: Optimization result: Pareto Front obtained with P = 100.

  Fig. A.1 in Appendix A where the number of products P is equal to 10.

Fig. 12 :

 12 Fig. 12: Optimization result: Pareto Front obtained with P = 10000.

Fig. A. 1 :

 1 Fig. A. 1: Optimization result: Pareto Front obtained with P = 10.

Table 1 :

 1 Example of assembly technique definition for a temporary joint between a component and a tool (fixed cost given in cost unit and duration given in time unit, N (0, 0.1) stands for a normal distribution with mean value of 0 and standard deviation value of 0.1).

	Technique name:	Tool
	Operations	Fixed cost Duration Resources	Quantities Distribution
	Docking	0	15	Tool; Operator	1;2	N (0, 0.1)
	Clamping	3	10	Tool; Operator	1;1
	Unclamping	3	10	Tool; Operator	1;1
	Splitting	0	15	Tool; Operator	1;2

Table 2 :

 2 Authorized assembly techniques for the joints of the product.

Table 3 :

 3 Probabilistic distributions and cost vs. tolerance model parameters of the components geometrical deviations (cu stands for cost unit).

	T lim,i a i	b i	m i	k i
	(mm) (cu) (cu/mm) (mm -1 )
	δl 1ab 0.01	0	200	1	1
	δl 2ab 0.01	0	200	1	1
	δl 3ab 0.01	0	50	1	1
	δl 4ab 0.01	0	50	1	1

Table 4 :

 4 Assembly techniques and tolerances for the square point in Fig.10.

	Joint	Assembly technique selected
	tj 1 and tj 2	1-Back-to-back positioning (temporary)
	tj3 to tj 7	4-Robot as tool
	J1 -m	6-Back-to-back positioning
	J1 -c and J2 7-Traditional bonding
	J 3 and J 5	9-Back-to-back positioning (with splitting)
	J 4 and J 6	10-Drilling with grid
	Geometrical characteristic	Tolerance allocated
	δl 1ab and δl 2ab		U(-0.40, 0.30)
	δl 3ab		U(-0.22, 0.21)
	δl 4ab

  .1, Table B.2, Table B.3. The characteristics of the resources required for the assembly operations of this library are given in Table B.4. Table B. 4: Resource list and characteristics.

	Resource Name	C r	C tu
	Operator	1000	0.833
	Non-contact measurement device 20000	0.067
	Glue gun	150	0.017
	Tool	5000	0.017
	Adjustable Tool	20000	0.033
	Robot	50000	0.7
	Specific Automated Station	1500000 0.8
	Grid	5000	0
	Numerical Grid	12000	0.01
	Drilling/Riveting End-effector	50000	0.05
	Specific Drilling/Riveting Station 1500000 0.8
	Automated Drilling Unit	2500	0.01

Preprint submitted to Journal of Manufacturing Systems December 12, 2014

Running on a single 2.53Ghz CPU.