# SPA-Based Modified Local Reachability Density Ratio wSVDD for Nonlinear Multimode Process Monitoring

Authors: Zhaojing Wang; Weidong Yang; Hong Zhang; Ying Zheng
Journal: Complexity
 (2021)

---

## Abstract

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).


---

## Body

## 1. Introduction

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. [4]. 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. [5] applied K-means to address the multimode clustering problem. Similarly, Zhang and Zhao [6] 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) [7], 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. [15] 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 [16] proposed the super principal component analysis (PCA) method, which can be considered as the first attempt based on this idea. Then, Zhang et al. [17] built the common and individual monitoring model by PCA-based and partial least squares- (PLS-) based scheme. Besides, Zhang et al. [17] utilized the common basis vectors to extract the common features among multimode data. For nearest-neighbor approaches, different modes should be standardized [18]. For example, Ma et al. [19] 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 [20], 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 [21]. Therefore, the SVDD method can deal with both nonlinear and non-Gaussian data [22]. 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. [23] 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 [24] 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 [25]. Li et al. [26] 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. [27] 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 [28], which is able to mine underlying statistic features of data and has been applied in process monitoring [29]. Zhu and Gu integrated SPA into local kernel principal component analysis (LKPCA) to enhance fault detection performance [30].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 Section2, 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.


## 2. Methodology

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 [28], 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 thatX∈ℝn×m is the original measured dataset with n samples and m variables, and a window of samples is denoted as(1)Xth=xt,1h,xt,2h,…,xt,mh=x1t−h+1x2t−h+1…xmt−h+1x1t−h+2x2t−h+2…xmt−h+2⋮⋮⋮⋮x1tx2t…xmt,where h is the window width and t 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 as(2)S≡μ|V|Ξ,where μ denotes the first-order statistics, that is, the mean, and its elements are calculated from the data in a window:(3)μit=1h∑q=0h−1xit−q,where V denotes the second-order statistics, which includes variance (vi), correlation (ri,j), autocorrelation (ri), and cross-correlation (ri,j) of different variables in the window as(4)vit=1h∑q=0h−1xit−q−μit2,(5)ri,jt=1h∑q=0h−1xit−q−μitxjt−q−μitvitvjt,(6)ridt=1h−d∑q=dh−1xit−q−μitxit+d−q−μitvit,(7)ri,jdt=1h−d∑q=dh−1xit−q−μitxjt+d−q−μitvitvjt,with d denoting the time lag between the variables; Ξ denotes the high-order statistics including skewness (γi) and kurtosis (κi) to measure the degree of nonlinearity and quantify the non-Gaussianity of the process variables [31] as(8)γit=1/h∑q=0h−1xit−q−μit3vit3/2,(9)κit=1/h∑q=0h−1xit−q−μit4vit2−3.The SP is constructed by stackingλ selected statistics by equations (1)‒(9) into a row vector. In the paper, three kinds of statistics μit, vit, and κit are selected from one-order, two-order, and high-order statistics to construct statistics matrix S∈ℝN×M with N=n−h+1, M=λ×m, and λ=3.


### 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 datasetS=s1T,s2T,…,sNTT∈ℝN×M is obtained according to equation (2), and it is normalized as S¯. In the procedure of SVDD modeling, S¯ is firstly mapped from the original space to a higher feature space by a nonlinear transformation function Φs¯i. 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:(10)minR,a,ξR2+C∑i=1Nwiξis.t.Φs¯i−a≤R2+ξi,ξi≥0,where a and R are the center and the radius of hypersphere, respectively; trade-off parameter C=1/p×N is introduced to make the hypersphere as small as possible while preventing the misclassification of the samples, in which p denotes the percentage of outliers permitted in the training set [26]; ξi is slack variable that allows outliers of the hypersphere; and wi is weight of the ith data, where smaller wi 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) [27] 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 [32]. 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 ofs¯ii=1,2,…,N are selected to construct a set KNNs¯i, which satisfies K=K1+K2, where K1 and K2 are the numbers of neighbors in space and time, respectively.Suppose thatK−distances¯i is the Euclidean distance from s¯ to its Kth-nearest neighbor. Then the distance ds¯i,s¯ik between s¯i and s¯ik (k = 1, 2, …, K) is computed, where s¯ik is the kth-nearest neighbor of s¯i. The reachability distance of s¯i with respect to s¯ik is defined as(11)reach−ds¯i,s¯ik=maxK−distances¯ik,ds¯i,s¯ik.Thus, mLRD is constructed as(12)mLRDs¯i=K∑k=1Kreach−ds¯i,s¯ik.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(13)ws¯i=mLRDRs¯i=mLRDs¯imaxs¯j∈KNNs¯imLRDs¯i,mLRDs¯j.It is obvious thatws¯i is kept within 0,1.According to equation (13), if s¯i is a normal point, its mLRD is close to its neighbors and the corresponding mLRDR approaches 1. If s¯i 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 [26], 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 α=α1,α1,…,αNT as follows:(14)minαi∑i=1NαiKers¯i,s¯i−∑i=1N∑j=1NαiαjKers¯i,s¯js.t.∑i=1Nαi=1,0≤αi≤wiC,where Kers¯i,s¯j=Φs¯i⋅Φs¯j=e−s¯i−s¯jTs¯i−s¯j/σ2, is kernel function and σ can be determined via detecting the “tightness” of the decision boundaries [33]. After solving equation (14), only objects s¯i with αi>0 are called the support vectors (SVs), and their subscripts set is SV=i|αi>0,i=1,2,…,N. Hence, the center and radius of the hypersphere are(15)a=∑i∈SVαiΦs¯i,(16)R=Kers¯sv,s¯sv−2∑i=1NαiKers¯sv,s¯i+∑i=1N∑j=1NKers¯i,s¯j,respectively.Whenh new samples are obtained for a window, SPs are calculated by equation (2) and normalized as s¯new, and the distance between s¯new and the center R can be calculated by(17)Dists¯new=Kers¯new,s¯new−2∑i=1NαiKers¯new,s¯i+∑i=1N∑j=1NKers¯i,s¯j.IfDists¯new≤R, 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 Figure1. The specific offline modeling and online monitoring steps are listed as follows:Offline modeling procedure:(i)
Step 1: adopt SPA on a multimode datasetX∈ℝN×M to obtain statistics matrix S≡μ|V|Ξ∈ℝN×M by equations (2)–(9).(ii)
Step 2: normalize statistics matrixS and obtain S¯.(iii)
Step 3: compute the mLRDR as a weight factorWs¯i by equations (11)–(13).(iv)
Step 4: set up weighted-SVDD model with the weightW, solve equation (14), and obtain the center and the radius of the hypersphere by equations (15) and (16), respectively.Online monitoring procedure:(i)
Step 1: for a current window datasetxnew, statistics vector snew∈ℝM is calculated by equations (2)–(9).(ii)
Step 2: normalize statistics vectorsnew to obtain s¯new.(iii)
Step 3: calculate the distanceDist between s¯new and the center of the hypersphere by (17fd17).(iv)
Step 4: ifDists¯new≤R, xnew is considered to be normal; otherwise, it is in a faulty state.Figure 1
Flowchart of SmLRDR-wSVDD scheme for multimode monitoring.In this paper, ifDists of 3successive windows are over the control limit R, 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.


## 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 [28], 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 thatX∈ℝn×m is the original measured dataset with n samples and m variables, and a window of samples is denoted as(1)Xth=xt,1h,xt,2h,…,xt,mh=x1t−h+1x2t−h+1…xmt−h+1x1t−h+2x2t−h+2…xmt−h+2⋮⋮⋮⋮x1tx2t…xmt,where h is the window width and t 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 as(2)S≡μ|V|Ξ,where μ denotes the first-order statistics, that is, the mean, and its elements are calculated from the data in a window:(3)μit=1h∑q=0h−1xit−q,where V denotes the second-order statistics, which includes variance (vi), correlation (ri,j), autocorrelation (ri), and cross-correlation (ri,j) of different variables in the window as(4)vit=1h∑q=0h−1xit−q−μit2,(5)ri,jt=1h∑q=0h−1xit−q−μitxjt−q−μitvitvjt,(6)ridt=1h−d∑q=dh−1xit−q−μitxit+d−q−μitvit,(7)ri,jdt=1h−d∑q=dh−1xit−q−μitxjt+d−q−μitvitvjt,with d denoting the time lag between the variables; Ξ denotes the high-order statistics including skewness (γi) and kurtosis (κi) to measure the degree of nonlinearity and quantify the non-Gaussianity of the process variables [31] as(8)γit=1/h∑q=0h−1xit−q−μit3vit3/2,(9)κit=1/h∑q=0h−1xit−q−μit4vit2−3.The SP is constructed by stackingλ selected statistics by equations (1)‒(9) into a row vector. In the paper, three kinds of statistics μit, vit, and κit are selected from one-order, two-order, and high-order statistics to construct statistics matrix S∈ℝN×M with N=n−h+1, M=λ×m, and λ=3.


## 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 datasetS=s1T,s2T,…,sNTT∈ℝN×M is obtained according to equation (2), and it is normalized as S¯. In the procedure of SVDD modeling, S¯ is firstly mapped from the original space to a higher feature space by a nonlinear transformation function Φs¯i. 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:(10)minR,a,ξR2+C∑i=1Nwiξis.t.Φs¯i−a≤R2+ξi,ξi≥0,where a and R are the center and the radius of hypersphere, respectively; trade-off parameter C=1/p×N is introduced to make the hypersphere as small as possible while preventing the misclassification of the samples, in which p denotes the percentage of outliers permitted in the training set [26]; ξi is slack variable that allows outliers of the hypersphere; and wi is weight of the ith data, where smaller wi 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) [27] 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 [32]. 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 ofs¯ii=1,2,…,N are selected to construct a set KNNs¯i, which satisfies K=K1+K2, where K1 and K2 are the numbers of neighbors in space and time, respectively.Suppose thatK−distances¯i is the Euclidean distance from s¯ to its Kth-nearest neighbor. Then the distance ds¯i,s¯ik between s¯i and s¯ik (k = 1, 2, …, K) is computed, where s¯ik is the kth-nearest neighbor of s¯i. The reachability distance of s¯i with respect to s¯ik is defined as(11)reach−ds¯i,s¯ik=maxK−distances¯ik,ds¯i,s¯ik.Thus, mLRD is constructed as(12)mLRDs¯i=K∑k=1Kreach−ds¯i,s¯ik.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(13)ws¯i=mLRDRs¯i=mLRDs¯imaxs¯j∈KNNs¯imLRDs¯i,mLRDs¯j.It is obvious thatws¯i is kept within 0,1.According to equation (13), if s¯i is a normal point, its mLRD is close to its neighbors and the corresponding mLRDR approaches 1. If s¯i 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 [26], 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 α=α1,α1,…,αNT as follows:(14)minαi∑i=1NαiKers¯i,s¯i−∑i=1N∑j=1NαiαjKers¯i,s¯js.t.∑i=1Nαi=1,0≤αi≤wiC,where Kers¯i,s¯j=Φs¯i⋅Φs¯j=e−s¯i−s¯jTs¯i−s¯j/σ2, is kernel function and σ can be determined via detecting the “tightness” of the decision boundaries [33]. After solving equation (14), only objects s¯i with αi>0 are called the support vectors (SVs), and their subscripts set is SV=i|αi>0,i=1,2,…,N. Hence, the center and radius of the hypersphere are(15)a=∑i∈SVαiΦs¯i,(16)R=Kers¯sv,s¯sv−2∑i=1NαiKers¯sv,s¯i+∑i=1N∑j=1NKers¯i,s¯j,respectively.Whenh new samples are obtained for a window, SPs are calculated by equation (2) and normalized as s¯new, and the distance between s¯new and the center R can be calculated by(17)Dists¯new=Kers¯new,s¯new−2∑i=1NαiKers¯new,s¯i+∑i=1N∑j=1NKers¯i,s¯j.IfDists¯new≤R, 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 Figure1. The specific offline modeling and online monitoring steps are listed as follows:Offline modeling procedure:(i)
Step 1: adopt SPA on a multimode datasetX∈ℝN×M to obtain statistics matrix S≡μ|V|Ξ∈ℝN×M by equations (2)–(9).(ii)
Step 2: normalize statistics matrixS and obtain S¯.(iii)
Step 3: compute the mLRDR as a weight factorWs¯i by equations (11)–(13).(iv)
Step 4: set up weighted-SVDD model with the weightW, solve equation (14), and obtain the center and the radius of the hypersphere by equations (15) and (16), respectively.Online monitoring procedure:(i)
Step 1: for a current window datasetxnew, statistics vector snew∈ℝM is calculated by equations (2)–(9).(ii)
Step 2: normalize statistics vectorsnew to obtain s¯new.(iii)
Step 3: calculate the distanceDist between s¯new and the center of the hypersphere by (17fd17).(iv)
Step 4: ifDists¯new≤R, xnew is considered to be normal; otherwise, it is in a faulty state.Figure 1
Flowchart of SmLRDR-wSVDD scheme for multimode monitoring.In this paper, ifDists of 3successive windows are over the control limit R, 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 [34], conventional SVDD, and LDR-wSVDD [26].


### 3.1. Tennessee-Eastman Benchmark Process

Tennessee-Eastman process (TE process) was put forward in 1993 by Downs and Vogel [35], which has been widely adopted for scientific research [36]. There are five units, reactor, condenser, compressor, stripper, and separator, in this nonlinear process. According to different G/H mass ratios (G and H are two products), 6 steady modes can be obtained as listed in Table 1.Table 1
The parameters of 6 modes.
ModeG (%)H (%)Product yield15050G,H: 7038kgh−121090G:140 8kgh−1, H:1266 9kgh−139010G:1000 0kgh−1, H:111 1kgh−145050Maximum productivity51090Maximum productivity69010Maximum productivityThe control scheme and normal operations were provided by Ricker [37] 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.Table 2
The parameters of 6 modes.
No.DescriptionType1A/C feed ratio, B composition constant (stream 4)Step2Composition, A/C ratio constant (stream 4)Step4Reactor cooling water inlet temperatureStep5Condenser cooling water inlet temperatureStep6A feed loss (stream 1)Step7C header pressure loss-reduced availability (stream 4)Step8A, B, C feed composition (stream 4)Step10C feed temperature (stream 4)Random variation11Reactor cooling water inlet temperatureRandom variation12Condenser cooling water inlet temperatureRandom variation13Reaction kineticsSlow drift14Reactor cooling water valveStickingFor 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 the201th sample. The confidence level of KPCA is chosen as 95%. For the other three algorithms, the Gaussian kernel width parameter σ is determined by the method in [33], and p=5% is selected for the penalty parameter C. In the proposed method, the window width h is assigned as 10, and the numbers of neighbors in space and time K1 and K2 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) [38] 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.Table 3
MAR for all test conditions on TE process.
FaultMode 1Mode 3KPCASVDDLDR-SVDDSmLRDR-SVDDKPCASVDDLDR-SVDDSmLRDR-SVDDT2SPET2SPE10.95130.0050.00250.002500.99250.01380.00880.00750.00520.87250.00750.0050.00250.00130.90830.270.03880.00750.036340.98630.97130000.98380.87500050.98630.97750.9650.90630.93130.97750.99750.08380.03250.002560.98540.57730000.98770.580200070.98380.93380000.97630.882500080.630.01880.01750.01750.01120.5650.01630.01630.0150.01100.96630.960.83380.53330.44130.95880.78630.5650.28750.105110.940.53720.16250.01500.94750.72130.2350.0150120.95880.81130.63630.48880.26880.60750.02130.01750.00750.0025130.46630.02380.02250.02250.02130.46750.070.05870.02620.0112140.98750.95380.70750.106200.9850.87250.60380.00750Table 4
FAR for all test conditions on TE process.
FaultMode 1Mode 3KPCASVDDLDR-SVDDSmLRDR-SVDDKPCASVDDLDR-SVDDSmLRDR-SVDDT2SPET2SPE10.030.020.010.020.010.010.0250.010.0250.0120.030.010.0050.0150.0050.0150.0250.0050.020.01540.020.0150.0050.0150.010.010.0250.0050.020.0150.020.0150.010.010.010.020.050.050.0150.01560.010.0150.0050.010.010.0150.050.0050.020.0270.020.020.0050.010.010.020.0550.0050.020.0280.0250.0150.010.020.0050.0150.050.010.0250.015100.0250.020.0050.020.0050.0150.040.0050.0250.01110.020.010.0050.0150.0050.010.040.0050.030.015120.020.0150.0050.0150.0050.020.040.010.0250.01130.020.0150.0050.020.010.010.0650.0050.0250.01140.020.010.010.020.0050.010.070.010.030.01Further, the detailed monitoring results of four methods for different faults of two modes are shown in Figures2–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.Figure 2
Monitoring results for fault 11 of Mode 1 by (a) KPCA, (b) conventional SVDD, (c) LDR-wSVDD, and (d) SmLRDR-wSVDD.
(a)(b)(c)(d)Figure 3
Monitoring results for fault 14 of Mode 1 by (a) KPCA, (b) conventional SVDD, (c) LDR-wSVDD, and (d) SmLRDR-wSVDD.
(a)(b)(c)(d)Figure 4
Monitoring results for fault 11 of Mode 3 by (a) KPCA, (b) conventional SVDD, (c) LDR-wSVDD, and (d) SmLRDR-wSVDD.
(a)(b)(c)(d)Figure 5
Monitoring results for fault 13 of Mode 3 by (a) KPCA, (b) conventional SVDD, (c) LDR-wSVDD, and (d) SmLRDR-wSVDD.
(a)(b)(c)(d)FARs of all methods are summarized in Table4. For all faults of Mode 1, FARs of KPCA and LDR-wSVDD are within 0.01,0.03 and 0.01,0.02, respectively, while those of conventional SVDD and SmLRDR-wSVDD are within 0,0.01. 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.Figure6 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.Figure 6
(a) LRD; the weight of each training data based on (b) LDR and (c) mLRDR.
(a)(b)(c)


### 3.2. Wastewater Treatment Process (WWTP)

In this subsection, benchmark simulation model no. 1 (BSM1) of WWTP [39] 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/.Figure 7
Process layout of the WWTP.By changing the flow rates and oxygen transfer coefficients in biological reactor, three different modes can be obtained [40]. 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.Table 5
FAR for all test conditions on TE process.
ParameterMode 1Mode 2Influent flow rate (m3/day)184469223Recycle flow rate (m3/day)184469223Internal recycle flow rate (m3/day)5533827669KLa-Tank 3 (h−1)107.5KLa-Tank 4 (h−1)107.5KLa-Tank 5 (h−1)3.53.5Table 6
Description and notation of the monitoring variables.
DescriptionNotationReadily biodegradable substrateSS,3 (Tank 3)SS,4 (Tank 4)SS,5 (Tank 5)OxygenSO,3 (Tank 3)SO,4 (Tank 4)SO,5 (Tank 5)Nitrate and nitrite nitrogenSNO,3 (Tank 3)SNO,4 (Tank 4)SNO,5 (Tank 5)NH4++NH3 nitrogenSNH,3 (Tank 3)SNH,4 (Tank 4)SNH,5 (Tank 5)AlkalinitySALK,3 (Tank 3)SALK,4 (Tank 4)SALK,5 (Tank 5)Total suspended solidTSS4 (Tank 4)Oxygen transfer coefficientKLa5 (Tank 4)Nitrate and nitrite nitrogenSettler underflowSNO,uPlant exitSNO,eThe advanced biological nitrogen removal processInfluent ammonium concentrationSNH,inThere are two disturbance types in the benchmark simulation [41]: 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 [41].In this example, the samples are collected every 15 minutes. 1344 samples of dry weather of Mode 1 and Mode 2 shown in Table7 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.Table 7
Six simulation cases in the BSM1 benchmark.
Disturbance typeSamplesSimulation conditionsNoneDry weather1st−672ndNormal Mode 1673rd−1344thNormal Mode 2ExternalCase 11st−800th,1001st−1344thNormal Mode 1801st−1000thRainInternalCase 21st−672ndNormal Mode 1673rd−1344thSpecific growth rate for autotrophs:μA from 0.5 to 0.3 each day in a linear fashionCase 31st−672ndNormal Mode 1673rd−1344thMaximum heterotrophic growth rateμH has a step change decreasing from 4.0 d−1 to 2.0 d−1Case 41st−672ndNormal Mode 1673rd−1344thDO controller setpoint: changed from 2 to 1 mg/L (operation)Case 51st−672ndNormal Mode 2673rd−1344thμH has a step change decreasing from 4.0 d−1 to 3.0 d−1, and maximum hetero-trophic decay rate (bH ) has a step change increasing from 0.3 d−1 to 0.1 d−1The 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 Table8. 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.Table 8
MAR under five cases on WWTP.
Cases no.KPCASVDDLDR-wSVDDSmLRDR-wSVDDT2SPECase 110.030.040.0350Case 210.21130.17410.12650.0164Case 30.97320.41370.07290.0030Case 40.75890000Case 50.96430.24110.15180.03130Specifically, the monitoring results for Case 1 are shown in Figure8. 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 672nd 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 672nd 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.Figure 8
Monitoring results for Case 1 by (a) KPCA, (b) conventional SVDD, (c) LDR-wSVDD, and (d) SmLRDR-wSVDD.
(a)(b)(c)(d)Figure 9
Monitoring results for Case 2 by (a) KPCA, (b) conventional SVDD, (c) LDR-wSVDD, and (d) SmLRDR-wSVDD.
(a)(b)(c)(d)As for FARs shown in Table9, 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.Table 9
FAR under five cases on WWTP.
Cases no.KPCASVDDLDR-wSVDDSmLRDR-wSVDDT2SPECase 10.00610.02270.02080.04110.0236Case 20.0060.01490.00150.01040.006Case 30.00450.01190.00450.01340.003Case 40.0060.01490.0030.01490.0015Case 50.00450.01040.00150.00150.0119Similarly, LRD, LDR, and mLRDR of each training sample are visualized in Figure10 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.Figure 10
(a) LRD; the weight of each piece of training data based on (b) LDR and (c) mLRDR.
(a)(b)(c)


## 3.1. Tennessee-Eastman Benchmark Process

Tennessee-Eastman process (TE process) was put forward in 1993 by Downs and Vogel [35], which has been widely adopted for scientific research [36]. There are five units, reactor, condenser, compressor, stripper, and separator, in this nonlinear process. According to different G/H mass ratios (G and H are two products), 6 steady modes can be obtained as listed in Table 1.Table 1
The parameters of 6 modes.
ModeG (%)H (%)Product yield15050G,H: 7038kgh−121090G:140 8kgh−1, H:1266 9kgh−139010G:1000 0kgh−1, H:111 1kgh−145050Maximum productivity51090Maximum productivity69010Maximum productivityThe control scheme and normal operations were provided by Ricker [37] 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.Table 2
The parameters of 6 modes.
No.DescriptionType1A/C feed ratio, B composition constant (stream 4)Step2Composition, A/C ratio constant (stream 4)Step4Reactor cooling water inlet temperatureStep5Condenser cooling water inlet temperatureStep6A feed loss (stream 1)Step7C header pressure loss-reduced availability (stream 4)Step8A, B, C feed composition (stream 4)Step10C feed temperature (stream 4)Random variation11Reactor cooling water inlet temperatureRandom variation12Condenser cooling water inlet temperatureRandom variation13Reaction kineticsSlow drift14Reactor cooling water valveStickingFor 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 the201th sample. The confidence level of KPCA is chosen as 95%. For the other three algorithms, the Gaussian kernel width parameter σ is determined by the method in [33], and p=5% is selected for the penalty parameter C. In the proposed method, the window width h is assigned as 10, and the numbers of neighbors in space and time K1 and K2 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) [38] 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.Table 3
MAR for all test conditions on TE process.
FaultMode 1Mode 3KPCASVDDLDR-SVDDSmLRDR-SVDDKPCASVDDLDR-SVDDSmLRDR-SVDDT2SPET2SPE10.95130.0050.00250.002500.99250.01380.00880.00750.00520.87250.00750.0050.00250.00130.90830.270.03880.00750.036340.98630.97130000.98380.87500050.98630.97750.9650.90630.93130.97750.99750.08380.03250.002560.98540.57730000.98770.580200070.98380.93380000.97630.882500080.630.01880.01750.01750.01120.5650.01630.01630.0150.01100.96630.960.83380.53330.44130.95880.78630.5650.28750.105110.940.53720.16250.01500.94750.72130.2350.0150120.95880.81130.63630.48880.26880.60750.02130.01750.00750.0025130.46630.02380.02250.02250.02130.46750.070.05870.02620.0112140.98750.95380.70750.106200.9850.87250.60380.00750Table 4
FAR for all test conditions on TE process.
FaultMode 1Mode 3KPCASVDDLDR-SVDDSmLRDR-SVDDKPCASVDDLDR-SVDDSmLRDR-SVDDT2SPET2SPE10.030.020.010.020.010.010.0250.010.0250.0120.030.010.0050.0150.0050.0150.0250.0050.020.01540.020.0150.0050.0150.010.010.0250.0050.020.0150.020.0150.010.010.010.020.050.050.0150.01560.010.0150.0050.010.010.0150.050.0050.020.0270.020.020.0050.010.010.020.0550.0050.020.0280.0250.0150.010.020.0050.0150.050.010.0250.015100.0250.020.0050.020.0050.0150.040.0050.0250.01110.020.010.0050.0150.0050.010.040.0050.030.015120.020.0150.0050.0150.0050.020.040.010.0250.01130.020.0150.0050.020.010.010.0650.0050.0250.01140.020.010.010.020.0050.010.070.010.030.01Further, the detailed monitoring results of four methods for different faults of two modes are shown in Figures2–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.Figure 2
Monitoring results for fault 11 of Mode 1 by (a) KPCA, (b) conventional SVDD, (c) LDR-wSVDD, and (d) SmLRDR-wSVDD.
(a)(b)(c)(d)Figure 3
Monitoring results for fault 14 of Mode 1 by (a) KPCA, (b) conventional SVDD, (c) LDR-wSVDD, and (d) SmLRDR-wSVDD.
(a)(b)(c)(d)Figure 4
Monitoring results for fault 11 of Mode 3 by (a) KPCA, (b) conventional SVDD, (c) LDR-wSVDD, and (d) SmLRDR-wSVDD.
(a)(b)(c)(d)Figure 5
Monitoring results for fault 13 of Mode 3 by (a) KPCA, (b) conventional SVDD, (c) LDR-wSVDD, and (d) SmLRDR-wSVDD.
(a)(b)(c)(d)FARs of all methods are summarized in Table4. For all faults of Mode 1, FARs of KPCA and LDR-wSVDD are within 0.01,0.03 and 0.01,0.02, respectively, while those of conventional SVDD and SmLRDR-wSVDD are within 0,0.01. 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.Figure6 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.Figure 6
(a) LRD; the weight of each training data based on (b) LDR and (c) mLRDR.
(a)(b)(c)


## 3.2. Wastewater Treatment Process (WWTP)

In this subsection, benchmark simulation model no. 1 (BSM1) of WWTP [39] 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/.Figure 7
Process layout of the WWTP.By changing the flow rates and oxygen transfer coefficients in biological reactor, three different modes can be obtained [40]. 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.Table 5
FAR for all test conditions on TE process.
ParameterMode 1Mode 2Influent flow rate (m3/day)184469223Recycle flow rate (m3/day)184469223Internal recycle flow rate (m3/day)5533827669KLa-Tank 3 (h−1)107.5KLa-Tank 4 (h−1)107.5KLa-Tank 5 (h−1)3.53.5Table 6
Description and notation of the monitoring variables.
DescriptionNotationReadily biodegradable substrateSS,3 (Tank 3)SS,4 (Tank 4)SS,5 (Tank 5)OxygenSO,3 (Tank 3)SO,4 (Tank 4)SO,5 (Tank 5)Nitrate and nitrite nitrogenSNO,3 (Tank 3)SNO,4 (Tank 4)SNO,5 (Tank 5)NH4++NH3 nitrogenSNH,3 (Tank 3)SNH,4 (Tank 4)SNH,5 (Tank 5)AlkalinitySALK,3 (Tank 3)SALK,4 (Tank 4)SALK,5 (Tank 5)Total suspended solidTSS4 (Tank 4)Oxygen transfer coefficientKLa5 (Tank 4)Nitrate and nitrite nitrogenSettler underflowSNO,uPlant exitSNO,eThe advanced biological nitrogen removal processInfluent ammonium concentrationSNH,inThere are two disturbance types in the benchmark simulation [41]: 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 [41].In this example, the samples are collected every 15 minutes. 1344 samples of dry weather of Mode 1 and Mode 2 shown in Table7 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.Table 7
Six simulation cases in the BSM1 benchmark.
Disturbance typeSamplesSimulation conditionsNoneDry weather1st−672ndNormal Mode 1673rd−1344thNormal Mode 2ExternalCase 11st−800th,1001st−1344thNormal Mode 1801st−1000thRainInternalCase 21st−672ndNormal Mode 1673rd−1344thSpecific growth rate for autotrophs:μA from 0.5 to 0.3 each day in a linear fashionCase 31st−672ndNormal Mode 1673rd−1344thMaximum heterotrophic growth rateμH has a step change decreasing from 4.0 d−1 to 2.0 d−1Case 41st−672ndNormal Mode 1673rd−1344thDO controller setpoint: changed from 2 to 1 mg/L (operation)Case 51st−672ndNormal Mode 2673rd−1344thμH has a step change decreasing from 4.0 d−1 to 3.0 d−1, and maximum hetero-trophic decay rate (bH ) has a step change increasing from 0.3 d−1 to 0.1 d−1The 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 Table8. 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.Table 8
MAR under five cases on WWTP.
Cases no.KPCASVDDLDR-wSVDDSmLRDR-wSVDDT2SPECase 110.030.040.0350Case 210.21130.17410.12650.0164Case 30.97320.41370.07290.0030Case 40.75890000Case 50.96430.24110.15180.03130Specifically, the monitoring results for Case 1 are shown in Figure8. 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 672nd 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 672nd 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.Figure 8
Monitoring results for Case 1 by (a) KPCA, (b) conventional SVDD, (c) LDR-wSVDD, and (d) SmLRDR-wSVDD.
(a)(b)(c)(d)Figure 9
Monitoring results for Case 2 by (a) KPCA, (b) conventional SVDD, (c) LDR-wSVDD, and (d) SmLRDR-wSVDD.
(a)(b)(c)(d)As for FARs shown in Table9, 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.Table 9
FAR under five cases on WWTP.
Cases no.KPCASVDDLDR-wSVDDSmLRDR-wSVDDT2SPECase 10.00610.02270.02080.04110.0236Case 20.0060.01490.00150.01040.006Case 30.00450.01190.00450.01340.003Case 40.0060.01490.0030.01490.0015Case 50.00450.01040.00150.00150.0119Similarly, LRD, LDR, and mLRDR of each training sample are visualized in Figure10 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.Figure 10
(a) LRD; the weight of each piece of training data based on (b) LDR and (c) mLRDR.
(a)(b)(c)


## 4. Conclusion

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.



---

*Source: 5517062-2021-04-21.xml*