**1. Introduction**

Under the innovative concept of Industry 4.0, automated and digitized systems in smart factories could enable the real-time integration and analysis of massive amounts of data by the use of electronics and information technologies. Moreover, this would finally result in a more flexible and optimized manufacturing process [1–3]. This theory points to the improvement of intelligent solutions in the tolerance design of mechanical products, including the replacement of empiricism with knowledge-intensive and data-based processes. The general practice of mechanical tolerance design is to treat tolerance design as a combination of tolerance modeling, assembly simulating, tolerance allocation, and optimization. The best combination of part tolerance is assigned through multi-dimensional objective constraints in order to meet environmental clearance, conform to functional requirements, improve product quality, and lower manufacturing costs in the late design stage [4–6]. To increase the accuracy and robustness of the designed tolerance scheme, the influence of processing feature degradation caused by machining precision deterioration is taken into consideration.

In the past few decades, researchers have presented several e fficient computer-aided methods for estimating geometry variation and product tolerance. Many models have been developed for geometrical feature representation, deviation accumulation, and tolerance zone estimation [7]. Conventional tolerance modeling methods are mainly categorized into those based on vector loops, Tolerance-Map (T-Map), small displacement torsor (SDT), a Homogeneous Transformation Matrix (HTM), and some other innovative models. Vector loops provide vectorial tolerance representation and analysis for the surfaces of di fferent components. However, only five standard surfaces are typically considered (plane, cylinder, sphere, cone, and torus) [8–11]. Consisting of a hypothetical volume of points, a T-Map contains all the possible locations for geometrical features, and is used to illustrate the accumulation of product variations [12–14]. The SDT model uses screws and constraints to establish the extreme limits of part tolerance zones [15–17]. To overcome the lack of analytical models between the target feature and each feature on the dimensional chain, HTM was combined with T-Map to obtain an expression of the decoupling pose of tolerance features [18]. Since the Jacobian matrices are quite suitable for deviation propagation [19], the unified Jacobian–Torsor model was developed for product precision performance estimation [20–22]. Furthermore, some innovative modeling methods, such as neural network [23], volumetric envelope [24], modal analysis [25,26] and graph theory [27,28] have also been applied to model the mating tolerance of both parts and assemblies. Some of the aforementioned models have contributed to the mainstream commercial computer-aided tolerancing software, such as 3DCS Variation Analyst, eM-TolMate, VisVSA, and CETOL 6 Sigma [29]. However, most of these tolerancing models excel at estimating the precision of a mechanical product at one state, while they cannot adapt to simulate a sequence of degradation trajectories. Moreover, the simplification of form deviation into a series of tiny dimensional, rotational, and translational deviations of geometric features in some of these models brings about inaccurate virtual representation in the modeling stage and unneglectable simulation errors in the synthesizing and optimizing stages.

As a response, the Skin Model Shape theory was introduced as a new paradigm for the modeling of product geometry considering shape variability. The theory of Skin Model Shapes is an integration of international geometric standards and the concept of the Skin Model [30–33]. Within the theoretical framework of GeoSpelling [34], the Skin Model Shape is commonly formed by discrete geometry schemes such as point clouds and surface meshes. It provides an approach for the employment of computational techniques on the Skin Model. Research studies have shown that the Skin Model Shape has grea<sup>t</sup> e fficiency for the representation of product geometry considering geometrical deviations and form errors [29,35,36]. Besides, since deviations in the Skin Model Shape are generated mainly by duplicating a sample pattern from the non-ideal manufacturing process, real part measurement, shape defect simulation, and statistical shape analysis are employed on Skin Model Shapes for convenience. For example, finite element analysis is introduced to help reduce local modeling drawbacks and improve model quality [37]. In addition, the discrete form of the Skin Model Shape has been proven to have grea<sup>t</sup> potential in further assembly simulation, performance evaluation, and tolerance analysis [38–40].

Significant e fforts have been made to improve the accuracy of product models when processing features are taken into consideration. Processing features, including systematic geometric deviations and microscopic form patterns, are essentially a ffected by and have been employed as an estimation for the performance of machine tools. The pattern of machining defects is detected in order to provide a comparison between the ideal model and the actual machining result, which is usually obtained by a coordinate measuring machine [41]. Most of the previous research has identified the machining state as a Boolean attribute: precision reliable or precision failure. If the measured degradation reaches a threshold, precision failure occurs, and the product fault time is defined. However—mainly due to the wear of the driving system, transmission gear, and cutting tool—a stochastic process, rather than binary classification, is more suitable for simulating and predicting the deterioration of machining precision and degradation of processing features. Sun et al. [42] emphasized that tool wear is a dynamic process extending from sharp to worn and possibly to breakage, and that multi-state classification could provide a more timely and accurate estimation. Dai et al. [43] pointed out that the degradation paths would be di fferent between one product and another. They conducted a comparison of Gaussian and logarithmic distribution in order to establish the degradation model. Ozcelik and Bayramoglu [44] verified the e ffect of tool wear on machining surface roughness and established a statistical model for prediction under various cutting conditions. Shu et al. [45] made the assumption that the next state of machine tool wear only depends on its current wear state. So, a non-homogeneous continuous-time Markov process model was used for modeling the total experience over the target life. After that, a linear mixed-e ffects model and maximum likelihood estimation was implemented to assess the wear evolution and lifetime of the tool [46]. Distribution of the residual life up to the wear threshold and

estimation of the wear level have also been researched. Moghaddass and Zuo [47] modeled the gradual degradation of a mechanical device as a continuous-time degradation process with multiple discrete states. A condition-monitored device based on an unsupervised parameter estimation method was developed with only incomplete information observable. A homogeneous continuous-time, finite-state semi-Markov model was established based on the past history of components by Cannarile et al. [48]. It was verified to be of grea<sup>t</sup> help in improving the diagnostic performance of an empirical classification system involving the degradation of mechanical systems. These applications of stochastic processes in manufacturing modeling have shown grea<sup>t</sup> potential in simulating the tolerance of a large branch of mechanical products in mass production based on measuring data.

The intention of this paper is to propose a solution to assembly tolerance design problems considering processing feature degradation. These problems are commonly encountered in the practical production of equipment manufacturing in industries such as high-precision computer numerical control machining, aeronautics, and astronautics. This paper focuses on assemblies containing observable parts (mainly outsourced parts, such as standard components or ultra high-precision parts) and predictable parts (mainly self-made parts, such as lathe beds, columns, and pillars). These two kinds of parts are treated in di fferent ways—sampling and measuring are usually implemented on the observable parts for precise and detailed modeling, while the predictable parts are modeled through mathematical simulations of the geometric tolerance, form errors, and assembly clearance. Some of the tolerance indices are given through measurements, while the to-be-designed tolerance indices are assumed and experimented repeatedly. Feature precision degradation is taken into consideration in tolerance modeling to enhance the reliability and robustness of a designed tolerance scheme. As a solution, a continuous-time variational multi-dimensional Markov process is introduced for modeling the degradation process, and the degraded surfaces are used in assembly simulations. On that basis, target assembly tolerance indices are synthesized through a series of numerical experiments. As a goal, a reduction in manufacturing cost and a guarantee in assembling probability is expected through the integration of these methods to the existing tolerance design framework, especially for products in mass production.

To achieve this, a tolerance design method considering processing feature degradation is proposed. The Skin Model Shapes form is adapted as a geometric basis. The nodal sampling point set on the machined surface is obtained by means of a high-precision coordinate measuring machine. A continuous-time multi-dimensional Markov process is trained to calculate the nodal displacement joint probability distribution on the machined surfaces. The degraded machined surfaces are predicted and applied in subsequent numerical experiments. To provide the precise assembly tolerance indices, assembly simulations and data analysis are conducted, and the assembly force constraints and assembly entity constraints are taken into consideration. Then, these tolerance indices are synthesized to provide tolerance schemes with low manufacturing costs and high assembly probabilities. Estimation and verification of the e fficiency of the proposed method are illustrated through an example of the transmission shaft on a five-axis high-precision machine tool (VTM200F).

#### **2. Predictive Machined Surface Modeling**

In this paper, machined surfaces are regarded as an integration of processing features and random geometric deviation. To a specific machine surface, the processing feature is determined by the basic geometric shape (plane, cylinder, sphere, etc.), the application and sequence of processing techniques (turning, milling, grinding, etc.), and the corresponding operating parameters (cutting speed, feed rate, cutting depth, etc.). Random geometric deviation depicts unpredictable and inevitable machining errors. Previous research has shown the rationality behind the idea that a stochastic process could represent the precision deterioration of a machining system. In the condition of mass production, the processing features are degrading due to the deterioration of the machining tool. Therefore, a multi-dimensional Markov model based on data training is introduced for modeling and prediction. The pattern of degradation becomes more explicit if the machined parts are from di fferent branches.

#### *2.1. Modeling of the Multi-Dimensional Markov Process*

The multi-dimensional Markov model is presented in this section, and a machined surface is used as an example. A uniform grid with a sampling interval Δ is placed on the nominal machined surface. The measuring process is calibrated through the grid nodes. A high-precision coordinate measuring machine is used to obtain the actual coordinates of sampling points on the real surface. The coordinate set is called the nodal sampling point set Ω. To a plane, the sampling interval Δ is a set distance along the x-direction and y-direction in the local machining coordinate system. To a spherical or cylindrical surface, the sampling interval Δ is a set angle around the central point or central axis.

First, we use *P*(*X<sup>t</sup> k* = *Lt k* ), *i* ∈ Ω to describe the probability that the actual coordinate of node *k* on the real surface is *Lt k* at time *t*. *Lt k* represents a discrete precision state by an interval of coordinate variation, as shown in Figure 1c. In a maintenance cycle of the machine tool, the shape error, dimension error, and location error of geometric features on real surfaces deteriorate over time. Therefore, the actual precision state of nodal points gradually moves away from the nominal positions:

$$P(X\_k^0 = X\_{k\\_norm}) > P(X\_k^1 = X\_{k\\_norm}) > \cdots > P(X\_k^n = X\_{k\\_norm}) > \cdots \tag{1}$$

Obviously, the actual precision state *Lt*+<sup>1</sup> *k* of node *k* at time *t* + 1 is directly related to *Lt k*, and has no relation to the precision states at and before time *t* − 1. Treated as a sampling on the time-axis of a continuous-time stochastic model of the machine surface, the coordinate change of a single nodal point over a period of time conforms to the discrete-time Markov property:

$$\begin{array}{l}P(\mathbf{X}\_{k}^{t}=L\_{k}^{t}) = P(\mathbf{X}\_{k}^{t}=L\_{k}^{t}|\mathbf{X}\_{k}^{t-1}=L\_{k}^{t-1}) \times P(\mathbf{X}\_{k}^{t-1}=L\_{k}^{t-1}|\mathbf{X}\_{k}^{t-2}=L\_{k}^{t-2}) \times \cdots\\\times P(\mathbf{X}\_{k}^{1}=L\_{k}^{1}|\mathbf{X}\_{k}^{0}=L\_{k}^{0}) \times P(\mathbf{X}\_{k}^{0}=L\_{k}^{0})\end{array} \tag{2}$$

In the machining process, a continuous cutting motion on the machining surface causes the coordinate change of relevant nodal points to conform to the Markov property. With the same form and parameters of the driving system and gear motion, the combination of motion between adjacent nodal points is fixed and constant. These nodal points are called relevant points. Inside the machining system, the same motion combination is conducted between pairs of relevant points. *Lt k* relates to the coordinates of all the neighboring points , and has no relation to those points {*k*1, *k*2, ··· , *km*}, and has no relation to those points outside the neighborhood, as shown in Equation (3):

$$P(X\_k^t = L\_k^t) = P(X\_k^t = L\_k^t | X\_{k\_1}^t = L\_{k\_1}^t, X\_{k\_2}^t = L\_{k\_2}^t, \dots, X\_{k\_m}^t = L\_{k\_m}^t), \ k\_1, k\_2, \dots, k\_m \in \mathbb{S}\_k \tag{3}$$

Combining Equations (1) and (2), the actual precision state of nodal points on the machined surface conform to the temporal and spatial Markov properties at the same time, as shown in Figure 1a.

**Figure 1.** Temporal and spatial Markov properties of planar machining surface (**a**), and an illustration of the precision state degradation in the production sequence (**b**,**<sup>c</sup>**).

#### *2.2. Calculation of Model Weight Parameters*

The degradation of the feature precision state of parts from the same batch {*T*0, *T*1, ··· , *Tn*} is a continuous-time process if the parts are manufactured on the same machine tool. Condition characteristics *m* and transmission characteristics *n* are unique properties on or between relevant points. Therefore, the weight parameters set {*w*} = {*<sup>w</sup>*1, *w*2, ··· , *wn*} is applied to the nodal sampling point set.

We use *x* to describe the condition of the rough surface before machining, and *y* to describe the machined surface at time *t*. As a description of feature precision degradation, the continuous-time nodal displacement joint probability distribution *<sup>P</sup>*(*y*|*<sup>x</sup>*, *t*) in the form of the temporal Markov process is formulated in Equation (4a), while *<sup>P</sup>*(*y*|*<sup>x</sup>*, *t*) in the form of the spatial Markov process is formulated in Equation (4b), as a description of feature precision specification at time *t*:

$$\begin{cases} P(y(t)|y(t-1), \mathbf{x}, t) = \frac{\exp\{\sum\_{i,j} \lambda\_i m\_i(t\_j, y, \mathbf{x}, i) + \sum\_{i,j} \mu\_{i,j} m\_j(t\_{j+1}, t\_j, y, \mathbf{x}, l)\}}{\sum\_{i,j} \exp\{\sum\_{i,j} \lambda\_i m\_i(t\_j, y, \mathbf{x}, i) + \sum\_{i,j} \mu\_{i,j} m\_j(t\_{j+1}, t\_j, y, \mathbf{x}, l)\}} (\mathbf{a}) \\\ P(y|\mathbf{x}, t) = \frac{\exp\{\sum\_{i,j} \lambda\_i m\_i(t\_j, y, \mathbf{i}) + \sum\_{h,j} \eta\_{h,j} m\_{h,j}(\mathbf{x}\_{h+1}, \mathbf{x}\_h, y, \mathbf{b})\}}{\sum\_{i,j} \exp\{\sum\_{i,k} \lambda\_i m\_k(\mathbf{x}\_{i-1}, \mathbf{x}\_i, y(t), \mathbf{i}) + \sum\_{i,j} \eta\_{h} n\_j(\mathbf{x}\_i, y, \mathbf{j})\}} \end{cases} (\mathbf{a})$$

where *i* is the number of nodal points, *j* represents the time sequence, and *l* and *h* represent the temporal and spatial relevant points. λ represents the weight parameter of *m*, and μ represents the weight parameter of *n*. The unified form of Equation (4) is shown in Equation (5):

$$P(y|\mathbf{x},t) = \frac{\exp\{\sum\_{k=1}^{K} w\_k f\_k(\mathbf{x}, y, t)\}}{\sum\_{\mathbf{x}, \mathbf{y}} \langle \exp\sum\_{k=1}^{K} w\_k f\_k(\mathbf{x}, y, t)\rangle} \tag{5}$$

Through the measured data from coordinate machining machine and the model shown in Equation (5), maximum likelihood estimation is applied to calculate the weight parameters. The value of *wk* indicates that the observed point set has the highest probability. We use a one-order step function as the object function. The overfitting problems in the training process caused by sampling data loss or over-training can possibly cause the verification error to grow larger as the calculation complexity grows and the training error decreases. To overcome this, an *L*2 regularization term weighted as 12δ2 is added to the object function, as shown in Equation (6):

$$g(w) = \frac{\partial \ln(p(y|\mathbf{x}, t))}{\partial w\_k} - \frac{\sum w^2}{2\delta^2}. \tag{6}$$

As a solution to Equation (6), a quasi-Newton iteration is implemented to calculate the parameter weights, which makes *g*(*w*) converge to zero:


$$\mathcal{M}^{(k+1)} = \mathcal{M}^{(k)} + \frac{y\_k y\_k^T}{y\_k^T \delta\_k} - \frac{\mathcal{M}^{(k)} \delta\_k \delta\_k^T \mathcal{M}^{(k)}}{\delta\_k^T \mathcal{M}^{(k)} \delta\_k} \tag{7}$$

(4) Use a line search function *h*(λ) = *f*(*w*(*k*) + <sup>λ</sup>*kpk*). Update *w*(*k*) when *h*(λ) reaches the global minimum.

(5) Output *w*(*k*) if *g*(*w*(*k*)) = 0. Otherwise, *k* = *k* + 1 and go back to step (2).

#### *2.3. Prediction of Degraded Surface*

When *<sup>P</sup>*(*y*|*<sup>x</sup>*, *t*) (as calculated in Section 2.2) reaches its maximum, the predictive machined surface is equal to *y*\*, which means that the greatest possible location distribution of sampling points at time *t* is obtained as the predictive surface:

$$y^\* = \underset{y}{\text{argmax}} \, P\_w(y|x, t) = \underset{y}{\text{argmax}} \sum\_{k=1}^K w\_k f\_k \tag{8}$$

Under a certain value of *<sup>P</sup>*(*y*|*<sup>x</sup>*,*<sup>t</sup>*), a backward path searching iteration is implemented after obtaining a specific sampling point's initial location probability <sup>δ</sup>1(*l*), so as to calculate the most probable location of the sampling point during the whole machining process:

$$\delta\_i(t) = \max\_{1 \le j \le m} \{ \delta\_i(t-1) + \sum\_{k=1}^K w\_k f\_k | y\_{i-1} = j, y\_i = t \}, t = 1, 2, \dots, N \tag{9}$$

The prediction of the feature precision state on the machined surface is achieved when the backward path searching reaches its end point, as shown in Equations (10) and (11):

$$\max\_{y} \sum\_{k=1}^{K} w\_k f\_k(y, \mathbf{x}) = \max\_{1 \le j \le m} \delta\_n(j), \tag{10}$$

$$y\_i^\* = \underset{1 \le j \le m}{\text{argmax}} \delta\_{\boldsymbol{\mu}}(j). \tag{11}$$

The calculated *<sup>P</sup>*(*y*|*<sup>x</sup>*,*<sup>t</sup>*) contains the information of crafts and steps. At time *t*, the predictive machined surface *y*\* (in Equation (12)) consists of the predicted coordinates of a nodal sampling point set under processing features. If time *t* = *tf* is large enough to cover the designated lifecycle of the to-be-designed assembly, the degradation of the part feature is considered by using the predictive machined surface *yf*\*:

$$y\_f^\* = \langle y\_{f,1'}^\* y\_{f,2'}^\* \cdots, y\_{f,N}^\* \rangle \tag{12}$$

#### **3. Constrained 3D Assembly Simulation and Tolerance Synthesis**

#### *3.1. Constrained Assembly Simulation*

Then, the tolerance models generated by the methods described in Section 2 are used in assembly simulations. Generally, the process of assembling two parts involves the transformation of one part from its initial condition to its assembling condition, while the other part is relatively static. In this paper, the former part is called assembly part *Q*, and the latter part is called reference part *P*. The assembling transformation ξ is in the form of a twist because of its splendid performance in representing the orientation, and its independence from the coordinate system transformation. Parts are considered under the geometric basis of Skin Model Shapes in this paper, which means that *Q* and *P* are in the form of meshes or point clouds. So, the assembly part *Q* is transformed through ξ from its initial condition *Q* to assembling condition *Q*, as shown in Equation (13):

$$Q\_{\prime} = \xi \cdot Q \tag{13}$$

Twist ξ is essentially an affine transformation. It can always be decomposed to a motion screw ξ*i* = (*si*, *vi*) in a three-dimensional (3D) velocity vector field, containing a rotation around its axis and a translation along the axis, as shown in Equation (14):

$$
\xi\_i = s\_i + p\_i \ltimes v\_i \tag{14}
$$

where *pi* is the vector from the spatial origin to point *i*. *vi* is the axial velocity vector at point *i*. *si* is the tangential velocity vector at point *i*. Assembly twists of all the points of the assembly part consist the assembly transformation: ξ = {ξ1, ξ2, ··· , ξ*<sup>i</sup>*, ···}.

Mathematically, *Q* is the global extremum of an object function containing *Q*, *P*, and twist ξ under assembly constraints, as shown in Equation (15):

$$
\xi = \operatorname\*{argmin}(f(P, Q')) \tag{15}
$$

Choose the root mean square Euclidean distance function in Equation (16) of the closest point pairs in point clouds *P* and *Q*:

$$f = RMSD(P, Q\nu) = \sqrt{\frac{1}{n} \sum\_{p\_i, q\_i \in CPP(i)} \|p\_i - q\_i\|^2} \tag{16}$$

Establish the iteration object function based on Equations (13)–(16). The end condition of the iteration is the appearance of three pairs of contact points in the assembly. A pair of contact points contains a point in *P* and a point in *Q*. Furthermore, the Euclidean distance between the two points is smaller than a threshold value. Usually, using the iterative closest points (ICP) is a common method under the circumstances to find the target. However, assembly entity constraints and assembly force constraints have restricted the direct use of ICP.

Assembly entity constraints restrict the intersection of part entity in assembly simulation. As shown in Figure 2a, the direct use of ICP causes an unreasonable intersection of part materials. However, if the rule of assembly entity constraints is considered, ICP would give a result that is much more rational (in Figure 2b). A signed distance field is introduced to supervise the compliance of the assembly entity constraints. The positive direction of the signed distance field is always against the assembling direction. The signed distance vector *dsdv i* (*i* = 1, 2, ··· , *n*) is from the points in *P* to the corresponding closest point in *Q*. When the angle between *dsdv i* and the positive direction is under 90 degrees, *dsdv i* is labeled positive. Otherwise, it is labeled negative. As shown in Figure 2c, all the signed distance vectors *dsdv i*(*i* = 1, 2, ··· , *n*) are positive during iteration before assembling. However, negative signed

distance vectors appear when the rule of assembly entity constraints is violated (in Figure 2d). So, the signals of *dsdv i*(*i* = 1, 2, ··· , *n*) are repeatedly checked during the iteration.

**Figure 2.** Illustration of assembly entity constraints and the application of a signed distance field in realizing the rule.

Assembly force constraints ensure the stability of assembling conditions through restricting the location of contact points according to the assembly force. Obviously, three points provide a stable support in 3D assembly simulation, so the three contact points should not be in a line. Additionally, the assembly force vector should cross the spatial triangle made up by the three contact points. Otherwise, the current condition is unstable, and additional rotations would happen. The compliance of assembly force constraints is only checked when the iteration ends, ensuring that the final assembling state is stable. The solution of the assembling condition is decomposed into four main steps:


**Figure 3.** Main steps in the computation of assembling condition when considering assembly force constraints.

#### *3.2. Static and Dynamic Tolerance Synthesis*

In this paper, tolerance synthesis assigns tolerance indices *X*, minimizes manufacturing cost, maximizes assembling probability, and meets the assembly functional requirement *Y* at the same time. The *X* tolerance index comprises the given indices *Xg* (mainly in outsourced parts) and to-be-designed indices *Xd* (mainly in self-made parts). *Xg* and *Xd* are independent from the perspective of manufacturing, but they become relevant because of the assembly functional requirement *Y*. Additional objectives including the manufacturing cost and assembling probability turn this issue into an optimization problem.

To increase the robustness of the design scheme, feature degradation during the machining process is considered when generating observable parts. Predictions of part surfaces at any time during the designed life are made using the methods in Section 2. Generally, surfaces in their late design age with seriously deteriorated geometric features are used in tolerance synthesis. The surfaces of predictable parts are simulated by a multivariate Gaussian process (MGP) based on given tolerance indices. Repetitive experiments and variated models are applied in the design process due to the stochastic characteristic of the MGP. Numerical experiments of assembly simulation provide assembly tolerance indices, which are usually categorized into static tolerance indices (such as verticality and parallelism) and dynamic tolerance indices (such as circular runout and end face runout). Static tolerance indices are certain indices based on measurements and the analysis of assembled Skin Model Shapes. These only vary when the tolerance models are varied (in Figure 4a). However, dynamic tolerance indices are obtained through monitoring a series of static tolerance indices by analyzing the possible motion and translational/rotational samplings. Measuring positions are sampled for experiments through a Monte Carlo process. As shown in Figure 4b, the upper end face parallelism of a sliding table on a high-precision guide is a typical dynamic tolerance index, which is calculated

through a series of assembly simulations after sampling measuring positions in its translational motion. Similarly, common runout errors of cylindrical parts (such as axle parts) can also be obtained by sampling measuring orientations in the rotational motion of those parts. A small sampling interval and large sampling quantity lead to a precise estimation.

**Figure 4.** Synthesis of static and dynamic assembly tolerance indices based on assembly simulation.

As shown in Figure 4, assembly tolerance indices are estimated through data analysis in each measuring position/orientation. The distribution of estimated assembly tolerance indices to the to-be-designed tolerance indices is summarized in the form of tables or figures. Then, an integration of assembly functional requirements and manufacturing cost estimation is introduced to outline the objectives. Due to the stochastic characteristic of tolerance models, not all the examples with the same assigned tolerance scheme can meet the assembly functional requirements. Generally, if the ratio of conformation reaches a certain threshold (such as 95%), the tolerance scheme is treated as reliable. On that basis, the manufacturing costs of all the reliable tolerance schemes are computed and compared based on industrial statistics and empirical formulation. That with the lowest costs is selected as the designed tolerance scheme.

## **4. Case Study**

#### *4.1. Description of the Tolerance Allocation Problem*

The proposed method is applied to solve an assembly tolerance design problem in the rotary feed system of a VTM200F5 five-axis turning and milling complex center. General tolerance indices of the rotary feed system have a significant influence on the cutting precision and are extremely sensitive to geometric deviation in the cutting process due to its direct link to the cutting tool. Considering assembly tolerance, the rotary feed system is simplified to two key parts: the hole part describing devices fixing rotary bearings, and a shaft part depicting the motor spindle component. The shaft part is an outsourced part made by Ti–Al thermostable alloy with its surface specially heat-treated. So, the key accuracy indices such as location accuracy, shape accuracy, and surface accuracy should not be modified in order to protecting the special surface coating. The hole part is a self-made part made by common bearing steel. Its key accuracy indices are to be designed, as shown in Figure 5. This is a typical component containing observable parts (with the same pattern of processing features and from di fferent batches) and predictable parts (to-be-designed parts).

**Figure 5.** An assembly tolerance design example in designing the rotary feed system of a five-axis turning and milling complex center (VTM200F5).

In this component, the runout accuracy of the distal end of the shaft part is set as the target accuracy index due to its direct influence on the cutting accuracy. In the shaft part, the nominal flatness of the distal end surface is 0.010 mm. The nominal total runout of the distal end surface is 0.020 mm. The nominal total runout of the proximal end surface is 0.010 mm. The nominal verticality between the distal end surface and axis C is 0.010 mm. The nominal flatness of the shoulder surface is 0.010 mm. In the hole part, surface A is the assembly positioning surface. The key accuracy indices to be designed in this tolerance allocation problem are listed in Table 1. The flatness of the shoulder surface is *FLF*. The total runout of inner hole is *TRD*. The verticality between axis D and surface A is *VDA*. The parallelism of the shoulder surface and surface A is *PLEA*.


**Table 1.** List of the target accuracy term of the assembly, the nominal value of the accuracy terms of the shaft part, and the to-be-designed accuracy terms of the hole part.

#### *4.2. Stochastic Process Training and Parameter Calculation*

Twenty outsourced parts were randomly chosen. Parallel lines in two directions were drawn on assembly surfaces with a uniform interval of 0.01 mm. Sampling points are defined as the intersection points of these lines. We used a non-contact coordinate measuring machine with an accuracy of at least 0.001 mm to measure the actual location of sampling points. The measured coordinates amounted to 20-point cloud sets {*T*0, *T*1, ··· , *T*20}. Take surface F as an example: the measured point clouds sets are {*TF*0, *TF*1, ··· , *TF*20}. Two of the three coordinates are determined by sampling points. The third coordinate is called the feature parameter. A first-order gradient object function of maximum likelihood estimation with weight decay is shown in Equation (17):

$$\log(w) = \sum\_{i=1}^{N} \sum\_{t=1}^{T} f\_k(\mathbf{x}\_t, y\_t, t) - \sum\_{i=1}^{N} \sum\_{y} \left[ p(y|\mathbf{x}, t) \times \sum\_{t=1}^{T} f\_k(\mathbf{x}\_t, y\_t, t) \right] - \frac{\sum w^2}{2\delta^2} \tag{17}$$

Set the initial parameter weights as *w*(0) = {0.1, 0.1, ··· , 0.1} and the initial positive definite matrix as *M*(0) = *diag*(0.1, 0.2, ··· , 2.0). The coefficient of the correct term is set to 0.5. The maximum number of iterations is 1000. The feature weights of *<sup>P</sup>*(*y*|*<sup>x</sup>*, *t*) corresponding to processing features are calculated according to Section 2.1. Some of the feature weights of surface F are shown in Table 2.


**Table 2.** Value of partial feature weights of surface F.

The prediction of degraded machined surface F is obtained using the methods in Section 2.3. Consider a particular sampling point (*i*, *j*), with the assumption that its feature parameter of the location on the predictive machined surface is *y*∗*i*,*j* and its feature parameter of the location on the real surface is measured as *y*(1) *i*,*j* . The geometric average of the overall location residuals *Er* is set as an estimation of the precision of the predictive machined surface model, as shown in Equation (18):

$$E\_r = \frac{\sum\_{N=1}^{20} \sqrt{\sum\_{i,j} \left( y\_{i,j}^{(N)} - \overline{y}\_{i,j}^{(N)} \right)^2}}{\sum\_{N=1}^{20} \sqrt{\sum\_{i,j} \left( y\_{i,j}^\* - y\_{i,j}^{(N)} \right)^2}} \tag{18}$$

The value of *Er* decreases when the difference between the predictive location and actual location of all the sampling points grows larger. Otherwise, the increase of *Er* indicates that the predictive machined surface model is in good agreemen<sup>t</sup> with the real surface. Considering all the values of surface F in the 20-part samples, a comparison of the predictive machined surface using a multivariate Gaussian process and using the methodology in this paper is shown in Figure 6. Compared with the commonly used multivariate Gaussian process, the multi-dimensional Markov process provides a predictive model with higher accuracy. When time *t* is applied to the designed life of the assembly, the degraded machined surface near the precision failure threshold is generated and used in the assembly simulation.

**Figure 6.** Comparison of overall location residuals of the predictive machined surface using the multivariate Gaussian process and using the method proposed in this paper.

#### *4.3. Tolerance Synthesis of Example Rotary Feed Component*

Numerical experiments were implemented with tolerance models of the shaft part and the hole part (in Figure 7a). Machined surface degradation prediction was combined as the model of the outsourced shaft part. Testing models for the hole part were generated based on MGP with given tolerance indices *e*1, *e*2, and *e*3. Using the method in Section 3.1, assembling conditions and contact points were estimated based on ICP when the assembly entity constraints and assembly force constraints were checked during iterations. It is worth noting that the location of contact points varies according to the assembling contact. There are mainly three basic patterns from the perspective of the

position distribution of contact points (illustrated in Figure 7a–c), which is the probability influenced by the scale of planar form errors, cylindrical form errors, and perpendicularity.

**Figure 7.** Numerical experiments (**a**) and three basic patterns of contact point positions (**b**–**d**).

The target assembly tolerance index, the runout error *ef* of surface F, was estimated by analyzing the transformed point cloud in assembling condition. A four-dimensional solution space was generated as an explanation of tolerance synthesis, constructed with *e*1, *e*2, and *e*3 as input and *ef* as an output. According to the functional requirement of the assembly, the nominal value of *ef* is 0.020 mm, which means the upper limit of *ef* is 0.020 mm in all the adopted tolerance schemes. However, an unnecessarily narrow tolerance interval arrangemen<sup>t</sup> leads to extremely high costs in fabrication and maintenance, while a combination of wide tolerance intervals causes a high ratio of assembling failure. As a result, the influence of manufacturing cost and assembly probability based on an empirical formula and experimental statistics is considered in tolerance synthesis. Assembly examples with acquired functional requirements are called "successful" assemblies. Assembling probability is estimated by the proportion of successful assemblies of all the experimental assemblies. According to industrial practice, a yield rate of 95% is accepted, which means that the assembling probability should be more than 95%. The region of interest in the solution space is generated following this rule.

The empirical relationship between the tolerance interval and manufacturing cost based on historical industrial practice is shown in Equation (19), where α1, α2, and α3 are manufacturing difficulty coefficients, and K is a scaling factor:

$$C = K(\frac{\alpha\_1}{e\_1^2} + \frac{\alpha\_2}{e\_2^2} + \frac{\alpha\_3}{e\_3^2}) \tag{19}$$

The influences of *e*1, *e*2, and *e*3 on the assembly probability and manufacturing cost are shown in Figure 8a. A Monte Carlo sampling method was applied in the four-dimensional solution space to calculate the output of all the sampled inputs. High-density sampling was conducted in the region of interest to ensure searching accuracy; the density of sampling points was lower, but even outside the region of interest, it avoided falling into the local optimum. The top 1% with the lowest manufacturing cost was approximately fitted into an ellipsoid limit surface, as shown in Figure 8b, where the axes of *e*1, *e*2, and *e*3 are even, and the axis of *ef* is uneven.

**Figure 8.** Distribution of assembly probability and manufacturing cost (**a**), and solution searching in a four-dimensional solution space constructed by even axes of *e*1, *e*2, and *e*3, and the uneven axis of *ef* (**b**).

Three examples with the lower manufacturing cost and satisfactory assembling probability constitute the designed tolerance scheme, as shown in Table 3. The designed scheme with an assembling probability larger than 95% and the lowest manufacturing cost is listed in bold. The designed value of all three to-be-designed tolerance terms *e*1, *e*2, and *e*3 were enlarged by 0.003, 0.004, and 0.004, respectively. The objective functional requirement was still met, and the assembling probability was satisfactory. However, the expansion of the tolerance interval caused a 35.6% reduction of the manufacturing cost. The other two listed example designed schemes can also be considered depending on the acceptable manufacturing cost considered.


**Table 3.** Design requirement of the object tolerance term and the nominal and designed value of the to-be-designed tolerance terms.
