Real Time Statistical Process Control for Autocorrelated Serial Data: A Simulation Approach

Computer measurement systems play an important role on process automation and Industry 4.0 implementation strategies. They can be easily integrated on modern production systems, enabling real time test and control of multiple product and process characteristics that need to be monitored. If for one side the big data provided by these systems is an important asset for production analytics and optimization, on the other hand, the high frequency data sampling, commonly used in these systems, can lead to autocorrelated data violating, this way, statistical independence requirements for statistical process control implementation. In this paper we present a simulation model, using digital recursive filters, to properly handle and deal with these issues. The model demonstrates how to eliminate the autocorrelation from data time series, creating and ensuring the conditions for statistical process control application through the application of real time control charts. A performance comparison between Shewhart of Residuals and Exponentially Weighted Moving Average (EWMA) of Individual Observations control charts is made for autocorrelated data time series with the presence of different mean shift amplitude perturbations.


I. INTRODUCTION
HE level of competitiveness experienced in modern factories are pushing even further the level of acuity and need for robust process control that ensure not only the quality of the end products but, simultaneously, the efficiency of the manufacturing processes.The Industry 4.0 [1] is a reflex of this environment where the digitalization and automation of manufacturing processes are part of the evolutionary path.If for one side the availability of the big data [2,3] enables the possibility of information generation to scrutinize processes, improve product quality and process efficiency and ultimately improve competitiveness, on the other side the way data collection is made leads to another type of problem, the serial data autocorrelation [4].Serial data autocorrelation can be found in many continuous and discrete manufacturing processes due to the high frequency data acquisition imposed by computer measurement systems and automated technology used for process sensing and measurement [5].In consequence of that, modern manufacturing industries have changed the effectiveness of traditional statistical process control (SPC) techniques, namely when process data violates the requirement of statistical data independence.The high frequency data acquisition not only leads to autocorrelated data, thus violating the assumptions of statistical data independence, but also makes it difficult to provide the real time statistical information.
Several authors have been suggesting methods to deal with autocorrelation, [6][7][8][9].All of them propose the same strategy to eliminate the serial autocorrelation from process data, which includes the fitting of an appropriate time series model to the process observations, followed by the SPC application to the sequence of obtained residuals.Despite of the effectiveness of time series modeling to eliminate autocorrelation, the lack of recursiveness of these models makes them fail when real time data control is required.
To properly handle these issues, an approach for real time process control is proposed.A simulation model was developed using the LabVIEW applying recursive digital filters, state space models and Z transform.Through the simulation model it is demonstrated how to extract serial data autocorrelation from process data using recursive digital filters, restoring, this way, the statistical data independence necessary for SPC application.The simulation model developed implements the SPC through the implementation of real time processing control charts, such as Shewhart of residuals and EWMA of individual observations.A comparison of effectiveness to detect process perturbations between both control charts is presented.

II. THEORETICAL BACKGROUND
This section presents a brief theoretical background used in the paper development and simulation model.

A. TIME SERIES MODELING
To model the serial autocorrelated process, it's proposed the use of difference equations [10] designated by autoregressive integrated moving average (ARIMA) models.An ARIMA model of (p,d,q) order for a process variable  can be represented by where B is the backward shift operator defined as   =  .∅() is the autoregressive (AR) operator of order p defined by ∅() = 1 − ∑ ∅  and () is the moving average (MA) operator of order q defined by () = 1 − ∑   ; where ∇=  −  defines the backward difference operator and ∇ = (1 − ) and  represents the order of differentiation.The process is said stationary for  = 0 and nonstationary when  ≥ 1;  represents a white noise variable where  ∼ (0,  ).Denoting the predicted value obtained from an appropriately fitted ARIMA model by  , the residuals  =  −  will behave like independent and identically distributed random variables.

B. THE Z TRANSFORM
The  transform is an analytical tool extremely useful in the analysis of discrete-time systems and computer modeling based on difference equations and transfer functions [11]).Taking () as a quality characteristic of a certain manufacturing process, the  transform of the discrete variable  is given by where  represents the transform operator and Z is the argument defined by a complex variable.On the other hand, the inverse  transform of () is given by The proof of the equivalence between the  transform and the back shift operator   =  is presented in [12] as follows As can be seen, the delay imposed by the backshift operator  to the variable  in the time domain, corresponds, in the  domain, to multiplying the  transform of the original function  by  .

C. RECURSIVE DIGITAL FILTERS
A digital filter () can be described as the implementation of an algorithm that computes a sequence of outputs () from a sequence of inputs ().The most common classes of digital filters are the Finite Impulse Response (FIR), also called moving average (MA) or all-zero filters, and the Infinite Impulse Response (IIR).The IIR filters can be split in autoregressive (AR) or all-pole filters and autoregressive moving average (ARMA) filters, which have both poles and zeros.ARMA digital filters can be represented using the Z transform, [13,14].

D. ARIMA DIGITAL FILTERS
ARMA digital filters can be represented by  transform using the following transfer function [13,14] In the previous equation, ∅ and  are, respectively, the n th and m th order parameters of AR(p) and MA(q) processes, () is the output variable and () is the input white noise stream.Therefore, the ARMA () digital filter will be of the form where the AR component of the digital filter H(z), called IIR, is given by and the MA component of the digital filter H(z), or FIR, will be

E. INVERSE DIGITAL FILETERS AND TIME SERIES MODELING
Recursive digital filters can be used effectively for modeling serial autocorrelated processes.On the other hand, when appropriate inverse recursive digital filters are applied to the autocorrelated observations  , the autocorrelation can be recursively removed and a sequence of normal and independently distributed errors  are obtained enabling, this way, the application of real time SPC.Considering the difference equation of an ARMA process defined as the respective recursive digital filter will be defined as with − = 1.
Considering now the difference equation of one-step-ahead prediction of an ARMA process [11] such as  = ∅  + ⋯ + ∅  −  −  , (11) the inverse digital filter becomes The prediction error () is found through the difference between the observed value and the one-step-ahead forecast as follows as − = 1 then The application of the inverse Z transform to (15) leads to that is,

F. ARIMA PROCESS REPRESENTED AS STATE-SPACE-MODEL
The space-state models (SSM) are based on recursive methods [15] and might be quite useful when recursive modeling and simulation of ARIMA processes are required.Due to their recursive properties, the SSM can be used in a real time applications once the model is updated each time a new observation becomes available.The SSM have been employed by several authors for simulation purposes [16,17].
According to [18], an ARMA (p,q) process can be represented by the following discrete state-space equations where [] is the state matrix of r x r dimension and ∅ are the appropriate coefficients of ∅(); [] defines an 1 x r vector containing the coefficients of () and [] defines a r x 1 vector of known constants; () is an r x 1 vector representing the state of process variables,  ∼ (0,  ) is the white-noise, and [] equals 0 or 1 depending on whether the process is or isn't subject to an additional white-noise term respectively.In our case we assume D = 0. Finally, , where p and q represent, respectively, the order of the autoregressive ∅() and moving average () polynomials.
SSM can also be used to represent ARIMA (p,d,q) processes [17].For that the ∅ , ∅ , … , ∅ parameters must be replaced by the outcome parameters of the representation (1 − ) ().The expanded transition and observation equations are defined respectively by (19) and (20).

G. SPC STRATEGIES
Traditional SPC [19] assumes that the observations a certain process characteristic  , obtained by the sample of order t are defined by where  defines the process mean and  represents the white noise error, where  ∼ (0,  )    ,  = 0 for any  ≠ 0, meaning that the sequence of errors  is statistically independent.This ensures that  follow as iid process (independent and identically distributed), with  ∼ (,  ) and   ,  = 0 for any for any  ≠ 0. This is a fundamental requirement for SPC application.
There is a close connection between control charts and hypothesis testing [19], null hypothesis against alternative hypothesis.Let´s assume that the  sample of the process characteristic has an average of  .If the value  is bounded by the control limits ( ± 3 ) we accept null hypothesis ( :  =  ).Under these conditions the process mean  is under control having only normal causes of variation.If the value  fall´s outside of the control limits ( ± 3 ) we reject null hypothesis ( :  =  ) and accept alternative hypothesis ( :  ≠  ).These conditions highlight the present of special causes of variation on the process mean  .
The presence of serial autocorrelation has a negative effect on the SPC application effectiveness as decreases the ARL VOLUME 22(2), 2023 (Average Run length) [6,7].This increases the level occurrences of Type I statistical error.Some authors, [6][7][8][9] demonstrated the effect of serial autocorrelation and presented appropriated strategies to eliminate undesired autocorrelation effects on control charts.

H. RESIDUALS CONTROL CHARTS
One of the strategies followed was to model the autocorrelative structure applying the appropriated ARIMA model to the autocorrelated time series.When an appropriate ARIMA model is fitted to an autocorrelated data time series, a sequence of independent and identically distributed residuals  =  −  is obtained.Once restored the iid conditions, the traditional Shewhart control charts can be applied to the residuals stream of data time series.These are called the residual control charts [19].Therefore, whenever a shift in the mean occurs, highlighting the presence of a special cause of variation, it will be transferred to the residuals and detected by the control charts [6].
The EWMA statistic [19] for a certain process characteristic  , obtained at order t, is defined by where t = 1, 2, …, and 0 <  < 1 is the weighted constant working as a filter to evidence recent or past observation.If the  data time series are uncorrelated, the control limits for the EWMA control chart for individual observations can be defined by the Lower and Upper Control Limit respectively, LCL and UCL, [8,19].
The EWMA can be used in situations where the data time series is autocorrelated, in specifically when the data follows a ARIMA (0,1,1) = IMA (1,1) process [8].
According to these authors the EWMA with  = 1 −  is the optimal one-step-ahead forecast for these processes, where  =  is the forecasted observation of instant t + 1 made at instant t.Therefore, we also have  =  .Using the above identities and replacing in (22)  ( Under these circumstances, the one-step-ahead prediction errors  =  −  , are independent and identically distributed with zero mean and standard deviation  .Control charts with control limits at ±3 can be applied to these one-step-ahead prediction errors.
According to [8], the EWMA can produce excellent onestep-ahead prediction for other ARIMA processes, if the observations are positively correlated, and the process mean does not drift too quickly.
Autocorrelated data requires a different approach for selecting  to achieve a certain average run length, when compared with uncorrelated data.The authors [19,27] proposed a method to select  based on the minimization of the sum of squares of the one-step-ahead prediction errors (SSE).According to these authors, the estimate of the variance of one-step-ahead errors  can be calculated by the quotient between the sum of squares of prediction errors for the optimal  and the  observations used for its determination.This method can be easily implemented in computer time series analysis.

III. SIMULATION MODEL
The Real Time simulation model presented in this work (Figure 1) was developed in LabVIEW, a system engineering graphical programing language widely used in virtual instrumentation and computer measurement applications.This simulation model can simulate any ARIMA process using the SSM approach from a white noise (0,1) base line.This way once created the desired autocorrelation data structure, it is also possible to eliminate the autocorrelation using digital recursive filters (FIR and IRR), reconverting the data into i.i.d.residuals.The model explores the application of individual observations control charts to the residuals and proper set of EWMA charts to autocorrelated data time series.The mean shifts, simulated through a step degree function, enable the evaluation of the two types of charts for different process perturbation amplitude shifts imposed at instant t = 2s.The sum of squares errors (SSE) algorithm was also implemented to support the selection of the optimum  to be used in the EWMA statistics.
and it has the following SSM representation CASE 2 -ARMA (1,1), autoregressive moving average model with parameters ∅ = 0.8 and  = 0.3 can be written by with the following SSM representation, In both cases the LabVIEW simulation model was properly configured with parameters of (32), (33), ( 35) and (36).The FIR and IIR digital filters were setup with the appropriated values to enable the autocorrelation elimination.Figure 2 and Figure 3 represents, respectively, the control chart for individual observations for the autocorrelated data time series Y(t) (blue line) and for the residuals e(t) (green line) for case 1 and case 2. As we can observe by the number of points from the autocorrelated time series Y(t) falling outside of control limits (UCL, LCL), the Type I statistical errors are significantly high on this case in comparison with the residuals e(t).
The parameter  of the EWMA control charts was selected, in both cases, using the minimum sum of square errors (SSE) of a(t), obtained from the one-step-ahead prediction errors  =  −  , according to (28) and (29).Figure 4 and Figure 5 illustrate the SSE function of the  for AR (1) and ARMA (1,1) cases, respectively.As it can be observed the values of  that minimizes the SSE on the AR (1) data time series is  = 0.7 while on ARMA (1,1) process is  = 0.45.

IV. DISCUSSION RESULTS
The data in Table 1 summarizes the results obtained during the simulations with the use of EWMA and Individual Observation of Residuals Control Charts for both processes, AR (1) and ARMA (1,1), for all the mean shifts.As expected, the EWMA Control Chart is more effective detecting small mean shifts in comparison with the Individual Observation of Residuals Control Chart for any of the processes AR (1) and ARMA (1,1) considered in the experience.Not only the detection is made earlier with the EWMA Control Chart, but it also provides, consistently, more points out of control limits after the mean shift perturbation.For a small mean shift of 0.25 the Individual Observation of Residuals Control Chart is not able to detect the perturbation, while the EWMA Control Chart detects the first point out of control limits at instant t = 3.38s and t = 4.01s respectively, for AR (1) and ARMA (1,1) processes.The effectiveness of the Residual Control Chart detecting the first point out of control improves with the increment of the mean shifts.Even though for mean shifts equal or bigger than 0.75 the Residuals and EWMA Control Charts shown a similar capability in detecting the first point out of control, EWMA shown more consistence in detecting points out of control after the mean shift.
While for the selected AR (1) and ARMA (1,1) processes the Residual Control Chart doesn't show differences in it performance as the residual series e(t) will be equal in both processes after autocorrelation filtering, the EWMA seams to react a bit slowly to the mean shifts on ARMA (1,1) in comparison with AR (1) process.Regarding the number of points out of control, the EWMA Control Chart show a good consistence on AR (1) and ARMA (1,1) processes, although the number of out-of-control points is bigger in the AR (1).Also, we can observe that EWMA tend to react more quickly in detecting points out of control as mean shifts perturbation amplitude increases.

TABLE 1. Summary of Simulation Results
Figure 8 and Figure 9 show simulated EWMA and Residual Control Charts for the AR (1) process with mean shift of 0.25, respectively.As referred, the Residual Control Chart is not able to perform any detection of points out of control after the perturbation at t = 2s, while EWMA Control Chart get the first point out of control at t = 3.38s.
The EWMA and Residual Control Charts for the AR (1) process with mean shift of 1 are shown on Figure 10 and Figure 11, respectively.As it can be seen and as per table 1, the detection of the mean shift in these two cases is made immediately after the mean shift perturbation.Note as EWMA reaction time to mean shift perturbation increased with shift amplitude (Figure 8, Figure 10).

V. CONCLUSIONS
A simulation model was developed to illustrate an approach for dealing with real time statistical process control with the presence of time series data autocorrelation very often found in technological advanced industrial environments.The model enables the generation of any ARIMA process.In this paper two autocorrelated processes were evaluated, a first order autoregressive AR (1) and a first order autoregressive moving average ARMA (1,1).The results show that appropriate digital recursive filters could be effectively used to remove autocorrelation in a real time environment enabling, that way, further application of statistical process control without violating the principle of statistical independence.Though the application of different mean shift amplitude perturbations, imposed at a preset time to each one of the autocorrelated processes, we evaluated the effectiveness of the application of both, Residuals of Individual Observations and EWMA control charts.This evaluation of performance consisted in the comparison of the time response achieved with the two control charts after the mean shift perturbation has been imposed to the time series data and by the consistence of points out of control after the occurrence of this perturbation.The results obtained are in line with theoretical expectations.For small amplitude mean shifts EWMA control chart shown a quicker time reaction to perturbation, exposing consistently points out of control limits after that point.Residuals Individual Observation control charts tend to improve reaction time to mean shift perturbations as its amplitude increases.They show similar time responses to EWMA for bigger amplitude mean shifts, in spite its consistence exposing points out of control is much smaller in comparison with EWMA control charts.

Figure 1 .
Figure 1.LabVIEW Simulation model for real time SPC.

Figure 2 .
Figure 2. Individual control chart for autocorrelated data time series Y(t) and residuals e(t) of AR(1) model.

Figure 3 .
Figure 3. Individual control chart for autocorrelated data time series Y(t) and residuals e(t) of ARMA (1,1) model.

Figure 8 .
Figure 8. EWMA Control Chart for AR (1) process with mean shift amplitude of 0.25.

Figure 9 .
Figure 9. Residual of Individual Observations Control Chart for AR (1) process with mean shift amplitude of 0.25.

Figure 10 .
Figure 10.EWMA Control Chart for AR (1) process with mean shift amplitude of 1.

Figure 11 .
Figure 11.Residual of Individual Observations Control Chart for AR (1) process with mean shift amplitude of 1.