An Improved Hidden Markov Model for Monitoring the Process with Autocorrelated Observations

Abstract

With the development of intelligent manufacturing, automated data acquisition techniques are widely used. The autocorrelations between data that are collected from production processes have become more common. Residual charts are a good approach to monitoring the process with data autocorrelation. An improved hidden Markov model (IHMM) for the prediction of autocorrelated observations and a new expectation maximization (EM) algorithm is proposed. A residual chart based on IHMM is employed to monitor the autocorrelated process. The numerical experiment shows that, in general, IHMMs outperform both conventional hidden Markov models (HMMs) and autoregressive (AR) models in quality shift diagnosis, decreasing the cost of missing alarms. Moreover, the times taken by IHMMs for training and prediction are found to be much less than those of HMMs.

1. Introduction

Process monitoring plays an essential role in intelligent manufacturing [1,2]. Statistical process control (SPC) is a quality control technique that uses statistical methods to monitor and control processes. The aim of SPC is to ensure that the process runs efficiently, producing more products that meet specifications, while reducing waste at the same time. Shewhart control charts are a key SPC tool used to determine whether a process is in control. If the observation value is within the upper and lower control limits, the process is in control; otherwise, the process is out of control. Shewhart charts often assume that the data collected from the process are independent. However, this assumption is not true in a variety of processes. For example, consecutive measurements on chemical and pharmaceutical processes or product characteristics are often highly autocorrelated, or the chronological measurements of every characteristic on every unit in automated test and inspection procedures are often autocorrelated. It has been shown in numerous studies that conventional charts do not work well, in the form of giving too many false alarms if the observations exhibit even a small amount of autocorrelation over time [3,4,5,6,7,8,9]. Clearly, better approaches are needed. In the following paragraphs, three common ways to solve the problems related to autocorrelation are introduced.

The first approach to solving an autocorrelation problem is simply to sample from the observation data stream at a lower frequency [10]. This seems to be an easy solution, although it has some shortcomings. Concretely speaking, it makes inefficient use of the available data. For example, if every tenth observation is sampled, approximately 90% of the information is discarded. In addition, since only every tenth datum is used, this approach may delay the discovery of a process shift.

The second approach is to re-estimate the real process variance, aiming to revise the upper and lower control limits. See, for example, [9,11,12,13,14,15,16,17,18,19].

The third approach is to use residual charts. The statistics in these charts are residuals; values are calculated by subtracting predicted values from observed values. In the implementation of residual charts, the key step is choosing a reasonable prediction model to obtain predicted values. Autoregressive integrated moving average (ARIMA) models are mostly used to model autocorrelated data. See, for example, [3,8,20,21,22,23,24,25,26]. In addition, multistage or multivariate autocorrelated processes are mostly studied by using Hotelling T² control charts, such as in [27,28,29,30]. Some exceptions include Pan et al., who proposed an overall run length (ORL) to replace T² charts [31], and S. Yang and C. Yang, who used a residual chart and cause-selecting control chart [32].

With the rapid development of artificial intelligence, machine-learning methods are becoming more and more popular in SPC in the case of observation autocorrelation [33,34,35,36,37,38]. Most of the related literature focuses on neural networks methods [39,40,41,42]. Some successful applications have also been introduced, for example, the failure diagnosis of wind turbine [43], and the prediction of the remaining useful life (RUL) of bearings [44]. Another machine-learning approach proposed in SPC is the hidden Markov model (HMM). Lee et al. proposed a modified HMM, combined with a Hotelling multivariate control chart to perform adaptive fault diagnosis [45]. The HMMs, whose training sets contain autocorrelated data, were employed to forecast observation values for residual charts in process monitoring [46]. However, although HMMs are supposed to deal with the autocorrelated processes, the essence of the model itself is inconsistent with the case of autocorrelations. This is because one key assumption in conventional HMMs is that the observations are independent of each other. Therefore, it is worth developing a modified HMM, by accounting for observation autocorrelations in models themselves, and observing whether it is better than a traditional HMM.

Therefore, to realize the goal of monitoring the autocorrelated process well, the residual chart based on an improved hidden Markov model (IHMM) with autocorrelated observations considered is developed. Due to the autocorrelations, the conventional expectation maximization (EM) algorithm for HMMs is not appropriate. A new EM algorithm is developed for the solutions. The Shewhart residual chart is employed in quality shift detection in conjunction with the IHMM. The residual is defined as the deviation of the predicted value by the IHMM and the current real observation value. Through the residual chart, we are able to see whether the process is in control. If the chart initiates an alarm, one running length is obtained. Thus, the average running lengths (ARLs) can be calculated with sufficient samples. The ARL is set as the comparison index for different models, including IHMM, HMM and AR.

The rest of this paper is organized as follows: Section 2 introduces the development of the IHMM model and its algorithm. In addition, the comparison of the prediction performances of different approaches is presented in Section 2. In Section 3, residual charts are introduced. In Section 4, numerical examples and interesting results are presented. The conclusions and possible areas for future research are given in Section 5.

2. Model Development

2.1. Hidden Markov Models

Denote a Markov chain with a finite state set $\left\{s_1,s_2,\cdots ,s_N\right\}$ by $\left\{S_n,n=1,2,\cdots \right\}$. Let $a_{ij}$ be the probabilities that the Markov chain enters state $s_j$ from state $s_i$ $\left(1\le i,j\le N\right)$ and ${\pi }_i=P\left\{S_1=s_i\right\},i=1,2,\cdots ,N$ be the initial state probabilities. Denote a finite set of signals by $\zeta$, and suppose a signal from $\zeta$ is sent each time the Markov chain enters a state. Suppose that when the Markov chain enters state $s_j$ independently of previous states and signals, the signal sent is $o_r$ with probability $p\left(o_r|s_j\right)$ that meets ${\sum }_{o_r\in \zeta }p\left(o_r|s_j\right)=1$. That is, if $O_n$ represents the $n$th signal, then it can be written by

$$P\left\{O_1=o_r|S_1=s_j\right\}=p\left(o_r|s_j\right),$$

(1)

$$P\left\{O_n=o_r|S_1,O_1,\cdots ,S_{n-1},O_{n-1},S_n=s_j\right\}=p\left(o_r|s_j\right).$$

(2)

Such a model, in which the signal sequence $O_1,O_2,\cdots$ can be observed while the underlying Markov chain state sequence $S_1,S_2,\cdots$ cannot be observed, is called a hidden Markov model [47].

Normally, an HMM contains the following elements [48]:

• A finite hidden state set $s=\left\{s_1,s_2,\cdots ,s_N\right\}$, where $N$ is the hidden state number.
• A set of possible observation values (signals) $q=\left\{q_1,\cdots ,q_K\right\}$, where $K$ is the number of possible observation values. Note that if the values of the observations are continuous, $K$ should be infinite.
• An observation sequence, $\mathit{o}=\left(o_1,\cdots ,o_T\right)$, where $T$ is the number of observations, $o_t\in q,1\le t\le T$.
• A distribution of state transition probabilities, $\mathit{A}=\left\{a_{ij}\right\}$, where

  $$a_{ij}=P\left\{S_{t+1}=s_j|S_t=s_i\right\},1\le i,j\le N,{\sum }_j{a}_{ij}=1.$$

  (3)
• A distribution of initial state probabilities, $\mathit{\pi }=\left\{{\pi }_i\right\}$, where

  $${\pi }_i=P\left\{S_1=s_i\right\},i=1,2,\cdots ,N,{\sum }_i{\pi }_i=1.$$

  (4)
• A conditional probability distribution of the observations, given $S_t=s_i$, $\mathit{B}=\left\{b_i\left(q_k\right)\right\}$, where

  $$b_i\left(q_k\right)=P\left\{o_t=q_k|S_t=s_i\right\},1\le i\le N,1\le t\le T,{\sum }_{q_k}{b}_i\left(q_k\right)=1.$$

  (5)

In general, an HMM contains three elements, denoted by $\lambda =\left\{\mathit{A},\mathit{B},\mathit{\pi }\right\}$.

2.2. The Improved HMM with Autocorrelated Observations

2.2.1. The Selection for the Order of Autocorrelation

Different from the traditional HMMs, in which current observations are assumed to be independent of previous states and observations, autocorrelated observations are considered. The observations are assumed to follow a Gaussian distribution, and the current observations are required to be dependent, not only on the current hidden state, but also on the previous observations.

The autocorrelated data from production may be multi-order, and it seems not cost-effective for judging this detailed order by using engineering experience and professional knowledge. Therefore, it is hoped to find a suitable order that will keep the implementation of the IHMM feasible, efficient and cost-effective.

The numerical experiment on the order selection for autocorrelated data in SPC application is designed as follows:

Step 1. Generate stationary autocorrelated data by a $d$th-order autoregressive model, $\mathrm{AR}\left(d\right),d\ge 2$;

Step 2. Use $\mathrm{AR}\left(p\right),p=d,d-1,\dots ,1$ models to fit the data, respectively, and obtain the parameters that need to be estimated by least-squares estimation;

Step 3. Generate 2000 stationary data sequences with a same length of 2000 under the processes with a mean shift magnitude of $\delta$;

Step 4. Predict each sequence by using $\mathrm{AR}\left(p\right),p=d,d-1,\dots ,1$ models, respectively, and calculate their residuals;

Step 5. Determine the central lines (CLs), upper control limits (UCLs) and lower control limits (LCLs) of $p$ residual charts;

Step 6. Generate stationary autocorrelated data by a $d$th-order autoregressive model;

Step 7. Calculate the average running lengths (ARLs) of $p$ residual charts.

In this study, twelve $\mathrm{AR}\left(2\right)$ and eight $\mathrm{AR}\left(3\right)$ models are used to generate data. The control limit coefficients of residual charts are 3. The shift magnitudes of the mean are set by 0, 1.5 and 3, respectively.

When data are generated by $\mathrm{AR}\left(2\right)$ models for an in-control process, the ARL of the residual charts is 371.3803 if the data are fitted by $\mathrm{AR}\left(2\right)$ models, and the ARL is 361.1688 if the data are fitted by $\mathrm{AR}\left(1\right)$ models. Thus, the probability of the type I error is increased by 2.83%. Similarly, when data are generated by $\mathrm{AR}\left(3\right)$ models, the ARL of the residual charts is 371.8250 if the data are fitted by $\mathrm{AR}\left(3\right)$ models, the ARL is 358.5401 if the data are fitted by $\mathrm{AR}\left(2\right)$ models, and the ARL is 345.8113 if the data are fitted by $\mathrm{AR}\left(1\right)$ models. Thus, the probabilities of the type I error are increased by 3.7% and 7.51%, respectively.

The ARLs of different situations are shown in Figure 1, in which the numbers on the x-axis represent the combinations of autocorrelation coefficients.

From Figure 1, it is seen that generally, AR(1) models outperformed other AR(d), $d\ge 2$ models in shift detection regardless of the autocorrelation order. Although AR(1) models lead to a slight increase in type I error, it seems to be insignificant compared with their good performances in detecting quality shifts. Therefore, only the first-order autocorrelation is considered in this study. As a result, there is no need to judge what order the autocorrelation is so that the modeling cost can be saved.

2.2.2. The Development of the IHMM

According to the analysis in Section 2.2.1, the construction of the IHMM, in which the current observation is related to its previous observation, is shown in Figure 2.

In Figure 2, arrows pointing to the right indicate the correlation between neighboring states or observations. It is intuitive that the IHMM will become a traditional HMM if there is no correlation between neighboring observations.

Let $o_i\left(t\right)$ be the observation value at time $t$ given state $s_i$; $o_i\left(t\right)$ can be fitted with the following function:

$$o_i\left(t\right)={\varsigma }_i+c_i{o}_{t-1}+{\epsilon }_i,1\le i\le N,1\le t\le T,$$

(6)

where $c_i$ is the first-order autocorrelation coefficient, ${\varsigma }_i$ is the constant term and ${\epsilon }_i$ is white noise, following a normal distribution with mean zero and variance ${\sigma }_i^2$.

Let ${\mu }_i\left(t\right)$ be the mean of $o_i\left(t\right)$ given state $s_i$, ${\mathit{x}}_t={\left(1,o_{t-1}\right)}^{\prime }$, then Equation (6) can be written as:

$${\mu }_i\left(t\right)={\mathit{C}}_i{\mathit{x}}_t,$$

(7)

where ${\mathit{C}}_i$ is a $1\times 2$ matrix consisting of $\left\{{\varsigma }_i\right\}$ and $\left\{c_i\right\}$.

Hence, we conclude that the states are a Markov process, and that the conditional observations given state $s_i$ follows a Gaussian distribution with mean ${\mathit{C}}_i{\mathit{x}}_t$ and variance ${\sigma }_i^2$, i.e.,

$$P\left\{S_{t+1}|S_1,\cdots ,S_t,o_1,\cdots ,o_t\right\}=p\left(S_{t+1}|S_t\right),$$

(8)

$$b_i\left(o_t|o_{t-1}\right)=P\left(o_t|o_{t-1},S_t=s_i\right)=\frac{1}{\sqrt{2\pi }{\sigma }_i}\exp \left(-\frac{{\left({\mathit{o}}_t-{\mathit{C}}_i{\mathit{x}}_t\right)}^2}{2{\sigma }_i^2}\right).$$

(9)

If the current observation is not related to its previous observations, ${\mu }_i\left(t\right)$ is the constant ${\varsigma }_i$. Thus, the IHMM becomes a traditional HMM.

2.3. Parameter Estimation

The aim of using an IHMM is to forecast observation values. Firstly, we estimate the parameters to obtain an optimal $\hat{\lambda }$ by maximizing the probability $P\left(\mathit{o}|\lambda \right)$.

The EM algorithm is a popular method to estimate the parameters for HMMs. However, since autocorrelation between observations is considered, the traditional algorithm could not be used directly. We redefine the parameters as $\lambda =\left\{\mathit{A},\mathit{C},{\mathit{\sigma }}^2,\mathit{\pi }\right\}$ where $\mathit{C}=\left\{{\mathit{C}}_i\right\},{\mathit{\sigma }}^2=\left\{{\sigma }_i^2\right\}$, and the definitions of ${\mathit{C}}_i$ and ${\sigma }_i^2$ are provided in Section 2.2.2. The flowchart demonstrating how to estimate the parameters in IHMM with an improved EM algorithm is shown in Figure 3.

Firstly, we introduce two types of probabilities: forward probability, $F_t\left(i\right)$, and backward probability, $B_t\left(i\right)$, defined as:

$$F_t\left(i\right)=P\left(o_1,\cdots ,o_t,S_t=s_i|\lambda \right),$$

(10)

$$B_t\left(i\right)=P\left(o_{t+1},\cdots ,o_T|S_t=s_i,\lambda \right).$$

(11)

The calculation of the two probabilities here is similar to that of traditional HMMs. Slightly different from traditional HMMs, the initial values are defined by $F_{d+1}\left(i\right)={\pi }_i{b}_i\left(o_{d+1}|o_1,\cdots ,o_d \right)$ and $B_T\left(i\right)=1$, respectively.

Based on the forward and backward probabilities, some intermediate probabilities are computed. Given $\lambda$ and $\mathit{o}$, denote the probability that the process is in state $s_i$ at time $t$ by ${\gamma }_t\left(i\right)$ and the probability that the process is in state $s_i$ at time $t$ and in state $s_j$ at time $t+1$ by ${\xi }_t\left(i,j\right)$. ${\gamma }_t\left(i\right)$ and ${\xi }_t\left(i,j\right)$ can be written as Equations (12) and (13). Please refer to reference [48] for the details of the derivations.

$${\gamma }_t\left(i\right)=\frac{F_t\left(i\right)B_t\left(i\right)}{{\sum }_j{F}_t\left(i\right)B_t\left(i\right)},$$

(12)

$${\xi }_t\left(i,j\right)=\frac{F_t\left(i\right)a_{ij}{b}_j\left(o_{t+1}\right)B_{t+1}\left(j\right)}{{\sum }_{i=1}^N{\sum }_{j=1}^N{F}_t\left(i\right)a_{ij}{b}_j\left(o_{t+1}\right)B_{t+1}\left(j\right)}.$$

(13)

Next, we develop an improved EM algorithm for the IHMM, step by step.

Step 1. Determine the log-likelihood function for complete data.

The complete data are $\left(o,S\right)=\left(o_1,\cdots ,o_T,s_1,\cdots ,s_T\right)$, and its log-likelihood function is $\log P\left(o,s|\lambda \right)$.

Step 2. E-step: determine $Q$ functions.

$$Q_1\left(\lambda ,{\lambda }^{\left(n\right)}\right)=E_s\left[\log P\left(\mathit{o},\mathit{s}|\lambda \right)|\mathit{o},{\lambda }^{\left(n\right)}\right]={\displaystyle \sum }_s\log P\left(\mathit{o},\mathit{s}|\lambda \right)P\left(\mathit{o},\mathit{s}|{\lambda }^{\left(n\right)}\right),$$

(14)

where $\lambda$ is the parameter that will maximize the $Q_1$ function, while ${\lambda }^{\left(n\right)}$ is the current parameter value.

With $P\left(\mathit{o},\mathit{s}|\lambda \right)={\pi }_{s_1}{b}_{s_1}\left(o_1\right)a_{s_1{s}_2}{b}_{s_2}\left(o_2\right)\cdots a_{s_{T-1}{s}_T}{b}_{s_T}\left(o_T\right)$, $Q_1\left(\lambda ,{\lambda }^{\left(n\right)}\right)$ can be written by

$$\begin{array}{c}{Q}_1\left(\lambda ,{\lambda }^{\left(n\right)}\right)={\displaystyle \sum }_s{\pi }_{s_1}\log P\left(\mathit{o},\mathit{s}|{\lambda }^{\left(n\right)}\right)\\ +{\displaystyle \sum }_s\left({\displaystyle \sum }_{t=1}^{T-1}\log a_{s_t{s}_{t+1}}\right)P\left(\mathit{o},\mathit{s}|{\lambda }^{\left(n\right)}\right)\\ +{\displaystyle \sum }_s\left({\displaystyle \sum }_{t=1}^T\log b_{s_t}\left(o_t\right)\right)P\left(\mathit{o},\mathit{s}|{\lambda }^{\left(n\right)}\right).\end{array}$$

(15)

For the re-estimation of $b_i\left(q_k\right)$, one more auxiliary function $Q_2$ is proposed by taking the conditional expectation of the log-likelihood of the observation sequence:

$$Q_2\left(\lambda ,{\lambda }^{\left(n\right)}\right)={\displaystyle \sum }_s{\displaystyle \sum }_{t=1}^T{\gamma }_t\left(i\right)\left(\ln \left(\frac{1}{\sqrt{2\pi }}\right)+\ln \frac{1}{{\sigma }_i}-\frac{{\left({\mathit{o}}_t-{\mathit{C}}_i{\mathit{x}}_t\right)}^2}{2{\sigma }_i^2}\right).$$

(16)

Step 3. M-step: re-estimation.

Re-estimate $\lambda$ that maximizes $Q_1\left(\lambda ,{\lambda }^{\left(n\right)}\right)$, that is

$${\lambda }^{\left(n+1\right)}=\arg \underset{\lambda }{\max }{Q}_1\left(\lambda ,{\lambda }^{\left(n\right)}\right).$$

(17)

The state transition probabilities are derived as:

$$a_{ij}{}^{\left(n+1\right)}=\frac{{\sum }_{t=1}^{T-1}{\xi }_t\left(i,j\right)}{{\sum }_{t=1}^{T-1}{\gamma }_t\left(i\right)}.$$

(18)

The initial state probabilities are derived as:

$${\pi }_i{}^{\left(n+1\right)}={\gamma }_1\left(i\right).$$

(19)

By re-estimating $\lambda$ that maximizes $Q_2\left(\lambda ,{\lambda }^{\left(n\right)}\right)$, $c_{i,\tau }{}^{\left(n+1\right)},{\varsigma }_i{}^{\left(n+1\right)},{\sigma }_i^2{}^{\left(n+1\right)}$ can be written as:

$$c_i{}^{\left(n+1\right)}=\frac{{\sum }_{t=1}^T{\gamma }_t\left(i\right)o_{t-1}\left(o_t-\frac{{\sum }_{t=1}^T{\gamma }_t\left(i\right)o_t}{{\sum }_{t=1}^T{\gamma }_t\left(i\right)}\right)}{{\sum }_{t=1}^T{\gamma }_t\left(i\right)o_{t-1}\left(o_{t-1}-\frac{{\sum }_{t=1}^T{\gamma }_t\left(i\right)o_{t-1}}{{\sum }_{t=1}^T{\gamma }_t\left(i\right)}\right)},$$

(20)

$${\varsigma }_i{}^{\left(n+1\right)}=\frac{{\sum }_{t=1}^T{\gamma }_t\left(i\right)o_t-{\sum }_{t=1}^T{\gamma }_t\left(i\right)o_{t-1}{c}_i{}^{\left(n+1\right)}}{{\sum }_{t=1}^T{\gamma }_t\left(i\right)},$$

(21)

$${\sigma }_i^2{}^{\left(n+1\right)}=\frac{{\sum }_{t=1}^T{\gamma }_t\left(i\right){\left(o_t-c_i{}^{\left(n+1\right)}{o}_{t-1}-{\varsigma }_i{}^{\left(n+1\right)}\right)}^2}{{\sum }_{t=1}^T{\gamma }_t\left(i\right)}.$$

(22)

Repeat Step 2 and Step 3 until the log-likelihood function converges.

2.4. Prediction

Once the parameter is determined by the improved EM algorithm, the model can be employed to forecast the expected value of the next observation.

The conditional probability distribution of the observation $o_{T+1}$ can be derived as

$$\begin{array}{c}P\left(o_{T+1}|o_T,\cdots ,o_1,\lambda \right)={\displaystyle \sum }_{i=1}^{N}P\left(o_{T+1}|S_T=s_i,o_T,\cdots ,o_1,\lambda \right)P\left(S_T=s_i|o_T,\cdots ,o_1,\lambda \right)\\ ={\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^{N}P\left(o_{T+1}|S_{T+1}=s_j,o_T,\cdots ,o_1,\lambda \right)P\left(S_{T+1}=s_j|S_T=s_i,\lambda \right){\gamma }_T\left(i\right)\\ ={\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^N{\gamma }_T\left(i\right)a_{ij}{b}_j\left(o_{T+1}|o_T,\cdots ,o_{T+1-d}\right).\end{array}$$

(23)

Therefore, the expectation of $o_{T+1}$ is computed by

$${\widehat{o}}_{T+1}=\int o_{T+1}P\left(o_{T+1}|o_T,\cdots ,o_1,\lambda \right)d{o}_{T+1}.$$

(24)

Given an observation sequence $\mathit{o}=\left(o_1,\cdots ,o_T\right)$, ${\widehat{o}}_t$ is predicted by

$${\widehat{o}}_t=\int o_{t}P\left(o_t|o_{t-1},\lambda \right)d{o}_t=\int o_t{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^N{\gamma }_{t-1}\left(i\right)b_j{a}_{ij}\left(o_t|o_{t-1}\right)d{o}_t,t\ge d.$$

(25)

Thus, the predicted values of the observations can be denoted by $\hat{\mathit{o}}=\left(o_1,{\widehat{o}}_2,\cdots ,{\widehat{o}}_T\right).$

2.5. Performance Comparison

The mean squared error (MSE), absolute mean error (AME) and mean absolute percentage error (MAPE) are common tools for measuring, fitting and predicting accuracy [49]. Both MSE and AME values determine the average deviation between fitting values and original values, while MAPE provides a measurement for testing the relevant difference between them. In this study, we use MSE as the criterion to evaluate the models. The equation for MSE is:

$$\mathrm{MSE}=\frac{1}{LT}{\sum }_{l=1}^L{\sum }_{t=1}^T{\left(o_t-{\widehat{o}}_t\right)}^2,$$

(26)

where $L$ is the number of predicted samples, and $T$ is the length of each sample.

We assume that a variable from a production process follows a normal distribution with mean 100 and variance 25 when the process is under control, and that the observations are first-order autocorrelated with a correlation coefficient of 0.6. We use IHMM, HMM and AR(1) methods to predict observation values, respectively. The MSEs for the three models are 15.1276, 16.1867 and 15.7861, respectively. Since these MSEs are very close, we conclude that the predicted performances of the three approaches are similar. The predicted results of an observation sequence with a length of 50 from the in-control process are shown in Figure 4, from which we can see that the three models have close performances. However, the time taken for prediction using the three models are quite different. For the prediction of an observation sequence with a length of 50, IHMM is 14.6523 s, HMM is 93.6521 s, and AR(1) is almost instantaneous under the environment of win10 OS (Microsoft, Redmond, WA, USA) with CPU of Intel(R) Core(TM) i7-7500U (Santa Clara, CA, USA).

Then, we suppose the process has a shift magnitude of 3. We still use the three methods to predict observation values, respectively. The predicted results of an observation sequence with a length of 50 from the out-of-control process are shown in Figure 5, from which we can see distinctly different performances of the three models. By observing the distances between the lines with different colors, obviously, if the MSEs are calculated, the MSE from AR(1) is much less than that from IHMM, and the MSE from IHMM is much less than that from HMM. The IHMM only results in a medium-level performance in the prediction for the autocorrelated process. However, it is very interesting that the performances of corresponding residual charts have the best performances in detecting quality shifts. This seems to suggest that the residual charts integrating IHMM can achieve a surprising effect. This result is verified in Section 3.

3. Statistical Process Control with Residual Chart

Residual control charts are an effective tool for online monitoring in the presence of autocorrelations. A residual chart called $e$ chart is developed in our study.

Residuals are obtained by subtracting the predicted values of observations from the original values, that is $\mathit{e}=\mathit{o}-\widehat{\mathit{o}}=\left(e_1,\cdots ,e_T\right).$ The control limits of the $e$ chart are given by

$$UCL={\mu }_e+k{\sigma }_e,$$

(27)

$$LCL={\mu }_e-k{\sigma }_e.$$

(28)

where $UCL$ represents the upper limit, while $LCL$ represents the lower limit,$k$ is the number of ${\sigma }_e$, ${\mu }_e$ represents the mean of $e$, and ${\sigma }_e$ represents the standard deviation of $e$. ${\mu }_e$ and ${\sigma }_e$ can be obtained by simulations based on sufficient samples.

If the value of $e_t$ drops within $UCL$ and $LCL$, the process is judged to be in control; otherwise, it is judged to be out of control.

4. Numerical Examples

We consider that the variable from a production process followed a normal distribution with mean 100 and variance 25 when the process is in control and that the observations were first-order autocorrelated. The correlation coefficient varies between −0.6 and 0.6 with increments of 0.3. Two shift magnitudes of 1.5 and 3 are considered. According to the definition of residual charts, the ARLs of in-control processes for all predicted methods are 370, so we focus our discussion on the out-of-control processes. By conducting multiple experiments, we find that it is appropriate to make the state number with 5 for both IHMM and HMM. The experimental results are shown in Figure 6 and Figure 7.

As shown in Figure 6 and Figure 7, when correlation coefficient changes from positive to negative, ARLs decrease dramatically, regardless of the approach used. Compared with positive correlations, the ARLs of negative correlations are relatively very small, and ARLs obtained by different models are very close to each other. Thus, the following discussions focus on positive correlations.

As pointed out in Section 3, although the predictions of both the IHMM and HMM are inferior to AR(1) models, the performances of residual charts from the former models are much better than the latter ones. As seen in Figure 4 and Figure 5, when the coefficients are larger than zero, the ARLs by IHMMs are shorter than those by HMMs, and by HMMs shorter than by AR(1) models.

As correlation coefficients increase, the ARLs generally increase, regardless of the approach used, especially as the shift magnitude decreases.

Generally speaking, when detecting quality shifts, the performances of IHMM, HMM and AR(1) models are ranked with IHMMs first, HMMs second and AR(1) last. Moreover, as pointed out in Section 2, the times taken by IHMMs are much shorter than HMMs under the same running environments.

5. Conclusions

In this paper, an IHMM with autocorrelated observations and a new EM algorithm are proposed. Residual charts in conjunction with the IHMM are employed for detecting quality shifts. The results demonstrates that: (1) the IHMM outperforms the HMM and AR(1) method with positive correlations; (2) the IHMM has similar performances with the HMM and AR(1) methods with negative correlations; (3) compared with positive correlations, the ARLs of the IHMM under negative correlations are relatively very small, as well as those of the HMMs and AR(1) models; (4) the IHMMs take a much shorter time than HMMs, for both training and prediction, but still longer than the AR(1) models.

Future research might focus on further experimental validations for the IHMM and its algorithm. The strict Gaussian distribution of observations could be extended to other probability distributions. Since multistage systems are commonplace in the manufacturing industry, it is worth extending this approach in this research direction.