Tolerance Analysis of Over-Constrained Assembly Considering Gravity Influence: Constraints of Multiple Planar Hole-Pin-Hole Pairs

Over-constrained assembly of rigid parts is widely adopted in aircraft assembly to yield higher stiffness and accuracy of assembly. Unfortunately, the quantitative tolerance analysis of over-constrained assembly is challenging, subject to the coupling effect of geometrical and physical factors. Especially, gravitywill affect the geometrical gaps inmechanical joints between different parts, and thus influence the deviations of assembled product. In the existing studies, the influence of gravity is not considered in the tolerance analysis of over-constrained assembly.This paper proposes a novel tolerance analysis method for over-constrained assembly of rigid parts, considering the gravity influence.This method is applied to a typical over-constrained assembly with constraints of multiple planar hole-pin-hole pairs. This type of constraints is non-linear, which makes the tolerance analysis more challenging. Firstly, the deviation propagation analysis of an over-constrained assembly is conducted. The feasibility of assembly is predicted, and for a feasible assembly, the assembly deviations are determined with the principle of minimum potential energy. Then, the statistical tolerance analysis is performed. The probabilities of assembly feasibility and quality feasibility are computed, and the distribution of assembly deviations is estimated. Two case studies are presented to show the applicability of the proposed method.


Introduction
Mechanical parts are always manufactured with more or less geometrical deviations. Tolerance analysis aims to predict the influence of these manufacturing deviations on the quality level of a mechanical system at the design stage. The gaps in mechanical joints between different parts, coupled with the manufacturing deviations, impact the functionality of an assembly system. The manufacturing deviations are assumed to be independent random variables defined by their tolerance zones and the probability distributions with respect to manufacturing process, while the gaps are characteristic of epistemic uncertainty which makes the tolerance analysis complex [1]. The gaps are seldom considered in iso-constrained assembly, but in over-constrained assembly they allow the parts to be assembled without deformation and interference, and cannot be ignored. Over-constrained assembly of rigid parts is widely used in products, such as aircrafts and automobiles, because this kind of assembly can increase stiffness and accuracy of mechanical structure [2]. However, the feasibility of assembly is weakened due to the introduction of redundant constraints into the assembly, and the deviations of assembled products are difficult to predict because of the uncertain gaps. This is one of the key challenges in assembly tolerance analysis.
Some studies have been conducted on tolerance analysis of over-constrained assembly. Compared to worst case method, statistical tolerance analysis is regarded to be better to understand the behavior of a mechanical system, and the over-constraints have to be considered for accuracy in the tolerance analysis of over-constrained assembly [3]. Dantan and Qureshi [4] introduced the quantifier notion, including existential quantifier "there exists" and universal quantifier "for all", to describe the condition corresponding 2 Mathematical Problems in Engineering to a functional requirement. The mathematical formulation of tolerance analysis was computed with quantified constraint satisfaction problem solvers and Monte Carlo simulation (MCS). With the mathematical formulation based on quantifier notion, Qureshi et al. [5] acquired the probability of the functional operation of an assembly based on MCS and optimization technique, which is used to find the worst case for gap configuration. This approach was generalized and applicable to both explicit and implicit response function, as well as both linear and non-linear constraints. For explicit function and linear constraints, MCS is time consuming and computationally expensive. Ballu et al. [3] solved probability of failure by using first order reliability method (FORM) and second-order reliability method (SORM) to determine the feasible extreme gap configuration. Beaucaire et al. [6,7] proposed the mathematical formulation from the structural reliability domain considering several contact point situations, and computed the defect probability by FORM. For non-linear constraints, Dumas et al. [8] investigated different linearization strategies and discussed the impact of these strategies on the predicted probability of failure. In addition, Dumas et al. [1] formulated the tolerance analysis problem by using the Lagrange dual form of optimization, where a simulation-based selective search algorithm, i.e., a genetic algorithm, was developed and then the probability of failure was estimated by FORM. The above studies [1,[3][4][5][6][7][8] focused on computing the probability that the assembled product can fulfill its intended functional requirement.
Another problem related to statistical tolerance analysis is the distributions of assembly deviations of functional points, where these deviations characterize the fulfillment of functional requirement of assembled product. In overconstrained assembly, the presence of gaps in mechanical joints makes the deviations of assembled product difficult to predict, because a large variety of possible gap configurations exist during the assembly process. Franciosa et al. [9] determined the best-fit assembly configuration by using optimization algorithm, where the optimization objective is to minimize the sum of the squares of the distances between object plane/axis and target plane/axis respectively. However, this optimization objective is not suitable for most of practical assemblies in production. Li et al. [10] defined two planar hole-pin pairs as 3-2-1 location mode, and expressed the ranges of location deviations by the intersection of inequations which is approximated to an elliptic. However, this study did not develop an approach to get the certain location deviations.
Furthermore, tolerance analysis of over-constrained assembly is subject to many nongeometrical factors involved in an assembly process, such as thermal flux, contact force, gravity, etc. Pierre et al. [11] integrated thermomechanical strains into variational tolerance analysis. Gouyou et al. [2] discussed the part deformation subject to contact forces which allow the parts to be assembled when the interference happens. However, although the influence of gravity on the assembly deviations of rigid parts has been well recognized [12], there is still limited research in tolerance analysis of overconstrained assembly integrating the influence of gravity. Specifically, in aircraft assembly, planar hole-pin-hole joints are widely used to join different parts and locate parts on the fixtures. Because many parts and fixtures in aircraft assembly are large and heavy, the gravity will significantly impact the configuration of gaps in hole-pin-hole joints, and thus the assembly accuracy. In view of the strict dimensional accuracy requirement in aircraft manufacturing, a reliable prediction of assembly deviations by considering the gravity influence is especially important to realize effective dimensional management in aircraft assembly. This paper proposes a method for tolerance analysis of over-constrained assembly incorporating the influence of gravity on assembly deviations. This method will be applied to a typical over-constrained assembly with constraints of multiple planar hole-pin-hole pairs. First of all, the deviation propagation throughout an over-constrained assembly is analyzed. Here, the feasibility of assembly is predicted; for a feasible assembly, the assembly deviations are determined considering the gravity influence. Then, the statistical tolerance analysis is performed by using MCS, and its output results include the probabilities of assembly feasibility and quality feasibility, as well as the distributions of assembly deviations. Especially, in the deviation propagation analysis, the judgment of assembly feasibility is formulated as a search problem. And for a feasible assembly, with the principle of minimum potential energy, the calculation of assembly deviations is modeled as an optimization problem. These search and optimization problems are inefficient and difficult to solve, because of the non-linear studied constraints of planar hole-pin-hole pairs. These problems are solved by using the search algorithms including sequential search and binary search, and the optimization algorithms including particle swarm optimization (PSO) and linear programming (LP).
The paper is organized as follows. Section 2 gives problem definition and mathematical representation. Sections 3 and 4 present deviation propagation analysis and statistical tolerance analysis of over-constrained assembly, respectively. Section 5 shows case studies. Section 6 draws conclusions.

Problem Definition and Mathematical Representation
. . Problem Description and Assumption. Tolerance analysis consists in predicting the functionality of an assembly system based on the consideration of the manufacturing deviations of parts. For over-constrained assembly of rigid parts, tolerance analysis consists of two key problems to be studied, including assembly feasibility prediction and assembly deviation calculation. In over-constrained assembly of rigid parts, the feasibility of assembly is weak because of redundant constraints, even with the gaps in mechanical joints between different parts. Furthermore, due to the uncertainty of gaps, the assembly deviations are difficult to predict, especially coupled with non-linear constraints and non-geometrical factors.
With considering the influence of gravity, the proposed method of tolerance analysis of over-constrained assembly is illustrated by a typical assembly with constraints of multiple Figure 1: Two parts to be assembled with constraints of multiple planar hole-pin-hole pairs. planar hole-pin-hole pairs. This type of constraints possesses the feature of non-linearity. As shown in Figure 1, each of the two parts, denoted as Part 1 and Part 2, has two or more pairs of mating holes on the mating plane. Part 1 is assumed to be fixed, and Part 2 will be assembled onto Part 1 with the planar hole-pin-hole joints. The centers of the mating holes of Part 1 are denoted as B 1 ∼B . Similarly, these centers of Part 2 are denoted as C 1 ∼C . N is the number of hole-pinhole joints. Furthermore, the quality requirement concerns the deviations of Part 2 with respect to (w.r.t.) Part 1. The functional points are on Part 2, denoted as K 1 ∼K . N is the number of functional points. The deviations of the functional points are used to represent the assembly deviations. In addition, G denotes the center of gravity of Part 2.
Notably, for simplifying the analysis and highlighting the research focus, some assumptions are made as follows: (a) The assembly tolerance analysis is performed only considering the degrees of freedom (DOFs) on plane, such that the mating planes are assumed to be perfect.
(b) The parts are treated as rigid bodies, and deformation and interference are not allowed.
(c) The direction of gravitational field is on the mating planes.
(d) The position of gravity center of Part 2 is not affected by its manufacturing deviations.
(e) The hole-pin-hole pairs adopt clearance fit.
(f) For the manufacturing deviations of hole-pin-hole pairs, dimensional deviations of diameters of holes and pins, and position deviations of centers of holes are considered, but the shape deviations are ignored.
. . Mathematical Representation . . . Definition of Coordinate Systems. In order to express the deviations as well as their relationships, four types of coordinate systems (CSs) are defined, including global coordinate system (GCS), part coordinate system (PCS), gravitational field coordinate system (GFCS), and functional point coordinate system (FPCS).
GCS is considered to be error free and unchangeable during the entire assembly process. PCS is associated with a part. Similarly to GCS, GFCS and FPCS are also unchanged. GFCS is constructed according to the direction of gravitational field, and the -axis is in the opposite direction to gravitational field. Each FPCS is used to represent the deviations of a functional point. The origin of a FPCS is in the nominal position of a functional point w.r.t. GCS, and the axes correspond to the direction required to be measured. Moreover, the spatial relationships among different CSs can be described by homogeneous transformation matrix (HTM) [13].
As shown in Figure 2, Part 2 without any deviation is in the nominal position w.r.t Part 1. XOY is the GCS, which is assumed to coincide with the PCS of Part 1. xoy is the PCS of Part 2 denoted as PCS 2 .
. . . Tolerance Representation. The dimensional and geometric tolerances indicate the deviation of an actual feature from its nominal one. The features are allowed to rotate and translate within their tolerance zone [10], and these rotations and translations can be described by the small displacement torsor (SDT) model [14]. In SDT model, tolerance zones and feature deviations are defined as vectorial representation [15], i.e., = [ ] T . = [ ] T denotes the rotation vector, which has three components , , and around x, y, and z axes, respectively. Similarly, T denotes the translation vector, which has three components u, v, and along x, y and z axes, respectively. In this paper, it is assumed that The assembly deviations are caused by deviation sources which are propagated and accumulated through the assembly process. This data stream of deviations can be analyzed by deviation propagation model [16][17][18]. For example, the deviation of feature i is = [ V ] T . An arbitrary point on feature i is A, and its coordinates are denoted as a = where the HTM T is ] . (2)

Deviation Propagation Analysis
Deviation propagation analysis studies the introduction, propagation and accumulation of data stream of deviations through an assembly process. For an over-constrained assembly with constraints of multiple planar hole-pin-hole pairs, the deviation propagation analysis is performed with considering the influence of gravity, which contains assembly feasibility judgment and assembly deviation determination.
. . Assembly Feasibility Judgment. An assembly with over constraints is usually disturbed by the infeasibility of assembly, even though normal joint gaps are reserved. Dantan and Qureshi [4] introduced the existential quantifier "there exists" to correctly formulate this problem concerning assembly with gaps, i.e., if there exits at least one gap configuration such that the assembly is feasible, then assembly feasibility of the mechanical structure is ensured. The mathematical expression is where M denotes the manufacturing deviations of parts, g denotes the gaps and g max is the allowable gap widths, f , = 1, 2, . . . , , are the functions of assembly requirement, and N is the number of the functions. All of these functions have to be non-positive to ensure the feasibility of assembly.
According to the above formulation, the judgment of assembly feasibility can be transformed to a search problem. Moreover, this search problem is inefficient to be solved due to the non-linear studied constraints of multiple planar holepin-hole pairs. Firstly, the search problem is modeled. Then in order to improve the efficiency, a solving method is given by reducing dimensionality of search object and narrowing down search space.
. . . Formulation of Assembly Requirement. In overconstrained assembly, the gaps in mechanical joints are correlated with each other and their worst configuration depends on the manufacturing deviations of parts [3]. For the constraints of multiple planar hole-pin-hole pairs, the correlation of gaps, coupled with the relationships between gap configurations and part manufacturing deviations, is complicated. In an assembly process, the deviation of PCS 2 w.r.t GCS is denoted as 2 = [ 2 2 V 2 ] T , and one 2 corresponds to a configuration of gaps. As shown in Figure 1, the coordinates of B , = 1, 2, . . . , , y' x' With the deviation 2 , the coordinates of C , where T 2 is the HTM of 2 . The functions of assembly requirement concern the planar hole-pin-hole joints. As shown in Figure 3, by taking the ith constraint of planar hole-pin-hole pair as an example, the function of assembly requirement is where = (1/2)( + − 2 ) and d , d , and d denote the diameters of the two mating holes and the pin respectively. It is definite that ≥ 0, due to clearance fit of hole-pin-hole pairs. According to the assembly requirement as seen in (3), the judgment of assembly feasibility is to judge whether f is non-positive, i.e., to judge whether | → B C |, is not larger than s , or whether | → B C | 2 is not larger than 2 . Therefore, from the view of judging assembly feasibility, (6) can be translated to (7) aŝ=

Mathematical Problems in Engineering
wherêis decorated with a topscript " ∧ ", in order to distinguish from the original assembly requirement function f . By referring to (5), (7) can be further arranged aŝ For coordinates and diameters above mentioned, the nominal values are decorated with a superscript "0, " whereas they represent the actual ones with manufacturing deviations of parts. In (8), coupled with the correspondence between 2 and configuration of gaps, the part manufacturing deviations and gap configuration are both incorporated.
According to (3) and (8), the assembly requirement under the constraints of N planar hole-pin-hole pairs is formulated as: if there exits at least one real vector 2 satisfying all the inequalitieŝ≤ 0, = 1, 2, . . . , , then assembly feasibility of the mechanical structure is ensured. The mathematical expression is where R 3×1 denotes three-dimensional real vector space.
. . . Solving Method of Search Problem. The judgment of assembly feasibility depends on solving the search problem modeled by (9). As a common search algorithm, sequential search algorithm is used. The search object is 2 = T with three dimensions. The search space is difficult to narrow, due to the non-linearity of the functionŝ, = 1, 2, . . . , . Therefore, the solving of this search problem is inefficient because of the multi-dimensional search object and the large search space. To improve the efficiency, the dimensionality of search object is reduced into one, and the search space is reasonably narrowed.
The feasible solution space of the search object 2 can be shown in a Cartesian coordinate system 2 − 2 V 2 2 . Referring to Figure 4, the feasible solution space for an inequalitŷ≤ 0 is denoted as D . When 2 is constant, D is a certain circle, and its radius always equals to . Besides, the locus of center of this circle with 2 is a right-hand circular helix. Then the feasible solution space under constraints of N planar hole-pin-hole pairs, denoted as D , can be expressed as the intersection of the spaces D , = 1, 2, . . . , . Figure 5 shows an example of non-empty D when N=2. When 2 is constant, D represents the intersection of N circles with the known centers and radii. The situation of this intersection denoted as ( 2 ) can be judged with a geometry algorithm [21] especially when N>2. It is assumed that when the intersection is empty, ( 2 ) < 0; when the intersection is only a point, ( 2 ) = 0; otherwise, ( 2 ) > 0. Then (9) can be redefined as follows: if there exits at least one real 2 such that the intersection of N circles is non-empty, then assembly feasibility of the mechanical structure is ensured. The mathematical expression is where D 2 = { 2 | 2 ∈ [ 2 min , 2 max ]} denotes the search space of 2 , which is determined as follows: (a) The ranges of 2 under constraints of each two planar hole-pin-hole pairs are successively calculated. For example, as shown in Figure 5, the range of 2 is denoted as 2 ∈ [ 2 12 min , 2 12 max ]. According to geometric characteristic of D 1 and D 2 , 2 12 min , and 2 12 min correspond to the situations when the two circleŝ1( 2 12 min,max , 2 , V 2 ) ≤ 0 and̂2( 2 12 min,max , 2 , V 2 ) ≤ 0 are externally tangent as If (11) has no real solutions, this range of 2 is empty.
(b) The intersection of these ranges of 2 in (a) is calculated as D 2 .
As shown in Figure 6, if D 2 is empty, it is sure that the assembly is infeasible; but if D 2 is non-empty, the feasibility of assembly is still not ensured and the search algorithm has to be performed, unless N=2. . . Assembly Deviation Determination. Assembly deviations are difficult to quantify in over-constrained mechanical structure due to uncertainty of gaps, when the assembly is feasible. Here, the influence of gravity is considered. According to the principle of minimum potential energy, the gravitational potential energy of Part 2 tends to be minimum. In other words, the center of gravity of Part 2 drops along the direction of gravitational field as large as possible, within the constraints of satisfying assembly requirement. Therefore this section involves an optimization problem. Furthermore, the non-linear studied constraints make the optimization difficult to solve. Firstly the optimization problem is modeled. Then similarly to Section 3.1, a solving method is given by reducing dimensionality of optimization object and narrowing search space, to improve both accuracy and efficiency.
Finally, the assembly deviations of functional points are calculated, and the quality feasibility is judged.
. . . Optimization Problem Modeling. As shown in Figure 1 respectively, which can be calculated similarly to (4)-(5), as  Then with the deviation 2 , the coordinates of G w.r.t.
Here, -axis of GFCS has the angle from Y-axis of GCS. Then (14) can be arranged as where is the coordinate of G against the direction of gravitational field.
When the influence of gravity is considered, is minimum within the constraints of satisfying assembly requirement, according to the principle of minimum potential energy. The optimization model is constructed as where 2 is the optimization object, ( 2 ) is the objective function, and 2 ∈ D is the constraint condition. Without considering the influence of gravity and other physical factors, the deviation 2 is randomly distributed within the feasible solution space D . When the influence of gravity on the assembly deviations is taken into account (the influence of other physical factors is not considered), the value of 2 corresponds to the optimal solution of (16). Figure 7 shows an example of the assembly under constraints of two planar hole-pin-hole pairs. With considering the gravity influence, the center of gravity G on Part 2 moves along the direction of gravitational field to the farthest position allowed by the assembly constraints. . . . Optimization Problem Solving. The assembly deviations correspond to the optimal solution of (16). As a selective search algorithm, PSO algorithm is used to solve the optimization problem. This optimization solving is difficult due to the non-linear constraint 2 ∈ D . Here, the dimensionality of search object is reduced, and the search space is calculated, to improve the accuracy and efficiency of solving.
As seen in (16), when 2 equals a real constant 2 con , the constraint 2 ∈ D ( 2 con , 2 , V 2 ) represents the intersection of N known circles by referring to Section 3.1, and the objective function ( 2 con , 2 , V 2 ) is regarded as a line. The slope of this line is unchangeable, which is only related to the direction of gravitational field. The line equation denoted as l(u 2 ,v 2 ) is where = −sin , = cos , = (− sin + cos )(cos 2 con − 1) + ( sin + cos ) sin 2 con . Therefore, the optimization modeled by (16) can be treated a LP problem when 2 = 2 con . Figure 8 gives a graphical schematic of constraint condition and objective function when N=3. Here, the minimum value of ( 2 con , 2 , V 2 ) 8 Mathematical Problems in Engineering is denoted as min ( 2 con , 2 , V 2 ). The coordinates of the point corresponding to min ( 2 con , 2 , V 2 ) are denoted as P min = [ 2 min V 2 min ] T . Obviously, this point belongs to intersection points of the arcs which comprise boundary of the domain D ( 2 con , 2 , V 2 ), or tangency points of the line l(u 2 ,v 2 ) and these arcs. Furthermore, it is certain that P min depends on 2 con , i.e., 2 min = 2 ( 2 con ), V 2 min = V 2 ( 2 con ). Hence, (16) can be remodeled with a nested LP as where the search space , which is calculated as follows: As seen in Section 3.1, a feasible solution of (10), denoted as 2 , has been found to ensure the assembly feasibility of the mechanical structure. It is obvious that 2 ∈ [ 2min , 2max ]. Moreover, for the search space of (10) denoted as D 2 = { 2 | 2 ∈ [ 2 min , 2 max ]}, it is definite that D 2 ⊇ D 2 , and then 2 min ≤ 2min ≤ 2 and 2 ≤ 2max ≤ 2 max . 2min and 2max can be approximately calculated by binary search algorithm. Especially, when N=2, 2min = 2 min and 2max = 2 max , thus the binary search algorithm is omitted.
The quality requirement is denoted as QR = , . An assembly is regarded to be qualified when all the assembly deviations of functional points meet the quality requirement, i.e., ∈ QR and ∈ QR , = 1, 2, . . . , .

Statistical Tolerance Analysis
Because as deviation sources, manufacturing deviations of parts are generally random deviations within their tolerance zones and statistical analysis is necessary rather than a single case analysis. As one of the most common methods for statistical analysis, the MCS is often computationally expensive and time-consuming, but quite comprehensive and easy to use even for implicit and non-linear problems. The MCS with deviation propagation analysis in Section 3 is used to perform statistical tolerance analysis. The flowchart of the statistical analysis is shown in Figure 9. Iteratively, a set of random deviations of the sources are generated within their tolerance zones; then the assembly feasibility is judged; for the feasible assembly, the assembly deviations of functional points are calculated and the quality feasibility is judged. The process of random generation, judgment and calculation is iterated until the samples are enough. As to random generation of deviations within their tolerance zones, it is assumed that the deviations are independent of one another and present normal distributions, and the deviations and their tolerances have a relationship as = 6 . After all the iterations, the probabilities of assembly feasibility and quality feasibility are estimated as [8] where N is the maximum number of iterations of the MCS and , (M) is an indicator function. The functions for assembly requirement and quality requirement are Furthermore, for feasible assemblies, the deviations distributions of functional points can be quantified, where the frequency distributions and the statistical characteristics are computed, seen in Section 5.
For the statistical tolerance analysis, the input is deviation sources, and the output contains probabilities of assembly feasibility and quality feasibility, as well as deviations distributions of functional points for feasible assemblies.

Case Studies
The proposed method is applied to two cases. Firstly stability and accuracy of the optimization seen in Section 3.2 are verified. Then the output results of statistical tolerance analysis for the gap configuration influenced by gravity are compared with the output results for random configuration and worst configuration of gaps.
The computational works are implemented by MATLAB. For sequential and binary search algorithms seen in Section 3, the precision is set as 10 −7 for 2 . In the statistical tolerance analysis with MC simulation seen in Section 4, the number of iterations in an experiment is set as N =10000.

. . Case Description
Case . This case is from a fixture with boards as shown in Figure 10, which is widely used for locating aircraft panels in the panel assembly process. The small working surfaces on a fixture board are applied to constraining the normal displacement of the panels. Each fixture board is mounted on fixture frame with two planar hole-pin-hole joints. This case is the assembly of a fixture board and a fixture frame, denoted as Part 1 and Part 2. Each small working surface of Part 2 is abstracted as a functional point.  . . Stability and Accuracy of Optimization (PSO). As seen in Section 3, PSO is used to determine assembly deviations, which further impacts the output results of statistical tolerance analysis. In this PSO algorithm, the swarm size is 10, and the number of iterations is 20. PSO is a heuristic optimization algorithm, with simple parameters and fast convergence, but it is easily trapped into local optimum. Therefore, it is necessary to verify stability and accuracy of this optimization algorithm.
Firstly, the experiment of statistical tolerance analysis is performed and repeated 10 times, and the 10 sets of output results are observed and compared. In order to eliminate the disturbance of different generated deviation sources, the same N sets of deviation sources are required for all experiments.
Part 2 Figure 11: CSs of a fixture frame and a fixture board.   Figure 13: CSs, Dimensions and tolerances of two fuselage sections.  Three iterative sub-processes are separated from the experiment process. The iterative sub-process of deviation sources generation is denoted as Sub-process 1, and N sets of generated deviation sources are stored in a matrix H . Then with these deviation sources, the iterative sub-process of assembly feasibility judgment is denoted as Sub-process 2, and N ( ≤ ) sets of deviation sources, which can make sure the feasibility of assembly, are selected and stored in a matrix H . Finally, with the deviation sources from H , the iterative sub-process of assembly deviation calculation and quality feasibility judgment is denoted as Sub-process 3. After the three sub-processes, the output results can be calculated.
Therefore, the 10 experiments are designed as: Subprocesses 1 and 2 are implemented in advance, and then Sub-process 3 is repeated 10 times. The 10 sets of output results with 4 figures after the decimal point are all the same, as shown in Table 1 for Case 1 and Table 2 for Case 2, where deviations of functional points (denoted as DFP) are classified into the deviations satisfying quality requirement (denoted as DFP QR) and that not satisfying quality requirement (denoted as DFP NQR).
For further verification, with a set of deviation sources arbitrarily selected from H , the PSO procedure is performed and repeated 10 times, and the 10 sets of optimization results are compared. Tables 3 and 4 show the selected deviation sources for Cases 1 and 2, respectively. The corresponding feasible solution space D is shown in Figure 14 for Case 1 and Figure 15 for Case 2. In the 10 PSO procedures, the global best particles 2 and the corresponding -coordinates are shown in Figure 16 for Case 1 and Figure 17 for Case 2. It is observed that the 10 optimization procedures for each case are all convergent, and the 10 sets of optimization results, including global best particles and objective function values, are all the same when these values are rounded to 6 decimal places.
Therefore, Tables 1 and 2 and Figures 16 and 17 indicate this PSO algorithm is quite robust and accurate, and thus the output results of tolerance analysis by using the proposed method are reliable.

. . Comparison with Random Configuration and Worst Configuration of Gaps.
The proposed method considers the influence of gravity on the gap configuration when determining the assembly deviations. In order to exhibit the significance of incorporating the influence of gravity to tolerance analysis, two comparison experiments with common gap configurations are given. The first comparison experiment is denoted as EX1, where the deviation 2 is randomly generated within the feasible solution space based on normal distribution. The second comparison experiment, denoted as EX2, uses the worst configuration of gaps, which makes the assembly deviations largest. Besides, the experiment considering the gravity influence is denoted as EX0.
In EX 2, the determination of the worst configuration of gaps is an optimization problem. The modeling and solving of this optimization are both similar to (16)- (19). For example, worst gap configuration corresponds to the largest assembly deviation of the jth functional point along Y -axis. The optimization is modeled as seen in (23), and solved by PSO with a nested LP.
In EX2 of Case 1, it is assumed that the worst configuration of gaps corresponds to the largest assembly deviations of K 2 along Y 2 -axis. In EX2 of Case 2, the worst configuration of gaps is assumed to concern the largest assembly deviations along Y-axis, to make a straight comparison between the output results of EX2 and EX0.
The output results of EX0, EX1 and EX2 are compared. In order to avoid the distractions of different generated deviation sources, the three experiments, denoted as a group of experiments, are designed as follows: Sub-processes 1 and  2 are implemented in advance, and then Sub-process 3 for EX0, EX1, and EX2 is performed, respectively. To enhance the reliability of the predicted output results, a group of experiments are repeated 10 times, and the 10 groups of output results are averaged, which are shown in Table 5 for Case 1 and Table 6 for Case 2. Furthermore, by taking the output results of a group of experiments as examples, the distributions of assembly deviations are shown in Figure 18 for Case 1 and Figure 19 for Case 2, and the distributions of the corresponding 2 are shown in Figure 20 for Case 1 and Figure 21 for Case 2.
From Tables 5 and 6 and Figures 18-21, the following is observed: (1) As to EX0, except for the output results along X-axis of Case 2, the mean of DFP is unequal to 0. The DFP NQR and the corresponding 2 are unsymmetrically distributed. This may be because of the following: (a) part 2 (PCS 2 ) drops along one constant direction to the boundary of feasible solution space D , where the constant direction is related toaxis; (b) Y -axis, = 1, 2, . . . , , is not perpendicular to -axis.
(2) As to EX1, P is the highest. The DFP NQR and the corresponding 2 are symmetrically distributed. This may be due to the normal distribution of 2 within D . The deviation 2 in the center of D has much more probabilities to be selected than the one around the boundary  of D . The former one is more likely to correspond to the smaller DFP than the latter one.
(3) As to EX2, the P is the lowest. For the DFP corresponding to    boundary of D , where the opposite directions are related to X -axis or Y -axis, = 1, 2, . . . , .
Certainly the output results of EX0 are significantly different from that of EX1 and EX2. For the feasible assembly of over-constrained mechanical structure, gravity has a remarkable impact on gap configuration and then assembly deviations. The gravity influence in tolerance analysis of overconstrained assembly cannot be ignored.

Conclusions
This paper presents a tolerance analysis method for overconstrained assembly considering the influence of gravity. The proposed method contains deviation propagation analysis and statistical tolerance analysis. The probabilities of assembly feasibility and quality feasibility can be computed, and the deviations distributions for feasible assemblies can be quantified.
The detailed tolerance analysis process is illustrated for a typical over-constrained mechanical structure with constraints of multiple planar hole-pin-hole pairs. This type of constraints has non-linearity, which increases complexity of assembly tolerance analysis. In the deviation propagation analysis, the judgment of assembly feasibility and the calculation of assembly deviations are formulated as a search problem and an optimization problem, respectively. Due to the non-linear constraints, it is inefficient and difficult to solve these problems. In order to improve accuracy and efficiency, according to characteristic of this type of constraints, the search and optimization problems are solved by reducing object dimensionality and narrowing search space, where the search algorithms, such as sequential search and binary   search, as well as the optimization algorithms, such as PSO and LP, are used. The results of the case studies show the following: (a) the PSO algorithm for the assembly deviation determination is demonstrated to be quite robust and accurate, thus the output results of the tolerance analysis are reliable and accurate enough; (b) for predicting the assembly deviations of the over-constrained mechanical structure, the influence of gravity is remarkable and cannot be ignored.
Our main contribution lies in developing a reliable and accurate method for assembly tolerance analysis with overconstraints of multiple planar hole-pin-hole pairs, considering the influence of gravity. Besides gravity, it can be used to incorporate the influence of other conservative forces into tolerance analysis of over-constrained assembly as well.
In the proposed method, Monte Carlo simulation (MCS) is used for the statistical tolerance analysis, due to its general applicability for implicit and non-linear problem. The accurate and reliable application of the MCS lightly depends on the number of iterations. The accuracy of prediction can be improved by increasing the number of iterations at the cost of computational time. To relieve the required computational effort, the first-order reliability method (FORM) method combining with a reasonable linearization strategy for the non-linear constraints will enable a reasonably accurate prediction of defect probabilities at a lower computational effort. However, it should be noted that the loss of accuracy of prediction on defect probability would be unavoidable.
Furthermore, the proposed method is applied to the twodimensional tolerance analysis of over-constrained assemblies, but the analytical framework of the proposed method may be extended to three dimension applications, i.e., for over-constrained mechanical structure, the assembly feasibility judgment and the assembly deviation determination considering the influence of gravity can be modeled as a search problem and an optimization problem, respectively. For the three-dimensional case, the DOFs in space instead of the ones on plane are required to be considered, and the analysis of assembly requirement is more challenging. Then, the search and optimization problems for judging assembly feasibility and determining assembly deviations are more difficult to be modeled and solved. The authors are already working on the extension of the method employed in this paper to the threedimensional tolerance analysis of over-constrained assembly. Besides, in future work, part deformation for eliminating the interference between different parts will also be integrated into the assembly tolerance analysis.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.