Behaviour Analysis, Complexity and Control of Networked Dynamical SystemsView this Special Issue
SPA-Based Modified Local Reachability Density Ratio wSVDD for Nonlinear Multimode Process Monitoring
Many industrial processes are operated in multiple modes due to different manufacturing strategies. Multimodality of process data is often accompanied with nonlinear and non-Gaussian characteristics, which makes data-driven monitoring more complicated. In this paper, statistics pattern analysis (SPA) is introduced to extract low- and high-order statistics from raw process data. Support vector data description (SVDD), which can deal with nonlinear and non-Gaussian problems, is applied to monitor multimode process in this paper. To improve detection performance of SVDD for training multimode data with outliers, modified local reachability density ratio (mLRDR) is proposed as a weight factor to be embedded in the weighted-SVDD (wSVDD) model, in which the local neighbors in terms of both space and time are considered. Finally, the effectiveness and superiority of our proposed method are demonstrated by the Tennessee-Eastman (TE) process and wastewater treatment process (WWTP).
The operation conditions of industrial processes will inevitably change with diverse customer requirements, set-points variation, and different intrinsic features, which leads to multiple modes [1‒3]. Because of the complexity of multimodal process, it is difficult to obtain satisfactory monitoring results. In the past decades, effective monitoring of multimode processes has attracted a lot of attention.
The data-driven monitoring methods of multimode process were recently reviewed extensively by Quinones-Grueiro et al. . They indicated that multimodel scheme is one of the widely applied monitoring methods, in which the identification of different modes should be performed in advance. For instance, Du et al.  applied K-means to address the multimode clustering problem. Similarly, Zhang and Zhao  clustered multimode data into different clusters through the kernel fuzzy c-means (KFCM) algorithm. In these methods, the cluster number should be confirmed in advance. As an alternative, the moving window-based methods, for example, recursive local outlier factor (LOF) , measure the similarity between windows with the spatial and temporal information of the features to identify the clusters. After mode identification, each mode is monitored by such traditional methods as principal component analysis/independent component analysis (PCA/ICA) and their extensions and so forth. [8‒10]. Recently, several efforts have been made on applying the multiple models with Bayesian fusion to multimode monitoring. Specifically, the Gaussian Mixture Model (GMM) was deployed to accomplish fault diagnosis by considering multiple models simultaneously [11, 12].
Compared with the multimodel scheme, single-model-based method simplifies the modeling and monitoring procedure. Hidden Markov model- (HMM-) based methods were proposed in [13, 14], whereas their parameter setting is still a complex task. Zhu et al.  proposed a recursive mixture factor analyzer for multimode time-variant process modeling and monitoring. Considering that there are interconnections among different modes, a strategy was developed to grasp common characteristics in them to construct a monitoring model. Hwang and Han  proposed the super principal component analysis (PCA) method, which can be considered as the first attempt based on this idea. Then, Zhang et al.  built the common and individual monitoring model by PCA-based and partial least squares- (PLS-) based scheme. Besides, Zhang et al.  utilized the common basis vectors to extract the common features among multimode data. For nearest-neighbor approaches, different modes should be standardized . For example, Ma et al.  presented a neighborhood-based global coordination (NBGC) method by aligning the local models into a global one for multimode process monitoring.
Among the single-model-based methods, SVDD model and its extensions are widely applied. Developed from support vector machine , SVDD is proposed based on the idea of one-class classification. SVDD projects the data into high-dimensional feature space without the requirement of Gaussian distribution . Therefore, the SVDD method can deal with both nonlinear and non-Gaussian data . However, SVDD is insensible to a very small fraction of outlier, and the trained hypersphere will sometimes tend to encircle outliers that are far away from the normal dataset. Besides, the normal samples with different density have no different effects on SVDD modeling. In the multimode training dataset, the sampling number and data density vary dramatically with different modes. Thus, SVDD is not suitable to be directly applied to multimode monitoring. The wSVDD method was developed to solve this problem, in which the weight factor is usually determined by nearest-neighbor information. To reduce the impact of outliers on modeling in single-mode process, Chen et al.  proposed robust-SVDD by introducing cutoff-distance-based local density of each data sample and the -insensitive loss function with negative samples. Wang and Lan  used the SD outlyingness to assign lower weight values to outliers. For multimode process, Zhao et al. utilized the weighted mean and standard deviation of each sample’s neighbors to standardize the dataset and applied weighted local standardization (WLS) strategy to wSVDD . Li et al.  proposed a local density ratio weighted support vector data description (LDR-wSVDD) to maintain the monitoring efficiency of a single hypersphere model. Nevertheless, the spatial information is still not fully mined, and the temporal information is not employed in the construction of the weight for wSVDD, which may lead to unsatisfactory performance of the existing methods. Consequently, wSVDD still needs to be further studied.
Recently, Breunig et al.  defined local reachability density (LRD) and indicated that the LRD of a sample is related to the local information of not only itself but also its nearest neighbors. By containing more comprehensive density information, LRD has the potential to be transformed into a new weight factor for wSVDD. However, LRD varies drastically with different modes; thus, it does not yet reflect their local information equally. In addition, conventional LRD only considers spatial information and does not take temporal information into account.
Moreover, the data in different modes may overlap each other in actual industrial process, so that the original data sometimes fails to embody the unique characteristics of certain modes. Statistics pattern analysis (SPA) was firstly introduced by Wang and He , which is able to mine underlying statistic features of data and has been applied in process monitoring . Zhu and Gu integrated SPA into local kernel principal component analysis (LKPCA) to enhance fault detection performance .
In this paper, SPA is applied to obtain statistical characteristics for fully mining the underlying information of multimode dataset firstly. Then, wSVDD is selected as a monitoring model, and modified local reachability density ratio (mLRDR) is proposed as a new weight for wSVDD to remove the multimodality of process and reduce the impact of outliers. The local information in terms of both space and time is considered in mLRDR calculation to further enhance monitoring performance. Ultimately, a novel scheme called SPA-based modified local reachability density ratio weighted support vector data description (SmLRDR-wSVDD) for nonlinear multimode monitoring is established.
The remainder of the paper is organized as follows: in Section 2, the proposed nonlinear multimode process monitoring scheme is elaborated in detail; in Section 3, the results and discussion of TE process and WWTP are presented; finally, conclusions are given in Section 4.
In this section, SmLRDR-wSVDD method is proposed to monitor nonlinear multimode process, in which statistics pattern dataset is formed by SPA, and a new weight factor named mLRDR is proposed based on LRD for wSVDD monitoring method.
2.1. Construction of Statistics Pattern
Although specific characteristics of different modes may not be embodied from the original data, they can be expressed in terms of the statistics. In the SPA framework , different statistics patterns for observed variables are selected to capture the dominant process characteristics such as dynamics and nonlinearity and are used as the modeling and monitoring object in the proposed method.
Suppose that is the original measured dataset with samples and variables, and a window of samples is denoted aswhere is the window width and is the time index. Generally, statistics pattern (SP) consists of three groups of statistics: first-order, second-order, and high-order statistics, which can be expressed aswhere denotes the first-order statistics, that is, the mean, and its elements are calculated from the data in a window:where denotes the second-order statistics, which includes variance (), correlation (), autocorrelation (), and cross-correlation () of different variables in the window aswith denoting the time lag between the variables; denotes the high-order statistics including skewness () and kurtosis () to measure the degree of nonlinearity and quantify the non-Gaussianity of the process variables  as
The SP is constructed by stacking selected statistics by equations (1)‒(9) into a row vector. In the paper, three kinds of statistics , , and are selected from one-order, two-order, and high-order statistics to construct statistics matrix with , , and .
2.2. mLRDR Weighted SVDD
The main idea of SVDD is to project data into high-dimensional space for constructing a minimum hypersphere. The training dataset is obtained according to equation (2), and it is normalized as . In the procedure of SVDD modeling, is firstly mapped from the original space to a higher feature space by a nonlinear transformation function . For conventional SVDD, all the training data has the same impact on the model construction, which makes it insensible to outliers and data density. Therefore, wSVDD is adopted by introducing a weight for each training data, and the corresponding hypersphere can be calculated by solving the following primal optimization problem:where and are the center and the radius of hypersphere, respectively; trade-off parameter is introduced to make the hypersphere as small as possible while preventing the misclassification of the samples, in which denotes the percentage of outliers permitted in the training set ; is slack variable that allows outliers of the hypersphere; and is weight of the data, where smaller indicates more possibility of an outlier.
For multimode processes, the weight factor of wSVDD should fully reflect the local distribution characteristics of each piece of training data. Density is an indicator reflecting the data distribution, which is widely used in the design of weight factor. Local reachability density (LRD)  of a data point depends on its distance to not only the point’s neighbors but also its neighbors’ neighbors.
Besides the space density, the current data is also related to the previous and the next points in time sequence, and monitoring performance will be improved by taking both the temporal information and spatial information into account . Consequently, the modified LRD is proposed in this paper by considering the local neighbors in both time and space to comprehensively extract the local features of the data.
Firstly, all K-nearest neighbors of are selected to construct a set , which satisfies , where and are the numbers of neighbors in space and time, respectively.
Suppose that is the Euclidean distance from to its -nearest neighbor. Then the distance between and (k = 1, 2, …, K) is computed, where is the -nearest neighbor of . The reachability distance of with respect to is defined as
Thus, mLRD is constructed as
The main property of weight factor for multimode process should equalize different mode and reduce the influence of outlier in modeling. According to equation (12), if distribution characteristics vary greatly with different modes, mLRD in different modes will change a lot; and some outliers will be hardly separated from normal data; that is, the local density of normal data with lower global density may be similar to that of the outliers around the data with higher global density. Accordingly, to weaken the multimodality feature and widen the gap between outlier and the normal data, a modified local reachability density ratio (mLRDR) is developed to be the weight factor of wSVDD as
It is obvious that is kept within .
According to equation (13), if is a normal point, its mLRD is close to its neighbors and the corresponding mLRDR approaches 1. If is an outlier, its mLRD is significantly smaller than that of its normal neighbors and mLRDR is close to 0, which indicates that outliers will have trivial influence on the determination of the hypersphere.
Compared with the existing weight factors such as LDR , mLRDR contains more comprehensive spatial and additional temporal information, and it can distinguish normal data with lower density and the outliers around the data points with higher density data; thus it will improve the monitoring performance of wSVDD.
For wSVDD modeling, equation (10) can be transformed into a dual-optimization problem by introducing Lagrange multipliers as follows:where is kernel function and can be determined via detecting the “tightness” of the decision boundaries . After solving equation (14), only objects with are called the support vectors (SVs), and their subscripts set is . Hence, the center and radius of the hypersphere arerespectively.
When new samples are obtained for a window, SPs are calculated by equation (2) and normalized as , and the distance between and the center can be calculated by
If , the system is considered in normal state. Otherwise, a fault occurs.
2.3. The Procedure of SmLRDR-wSVDD Scheme
In this subsection, the proposed SmLRDR-wSVDD scheme is illustrated in Figure 1. The specific offline modeling and online monitoring steps are listed as follows: (i)Step 1: adopt SPA on a multimode dataset to obtain statistics matrix by equations (2)–(9).(ii)Step 2: normalize statistics matrix and obtain .(iii)Step 3: compute the mLRDR as a weight factor by equations (11)–(13).(iv)Step 4: set up weighted-SVDD model with the weight , solve equation (14), and obtain the center and the radius of the hypersphere by equations (15) and (16), respectively. (i)Step 1: for a current window dataset , statistics vector is calculated by equations (2)–(9).(ii)Step 2: normalize statistics vector to obtain .(iii)Step 3: calculate the distance between and the center of the hypersphere by (17fd17).(iv)Step 4: if , is considered to be normal; otherwise, it is in a faulty state.
In this paper, if s of 3successive windows are over the control limit , the system is considered in faulty condition, and the last sample in the first faulty window is determined as the beginning of the fault.
It must be noted that the identification of new normal modes is not considered in the proposed method. Although a new normal mode will be detected as a fault by SmLRDR-wSVDD, it may be redefined by analyzing the stationarity of sampling data after enough data are collected.
In summary, the proposed method adopts SPA to mine more features from different-order statistics of the original data compared with other existing SVDD-based methods. Besides, a new weight factor is proposed for wSVDD modeling, in which the additional neighbors’ local spatial information and temporal information are considered, which makes the proposed method more sensitive to outliers and density compared with other methods.
3. Case Study
In this section, the proposed SmLRDR-wSVDD method is applied on Tennessee-Eastman (TE) benchmark process and wastewater treatment process (WWTP). In both cases, the proposed method is compared with the existing methods including KPCA , conventional SVDD, and LDR-wSVDD .
3.1. Tennessee-Eastman Benchmark Process
Tennessee-Eastman process (TE process) was put forward in 1993 by Downs and Vogel , which has been widely adopted for scientific research . There are five units, reactor, condenser, compressor, stripper, and separator, in this nonlinear process. According to different mass ratios ( and are two products), 6 steady modes can be obtained as listed in Table 1.
The control scheme and normal operations were provided by Ricker  and the corresponding simulation platform can be downloaded from http://depts.washington.edu/control/LARRY/TE/download.html. Mode 1 and Mode 3 are simulated to generate modeling and testing datasets. There are 53 variables: 12 manipulated variables, 22 continuous process variables, and 19 composition measurements. Since the steady-state values of the recycle value and steam value in Mode 1 and the agitator rates in both modes do not change, they are not included in the monitored variables, and the remaining 9 manipulated variables and 22 easily measured continuous variables are chosen in this case. In addition, TE process simulation platform contains 20 faults. Excluding faults 3, 9, and 15, which are slight disturbances, and unknown faults 16–20, 12 faults given in Table 2 are employed; and the sampling interval is set as 3 minutes.
For model training, 800 normal samples from Mode 1 and 200 normal samples from Mode 3 are collected, and 20 outliers are mixed with them. For online monitoring, each test dataset contains 1000 samples, in which the first 200 samples are normal, and the fault occurs at the sample. The confidence level of KPCA is chosen as . For the other three algorithms, the Gaussian kernel width parameter is determined by the method in , and is selected for the penalty parameter . In the proposed method, the window width is assigned as 10, and the numbers of neighbors in space and time and are chosen as 8 and 2, respectively. The above-mentioned parameters are the same in the two cases.
The four methods, that is, KPCA, SVDD, LDR-wSVDD, and the proposed SmLRDR-wSVDD, are applied and compared on the TE process. The missed alarm rates (MARs)  for each fault type of Mode 1 and Mode 3 are shown in Table 3, in which the best results are marked in bold, and the false alarm rates (FAR) are shown in Table 4. Specifically, MARs of KPCA on faults 4, 5, 6, 7, 10, 11, 12, and 14 of Mode 1 and faults 4, 5, 6, 7, 10, 11, and 14 of Mode 3 are all more than 0.5, indicating that a single KPCA model is not suitable to monitor multimode process. Among the four methods, the proposed SmLRDR-wSVDD method achieves the lowest MARs for almost all the faults except for fault 5 of Mode 1 and fault 2 of Mode 3. For example, the proposed method greatly decreases the MARs of fault 10 for both Mode 1 and Mode 3 and fault 12 for Mode 1 by over 9% compared to LDR-wSVDD. It must be noted that all methods have high MAR for fault 5 of Mode 1, since it is actually a slight disturbance.
Further, the detailed monitoring results of four methods for different faults of two modes are shown in Figures 2–5. For fault 11 of Mode 1, the SPE of KPCA and Dist of SVDD fluctuate violently and some of them are under the limits shown in Figures 2(a) and 2(b); as illustrated in Figure 2(c), Dist fluctuates slightly, and some fall below the limit by LDR-wSVDD; meanwhile, with the proposed method shown in Figure 2(d), all the statistics are far above the control limit. Similarly, the proposed SmLRDR-wSVDD can separate normal condition and fault more obviously than the other 3 methods as shown in Figures 3 and 4. The monitoring result of fault 13 for Mode 3 is shown in Figure 5. For this fault, there are 48, 47, 27, and 9 samples of detection delay by KPCA, SVDD, LDR-wSVDD, and SmLRDR-wSVDD. Thus, one can infer that SmLRDR-wSVDD is more sensitive to the four faults than the other methods.
FARs of all methods are summarized in Table 4. For all faults of Mode 1, FARs of KPCA and LDR-wSVDD are within and , respectively, while those of conventional SVDD and SmLRDR-wSVDD are within . In Mode 3, FARs of SPE in KPCA are the largest for all faults, sometimes even larger than 0.05, and those of LDR-wSVDD for all faults range from 0.015 to 0.035. Note that all FARs of the proposed SmLRDR-wSVDD method are less than 0.02, though most of them are higher than those of conventional SVDD within 0.01.
Figure 6 shows LRD and the weight of each training sample based on LDR and the proposed mLRDR, respectively, in which blue, green, and red dots are from Mode 1, Mode 3, and outliers, respectively. In Figure 6(a), it is easy to find that the normal modes vary greatly in terms of LRD, and there is no distinctive gap between Mode 3 and outliers. As shown in Figures 6(b) and 6(c), LDRs and mLRDRs of samples in two normal modes distribute around 0.6 and 0.9 evenly. Compared with LDRs, mLRDRs of outliers are quite lower than those of normal data. Therefore, the decision boundary is not susceptible to outliers by the proposed method, which will improve its monitoring performance compared with other methods.
3.2. Wastewater Treatment Process (WWTP)
In this subsection, benchmark simulation model no. 1 (BSM1) of WWTP  is considered. This benchmark is a complicated nonlinear system, which consists of five compartments (two anoxic and three aerobic) in the biological reactor and a ten-layer clarifier. The layout of the plant is shown in Figure 7, and more detailed information about BSM1 can be found at http://www.benchmarkwwtp.org/.
By changing the flow rates and oxygen transfer coefficients in biological reactor, three different modes can be obtained . Without loss of generality, Mode 1 and Mode 2 are chosen in this case, and the corresponding parameters are shown in Table 5. The process contains 223 variables, 20 of them is related to the biological phenomena, and they are selected in this example as shown in Table 6.
There are two disturbance types in the benchmark simulation : external disturbances and internal disturbances. Specifically, dry weather, storm events, and prolonged rain are three weather situations, among which the first one is regarded as normal and the latter two are external disturbances, which can be detected by the influent characteristics. As for the internal disturbances, decreasing nitrification, decreasing settling velocity, nitrate sensor failure, and setpoint change of DO controller were mentioned in .
In this example, the samples are collected every 15 minutes. 1344 samples of dry weather of Mode 1 and Mode 2 shown in Table 7 are mixed with 20 additional outliers as training dataset. 5 simulation cases (Case 1–Case 5) are designed for the purpose of testing [40, 41], which are also illustrated in Table 7.
The four above-mentioned algorithms are performed for the 5 cases, and the parameters of them are the same as those in the TE process study. The MARs for all faults are shown in Table 8. For each case, the best result is marked in bold. All the four methods have low MARs in Cases 1 and 4. Meanwhile the MARs of KPCA are beyond 0.2 in Cases 2, 3, and 5, and those of conventional SVDD are more than 0.1 in Cases 2 and 5. Most MARs of LDR-wSVDD are below 0.05, except 0.1265 in Case 2. It must be noted that the MARs of the proposed method in Cases 1, 3, 4, and 5 are 0 and only 0.0164 in Case 2, all of which are the lowest.
Specifically, the monitoring results for Case 1 are shown in Figure 8. As shown in Figures 8(a)–8(c), the fault is detected with 6, 8, and 7 points’ delay by SPE of KPCA, the conventional SVDD, and LDR-wSVDD, respectively; meanwhile the proposed method can detect the fault immediately after its occurrence, as shown in Figure 8(d). For Case 2 shown in Figure 9, the fault occurs after the sample, and it changes slowly away from normal. From Figures 9(a) and 9(b), there are a lot of statistics falling below the control limits after the sample under KPCA and the conventional SVDD. As shown in Figure 9(c), though the weight factor, that is, LDR, is applied in wSVDD, many statistics are still under the limit at the beginning of the fault. As shown in Figure 9(d), the fault is detected quickly by the proposed method, and the following statistics are far above the control limit. Compared with other methods, SmLRDR-wSVDD is more effective for the faults with slow drift in the budding stage.
As for FARs shown in Table 9, the proposed SmLRDR-wSVDD method has the lowest FARs in Cases 2, 3, and 4. Although its FARs are slightly higher than those of conventional SVDD for Cases 1 and 5, both of them are fewer than 0.03, which is too trivial to be considered.
Similarly, LRD, LDR, and mLRDR of each training sample are visualized in Figure 10 to explain the superiority of our proposed method. Blue, green, and red dots are from Mode 1, Mode 2, and outliers, respectively. Obviously, LRD changes with different modes, and that of outliers is close to normal in Figure 10(a). Although LDRs can remove multimodality as shown in Figure 10(b), the LDRs of outliers and normal samples are not distinguished. As a comparison, mLRDRs of outliers are far away from those of normal as shown in Figure 10(c). It is demonstrated once again that the proposed method can deal with multiple modes evenly and reduce the influence of outliers on modeling.
In this paper, a nonlinear multimode process monitoring method named SmLRDR-wSVDD is proposed to further improve the monitoring performance. Specifically, the SPA is used first to extract low- and high-order statistics; then, weight-SVDD is adopted to build a monitoring model, in which a new weight factor mLRDR is proposed based on LRD. Compared with the traditional methods, more features that are difficult to display are mined from raw data by the proposed method, and more comprehensive local space and time information is integrated into the weight factor, which can remove the multimodality uniformly and greatly weaken the influence of outliers on the boundary of the hypersphere. The applications in the TE process and WWTP have demonstrated that the proposed method is effective in nonlinear multimode process monitoring with lower MAR/FAR and higher fault sensitivity compared with the other three existing methods.
The data in the manuscript are generated according to http://depts.washington.edu/control/LARRY/TE/download.html and http://www.benchmarkwwtp.org/.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
This work was supported by the National Natural Science Foundation of China under Grant 61873102.
D. Ha, U. Ahmed, H. Pyun, C.-J. Lee, K. H. Baek, and C. Han, “Multi-mode operation of principal component analysis with k-nearest neighbor algorithm to monitor compressors for liquefied natural gas mixed refrigerant processes,” Computers & Chemical Engineering, vol. 106, pp. 96–105, 2017.View at: Publisher Site | Google Scholar
J. Alex, L. Benedetti, J. Copp, et al., Benchmark Simulation Model no.1 (BSM1), pp. 1–62, 2008.