Nonlinear Profile Monitoring Using Spline Functions

Abstract

In this study, two new integrated control charts, named $T^2$-MAE chart and MS-MAE chart, are introduced for monitoring the quality of a process when the mathematical form of nonlinear profile model for quality measure is complicated and unable to be specified. The $T^2$-MAE chart is composed of two memoryless-type control charts and the MS-MAE chart is composed of one memory-type and one memoryless-type control charts. The normality assumption of error terms in the nonlinear profile model for both proposed control charts are extended to a generalized model. An intensive simulation study is conducted to evaluate the performance of the $T^2$-MAE and MS-MAE charts. Simulation results show that the MS-MAE chart outperforms the $T^2$-MAE chart with less false alarms during the Phase I monitoring. Moreover, the MS-MAE chart is sensitive to different shifts on the model parameters and profile shape during the Phase II monitoring. An example about the vertical density profile is used for illustration.

1. Introduction

Statistical process control (SPC) is the methodology for quality control using statistical methods for monitoring the quality of a process. Because it can maintain the process to be operated efficiently, produce more confirming products and reduce rework or scrap, SPC has been widely used to monitor the quality of production processes to decline the variability among products and keep their quality close to the desired level. It should also mention that the major tool, control charts, can be treated as a visualization of statistical hypothesis testing. Engineers can use control charts to process data and identify common-cause or assignable-cause variations during the production and distinguish between the process in control and the problematic variation. If the variability is due to an assignable-cause, an out-of-control signal will be flagged before large amount of nonconforming products produced. Therefore, engineers can efficiently maintain the process stable through using SPC. Comprehensive demonstrations about the fundamental to establish control charts and the applications of using control charts can be found in Montgomery [1]. In many instances, the quality of a product can be expressed as a functional form of explanatory variables. The SPC for monitoring functional form data is also named profile monitoring. Recent advances in manufacturing have been proposed, the functional form under monitoring could be complicated and contains a lot of unknown parameters. Multivariate control chart methods are efficient to simultaneously monitor the shift on the model parameters. Two concerns have been raised for using a multivariate control chart to monitor a process. Multivariate control chart methods are usually less sensitive to catch process shifts in Phase II monitoring if the data dimension is 10 or larger. Moreover, parametric multivariate control charts, for example, the Hotelling $T^2$ chart and multivariate exponentially weighted moving average (MEWMA) chart, are sensitive to the violation of normality assumption. If the random vectors of dependent variables do not follow a multivariate normal distribution, the false alarm rate (FAR) is usually inflated and over the nominal level. The FAR in the Phase I monitoring of quality control is also known as the type I error in statistical hypothesis testing methods. The EWMA-type control chart is more sensitive to a small to moderate shift on the process parameter due to using memory-type statistics to establish the control chart.

In profile monitoring, the parameters of interest are often the relationship between the response and explanatory variables and the nature of the variance between and within the profiles. Multivariate control charts are often unable to well monitor the relationship between the response and explanatory variables and the autocorrelation between the observations. In many occasions, more than ten sampling points per profile are measured as data and this fact makes the design for a multivariate control chart method cumbersome. Most of the existing profile monitoring methods assume that the error terms follow a normal distribution to develop the monitoring scheme. How to expand the normality assumption of error terms to a generalized model is also an important issue for profile monitoring.

B-spline approximation is a nonparametric method, which has been used to characterize the functional relation between the response variable and explanatory variables without subjective setting a nonlinear functional form. In the the area of numerical analysis, a B-spline or basis spline is a spline function, which has a minimal support with respect to a given degree, smoothness, and domain partition. Any fixed-degree spline function can be expressed as a linear combination of the B-spline with that degree. The B-spline curve can be used for curve-fitting for a nonlinear data set. The advantage of using a B-spline curve is given as follows: B-spline curve is a generalization of the Bézier curve. The B-spline uses polynomial curve for establishment and requires the information regarding the degree of the curve, a knot vector and a more complex theory than Bézier curves. A B-spline curve can be a Bézier curve. Moreover, the B-spline curve satisfies all important properties that the Bézier curves has. The B-spline curve has more control flexibility than the Bézier curve for modeling. For example, the user can use a lower degree B-spline curve to maintain a large number of control points. The user can change the position of a control point without globally changing the shape of the whole curve based on the local modification property. The B-spline curve can provide a finer shape control due to it satisfies the strong convex hull property. Many techniques for designing and editing the shape of a curve can be done such as changing the knots in a B-spline curve.

Chang and Yadama [2] used the B-spline approximation method to characterize the shape of nonlinear profile and split the nonlinear profile into several sets to reduce the dimension of input variables. After using dimension reduction techniques and doing data transformation, Chang and Yadama [2] suggested using traditional control chart methods for monitoring nonlinear profile processes. Chuang et al. [3] combined the B-spline approximation and control chart methods to develop a nonparametric SPC method for monitoring profile data. Their monitoring method includes five steps: data cleaning, fitting B-spline models, resampling for response data using block bootstrap method, constructing the confidence band based on bootstrap curve depths, and monitoring profiles online based on curve matching. Winistorfer et al. [4] used B-spline approximation and Wavelet transformation methods to construct control charts for monitoring vertical density profile (VDP) data. The VDP is the density distribution through the panel thickness to measure the quality of a wood-based composite panel about its well with strength and physical properties. Woodall et al. [5] discussed some of the general profile monitoring issues via using control chart methods. They also gave a comprehensive review for related SPC literature about profile monitoring. Chang and Chou [6] provided a general discussion about using the Wavelet transformation and B-spline approximation methods for monitoring nonlinear profile models. Hadidoust et al. [7] proposed a Phase II B-spline approximation monitoring method to detect the change-point in a nonlinear profile process.

Fan et al. [8] suggested a Phase II nonlinear profile monitoring procedure via using the transformation of the sum of sine functions. Their method can tackle complicated nonlinear functional form of profiles through heavy numerical computation. Shiau et al. [9] aimed at constructing monitoring schemes for the nonlinear profile model of random effects. They proposed to use principal component method to obtain principal component scores and then use the principal component scores to construct Hotelling $T^2$ chart. Vaghefi et al. [10] proposed two approaches for monitoring nonlinear profile processes: the first approach was to use control chart method with parametric estimates from a regression model as inputs and the second approach was to use metrics to measure the deviation from a reference curve to escape the complexity problems from the coefficient estimation of nonlinear profiles. Williams et al. [11] used Hotelling $T^2$ chart to monitor the coefficients of a parametric nonlinear regression model for a profile process. They used three general formulation methods for the $T^2$ statistics and suggested statistical methods to determine the Phase I control limits. Boullosa-Falces et al. [12] studied the methods to detect small and sudden deviations in the fuel oil process of marine diesel engine when the correlations between input variables were low. They proposed a monitoring method to combine the Hotelling $T^2$ and cumulative sum (CUSUM) charts. The CUSUM chart is also a memory-type control chart, which is developed based on the cumulative sum statistics.

Three concerns in current existing works on nonlinear profile monitoring in literature are stated as: (1) the nonlinear profile model is complicated and difficult to be expressed as an explicit mathematical form in some instances; (2) the existing control chart methods often assume that the error terms in a nonlinear profile model are independent and follow a normal distribution; or (3) the current existing nonlinear profile control chart could not detect process shift when the process shift is due to both the model parameters and profile shape shift together. To enhance the applicability of the control chart methods, expanding the normality assumption of error terms to a generalized distribution family is imperative. Moreover, the proposed nonlinear profile control chart methods are expected more sensitive to alarm the process shift on the model parameters or profile shape during the Phase II monitoring.

B-spline approximation methods are easy for implementation and can be used to fit a lot of nonlinear profile models. In this study, the cubic B-spline curve is used to model nonlinear profile data and then two integrated control charts are proposed for process monitoring. The first integrated chart is to combine the Hotelling $T^2$ chart and the mean absolute error (MAE) chart and named $T^2$-MAE chart. The second integrated chart is to combine the multivariate sign exponentially weighted moving average (MSEWMA) chart and the MAE chart and named MS-MAE chart. Because both Hotelling $T^2$ and MAE charts are memoryless-type control charts and the MSEWMA chart is a memory-type control chart. Hence, the $T^2$-MAE method is composed of two memoryless-type control charts and the MS-MAE chart is composed of one memory-type and one memoryless-type control charts. The Hotelling $T^2$ chart is a parametric control chart and the MSEWMA and MAE charts are nonparametric control charts.

The MSEWMA chart was originally proposed by Zou and Tsung [13]. They have shown the MSEWMA chart has some appealing properties. The MSEWMA statistic is fast to compute and is affine invariant. Moreover, for the distributions with elliptical directions, the charting sequence is a Markov chain process. Then, the in-control average run length (ARL) can be evaluated. The ARL is the expected run length of the interval between out-of-control events. The ARL in Phase I monitoring is a reference measure to design control charts. Compared with the MEWMA chart, the MSEWMA chart is more robust for in-control performance and more sensitive to the small and moderate shifts in location parameters for skewed and heavy-tailed multivariate observations. Zi et al. [14] proposed a distribution-free robust procedure based on using rank-based regression method to reduce the impact from the violation of multivariate normal distribution. Their proposed method is efficient to monitor linear profile processes. In this study, we use Monte Carlo simulations to evaluate the performance of the proposed $T^2$-MAE chart and MS-MAE chart, and then the best control chart method is recommended to monitor nonlinear profile data.

The rest of this manuscript is organized as follows: The construction of the $T^2$-MAE chart and MS-MAE chart via using the cubic B-spline approximation method are studied in Section 2. In Section 3, Monte Carlo simulations are conducted to verify whether the $T^2$-MAE chart and MS-MAE chart can be operated close to the nominal level of FAR during the Phase I monitoring. The performance of the recommended control chart method for out-of-control Phase II monitoring is also evaluated using simulations. The strengths and weaknesses of using the proposed control chart method and the deep learning method of Autoencoder to monitor the quality of a nonlinear profile process are discussed, too. Section 4 demonstrates the application of the recommended control chart method with a re-generated data set for a real example regarding the VDP quality of particleboards. Some concluding remarks are given in Section 5.

2. The B-Spline Approximation and Control Charts

2.1. B-Spline Approximation

To trace the mathematical notations used in the proposed control chart method, a list of mathematical notations is given at the end of the paper for reference. Let Y and X denote the response and explanatory variables, respectively. The realizations of data are denoted by $(x_{ji},y_{ji}),i=1,2,\cdots$,$n_j$ in the jth profile for $j=1,2,\cdots$. The response variable $y_{ji}$ can be expressed as a functional form of $x_{ji}$ through the following model:

$$y_{ji}=S\left(x_{ji}\right)+{ϵ}_{ji},i=1,2,\cdots ,n_j,j=1,2,\cdots ,$$

(1)

where $S(·)$ is a nonlinear function, and ${ϵ}_{ji}$ is an random error. The explicit mathematical function form of $S(·)$ could be complicated and unable to be specified in practical applications. In this study, the error terms are assumed to follow a skew-normal distribution (SND), its probability density function (PDF) can be defined by

$$f(ϵ;{\mu }_{ϵ},\sigma ,\lambda )=\frac{2}{\sigma }\varphi \left(\frac{ϵ-{\mu }_{ϵ}}{\sigma }\right)\Phi \left(\lambda \frac{ϵ-{\mu }_{ϵ}}{\sigma }\right),ϵ\in R,$$

(2)

where ${\mu }_{ϵ}\in R$ and $\sigma >0$ are the location and scale parameters, respectively, and $\lambda$ is the skewness parameter, which controls the skewness of the SND, see Azzalini [15] and Su et al. [16]. The SPC applications via using the SND for Shewhart control chart, economic control chart and control chart for monitoring the quality of linear profile process can be found in Tsai [17], Li et al. [18], Su et al. [19], and Li and Tsai [20]. Denote the SND by SN$({\mu }_{ϵ},\sigma ,\lambda )$ here and after. When $\lambda =0$, SN$({\mu }_{ϵ},\sigma ,0)$ reduces to a normal distribution. The mean and variance of $ϵ$ can be obtained by

$$E\left(ϵ\right)={\mu }_{ϵ}+\sigma \sqrt{\frac{2}{\pi }}\frac{\lambda }{\sqrt{1+{\lambda }^2}}$$

(3)

and

$$Var\left(ϵ\right)=\left(1-\frac{2{\lambda }^2}{\pi (1+{\lambda }^2)}\right){\sigma }^2,$$

(4)

respectively. Let $E\left(ϵ\right)=0$, it can be shown that ${\mu }_{ϵ}=-\sigma \sqrt{\frac{2}{\pi }}\frac{\lambda }{\sqrt{1+{\lambda }^2}}$. The error terms are independent and follow the SN$\left({\mu }_{ϵ},\sigma ,\lambda \right)$.

In practice, an order-d B-spline function can be used to approximate the nonlinear function $S\left(t\right)$. Let ${\mathit{p}}_0,{\mathit{p}}_1,\cdots ,{\mathit{p}}_N$ denote $(N+1)$ control points, $t_0,t_1,\cdots ,t_{N+d}$ be $N+d+1$ non-decreasing knots in the range of t, and $\left\{B_{i,d}\left(t\right),i=0,1,2,\cdots ,N\right\}$ be a set of basis functions of order d, The nonlinear function $S\left(t\right)$ can be represented by

$$S\left(t\right)=\sum _{i=0}^N{\mathit{p}}_i{B}_{i,d}\left(t\right),N\ge (d-1),t\in \left[t_{d-1},t_{N+1}\right],$$

(5)

where

$$B_{i,0}\left(t\right)=\left\{\begin{array}{cc}1,\hfill & t_i<t<t_{i+1}\hfill \\ 0,\hfill & otherwise\hfill \end{array}\right.$$

(6)

$$B_{i,d}\left(t\right)=\frac{t-t_i}{t_{i+d}-t_i}{B}_{i,d-1}\left(t\right)+\frac{t_{i+d+1}-x}{t_{i+d+1}-t_{i+1}}{B}_{i+1,d-1}\left(t\right),$$

(7)

and $\frac{0}{0}=0$ for the coefficients of Equation (7). The basis functions $B_{i,d}\left(x\right)$,$i=0,1,2,\cdots ,N$ have the following properties:

(i)

Partition of unity: $\sum _{i=0}^N{B}_{i,d}\left(t\right)=1$.

(ii)

Positivity: $B_{i,d}\left(t\right)\ge 0$.

(iii)

Local support: $B_{i.d}\left(t\right)=0$ if $x\notin \left[t_i,t_{i+d+1}\right]$.

(iv)

Continuity: the $(d-2)$th derivation of $B_{i,d}\left(t\right)$ is continuous.

Each knot span $t_i\le t\le t_{i+1}$ is mapped onto a polynomial between $S\left(t_i\right)$ and $S\left(t_{i+1}\right)$. Define knots in the interval $[a,b]$. The first and last d knots are set to be $t_0=t_1=\cdots =t_{d-1}=a$ and $t_{N+1}=t_{N+2}=\cdots =t_{N+d}=b$, and the middle $N-d+1$ knots can be chosen to be equidistant. The normalization of the knots resulting in $a=0$ and $b=1$ is helpful to improve numerical accuracy in floating point arithmetic computation. Because a cubic B-spline with $d=4$ is sufficient to approximate most nonlinear curves in realistic applications. In this study, the cubic B-spline method is used to approximate the nonlinear function $S\left(t\right)$. A comprehensive review for using R codes to implement cubic B-spline approximation can be found in de Boor [21] and Perperoglou et al. [22]. Let m in-control profile data sets be collected as Phase I samples. A multivariate control chart for monitoring nonlinear profile processes can be obtained through the following two steps:

Step 1:

Use a cubic B-spline to approximate $S(·)$ in Model (1) and obtain the estimates of spline coefficients and the predicted value of $y_{ji}$, ${\widehat{y}}_{ji}$ for $i=1,2,\cdots ,n_j$, $j=1,2,\cdots ,m$.

Step 2:

The residuals can be obtained by ${\widehat{ϵ}}_{ji}=y_{ji}-{\widehat{y}}_{ji}$, $i=1,2,\cdots ,n_j,j=1,2,\cdots ,m$. Treat the residuals ${\widehat{ϵ}}_i$, $i=1,2,\cdots ,n_j$ from the jth profile as a random sample from $SN\left(-\sigma \sqrt{\frac{2}{\pi }}\frac{\lambda }{\sqrt{1+{\lambda }^2}},\sigma ,\lambda \right)$ and find the maximum likelihood estimates (MLEs) of $\sigma$ and $\lambda$ for each profile and denoted them by ${\widehat{\sigma }}_j$ and ${\widehat{\lambda }}_j$ , respectively, for $j=1,2,,\cdots ,m$. Let ${\mathbf{z}}_j^T=\left({\widehat{\sigma }}_j,{\widehat{\lambda }}_j\right)$, $j=1,2,\cdots ,m$. Then a multivariate control chart can be constructed based on the Phase I sample ${\mathbf{z}}_j$, $j=1,2,\cdots ,m$.

The Hotelling $T^2$ chart and MSEWMA chart are used to develop the proposed monitoring methods. The Hotelling $T^2$ chart is a multivariate extension of the Shewhart-type $\overline{X}$ chart, which is a memoryless-type control chart. The MSEWMA chart is a nonparametric version of the MEWMA chart, which is a memory-type control chart. In general, a memory-type control chart can alarm quicker for out-of-control in Phase II than a memoryless-type control chart, and a parametric control chart is easier to establish than a nonparametric control chart. However, if the normality assumption of the quality variables is violated, then the FAR of the parametric control chart could be inflated and over the nominal level.

2.2. The New Proposed Methods

Let the Phase I sample be ${\mathbf{z}}_j^T=\left({\widehat{\sigma }}_j,{\widehat{\lambda }}_j\right)$, $j=1,2,\cdots ,m$,

$$\overline{\mathbf{z}}=\sum _{j=1}^m{\mathbf{z}}_j/m$$

and

$$S_{\mathbf{z}}=\frac{1}{m}\sum _{j=1}^m\left({\mathbf{z}}_j-\overline{\mathbf{z}}\right){\left({\mathbf{z}}_j-\overline{\mathbf{z}}\right)}^T.$$

The Hotelling $T^2$ statistics can be represented by

$$T_i^2={\left({\mathbf{z}}_i-\overline{\mathbf{z}}\right)}^T{S}_{\mathbf{z}}^{-1}\left({\mathbf{z}}_i-\overline{\mathbf{z}}\right),i=1,2,\cdots .$$

(8)

The control limit in Phase I is

$$L_H=\frac{{(m-1)}^2}{m}{\kappa }_{\alpha ,p/2,(m-p-1)/2},$$

(9)

where ${\kappa }_{\alpha ,p/2,(m-p-1)/2}$ is the upper $\alpha$th quantile of the beta distribution with parameters $p/2$ and $(m-p-1)/2$. The control limit in Phase II is

$$L_H=\frac{p(m+1)(m-1)}{m^2-mp}{F}_{\alpha ,p,m-p},$$

(10)

where $F_{\alpha ,p,m-p}$ is the upper $\alpha$th quantile of the F distribution with degrees of freedom p and $m-p$. If m is large, $L_H\cong {\chi }_{\alpha ,p}^2$, where ${\chi }_{\alpha ,p}^2$ is the upper $\alpha$th quantile of the chi-square distribution with degrees of freedom p.

It is known that a memoryless-type control chart is insensitive to a small or moderate shift on the process parameter. Among all memory-type control charts, the EWMA-type chart is easy to implement and can quickly alarm for small to moderate shifts of the process parameter. The EWMA chart for univariate process is insensitive to the violation of normality assumption; but the MEWMA chart, which is a multivariate version of the EWMA control chart, is very sensitive to the violation of multivariate normal distribution assumption. It is known that the FAR for the MEWMA chart will be inflated and over the nominal level if the random vectors do not follow a multivariate normal distribution.

Because the domain of $\sigma$ is $R^{+}=\left\{\sigma \mid \sigma >0\right\}$, the vector of ${\mathbf{z}}_j^T=\left({\widehat{\sigma }}_j,{\widehat{\lambda }}_j\right)$ does not follow a bivariate normal distribution. In this case, the MSEWMA chart proposed by Zou and Tsung [13] can be considered to implement the process monitoring. The construction of the MSEWMA control chart is addressed as follows: Let ${\delta }_i={\mathbf{z}}_i-\overline{\mathbf{z}},i=1,2,\cdots ,m$. Two unknown parameters, the ${\theta }_0$ and $A_0$, are needed to construct MSEWMA statistics. They can be estimated based on the Phase I sample $\{{\delta }_i,i=1,2,\cdots ,m\}$. Let

$$v_i=\frac{A_0({\delta }_i-{\theta }_0)}{∥A_0({\delta }_i-{\theta }_0)∥},i=1,2,\cdots ,m.$$

(11)

The MSEWMA series can be defined by

$$w_i=(1-\xi )w_{i-1}+\xi v_i,i=1,2,\cdots ,$$

(12)

where $v_0=0$ is the initial value and $0<\xi <1$ is a constant. In general, $\xi =0.1$ and 0.2 are two widely used parameters for implementing the MSEWMA chart. Based on the MSEWMA series in Equation (12), the MSEWMA statistic can be defined by

$$Q_i=\frac{2-\xi }{\xi }p{w}_i^T{w}_i,i=1,2,\cdots$$

(13)

The process is claimed as out-of-control if $Q_i>L_M$. The $L_M$ is a positive control limit. Zou and Tsung [13] have provided tables to report the reference values of $L_M$, those reference values were obtained by simulations.

Let ${\widehat{\theta }}_0$ and ${\widehat{A}}_0$ be the estimates of ${\theta }_0$ and $A_0$, respectively. The values of ${\widehat{\theta }}_0$ and ${\widehat{A}}_0$ can be searched based on the iterations proposed by Hettmansperger and Randles [23].

Step 1:

Obtain the initial value of $\widehat{\theta }$ by

$${\widehat{\theta }}_0^{\left(0\right)}=arg\underset{\theta }{min}\sum _{i=1}^m∥{\delta }_i-\theta ∥.$$

(14)

Step 2:

For iteration $l=0,1,2,\cdots$, search the value of ${\widehat{A}}_0^{\left(l\right)}$ that satisfies Equation (15), given ${\widehat{\theta }}_0^{\left(l\right)}$:

$$\frac{1}{m}\sum _{i=1}^m\left(\frac{A_0^{\left(l\right)}({\delta }_i-{\widehat{\theta }}_0^{\left(l\right)}){({\delta }_i-{\widehat{\theta }}_0^{\left(l\right)})}^T{A}_0^{T\left(l\right)}}{{∥A_0^{\left(l\right)}({\delta }_i-{\widehat{\theta }}_0^{\left(l\right)})∥}^2}\right)=\frac{1}{p}{I}_p.$$

(15)

In practical applications, the value of ${\widehat{A}}_0^{\left(l\right)}$ can be searched based on the Step(a)-Step(e), see Tyler [24]:

Step (a):

Let initial value of $\Omega$ be $\Omega =I_p$.

Step (b):

Evaluate the value of ${\Omega }_x$ by ${\Omega }_x=\left[\frac{p}{trace(\Omega )}\right]\Omega$.

Step (c):

Choose the value of $A_{\Omega }$ such that the condition $A_{\Omega }^T{A}_{\Omega }={\Omega }_x^{-1}$ is true by using Cholesky decomposition method for the positive definite matrix ${\Omega }_x^{-1}$.

Stpe (d):

Update $\Omega$ by

$$\Omega \leftarrow p{\Omega }^{\frac{1}{2}}\frac{1}{m}\sum _{i=1}^m\left(\frac{A_{\Omega }\left({\delta }_i-{\widehat{\theta }}_0^{\left(l\right)}\right)}{∥A_{\Omega }\left({\delta }_i-{\widehat{\theta }}_0^{\left(l\right)}\right)∥}\right){\left(\frac{A_{\Omega }\left({\delta }_i-{\widehat{\theta }}_0^{\left(l\right)}\right)}{∥A_{\Omega }\left({\delta }_i-{\widehat{\theta }}_0^{\left(l\right)}\right)∥}\right)}^T{\Omega }^{\frac{1}{2}}.$$

(16)

Stpe (e):

Repeat the Steps (b) to (d) until convergence.

Step 3:

Update ${\widehat{\theta }}_0$ by ${\theta }_0^{(l+1)}={\left[A_0^{\left(l\right)}\right]}^{-1}{\theta }_{\delta }$,

$${\theta }_{\delta }=arg\underset{\theta }{min}\sum _{i=1}^m∥D_i-\theta ∥,$$

(17)

where $D_i=A_0^{\left(l\right)}{\delta }_i$, $i=1,2,\cdots ,m$.

Step 4:

Let $l=l+1$. Repeat Step 2 to Step 3 until

$$\frac{∥{\widehat{\theta }}_0^{\left(l\right)}-{\widehat{\theta }}_0^{(}l-1)∥}{∥{\widehat{\theta }}_0^{(l-1)}∥}\le e,$$

(18)

where e is a predetermined small positive number. The vales of ${\widehat{\theta }}_0^{\left(l\right)}$ and ${\widehat{A}}_0^{\left(l\right)}$ at the final iteration when the algorithm is stopped can be the ${\widehat{\theta }}_0$ and ${\widehat{A}}_0$, respectively.

The MSEWMA chart and Hotelling $T^2$ chart are insensitive to the shift on the profile shape. Hence, we suggest to integrate the MSEWMA chart (or Hotelling $T^2$ chart) with the MAE chart, and name the new integrated chart by MS-MAE chart (or $T^2$-MAE chart). The MAE for profile j is defined by

$${\widehat{\eta }}_j=\frac{{\sum }_{i=1}^n\left|{\widehat{ϵ}}_{ji}\right|}{n},j=1,2,\cdots ,m.$$

(19)

Let the upper control limit of the MAE chart be $q_{\alpha }$, which is the upper $\alpha$th quantile such that $P(\widehat{\eta }>q_{\alpha })=\alpha$. In practice, we can use the empirical distribution of $\left\{{\widehat{\eta }}_j,j=1,2,\cdots ,m\right\}$ to obtain an approximated value of $q_{\alpha }$ when m is large. Let the FARs of the MSEWMA (or Hotelling $T^2$) and MAE charts be ${\alpha }_1$ and ${\alpha }_2$, respectively. Then, the overall FAR can be $\alpha =(1-{\alpha }_1)(1-{\alpha }_2)$; that is, we can select the values of ${\alpha }_1$ and ${\alpha }_2$ to maintain the overall FAR at the desired level. For example, if we can take ${\alpha }_1={\alpha }_2=0.0027$ to construct the proposed MS-MAE chart or $T^2$-MAE chart, the overall FAR of the integrated chart is $\alpha =0.0054$. The corresponding ARL for the in-control Phase I process is about 185. The two integrated control charts can be implemented through using the following Step 1 to Step 4:

Step 1:

Obtain the estimates of $\{{\mathbf{z}}_j^T=({\widehat{\sigma }}_j,{\widehat{\lambda }}_j),j=1,2,...,m\}$ and the residuals $\{{\widehat{ϵ}}_{ij},i=1,2,...,n,j=1,2,...,m\}$ based on m in-control Phase I subgroups by using the proposed estimation method in Section 2.1.

Step 2:

Obtain the Hotelling $T^2$ series based on Equation (8) and obtain the mean residuals, $\{{\widehat{\eta }}_j,j=1,2,...,m\}$, based on Equation (19).

Step 3:

The first control chart in the integrated control chart can be the Hotelling $T^2$ chart or MSEWMA chart. The second control chart in the integrated control chart is the MAE chart. Let the FAR in the first and second control charts are ${\alpha }_1$ and ${\alpha }_2$, respectively, and the overall FAR is $\alpha =(1-{\alpha }_1)(1-{\alpha }_2)$.

Step 3.1:

The control limits in Equations (9) and (10) with $FAR={\alpha }_1$ can be used to establish the Hotelling $T^2$.

Step 3.2:

The control limit of the MSEWMA chart with $FAR={\alpha }_1$ can be selected from the tables proposed by Zou and Tsung [13] or we can search the control limit via using simulation methods.

Step 3.3:

The quantile, $q_{{\alpha }_2}$, based on the empirical distribution of $\widehat{\eta }$ can be obtained as the control limit of MAE chart.

Step 4:

The $T^2$-MAE chart is to integrate the Hotelling $T^2$ and MAE charts in Step 3 and the MS-MAE chart is to integrate the MSEWMA and MAE charts in Step 3.

3. Monte Carlo Simulations and Discussions

3.1. Performance Evaluation

The R language codes with the packages of “sn”, “splines” and “optim” are used to generate random samples from the SND, implement the B-spline approximation and obtain the estimates of the model parameters, respectively. After the proposed two integrated control charts are established, a simulation study is conducted with has two goals. The Goal 1 is to study whether the proposed integrated control charts of the $T^2$-MAE chart and MS-MAE chart can be maintained at the nominal FAR level for the in-control Phase I monitoring. The Goal 2 is to study the speed to alarm of the recommended integrated control chart for parameter shifts in Phase II monitoring when the process is out of control. Nonlinear profile data are simulated using Equation (20), which was proposed in Noorossana et al. [25] for the VDP data:

$$S\left(x_{ji}\right)=\left\{\begin{array}{cc}{\theta }_1{(x_{ji}-c)}^{{\theta }_3}+{\theta }_6,\hfill & x_{ji}>{\theta }_5,\hfill \\ {\theta }_2{(-x_{ji}+c)}^{{\theta }_4}+{\theta }_6,\hfill & x_{ji}\le {\theta }_5,\hfill \end{array}\right.$$

(20)

for $i=1,2,\cdots ,n$ and $j=1,2,\cdots ,m$. The total 24 sets of $({\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4,{\theta }_5,{\theta }_6)$ given by Noorossana et al. [25] are used to generate nonlinear profiles, we take the median of each set of ${\theta }_h$ for $h=1,2,\cdots ,6$ to generate the reference nonlinear profile for simulations in this section, and the parameters used to generate the reference nonlinear profile are ${\theta }_1=3506$, ${\theta }_2=2561$, ${\theta }_3=4.855$, ${\theta }_4=4.125$, ${\theta }_5=0.3$ and ${\theta }_6=44.455$. Following the suggestions of Noorossana et al. [25], Equations (1) and (2) with $x=0,0.02,0.04,0.006,\cdots ,0.626$, $\sigma =2$ and $\lambda =3$ are used to generate total $n_j=314$ pairs of $(x_{ij},y_{ij})$ for profile j. Each generated data set can be used to obtain one vector of ${\mathbf{z}}^T=(\widehat{\sigma },\widehat{\lambda })$ and one $\widehat{\eta }$. All ${\mathbf{z}}_j$ and ${\widehat{\eta }}_j$, $j=1,\cdots ,m$, are used to construct the $T^2$-MAE and MS-MAE charts with the FARs ${\alpha }_1={\alpha }_2=0.0027$. Hence, the overall FAR for two integrated control charts are $\alpha =0.0054$. The nominal in-control Phase I ARL is about 185, denoted by ARL${}_0$ = 185.

Figure 1 reports the simulated values of ARL${}_0$ of the MS-MAE chart with $\xi$ = 0.2 and the $T^2$-MAE chart. Figure 1 shows that the simulated ARL${}_0$ of the MS-MAE chart is very closed to the nominal value when the number of samples $m\ge 800$. However, the simulated ARL${}_0$ of the $T^2$-MAE chart always seriously underestimates the nominal ARL${}_0$ even the number of samples increases to 2000. It is found that the Hotelling $T^2$ chart has more false alarms than the MAE chart. Because the input vectors of $z_j$ do not follow a bivariate normal distribution, the $T^2$-MAE chart causes more false alarms for monitoring the parameter estimates than the expectation. Hence, the $T^2$-MAE chart is not a suitable tool for monitoring nonlinear profile data in this study. We recommend using the MS-MAE chart to monitor nonlinear profile processes. Hence, only the ability of the MS-MAE chart for detecting process shifts in Phase II is studied in the following simulations.

For generating out-of-control profile data to evaluate the performance of the MS-MAE chart to capture the shift of the curve at different profile locations, 10,000 data sets, $(x_{ji},y_{ji}),$$i=1,2,\cdots ,1000,$$j=1,2,\cdots ,$ 10,000, are generated based on Equation (20). Then all generated data set are fitted by Model (21) via using the package splines of the software R, where ${\beta }_0$, ${\beta }_1$, ⋯, ${\beta }_7$ are parameters and $B_h\left(x\right)$, $h=1,2,\cdots ,7$ are basis functions in the cubic B-spline approximation method.

$$Y_{ji}={\beta }_0+{\beta }_1{B}_1\left(x\right)+{\beta }_2{B}_2\left(x\right)+\cdots +{\beta }_7{B}_7\left(x\right)+ϵ$$

(21)

Let ${\mathit{\beta }}^T=({\beta }_0,{\beta }_1,\cdots ,{\beta }_7)$ and the estimate of $\mathit{\beta }$ based on the jth data set $\left\{(x_{ji},y_{ji}),i=1,2,\cdots ,1000\right\}$ be denoted by ${\widehat{\mathit{\beta }}}_j$, $j=1,2,\cdots ,$ 10,000. In the following simulation study, we use the mean vector, ${\gamma }^T=\frac{1}{10,000}{\sum }_{j=1}^{10,000}{\widehat{{\mathit{\beta }}_{\mathbf{j}}}}^T=({\gamma }_0,{\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4,{\gamma }_5,{\gamma }_6,{\gamma }_7)$, as the in-control vector of nominal parameters for the nonlinear profile model. Ten scenarios of profile shifts, labeled by Shift 1 to Shift 10 and given as follows, are considered for Phase II monitoring:

Shift 1:

${\gamma }^{\left(1\right)}=({\gamma }_0,0.98\times {\gamma }_1,0.98\times {\gamma }_2,{\gamma }_3,{\gamma }_4,{\gamma }_5,{\gamma }_6,{\gamma }_7)$

Shift 2:

${\gamma }^{\left(2\right)}=({\gamma }_0,{\gamma }_1,0.98\times {\gamma }_2,0.98\times {\gamma }_3,{\gamma }_4,{\gamma }_5,{\gamma }_6,{\gamma }_7)$

Shift 3:

${\gamma }^{\left(3\right)}=({\gamma }_0,{\gamma }_1,{\gamma }_2,0.98\times {\gamma }_3,0.98\times {\gamma }_4,{\gamma }_5,{\gamma }_6,{\gamma }_7)$

Shift 4:

${\gamma }^{\left(4\right)}=({\gamma }_0,{\gamma }_1,{\gamma }_2,{\gamma }_3,0.98\times {\gamma }_4,0.98\times {\gamma }_5,{\gamma }_6,{\gamma }_7)$

Shift 5:

${\gamma }^{\left(5\right)}=({\gamma }_0,{\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4,0.98\times {\gamma }_5,0.98\times {\gamma }_6,{\gamma }_7)$

Shift 6:

${\gamma }^{\left(6\right)}=({\gamma }_0,1.02\times {\gamma }_1,1.02\times {\gamma }_2,{\gamma }_3,{\gamma }_4,{\gamma }_5,{\gamma }_6,{\gamma }_7)$

Shift 7:

${\gamma }^{\left(7\right)}=({\gamma }_0,{\gamma }_1,1.02\times {\gamma }_2,1.02\times {\gamma }_3,{\gamma }_4,{\gamma }_5,{\gamma }_6,{\gamma }_7)$

Shift 8:

${\gamma }^{\left(8\right)}=({\gamma }_0,{\gamma }_1,{\gamma }_2,1.02\times {\gamma }_3,1.02\times {\gamma }_4,{\gamma }_5,{\gamma }_6,{\gamma }_7)$

Shift 9:

${\gamma }^{\left(9\right)}=({\gamma }_0,{\gamma }_1,{\gamma }_2,{\gamma }_3,1.02\times {\gamma }_4,1.02\times {\gamma }_5,{\gamma }_6,{\gamma }_7)$

Shift 10:

${\gamma }^{\left(10\right)}=({\gamma }_0,{\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4,1.02\times {\gamma }_5,1.02\times {\gamma }_6,{\gamma }_7)$

Figure 2 shows the positions of profile shape shift. Shift 1 to Shift 5 indicate that the curve shifts to its inside at different locations of the reference profile, and Shift 6 to Shift 10 indicate that the curve shifts to its outside at different locations of the reference profile. Shift 1 and Shift 6 indicate that the curve shifts at the upper-side arm of the reference profile; Shift 2 and Shift 7 indicate that the curve shifts at the middle bottom of the reference profile; Shift 3 and Shift 9 indicate that the curve shifts at the right bottom the reference profile; and Shift 4 and Shift 10 indicate that the curve shifts at the left bottom of the reference profile.

The control limits of the MS-MAE chart is obtained with the FARs ${\alpha }_1$ = ${\alpha }_2=0.0027$ and the overall FAR can be maintained at $\alpha =0.0054$. The ARL in Phase II to alarm an out-of-control case is labeled as ARL${}_1$, which is evaluated based on 10,000 runs of simulation for each of the combinations of the SND parameters $(\sigma ,\lambda )=$(1,1),(1,3), (3,1), (3,3), (5,1) and (5,3), respectively. The smaller the ARL${}_1$ is, the better of the ability of the control chart to alarm an out-of-control signal. Moreover, the corresponding standard deviation of run length (SDRL) for each shift, from Shift 1 to Shift 10, is respectively evaluated based on 10,000 simulation runs, too. All simulation results are reported in

Table 1. The ARL${}_1$ and SDRL of the MS-MAE chart for different combinations of $\sigma$, $\lambda$ and $\xi =0.1$.

		$(\mathit{\sigma },\mathit{\lambda })$
		(1,1)	(1,3)	(3,1)	(3,3)	(5,1)	(5,3)
Shift 1	ARL1	3.8412	2.3802	7.7024	13.1864	12.2773	40.5034
	SDRL	2.2930	0.7246	1.3110	4.3776	4.9731	30.9255
Shift 2	ARL1	2.0835	2.0009	7.4032	10.7722	11.7010	29.1175
	SDRL	0.3009	0.0300	1.6210	4.9118	5.0599	23.0682
Shift 3	ARL1	2.0289	2.0005	7.2282	9.5635	11.4067	26.2756
	SDRL	0.1757	0.0224	1.7408	4.8439	5.0496	20.6749
Shift 4	ARL1	2.0840	2.0017	7.4103	10.6689	11.5696	29.1012
	SDRL	0.3046	0.0412	1.6292	4.8818	5.0152	23.1551
Shift 5	ARL1	3.3459	2.2565	7.6737	13.0108	12.2735	39.3576
	SDRL	1.7959	0.5686	1.3659	4.4637	5.0878	30.4337
Shift 6	ARL1	3.8125	2.3554	7.6941	13.2833	12.3074	40.7550
	SDRL	2.2551	0.6813	1.2853	4.3092	4.9909	31.5401
Shift 7	ARL1	2.0853	2.0018	7.3737	10.7119	11.7145	30.4102
	SDRL	0.3020	0.0424	1.6213	4.8313	5.1583	23.8104
Shift 8	ARL1	2.0331	2.0003	7.2480	9.5152	11.4173	27.1037
	SDRL	0.1866	0.0173	1.7236	4.7609	5.0859	21.2876
Shift 9	ARL1	2.0797	2.0023	7.3797	10.7447	11.6904	30.7924
	SDRL	0.2939	0.0479	1.6091	4.8852	5.0775	23.7903
Shift 10	ARL1	3.3841	2.2383	7.6936	13.2156	12.2797	40.9350
	SDRL	1.7945	0.5373	1.3583	4.4369	5.0288	31.8419

for $\xi =0.1$ and in

Table 2. The ARL${}_1$ and SDRL of the MS-MAE chart for different combinations of $\sigma$, $\lambda$ and $\xi =0.2$.

		$(\mathit{\sigma },\mathit{\lambda })$
		(1,1)	(1,3)	(3,1)	(3,3)	(5,1)	(5,3)
Shift 1	ARL1	4.1565	2.6312	6.4459	11.9355	15.2565	44.0534
	SDRL	2.6307	1.0141	1.4416	5.7643	8.7400	38.9507
Shift 2	ARL1	2.1123	2.0045	6.0158	8.5698	14.5249	28.9411
	SDRL	0.3520	0.0669	1.7091	4.8396	8.3630	25.5944
Shift 3	ARL1	2.0360	2.0004	5.7817	7.3016	14.2728	24.3020
	SDRL	0.1932	0.0200	1.8407	4.4238	8.4304	21.2048
Shift 4	ARL1	2.1066	2.0047	6.0167	8.5518	14.5316	28.0087
	SDRL	0.3462	0.0684	1.7050	4.8475	8.2793	24.8782
Shift 5	ARL1	3.7319	2.4397	6.4609	11.7575	15.1310	43.0275
	SDRL	2.1750	0.7979	1.4777	5.7632	8.8302	38.6364
Shift 6	ARL1	4.2259	2.6179	6.4738	12.0545	15.2113	45.7271
	SDRL	2.6815	1.0208	1.4308	5.8186	8.6970	40.0272
Shift 7	ARL1	2.1157	2.0058	6.0187	8.5165	14.5941	29.7988
	SDRL	0.3607	0.0759	1.7339	4.8585	8.3315	26.3605
Shift 8	ARL1	2.0373	2.0006	5.7978	7.2286	14.3202	25.5466
	SDRL	0.1957	0.0245	1.8242	4.2111	8.1342	22.2157
Shift 9	ARL1	2.1107	2.0046	6.0441	8.5784	14.7236	29.8490
	SDRL	0.3476	0.0677	1.7032	4.7884	8.5263	26.1232
Shift 10	ARL1	3.7231	2.4397	6.4618	11.7864	15.1848	43.5996
	SDRL	2.1656	0.7994	1.4660	5.6724	8.6189	39.0241

for $\xi =0.2$.

In view of Table 1 and Table 2, it can be noticed that the proposed MS-MAE chart is sensitive to alarm the shift at different locations of the reference curve with a short ARL${}_1$. When the scale parameter of the SND is large or the SND is highly skewed, the ability of the MS-MAE chart to alarm out-of-control is reduced. When $\sigma$ or $\lambda$ in the SND become large, the MS-MAE chart needs more time or a longer ARL${}_1$ to alarm out of control in Phase II. Moreover, Table 1 and Table 2 also show that the MS-MAE chart with $\xi$ = 0.1 and MS-MAE chart with $\xi$ = 0.2 are competitive. Both charts have two close values of ARL${}_1$ to alarm process shifts for the same parameter combination of $\sigma$ and $\lambda$. The SDRL in Table 1 and Table 2 are also close for the same parameter combination of $\sigma$ and $\lambda$. Overall, the MS-MAE chart with $\xi =0.1$ can alarm out of control a little quicker than the MS-MAE chart with $\xi =0.2$.

3.2. Discussions

Today, deep learning methods have earned more attention in machine learning applications. Deep learning refers to multi-layer neural networks, which can learn complex patterns for pattern recognition. If we treat the control chart method as an issue of pattern recognition, the deep leaning methods are competitive with the proposed control chart method. The Autoencoder is one popular deep learning method in an unsupervised manner. The Autoencoder can be used to learn a representation of encoding for a data set via training the network to ignore the signal of noise for dimensionality reduction. A reconstructing side is learnt from the reduction side, in which the Autoencoder tries to copy its input to its output to maintain the representation as close as possible to its input. An internal or hidden layer in an Autoencoder can be used to describe a code for representing the input, and it is composed of two main parts: an encoder and decoder for mapping the input into the code and mapping the code into a reconstruction of the input, respectively. The idea of Autoencoder has been widely applied in neural networks for deep learning in the past few decades.

As other deep learning algorithms, the Autoencoder requires large amounts of data for training via intensive computation. Moreover, the Autoencoder requires much more expertise to set the architecture and hyper-parameters. Hence, the Autoencoder could be unsuitable as a general-purpose algorithm. In many occasions of quality control, it could be difficult to offer large amounts of data for training model in a deep learning algorithm due to the considerations of sampling cost, destructive testing or other reasons. The proposed control chart method uses a cubic B-spline curve to characterize nonlinear profile data and monitor the quality of a process via using an integrated control chart. The proposed control chart method could be less effective than the Autoencoder algorithm to process a broad nonlinear functions. However, the proposed control chart method requires less sample resource to obtain reliable monitoring results than using a machine learning method. As we have mentioned in the literature review in Section 1, the cubic B-spline function has been proved to well approximate a lot of nonlinear functions in quality control applications. When the sample resource is limited for training model, the proposed control chart method is competitive for monitoring the quality of a nonlinear profile process. If large amounts of data can be offered by a nonlinear profile process for training model, the Autoencoder algorithm can be considered to replace the proposed control chart method for process monitoring. By the way, how to conduct a diagnostic process for identifying the causes of out of control via using a deep learning algorithm for process monitoring is still an open question in SPC.

4. An Application

In this section, an example regarding monitoring the quality of particle boards is used for illustration. The VDP through the panel thickness has been identified as an important panel characteristic, which is related to the strength and physical properties of particleboards. A profilometer is used to monitor the density of the finished particle boards over time at fixed depths across the thickness of the particle board during the manufacturing of particle boards. Density measurements for the VDP data were taken at depths $x_i=0.002\times i$, $i=0,1,\cdots ,313$ to give a sample of $(x_{ji},y_{ji})$, $i=1,2,\cdots ,n$, $j=1,2,\cdots ,m$.

Because we do not have the original data set, the process proposed by Noorossana et al. [25] is followed to generate VDP data sets for illustration. Equation (21) with the obtained $\gamma$ in Section 3, based on in-control Phase I samples with $\sigma$ = 1 and $\lambda$ = 3, is used to generate 50 in-control Phase I nonlinear profiles, and we use ${\gamma }^{*}=({\gamma }_0,1.05\times {\gamma }_1,1.02\times {\gamma }_2,{\gamma }_3,{\gamma }_4,{\gamma }_5,{\gamma }_6,{\gamma }_7)$ to generate additional 25 out-of-control profiles; that is, the SND parameters of $\sigma$ and $\lambda$ do not shift, but the profile shape shifts at its right side as that in Figure 2g. Using the FARs ${\alpha }_1={\alpha }_2=0.0027$ and $\alpha =0.0054$ for the in-control Phase I monitoring, the numerical control limits $L_M$ = 8.567 and $L_A=0.144$ can be obtained via using simulations for the MS-MAE chart with $\xi =0.2$. The MS-MAE chart for this data set is given in Figure 3. The MSEWMA control chart in the MS-MAE method is used to monitor the parameters $\sigma$ and $\lambda$, and the MAE control chart is used to monitor the shift on the profile shape.The MS-MAE chart triggers a signal for out of control at profile 52, and then several out-of-control signals follow. In scan of Figure 3, we find that the MSEWMA chart does not alarm for out of control but the MAE chart triggers out-of-control signals at profile 52 and after. Therefore, the MS-MAE chart can quickly alarm out-of-control in this example.

5. Concluding Remarks

In this study, two new integrated control charts, the $T^2$-MAE chart and MS-MAE chart, are proposed to monitor nonlinear profile data when the mathematical form of the nonlinear profile model is complicated and unable to be specified. Moreover, the error terms in the nonlinear profile model are assumed to follow a SND, which is a generalized version of the normal distribution. The cubic B-spline approximation method is used for modeling the profile data, and the residuals are used to obtain the MLEs of the SND parameters.

Monte Carlo simulations are conducted to verify the performance of the proposed $T^2$-MAE chart and MS-MAE chart. The simulation results have shown that the $T^2$-MAE chart was unable to perform well due to the input vector of the MLEs of the SND parameters does not follow a bivariate normal distribution. The FAR of the $T^2$-MAE chart highly overestimates the nominal value of the FAR in the Phase I monitoring. The MS-MAE chart is more reliable than the $T^2$-MAE chart. Moreover, the MS-MAE chart is sensitive to alarm the process shifts on the distribution parameters or profile shape. The re-generated VDP data sets of particle boards have been used to show that the proposed MS-MAE chart can quickly alarm out of control in Phase II monitoring. If large amounts of data can be offered by a production process for training model, the deep learning method via using the Autoencoder algorithm can be considered to replace the proposed control chart method for process monitoring.

This study focuses on providing a general control chart method to monitor the quality of a nonlinear profile process. The VDP data is used for illustration. However, the proposed method can be used to monitor any nonlinear profile process with data that can be characterized via using a cubic B-spline approximation. In view of the simulation results, the parameter estimation for the SND was found not stable if the sample size in each profile is small. How to obtain reliable parameter estimation for the SND via using other algorithms is an open question. It also is important to reduce the sample resource to establish the MSEWMA chart in the proposed control chart method. Some machine learning methods are also competitive for the process monitoring purpose. However, the sample resource limitation is a problem for many production processes. How to conduct a efficient diagnostic process for identifying the out-of-control causes via using a deep learning algorithm for process monitoring is still an open question in SPC. It is also important to design a new control chart based on the economical purposes or conditions for monitoring the quality of a process. These topics will be studied in the future.