Manufacturing variation models in multi-station machining systems

In product design and quality improvement fields, the development of reliable 3D machining variation models for multi-station machining processes is a key issue to estimate the resulting geometrical and dimensional quality of manufactured parts, generate robust process plans, eliminate downstream manufacturing problems, and reduce ramp-up times. In the literature, two main 3D machining variation models have been studied: the stream of variation model, oriented to product quality improvement (fault diagnosis, process planning evaluation and selection, etc.), and the model of the manufactured part, oriented to product and manufacturing design activities (manufacturing and product tolerance analysis and synthesis). This paper reviews the fundamentals of each model and describes step by step how to derive them using a simple case study. The paper analyzes both models and compares their main characteristics and applications. A discussion about the drawbacks and limitations of each model and some potential research lines in this field are also presented.


Introduction
Traditionally, product design has been separated from manufacturing process design along the product development cycle increasing ramp-up times, product changes costs, and variability of product quality. This product-oriented approach, named as over-the-wall design due to the sequential nature of the design activities, prevents the integration of design and manufacturing activities to improve product development [1]. In order to overcome this limitation, manufacturers have begun to investigate the ways to simultaneously evaluate product designs and manufacturing processes in an attempt to eliminate downstream manufacturing problems and reduce ramp-up times. For this purpose, product design requires the application of processoriented approaches through 3D manufacturing variation models to integrate product and manufacturing process information [2]. However, the application of 3D manufacturing variation models is nowadays limited, specially in multi-station machining processes (MMPs) where a large number of machining operations are conducted in different stations with different fixture devices. In fact, in product design and quality improvement fields, the development of reliable 3D variation models of MMPs is a key issue to estimate the resulting geometrical and dimensional quality of manufactured parts.
In the literature, one can distinguish two important group of researchers dealing with the development of 3D variation models for MMPs: (a) a first group of researchers more focused on the quality improvement field, mostly universities from USA, and (b) a second group of researchers focused on product design, mostly universities from France and Canada. The research conducted by the first group that we name as Industrial Operations school ("IO school") is focused on modeling the dimensional and geometrical variation of manufactured parts by the state-space model which is commonly applied in control theory [3]. Through this model, a large number of quality improvement activities can be conducted such as part quality control, optimal placement of inspection stations, robust process planning, process-oriented tolerancing, etc. The research conducted by the second group that we name as Product Design school ("PD school") is focused on modeling dimensional and geometrical part quality variations considering the system workpiece/fixture/machine tool as a mechanical assembly. By this consideration, well-known approaches for the analysis of mechanical assemblies can be applied. Specifically, the PD school applies the concept of small displacement torsors (SDTs) to describe and propagate the surface deviations along the different machining stations. The application of the PD school investigations are mainly focused on tolerance analysis and synthesis.
Both schools have been very active during last decade, and their 3D manufacturing variation models have been used in a large number of applications. Surprisingly enough, these schools have been working in parallel, without using potential synergies or adapting the advances from one school to the other. Moreover, despite the success of both schools and the large number of publications related to them, there is no review in the literature that analyzes and compares both 3D manufacturing variation models. This paper contributes to fill this literature gap, providing a critical comparison of both models, paying special attention to (a) fixtures and processes supported, (b) limitations of virtual inspection, (c) GD&T conformance, (d) applicability, (e) modeling accuracy, and (f) modeling complexity. Furthermore, the review also exposes future lines of research that should be addressed by both schools. The paper is organized as follows: Section 2 describes the fundamentals of the 3D manufacturing variation model applied by the IO school and its main applications. Similarly, Section 3 describes the 3D manufacturing variation model applied by the PD school together with its main applications. Section 4 presents two modeling 2D examples (extensible to any 3D case study) to show, step by step, the implementation of both models, analyzing them symbolically and numerically. Section 5 discusses both 3D variation models, comparing their main drawbacks and advantages. Finally, Section 6 concludes the paper and outlines some potential research lines in the field.

Fundamentals
Manufacturing variability in MMPs and its impact on part quality can be modeled by capturing the mathematical relationship between the sources of variation of manufacturing process variables that are critical to part quality (named key control characteristics (KCCs) of the process) and the deviations of the functional features of the product (named key product characteristics (KPCs)). This relationship is modeled through a non-linear function where y 1 is a KPC and u 1 , u 2 , . . . , u n are the KCCs in the MMP. By the assumption of small variations, the non-linear model can be linearized through a Taylor series expansion, and the part quality variation can be defined as where u = [u 1 , u 2 , . . . , u n ] T , ε 1 contains the high-order non-linear residuals of the linearization, and the linearization point is defined byū = [ū 1 ,ū 2 , . . . ,ū n ] T . This linear approximation can be considered good enough for many MMPs [4]. Considering that there are M KPCs in the part which are stacked in the vector Y = [ y 1 , . . . , y M ] T , Eq. 1 can be re-written in a matrix form as where U = [ u 1 , . . . , u n ] T and u j = u j −ū j for j = 1, . . . , n defines the small variations of the KCCs in a MMP; is the matrix ; and ε is the stacked vector of the high-order non-linear residuals. In MMPs, the derivation of Eq. 2 is a challenging task. Researchers from USA universities have proposed the adoption of the well-known state-space model from control theory [3] to represent mathematically the relationship between the sources of variation of a MMP and the deviation of the machined surfaces at each station, including how the deviation of previous machined surfaces influences at current station when these surfaces are used as locating datums. In this representation, dimensional deviations of part surfaces from nominal values are defined by 6×1 vectors named state vectors, in the form of T is the small orientation deviation of the local CS of the ith part surface. The notation • i refers to the nominal CS and i to the current CS. The deviation of all part surfaces at the station k are stacked in the state vector Note that each feature deviation, x k,i , is expressed in its own CS.
In a machining station, three main sources of variation can be distinguished: datum-induced deviations, fixture-induced deviations, and machining-induced deviations. The state-space model defines analytically how these three main sources of error influence on the final part quality deviation. To illustrate how these three main sources of variation influence on part quality, we consider an N-station machining process shown in Fig. 1 and the kth machining station with the workpiece and the fixture device shown in Fig. 2. At this kth station, the sources of variation below exist. First, the deviations of the datum surfaces used for locating the workpiece deviate the workpiece location from its nominal value. This term can be estimated as is the vector of part surface deviations from upstream machining stations and A k linearly relates the datum deviations with the machined surface deviation due to the locating deviation of the workpiece.
Secondly, the fixture-induced deviations deviate the workpiece location on the machine tool table and produce a machined surface deviation. This term can be estimated as k is the vector of locator deviations and B f k is a matrix that linearly relates locator deviations with the machined surface deviation.
Thirdly, the operation or machining deviations such as those due to geometrical and kinematic errors, toolwear errors, etc. deviate the cutting tool tip during machining, and thus, the machined surface is deviated from its nominal value. This term is modeled as x m k = B m k · u m k , where u m k is the vector that defines the KCCs related to operation or machining deviations and B m k is a matrix that linearly relates these KCCs with the machined surface deviations. Therefore, for an N-station machining process, the derivation of the state-space model can be defined in a generic form as where B k · u k represents the deviations introduced within station k due to the KCCs (related to fixturing and machining) and it is defined as T and w k is the unmodeled system noise and linearization errors.
The derivation of the state-space model in MMPs, named the stream of variation (SoV) model, was firstly presented by Huang et al. [5]. Djurjanovic et al. [6] expanded Huang's work in order to explicitly derive the linear equations that model the relationships between fixtures, locating datum, and measurement datum features. In their research work, a complex mathematical derivation was required, and the methodology proposed was not straightforward to be applied. The SoV model derivation was improved by [4] who applied the differential motion vector (DMV) concept from robotics to represent the geometric deviation of each machined feature. In their work, a step-by-step methodology was proposed in order to derive the matrices A k and B k at each station using product and process information (i.e., part geometry and fixture layouts). Previous works were limited to 3-2-1 orthogonal fixture layouts based on locators and generic cutting tool path deviations without explicitly including specific machining-induced errors. Loose et al. [7] extended the state-space model formulation by including general non-orthogonal fixture layouts based on locators. More recently, Abellan-Nebot et al. [8] expanded the formulation of matrix B k in order to include common machining sources of variation such as those due to tool wear, thermal expansions, cutting tool deflections, and geometric-kinematic machine tool errors.

Virtual part quality inspection and verification
Somewhere along the MMP, an inspection station can be placed in order to inspect the KPCs and verify whether the workpiece/part is within specifications. Following the state-space model formulation from control theory [3], a virtual inspection after the kth machining station can be conducted using the expression: where y k represents the deviations of the inspected KPCs, C k · x k are the deviations of the KPCs that are defined as a linear combination of the deviations of workpiece features at the kth station, and v k is the measurement noise of the inspection process.
In a similar way to x k , vector y k is defined as [y T k,1 , . . . , y T k,q , . . . , y T k,M ] T , where y k,q is the inspected deviation of the qth KPC (denoted as S q ) defined by the vector y k,q = [(d S m S q ) T , (θ S m S q ) T ] T , where S m is the measurement datum surface and M is the number of KPCs inspected. In Eq. 4, matrix C k depends on what KPCs they are and which measurement datums are used to locate the part in the inspection station. How to derive matrix C k is explained in detail in [4].
In order to express the part quality measurements by an explicit linear function of the KCCs presented along the MMP and considering that the inspection station is placed at the end of the MMP (after machining station N), Eqs. 3 and 4 can be combined and rewritten in the input-output form as: where the vectors Y and U are the stacking quality vectors after inspection and the vectors of sources of error, respectively, from the stations k = 1, 2, . . . , N.
. . , u T N ] T , and matrices and ε are defined as: where Using Eq. 5, the deviation of the local CS of an inspected KPC S q w.r.t. the measurement datum S m can be estimated. However, if one wants to verify if an inspected KPC is within its tolerance zone according to geometric dimensioning and tolerancing (GD&T) practices, the deviation of the boundary points of the inspected KPC w.r.t. the measurement datum should be evaluated. For this purpose, the deviation of the rth boundary point P r of the inspected qth KPC w.r.t. S m can be evaluated according to the expression: where y N,S q is the deviation of the qth KPC obtained from Eq. 5, I 3×3 is the 3 × 3 identity matrix, 0 3×3 is the 3 × 3 zero matrix, andt S q P r is the skew matrix of the nominal position vector t S q P r which describes the position of the point P r w.r.t. the S q . The resulting deviation of point P r w.r.t. S m is then defined by a 6 × 1 deviation vector in the form of y N,P r = [(d S m P r ) T , (θ S m P r ) T ] T . The deviation of the rth point of the toleranced surface following the direction of part verification, defined by the vector n = [n x , n y , n z ] T , is evaluated by the expression: The tolerance zone where the variability of the manufactured feature will lie can be obtained by analyzing the deviation from nominal values of all boundary points of the KPC. If the tolerance of a KPC is defined from design specifications, then one can be interested in verifying the part according to a given MMP. For this purpose, the distance between each deviated boundary point and the specified tolerance zone from design should be evaluated. For the point P r , this distance is defined as the gap distance Gap r , and it is formulated as: where τ is the maximum deviation of the rth point according to the tolerance size (e.g., for a positional tolerance, τ is t/2, where t is the tolerance value). The rth point of the inspected surface will be within the tolerance zone if Gap r remains positive or null (see Fig. 3). Analyzing the deviation of all boundary points of the KPC, the verification of the GD&T tolerances applied to the KPC can be conducted. It should be remarked that the virtual inspection and verification can be conducted following two main approaches: the worst-case approach and the statistical approach. Depending on which type of approach is applied, the estimation of the KPC deviation and, thus, the estimation of the deviation of its boundary points in order to analyze a functional specification will be more or less conservative. For each approach, the resulting estimation is derived as follows: -Worst-case approach: The estimated deviation of the qth KPC will be the maximum according to the expected sources of variation. According to Eq. 5, the worst-case analysis can be conducted assuming that all coefficients from matrix and vector U are positive and the measurement error also increases the expected deviation. The worst-case deviation of the S q CS is: and the worst possible part quality considering the point boundary deviation is defined as: Gap wc = min Gap r wc ∀r ∈ boundary points, where Gap r wc is evaluated by Eqs. 10-12 considering y N,S q−wc instead of y N,S q . -Statistical approach: The worst-case analysis produces an estimation that is highly improbable, specially for a large number of sources of variation due to the randomness of the sources of variation To estimate a more probable scenario, the statistical analysis is commonly applied. In this analysis, the sources of variation are assumed to be independent to each other and normally distributed. Under these assumptions, the covariance of the S q CS can be estimated as: where • is the covariance matrix of •. Therefore, the deviation of the KPC is estimated, according to 6σ , as: The part quality considering the point boundary deviations is defined as: where Gap r st is evaluated by Eqs. 10-12 considering y N,S q−st instead of y N,S q .

Part quality estimation and process planning
The straightforward application of the SoV model is part quality estimation (i.e., tolerance analysis) which leads the designer to estimate if the MMP is able to manufacture parts within specifications. By this analysis, the process planner can search the most robust MMP to manufacturing disturbances from a group of candidates and conduct specific modifications to improve the manufacturing process. Zhang et al. [9] presented a sensitivity analysis based on the SoV model to assess how sensitive the KPCs are to certain fixture-induced variations in an MMP. Through the sensitivity indices, the robustness of each process plan candidate can be evaluated, and the critical stations and fixture components of each MMP can be detected and modified. Liu et al. [10] proposed a quality-assured setup planning to select the optimal process plan from a group of process plan candidates with different fixture layouts. The optimal process plan was referred to as the process plan candidate that minimizes the cost related to process precision and satisfies the quality specifications. Abellan-Nebot et al. [11] proposed the use of historical shop-floor quality data from existing MMPs to extract manufacturing operation capabilities in order to conduct a more accurate process planning. The process plan selection procedure was based on three components: (a) inference on the process capabilities from shop floor data, (b) sensitivity analysis of candidate process plans to identify critical fixtures and manufacturing stations/operations, and (c) optimal selection of candidate process plans evaluating a multivariate capability ratio.

Fault cause identif ication
The issue of diagnosability refers to the problem of whether the measurements of the KPCs contain enough information for the diagnosis of critical process faults [14]. For instance, knowing the SoV model defined by Eq. 5 and measuring different KPCs (vector Y), it may be possible to infer the sources of variation (vector U). However, the MMPs are usually not diagnosable due to the inherent dimensional coupling between cutting tool deviations and fixture deviations at each machining station. That is, fixture-induced deviations and machining-induced deviations may produce the same pattern deviation of KPCs. Consequently, it is difficult to distinguish error sources at each operation.
To overcome this limitation, Wang et al. [12] applied the SoV model and proposed the equivalent fixture error concept. With this concept, datum-induced and machining-induced errors are transformed to equivalent fixture-induced errors at each operation. Using this approach, a sequential root cause identification can be conducted minimizing the number of measurements required, isolating firstly the faulty station. Assuming measurement and un-modeled noises to be negligible, Ding et al. [13] studied the diagnosability of an MMP through the definition of a diagnosability matrix. According to this matrix, three different types of diagnosability were defined: (a) diagnosability within MMP, (b) diagnosability within station, and (c) diagnosability between stations. Zhou et al. [14] extended the diagnosability analysis of MMPs when measurement and un-modeled noises are not negligible. Besides analyzing the diagnosis capability of the MMP, other research works have also studied how to identify a specific root fault cause when it is diagnosable. For this purpose, pattern recognition techniques [16] and direct estimation methods [15] have been tested. The definition of at which station an inspection of part/workpiece quality should be conducted and which features should be inspected is crucial for a successful identification of the root fault causes and process improvement. Djurdjanovic and Ni [17] proposed a Bayesian-based method to analyze the measurement schemes (i.e., placement of the inspection station and features to be inspected) in a MMP. Later, the same authors presented in [20] other non-Bayesian methods for analyzing different measurement schemes when statistical characteristics of the sensor noise term ε are known. Other research works such as [18,19] tackle the synthesis problem to define which is the optimal placement of the inspection stations for a given MMP.

Dimensional quality control
As an extension of diagnosis methodologies, some researchers have developed an in-line process adjustment technique to reduce variability in MMPs. The basic idea is to control the product quality through in-line adjustments of certain process parameters such as the fixture locations or the cutting tool path itself. The SoV model is applied to estimate the impact that those potential control actions will have on the quality of the final product. Active control for variation reduction requires two enablers [21]: in-line dimensional measurement sensors to measure actual part deviation and real-time actuators such as CNC machining stations or flexible tooling [29] to act on the manufacturing process. By these enablers, dimensional quality control can be based on feed-back control or feed-forward control [21]. Feed-back control implies that the control actions (corrections) are determined using downstream measurements obtained at the end of the process or in certain intermediate stations. This dimensional quality control can only be used to compensate mean shifts, but not to reduce variability. On the other hand, feedforward control uses in-line measurements to determine the current deviation of the workpiece in order to subsequently apply control actions to minimize the effect of this deviation in the final part quality. In this way, feed-forward control compensates current deviations instead of reacting to past deviations as feed-back control does [21].
The first work in the field of active control for variation reduction was conducted by Djurdjanovic and Zhu [22]. The feed-back and feed-forward control strategies for the placement of stations with dimensional adjustment capability were proposed. Innovatively addressing the dimension compensation problem, this work considers only the deterministic effects, neglecting the noise due to the linearization, unmodeled effects, process noise, and sensor imperfec-tion. Furthermore, the concept of compensability was introduced to quantitatively evaluate the capability of variation compensation in a specific system. Izquierdo et al. [21] extended the feed-forward control strategy to include parts/process requirements and specific engineering constraints on the magnitudes of control actions, such as physical limits and inaccuracy of tooling adjustments. These works were focused on the study of feed-forward control with full control over all tooling elements. This assumption may not be realistic, since tooling adjustments through flexible fixtures or CNC machine tools may only be assigned to selected stations in the system due to their high costs. Thus, Djurdjanovic and Ni [23] proposed a feed-forward control strategy with distributed actuation capabilities, taking into consideration the actuation accuracy and noise. However, they only select the best placement from the potential and distributed tooling adjustments, without considering the interaction of multiple tooling adjustments. Metaheuristic optimization approaches were used in [24,25], where the research work in [23] was extended to deal with variation reduction considering multiple tooling adjustments. More recently, Abellan-Nebot et al. [26] proposed a methodology to implement sensor-based fixtures in MMPs analyzing at which stations the sensor-based fixture should be installed to produce the higher compensability rate, increasing the final product quality. The paper also deals with the optimal sensor distribution within the fixture to increase the compensation capability.

Process-oriented tolerancing
Process-oriented tolerancing approach is a new tolerance approach that tries to overcome the traditional limitations of the product-oriented tolerancing approach. Unlike product-oriented tolerancing, where part tolerances are optimally allocated only considering an associated manufacturing cost from very generic process planning guidelines, process-oriented tolerancing optimally allocates tolerances of manufacturing process variables considering explicitly their associated manufacturing costs and their relationship with product quality. Basically, process-oriented tolerancing is essentially a tolerance synthesis problem where the quality specification of the final product is ensured by allocating tolerances of manufacturing process variables (such as locator tolerances) for a minimum cost. The framework of process-oriented tolerance synthesis was firstly proposed in [28,30]. In these research works, the tolerances of process variables in a MMP are optimally allocated by solving a non-linear constrained optimization problem defined by cost functions, the SoV model, a process degradation model of fixture components, a tolerance accumulation model, and several constraints related to part specifications (tolerances). Chen et al. [27] expanded the work in [28] to integrate process-oriented tolerancing with maintenance planning in multi-station assembly processes. They incorporated tool fabrication cost, fixture maintenance cost, and quality loss functions to optimize the tolerance allocation of manufacturing process variables and the frequency of fixture maintenance operations.

Fundamentals
The PD school deals with the manufacturing variation analysis in MMPs applying some of the concepts used in analyzing the geometrical variations of mechanical assemblies due to the imperfections of their components. The main idea of the PD school is to consider the manufacturing setup in a machining station as a mechanism, so the knowledge related to dimensional and geometrical variation analysis in mechanisms can be applied. In the literature, the study of the variation propagation in mechanisms can be conducted through different modeling approaches according to the nature of the model. The most common approaches applied are (a) kinematic models such as models based on SDTs [31] or vector-loop based models [32,33] and (b) degree of freedom models such as the tolerance maps models (T-Maps) [34]. The PD school has mainly applied the SDT approach to model and propagate the surface variations from parts, fixtures, and cutting tools, deriving the so-called model of the manufactured part (MoMP) [35]. The aim of the MoMP is to simulate the deviations generated in the manufacturing process considering two independent phenomena: positioning deviations and machining deviations. These deviations are accumulated over successive setups propagating the manufacturing variability. Positioning deviations are caused by fixture surface deviations and locating datum surface deviations which have been generated in previous setups. Machining deviations are caused by multiple sources of error such as geometric and kinematic errors, thermal errors, cutting force-induced errors, cutting tool wear errors, etc. The positioning and machining deviations of part surfaces are modeled by the SDT approach assuming that the expected manufacturing variations are small and parts behave as solid rigid. By this approach, dimensional and geometrical variations of manufactured parts are obtained by propagating the deviations of the elements that take part in the manufacturing process (e.g., part-holder surfaces and workpiece surfaces), modeled as a chain of small displacement torsors.
The dimensional and geometrical deviation of each element, described by an SDT, depends on the types of surfaces and tolerances involved. An SDT of a surface is composed of the small translation and orientation deviations that define the deviation between the nominal surface and the substitute surface, which is an ideal representation of the real one. For instance, a surface with a planar geometry can only present translation variations in the Z -axis and orientation variations around the X-and Y-axes, considering the normal vector of the planar surface in Z direction of its local CS. Other variations (translation in the X-and Y-axes and rotations around the Z -axis) keep the surface invariant, and thus, these deviations are considered to be undetermined. The SDT that describes the deviation between the substitute plane S i and the nominal plane N i , denoted as T S i ,N i , is thus defined as a translation deviation vector D = U · x + U · y + w · z and an orientation deviation vector = α · x + β · y + U · z. This SDT is defined as: where U is an undetermined component, w is the translation deviation in the Z -axis, and α and β are the orientation deviations around the X-and Y-axes, respectively. Similar SDTs have been defined in the literature [35][36][37] for other types of surfaces, and some of them are shown in Table 1. Note that the torsor components are constrained to keep the surface within the tolerance range. In addition to the surface torsors, it is also defined the link torsor and the part torsor. The link torsor represents the link between two substitute surfaces from different parts and shows the degrees of freedom constrained by the link (joint). The part torsor represents the part's displacement within the assembly in relation with geometric errors, joints, and its the nominal position of the part. For each part, a part torsor is defined, and for each contact between parts, a link torsor is defined.
In assemblies, the resulting deviation of a part of an assembly can be directly computed from the summation of torsors of the assembled parts that define the position of the part analyzed. For instance, consider an

Cylinder
Cylindrical annulus Generic surface Offsetting of a surface assembly of two parts, A and B. The computation of the part torsor T R,B (SDT of part B w.r.t. frame R) is evaluated for any set of joints between parts A and B (i.e., for any set of two interacting surfaces A i ,B j ) as follows [35]: where T A,A i and T B,B j are the surface torsors that represent the deviation of surface A i and B j of part A and B, respectively, in the reference frame R, T A i ,B j is the link torsor that represents the deviation of the link between surface A i and surface B j in the reference frame R, and T R,A is the SDT of part A w.r.t. frame R.
Note that T R,B is identical whatever the contact surface is considered, i.e., if there are n interacting surfaces A i − B j that compose the joints in the assembly of parts A and B; thus, there are n equations defined as Eq. 19 that result in the same value of T R,B . In the system of n linear equations, the unknown parameters are the link torsors. Obviously, different joints suppress different degrees of mobility, so for isostatic assemblies, the resulting system of linear equations can be solved.
In the MoMP, the resulting variation in part quality at the end of the MMP is obtained by evaluating the chain of torsors that influence the manufacturing performance at each station, expressed all torsors in the same CS. This chain of torsors considers the positioning and the machining deviation at each station. Following the research work in [35], the position and machining deviation can be evaluated as follows: -Positioning deviation: positioning deviation, expressed as the torsor T F k ,D , is the deviation of the nominal design reference (D) at the nominal fixture setup (F k ), and it is formulated as: In this expression, T D,S l is the SDT of locating datum surfaces (S l ) at D, and it is obtained from previous stations. Torsor T F k ,H k, j indicates the jth part-holder surface deviation at the k station, and the maximum values of the parameters of this torsor represent the part holder precision (maximum deviations expected from part-holder surfaces). Torsor T H k, j ,S l represents the relative position between the locating datum surface S l and the jth part-holder surface which depends on how both surfaces contact (joint type). The torsor T H k, j ,S l is also called the link torsor, and its parameters are called link parameters. Depending on the joint type of each pair of mating surfaces in the workpiece/part-holder assembly, different link torsors are defined [37]. Some of these torsors are shown in Table 2. The methodology to obtain the values of the link parameters for a generic workpiece/part-holder assembly (hyperstatic or isostatic assembly) is shown in detail in [38]. Note that Eq. 20 should be accomplished for any pair of locating datum and part-holder surfaces S l and H k, j that define the workpiece/part-holder assembly. -Machining deviation: machining deviation, expressed as the SDT T F k ,S i for the machining station k, is the deviation of the surface machining operation (M k,o i that generates surface S i ) at the nominal machine tool setup (i.e., the F k ) and it is formulated as: where T F k ,M k,o is the SDT of the oth machining operation due to geometrical-kinematic variations and thermal distortions of the machine tool; T M o ,M k,o i is the SDT of the machining operation due to cutting tool wear or cutting force-induced deviations when machining the surface S i . It is considered that T F k ,S i is equal to T F k ,M k,o i , that is to say, that there is an identity between the surface generated by the machining operation and the surface machined on the part. Note that the parameters of the torsors T F k ,M k,o and T M k,o ,M k,o i and their constraints (maximum values) represent the machine tools and tooling capabilities (i.e., maximum expected deviations of the cutting tool path due to machining inaccuracies). The torsor T F k ,M k,o i will be defined according to the type of the surface geometry generated and the capability of the manufacturing process. For instance, for a face milling operation that generates a planar surface parallel to the machine tool table, the following torsor will be defined: where w k M i , α k M i and β k M i are the machining deviations (translation and orientation deviations).
According to the positioning deviation and the machining deviation, the deviation of the actual part surface at the D CS in a single set-up, defined by the SDT T D,S i , can be evaluated as follows: Note that the resulting torsors T D,S i , T D,S i+1 , . . . related to the ith, (i + 1)th, . . . manufacturing features at one station may be used as inputs in subsequent stations in case that these features are used as locating datums. Thus, the deviation of machined surfaces are propagated in the subsequent stations through the SDT T F k ,D . As shown in Fig. 4, the station-by-station evaluation of all torsors defines the MoMP. The derivation of the MoMP was firstly proposed by Villeneuve et al. [35] for milling processes, considering both isostatic and hyperstatic fixtures, the later with a specific hierarchy of positioning surfaces (primary, secondary, and tertiary positioning surfaces). Vignat and Villeneuve [39] extended the formulation of MoMP for modeling manufacturing variation in turning processes considering negligible the vibration effects and the rotation defect of the lathe. The generic resolution of the positioning problem between workpiece and part holder was studied by Villeneuve and Vignat [38] providing a straightforward methodology to obtain the values of the link parameters. More recently, Nejad et al. [40] proposed the combination of the MoMP and the unified Jacobian-torsor model developed by Ghie et al. [41] for modeling mechanical assemblies. They used the interval arithmetic to associate bounds on the components of small displacement torsors in order to determine the lower an upper limit between which the actual surfaces must lie.

Virtual part quality inspection and verification
The formulation of the MoMP ends with the inclusion of the virtual inspection of the part using a virtual The gauge and the resulting part from the MoMP are assembled, and similar to the assembly process of the workpiece/part holder in the machining setup, a positioning deviation can be defined when the virtual inspection is conducted (see Fig. 4). For the virtual gauge, the gauge CS G and gauge surfaces are defined, the pth gauge surface being defined as G p . The gauge positioning deviation is defined as: where torsor T G,G p indicates the deviation of the positioning surface p of the gauge and the maximum values of the parameters of this torsor represent the gauge precision (if we assume the inaccuracy of the gauge to be negligible, T G,G p = {0 3×1 0 3×1 }), torsor T D,S m represents the deviation of the measurement datum surface S m , and the link torsor T S m ,G p represents the relative position between S m and the positioning gauge surface G p , which depends on the part/gauge assembly condition (joint type). After assembling the gauge and the final part, the functional tolerance compliance is verified by measuring the signed distance between the virtual gauge and the inspected surface S q . This distance is evaluated at the boundary points of the toleranced surface projected along the inspection direction, obtaining the distance Gap r for each rth boundary point. To calculate the Gap r distance, first, the deviation between the inspected surface and its tolerance zone (defined as T Z q ) should be calculated. This deviation is expressed with the SDT T T Z q ,S q , and it is evaluated following Eq. 25 (see [42]): where T D,S q represents the deviation of the surface to be inspected S q and T G,T Z q is the deviation of the tolerance zone of the inspected surface w.r.t. the gauge CS, which is assumed to be {0 3×1 0 3×1 } if gauge errors are negligible. Considering the SDT T T Z q ,S q as follows: the SDT that defines the deviation of the rth boundary point of the toleranced surface S q is expressed as: where t S q P r is the translation vector from S q to P r and × is the cross product operator. By knowing the SDT T T Z q ,P r , the deviation of the rth point of the toleranced surface along the direction of part verification, defined by the vector n = [n x , n y , n z ] T , is evaluated by the expression: Analyzing the distance between the point deviation and the tolerance zone, the gap distance defined as Gap r is formulated as: where τ is the maximum deviation of the rth boundary point according to the tolerance value. The rth boundary point of the inspected surface will be within the tolerance zone if Gap r remains positive or null.
In a similar way to the virtual verification by the SoV model, the virtual measurement and verification by the MoMP can be conducted following two main approaches: the worst-case approach and the statistical approach. For each approach, the resulting estimation is derived as follows: -Worst-case approach: Villeneuve and Vignat [43] reported that for a worst-case analysis, the tolerance compliance is conducted by solving an optimization problem in which the minimum gap distance from Eq. 29 at all boundary points of the toleranced surface is evaluated. This optimization problem is defined as: where Gap min = min Gap r , ∀r ∈ boundary points.
In Eq. 30, the term Gap min is the minimum distance between the virtual gauge and the toleranced surface inspected after measuring the distance at all boundary points. The expression max CGP LGP Gap min defines the inspection process, where the gauge is assembled with the part according to the standard ISO or ASME tolerance specifications shown in the design drawings. The resulting assembly depends on how the part is inspected (defined by the link parameters denoted as LGP, which are the parameters of the torsor T S m ,G p ) and the positioning limits defined by the constraints of the positioning algorithm (denoted as CGP), as explained in [38]. Within the limits of these displacements, the most favorable position for the virtual gauge relative to the part can be found by maximizing the Gap min value. In Eq. 30, material condition modifiers or incomplete datum frames can be considered in the tolerance verification, since they are related to the link parameters LGP of the gauge/part assembly and the positioning constraints CGP. The term min CM,CH,CGP DM,DH,LHP (·) is the search expression of the worst-case combination of the manufacturing defects DM, DH, LHP (machining, part holder, and workpiece-fixture assembly deviations, respectively) within the estimated range of variations expressed by the constraints CM, CH, CHP (machining, part holder, and workpiecefixture assembly constraints, respectively, which are related to machine tool and fixture capabil-ities and workpiece-fixture configurations). According to this worst-case analysis, a process plan will be considered able to satisfy the functional tolerance if Gap wc remains positive or null, which means that the deviation of the inspected surface is within the tolerance zone defined in the part drawing.
-Statistical approach: As reported above, the worstcase search is defined in Eq. 30 by the term min CM,CH,CGP DM,DH,LHP (·). For a statistical analysis, instead of conducting a search for the worst-case combination, a large number of simulations are conducted in which the sources of variation DM, DH, LHP (machining, part holder, and workpiece-fixture assembly deviations, respectively) are simulated following a specific probability distribution function. For each simulation, the max CGP LGP (Gap min ) is evaluated. After running thousands of simulations, the resulting probability distribution of the variable max CGP LGP (Gap min ) defines whether, statistically, the parts comply with the functional tolerances [44].

Tolerance analysis
The purpose of a tolerance analysis is to verify whether the design tolerance requirements can be met for a given process plan with specified manufacturing deviations. In tolerance analysis, the cumulative effect of individual variations with respect to the specified functional tolerance in all machining operations is studied in order to check a product's functionality compared with its design requirements. This is also referred to as error propagation, error stack-up, and tolerance stack-up in a MMP.
In this field, Ayadi et al. [45] applied the MoMP for MMPs based on 3-2-1 generic fixtures. They proposed an ascendant transfer method to simplify the resolution of the virtual inspection equation (e.g., the worstcase analysis given by Eq. 30) and make possible to run the tolerance analysis more straightforward. Louati et al. [46] applied the MoMP to quantify the part quality variation due to different setups in order to select the best setting solution. Tichadou et al. [47] compared the use of the MoMP and an integrated CAD/CAM system for tolerance analysis in MMPs remarking the acceptable accuracy of the MoMP according to common requirements in mechanics and the difficulty of translating the GD&T specifications to be applied within the MoMP. Nejad et al. [42] proposed a detailed mathematical formulation of tolerance analysis based on searching the worst case using two different optimization methods such as genetic algorithms and sequential quadratic programming. To simplify the resolution of Eq. 30, the optimization problem was broken down into two subproblems: the worst possible part produced according to the MMPs and the optimal inspection of one individual part. The resolution of the first subproblem is considered the input for the second subproblem. The tolerance analysis presented also considered different strategies to estimate the constraints related to the deviation torsor parameters of fixture's and machine tool's capabilities. Nejad et al. [40] studied the tolerance analysis problem using a combined approach of the MoMP and the Jacobiantorsor model. This work applied the interval arithmetic so all torsors were expressed by interval ranges. The worst-case analysis is obtained studying the error stackup on the functional elements expressed in interval ranges and makes the resolution quite rapid compared with previous methods. However, the approach is limited by the fact that it considers the torsor parameters independently, so they can reach their extreme values simultaneously which is not realistic. A detail comparison of both solution techniques presented in [40,42] was reported in [48], considering both worst-case and statistical analysis. The comparison showed that the interval approach performance is faster and accurate than the other approaches when the parameters related to the sources of error are assumed independent. Furthermore, the optimization approaches were recommended when using realistic quality constraints (parameters not independent) although they are time-consuming. Additionally, the comparison was also conducted applying different strategies to simulate the torsors related to the performance of fixtures and machine tools (constraints parameters that define the machine tool capability and fixture accuracy). These strategies were previously studied in other research works [51,52].

Tolerance synthesis
From the functional requirements of a part, the manufacturing requirements at each setup can be derived using the MoMP. In fact, from the functional inequalities system, the corresponding manufacturing inequalities are determined. These inequalities limit the defect allowed for each setup to produce a final part conform with the functional tolerance. The general tolerance synthesis problem is presented in [43]. Anselmetti and Louati [49] described in detail the tolerance synthesis considering the ISO standard. A simple algorithm is used to directly provide a complete set of manufacturing specifications in compliance with ISO standards, with the orientation and location specifications and datum reference frame. The method is applied for each functional requirements, deriving for each case the corresponding manufacturing specifications. Similarly, Vignat and Villeuve [50] studied the derivation of ISO manufacturing tolerances for each station. By their method, it can be determined if the proposed set of manufacturing tolerances is complete and if there are unnecessary manufacturing specifications which can be eliminated from the optimization problem.

Modeling example: SoV model
For illustrative purposes of the SoV model derivation, consider the part design and its associated raw material shown in Fig. 5 and the two-station machining process used in the manufacturing process shown in Fig. 6. In order to evaluate the final part variability due to both fixture-and machining-induced deviations, the SoV methodology was applied. The methodology to derive the SoV model can be summarized in the following steps: -Step 1: Define the CSs of fixtures and part surfaces.

Symbolic resolution
Applying the SoV model for the 2D case study, the resulting deviation of each KPC due to fixture-and machining-induced deviations is defined as follows: -KPC 1 : the deviation of the KPC 1 is defined as where these deviations are related to fixture and machining deviations as: Note that most of the deviations in the first station are propagated downstream, affecting the KPC 1 . However, note that locator deviations l 1 3x and l 2 3x and machining deviations along X direction (u 1 M 5 and u 2 M 6 ) do not influence on KPC 1 . -KPC 2 : the deviation of the KPC 2 is defined as where y 4,P 7A y = +22.5γ 1 M 5 + 2.5γ 2 M 7 + 0.9 l 1 1y − 0.9 l 1 2y −0.9 l 2 2y + 0.9 l 2 1y + l 2 3x − u 2 M 7 , (36) In this case, note that the locator deviations l 1 3x do not influence on this KPC; however, the same locator deviation at station 2 ( l 2 3x ) does influence. Furthermore, only machining deviations when milling surface S 7 influence on the KPC (except the deviation v 2 M 7 ), and the machining deviations when milling surface S 6 and S 5 do not influence at all. -KPC 3 : the deviation of KPC 3 is defined as where y 3,P 6A y and y 3,P 6B y are defined by Eqs. 33 and 34, respectively. By substituting, the deviation of KPC 3 is defined as |10γ 1 Note that as this KPC is a parallelism relationship, locator deviations l 1 3x and l 2 3x and translational machining deviations at any station do not influence.

Numerical resolution
The case study is numerically solved analyzing the worst case and the statistical approach. The expected variability range for each manufacturing process variable is shown in Table 5, and all manufacturing process variables are assumed to be independent to each other. For the statistical analysis, the manufacturing process variables are assumed to be normally distributed, and the ranges shown in Table 5 cover the 6σ interval. In order to validate the model, the worst-case analysis was also conducted using Pro/Engineer Wildfire 5.0. Using this software, the assembly of the workpiece and the fixture at each machining station was generated, and a material removal operation to simulate the machining operation was also defined. As a result, the final machined part was obtained as a function of the previous assemblies and operations. The worst-case analysis was then conducted by analyzing the resulting part at the extreme values of the manufacturing sources of variation and measuring all KPCs. Table 6 shows the KPC variations according to the type of analysis conducted.

±0.001
Dimensional deviations in millimeters and angular deviations in radian  Fig. 5 is used to describe the use of the MoMP. Unlike the previous example, the current two-station manufacturing process applies fixture surfaces instead of fixture locators as it is shown in Fig. 7, since this modeling approach lets model surface-based fixtures. The MoMP is built applying the methodology composed of the following steps:   Dimensions in millimeters and radian expressed in the frame D. Note that torsors should be expressed in the same frame to be summed. -Step 4: For the next kth machining station (starting from station 1), derive the following torsors:   where: The numerical resolution for worst-case analysis requires a search algorithm. In this case study, it is applied two algorithms sequentially. Firstly, a genetic algorithm (GA) in order to find a region close to the optimal solution. Secondly, the solution provided by the GA is used as the initial point in a mesh adaptative direct search algorithm for tuning the optimal result. The resolution was repeated ten times to ensure the optimization convergence, so the global minimum reached provides the worst-case solution. The numerical resolution for the statistical analysis requires the evaluation of thousands of samples through Monte Carlo simulations. For this case study, 5,000 Monte Carlo simulations were run assuming that all torsor parameters are normally distributed with the 6σ ranges shown in Table 9 subjected to their tolerance constraints. In order to validate the model, the same worst-case analysis using Pro/Engineer Wildfire 5.0 as explained in Section 4.1.2 was conducted. Table 10 shows the KPC variations according to the type of analysis conducted.

Fixture and processes supported
The SoV model is based on a linear system of equations described in a matricial form in a similar way as the well-know state-space model from control theory. For this reason, the SoV model has been used for modeling isostatic fixtures based on punctual locators according to the 3-2-1 layout. Hyperstatic fixtures based on locators or other workholding devices such as vises or four-jaw chucks have not been considered. The MoMP overcomes this limitation since the formulation lets consider non-punctual contacts for modeling the fixture assembly. Therefore, hyperstatic fixtures with a defined hierarchy of contacts and industrial fixtures such as vises or chucks can be derived.
In terms of machining systems, the SoV model deals only with milling operations since in this process, all degrees of freedom are constrained, and the application of the SoV model in its matricial form for diagnosis and process improvement is straightforward. The MoMP has been derived for both milling and turning processes, making it more applicable in real MMPs.

Limitations of virtual inspection
At the virtual inspection station, the SoV model conducts a point-based inspection similar to that conducted by a coordinate measuring machine (CMM) including in the model the measurement error in the inspection process. The virtual inspection and verification can be conducted straightforward for the worst-case and statistical analysis by evaluating one single expression (without conducting a large number of simulations or using an optimization algorithm). In the case of the MoMP, the deviation of the inspected surfaces is analyzed through a virtual gauge, so the inspection is a surface-based process where the precision of the gauge system can be considered. According to the tolerance accumulation, the MoMP can conduct both worst-case and statistical analysis. However, the worst-case analysis requires the resolution of a complex optimization problem where the global solution (and thus, the real worst-case solution) may not be reached. On the other hand, the statistical analysis requires time-consuming Monte Carlo simulations to infer, from thousands of simulated samples, the statistical results of the virtual inspection. Thus, the resolution of both worst-case and statistical analysis are quite complex in comparison with the same analysis conducted under the SoV model. This is because the SoV model includes the individual source of errors in the manufacturing process (e.g., individual locator errors in isostatic 3-2-1 fixture devices) which are independent to each other, whereas in the MoMP models, the sources of error are described by the parameters of the SDTs of surfaces which may be not independent. In other words, the SoV model is applied only in simple isostatic fixture configurations based on locators that greatly facilitates both worstcase and statistical analysis. For more complex fixtures based on surfaces, the SoV cannot be applied and its resolution should be conducted applying the MoMP solving a complex optimization problem for both worstcase and statistical cases.

GD&T conformance
The SoV model virtually inspects the parts as they were inspected in a CMM. Therefore, the SoV model is oriented to vectorial dimensioning and tolerancing (VD&T) specifications, and common product design specifications based on GD&T should be translated in order to be applied correctly. The mathematical translation of some GD&T specifications such as parallelism, angularity, and perpendicularity to be used by the SoV was investigated in [53]. Kong et al. [54] also included into the SoV model other GD&T specifications such as material condition modifiers.
Unlike SoV model, the MoMP describes surface deviations using SDTs and uses gauges for inspection. These two characteristics make the MoMP more friendly to deal with GD&T specifications than the SoV model. For instance, the MoMP has been successfully applied for manufacturing tolerancing according to ISO standards and GD&T specifications [49]. The MoMP has been also applied to deal with material condition modifiers and incomplete datum frames [39,42].

Applicability
The state-space model formulation in the SoV model makes it to be preferred to the MoMP model when the applicability of the model is mainly focused on part quality improvement activities such as fault diagnosis, process diagnosability analysis, dimensional quality control, and fixture maintenance planning. However, this applicability can be only conducted in milling processes with isostatic fixtures. In this field, the formulation of the MoMP makes it more difficult to be applied since there is no single matrix equation to be straightforward analyzed as there is in the SoV model. Unlike the SoV model, the MoMP is oriented to product design activities. This is mainly because of the ability of SDTs to model surface tolerances and because of the contact between surfaces is not only based on punctual contacts but also on surface contacts.

Model accuracy
Considering the modeling mechanism, both SoV and MoMP models are based on kinematic chains and homogeneous transformation matrices. Both assume solid rigid parts and they cannot deal with form errors (intrinsic part tolerances), so they are focused on position and orientation errors (extrinsic part tolerances). In manufacturing context, angular deviations due to manufacturing error remain very small, so the general matrix approach can be simplified using SDT or DMV [55]. Simulations conducted with Pro/Engineer are shown in Tables 6 and 10, and they reveal that, for the worst-case analysis in the case study, the estimation error of KPCs is below 1 μm for both approaches.
It should be noted that both models require a comprehensive information of the product and the manufacturing process such as geometry of nominal part surfaces, manufacturing and fixture capabili-ties, geometrical information of fixture configurations, placement of the inspection stations, capability of the inspection systems, etc. The large quantity of a priori engineering knowledge required and its level of accuracy has an important impact on the final accuracy of the model, which might prevent their application in large-scale MMPs [56]. However, for low and mediumscale systems, the application of SoV and MoMP models is highly recommended to eliminate downstream manufacturing problems and reduce ramp-up times.

Modeling complexity
The derivation of the SoV model is, in authors' opinion, more oriented to MMPs than the MoMP due to the adoption of the state-space model of control theory, and it is more easily to be automated due to its matrix formulation. The multi-station representation by the SoV model is clearer than that from the MoMP model, and the inclusion of the inspection stations is also straightforward. Furthermore, the resolution of worstcase and statistical analysis is also easier and faster since it requires only the evaluation of a system of equations in matrix form. However, it should be noted that there are more limitations applied in the SoV model, so complex MMPs based on turning-milling processes and non-punctual isostatic fixtures should require the application of the MoMP.

Conclusion and future work
The development of 3D manufacturing variation models for integrating product and manufacturing process design is essential to improve product development, eliminate downstream manufacturing problems, and reduce ramp-up times. In the literature, two main 3D manufacturing variation models in MMPs have been studied: the SoV model and the MoMP. This paper has presented both SoV and MoMP models, analyzing and comparing their main characteristics and applications. Furthermore, the methodology to derive both models has been described step by step by a simple 2D case study (extensible for any 3D part). As future work, some potential improvements on manufacturing variation propagation models in MMPs should be conducted. As some ideas, the authors suggest: -Both SoV and MoMP models are restricted to solid rigid components. Industrial machining operations also require the analysis of manufacturing variability on complaint parts due to deformations induced by clamp forces or cutting tool forces.
-Both SoV and MoMP models do not include form errors. Although form errors are considered to be not as critical as dimension or orientation errors for 3D manufacturing variation, research efforts should be conducted to include in somehow these geometrical variations. -SoV model deals with fixtures composed of locators, but their extension to cover common industrial fixtures such as vises or chucks should be addressed. -SoV model deals with non-orthogonal 3-2-1 fixture layouts, but common hyperstatic fixture schemes in industry such as N-2-1 layouts are still unaddressed. -The derivation of the MoMP includes generic machine tools capabilities by defining different torsor components. However, machine tool capabilities can be expressed by their sources of variation themselves such as cutting tool deflections, tooling wear, etc. The inclusion of additional chain of torsors to explicitly indicate the effect of these machining variations could be interesting for future applications in manufacturing tolerancing. -The application of the MoMP for process tolerance allocation is still unaddressed. This process tolerance allocation should be analyzed by integrating the different cost functions of the manufacturing process such as fixture maintenance operations or cutting tool replacement policies among others.