Multioutput least square SVR-based multivariate EWMA control chart: The performance evaluation and application

Autocorrelation leads to a bias estimator of standard control charts. It is important to develop control chart that allows autocorrelation and to evaluate its performance. The objective of this paper is to evaluate the performance of multioutput least square support vector regression (MLS-SVR)-based multivariate exponentially weighted moving average (MEWMA) control chart for monitoring multivariate autocorrelated data. For first order of vector autoregressive (VAR) and first order of vector moving average data, the proposed control chart tends to yield stable in-control average run length at about 200. The proposed control chart becomes more insensitive due to the increase of MEWMA smoothing parameter. In the real application, the proposed method is successfully applied to monitor water turbidity and chlorine residual data in the drinking water manufacturing process. Subjects: SPC/Reliability/Quality Control; Multivariate Statistics; Quality Control & Reliability


Introduction
Traditional control chart usually assumes that variables are statistically independent over time. However, data collected in time often show serial dependency and always affect false alarm rate and shift detection power. Consequently, a control chart developed under the assumption of independent observation would have poor performance when it is applied to monitor Muhammad Mashuri ABOUT THE AUTHORS The corresponding author, Dr. Muhammad Mashuri, is an associate professor in the Department of Statistics, Institut Teknologi Sepuluh Nopember, Surabaya, Indonesia. His research interest includes statistical quality control and multivariate analysis, especially in industrial field. Dr. Suhartono and Dr. Dedy Dwi Prastyo are senior lecturers whose research interest includes time series analysis as well as machine learning. Hidayatul Khusna and Muhammad Ahsan are PhD students through Master Program of Education Leading to Doctoral Degree for Excellent Graduates (PMDSU), an accelerated programme for undergraduate prepared to be candidate lecturers or researchers with doctoral degrees.

PUBLIC INTEREST STATEMENT
A control chart is one of the improvement tools to monitor how a process changes over time. A good control chart will provide precise conclusion if the special causes of variation affect the actual process. The multivariate control chart for monitoring variables with time dependency has been developed. In this study, one of the multivariate control charts for monitoring autocorrelated data named MLS-SVR-based MEWMA control chart is evaluated based on two types of linear time series data. Furthermore, the proposed control chart draws a correct decision when the drinking water manufacturing process needs to be improved. autocorrelated processes. To overcome this problem, researchers developed two general approaches to deal with autocorrelation in the process. First, fitting time series models to the autocorrelated data then monitoring the residuals using conventional control chart (Alwan & Roberts, 1988;Montgomery & Mastrangelo, 1991). Second, traditional control charts with modified control limits are used to monitor autocorrelated data (Lu & Reynolds, 1999;Vanbrackle & Reynolds, 1997). First approach seems to work better in high level of autocorrelation while second approach will work better in low-to-moderate level of autocorrelation (Yashchin, 1993).
Some researchers investigated the effect of autocorrelation on the control chart performance. Serially dependent data may lead to incorrect out-of-control signal and reduce the effectiveness of a control chart (Noorossana & Vaghefi, 2006). If data follow autoregressive moving average model with order (1,1), exponentially weighted moving average (EWMA) control chart will detect shift more quickly than Shewhart control chart (Wardell, Moskowitz, & Plante, 1994). Jarrett and Pan (2007) evaluated mean shift effect on the performance of multivariate control chart based on the residual of vector autoregressive (VAR) model. For small autocorrelation coefficient, Hotelling's T 2 control chart based on the residual of VAR model is effective to detect small shift of mean .
One of the well-known performance indicators to evaluate and to compare the effectiveness of various control charts is average run length (ARL). However, the distribution of run length is affected by serial correlation (Bagshaw & Johnson, 1975;Johnson & Bagshaw, 1974). Böhm and Hackl (1996) derived close form approximation of in-control ARL to evaluate the effect of autocorrelation on the performance of cumulative sum (CUSUM) control chart. Both Shewhart and CUSUM charts become less sensitive to detect small shift in the mean when applied to the residual of autocorrelated process (Runger, Wlllemain, & Prabhu, 1995). Kramer and Schmid (1997) evaluated the performance of residual-based multivariate exponentially weighted moving average (MEWMA) control chart for first order of bivariate VAR model. Śliwa and Schmid (2005) found that for monitoring cross-covariance purpose, residual-based MEWMA control chart has better performance than modified MEWMA control chart.
Recently, some researchers also use machine learning to improve the performance of residualbased multivariate control chart. Khediri, Weihs, and Limam (2010) studied the performances of residual-based T 2 , CUSUM of T 2 (COT), and T 2 of EWMA (EWMAT) control chart. Khediri et al. (2010) proved that multivariate control charts based on the residual of support vector regression (SVR) are more effective than multivariate control charts based on the residual of Artificial Neural Network (ANN). The SVR-based MCUSUM control chart has better performance than MCUSUM control charts based on the residual of ANN and VAR (Issam & Mohamed, 2008).
SVR is one of the machine learnings that provides global optimal solution. Meanwhile, least square SVR (LS-SVR) is least square version of standard SVR which replaces quadratic programming problem with linear programming problem (Vapnik, 1998(Vapnik, , 1995 and changes the inequality constraints formula by the equality ones (Suykens, Gestel, Brabanter, Moor, & Vandewalle, 2002;Suykens & Vandewalle, 1999). The linear equation in LS-SVR is simple to solve and good in computational time saving. When the variables of interest are more than one, the issue becomes multioutput case. Xu, An, Qiao, Zhu, and Li (2013) introduced multioutput LS-SVR (MLS-SVR) in order to cover the Hierarchical Bayes (Heskes, 2000) intuition, where each output is permitted to have different slope function. The MLS-SVR has ability to map the multivariate input space to the multivariate output space.  proposed MLS-SVR-based MEWMA control chart and evaluated its ability to detect the presence of both additive and innovative outliers in the simulation data which follow vector autoregressive and moving average model with order (1,1). However, the performance of MLS-SVR-based MEWMA control chart based on ARL criteria has not been investigated. Therefore, the aim of this research is to evaluate the performance of MLS-SVR-based MEWMA control chart  using ARL criteria. This research calculates the ARLs using simulation studies for several combinations of smoothing parameter and multivariate linear time series data, including VAR (1) and VMA (1) model. Moreover, this research is also aimed to provide an illustrative example of MLS-SVR-based MEWMA control chart for monitoring real autocorrelated data.
The rest of this paper is organized as follows. Section 2 describes MLS-SVR-based MEWMA control chart and the way the chart is constructed. ARL algorithm is displayed at Section 3. In addition, Section 4 presents the performance of MLS-SVR-based MEWMA control chart for VAR (1) and VMA (1) data. Application of the proposed method to monitor water quality data is displayed at Section 5. Finally, some concluding remarks and future research opportunity are given in Section 6.

MLS-SVR algorithm
The standard formula of LS-SVR trains the different output in multioutput case separately. Xu et al. (2013) proposed MLS-SVR algorithm to cope the potential cross-relatedness among different output such that the relation between outputs can be captured by learning all outputs simultaneously. Let Y ¼ ½y ij 2 R nÂ, are the observable outputs, where i ¼ 1; 2; . . . ; n, j ¼ 1; 2; . . . ; ,, i denotes number of samples and j for number of outputs. Suppose ðx 1 ; y 1 Þ; ðx 2 ; y 2 Þ; . . . ; ðx n ; y n Þ È É are the specific independent and identically distributed samples, where x i 2 R d , d is dimension of input, and y i 2 R , .
In order to cover the Hierarchical Bayes intuition, the parameter w j can be written as w j ¼ w 0 þ v j . The parameter w 0 describes mean vector whereas the small values of vector v j indicate that different outputs are similar to each other. In the other words, mean vector w 0 carries commonality information while vector v j carries specialty information. Let φ : R d ! R h is a mapping function to higher dimensional Hilbert space with h dimension. All MLS-SVR parameters are assumed to be associated with φðxÞ so that parameter w j 2 R h ðj 2 N , Þ, vector v j 2 R h ðj 2 N , Þ, and mean vector w 0 2 R h .
The objective of MLS-SVR algorithm is to find the function f ðxÞ ¼ φðxÞ T W þ ðbÞ T which has most deviation from the actual observable outputs and at the same time is as flat as possible. On the other words, the errors are acceptable as long as the errors are less than or equal to the deviation. The flatness of function f ðxÞ can be obtained by looking for the small value of parameter W. This condition is equivalent with minimizing parameter w T 0 w 0 as well as trace V T V . Therefore, estimating vector Þ 2 R , simultaneously could be attained by minimizing the following objective function with constraints (Xu et al., 2013): (1) meter, and γ 0 ; γ 00 2 R þ are regularized parameter.
The optimization problem has the assumption that it should be feasible, the function f ðxÞ actually exists for all approximates samples ðx 1 ; y 1 Þ; ðx 2 ; y 2 Þ; . . . ; ðx n ; y n Þ È É with precision equal to an acceptable deviation. A matrix Ξ containing slack variable is involved in optimization problem (1) in order to cope the infeasible constraints. In addition, the regularized parameters γ 0 and γ 00 describe the trade-off between the flatness of function f ðxÞ and the amount up to which errors larger than deviation are tolerated.
A set of dual variables is introduced in order to construct a Lagrange function from the objective function along with its constraints. The Lagrange function from the optimization problem (1) can be written as follows: where A ¼ α 1 ; α 2 ; . . . ; α , À Á T 2 R nÂ, contains Lagrange multiplier. According to Karush-Kuhn-Tucker method, a set of linear equation in Equation (2) yields an optimal condition.
Following the idea of LS-SVR, the linear equation system (3) is obtained by eliminating WandΞ matrices from Equation (2).
Lagrange multiplier, and output variables are defined However, solving linear equation system (3) directly is difficult because it is not positively definite. For this reason, linear equation system (3) is reformulated into the following one: It can be approved that G is a positive definite matrix. Thus, the solutions of linear equation system (4) can be obtained using following steps: (1) Solve # and ν from M# ¼ N and Mν ¼ y.
The solutions of the linear equation system (4) can be obtained by solving two linear equations which have same positive definite matrix M. Singular value decomposition is an alternative method to solve those linear equations. Supposed thatα ¼α andb are the solutions of linear equation system (4), then the decision function of MLS-SVR can be written as follows: Grid search method (Hsu, Chang, & Lin, 2016) is used to identify the proper hyper-parameters of MLS-SVR. For all possible combinations of kernel function parameter, σ 2 2 À15 ; 2 À13 ; . . . ; 2 3 È É , as well as regularized parameters, γ 00 2 2 À10 ; 2 À8 ; . . . ; 2 10 È É and γ 0 2 2 À5 ; 2 À3 ; . . . ; 2 15 È É ; the optimal pair of hyper-parameters σ; γ 0 ; and γ " is selected based on minimum criterion from average of mean square error (MSE) over each output. An evolutionary algorithm as proposed by Härdle, Prastyo, and Hafner (2014) can also be employed to optimize SVR parameters. This might become useful for future work.
Lowry, Woodall, Champ, and Rigdon (1992) introduced MEWMA control chart to monitor small change in the mean vector of independent observation. If the observation violates the independent assumption, then residual-based MEWMA control chart is employed. Supposed that e i ¼ ðe i1 ; e i2 ; . . . ; e ij ; . . . ; e i, Þ is vector of residual from ith observation. This vector of residual is transformed to Z i as follows: where λ is smoothing parameter with 0 < λ < 1 and initial value Z 0 ¼ 0. The transformed value Z i in Equation (7) is used to construct following MLS-SVR-based MEWMA statistic: where inverse covariance matrix of Z i is calculated by In addition, AE is covariance matrix of observation whereas T 2 i in Equation (8) denotes MLS-SVR-based MEWMA statistic. The process is said to be in-control if T 2 i statistic is not greater than the upper control limit (UCL), H, which refers to Prabhu and Runger (1997). Under in-control condition, the residuals of MLS-SVR model using the proper inputs selected using Equation (6) and the optimal parameters obtained as in Equation (5) are assumed to follow independent observation. Thus, the UCL of MLS-SVR residual-based MEWMA control chart is equal to the UCL of MEWMA control chart.

ARL algorithm
ARL is one of the well-known methods to measure and to compare the control chart performance. Run length (RL) is defined as the number of observations until the first observation is detected outside the control limits. The control chart performance could be assessed by two types of ARLs that are ARL 0 and ARL 1 . The ARL 0 or in-control ARL is defined as the expected number of observations before an out-of control observation is detected when the process is actually in-control. The ARL 1 or out-of control ARL is the expected number of samples before an out of control signal is received when the process is actually shifted to an out-of-control state (Montgomery, 2009). For a fixed value of ARL 0 , a chart is considered to be more effective than other charts if it has a smaller ARL 1 .
There are some general methods to calculate the ARL of a control chart such as simulation, integral equation, and Markov chain approximation. Those methods had been developed by several researchers to calculate the ARL of MEWMA control chart. Lowry et al. (1992) and Love (2000a, 2000b) calculated the ARL of MEWMA control chart using simulation. Rigdon (1995) and Bodden and Rigdon (1999) acquainted integral equation approximation to calculate the ARL of MEWMA chart. Furthermore, a Markov chain approximation for calculating the ARL of MEWMA chart is proposed by Runger and Prabhu (1996) as well as Molnau et al. (2001).
The ARLs of a control chart for time series data have a complex calculation. Śliwa and Schmid (2005) confirmed that even for univariate case, there is no explicit formula to calculate the ARL of a control chart for time series data. Kramer and Schmid (1997) and Śliwa and Schmid (2005) utilized extensive Monte Carlo simulation to determine the ARL of MEWMA control chart for time series data. In this research, the ARL of MLS-SVR-based MEWMA control chart is calculated using simulation study based on Algorithm 1.
Algorithm 1. The ARL of MLS-SVR-based MEWMA Control Chart 1. Specify the number of variables , and smoothing parameter λ. 2. Define the UCL H for specified parameters given in step 1. 3. Specify the parameter of multivariate linear time series data Φ p or Θ q . 4. Generate multivariate linear time series data using parameter Φ p or Θ q given in step 3 and the residuals which satisfy multivariate normal distribution with specified mean vector μ a ¼ 0 and covariance matrix AE ¼ I. 5. Train generated data resulted from step 4 using MLS-SVR algorithm, where the inputs are selected using Equation (6) and the optimal parameters along with the optimal hyper-parameters are obtained as in Equation (5). 6. For C 0 repetitions, follow these steps: a. Generate n samples of multivariate time series data using parameter Φ p or Θ q given in step 3 and the residuals which have mean vector μ a ¼ 0 and covariance matrix AE ¼ I. b. Train the data which are generated in step 6(a) using inputs, parameters, and hyperparameters obtained from step 5. c. Calculate statistics T 2 i ; i ¼ 1; 2; . . . ; n using Equation (8). d. Calculate the RL, number of samples until finding the first statistic T 2 i that greater than H. 7. Calculate ARL 0 from the average of RL over C 0 repetitions. 8. For C 00 repetitions, where μ 0 a ¼ μ a þ mean vector shift, follow these steps: a. Repeat step 6 which multivariate time series data are generated with residuals that have mean vector μ 0 a and covariance matrix AE ¼ I: b. Calculate ARL 1 from the average of RL over C 00 repetitions. 9. Repeat step 1 until 7 for different ,, λ and the associated H. 10. Plot the ARL 0 and ARL 1 over the mean shift.

Performance evaluation of the proposed control chart
This section is aimed to investigate the performance of the proposed control chart using ARL. For this purpose, simulation studies are designed to generate multivariate linear time series data which follow VAR (1) or VMA (1) model using this mathematical expression: VMAð1Þ where a i is white noise residual which follows multivariate normal distribution with zero mean vector, μ a ¼ 0 , , and covariance matrix AE ¼ I , . The performance of the proposed control chart is evaluated for , ¼ 2 and , ¼ 3 quality characteristics. For , ¼ 2 quality characteristics, the mean vector is set to be μ ¼ μ 1 μ 2 ½ T ¼ 5 10 ½ T . Otherwise, the mean vector is equal to The parameters of VAR (1) and VMA (1) used in the simulation studies are shown in Table 1.
The simulation studies to calculate the ARL of the proposed control chart are carried out based on Algorithm 1 with n ¼ 1000 samples, C 0 ¼ 1,000 repetitions and C 00 ¼ 14 repetitions. These simulation studies utilize the UCL which is chosen for in-control ARL equal to 200. The in-control ARL of the proposed control chart is calculated using generated in-control data which follow VAR (1) or VMA (1) model, whereas the out-of-control ARL is obtained by adding mean shift 0.2 to each residual of in-control data, i.e. μ 0 a ¼ μ a þ mean vector shift. The multivariate linear time series data which follow VAR (1) model are modelled using MLS-SVR algorithm, where Y 1ðiÀ1Þ and Y 2ðiÀ1Þ are selected as inputs, whereas the generated data which follow VMA (1) model are modelled using Y 1ðiÀ1Þ ; Y 1ðiÀ2Þ ; . . . ; Y 1ðiÀ10Þ and Y 1ðiÀ1Þ ; Y 2ðiÀ2Þ ; . . . ; Y 2ðiÀ10Þ as the inputs in MLS-SVR algorithm.
The selection of the input is based on the significance lag of partial autocorrelation function. The radial basis function (RBF) kernel function is employed in MLS-SVR modelling.
Tables 2 and 3 present the ARL of the proposed control chart for VAR (1) data with , ¼ 2 and , ¼ 3, respectively. For various smoothing parameters λ; MLS-SVR-based MEWMA control chart yields stable ARL 0 at about 200. The ARL 1 is increasing when the magnitude of smoothing parameter λ is rising. In other words, the probability of the proposed control chart to detect a specified shift in the mean is decreasing as the increasing of the magnitude of smoothing parameter λ. Moreover, the ARL 1 is decreasing as the increasing of the mean shift. This indicates that the larger the mean shift, the faster ability of the proposed control chart to detect those actual shift.   The ARLs of the proposed control chart for VMA (1) data with , ¼ 2 and , ¼ 3 are shown by Tables 4 and 5, respectively. The proposed control chart tends to yield ARL 0 at about 200. The ARLs of the proposed control chart for VMA (1) have similar pattern with those for VAR (1) data. If the mean of each variable has shifted from 0 to 0.2, then the ARL of the proposed control chart with smoothing parameter λ = 0.2 for monitoring VMA (1) with , = 2 variables decreases from 200 to about 150. However, the ARL of the proposed control chart for monitoring VMA (1) data with , = 3 variables diminishes from 196 to about 139 to detect that shift.
The ARLs comparison for VAR (1) and VMA (1) data is exhibited at Figure 1. The smaller the smoothing parameter λ, the faster ability of the proposed control chart to detect an actual mean shift of a process. If the mean shift of each variable has shifted from 0 to 0.2, then the ARL of the proposed control chart for monitoring VAR (1) data with λ = 0.05 diminishes to about 50. Otherwise, the ARL of the proposed control chart for monitoring VMA (1) data with λ = 0.05 decreases to about 100 to detect the same mean shift of a process. This indicates that the proposed control chart for monitoring VAR (1) data is more sensitive than that for monitoring VMA (1) data. Moreover, the plots of ARL for monitoring VMA (1) data with λ = 0.05 are significantly different with those with λ = 0.8. Hence, the performance of the proposed control chart to detect mean shift of VMA (1) data depends much on smoothing parameter λ. On the contrary, the smoothing parameter λ does not play a significant role to the ARL of the proposed control chart for VAR (1) data.

Application
In this section, the proposed control chart is applied to monitor water quality data in Surabaya, Indonesia. Two important quality characteristics in drinking water manufacturing process are water turbidity and chlorine residual. Turbidity is measured by the concentration of dissolution and the presence of particles in a liquid using Nephelometric Turbidity Units (NTU). Chlorination or affixing chlorine in to contaminated water is intended primarily for microbes killing, where chlorine residual is measured in ppm unit. Drinking water would be safe from bacteria if it has minimum chlorine residual of 0.2 ppm. However, excessive chlorine addition would affect the smell and the  Khusna et al., Cogent Engineering (2018), 5: 1531456 https://doi.org/10.1080/23311916.2018 taste of water. Figure 2 displays time series plots of water turbidity and chlorine residual for Phase I monitoring process. Those variables are measured hourly from 19 August 2016 until 29 August 2016.  Water quality data displayed at Figure 2 are modelled using MLS-SVR algorithm, where Y 1ðiÀ1Þ ; Y 1ðiÀ2Þ ; Y 1ðiÀ7Þ ; Y 1ðiÀ8Þ and Y 2ðiÀ1Þ ; Y 2ðiÀ3Þ ; Y 2ðiÀ9Þ are selected as input variables. These inputs are chosen by considering the significant lags of partial autocorrelation function. MLS-SVR modelling using RBF kernel function yields the optimal combination of hyper parameters γ 0 ¼ 2 À5 , γ 00 ¼ 2 À8 , and σ ¼ 2 3 . The residuals of MLS-SVR model satisfy multivariate normal distribution and white noise condition with minimum MSE equal to 0.000605. Furthermore, these residuals are monitored using MEWMA control chart with significance level α ¼ 0:005 and smoothing parameter λ ¼ 0:2 as displayed at Figure 3. It can be shown that MLS-SVR-based MEWMA control chart for Phase I monitoring process fulfils incontrol condition. Thus, the optimal combination of hyper parameters and the control limit can be applied in Phase II monitoring process.   Time series plots of water turbidity and chlorine residual data for Phase II monitoring process are displayed at Figure 4. Both variables are measured in 5 days starting from 30 August 2016. It can be known that the time series plots for first four days in Phase II monitoring process follow the same pattern as the time series plots of water quality data in Phase I. However, time series plots for last day in Phase II monitoring process, starting from 10.00 AM, show unusual pattern. Time series plot of water turbidity indicates increasing pattern while time series plot of chlorine residual exhibits steep shift pattern.
The water quality data displayed at Figure 4 are then modelled using MLS-SVR algorithm using the optimal combination of hyper parameters obtained from Phase I monitoring process. The MLS-SVRbased MEWMA control chart for Phase II monitoring process is shown at Figure 5. The residuals resulted from Phase II monitoring process do not satisfy white noise condition. Consequently, it leads to many false alarms and tight control limit. The increasing pattern far away from control limit in last eight observations indicates a serious problem in water quality data. Hence, monitoring water quality data using MLS-SVR-based MEWMA control chart becomes an early warning such that water manufacturing process needs to be improved. Looking for the assignable causes is needed in order to improve the manufacturing process. In 3 September 2016, leak pipeline is found in one of the water distribution area. Moreover, pipeline maintenance process and flow meter installation are started at 08.00 PM. As a result, both water turbidity and chlorine residual plots indicate unusual pattern.
The MLS-SVR-based MEWMA control chart for water quality data needs to be updated in order to fit the next Phase I monitoring process. This step could be solved by replacing original unusual pattern with predicted value as displayed at Figure 6. Data displayed at Figure 6 are then modelled using MLS-SVR algorithm as testing data. Furthermore, white noise residuals resulted from this  process are monitored using MEWMA control chart. Figure 7 shows MLS-SVR-based MEWMA control chart for updated water quality data. All of the observations fulfil in-control condition such that this control chart could be fitted as the next Phase I monitoring process.

Conclusions and future research
The performance of MLS-SVR-based MEWMA control chart for monitoring multivariate linear time series data which follow VAR (1) and VMA (1) model is investigated using ARL criteria. The smaller the smoothing parameter, the more sensitive of the proposed control chart to detect an actual mean shift of a process. For a certain value of smoothing parameter, the proposed control chart for monitoring VAR (1) data is more sensitive than that for monitoring VMA (1) data. The ARLs of the proposed control chart for monitoring VMA (1) data depend much on the smoothing parameter. In addition, monitoring hourly water quality data using MLS-SVR-based MEWMA control chart becomes an early warning to improve drinking water manufacturing process if the actual assignable causes are occurred. Developing the theoretical ARL formula for MLS-SVR-based MEWMA control chart is one of the open future research topics. Bootstrap resampling method  and kernel density estimation (Ahsan, Mashuri, Kuswanto, Prastyo, & Khusna, 2018a, 2018b are two computational techniques which can be used to develop the more accurate control limit of the proposed chart. Furthermore, the mixed monitoring scheme as demonstrated by Ahsan, Mashuri, Kuswanto, Prastyo, and Khusna (2018b) can also be considered for future development of this proposed chart.