**2. Process Description**

The process considered in this work regards a surface mount industrial unit (Bosch Car Multimedia Portugal), composed by several lines where electronic components are attached to PCBs through reflow soldering. In this process, SPDs are accurately placed in designated points of copper PCBs that will be used to fix and functionalize electronic components further ahead in the assembly line, through reflow soldering under a predefined temperature profile. Several thousands of SPDs are set in place, with a cycle time that can achieve approximately 20s. The resulting PCBs with printed SPDs are subject to 100% 3D Solder Past Inspection (3DSPI). This technology has been a pillar for process monitoring and improvement in this plant, and will be further explored in this work. 3DSPI generates, for each single SPD printed in the pad, a full 3D profile, retaining 5 features for analysis: area (*a*), height (*h*), volume (*v*), o ffset in the X-direction (*x*) and o ffset in the Y-direction (*y*).

In this work, we consider the production of one specific product from which 3507 SPDs (we will also use the designation of "pads" to refer to such SPDs) are monitored, totaling 3507 × 5 = 17,535 variables to be handled. Three datasets are available and will be used in this study to test the proposed approach: CS1 contains information for 337 NOC PCBs and is used to validate the AGV methodology and confirm/tune some of its parameters; CS2 and CS3, contain 2080 and 1330 PCBs, respectively, mostly regarding good products, but also containing some PCBs with a few faulty pads. CS2 and CS3 will be used to further test the proposed methodology.

#### **3. Artificial Generation of (Common Cause) Variability (AGV)**

The proposed methodology for Artificial Generation of (common cause) Variability (AGV) aims at completing the reference dataset with the long term dominant sources of common cause variability that, from engineering knowledge, are known to exist. This is done by simulating the printing of solder paste deposits in Printed Circuit Boards (PCB) under situations where the dominant sources of common cause variation span their normal range (as extracted from the analysis of historic data from production runs of related products). The simulations are based on accurate mechanistic/first principle models and reproduce the e ffects of underlying common causes on the 3DSPI measurements of solder paste deposits: area (*a*), height (*h*), volume (*v*), o ffset-X (*x*) and o ffset-Y (*y*). As part of the common cause variation is introduced when di fferent lots of the same product are produced while others are active during the production of a given lot, the AGV module also accounts for both inter-lot, intra-lot and pad specific variability.

During the process simulation, the AGV module generates data for *L* lots with *B* PCBs each, where each PCB is composed by *P* pads (a pad is a small part of the PCB where SPD should be placed). A more comprehensive description of the simulation settings for each of the variability sources are provided in the following subsections. A schematic representation of the AGV module is presented in Figure 2.

**Figure 2.** Schematic representation of the AGV module.

The basic inputs of the AGV module are the nominal values and tolerance limits for the measured parameters (*a*, *h*, *v*, *x*, *y*) and the geometrical coordinates of each pad (Table 1). Based on this information, the allowed variability given the existing tolerance limits is estimated by,

$$\sigma\_{\text{lul}}\text{(parameter, pad)} = \frac{\text{UTL}\{\text{parameter, pad}\} - \text{TTL}\{\text{parameter, pad}\}}{6}.\tag{1}$$

This parameter is analogous to the standard deviation and is based on the fact that, for normally distributed data, 99.73% of the data falls within the ±3 sigma interval.



Along with the PCB's tolerance specifications, the user can also define a set of tuning parameters (Table 2) that allow for the full configuration of the simulated conditions of the printing process. The relationship between the tuning parameters and the process phenomena is discussed in the following subsections. It is expected that, in future applications for the same lines and for related products, the settings for most of these parameters will not change. However, the need for some tuning can be readily assessed after collecting the 3DSPI measurements from the first few products.

**Table 2.** Tuning parameters in the AGV module and the values used in this work.



**Table 2.** *Cont.*

(1), (2), (3), (4) These parameters are between 0 and 1 and the sum of their squares is equal to 1. (5) These values were selected based on engineering knowledge and to achieve a realistic representation of the printing process discussed in Section 5.

#### *3.1. Common Cause Variation A*ff*ecting O*ff*set-X and O*ff*set-Y*

The relative location of the pads can be affected by rigid body translations and rotations of the PCB. Furthermore, offsets in the *y*-coordinate (offset-Y) can also be attributed to the printing direction of the squeegee.

The translation effects are identical for both offset-X and offset-Y. For the sake of brevity, only the case of offset-X is described. To simulate the rigid body translation effect (Figure 3), we considered that it can be split into three components. These are relative to:


In these simulations, *z*(*inter*) *<sup>x</sup>*,*l*,:,: , *z*(*intra*) *<sup>x</sup>*,*l*,*b*,: and *z*(*pad*) *<sup>x</sup>*,*l*,*b*,*p* are generated by taking random draws from the standard normal distribution. The fraction of total variation attributed to each component is determined by the scaling factors α(*inter*) *trans* , α(*intra*) *trans* and α(*pad*) *trans* . These scaling factors must fall between 0 and 1 and are further restricted to comply with:

$$
\left(\alpha\_{trans}^{(intra)}\right)^2 + \left(\alpha\_{trans}^{(intra)}\right)^2 + \left(\alpha\_{trans}^{(pad)}\right)^2 = 1.\tag{2}
$$

Thus, the sum of all components amounts to one unit of allowed variance. Finally, the components are scaled by the allowed variance <sup>σ</sup>2*tolx*, *<sup>p</sup>*ϕ<sup>2</sup>*x*, where <sup>σ</sup>*tolx*, *p* is the allowed standard deviation for the offset-X of the *p*-th pad and ϕ*x* is a scaling factor.

Following the above definitions, the translation effect on the offset in the *x*-coordinate for the *p*-th pad of the *b*-th PCB of the *l*-th lot (*x*(*trans*) *l*,*b*,*p*) is given by (Table 3):

$$x\_{l,b,p}^{(trms)} = \left(\alpha\_{t\text{mus}}^{(inter)} z\_{\text{x},l,:,:}^{(intr)} + \alpha\_{t\text{mus}}^{(intr)} z\_{\text{x},l,b,:}^{(intr)} + \alpha\_{t\text{mus}}^{(pad)} z\_{\text{x},l,b,p}^{(pad)}\right) \mathbf{o}\_{tdl}\{\mathbf{x},p\} \varphi\_{\text{x}}.\tag{3}$$

Likewise, the translation effect on the offsets in the *y*-coordinate (*y*(*trans*) *l*,*b*,*p*) are simulated by (Table 3):

$$\mathbf{x}\_{l,b,p}^{(trans)} = \left(\alpha\_{\text{trans}}^{(inter)} z\_{y,l,r;}^{(intra)} + \alpha\_{\text{trans}}^{(intra)} z\_{y,l,b,r}^{(intra)} + \alpha\_{\text{trans}}^{(pad)} z\_{y,l,b,p}^{(pad)}\right) \mathbf{e}\_{l,l} \{y,p\} \mathbf{q}\_{y}. \tag{4}$$

**Table 3.** Pseudo-code for generating the translation effect on offset-X and offset-Y.

	- 1.1. Generate inter-lot variability: *z*(*inter*) *<sup>x</sup>*,*l*,:,: ∼ *N*(0, <sup>1</sup>), *z*(*inter*) *y*,*l*,:,: ∼ *N*(0, <sup>1</sup>); 1.2. For each PCB, *b* = 1, ... , *B*: 1.2.1. Generate intra-lot variability: *z*(*intra*) *<sup>x</sup>*,*l*,*b*,: ∼ *N*(0, <sup>1</sup>), *z*(*intra*) *y*,*l*,*b*,: ∼ *N*(0, <sup>1</sup>); 1.2.2. For each pad, *p* = 1, ... , *P*: 1.2.2.1.Generate random variability: *z*(*pad*) *<sup>x</sup>*,*l*,*b*,*p* ∼ *N*(0, <sup>1</sup>), *z*(*pad*) *y*,*l*,*b*,*p* ∼ *N*(0, <sup>1</sup>); 1.2.2.2.Compute translation effect on offset-X: *x*(*trans*) *l*,*b*,*p* = α(*inter*) *trans z*(*inter*) *<sup>x</sup>*,*l*,:,: + α(*intra*) *trans z*(*intra*) *<sup>x</sup>*,*l*,*b*,: + α(*pad*) *trans z*(*pad*) *<sup>x</sup>*,*l*,*b*,*<sup>p</sup>*<sup>σ</sup>*tolx*, *p*ϕ*<sup>x</sup>*. 1.2.2.3.Compute translation effect on offset-Y: *y*(*trans*) *l*,*b*,*p* = α(*inter*) *trans z*(*inter*) *y*,*l*,:,: + α(*intra*) *trans z*(*intra*) *y*,*l*,*b*,: + α(*pad*) *trans z*(*pad*) *<sup>y</sup>*,*l*,*b*,*<sup>p</sup>*<sup>σ</sup>*toly*, *<sup>p</sup>*ϕ*y*. >Žƚ >Žƚ W W >Žƚ W

**Figure 3.** Schematic representation of the inter-lot, intra-lot and pad specific components used to generate the translation effect on offset-X.

The rigid body rotation effects (Figure 4), are taking into account by simulating rotations in the PCB with different centers and magnitudes. This rotation can be further decomposed into:


The *z*(*inter*) *l*,:,: and *z*(*intra*) *l*,*b*,: variables are generated from the standard normal distribution and scaled by the allowance rotation angle θ, divided by 3 to ensure that the rotation angle falls, with high probability, between −θ θ . Furthermore, the amount of variability attributed to the inter-lot and intra-lot components is defined by the scale factors α(*inter*) *rot* and α(*intra*) *rot* , which are between 0 and 1 and comply to the condition:

$$
\left(\alpha\_{r\text{rot}}^{(intra)}\right)^2 + \left(\alpha\_{r\text{rot}}^{(intra)}\right)^2 = 1.\tag{5}
$$

Based on this, the rotation angle of the *b*-th PCB of the *l*-th lot (*tl*,*b*,:) is computed as,

$$t\_{l,b,:} = \frac{\alpha\_{rot}^{(inter)} z\_{l,::}^{(inter)} \theta + \alpha\_{rot}^{(intra)} z\_{l,b,:}^{(intra)} \theta}{3}. \tag{6}$$

The rotation matrix is subsequently obtained by,

$$R(t\_{l,b,:}) = \begin{bmatrix} \cos(t\_{l,b,:}) & -\sin(t\_{l,b,:}) \\ \sin(t\_{l,b,:}) & \cos(t\_{l,b,:}) \end{bmatrix}. \tag{7}$$

The center of rotation is randomly generated for each PCB using the uniform distribution. More specifically, the *x*-coordinate for the center of rotation (*rx*) is taken from the uniform distribution on the interval min(**<sup>c</sup>***x*) max(**<sup>c</sup>***x*) , where **c***x* is a vector with the *x*-coordinates of each pad. Similarly, the *y*-coordinate of the center of rotation (*ry*) is generated from the uniform distribution in the interval min(**<sup>c</sup>***y*) max(**<sup>c</sup>***y*) , where **c***y* is a vector with the *y*-coordinates of each pad.

Using the above information, the rotation effects on the offset-X, *x*(*rot*) *l*,*b*,*p* , and offset-Y, *y*(*rot*) *l*,*b*,*p* , for the *p*-th pad of the *b*-th PCB of the *l*-th lot, are determined by (see Table 4.):

$$
\begin{bmatrix} x\_{l,l;p}^{(nt)} & y\_{l,l;p}^{(nt)} \end{bmatrix} = \begin{pmatrix} \begin{bmatrix} \mathbf{c}\_{\mathbf{x}}\langle p \rangle & \mathbf{c}\_{\mathbf{y}}\langle p \rangle \end{bmatrix} - \begin{bmatrix} \begin{array}{cc} r\_{\mathbf{x}} & r\_{\mathbf{y}} \end{bmatrix} \end{pmatrix} \mathbb{R}^{\mathsf{T}} + \begin{bmatrix} \begin{array}{cc} r\_{\mathbf{x}} & r\_{\mathbf{y}} \end{bmatrix} - \begin{bmatrix} \mathbf{c}\_{\mathbf{x}}\langle p \rangle & \mathbf{c}\_{\mathbf{y}}\langle p \rangle \end{bmatrix} \end{bmatrix}. \tag{8}
$$

**Table 4.** Pseudo-code for generating the rotation effect on offset-X and offset-Y.

	- 1.1. Generate inter-lot rotation angle: *z*(*inter*) *l*,:,:∼ *N*(0, <sup>1</sup>);
		- 1.2. For each PCB, *b* = 1, ... , *B*:
			- 1.2.1. Generate intra-lot rotation angle: *z*(*intra*) *l*,*b*,:∼ *N*(0, <sup>1</sup>);
			- 1.2.2. Compute overall rotation angle: *tl*,*b*,:= α(*inter*) *rotz*(*inter*) *l*,:,:θ + α(*intra*) *rotz*(*intra*) *l*,*b*,:θ/3;


 ∼

$$\begin{aligned} r\_{\mathcal{X}} &\sim \mathcal{U}(\min(\mathbf{c}\_{\mathcal{X}}), \max(\mathbf{c}\_{\mathcal{X}})), \\ r\_{\mathcal{Y}} &\sim \mathcal{U}(\min(\mathbf{c}\_{\mathcal{Y}}), \max(\mathbf{c}\_{\mathcal{Y}})). \end{aligned}$$

$$1.2.4.1 \text{ For each pad, } p = 1, \dots, \stackrel{\cdot}{P};$$

 ... 1.2.4.2.Compute rotation effect on offset-X and offset-Y:

,

$$\begin{aligned} \left[ \begin{array}{cccc} \mathbf{x}\_{l,b,p}^{(rot)} & y\_{l,b,p}^{(rot)} & \end{array} \right]\_{\mathbf{c}} &= \left( \begin{bmatrix} \mathbf{c}\_{\mathbf{x}}\{p\} & \mathbf{c}\_{\mathbf{y}}\{p\} \end{bmatrix} - \begin{bmatrix} r\_{\mathbf{x}} & r\_{\mathbf{y}} & \end{bmatrix} \right) \mathbf{R}^{T} + \dots \\ & + \begin{bmatrix} r\_{\mathbf{x}} & r\_{\mathbf{y}} \end{bmatrix} - \begin{bmatrix} \mathbf{c}\_{\mathbf{x}}\{p\} & \mathbf{c}\_{\mathbf{y}}\{p\} \end{bmatrix} \end{aligned}$$

**Figure 4.** Schematic representation of the inter-lot and intra-lot components used to generate the rotation effect on offset-X and offset-Y. The cross represents the specific center of rotation of each PCB.

Along with the translation and rotation effects, the offset-Y is also affected by an additional systematic bias caused by the printing direction of the squeegee (the squeegee is the device that pushes the solder past deposits placed on the top side of the stencil towards the existing apertures in the stencil). This phenomenon leads to a positive deviation in the offset-Y when the printing is performed in one direction and a negative deviation when the printing direction is reversed. Without loss of generality it is assumed that PCBs with an odd index *b* are printed in the direction of increasing *y*-coordinates. Conversely, PCBs with an even index *b* are assumed to be printed in the opposite direction of decreasing magnitude *y*-coordinates. To simulate the squeegee effect on offset-Y (*y*(*sqee*) *l*,*b*,*p* ; Table 5), the allowance deviation that can be caused by the squeegee is defined by <sup>Δ</sup>*y* and its impact is alternated with the PCB index. Furthermore, to account for local random variation, for each PCB, the allowance deviation is multiplied by a random variable with uniform distribution in the interval [0,1].

**Table 5.** Pseudo-code for generating the squeegee effect on offset-Y.


Finally, the offset-X (*x*(*agv*) *l*,*b*,*p* ) and offset-Y (*y*(*agv*) *l*,*b*,*p* ) are determined by summing all effects:

$$\mathbf{x}\_{l,b,p}^{(a\overline{\chi}v)} = \mathbf{N}\{\mathbf{x}, p\} + \mathbf{x}\_{l,b,p}^{(trans)} + \mathbf{x}\_{l,b,p}^{(rot)} \tag{9}$$

$$y\_{l,b,p}^{(agv)} = \mathbf{N}\{y, p\} + y\_{l,b,p}^{(trans)} + y\_{l,b,p}^{(nt)} + y\_{l,b,p}^{(spc)},\tag{10}$$

where **N***x*, *p* and **N***y*, *p* are the nominal value for the offset-X and offset-Y of the *p*-th pad.

#### *3.2. Common Cause Variation A*ff*ecting Height*

The height of solder paste in each pad is affected by the solder mask and the squeegee's action. The effect of the solder mask translates into systematic deviations across all PCBs in the same lot (because the same mask is used in each production lot). The allowed deviation due to the solder mask is here defined by <sup>Δ</sup>*h*,*sold*. Furthermore, its impact on the height is simulated through (Figure 5):


For both cases, *z*(*inter*) *l*,:,: and *z*(*intra*) *l*,*b*,: are generated from the standard normal distribution. The α(*inter*) *h* and α(*intra*) *h* factors scale the contribution of each component. These factors are between 0 and 1 and conditioned to:

$$\left(a\_h^{(intra)}\right)^2 + \left(a\_h^{(intra)}\right)^2 = 1.\tag{11}$$

As the total variability of the height is greater than that caused by the solder mask, the remaining variability is added as a specific component for each pad, *z*(*pad*) *l*,*b*,*p* σ2*tolh*, *<sup>p</sup>*ϕ<sup>2</sup>*h* − <sup>Δ</sup>2*h*,*sold*. Following these considerations, the solder mark effect on the height of the *p*-th pad of the *b*-th PCB of the *l*-th lot (*h*(*sold*) *l*,*b*,*p* ) is computed as follows (see also Table 6):

$$h\_{l,b,p}^{(\text{sold})} = \left(a\_h^{(\text{intra})} z\_{l,:,:}^{(\text{inter})} + a\_h^{(\text{intra})} z\_{l,b,:}^{(\text{intra})}\right) \Delta\_{h,\text{sold}} + z\_{l,b,p}^{(\text{pad})} \sqrt{\mathbf{o}\_{\text{td}}^2 \{h,p\} q\_{\text{ft}}^2 - \Delta\_{h,\text{sold}}^2} \tag{12}$$

**Table 6.** Pseudo-code for generating the solder mask effect on height.

**Figure 5.** Schematic representation of inter-lot, intra-lot and pad specific components used to generate the solder mask effect on height.

As in the case of the offset-Y, the height is also affected by the printing direction of the squeegee. This effect relates to a progressive increase in the pressure of the squeegee. At the beginning of the printing process, the pressure is low and thus the heights are typically lower than intended. As the printing progresses, the heights reach the intended nominal values. The same applies when the printing direction is reversed, but now affecting the pads on the other extreme of the PCB (Figure 6). To be consistent with the squeegee effect on the offset-Y (see Section 3.1), it is considered that PCBs with an odd index *b* are printed in the direction of increasing *y*-coordinates, while PCBs with an even index are printed in the reversed direction.

To systematically address both printing directions, the original *y*-coordinate of the pads is transformed into an equivalent distance, related to the point where the printing beginnings, by:

$$d = (-1)^{b+1} \left( \mathbf{c}\_y \{ p \} - \frac{\max(\mathbf{c}\_y) + \min(\mathbf{c}\_y)}{2} \right) + \frac{\max(\mathbf{c}\_y) - \min(\mathbf{c}\_y)}{2},\tag{13}$$

where **c***y* is a vector with the *y*-coordinates of each pad, *p* is the pad of interest and *b* is the PCB index. Based on this, the influence of the squeegee in the height (*h*(*squee*) *l*,*b*,*p* ) is modeled through an exponential decay function (Table 7):

$$\mathcal{H}\_{l,b,p}^{(\text{square})} = -\Delta\_{b,\text{square}} \mathcal{U}\_{l,b,:} \mathcal{e}^{-\frac{d}{\overline{\tau}}},\tag{14}$$

where <sup>Δ</sup>*h*,*squee* is the maximum allowance deviation, τ is a length constant and *ul*,*b*,: is a random variable following *U*(0, 1). Process knowledge indicates that the deviation in *h*(*squee*) *l*,*b*,*p* is more pronounced in the first third of the PCB (i.e., until *d* = max(**<sup>c</sup>***y*) − min(**<sup>c</sup>***y*) /3). Thus, by considering that for *d* = 2τ the deviation in *h*(*squee*) *l*,*b*,*p* is reduced to 13.5% of the allowance deviation, the length constant is set to:

$$
\pi = \frac{1}{6} (\max(\mathbf{c}\_{\mathcal{Y}}) - \min(\mathbf{c}\_{\mathcal{Y}})).\tag{15}
$$

**Table 7.** Pseudo-code for generating the squeegee e ffect on height.

**Figure 6.** Schematic representation of the squeegee e ffect on height. Dimensions are not to scale.

The final value for the height of each pad (*h* (*agv*) *l*,*b*,*p* ) is determined by summing the solder mask and squeegee e ffects:

$$h\_{l,b,p}^{(agv)} = \mathbf{N}\{h, p\} + h\_{l,b,p}^{(adsd)} + h\_{l,b,p}^{(space)}.\tag{16}$$

#### *3.3. Common Cause Variation A*ff*ecting Area*

The area common cause variation can be conceptually subdivided into three components representing (the structure of each component is analogous to Figure 3):

(i). Inter-lot deviation that impacts all pads of all PCBs in the same lot—inter-lot variation: α(*inter*) *az* (*inter*) *l*,:,:σ*tol* - *a*, *p* ϕ*a*;

(ii). Intra-lot deviation that a ffects all pads in each PCB—intra-lot variation: α(*intra*) *a z* (*intra*) *l*,*b*,: σ*tol* - *a*, *p* ϕ*a*;

(iii). Local variability that a ffects each pad—pad specific variability: α(*pad*) *a z* (*pad*) *l*,*b*,*p* σ*tol* - *a*, *p* ϕ*<sup>a</sup>*.

The random part of each component (*z* (*inter*) *l*,:,: , *z* (*intra*) *l*,*b*,: and *z* (*pad*) *l*,*b*,*p* ) is drawn from the standard normal distribution and scaled by the allowed standard deviation of each pad's area ( σ*tol* - *a*, *p* ) using a scaling factor for the total variance ( ϕ*a*). Furthermore, the contribution of each component is codified by the scaling factors α(*inter*) *a* , α(*intra*) *a* and α(*pad*) *a* , which are between 0 and 1 and are subject to:

$$
\left(\alpha\_a^{(inter)}\right)^2 + \left(\alpha\_a^{(intra)}\right)^2 + \left(\alpha\_a^{(nd)}\right)^2 = 1. \tag{17}
$$

Therefore, the final variance for the area of each pad is <sup>σ</sup>2*tol* - *a*, *p* ϕ2 *a* . Following these definitions, the values for the area (*a* (*agv*) *l*,*b*,*p*) are obtained by (Table 8):

$$a\_{l,b;p}^{(agv)} = \mathbf{N}\{a,p\} + \left(a\_{\mathbf{a}}^{(intra)}z\_{l,:,:}^{(inter)} + a\_{\mathbf{a}}^{(intra)}z\_{l,b;:}^{(intra)} + a\_{\mathbf{a}}^{(pad)}z\_{l,b;p}^{(pad)}\right)\mathbf{o}\_{l:b}\{a,p\}\mathbf{q}\_{\mathbf{f}}.\tag{18}$$

**Table 8.** Pseudo-code for generating the area.


#### *3.4. Common Cause Variation A*ff*ecting Volume*

The volume of solder paste deposit in each pad (*v*(*agv*) *l*,*b*,*p* ) is considered to be proportional to the product of area (*a*(*agv*) *l*,*b*,*p* ) and height (*h*(*agv*) *l*,*b*,*p* ):

$$
\psi\_{l,b,p}^{(a\emptyset v)} = a\_{l,b,p}^{(a\emptyset v)} h\_{l,b,p}^{(a\emptyset v)} \mathbf{f} \{ p \},\tag{19}
$$

where **f***p* is a correction constant (Table 9). For the case studies considered in this article, three different values for **f***p* were observed. However, it was not possible to establish a relationship between this parameter and the pads' location or geometry. Therefore, in the absence of further prior information, it is suggested that **f***p* is estimated using the nominal values for volume, area and height such that:

$$\mathbf{f}\{p\} = \frac{\mathbf{N}\{\upsilon, p\}}{\mathbf{N}\{a, p\}\mathbf{N}\{h, p\}}.\tag{20}$$

**Table 9.** Pseudo-code for generating the volume.

```
1. For each lot, l = 1, ... , L:
       1.1. For each PCB, b = 1, ... , B:
                 1.1.1. For each pad, p = 1, ... , P:
                           1.1.1.1.Compute volume:
                                   v(agv) l,b,p = a(agv) l,b,p h(agv) l,b,p N{v,p} N{a,p}N{h,p}.
```
#### **4. Multivariate Statistical Process Monitoring Based on Principal Component Analysis (MSPM-PCA)**

In the previous section, details were given regarding the artificial generation of common cause variability for the SPD printing process in a Surface Mount Technology (SMT) production line. This is done in the AGV module, and the simulation outputs for area (*a*), height (*h*), volume (*v*), offset-X (*x*) and offset-Y (*y*), are used to augmen<sup>t</sup> the reference NOC dataset, enriching it with a wide overage of common cause variation modes. The augmented NOC dataset can then be used to set up a multivariate statistical process monitoring scheme, which in the present case is based on principal component analysis (PCA).

PCA is a latent variable methodology that decomposes the original data matrix, **X**, with *n* observations and *m* variables, as:

$$\mathbf{X} = \mathbf{T}\mathbf{P}^{\mathrm{T}} + \mathbf{E},\tag{21}$$

where **T** is a (*n* × *k*) matrix of PCA scores, **P** is a (*m* × *k*) matrix with the PCA loadings, **E** is a (*n* × *m*) matrix of residuals and *k* is the number of retained principal components (PC). As PCA is scale-dependent, the data matrix **X** is usually standardized or "autoscaled" to zero mean and unit variance.

The principal components of PCA are monitored by the Hotelling's *T*<sup>2</sup> as proposed by [52,53],

$$T^2 = \sum\_{i=1}^k \frac{t\_i^2}{\lambda\_i} = \mathbf{x}^T \mathbf{P} \mathbf{A}\_k^{-1} \mathbf{P}^T \mathbf{x}\_\prime \tag{22}$$

where **x** is a (*m* × 1) vector with the current observations, **Λ***k* is a diagonal matrix with the first *k* eigenvalues in the main diagonal. If the process follows a multivariate normal distribution, then the upper control limit (UCL) for *T*<sup>2</sup> can be computed as [54,55]:

$$MLCL\_{T^2} = \frac{k(n-1)(n+1)}{n^2 - nk} F\_{a,k,n-k\prime} \tag{23}$$

where *<sup>F</sup>*<sup>α</sup>,*k*,*n*−*<sup>k</sup>* is the upper α-percentile of the *F* distribution, with *k* and *n* − *k* degrees of freedom.

The complementary residuals' subspace is monitored through the *Q*-statistic based on the squared prediction error (SPE) of the residuals, **e** (*m* × 1) [53]:

$$\mathbf{Q} = \mathbf{e}^{\mathsf{T}} \mathbf{e} = \left(\mathbf{x} - \hat{\mathbf{x}}\right)^{\mathsf{T}} \left(\mathbf{x} - \hat{\mathbf{x}}\right) = \mathbf{x}^{\mathsf{T}} \left(\mathbf{I} - \mathbf{P} \mathbf{P}^{\mathsf{T}}\right) \mathbf{x},\tag{24}$$

where **x** (*m* × 1) is the projection of **x** (*m* × 1) onto the PCA subspace. The UCL for this statistic is usually determined by [53,56]:

$$\text{LCL}\_Q = \theta\_1 \left( \frac{z\_a \sqrt{2\theta\_2 h\_0^2}}{\theta\_1} + 1 + \frac{\theta\_2 h\_0 (h\_0 - 1)}{\theta\_1^2} \right)^{1/h\_0},\tag{25}$$

with, **^**

$$\partial\_{i} = \sum\_{j=k+1}^{m} \lambda\_{j\prime}^{i} \text{ i } i = 1 \text{ 2, 3,} \tag{26}$$

$$h\_0 = 1 - \frac{2\theta\_1 \theta\_3}{3\theta\_2^2},\tag{27}$$

where, *z*α is the upper 1 − α percentile of the standard normal distribution.

Due to deviations in the model assumptions, the theoretical control limits for the PCA monitoring statistics are often found to be inaccurate. Therefore, it is often preferable to set the control limits by adjusting a scaled χ2 distribution to the empirical distributions of *T*<sup>2</sup> and *Q* [57–59], leading to:

$$\text{UCL}\_{T^2} = \mathcal{g}\_{T^2} \cdot \chi^2 \Big(\alpha\_\prime h\_{T^2}\Big) \tag{28}$$

$$\text{LCL}\_Q = \text{g}\_Q \cdot \chi^2 \Big(a, h\_Q\big) \tag{29}$$

where *g* is a weighting factor and *h* is the effective number of degrees of freedom for the χ2 distribution. Both values are obtained by matching the moments of the χ2 distribution with those of the empirical distributions obtained from NOC data, leading to *g* = *v*/(<sup>2</sup>*u*) and *h* = <sup>2</sup>*u*2/*<sup>v</sup>*, where *v* is the sample variance and *u* the sample mean of the monitoring statistic.

Under the PCA framework, fault diagnosis can be performed by resort to contribution plots [60–63]. In this context, the contribution of the *i*-th variable for the *T*<sup>2</sup> and *Q* statistics are computed as [62,64],

$$T\_i^2 = \left\| \mathbf{A}^{-1/2} \mathbf{P}\_i^T \mathbf{x}\_i \right\|^2 \tag{30}$$

$$Q\_i = \left\| \mathbf{x}\_i - \mathbf{\hat{x}}\_i \right\|^2 \tag{31}$$

where **P***i* is the *i*-th row of **P** and **x***i* is the *i*-th element of **x**. The theoretical control limits for each *T*2*i* and *Qi* can be found in Refs. [62,64]. Alternatively, the scaled χ2 approximation described before can also be used.

#### **5. Results for the Implementation of AGV in an Industrial SMT Process**

The results presented next regard the printing of solder paste deposits in PCBs composed by 3517 pads. For each pad, measurements for five parameters are collected: (i) area; (ii) height; (iii) volume; (iv) offset-X; (v) and offset-Y. For monitoring purposes, only 3507 pads are considered since the remaining pads are associated with atypical shapes for which the accumulated historic data is not ye<sup>t</sup> enough to set up the AGV module parameters; they are also known to produce unusual offsets in *x*- and *y*-coordinates even in NOC. Consequently, the total number of variables under monitoring is *m* = 3507 *pads* × 5 *parameters* = 17,535 *variables*.

To analyze the characteristics of the proposed AGV module, three real datasets are considered in this study. The first dataset (CS1; Section 5.1) corresponds to NOC data and is composed by 337 PCBs. The second dataset (CS2; Section 5.2) contains 2080 PCBs taken from a different production run. Similarly, the third dataset (CS3; Section 5.3) has 1330 PCBs taken from a third production run. The CS2 and CS3 datasets are mostly composed by normal PCBs, but also contain some faulty PCBs. Furthermore, due to structured common cause variation in the process, the operational state of the three case studies is typically distinct.

In this study, the PCA model is trained using only simulated data generated by the AGV module with the values for the tuning parameters presented in Table 2. No data from the process was employed for building the PCA model or setting up the monitoring charts (they can be also considered, for this reason, pre-control charts). 20 lots with 300 PCBs in each lot were generated in the AGV module. Ten of these lots were used to build the PCA model (training dataset). The other ten lots were used to set the control limits by adjusting the empirical distributions of the monitoring statistics (validation dataset). The control limits were set to a false alarm rate of 0.01. Although we only present here the results for one set of the 20 lots, several replicates were analyzed leading to similar results.

From analysis of the scree plot in Figure 7, the number of retained principal components was set to 5. By further inspection of the PCA loadings (Figure 8) it is verified that the first five components have a close relationship with known phenomena. PC1 describes the relationship between heights and volumes. This PC also shows that the variation in the SPD's height and volume is related to the pads' location in the PCB. PC2 and PC3 explain the rigid body translation effects in the PCBs, while PC4 describes the rotation effects. PC5 is related to the squeegee effect in the heights and its propagation to the volumes. For this reason, PC5 is also related to the printing direction. As for the remaining PCs, they are mostly associated with relationships between areas and volumes of different SPDs. However, since the variation in the area is assumed to be independent from pad to pad, each one of them tend to generates a different PC. Therefore, many PCs are generated by this effect (note that there are over three thousands of pads involved in this mechanism). This justifies why the first 5 PCs explaining the most relevant structured phenomena for process monitoring contain less than 5% of the total variance in the simulated data. However, the monitoring procedure still comprises all the process variation: both the variation captured by the first 5 PCs (through the Hotelling's *T*<sup>2</sup> statistic for the selected PC scores) and the complementary residual variation (using the *Q*-statistic for the PCA residuals). Therefore, this asymmetrical distribution of common cause variability is consistent with the physics of the process and does not pose any problem or limitation for the monitoring scheme.

**Figure 7.** Scree plot of the first 20 eigenvalues of the PCA model.

**Figure 8.** Loadings of the first six principal components. Each block of loadings corresponds to (from left to right) areas, heights, volumes, offset-X and offset-Y.

#### *5.1. Analysis of Dataset CS1: NOC Data*

The CS1 is composed by NOC data, and therefore it will be used to validate the AGV simulation module. By applying the PCA model to this dataset the score plots in Figure 9 are obtained. The results show that the scores from CS1 do not fully overlap the scores from the simulated lots. This is visible in the plot of PC3 vs. PC4 (see Figure 9c), and is an expected behavior since the inter-lot variability causes the lot to fall in a different region of the scores space. In this case, CS1 has lower scores in PC3 than the simulated lots. PC3 is related to translation effects in the PCBs, which implies that the PCBs in CS1 have a systematic translation in a given direction. Note however that some of the simulated lots have translations with similar magnitudes but on the opposite direction. As these inter-lot translations occur randomly and in either direction, we consider that the simulated data is still representative of real operation conditions. The different behavior across the lots is also visible for PC1 and PC4. Since both PCs are related with changes in offset-X and offset-Y, this indicates that the alignment of the PCBs in each lot is a critical factor in the final quality of the products. This information will be further explored in Sections 5.2 and 5.3 to identify possible problems in the different datasets under analysis. As for the intra-lot variability, it is verified that the dispersion of the real data is similar to that obtained with the AGV module.

**Figure 9.** Relation between the principal components of the simulated validation lots and the data Figure 1. (**a**) first and second principal components; (**b**) first and fourth principal components; (**c**) third and fourth principal components; (**d**) fourth and fifth principal components.

As for the control charts produced by the PCA model (Figure 10), it is noticeable that the monitoring statistics are relatively low compared to the control limits. This happens because the PCA model was built using multiple lots, while the monitored data respects to a single lot operating in a well-defined and localized region, spanning only a small fraction of the operational NOC space.

One word of caution should be referred about the proposed approach, resulting from our accumulated experience on its application to real data. Although the proposed methodology decreases the chance of declaring a false alarm due to a change to another lot and solves the issue of false alarms that render classical monitoring methods useless, it also decreases the sensitivity to detect small deviations within a lot since a slightly deviating PCB may be treated as belonging to a different lot rather than being abnormal. This is a consequence of incorporating inter-lot variability in the AGV simulation scenarios.

**Figure 10.** Control charts for the PCBs in CS1: (**a**) *T2*-statisctic; (**b**) *Q*-statistic.

Related to this issue, it was also noted that the *Q*-statistic is affected by the large number of monitored variables (*m*). In this case, there are 17,535 variables, but only 5 PCs are retained in the PCA model. These 5 PCs also explain a small fraction of the total variance (less than 5%), meaning that most of the variability of the data is monitored by the *Q*-statistic. In other words, the *Q*-statistic cumulates a large number of squared residuals with rather large variances. Thus, when a fault occurs, even with a large magnitude, it can suffer from dilution effect from the contribution of all residuals.

From the analysis of the dataset CS1, we can confirm that real NOC data fall inside the expected NOC envelop obtained from the AGV module. The dispersion characteristics of the CS1 lot in the PCA subspace are also similar to the simulated lots, which occupy different regions in the PCA subspace due to the expected (and simulated) inter-lot variability. These observations consubstantiate the rigor placed on the model development process, and provide confidence to the analysis of the remaining datasets using the PCA model developed in this section based on the simulated AGV data.

#### *5.2. Analysis of Dataset CS2: Test Data (NOC and Faulty)*

From analysis of the scores obtained for this case study (Figures 11 and 12) it is apparent that this dataset is composed by at least five distinct operation periods. By comparison with the previously simulated validation lots, it is verified that all periods still fall within the typical operation conditions described by the PCA model applied to AGV data. It is also visible that each period resembles individual lots with specific characteristics. In this case study, a clear distinction in the printing direction is observed for some of the production periods. This phenomenon is more noticeable when PC5 is plotted against PC2 (Figure 13).

By inspecting the control charts for this case study (Figure 14), it is verified that all PCBs are marked as in-control. Nevertheless, it is noticeable that PCBs #159 and #1596 have atypical *Q*-statistics. While for PCB #1596 most of the variables are within their tolerance specifications, PCB #159 present significant deviations in 151 variables. Therefore, PCB #159 is clearly a false negative. The deviations in the measurements of PCB #159 are passed to the residual subspace of PCA. However, since the critical residuals are only a small part of the sum of squared residual over all variables, they are not large enough to trigger an alarm. In spite of this, the contribution plots of PCA are able to identify the abnormal pads and parameters (Figure 15). This result suggests that a monitoring statistic focused solely on the relevant residuals (i.e., those with a significant contribution in the residual subspace) could lead to a more sensitive fault detection.

**Figure 11.** Principal components of CS2: (**a**) first principal component; (**b**) fourth principal component.

**Figure 12.** Relation between the principal components of the simulated validation lots and the data from CS2: (**a**) first and second principal components; (**b**) first and fourth principal components; (**c**) third and fourth principal components; (**d**) fourth and fifth principal components.

**Figure 13.** Second and fifth principal components of CS2 stratified by printing direction: (**a**) PCBs #435 to #834. (**b**) PCBs #835 to #1595.

**Figure 14.** Control charts for the PCBs in CS2: (**a**) *T2*-statisctic; (**b**) *Q*-statistic.

**Figure 15.** Contribution of each pads' parameter to the *Q*-statistic of PCB #159 of CS2: (**a**) contributions for the pads' area; (**b**) contributions for the pads' height.

#### *5.3. Analysis of Dataset CS3: Test Data (NOC and Faulty)*

This dataset comprises at least six distinct operation periods (Figure 16). Once again, these periods are mostly identified by PC1 and PC4, meaning that they are related to lots with di fferent characteristic rotations in the PCBs during the SPD printing process. This result also highlights that the rotation effects on the o ffsets is mostly driven by inter-lot variation. From Figure 17 it is also verified that the selected values for the tuning parameters of the AGV module led to simulated lots resembling real data.

Although the AGV module generates data with realistic features, the PCA-based monitoring methodology presents limitations in its detection capabilities due to the incorporation of inter-lot variability (Figure 18). For instance, from the *Q*-statistic, it is observed that PCBs #248, #401 and #1191 have larger deviations than the others. Nevertheless, these deviations are not large enough to overcome the sum of squared residuals over a large number of variables. By inspecting the contribution plots for these PCBs, it is verified that some of the pads have unusual heights (Figure 19) and therefore these PCBs should have been signaled as potentially faulty. Once again, this issue may be solved in the future if the squared residuals below a given threshold are not included in the *Q*-statistic, which will cause the monitoring scheme to be more sensitive to critical deviations.

**Figure 16.** Principal components of CS3: (**a**) first principal component; (**b**) fourth principal component.

**Figure 17.** Relation between the principal components of the simulated validation lots and the data from CS3: (**a**) first and second principal components; (**b**) first and fourth principal components; (**c**) third and fourth principal components; (**d**) fourth and fifth principal components.

**Figure 18.** Control charts for the PCBs in CS3: (**a**) *T2*-statisctic; (**b**) *Q*-statistic.

**Figure 19.** Contribution of the pads' height to the *Q*-statistic of PCB #248 of CS3.
