# Addressing Epistemic Uncertainty in Statistical Tolerance Analysis: Worst-Case Non-Conformity Rate Assessment

Sonia Carolina Garcia Gomez, Pierre Beaurepaire, Nicolas Gayton, Chahrazed Amrane

► To cite this version:

Sonia Carolina Garcia Gomez, Pierre Beaurepaire, Nicolas Gayton, Chahrazed Amrane. Addressing Epistemic Uncertainty in Statistical Tolerance Analysis: Worst-Case Non-Conformity Rate Assessment. 26e Congrès Français de Mécanique, Laboratoire d'Etude des Microstructures et de Mécanique des Matériaux (LEM3 - UMR CNRS 7239), Aug 2025, Metz, France. hal-05312883

HAL Id: hal-05312883

https://hal.science/hal-05312883v1

Submitted on 14 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Addressing Epistemic Uncertainty in Statistical Tolerance Analysis : Worst-Case Non-Conformity Rate Assessment

Sonia C. GARCIAª, Pierre BEAUREPAIREᵇ, Nicolas GAYTON, Chahrazed AMRANE

a. Université Clermont Auvergne, Clermont Auvergne INP, CNRS, Institut Pascal, F-63000 Clermont-Ferrand, France – France

e-mail : sonia.garcia gomez,pierre.beaurepaire,nicolas.gayton@sigma-clermont.fr

b. Phimeca Engineering SA, Phimeca Engineering SA, France
e-mail : amrane@phimeca.com

## Abstract :

Deviations in component dimensions are unavoidable due to manufacturing and measurement imperfections. Tolerance analysis evaluates how these variations impact the functionality of a mechanism by ensuring the stack-up of deviations remains within acceptable limits. The tolerance zone defines the permissible range of variations concerning the nominal dimension of the components, establishing boundaries for their acceptable deviations. In statistical tolerance analysis, the variability of components is represented using statistical distributions that model the spread and likelihood of component dimensions within a defined tolerance zone. The choice of distribution is crucial, as there may be an infinite number of probabilistic distributions that could fit the tolerance information. Traditional statistical tolerance methods often rely on precise probabilistic models that do not account for epistemic uncertainty due to a lack of knowledge about the exact distribution of deviations. Making arbitrary choices concerning the distribution of the deviations can affect the accuracy and the validity of the tolerance analysis. To address this issue, methods such as the APTA method [5] have been developed to incorporate ranges of variation for the mean and/or standard deviation. This paper aims to determine the worst-case non-conformity rate by accounting epistemic uncertainty to build the probability box (p-box) that represents all the possible distributions of the cumulative stack-up of deviations for serial and parallel systems. The method is exemplified using simple case studies.

## Mots clefs : Statistical tolerance analysis, non-conformity rate, epistemic uncertainty, imprecise probability

## 1 Introduction

Dimensional and geometric deviations introduced during manufacturing are inevitable and directly impact the functional performance of assembled products. To control unintended part deviations, the designer must specify and allocate tolerances to guarantee that specific quality requirements are met [9].
Nevertheless, defining tolerances is a complex task because the assigned values not only affect the functionality and manufacturing cost of the design, but they are often interrelated and contribute to a given resultant tolerance [15].

The problem of tolerance design is traditionally approached in two different ways : tolerance analysis and tolerance synthesis [10]. Tolerance synthesis is often viewed as an optimization problem where the goal is to determine the set of tolerances that either maximizes or minimizes a cost, quality, or a combined cost-quality function. In contrast, tolerance analysis focuses on evaluating the compliance of a given set of tolerances with the specified functional requirements of the products [11].

Tolerance analysis, and tolerance design in general, can be classified into two main approaches : statistical tolerance analysis and worst-case analysis [8]. The worst-case analysis assumes that all dimensional deviations will occur at their maximum allowable limits, ensuring that even under the worst combination of deviations, the assembled product will still meet its functional requirements [3]. The primary benefit of this approach is its robustness–providing deterministic outcomes and guaranteeing product performance. However, it often leads to unnecessarily tight tolerances, as it does not account for the probabilistic nature of manufacturing variations. This can increase manufacturing costs, as tighter tolerances typically require more precise, and therefore more expensive, production processes.

In contrast, statistical tolerance analysis considers a probabilistic approach, treating dimensional variations as random variables that follow defined distributions (often normal). This method estimates the likelihood that an assembly will meet functional requirements by considering the combined effects of part deviations [16]. One of its main advantages is the ability to relax tolerances, compared to worst-case analysis, which can lead to cost savings in manufacturing. It also more accurately reflects the variation typically found in production environments.

However, statistical tolerance analysis has some limitations. It does not guarantee that every part will meet specifications, making it less suitable for high-reliability applications requiring strict functional compliance [14]. Additionally, its effectiveness depends on the accuracy of statistical data, such as part variation distributions and correlations, which can be difficult to define due to limited production data. Traditional methods, like the Root Sum Squared (RSS) model [6] and second-order tolerance-analysis techniques [1], assume precise defect distribution models, which may not always be available. To overcome these challenges, more advanced methods have been developed, focusing on defining the admissible conformity space using tolerance intervals and statistical characteristics. These approaches enable the simulation of a broader range of statistical cases, offering more robust solutions by better-incorporating uncertainty in process variability [5].

Managing uncertainty in process variability inherently requires addressing epistemic uncertainty. In this context, methods based on imprecise probabilities - such as interval theory, P-Boxes, and fuzzy set theory - have gained popularity [4]. These approaches incorporate the lack of precise knowledge about the system into the analysis, ensuring that uncertainty is systematically propagated to the results. These methods allow engineers to estimate upper and lower bounds on defect probabilities, leading to more robust and realistic assessments of design tolerances.

This paper presents an approach integrating epistemic uncertainty into statistical tolerance analysis using probability boxes. The proposed methodology estimates the worst-case non-conformity rate of a system, without requiring precise assumptions about defect distributions, making it more adaptable to real-world manufacturing variability. To demonstrate its practical benefits in mechanical design, the approach is
applied to two real-world case studies, highlighting its effectiveness in addressing uncertainty and improving reliability.

The article is subdivided into three major parts. Firstly, Section 2 presents the framework of the statistical tolerance analysis problem, focusing on the propagation of variability through the stack-up (assembly response) function and the general calculation of the non-conformity rate. Then, Section 3 introduces imprecise probability analysis for obtaining the worst-case non-conformity rate for each single component and, in general, for an assembly response function. Section 4 illustrates the use of imprecise probability for calculating the worst-case non-conformity rate in two case studies and compares the results with those presented in [5]. Finally, the discussion and conclusions are presented. It should be noted that, unless otherwise specified, all dimensions are assumed to be in millimeters (mm).

## 2 Framework of Statistical Tolerance Analysis

From a mathematical point of view, the functional relationship between the tolerances of the individual components and the behaviour of the final assembly can be expressed as :

$$Y = f(X_1, X_2, \dots, X_n) = f(\mathbf{X}) \qquad (1)$$

where, Y is the assembly response function (also known as the performance characteristic), and X are random variables representing the individual dimensions or tolerances of the components. The values of each random variable $X_i$ are typically derived from the drawing dimensions and lie within a specified range, bounded by a lower specification limit $LSL_i$ and an upper specification limit $USL_i$. The assembly quality is determined by whether Y meets the functional condition, usually represented by functional variable bounds ($LSL_Y, USL_Y$). The non-conformity rate $P_{nc}^Y$ represents the probability that Y falls outside its acceptable range, which depends on the variability of the random variables. To compute this non-conformity rate, it is necessary to know the probability density functions (PDFs) of each $X_i$. Since these distributions are not always readily available, certain assumptions must be made [13]. Some models that can be considered are : uniformly distributed statistical model, non-shifted distribution, shifted distribution and dynamic shifted distribution. These models are briefly explained in [5].

In most of the cases, each variable $X_i$ is assumed to follow a normal distribution. Under the non-shifted distribution model, the nominal value corresponds to the mean, and the standard deviation is defined in terms of the component tolerance $T_i$ as :

$$\sigma_i = \frac{T_i}{k} \qquad (2)$$

where k is a constant associated with the quality, and its value influences the percentage of realizations within the specified tolerance zone. For instance, if it is set to 6, approximately 99.9% of the values of $X_i$ will fall within the tolerance limits.

Manufacturing processes inherently exhibit natural variability, which significantly impacts the quality of a batch of components in mass production [12]. The natural process tolerance is generally defined as ±3 times the standard deviation of the manufacturing process [13]. Consequently, the specified tolerance zone $T_i$ for each component $X_i$ must be smaller than the natural tolerance $t_i$ associated with its respective manufacturing process, hence the capability index of the manufacturing process $C_{p,i}$ can be calculated as follows :
$$C_{p,i} = \frac{\text{Tolerance width}}{\text{Process variation}} = \frac{T_i}{t_i} = \frac{T_i}{6\sigma_i} \qquad (3)$$

where $\sigma_i$ is the standard deviation of a batch of $X_i$. Due to the natural process variability, the mean of the batch can be shifted. This mean shift can be modelled through the constant $C_{pk,i}$ which also includes the effect of the absolute mean shift $\delta_i$ from the nominal value [7] :

$$C_{pk,i} = \frac{T_i/2 - |\delta_i|}{3\sigma_i} \qquad (4)$$

Since the capability of a manufacturing process directly influences the variability of $X_i$, the minimum capability indices $C_{p,i}^{(r)}$ and $C_{pk,i}^{(r)}$ must be established as requirements. Using Equations (3)-(4), the conformity region for any $X_i$ can be determined. Figure 1a illustrates the conformity region for a given variable $X_i$ when the feasibility constraints are :

$$C_{p,i}^{(r)} = C_{pk,i}^{(r)}, \quad C_{p,i}(\sigma_i) = C_{p,i}^{(r)}, \quad C_{pk,i}(\sigma_i, \delta_i) = C_{pk,i}^{(r)}, \quad \text{and} \quad C_{p,i}^{(r)} < C_p^0 < \infty \qquad (5)$$

It is important to note that while $C_{pk,i}^{(r)}$ is not required to be equal to $C_{p,i}^{(r)}$, it must not be greater than it. When $C_{pk,i}^{(r)} < C_{p,i}^{(r)}$, the representation of the conformity region is no longer a triangle but a trapezoid, see Figure 1b. Furthermore, as $C_p^0$ tends to infinity, the maximum mean shift $\delta_{max}(\sigma_{min})$ approaches $t/2$, while $\sigma_{min}$ tends to zero. A batch of $X_i$ is considered compliant if its $\sigma_i$ and $\delta_i$ associated are inside the conformity region.

## 2.1 Propagation of variability and the non-conformity rate of Y

Once the conformity region for each of the individual process variables is determined, it is important to understand how the variability in these input variables propagates through the functional relationship. The distribution of Y depends on the nature of the function $f(X)$ and the statistical properties of the individual input variables.

For a linear relationship, where $Y = \sum_{i=1}^{n} X_i$, the resulting distribution of Y will also follow a normal distribution with the mean and variance determined by the mean and variance of each $X_i$ :

$$\mu_Y = \sum_{i=1}^{n} (\mu_i + \delta_i) \\ \sigma_Y^2 = \sum_{i=1}^{n} \sigma_i^2 \qquad (6)$$

In the case of having a non-linear function, the distribution of Y cannot be easily derived, and methods such as Monte-Carlo simulation are often employed to estimate it.

The non-conformity rate $P_{nc}^Y$ of Y is then defined as the probability that it falls outside the specified tolerance limits for the final product. Specifically, given the functional limits $LSL_Y$ and $USL_Y$ the non-conformity rate is calculated as :

$$P_{nc}^Y = P(Y \leq \text{LSL}_Y) + P(Y \geq \text{USL}_Y) \qquad (7)$$
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Conformity region or Variability Domain (VD).</figcaption></figure>

To calculate this, the probability distribution of Y must be known. If Y is normally distributed, the non-conformity rate is given by :

$$ 
\begin{align} 
P(Y \le LSL_Y) &= \Phi \left( \frac{LSL_Y - \mu_Y}{\sigma_Y} \right) \tag{8} \\ 
P(Y \ge USL_Y) &= 1 - \Phi \left( \frac{USL_Y - \mu_Y}{\sigma_Y} \right) 
\end{align} 
$$

where Φ is the Cumulative Density Function (CDF) of the normal distribution. These calculations allow for an estimation of the likelihood that the assembly will fail to meet the required specifications.

When Equation (1) takes a non-linear form, the distribution of Y is not normal. However, the non-conformity probability $P_{nc}^Y$ can always be estimated using Monte Carlo simulations. While this method provides accurate results, it is computationally expensive. To address this, some authors suggest alternative approaches, such as the First-Order Reliability Method (FORM) [2] or the Second-Order Taylor Approximation (SOTA) method. The latter is based on a second-order Taylor expansion of Y and approximates its distribution using the generalized lambda distribution. Then, knowing the mean and standard deviation of each component, the calculation of $P_{nc}^Y$ is performed using the cumulative density
function of the generalized lambda distribution [1].

In [5], the authors propose a methodology to take into account the process variability in order to calculate the non-conformity rate $P_{nc}^Y$:

$$P_{nc}^Y = \int_{VD} P_{nc|\delta,\sigma}^Y(\sigma, \delta) \prod_{i=1}^{n} h_{\delta,\sigma}(\delta_i, \sigma_i) d\delta_i d\sigma_i \quad (9)$$

where $VD$ represents the variability domain within the conformity region (Figure 1), $n$ is the number of variables $X$, and $h_{\delta,\sigma}(\delta_i, \sigma_i)$ denotes the joint probability density function of the mean shifts $\delta_i$ and standard deviations $\sigma_i$.

When there is not explicit knowledge of $h_{\delta,\sigma}(\delta_i, \sigma_i)$, due to dependencies on process conditions and potential correlations, it is still possible to calculate the upper bound of the non-conformity rate:

$$P_{nc}^{Y(U)} = \max_{(\delta, \sigma) \in VD} P_{nc|\delta,\sigma}^{Y}(\sigma, \delta) \quad (10)$$

The former equation particularly useful, as it allows to estimate the worst-case non-conformity rate without the need for precise modeling.

# 3 Imprecise Probability Analysis

As already stated, due to the natural variability in manufacturing processes, the parameters that define each $X_i$ cannot be known exactly, moreover, they can fluctuate. When the requirements for the capability indices, $C_{pk,i}^{(r)}$ and $C_{p,i}^{(r)}$, are set, and the maximum capability of the process, $C_{p,i}^{(0)}$, is known, it becomes possible to calculate the ranges of variation for each $\sigma_i$ and $\delta_i$. With this information it is possible to build the credal set that includes the unknown CDF of each $X_i$. In the same way, knowing the interaction between the $X_i$ it is then possible to calculate the credal set of $Y$. P-Boxes are commonly used in engineering applications to represent a set of possible probability distributions for $F(X_i)$ bounded by a lower $F(X_i)^L$ and an upper $F(X_i)^U$ CDF. P-Boxes are computationally efficient, easy to build, and provide a simple graphical representation [4].

## 3.1 Modelling a single component

Consider a component $X_i$ of dimension 6. $\pm 0.3$, with $C_p^{(r)} = 0.5$, $C_{pk}^{(r)} = 0.2$ and $C_p^{(0)} = 1.5$. Since $C_p \ge C_{pk}$, the representation of the variability domain is a trapezoid, see Figure 1b. The P-Box that represents the credal set of $X_i$, Figure 2a, is bounded by the CDFs associated with the mean shift and standard deviations of the vertices, see Figure 3, in other terms, the upper bound and the lower bound of the P-Box of $X_i$ can be calculated as:

$$ 
\begin{aligned}
F(X_i)^U &= \begin{cases} \Phi\left(\frac{x-(\mu_i-\delta_{\min,i}(\sigma_{\max,i}))}{\sigma_{\max,i}}\right), & \text{if } x < \mu_i - 3C_{p,i}^{(r)}\sigma_{\max,i} \\ \Phi\left(\frac{x-(\mu_i-\delta_{\max,i}(\sigma_{\min,i}))}{\sigma_{\min,i}}\right), & \text{if } x \ge \mu_i - 3C_{p,i}^{(r)}\sigma_{\max,i} \end{cases} \\
F(X_i)^L &= \begin{cases} \Phi\left(\frac{x-(\mu_i+\delta_{\max,i}(\sigma_{\min,i}))}{\sigma_{\min,i}}\right), & \text{if } x < \mu_i + 3C_{p,i}^{(r)}\sigma_{\max,i} \\ \Phi\left(\frac{x-(\mu_i+\delta_{\min,i}(\sigma_{\max,i}))}{\sigma_{\max,i}}\right), & \text{if } x \ge \mu_i + 3C_{p,i}^{(r)}\sigma_{\max,i} \end{cases}
\end{aligned} 
\quad (11) $$
<figure><img src="image_3.png" /><figcaption>FIGURE 2 – P-Box for a variable $X_i$ representing all the possible CDF within the variability domain for $C_{p,i}^{(r)} = 0.5$, $C_{pk,i}^{(r)} = 0.2$ and $C_{p,i}^{(0)} = 1.5$. In (b), the square markers represent the results of the maximization of the CDF, while the triangle markers represent the results of the minimization.</figcaption></figure>

when $C_{pk,i}^{(r)} = C_{p,i}^{(r)}$ the former set of equations is simplified because $\delta_{min,i}(\sigma_{max,i}) = 0$.

When the calculation of the P-Box is not straightforward, numerical methods for propagating P-Boxes can be used [4]. In our case, we use the double-loop approach, in which $F_{X_i}^U$ and $F_{X_i}^L$ are obtained by solving an optimization problem that maximizes and minimizes the CDF of $X_i$, respectively :

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Association of the CDFs that constitute the P-Box for a variable $X_i$, with their respective mean shifts and standard deviations from the Variability Domain (the square and triangular markers are associated with the upper bound and the lower bound, respectively).</figcaption></figure>
$$ F_{X_i}^{U,\text{opt}} = \max F_{X_i}(\sigma_i, \delta_i) \\ \text{subject to :} \quad \sigma_{\min,i} \le \sigma_i \le \sigma_{\max,i}, \\ \qquad -\delta_{\max,i}(\sigma_i) \le \delta_i \le \delta_{\max,i}(\sigma_i) \tag{12} $$

$$ F_{X_i}^{L,\text{opt}} = \min F_{X_i}(\sigma_i, \delta_i) \\ \text{subject to :} \quad \sigma_{\min,i} \le \sigma_i \le \sigma_{\max,i}, \\ \qquad -\delta_{\max,i}(\sigma_i) \le \delta_i \le \delta_{\max,i}(\sigma_i) \tag{13} $$

When calculating the P-Box via optimization for the previous case, we obtain a good approximation to the analytical solution (see Figure 2b).

## 3.2 Illustrative example : Two-part assembly

The computation of the P-Box of Y for any assembly response function can be done using the double loop approach. To illustrate the process, consider the mechanical assembly of two parts shown in Figure 4. In this case, the stack-up function of the system can be written as :

$$ Y = X_1 + X_2 \tag{14} $$

where the nominal dimensions of the components $X_1$ and $X_2$ are 6 and 4, respectively. The tolerances of the components are equal, with $t_1 = t_2 = 0.6$. Their capability requirements are also equal, set as $C_{p1}^{(r)} = C_{p2}^{(r)} = 1$, $C_{pk1}^{(r)} = C_{pk2}^{(r)} = 1$ and $C_{p1}^{0} = C_{p2}^{0} = 1.5$. The mechanical requirement for the assembly can be written as :

$$ Y = X_1 + X_2 \in [9.5, 10.5] \tag{15} $$

where 9.5 and 10.5 represent the functional variable bounds, LSLY and USLY, respectively. Figure 5a shows the calculated P-Box for Y given the tolerances and the capability requirements. The worst-case non-conformity rate can be calculated by means of an optimization process

$$ P_{\text{nc}}^{Y(U)} = \max(F_Y(\sigma_i, \delta_i, \text{LSL}_Y) + (1 - F_Y(\sigma_i, \delta_i, \text{USL}_Y))) \\ \text{subject to :} \quad \sigma_{\min,i} \le \sigma_i \le \sigma_{\max,i}, \\ \qquad -\delta_{\max,i}(\sigma_i) \le \delta_i \le \delta_{\max,i}(\sigma_i) \tag{16} $$

It is important to note that, in this case, the previous optimization process, Equation (16) differs from the optimization used to generate the P-Box of the stack-up function (Equations (12)-(13)). As seen in

<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Simple mechanical assembly of two parts in a serial configuration</figcaption></figure>
Equation (15), the functional condition includes both a lower and an upper bound, meaning the worst-case non-conformity rate is influenced by both boundaries of the P-Box.

For the functional condition, Equation (15), the worst-case non-conformity rate $P_{nc}^{Y(U)}$ is approximately 731 ppm. In this case, the mean shift and standard deviation for both variables associated with the calculation of $P_{nc}^{Y(U)}$ are presented in Figure 6.

Figure 5b illustrates the evolution of the maximum non-conformity rate as the functional variable bounds, $LSL_Y$ and $USL_Y$, vary—expanding the range of the functional condition. Specifically, $LSL_Y$ ranges from [9.6, ..., 9.2] and $USL_Y$ from [10.4, ..., 10.8], resulting in a tolerance zone for $Y$ of [0.8, ..., 1.6], computed as $S = USL_Y - LSL_Y$.

In Figure 6, we observe how the mean shift and standard deviation associated with the maximum non-conformity rate evolve, shifting toward the upper vertex of the triangle when the tolerance zone is less than 1.1. This is noteworthy because it suggests that the worst-case non-conformity rate is not always linked to the minimum standard deviation and an extreme mean shift. Instead, an alternative configuration may yield an even higher non-conformity rate, indicating that the optimization problem becomes nontrivial in more complex systems.

# 4 Case studies

Let's consider the two industrial applications presented in [5]. Both industrial cases are provided by RADIALL; the first one is a linear application, while the second one is non-linear.

## 4.1 Linear industrial application

This example consists of an electrical pin whose functionality depends on a variable $Y$ that represents the contact length of the pin in its socket. $Y$ depends on 17 dimensions linked together by the following stack-up function :

$$Y = f(X_i) = \sum_{i=1}^{17} k_i X_i, \quad k_i = \pm 1 \qquad (17)$$

<figure><img src="image_6.png" /><figcaption>FIGURE 5 – (a) P-Box of $Y$ for the simple mechanical assembly of two parts in a serial configuration, and (b) sensitivity of the non-conformity rate regarding the functional condition $S$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIGURE 6 – Values of mean shift and standard deviation of each variable that generates the $P_{nc}^{Y(U)}$ for different functional conditions.</figcaption></figure>

The functional requirement of the system is to ensure that the contact length of the pin (Y) is greater than or equal to a given threshold S, $Y \ge S$. This requirement is considered the minimum contact length needed to ensure good conductivity in an electrical circuit. The values of $k_i$, the nominal dimension of each X, their tolerances, and the capability requirements are given in Table 1. The parts $X_i$ for $i = [4, 5, 12, 13]$ are considered difficult to manufacture within tolerance; hence, each dimension is verified after fabrication. This consideration is taken into account in the optimization process by setting

<table><caption>TABLE 1 – RADIALL Linear application : Dimension details, from [5].</caption><thead><tr><th>Dim.</th><th>$k_i$</th><th>Nominal Dimension</th><th>$t_i$</th><th>$C_{p,i}^{(r)} = C_{pk,i}^{(r)}$</th><th>$C_{p,i}^{0}$</th></tr></thead><tbody><tr><td>1</td><td>-1</td><td>0.750</td><td>0.20</td><td>1.8</td><td>2.6</td></tr><tr><td>2</td><td>1</td><td>10.710</td><td>0.16</td><td>1.8</td><td>2.6</td></tr><tr><td>3</td><td>1</td><td>3.750</td><td>0.10</td><td>1.8</td><td>2.6</td></tr><tr><td>4</td><td>-1</td><td>5.125</td><td>0.05</td><td>n.a</td><td>n.a</td></tr><tr><td>5</td><td>-1</td><td>3.125</td><td>0.05</td><td>n.a</td><td>n.a</td></tr><tr><td>6</td><td>1</td><td>4.505</td><td>0.03</td><td>1.7</td><td>1.94</td></tr><tr><td>7</td><td>1</td><td>5.985</td><td>0.09</td><td>1.8</td><td>2.8</td></tr><tr><td>8</td><td>-1</td><td>21.535</td><td>0.11</td><td>0.38</td><td>1.32</td></tr><tr><td>9</td><td>-1</td><td>8.195</td><td>0.11</td><td>0.38</td><td>1.32</td></tr><tr><td>10</td><td>1</td><td>5.985</td><td>0.09</td><td>1.8</td><td>2.8</td></tr><tr><td>11</td><td>1</td><td>12.305</td><td>0.03</td><td>1.7</td><td>1.94</td></tr><tr><td>12</td><td>-1</td><td>10.325</td><td>0.05</td><td>n.a</td><td>n.a</td></tr><tr><td>13</td><td>-1</td><td>5.125</td><td>0.05</td><td>n.a</td><td>n.a</td></tr><tr><td>14</td><td>1</td><td>3.750</td><td>0.10</td><td>1.8</td><td>2.6</td></tr><tr><td>15</td><td>1</td><td>0.825</td><td>0.05</td><td>1.05</td><td>1.85</td></tr><tr><td>16</td><td>1</td><td>8.800</td><td>0.10</td><td>1.8</td><td>2.6</td></tr><tr><td>17</td><td>-1</td><td>0.400</td><td>0.20</td><td>1.8</td><td>2.6</td></tr></tbody></table>
$\sigma_i = 0$ and $\delta_i \in [-t_i/2, t_i/2]$ for these four variables.

The P-Box of Y can be calculated through the double-loop optimization, as shown in Figure 7a. Since the functional condition limits only the lower bound of Y, we can calculate the maximum non-conformity rate by considering only the upper bound of the P-Box:

$$P_{\text{nc}}^{Y(U)} = F_{Y}^{U,\text{Opt}}(S) \qquad (18)$$

To perform this optimization, we employ a local optimization algorithm initialized at 100 different random starting points. This approach helps mitigate the risk of converging to a local minimum. For $S = 1.75$, the maximum non-conformity rate is approximately 98.7%. This corresponds to the mean shift and standard deviations presented in Figure 8a. These results are consistent with those presented in [5], with slight differences potentially arising from the optimization process.

Additionally, Figure 7b illustrates the evolution of the maximum non-conformity rate as a function of S, highlighting its sensitivity to this parameter. For $S < 1.65$, the maximum non-conformity rate values are higher than those reported in [5]. This discrepancy arises because that study assumes the worst-case scenario always occurs at the lower vertex of the variability domain (smallest standard deviation and for the mean shift $-k_i\delta_{\text{max},i}(\sigma_{\text{min},i})$). However, our results demonstrate that this is not always the case. As shown in Figure 8b, for $S < 1.65$, the mean shift and standard deviation that yield the worst-case value can sometimes be associated with the upper vertex of the variability domain (components $X_8, X_9$).

## 4.2 Non-linear industrial application

This application involves modeling the amplitude of a pin's extremity, which can lead to assembly issues in an electrical system, see Figure 9. The function describing this variation depends on 14 dimensions and follows the relationships provided below:

<figure><img src="image_8.png" /><figcaption>FIGURE 7 – (a) P-Box of Y for the linear case representing all the possible CDF taking into account the variability of all the $X_i$, and (b) the sensitivity of the non-conformity rate regarding the value of S.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIGURE 8 – Mean shift and standard deviation values that generate the maximum non-conformity rate ($P_{nc}^{Y(U)}$) for different values of $S$.</figcaption></figure>
$$c = \frac{cm_4 + cm_8 + cm_7}{2}$$

$$i = \frac{ima_{17} + ima_{19} + ima_{20}}{2}$$

$$h = ima_2 + ima_3$$

$$\alpha = \arccos\left(\frac{c}{\sqrt{i^2 + h^2}}\right) - \arccos\left(\frac{i}{\sqrt{i^2 + h^2}}\right)$$

$$r2 = \frac{1}{\tan(\alpha)} \left( \frac{ima_{19} + ima_{20}}{2} - \frac{cm_8 + cm_7}{2 \cos(\alpha)} \right)$$

$$z = \frac{r_2}{\cos(\alpha)} + \left( \frac{cm_9 + cm_{10}}{2} + \frac{cm_7}{4} \right) \tan(\alpha)$$

$$J1 = (cm_1 - cm_3 - z) \sin(\alpha)$$

$$J2 = \frac{cm_7}{4} \cos(\alpha)$$

$$J3 = \frac{cm_5 + cm_6}{2} \cos(\alpha)$$

$$Y = J1 + J2 + J3$$

The information regarding component dimensions, tolerances, and process capability indices is presented in Table 2. The functional requirement of the system is to ensure that the amplitude of the pin's extremity is below a threshold S, Y < S.

Since the stack-up function of the system is nonlinear, we will employ the SOTA methodology to construct the P-Box and compute the non-conformity rate. Figure 10a illustrates the P-Box of the system. Given that the controlled bound in the functional condition is the upper bound of Y, the worst-case non-conformity rate must be calculated using the lower bound of the P-Box. Thus:

$$P_{nc}^{Y(U)} = 1 - F_{Y}^{L, Opt}(S) \qquad (19)$$

<figure><img src="image_10.png" /><figcaption>FIGURE 9 – Non-linear industrial case, reproduced from [5].</figcaption></figure>
<table><caption>TABLE 2 – RADIALL Non-Linear application : component dimensions and process capability indices, from [5].</caption><thead><tr><th>Dim.</th><th>Nominal Dimension</th><th>t<sub>i</sub></th><th>C<sup>(r)</sup><sub>p,i</sub> = C<sup>(r)</sup><sub>pk,i</sub></th><th>C<sup>0</sup><sub>p,i</sub></th></tr></thead><tbody><tr><td>cm1</td><td>10.53</td><td>0.2</td><td>1.1</td><td>1.6</td></tr><tr><td>cm3</td><td>0.75</td><td>0.04</td><td>1.1</td><td>1.6</td></tr><tr><td>cm4</td><td>0.643</td><td>0.015</td><td>1.1</td><td>1.6</td></tr><tr><td>cm5</td><td>0.1</td><td>0.2</td><td>1.1</td><td>1.6</td></tr><tr><td>cm6</td><td>0</td><td>0.06</td><td>1.1</td><td>1.6</td></tr><tr><td>cm7</td><td>0</td><td>0.2</td><td>1.1</td><td>1.6</td></tr><tr><td>cm8</td><td>0.72</td><td>0.04</td><td>1.1</td><td>1.6</td></tr><tr><td>cm9</td><td>1.325</td><td>0.05</td><td>1.1</td><td>1.6</td></tr><tr><td>cm10</td><td>0</td><td>0.04</td><td>1.1</td><td>1.6</td></tr><tr><td>ima2</td><td>3.02</td><td>0.06</td><td>0.86</td><td>1.36</td></tr><tr><td>ima3</td><td>0.4</td><td>0.06</td><td>0.86</td><td>1.36</td></tr><tr><td>ima17</td><td>0.72</td><td>0.04</td><td>0.86</td><td>1.36</td></tr><tr><td>ima19</td><td>0.97</td><td>0.04</td><td>0.86</td><td>1.36</td></tr><tr><td>ima20</td><td>0</td><td>0.04</td><td>0.86</td><td>1.36</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>FIGURE 10 – (a) P-Box of Y for the non-linear case representing all the possible CDF taking into account the variability of all the $X_i$, and (b) the sensitivity of the non-conformity rate regarding the value of S.</figcaption></figure>

For $S = 0.52$, the maximum non-conformity rate is approximately 4.1%, corresponding to the mean shifts and standard deviations presented in Figure 11. Figure 10b illustrates the evolution of the maximum non-conformity rate as a function of $S$, showing that it decreases as $S$ increases. Since this is a non-linear case, the cumulative distribution function (CDF) is computed through the SOTA methodology, using a generalized lambda distribution (GLD). The resulting CDF serves as input for our optimization algorithm—a local optimization approach initialized at 100 different random starting points. However, this method becomes computationally expensive, primarily due to the number of points required for GLD estimation. In our case, $2^{14}$ samples are generated using a Sobol sequence for the 14 variables. Additionally, as we approach the tail of the distribution, the optimization algorithm tends to get trapped in local minima, making it increasingly difficult to determine the worst-case non-conformity rate beyond $S = 0.55$.
<figure><img src="image_12.png" /><figcaption>FIGURE 11 – Values of mean shift and standard deviation of each variable that generates the $P_{nc}^{Y(U)}$ for different values of $S$ in the non-linear case.</figcaption></figure>

If we assume the worst-case configuration to be the one presented in Figure 11 for all values of $S$, we obtain results that align with our optimization outcomes (up to $S = 0.55$) and with those reported in [5]. Slight differences from [5] may arise because we employ the SOTA methodology, whereas the authors in [5] used the FORM approach to calculate the non-conformity rate, along with different numerical approximations.

# 5 Discussion

From the illustrative example of the two-part assembly, we observe that determining the worst-case non-conformity rate of a system is not necessarily straightforward. Unlike previous studies, which assumed a fixed relationship between standard deviations and mean shifts, our results demonstrate that these parameters can vary depending on the system's functional condition, as shown in Figure 6.

This observation is further reinforced in the linear industrial application. As we move further into the tail of the stack-up function's distribution, the mean shift and standard deviation of certain parts transition from the lower-left vertex to the upper vertex of the variability domain, as shown in Figure 8. This shift highlights the complexity of accurately identifying worst-case scenarios.

In the non-linear case study, we successfully computed the worst-case non-conformity rate for certain values of $S$. However, as we progress further into the tail of the distribution, the optimization algorithm encounters difficulties in identifying the optimal solution, Figure 10b. This occurs because the algorithm is more likely to become trapped in local minima.

To validate our approach, we compared our results with those reported in [5] for the same case studies. These comparisons indicate that our methodology has significant potential. However, to improve convergence in the non-linear case, it would be necessary to implement a global search algorithm—such as simulated annealing, genetic algorithms, or other heuristic methods—to enhance the optimization process
and address these challenges effectively.

# 6 Conclusions

In this work, we introduced an approach that integrates epistemic uncertainty into statistical tolerance analysis using P-Boxes. The methodology enables the estimation of worst-case non-conformity rates without requiring precise assumptions about defect distributions.

As an initial step, we applied the approach to a basic serial system consisting of two parts. We computed the corresponding P-Box and analyzed the maximum non-conformity rate under different functional conditions. This analysis revealed that the worst-case non-conformity rate does not always correspond to the same combination of standard deviation and mean shift for the components. This insight suggests that, in more complex systems, the optimization process will play an even more critical role in accurately identifying worst-case scenarios.

Furthermore, in both the linear and non-linear case studies reproduced from [5], we observed some differences in the maximum non-conformity rate compared to their results. This may be due to the fact that, when calculating the maximum non-conformity rate, the previous study considered the worst-case to be the smallest standard deviation with the largest mean shift. However, our analysis shows that this is not always the case.

The results from the case studies show that the pairs of mean shifts and standard deviations associated with each component $X_i$, that generate the maximum non-conformity rate are located at a vertex of the variability domain for each component. However, we cannot ensure that this holds for more complex systems, particularly parallel or over-constrained systems, where the worst-case configurations may follow a different pattern. Investigating this further would be an interesting topic for future research, as it could provide deeper insights into worst-case non-conformity rate behavior in other architectures.

# Acknowledgments

This work was supported by the French National Research Agency (Agence Nationale de la Recherche, ANR) under project number ANR-21-CE46-0009.

# Références

[1] A Second-Order Method for Assembly Tolerance Analysis, volume Volume 1 : 25th Design Automation Conference of International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, 09 1999.

[2] Karl Breitung. Asymptotic approximations for multinormal integrals. *Journal of Engineering Mechanics*, 110(3):357-366, 1984.

[3] Vincent Delos, Denis Teissandier, Alexander Malyshev, Grégory Nuel, and Sonia García. Polyhedral-based Modelling and Algorithms for Tolerancing Analysis. *Computer-Aided Design*, page 103071, Jul 2021.

[4] Matthias G.R. Faes, Marco Daub, Stefano Marelli, Edoardo Patelli, and Michael Beer. Engineering analysis with probability boxes : A review on computational methods. *Structural Safety*, 93:102092, 2021.
[5] Gayton, Nicolas, Beaucaire, Paul, Bourinet, Jean-Marc, Duc, Emmanuel, Lemaire, Maurice, and Gauvrit, Laurent. Apta : advanced probability-based tolerance analysis of products. *Mécanique & Industries*, 12(2):71-85, 2011.

[6] W. H. Greenwood and K. W. Chase. Root sum squares tolerance analysis with nonlinear problems. *Journal of Engineering for Industry*, 112(4):382-384, 11 1990.

[7] Mikel J Harry and Reigle Stewart. *Six sigma mechanical design tolerancing*. Motorola University Press Schaumburg, IL, 1988.

[8] Y.S. Hong and T.C. Chang. A comprehensive review of tolerancing research. *International Journal of Production Research*, 40(11):2425-2459, 2002.

[9] Zbigniew Humienny. State of art in standardization in the geometrical product specification area a decade later. *CIRP Journal of Manufacturing Science and Technology*, 33:42-51, 2021.

[10] Sandipan Karmakar and Jhareswar Maiti. A review on dimensional tolerance synthesis : Paradigm shift from product to process. *Assembly Automation*, 32:373-388, 2012.

[11] R. Sampath Kumar, V. Soundararajan, and N. Alagumurthi. Review of tolerance analysis, al- location and constraints in manufacturing. *Journal for Manufacturing Science and Production*, 10(1):1-16, 2009.

[12] E. M. Mansoor. The application of probability to tolerances used in engineering designs. *Proceed- ings of the Institution of Mechanical Engineers*, 178(1):29-39, 1963.

[13] Swami D Nigam and Joshua U Turner. Review of statistical approaches to tolerance analysis. *Computer-Aided Design*, 27(1):6-15, 1995.

[14] Jami J. Shah, Gaurav Ameta, Zhengshu Shen, and Joseph Davidson. Navigating the tolerance analysis maze. *Computer-Aided Design and Applications*, 4(5):705-718, 2007.

[15] Michael Walter and Sandro Wartzack. Statistical tolerance-cost-optimization of systems in motion taking into account different kinds of deviations. In Michael Abramovici and Rainer Stark, editors, *Smart Product Engineering*, pages 705-714, Berlin, Heidelberg, 2013. Springer Berlin Heidelberg.

[16] D.E. Whitney, O.L. Gilbert, and M. Jastrzebski. Representation of geometric variations using matrix transforms for statistical tolerance analysis in assemblies. *Research in Engineering Design*, 6:191-210, 1994.