Edoh Goka

Pierre Beaurepaire

Lazhar Homri

Jean-Yves Dantan

Probabilistic-based approach using Kernel Density Estimation for gap modeling in a statistical tolerance analysis

The statistical tolerance analysis has become a key element used in the design stage to reduce the manufacturing cost, the rejection rate and to have high quality products. One of the frequently used methods is the Monte Carlo simulation, employed to compute the non-conformity rate due to its efficiency in handling the tolerance analysis of over-constrained mechanical systems. However, this simulation technique requires excessive numerical efforts. The goal of this paper is to improve this method by proposing a probabilistic model of gaps in fixed and sliding contacts and involved in the tolerance analysis of an assembly. The probabilistic model is carried out on the clearance components of the sliding and fixed contacts for their assembly feasibility considering all the imperfections on the surfaces. The kernel density estimation method is used to deal with the probabilistic model. The proposed method is applied to an over-constrained mechanical system and compared to the classical method regarding their computation time.

1.

Introduction

In a competitive industrial area, companies work continuously to improve the quality of their products, reduce the manufacturing costs and decrease the number of defective products. Tolerance analysis is a key tool integrated into the design phase of the manufacturing products, which allows designers to achieve such goals [START_REF] Chen | A comprehensive study if three dimensional tolerance analysis methods[END_REF][START_REF] Schleich | A discrete geometry approach for tolerance analysis of mechanism[END_REF]. Tolerance analysis aims to verify if the individual tolerance specified to each part in the product allows the designer to reach the assembly feasibility and the functional requirements. Two methods for the tolerance analysis of mechanical assemblies are usually used: the worst-case and the statistical approaches. The worst-case tolerancing often leads to tight tolerances and high production costs. The statistical tolerancing is widely used because it leads to large tolerances and allows practitioner to get a desire yield and low production cost [START_REF] Hong | A comprehensive review of tolerancing research[END_REF]. The statistical tolerance analysis can be achieved based on several methods [START_REF] Hong | A comprehensive review of tolerancing research[END_REF][START_REF] Nigam | Review of statistical approaches of tolerance analysis[END_REF][START_REF] Skowronski | Using Monte-Carlo variance reduction in statistical tolerance synthesis[END_REF]: the root-sum-squares method, Croft's method, extended Taylor series approximation, Hasofer-Lind index method, approximation by numerical integration or quadrature technique, Taguchi's method and the Monte Carlo simulation. One of the current challenges in tolerancing is the consideration of the over-constrained mechanical systems with gaps. The assembly of such mechanical systems is possible because of the gaps [START_REF] Dumas | A new system formulation for the tolerance analysis of overconstrained mechanisms[END_REF], even though the parts exhibit geometric imperfections. Consequently, the study of over-constrained mechanical systems is more complex than the iso-constrained mechanical systems. Several methods to study over-constrained mechanical systems have been proposed. The combination of the set point variation, homogeneous transform matrix and Monte Carlo simulation has been proposed in [START_REF] Yang | Statistical Tolerancing based on Variation of Point-set[END_REF]. An alternative strategy consists of identifying the contact configurations of the mechanical assembly [START_REF] Beaucaire | Statistical tolerance analysis of overconstrained mechanisms with gaps using system reliability methods[END_REF][START_REF] Ballu | Geometrical reliability of overconstrained mechanisms with gaps[END_REF]. Thus, some system reliability methods are used to compute the non-conformity rate of the overconstrained mechanical assemblies. However, for these system reliability methods, even though the gap is considered, it is up to the designer to define the possible contact configurations of the mechanical assembly, which may be a challenging task. Moreover, these methods are not able to consider parts' form defects and are globally applied on 2D cases. Furthermore, to overcome these previous limitations, a new method based on the Lagrange dual form of the optimization problem and FORM method has been proposed to compute the probability of failure with reduced computational efforts and carrying out a preliminary selective search algorithm to identify the configurations which contribute mostly to this probability of failure [START_REF] Dumas | A new system formulation for the tolerance analysis of overconstrained mechanisms[END_REF]. Though this work considers 3D over-constrained mechanical systems with multi-configuration of gaps, it has not considered the form defects nor the different types of contacts in the tolerance analysis procedure. Therefore, this paper proposes a probabilistic model of gaps in fixed and sliding contacts, applied in a global statistical tolerance analysis method. The statistical tolerance analysis is carried out based on constrained optimization and Monte Carlo Simulation (MCS). The probabilistic model is proposed to determine the gap components in fixed and sliding contacts between non-ideal surfaces, instead of computing these components using the MCS combined with an optimization procedure. The Kernel Density Estimation (KDE) method is used to define the probabilistic model. A first optimization is carried out on the gap components of each sliding or fixed contact. Then the probabilistic model is identified based on these values gathered from this first optimization to generate new samples in good agreement with the distribution and the correlation of the original gap-values. For the floating contacts a set of constraints is developed to prevent the inferences of the mating surfaces. Then, a general optimization is carried out to determine the values of their gap components. Finally, a Monte Carlo simulation is carried out based on all the previous developments to determine the probabilities of assembly failure and functional failure. The paper is structured as follows: Section 2 presents the main algorithm currently used to deal with the statistical tolerance analysis and highlights the step at which the probabilistic model of gap is integrated. Section 3 presents the KDE method used to determine the probabilistic model of gaps, Section 4 describes an application of the proposed method to an over-constrained mechanical system, and finally, discussions and conclusion are given at the end.

2.

Tolerance analysis of over-constrained mechanical systems

Classical tolerance analysis approach

The statistical tolerance analysis method aims to calculate the probabilities of failure of the assembly and functioning requirements of a mechanical system. Generally, three main issues are considered in the tolerance analysis procedure: the geometrical deviations modeling, the geometrical behavior modeling of the mechanical systems and the mathematical formulation of the tolerance analysis techniques [START_REF] Qureshi | A statistical tolerance analysis approach for over-constrained mechanism based on optimisation and Monte-Carlo simulation[END_REF]. The main contribution of this paper concerns the geometrical behavior modeling of assemblies involving form defects, where a probabilistic model is proposed for gaps modeling in fixed and sliding contacts. This section presents a state of art of the statistical tolerance analysis approach. It is necessary to use a realistic model of the surfaces of a mechanical system (i.e. modeled imperfections in good agreement with the imperfections observed on the manufactured components) to guarantee accurate results of the tolerance analysis. Globally, three models of surfaces are applicable: the nominal surface [START_REF] Dumas | Impact of a behavior model linearization strategy on the tolerance analysis of over-constrained mechanisms[END_REF], the substitute surface [START_REF] Dumas | Impact of a behavior model linearization strategy on the tolerance analysis of over-constrained mechanisms[END_REF], and the non-ideal surface represented by skin model concept [START_REF] Schleich | Skin Model Shapes: A new paradigm shift for geometric variations modelling in mechanical engineering[END_REF]. Both, substitute and non-ideal surface, have an important role in the tolerance analysis method. Several models have been proposed to characterize the geometrical deviations of surfaces [START_REF] Qureshi | A statistical tolerance analysis approach for over-constrained mechanism based on optimisation and Monte-Carlo simulation[END_REF][START_REF] Wittwer | The direct linearization method applied to position error in Kinematic linkages[END_REF][START_REF] Arroyave-Tobón | Applying screw theory for summing sets of constraints in geometric tolerancing[END_REF]. In this paper, the SDT approach [START_REF] Bourdet | The concept of the small displacement torsor in metrology[END_REF] is used to model surface deviations and gaps between the mating surfaces. Thus, each surface is characterized by three components of translation (u, v, w) and three components of rotation (α, β, γ). The set of dimension, position and orientation deviations of the mechanical system is denoted by X and is randomly generated. González and Sánchez [START_REF] González | Statistical tolerance synthesis with correlated variables[END_REF] proposed a method for the optimal tolerance allocation using the statistical approach with dependent (e.g. correlated) variables; the structure of dependence can be associated with the manufacturing process. The non-ideal surface is built by adding form defects to the substitute surface. The method used to build the non-ideal surface has been explained in [START_REF] Homri | Tolerance analysis -form defects modeling and simulation by modal decomposition and optimization[END_REF][START_REF] Goka | Geometrical variation simulation for assembly with form defects[END_REF]. The vector F used throughout this paper, represents the set of form defects amplitudes of the surfaces defined in correlation with the discretization of the surfaces. The geometrical behavior modeling of an over-constrained mechanical system is obtained through the development of a set of constraints [START_REF] Qureshi | A statistical tolerance analysis approach for over-constrained mechanism based on optimisation and Monte-Carlo simulation[END_REF][START_REF] Dumas | Impact of a behavior model linearization strategy on the tolerance analysis of over-constrained mechanisms[END_REF][START_REF] Goka | Geometrical variation simulation for assembly with form defects[END_REF][START_REF] Dantan | Worst case and statistical tolerance analysis based on quantified constraint satisfaction problems and monte carlo simulation[END_REF]: compatibility equations (denoted by

  ,0 c C  XG ),
interface constraints (denoted by

  0 i C  X,F,G
), and the functional conditions (denoted by

  0 f C  X,F,G
).

These constraints have been detailed in [START_REF] Goka | Statistical tolerance analysis of over-constrained mechanical assemblies with form defects considering contact types[END_REF] The assembly and functioning problems must be formulated to perform the tolerance analysis of an overconstrained mechanical system. The formulation of the assembly and functioning problems have been detailed in [START_REF] Qureshi | A statistical tolerance analysis approach for over-constrained mechanism based on optimisation and Monte-Carlo simulation[END_REF][START_REF] Dantan | Tolerance synthesis: quantifier notion and virtual boundary[END_REF] and the mathematical formulations of the probabilities of assembly and functioning failure are given in [START_REF] Goka | Statistical tolerance analysis of over-constrained mechanical assemblies with form defects considering contact types[END_REF]. To perform the computation of these two probabilities, a Monte Carlo simulation combined with a constrained optimization is employed. This method has been introduced by Qureshi et al. [START_REF] Qureshi | A statistical tolerance analysis approach for over-constrained mechanism based on optimisation and Monte-Carlo simulation[END_REF] considering only size, position and orientation defects; and further extended to integrate parts' form defects [START_REF] Homri | Tolerance analysis -form defects modeling and simulation by modal decomposition and optimization[END_REF][START_REF] Goka | Geometrical variation simulation for assembly with form defects[END_REF] and to consider different contact types [START_REF] Goka | Statistical tolerance analysis of over-constrained mechanical assemblies with form defects considering contact types[END_REF]. The overall tolerance analysis technique of an over-constrained mechanical system is explained in Algorithm 1.

ALGORITHM 1:

Computation of the probability of assembly failure ( na P ) and the probability of functional failure ( nf P ) Require: Total number of Monte Carlo Samples (N mc ) Require: Total number of fixed and sliding contacts in the over-constrained mechanical system (N ct ) Ensure: total number of Monte Carlo populations and compute the probabilities na P , nf P

1: Constraints development:   , c C XG ,   ,, i C X G F ,   f C X,G,F
// min f(G') subject to   , c C XG ,   , , 0 i C  ' X G F ,   0 f C  X,G,F 11:
End for 12: Estimation of na P and nf P , by using the following equations [22]:

  ( ) ( ) na 1 1 P mc na N ii D i mc I N    X , F ,   ( ) ( ) nf 1 1 P mc nf N ii D i mc I N    X , F ,
  1 if no solution is found 0 if a solution is found D I     X, F
and for the functional requirements, it is represented by:

      1 if 0 0 if 0 f D f C I C         X,G, F X, F X,G, F
Algorithm 1 sums up the local simulation to characterize the fixed and sliding contacts and the general optimization carried out on the floating contacts with form defects considered [START_REF] Goka | Statistical tolerance analysis of over-constrained mechanical assemblies with form defects considering contact types[END_REF]. The steps highlighted in Algorithm 1 (framed steps) represent the focus of this paper, and a novel procedure is proposed. This local simulation on the fixed and sliding contacts is discussed in detail in the upcoming subsection.

Mathematical formulation of the assembly feasibility of the fixed and sliding contacts of parts with form defects

The goal of the contact study in this section is to determine the clearance components which allow an admissible assembly configuration in case of fixed or sliding contacts. Several methods to deal with different contacts fitting (fixed, sliding and floating) have been proposed. Some researchers have proposed numerical optimization approaches with several objective functions (convex hull volume, Euclidean distance, weighted unsigned distance, weighted signed distance) [START_REF] Schleich | Contact and mobility simulation for mechanical assemblies based on skin model shapes[END_REF], experimental and theoretical investigations using the SDT and clearance domain (gap hull) models [START_REF] Lê | Experimental and theoretical investigations of mechanical joints with form defects[END_REF], the Fast Fourier Transform and the Conjugate Gradient Method in case of elastic contact of rough surfaces [START_REF] Hu | Numerical Analysis for the Elastic Contact of Real Rough Surfaces[END_REF][START_REF] Wang | Numerical solution of dry contact problem based on fast fourier transform and conjugate gradient method[END_REF], development of an optimization algorithm and geometric interference relationships [START_REF] Hummel | Conguration design and optimization of universal joints with manufacturing tolerances[END_REF] and some other methods in [START_REF] Putignano | A new efficient numerical method for contact mechanics of rough surfaces[END_REF][START_REF] Adragna | Statiscal assemblies with form errors -a 2D example[END_REF][START_REF] Grandjean | On the role of form defects in assemblies subject to local deformations and mechanical loads[END_REF]. Moreover, some others researchers proposed to study the floating contacts by means of a regularized closure function and a constrained optimization problem [START_REF] Rameau | Clearance vs. tolerance for rigid overconstrained assemblies[END_REF], Virtual Work Principle and a contact model of joints to determine the displacements at each node and relative displacements/rotations of the kinematic pairs [START_REF] Cammarata | A novel method to determine position and orientation errors in clearance-affected overconstrained mechanisms[END_REF], constrained registration techniques and difference surface in a sequential positioning problem [START_REF] Schleich | Approaches for the assembly simulation of skin model shapes[END_REF], Tooth Contact Analysis algorithm [START_REF] Schleich | A discrete geometry approach for tolerance analysis of mechanism[END_REF][START_REF] Lin | Kinematic error analysis and tolerance allocation of cycloidal gear reducers[END_REF], the Gap Hull estimation [START_REF] Schleich | Gap hull estimation for rigid mechanical joints considering form deviations and multiple pairs of mating surfaces[END_REF]. Among all the previous methods, the Signed Distance Approach [START_REF] Homri | Tolerance analysis -form defects modeling and simulation by modal decomposition and optimization[END_REF][START_REF] Goka | Geometrical variation simulation for assembly with form defects[END_REF][START_REF] Goka | Statistical tolerance analysis of over-constrained mechanical assemblies with form defects considering contact types[END_REF][START_REF] Schleich | Contact and mobility simulation for mechanical assemblies based on skin model shapes[END_REF] combined with an optimization technique is selected to characterize the fixed and sliding contacts assembly in the mechanical system. Thus, this subsection presents this local simulation on the gap components of a fixed or sliding contact added to the classical tolerance analysis approach. Based on the SDT approach [START_REF] Bourdet | The concept of the small displacement torsor in metrology[END_REF][START_REF] Homri | Tolerance analysis -form defects modeling and simulation by modal decomposition and optimization[END_REF][START_REF] Goka | Geometrical variation simulation for assembly with form defects[END_REF], the clearance torsor g of a contact can be written as follows:

M u v w          g (1)
where u, v and w represent the translations in the direction of the x-, y-and z-axes, respectively; α, β and γ represent the rotations around the x-, y-and z-axes respectively. Problems involving form defects in the sliding and fixed contacts are more complex than problems without form defects and specific procedures are required to perform the tolerance analysis.

Thus, two cases may be considered when a tolerance analysis problem involves fixed or sliding contacts: (1) form defects are not considered and only size, position and orientation defects are involved; (2) form defects are considered. In the first case, the representation of the clearance torsor depending on the different nature of contact is shown in Fig. 1a. The focus of this paper concerns the second case.

Fig. 1. Gap torsor for: (a) perfect planar contact, (b) contact of non-ideal planar mating surfaces

The second case concerns the contact of planar mating surfaces with form defects (See Fig. 1b). The clearance components depend on the form defects amplitudes and the contact configuration becomes non-trivial. The values of the clearance components are set to a specific value associated with the contact configuration.

For instance α=α*, β=β* and w=w* in the case of a sliding planar contact, where α*, β* and w* are non-zero even though the two surfaces are in contact. The optimization strategy of g* consists of minimizing the clearance components such that the signed distances between the two non-ideal surfaces are respected. This optimization technique is well detailed in [START_REF] Homri | Tolerance analysis -form defects modeling and simulation by modal decomposition and optimization[END_REF][START_REF] Goka | Geometrical variation simulation for assembly with form defects[END_REF][START_REF] Goka | Statistical tolerance analysis of over-constrained mechanical assemblies with form defects considering contact types[END_REF].

For the fixed contact, the six components of the clearance components are determined during the optimization procedure to fully characterize the contact. Regarding the sliding contacts, only three components are identified during the optimization (the translation component w, and the two rotation angles α and β) for a planar contact.

The three remaining components of the clearance torsor are not determined at this stage of the procedure. They may be interpreted as "free variables", i.e. these parameters are computed in the general optimization of the tolerance analysis which is associated with the floating contacts (as discussed in section 2.1). The components of g* depend on the nature of the contact and they are shown in Fig. 1b. More information about this clearance g* and the mathematical formalization of its optimization can be found in [START_REF] Goka | Statistical tolerance analysis of over-constrained mechanical assemblies with form defects considering contact types[END_REF].

The tolerance analysis algorithm combines a MCS and separate optimization methods for each sample. Both procedures are numerically demanding and combining them may lead to excessive computational efforts. Hence, the probabilistic models are developed to estimate the distributions of the clearance components for fixed and sliding contacts.

Probabilistic model for gap modeling in fixed and sliding contacts

In this section, a probabilistic model is proposed for gap components in fixed and sliding contacts to determine their relative positioning without interferences. In the statistical tolerance analysis procedure described above, the optimizations for fixed and sliding contacts are carried out for each sample generated during the Monte Carlo simulation and that is very time-consuming. Therefore, the joint probability density functions associated with the clearance components are identified and subsequently used to generate samples instead of using the optimizations carried out on these contacts. The Kernel Density Estimation (KDE) approach is selected to calculate these joint probability density functions. The main goal of using the KDE approach is to generate new samples from the original gap components optimized in such manner that these new samples have the same distribution and features (mean, standard deviation, etc.) as the original values. In Section 3.1, the mathematical definition of the KDE approach is given. Section 3.2 presents the method of sampling using the original data from the classical optimizations and the Section 3.3 describes how these samples from the probabilistic model are used in the tolerance analysis.

3.1

Kernel Density Estimation approach.

The Kernel Density Estimation (KDE) is a statistical non-parametric method used to estimate the Probability Density Function (PDF) of a random variable [START_REF] Tai | Resampling of Data for Offshore Grid Design Based on Kernel Density Estimation and Genetic Algorithm[END_REF]. This method can be used for univariate and multivariate continuous variables. The KDE may be used to estimate a probability density function [START_REF] Silverman | Density estimation for statistics and data analysis, Chapman&Hall/crc[END_REF] like some other methods: histograms, naïve estimator, nearest neighbor method, variable kernel estimator, orthogonal series estimators and maximum penalized likelihood estimators. The KDE approach is applied on the realization of random variables of the clearance components to approximate their probability density. The kernel density estimator of a univariate continuous variable y is defined as:

  1 1 n i i yy f y K nh h        (2)
where K(.) represents the kernel function centered on the observation i y , h is the bandwidth or the window width. Several kernels are available in the literature and can be used in this kernel estimation [START_REF] Silverman | Density estimation for statistics and data analysis, Chapman&Hall/crc[END_REF]: Epanechnikov kernel, Biweight kernel, Triangular kernel, Gaussian kernel, rectangular kernel. In our work, we chose to use the Gaussian kernel defined as follows:

  2 11 exp 2 2 K u u      (3)
where

  =/ i u y y h  .
Here, the clearance components are multivariate continuous variables which can be correlated. Therefore, the Multivariate Kernel Density Estimation (MKDE) is used. Thus, the MKDE is defined in general by the same function [START_REF] Rajagopalan | Multivariate nonparametric resampling scheme for generation of daily weather variables[END_REF] as in Eq. ( 3) (generalized for the multivariate case) but

    1 2 = T yy u h   ii Y S Y ,   . K is the multivariate Gaussian kernel function,   12 , ,..., d y Y Y Y 
represents the d dimensional random vector whose density is being estimated,   12 , ,...,

T i i di Y Y Y  i Y
i=1 to n (n representing the total number of sample available to construct the KDE estimation) represents the sample value of y and S represents the sample covariance matrix. The multivariate Gaussian kernel function [START_REF] Rajagopalan | Multivariate nonparametric resampling scheme for generation of daily weather variables[END_REF] is defined as:

        / 2 1/ 2 1 exp / 2 2 det d d K u u h   S (4)
An important parameter of the KDE approach is the bandwidth h which has a huge impact on the KDE model. The bandwidth h is used to control the amount of smoothing of the data in the density estimate [START_REF] Rajagopalan | Multivariate nonparametric resampling scheme for generation of daily weather variables[END_REF]. There are several techniques used to determine the best h: the cross-validation method, the plug-in method, the Scott's Rule and the Silverman Rule (Rule of thumb) [START_REF] Scott | Multivariate density estimation: theory, practice, and visualization[END_REF][START_REF] Heidenreich | Bandwidth selection for kernel density estimation: a review of fully automatic selectors[END_REF]. The cross-validation method and the plug-in method are based on an optimization procedure to find the best h value, and require considerable numerical efforts [START_REF] Rajagopalan | Multivariate nonparametric resampling scheme for generation of daily weather variables[END_REF]. Thus, in this paper, Silverman appropriate h proposed in [START_REF] Silverman | Density estimation for statistics and data analysis, Chapman&Hall/crc[END_REF] is used for the Gaussian kernel as the optimal bandwidth. The bandwidth h proposed by Silverman is defined as:

          1/ 4 1/ 4 4 / 2 1 . d d h d n    (5)
where n represents the number of the data points, d the dimension of the set of random variables. After describing the KDE approach, the resampling method is detailed in the next section.

Generation of new KDE based variables using the original data.

The KDE approach is introduced in this article to determine firstly, the probability density of the variables and secondly, to generate new samples or realizations based on the KDE in order to reduce the numerical efforts. The process to generate the new samples which are going to be used to replace the repetitive optimization is detailed in this section. In the classical tolerance analysis method presented above, a sample of variables X is generated in each iteration of the Monte Carlo simulation. For each sample X, an optimization is carried out on each fixed and sliding contact to determine its corresponding g* (see Section 2). To introduce the KDE simulation method, a preliminary Monte Carlo simulation is launched on the optimization of the components of g*. Realizing the preliminary Monte Carlo simulation on the optimization of g* from a fixed or sliding contact, allows one to get several realizations of each component of g*. These realizations are considered as the original data on which the KDE simulation method is applied. Thereafter, the objective is to approximate the joint distribution of these realizations of the components from g* by using the KDE simulation method. Several experiments on the optimization of g* from a fixed or a sliding contacts have shown that, there is a great variability on the marginal distributions of the components from g* and that the structure of dependence between them is complex. The standard methods of approximation of marginal distributions which on one hand include the use of simple predefined distributions (Gaussian, uniform, beta, gamma, etc.) and on the other hand, the use of usual copulas for the dependence structure (Gaussian copula, Archimedean copula …) are not flexible enough to describe with accuracy the distribution of each component from g* and their structure of dependence. Hence, it is necessary to find an approach that allows designers to generate new samples with various variables' distribution and which respects also to the dependence structure between the new simulated variables. The KDE simulation method has been chosen because of its simplicity, flexibility and accuracy to approximate the marginal distribution of each component of g* and to model the structure of dependence between these components of g*. An algorithm is developed to generate new samples of the variables based on the KDE approach. The method to generate the KDE simulation data is explained as follows. Let us consider N mc as the total number of Monte Carlo simulation, k the number of components of the clearance parameters and KDE N the total number of the preliminary MCS. The local simulation is performed N KDE times and the results are kept in a matrix. After the last iteration, the completed matrix is considered as the original values of the clearance components. This original data is denoted by D and represents a N KDE -by-k matrix. Then, the probabilistic model based on the Kernel density Estimation is launched by using the matrix D. The sampling method is described by Rajagopalan et al [START_REF] Rajagopalan | Multivariate nonparametric resampling scheme for generation of daily weather variables[END_REF] and summarized in Algorithm 2. 8: Generate a multivariate gaussian variables N centered on 0 and has h 2 xS as the standard deviation, N is a matrix of size N mc x k 9: Compute P = M + N. Each column from P represents the probabilistic model of each components from the clearance torsor Algorithm 3 is used to get the probabilistic model for each component from the clearance of each fixed and sliding contacts. Therefore, by using this algorithm 3, the distributions of the new samples of the clearance components are in good agreement with the original data as well.

Integration of the KDE based variables in the statistical tolerance analysis approach

The classical approach for tolerance analysis of over-constrained mechanisms involves a probabilistic model of the geometry of the parts, i.e. the dimension, orientation, shape or form defects are characterized by a set of random variables. Monte Carlo simulation is used to generate samples of the geometry and for each sample, the components of the gap torsor are analyzed. This task is non-trivial and requires to perform (1) an optimization of the gap components for each fixed and sliding contact to compute the value of the limited degrees of freedom;

(2) optimizations on the free degrees of freedom of the contacts, which have not been previously identified, in order to determine whether the mechanism can be assembled and whether it is functional. This method allows us to compute the probability of non-functionality and the probability of non-assembly.

Fig. 2. Frameworks: (a) classical method, (b) using the proposed method

The proposed method in this paper is focused on the first step described above. It is suggested to investigate each sliding and fixed contact individually. Samples of the imperfect geometry are generated using Monte Carlo simulation and the corresponding components of the gap torsor are computed. This data is subsequently used to identify the joint probability density function of the limited displacements components of the gap torsor using the KDE. Once this probability density function is characterized, it is possible to use it in order to directly generate samples of the components of the gap torsor without using any numerically demanding optimization procedure. The interface constraints of the geometrical behavior of a mechanism, involve solely the floating contacts. Therefore, the samples of the proposed KDE based method are not involved in these interface constraints. The compatibility equations and the functional conditions are used in the step (2) detailed in the geometrical behavior modeling. They represent the constraints of a general optimization procedure. The variables in P are involved as inputs of the compatibility equations and the functional requirements. The compatibility equations and the functional requirements are concerned about the fixed, the sliding and the floating contacts. Thus, at each iteration of the MCS, a sample of gap components from P are introduced in the compatibility equations and in the functional requirements inequalities. The remaining degrees of freedom of the fixed and sliding contacts which are not limited, are considered as free variables in the step (2) described above, i.e. they are adjustable parameters of the optimization performed at this step.

The whole process for the tolerance analysis with the probabilistic models is depicted in Fig. 2b.

Application

In this section, the proposed probabilistic model approach is implemented to model the distribution of the gap components of a sliding contact in an over-constrained mechanical system, chosen as case study. In section 4.1, the mechanical system is presented and then in Section 4.2, the results of the probabilistic modeling are given as well as the results of the tolerance analysis. A comparison is addressed between the classical tolerance analysis method and the proposed method.

Description of the mechanical system

The case study is a part of an electric engine, see Fig. 3. The mechanical system is composed of a shaft, a housing and a body.

Fig. 3. Studied mechanical system

The mechanical system is composed of 4 floating contacts (1b/2b, 2c/3c, 1d/3d, and 1e/3e) and one sliding contact (1a/2a). The sliding contact is the most important for the approach proposed in this paper because the probabilistic models are established on its clearance parameters. The functional requirement characterizes the floating contact 1d/3d, where a minimum clearance threshold should be respected between the shaft and the housing. The tolerance analysis of this mechanical system is carried out based on the tolerances specifications taken from [START_REF] Goka | Statistical tolerance analysis of over-constrained mechanical assemblies with form defects considering contact types[END_REF][START_REF] Anselmetti | Generation of functional tolerancing based on positioning features[END_REF]. The tolerance specifications of the different parts of the mechanical system used in these studies and more information about the other tolerances used in the overall tolerance analysis method can be found in our previous work [START_REF] Goka | Statistical tolerance analysis of over-constrained mechanical assemblies with form defects considering contact types[END_REF]. Firstly, a preliminary MCS is launched on the optimization of the gap components in contact 1a/2a. Then, the KDE method is used to get new samples of these realizations, which are after introduced in the overall tolerance analysis method to replace the highlighted steps from algorithm 1. The contact 1a/2a is composed of two circular planar surfaces. The form defects are generated according to the z axis. The contact 1a/2a is sliding because body is bolted to the shaft and this assembly technique allows slight displacement in the direction of the contact plane. Thus, three clearance parameters are identified as shown in the following clearance g*:

M w           * g (6)
where α, β are components of the rotations around the x-and y-axes, and w is the translation in the z-axes. These components of g* are not trivial and are dependent on form defects. An optimization procedure is used to identify them. The probabilistic model method proposed in Section 3 is implemented on these three components. The rest of g* components (γ, u, v) are not determined during this optimization or probabilistic model method.

They are considered as free variables and are included in the optimization of G'. The results for the components w, α, β are presented and interpreted in the next section.

Results

In this subsection, two categories of simulations have beeen launched to apprehend the different aspects of the use of the KDE method in the tolérance analysis such as the reduction of simulation time, the suitable meshing to modelize corectly the geometrical behavior of the mechanical system. The first categorie concerns the application of two differents distributions of the mating surfaces of 1a/2a in order to dertemine the appropriate meshing that modelize properly the geometrical behavior of the mechanical system with steady calculated probabilities. The computing time is not the important aspect of the simulations of this first categorie. The Gaussian distribution and the uniform distribution are considered here to generate the amplitudes of the form defects of the surfaces involved in the sliding contact (contact 1a/2a). These distributions are employed with different standard deviations and supports to test the robutness of the proposed method. The Gaussian distribution is composed of three different standard deviations with the same mean. The uniform distribution is also composed of three different supports. These two distributions are chosen in order to determine whether the distributions of the amplitudes of the form defects have an impact on the distribution of the components of g* and if the KDE is capable of approximating the marginal distributions of the clearance components and keeping their structure of dependence. Three levels of meshing refinement of the mating surfaces are considered: the coarse, the standard and the refined mesh with respectively 100, 2500, and 10000 discretization points. The consideration of multiple levels of refinement allows us to determine whether the mesh is sufficiently fine to capture the behavior of the components of g*, and to detemine whether the probabilities (of assembly failure and functional failure) have converged towards their reference value. Table 1 shows the different strategies considered to model the form defects. Uniform 100 2500 10000

[2x10 -3 , 3x10 -3 ] [1.5x10 -3 , 3.5x10 -3 ] [-7.5x10 -3 , 1x10 -2 ]
Based on Table 1, all the combinations possible of parameters shown are studied. Two distributions (Gaussian and uniform) are available to generate the amplitudes of form defects. For each distribution, the three levels of refinement of the mesh are used and for each level of refinement, the three standard deviations (for the Gaussian distribution) or the three supports (for the uniform distribution) are used to generate the form defects amplitudes.

In total, considering the combination between the types of meshing and the standard deviation or support, 18 types of form defects amplitudes are generated (nine for the Gaussian distribution and nine for the uniform distribution). For each distribution of the considered form defects amplitudes the probabilistic model of g* of the sliding contact is established and the tolerance analysis of the whole over-constrained mechanical system is thereafter carried out. As an example of simulation, the form defects may be generated using 100 discretization points and amplitudes following a Gaussian distribution with a mean equal to 2.5x10 -3 mm and a standard deviation equal to 10 -6 . After every simulation, some comparisons are carried out between the KDE simulation variables and the original variables from the local simulation (i.e. obtained using an optimization procedure). These comparisons aim to verify whether the simulated data from the KDE correspond to the original data and if they have the same features. They involve the components w, α, β of the clearance parameters, their histograms, and the Probability Density Function (PDF). As all the comparison figures of the 18 simulations cannot be shown together, just figures from the Gaussian distribution N(2.5x10 -3 , 10 -4 ) with 100 discretization points are presented. Nevertheless, the results concerning the probabilities (P nf , P na ) and the time of simulations assessed during the tolearance analysis carried out after, are presented for all the 18 simulations (see the Appendix). For each simulation, the matrix P defined in the Algorithm 3 can be written as: 1. The total number samples used in the Monte Carlo Simulation is N mc = 100000.  3 components of g* (w, α, β) for the first local simulation (optimization) and then after the probabilistic model of the data are determined.  The number of the optimization constraints (signed distances) depends on the mesh discretization.

  1 1 1 2 2 2
The reference simulation involves Gaussian distribution for the amplitudes of form defects, N(2.5x10 -3 , 10 -4 ) with 100 as meshing type. Based on the simulation from this example, the data from the local simulation (original data) and the data from the KDE simulation are shown respectively in Fig. 4a and Fig. 4b. In these figures, the histogram of each component is shown as well as the repartition between the components. It should be noticed that the original data have the same histograms as the KDE simulation data. Furthermore, the structures of dependence are in good agreement.

Fig. 4. Data scatterplot matrices: (a) Original data (b) KDE simulation data

Fig. 4 and some other figures presented in the appendix (Fig. 5 and Fig. 6) compare the data obtained using the classical approach (MCS combined with an optimization of the gap torsor components) and the proposed KDE method in the particular case where the amplitude of the form defects exhibits a Gaussian distribution with a mean equal to 2.5x10 -3 , a standard deviation equal to 10 -4 and the mesh discretization involves 100 points. The same comparison about the distribution of data points was performed for the other configurations described in Table 1, and we observed that globally, the KDE samples are in good agreement with the Monte Carlo data; however, these comparisons are not shown here and neither in the appendix to keep the paper short and concise.

The KDE simulation data are collected for each simulation and introduced mainly in the compatibility equations already developed. For each kind of distribution of the form defects amplitudes, two tolerance analysis simulations are run. The first tolerance analysis simulation is based on the technique explained in Section 2.1 i.e. based on the direct optimization of the clearance parameters of the sliding contact at each iteration (local simulation) of the Monte Carlo simulation [START_REF] Goka | Statistical tolerance analysis of over-constrained mechanical assemblies with form defects considering contact types[END_REF]. The second tolerance analysis simulation uses the KDE simulation variables as described in Section 3.2. A comparison is then established between these two tolerance analysis simulations as shown in the Tab. 2 for one type of simulation (F ⁓ N(2.5x10 -3 , 10 -4 )). For the other cases of simulation their results are presented in the tables presented in the appendix (Table 4-Table 8). It is important to note that the total simulations time presented in these tables in case of tolerance analysis using the KDE simulation samples, includes the time of the preliminary MCS applied to obtain the original data, the time to get the joint PDF of the clearance components and the new samples and finally the time to carry out the whole tolerance analysis simulation. Thus, 1

T represents the addition of the time to carry out the preliminary MCS to get the original data and the time to generate new samples (matrix P) of gap components based on the KDE method; 2 T represents the time to carry out the tolerance analysis of the mechanical system using the KDE simulation samples. This particular example does not imply multiple fixed, sliding and floating contacts during simulations. The gain in computation time can be much greater on mechanical systems with more fixed and sliding contacts. P nf = 1.64x10 -2 P na = 1.00x10 -5 P nf = 1.64x10 -2 P na = 1.00x10 -5

T 1 13min 47min 3h36min

T 2 9h13min 9h39min 9h16min
Total time (T 1 + T 2 ) 9h26min 10h26min 12h52min
The first results shown in Tab. 2 and the others from the appendix imply the reduction of the simulation time for the tolerance analysis based on the KDE simulation data. The speedup obtained using the KDE increases as the mesh is refined. The simulations reduction time is on average 14.50% for the mesh with 100 discretization points, 33.55% for the mesh with 2500 discretization points and 73.24% for the mesh with 10000 discretization points. Consequently, the tolerance analysis method based on the probabilistic models allows a gain of simulation time.

The second finding concerns the probabilities assessed using the different types of meshing. It is observed that the meshes with 2500 and 10000 discretization points have similar values of the probability of assembly failure (P na ), whereas its value is significantly different if the mesh with 100 discretization points is used. A similar conclusion can be drawn regarding the probability of functional failure (P nf ). Gradually as the mesh refinement is increased, the probabilities converge towards their reference value. The choice of the number of the discretization points for this tolerance analysis problem is solved in this manner. Based on the results and considering the constraints in terms of simulation time and accuracy of the results, the tolerance analysis using the meshing refinement of 2500 points can be taken as the best compromise solely for these configurations of meshing refinement. The third finding concerns the impact of the distributions on the probabilities. For both distributions (Gaussian and uniform) chosen to generate the amplitudes of the form defects on the mating surfaces, and for all the types of meshing refinement, the probabilities of assembly failure are increasing gradually as the amplitudes are increasing. This result is in accordance with engineering experiences. Indeed, as the amplitudes of the form defects are increasing, the mating surfaces are more subject to interferences and that can prevent the assembly to occur. It is observed then that the probabilities of assembly failure are very sensitive to the choice of the distributions, whereas the distributions have a reduced impact on the probabilities of functional failure. The final finding concerns the simulation times, which are increasing gradually as the meshing density increases. Indeed, the total number of constraints involved in the tolerance analysis increases as the mesh is refined, which increases the numerical efforts.

NN 

) of samples which are going to be used in the tolerance analysis.

mc N samples are sufficient to allow the tolerance analysis process to evaluate steady and good probabilities. The numerical efforts are reduced because the algorithm of this proposed method and also because a small number of samples is used in the preliminary MCS to generate the original data and identify the joint PDF associated with it. Thus, to show the advantage of the reduction of samples number in the preliminary MCS and the gain in time associated with it, simulations have been launched on the distribution F ⁓ U(1.5x10 -3 , 3.5x10 -3 ) with a total number of MCS 3, where it is observed that the probabilities can be fairly estimated by using just 10000 samples to generate the KDE data. These results show that it is not necessary to have the number of the preliminary MCS equal to the second MCS where the probabilities are calculated. By doing this, there is a gain in computing time about 6.30% based on the example shown in Table 3, between the simulation in which To sum up all the results, the first simulations based on the variations of the meshing and the distribution of form defects amplitudes, have helped to determine the best meshing choice resulting on steady probabilities of failure.

The second batch of simulations on the reduction of the number KDE N have allowed to determine that it is not necessary to have the same number of samples in the preliminary MCS as in the second MCS. And that is very important to reduce the simulation time. Furthermore, based on the results of Tab 3 and Tab 7 (in the Appendix), the simulations based on the computer formatting and development of the algorithm from the proposed method are faster than the simulations of the algorithm from the classical method. This explains the important reduction of time even if the N KDE is equal to N mc .

5.

Discussions

The application of the proposed method has shown in the previous section some of its advantages. However, some others can be added.

 Some mechanical systems are composed of several identic contacts. In this case, the KDE method is very suitable because it allows designers to launch just one time, a simulation to determine the original data of a contact. Thereafter, this original data can be used in the evaluation of the joint PDF of its contact and also for the rest of the identic contacts. This helps the designer to save computational time.

In fact, the identic contacts have the same elementary defects and the same characteristics. Therefore, they will have the same components of clearance for the best configuration allowing the assembly of the mechanical system. This gain in time is not possible with the classical method because a local simulation should be carried out on each of the identic contacts and at each iteration of the MCS. These repetitive simulations would increase the computational time. Furthermore, additional reduction of computing time can be obtained by using a reduced number of samples ( KDE N ) during the preliminary MCS as shown in Table 3.  Moreover, the proposed method allows to study each of the fixed or sliding contacts separately.

Therefore, a broad spectrum of these fixed and sliding contacts can be studied in advance to determine their original data. This will be consuming in time but it will help the designers to earn computational time once the design starts. During the design, studies must be carried out very fast to be competitive on the market. Thus, the joint PDF of the original data and the new samples will be obtained very fast to proceed to the tolerance analysis of the mechanical system. This computational time gain is possible because the original data of the contacts already exist and also because the KDE method is faster.  The use of the KDE method enables one to use another method than the MCS to study the tolerance analysis of over-constrained mechanical systems. This other method can be the First Order Reliability Method (FORM) for systems applied in [START_REF] Dumas | A new system formulation for the tolerance analysis of overconstrained mechanisms[END_REF] to compute the probabilities of failure. The use of this FORM method in not compatible with the classical statistical tolerance analysis method presented in Section 2.

Moreover, some probabilities calculated in the tables are slightly different. Nonetheless, these differences of probabilities between the local simulation and the proposed method are not sufficiently low to consider that the proposed method converges. Indeed, the probabilities calculated are of the same order of magnitude. It is reminded that the results are obtained using MCS which converges slowly. This method relies on the generation of pseudo-random numbers, and therefore two independent runs of the simulation lead to slightly different results.

Conclusion and future work

In this paper, a probabilistic-based model is proposed to deal with the gap components in fixed and sliding contacts. The model is then introduced in the global statistical tolerance analysis method of overconstrained mechanical system with form defects. This statistical tolerance analysis method is based on the Monte Carlo Simulation coupled with multi-scale optimization methods. In this new approach, the contacts (fixed, sliding and floating) are considered differently. The probabilistic models are developed to characterize the clearance parameters from the sliding and fixed contacts. The kernel Density Estimation (KDE) method is used to establish these probabilistic models of the clearance parameters. The main drawback of the Monte Carlo simulation is the considerable numerical efforts associated with it. These numerical efforts are reduced by proposing this new algorithm and by using reduced number of samples in a preliminary MCS which data are employed to get new samples through the KDE simulation method. Thereafter, these new samples are introduced in the general tolerance analysis method. Furthermore, simulations have been proposed to get best choice of the meshing to get steady probabilities of failure. The amplitudes of form defects have an impact on the probability of assembly failure just as the meshing points have an impact on the simulation time.

The proposed method cannot handle yet the deformations which occur on the different parts due the forces applied on them. These deformations can play a significant role in the assembly of the mechanical systems as shown in [START_REF] Liu | Integrating form errors and local surface deformations into tolerance analysis based on skin model shapes and a boundary element method[END_REF]. A possible outlook can be to see in which manner the studies carried out in [START_REF] Liu | Integrating form errors and local surface deformations into tolerance analysis based on skin model shapes and a boundary element method[END_REF] can be integrated in the tolerance analysis method of the present work or in the method addressed in [START_REF] Goka | Statistical tolerance analysis of over-constrained mechanical assemblies with form defects considering contact types[END_REF]. To show the similarity between the original data and the KDE simulation data in more detail, the distributions points between the clearance components are presented in the Fig. 5 for the case of simulation F ⁓ N(2.5x10 -3 , 10 -4 ). This Fig. 5 shows that the KDE based data globally correspond to the original data.

Nomenclature

Fig. 5. Comparison of the distribution between the clearance parameters (w, α, β)

Furthermore, the joint Probability Density Function (PDF) between each pair of KDE variables show the match between the original data points and the PDF contour plot. As shown in Fig. 6, for every combination of the clearance parameters, the density of the original data points is corresponding to the color of the PDF contour plot. P nf = 1.63x10 -2 P na = 1.00x10 -5 P nf = 1.62x10 -2 P na = 1.00x10 -5

T 1 15min 46min 3h28min

T 2 9h37min 9h39min 9h53min

Total time (T 1 + T 2 ) 9h52min 10h25min 13h21min

  F represents the indicator function; for assembly conditions, it is represented by:

ALGORITHM 2 : 1 :gStep 2 : 3 : 3 : 6

 212336 Sampling of new KDE-based variablesStep Let us consider D= * * * 12 , ,..., is of size 1 x k, i=1 to N KDE , obtained from the N KDE time optimization of the clearance parameters of a fixed or a sliding contact. Pick uniformly * matrix of size N mc x k. Simulate new variables from the kernel i.e. simulate multivariate gaussian variables centered on * i g of the new matrix and have as standard deviation 2 * h S . This simulation N ( * i g , 2 x h S ) is a matrix with dimension N mc x k. h represents the bandwidth whose calculation is explained in the previous section and S represents the sample covariance matrix.Step These new variables obtained, represent the probabilistic model of the clearance components, which are then introduced in the constraints (compatibility equations, interface constraints) developed.Complete algorithm of the overall method for the probabilistic model of the clearance components is summarized in the following ALGORITHM Probabilistic model of the clearance components of a fixed or a sliding contact Require: Set the number of Monte Carlo Samples (N mc ) Require: Set the total number of the preliminary MCS (N KDE ) Require: Ideal mating surfaces of a fixed or sliding contact Require: An empty matrix D of size N KDE x k Ensure: number of the optimization iterations at a contact 1: Generate the form defects amplitudes for each surface of the contact, meshing points on each surface 2: Apply the amplitudes to the perfect discretized mating surfaces to get the real mating surfaces i.e. surfaces with form defects 3: for each i N KDE do 4: min * () f g subject to Signed distances (S 1 , S 2 ) ≥ 0 (S 1 , S 2 : the non-ideal mating surfaces) 5:// Change the i th row of D by the components of * g matrix M of size N mc x k.

  summary for each simulation the following mathematical characteristics are assumed: N mc random variables following the Gaussian distribution X ⁓ N(µ, σ i ) and the uniform distribution U (a i , b i ) where µ=2.5x10-3 ,

Table 2 :

 2 Results of the simulation F ⁓ N(2.5x10 -3 , 10 -4 ), N KDE = 100000, Nmc = 100000

  with discretization type of 2500 and different numbers of KDE N . The results are shown in Table

  10000 and the one where it is equal to 100000.

Appendix A. 1

 1 Figures to show the similarity between the original data and the KDE simulation data

Fig. 6 .Table 4 :

 64 Fig. 6. Contour plot of the Probability Density Function of the KDE data compared to the original data points

Table 1 :

 1 The local simulation parameters

	Distribution	Meshing	Standard deviations/ Supports
		100	10 -6
	Gaussian (µ=2.5x10 -3 )	2500	10 -4
		10000	10 -2

  The second category of simulation concerns the use of reduced number of samples during the preliminary MCS in order to gain in computing time. During the first MCS to generate the original data, a small number KDE N of samples can be used and it is considered that this number KDE N of simulation, is sufficient to evaluate a steady joint PDF of the original data. This latter can be used to generate a large number

		N	(	NN 	or
		mc		mc	KDE
	mc	KDE

Table 3 :

 3 Results of the simulation F ⁓ U(1.5x10-3

			, 3.5x10 -3 ), N mc =100000, with different	N	KDE
	Simulations			Results
				2500
	N KDE	10000	20000	30000	50000	100000
	Probabilities	P nf = 1.63x10 -2 P na = 3.00x10 -5	P nf = 1.64x10 -2 P na = 3.00x10 -5	P nf = 1.63x10 -2 P na = 3.00x10 -5	P nf = 1.63x10 -2 P na = 3.00x10 -5	P nf = 1.63x10 -2 P na = 3.00x10 -5
	T 1	10min	13min	19min	24min	50min
	T 2	9h43min	9h44min	9h42min	9h45min	9h43min
	Total time (T 1 +	9h53min	9h57min	10h01min	10h09min	10h33min
	T 2 )

Table 5 :

 5 Results of the simulation F ⁓ N(2.5x10 -3 , 10 -2 ), N KDE = 100000, Nmc = 100000

	Simulation		Meshing
		100	2500	10000
	Classical tolerance analysis	P nf = 4.45x10 -3 P na = 8.41x10 -1	P nf = 4.39x10 -3 P na = 8.81x10 -1	P nf = 4.38x10 -3 P na = 8.81x10 -1
	Total time	12h36min	21h22min	66h23min
	Tolerance analysis based on KDE	P nf = 4.54x10 -3 P na = 8.49x10 -1	P nf = 4.39x10 -3 P na = 8.80x10 -1	P nf = 4.39x10 -3 P na = 8.81x10 -1
	simulation
	T 1	12min	44min	3h30min
	T 2	9h20min	9h19min	9h24min

Table 6 :

 6 Results of the simulation F ⁓ U(2x10 -3 , 3x10 -3 ), N KDE = 100000, Nmc = 100000

	Simulation		Meshing
		100	2500	10000
	Classical tolerance analysis	P nf = 1.72x10 -2 P na = 1.00x10 -5	P nf = 1.63x10 -2 P na = 1.00x10 -5	P nf = 1.63x10 -2 P na = 1.00x10 -5
	Total time	10h33min	13h59min	46h42min
	Tolerance analysis based on KDE	P nf = 1.71x10 -2 P na = 1.00x10 -5	P nf = 1.64x10 -2 P na = 1.00x10 -5	P nf = 1.64x10 -2 P na = 1.00x10 -5
	simulation
	T 1	13min	39min	3h44min
	T 2	9h42min	9h40min	9h42min
	Total time (T 1 + T 2 )	9h55min	10h19min	13h26min

Table 7 :

 7 Results of the simulation F ⁓ U(1.5x10 -3 , 3.5x10 -3 ), N KDE = 100000, Nmc = 100000

	Simulation		Meshing
		100	2500	10000
	Classical tolerance analysis	P nf = 1.73x10 -2 P na = 5.00x10 -5	P nf = 1.64x10 -2 P na = 1.00x10 -5	P nf = 1.64x10 -2 P na = 1.00x10 -5
	Total time	10h36min	14h07min	46h41min
	Tolerance analysis based on KDE	P nf = 1.71x10 -2 P na = 1.00x10 -5	P nf = 1.63x10 -2 P na = 3.00x10 -5	P nf = 1.64x10 -2 P na = 2.00x10 -5
	simulation
	T 1	12min	50min	3h32min
	T 2	9h33min	9h44min	9h36min
	Total time (T 1 + T 2 )	9h45min	10h34min	13h08min

Table 8 :

 8 Results of the simulation F ⁓ U(-7.5x10 -3 , 1x10 -2 ), N KDE = 100000, Nmc = 100000

	Simulation		Meshing
		100	2500	10000
	Classical tolerance analysis	P nf = 1.71x10 -2 P na = 2.00x10 -1	P nf = 1.69x10 -2 P na = 2.15x10 -1	P nf = 1.69x10 -2 P na = 2.15x10 -1
	Total time	10h39min	14h03min	46h51min
	Tolerance analysis based on KDE	P nf = 1.71x10 -2 P na = 2.09x10 -1	P nf = 1.67x10 -2 P na = 2.23x10 -1	P nf = 1.70x10 -2 P na = 2.24x10 -1
	simulation
	T 1	15min	49min	3h34min
	T 2	9h31min	9h47min	9h35min
	Total time (T 1 + T 2 )	9h46min	10h36min	13h09min

Acknowledgements

The authors would like to acknowledge the support of "AHTOLAND" project.


Maria Vakalopoulou

Guillaume Chassagnon

Norbert Bus

Rafael Marini Silva

Evangelia I Zacharaki

Marie-Pierre Revel

Nikos Paragios

R Marini

M.-P Revel

AtlasNet: Multi

Keywords: atlas Non-linear Deep Networks for Medical Image Segmentation Encoder-decoder, pixel-wise classification, deformable registration, interstitial lung disease

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

Introduction

Image segmentation is one of the most well studied problems in medical image analysis [START_REF] Paragios | hyper)-graphical models in biomedical image analysis[END_REF][START_REF] Iglesias | Multi-atlas segmentation of biomedical images: A survey[END_REF]. Segmentation seeks to group together voxels corresponding to the same organ, or to the same tissue type (healthy or pathological). Existing literature can be classified into two distinct categories, model-free and modelbased methods. Model-based methods assume the manifold of the solution space can be expressed in the form of a prior distribution, with sub-space approaches (e.g. active shapes), probabilistic or graphical models and atlas-based approaches being some representatives in this category. Model-free approaches on the other hand rely purely on the observation space combining image likelihoods with different classification techniques.

The emergence of deep learning as disruptive innovation method in the field of computer vision has impacted significantly the medical imaging community [START_REF] Zhou | Deep Learning for Medical Image Analysis[END_REF].

Numerous architectures have been proposed to address task-specific segmentation problems with the currently most successful technique being the Fully Convolutional Network (FCN) [START_REF] Shelhamer | Fully convolutional networks for semantic segmentation[END_REF]. Additionally, FCNs have been combined with upsampling layers, creating a variety of networks [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF][START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for image segmentation[END_REF], and have been extended to 3D [START_REF] Milletari | V-net: Fully convolutional neural networks for volumetric medical image segmentation[END_REF], boosting even more the accuracy of semantic segmentation.

The main challenges for deep learning in medical imaging arise from the limited availability of training samples -that is amplified when targeting 3D architectures -, the lack of discriminant visual properties and the three-dimensional nature of observations (high dimensional data). In this paper, we propose a novel multi-network architecture that copes with the above limitations. The central idea is to train multiple redundant networks fusing training samples mapped to various anatomical configurations. These configurations correspond to a representative set of observed anatomies and are used as reference spaces (frequently referred to as atlases). The mapping corresponds to a non-linear transformer. Elastic registration based on a robust, multi-metric, multi-modal graph-based framework is used within the non-linear transformer of the network. Training is performed on the sub-space and back-projected to the original space through a de-transformer that applies an inverse nonlinear mapping. The responses of the redundant networks are then combined to determine the optimal response to the problem.

The proposed framework relates also to the multitask learning (MTL) paradigm, where disparate sources of experimental data across multiple targets are combined in order to increase predictive power. The idea behind this paradigm is that by sharing representations between related tasks, we can improve generalization. Even though an inductive bias is plausible in such paradigms, the implicit data augmentation helps reducing the effect of the data-dependent noise. The idea of MTL for image segmentation has been incorporated before, such as in deep networks [START_REF] Ruder | An overview of multi-task learning in deep neural networks[END_REF] where soft or hard parameter sharing of hidden layers is performed, or in multi-atlas segmentation [START_REF] Iglesias | Multi-atlas segmentation of biomedical images: A survey[END_REF], where multiple pre-segmented atlases are utilized in order to better capture anatomical variation. As in most ensemble methods, the concept is that the combination of solutions by probabilistic inference procedures can offer superior segmentation accuracy.

AtlasNet differs from previous methods in respect to both scope and applicability. In (single or multi) atlas segmentation, the aim is to map a pre-segmented region of interest from a reference image to the test image, therefore applicability is limited to normal structures (e.g. organs of the body or healthy tissue) that exist in both images. Exploitability is even more reduced in the case of multi-atlas segmentation due to the rareness of multiple atlases. The proposed strategy on the contrary is suitable also for semantic labelling of voxels (as part of healthy or pathological tissue) without the requirement of spatial correspondence between those voxels in atlas and test image.

AtlasNet uses multiple forward non-linear transformers that map all training images to common subspaces to reduce biological variability and a backward de-transformer to relax the effect of possible artificial local deformations. In fact, due to the ill-posedness of inter-subject image registration, regularization constraints are applied to derive smooth solutions and maintain topological relationships among anatomical structures. Consequently, image registration does not always produce a perfectly conforming diffeomorphism due to the nonexistence of a single atlas that matches all anatomies. The use of multiple spaces comes to reduce the atlas selection bias, while the backward transformation aims to balance the effect of possible alterations in local image texture due to the nonlinearity in the transformation. Highly promising results comparable to human accuracy on the challenging task of interstitial lung disease (ILD) segmentation demonstrate the potential benefits of our approach. Furthermore, the obtained performance outreached redundant conventional networks. Finally, the proposed approach addresses most of the limitations of existing neural network approaches. First, it requires fairly small number of training examples due to the reduced diversity of observations once mapped to a common anatomy. Second, it performs data augmentation in a natural manner thanks to the elastic mapping between observations and representative anatomies. Third, it inherits robustness, stability and better generalization properties for two reasons: the limited complexity of observations after mapping, and the "anatomically" consistent redundancy of the networks.

Methodology

The method consists of two main parts, a transformer and a de-transformer part. The former maps a sample S to N different atlases A i , i ∈ {1, . . . , N }, constructs their warped versions, and trains N different networks, while the latter projects back the N predictions to the initial space. These projections are then combined to obtain the final segmentation. The transformer part consists of a non-linear deformable operator (transformer T i ) and a segmentation network C i while the de-transformer part uses the inverse deformable operator (de-transformer T -1 i ) to map everything back to the initial space of a sample S. The framework is flexible, enables any suitable transformation operator (with an existing inverse) to be coupled with a classifier. The inference framework is illustrated in Fig. 1.

Multimetric Deformable Operator

The multimetric deformable operator, responsible for mapping samples to different anatomies (atlases) therefore reducing variance and producing anatomically meaningful results, is an elastic image registration method that follows a context-driven metric aggregation approach [START_REF] Ferrante | Deformable Registration Through Learning of Context-Specific Metric Aggregation[END_REF] which aims to find the optimal combination of different similarity metrics. The elastic operator is implemented using a deformable mapping from a source image S to a given atlas A i . Let us consider that a number of metric functions ρ j , j ∈ {1, . . . , k}, can be used to compare the deformed source image and the target anatomy A i . The non-linear transformer T corresponds to the operator that optimizes in the domain Ω the following energy:

E( T ; S, A i ) = Ω k j=1 w j ρ j (S • T , A i )dΩ + α Ω ψ( T )dΩ
where w j are linear constraints factorizing the importance of the different metric functions, and ψ() is a penalty function acting on the spatial derivatives of the transformation as regularization to impose smoothness. Such a formalism can be considered either in the continuous setting that requires differentiable functions with respect to the metric functions ρ j or in a discrete setting. The advantage of a discrete variant is that it can integrate an arbitrary number and nature of metric functions as well as regularizers while offering good guarantees concerning the optimality properties of the obtained objective function. Inspired by the work done in [START_REF] Glocker | Dense image registration through mrfs and efficient linear programming[END_REF] we express the non-linear operator as a discrete optimization problem acting on a quantized version of the deformation space.

We used free form deformations as an interpolation strategy, invariant to intensity image metrics, pyramidal implementation approach for the optimization and belief propagation for the estimation of the optimal displacement field in the discrete setting. Details on the implementation can be found in [START_REF] Glocker | Dense image registration through mrfs and efficient linear programming[END_REF].

Segmentation Networks

The segmentation networks C i operate on the mapped image, T i (S), to produce a segmentation map and can be the same or different depending on the task and the application and are completely independent of the exact classifier. After defining the optimal deformations T i , i = 1..N , between the source image and the different atlases in the transformer part, AtlasNet uses the inverse transformations to project back to the initial space of the source image S the predicted segmentation maps:

S seg i = T -1 i (C i (T i (S)
)) In this work, motivated by the state-of-the-art performance of F-CNNs in several problems we adapted them for dense labeling. We use the SegNet deep learning network [START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for image segmentation[END_REF] which performs pixelwise classification and is composed of an encoder and a decoder architecture and follows the example of U-net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF]. It consists of 13 layer groups, similar to the ones of the VGG16 network.The entire architecture consists of repetitive blocks of convolutional, batch normalization, rectified-linear units (ReLU) and indexed max-pooling layers. For more details we refer to the original publication.

Different fusion strategies can be used for the combination of the segmentations. We used the probabilistic output of the classifiers (before hard decision) and fused the output of the different networks based on majority voting.

Implementation Details

For the registration, we used the same parameters for all images and all atlases. Three different similarity metrics have been used, namely, mutual information, normalized cross correlation and discrete wavelet metric. For the mutual information 16 bins were used, in the range of -900 to 100.

We used the same parameters for training all SegNet networks (initial learning rate = 0.01, decrease of learning rate = 2.5•10 -3 every 10 epochs, momen-tum=0.9, weight decay =5•10 -4 ). The training of a single network required around 16 hours on a GeForce GTX 1080 GPUs, while the prediction for a single (volumetric) dataset lasted only a few seconds.For data augmentation we performed only random rotations (between -10 and 10 degrees) and translations (between 0 and 20 pixels per axis) avoiding local deformations since the anatomy should not artificially change. Moreover, for training, we performed median frequency balancing [START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for image segmentation[END_REF] to balance the data, as the samples with disease are considerably fewer than the rest of the samples.

Experimental Results & Dataset

We used as case study to evaluate our method the ILD segmentation in CT images because it is a challenging problem; boundaries are difficult to detect and delineation suffers from poor-to-moderate interobserver agreement [START_REF] Camiciottoli | Lung ct densitometry in systemic sclerosis: Correlation with lung function, exercise testing, and quality of life[END_REF]. Moreover, although several visual scoring systems have been proposed to quantify the disease, they only allow basic quantification of ILD severity. The dataset includes 17 (volumetric) CT images consisting of 6000 slices in total, each being of 512 × 512 dimension, and annotations of lung and disease. The ILD annotation was performed by a medical expert by tracing the disease boundaries in axial view over all slices and used for training the classification model. Assessment of the method was performed on images from 29 additional patients being fully annotated only on selected CT slices (n = 20) by three different observers. Note that the data were multi-vendor (GE & Siemens) and correspond to the same moment of the respiratory cycle.

For all experiments we used six different atlases and registered both training and testing images to them. The choice of atlases was made by a radiologist towards integrating important variability of the considered anatomies. Our experimental evaluation has two objectives: The first one is to show that AtlasNet provides more robust and accurate solutions compared to conventional networks and the second is to examine whether the proposed methodology can truly be trained with fewer examples while leading to good performance. We used five metrics, namely sensitivity, precision, Hausdorff, average contour distance and dice coefficient (over the number of epochs), to evaluate the performance of the proposed method.

On the number of atlases: Fig. 2a presents the behavior of our method using different number of atlases. It can be observed that the dice initially increases and tends to stabilize for more than 5 templates.

Note that, even with the use of only one atlas the deformable operator of AtlasNet helps to increase the dice coefficient (from 0.533 to 0.604), as indicated by Fig. 2b and achieves the highest values of dice compared to conventional networks and usual data augmentation techniques.

On the number of training samples: To evaluate the performance of our architecture with less samples we used a reduced number of samples (30%, 50% and 70% respectively) for the same number of epochs (18) and compare the performance with the one in [START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for image segmentation[END_REF]. The obtained mean dice coefficient values in [START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for image segmentation[END_REF] were 0.434, 0.462, 0.487, while for AtlasNet were 0.613, 0.646 and 0.672 respec- tively, indicating the robustness of AtlasNet with a significantly lower number of samples. In simple words, the proposed architecture produces better or similar results with 30% of the samples compared to the state-of-the-art architecture [START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for image segmentation[END_REF] with and without data augmentation.

Comparison with the state-of-art: Although results on different datasets are not directly comparable, we compare our method with works related to ILD segmentation. Anthimopoulos et al. [START_REF] Anthimopoulos | Lung Pattern Classification for Interstitial Lung Diseases Using a Deep Convolutional Neural Network[END_REF] classified CT image patches with ILD patterns using a CNN and obtained 0.856 accuracy for 6 disease classes. By extracting patches on our data (where different patterns are annotated as a single class) in the same way as in [START_REF] Anthimopoulos | Lung Pattern Classification for Interstitial Lung Diseases Using a Deep Convolutional Neural Network[END_REF] we obtained 0.916 mean accuracy. In [START_REF] Vakalopoulou | Deep patch-based priors under a fully convolutional encoder-decoder architecture for interstitial lung disease segmentation[END_REF] a patch-based CNN was augmented with a deep encoder-decoder to exploit partial annotations. By applying AtlasNet on the same dataset as in [START_REF] Vakalopoulou | Deep patch-based priors under a fully convolutional encoder-decoder architecture for interstitial lung disease segmentation[END_REF], we increased the mean dice from 0.671 to 0.725.

Moreover we compared AtlasNet in respect to disease segmentation with standard frameworks (without registration and with or without data augmentation) for the same number of epochs (18) and illustrate results in Table 1 and Fig. 2b. For equal comparison, we assessed accuracies using the same classification strategy [START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for image segmentation[END_REF] trained on the initial CT slices, and after performing data augmentation as described earlier. The proposed method reports the best accuracy with respect to Hausdorff distance, average contour distance and dice, indicating that the disease segmentation is much more accurate than by the conventional frameworks with or without data augmentation. This can be inferred also from Fig. 3 where axial slices of two different subjects are depicted. It is clear that the proposed approach segments accurately the boundaries of the disease.

For a more complete evaluation, we compare AtlasNet also with inter-observer agreement using the annotations of three different medical experts. In particular, the annotations of one observer have been used as ground truth to evaluate the rest. From Table 1 and Fig. 3, it can be observed that AtlasNet demonstrates more robust performance than manual segmentation. Finally, it is worth mentioning that even if the network operates on 2D slices, without accounting for out-of-slice connections, the fusion of the different atlases' predictions makes the final segmentation smooth across all three axes.

Concerning the computational resources, we use a single segmentation network [START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for image segmentation[END_REF] for each of the N atlases, therefore the time and memory usage for one atlas is that of the CNN, while we also showed that a small N (such as 6) suffices. For segmentation of one volumetric CT on a single GPU the total testing time (using 6 atlases) is 3-4min, which includes the registration step. The registration cost is negligible since a graph-based GPU algorithm is used taking 3-5sec per subject. This cost drops linearly with the number and computing power of GPUs. Thus we believe that the additional complexity of AtlasNet is fully justified, since it improves performance by more than 20% and also maintains it stable with only 30% of the training data compared to conventional single networks.

Conclusion

In this paper, we present a novel multi-network architecture for (healthy or pathological) tissue or organ segmentation that maximizes consistency by exploiting diversity. Evaluation of the method on interstitial lung disease segmentation highlighted its advantages over previous competing approaches as well as inter-observer agreement. The investigation of techniques for soft parameter sharing of hidden layers, and information transfer between the different networks and atlases is our direction for future work. Finally, the extension to multi-organ segmentation including multiple classes loss functions is one of the potential directions of our method.

Fig. 1 :

 1 Fig. 1: The proposed AtlasNet framework. A i indicates atlas i.

Fig. 2 :

 2 Fig. 2: Quantitative evaluation of the dice coefficient for the proposed method and different number of atlases.

( a )

 a Expert 1 (b) Expert 2 (c) SegNet [2] (d) Augm. & [2] (e) AtlasNet

Fig. 3 :

 3 Fig. 3: Interstitial lung disease segmentation (depicted with red color) on two testing subjects using the different employed strategies.

Table 1 :

 1 Different evaluation metrics for the testing dataset for scleroderma disease.

	Method	Sensitivity Precision Hausdorff Dist. Average Dist. Dice
	SegNet [2]	0.348	0.623	4.984	1.891	0.533
	Augmentation & SegNet [2]	0.534	0.567	4.077	1.309	0.619
	Inter-observer	0.693	0.522	4.005	1.317	0.662
	AtlasNet	0.682	0.545	3.981	1.274	0.677

Acknowledgements

This work has been partially supported by the european project ERC-PoC 737604 TheraPanacea.