# Modifiedt-Distribution Stochastic Neighbor Embedding Using Augmented Kernel Mahalanobis-Distance for Dynamic Multimode Chemical Process Monitoring

Authors: Haoyu Gu; Li Wang
Journal: International Journal of Chemical Engineering
 (2022)

---

## Abstract

The traditional data-driven process monitoring methods may not be applicable for the system which has dynamic and multimode characteristics. In this paper, a novel scheme namedmodified t-distribution stochastic neighbor embedding using augmented Mahalanobis-distance for dynamic multimode chemical process monitoring (AKMD-t-SNE) is proposed to realize dynamic multimodal process monitoring. First, the augmented matrix strategy is utilized to ensure the sample contains the autocorrelation of the process. Moreover, AKMD-t-SNE method eliminates the scale and spatial distribution differences among multiple modes by calculating the kernel Mahalanobis distance between the samples to establish the global model. The features extracted via the proposed method are obviously non-Gaussian, and there will be a deviation in the construction of traditional statistics. Then, the support vector data description (SVDD) method is used to construct statistics to deal with this problem. In addition, a hybrid correlation coefficient method (HCC) is proposed to achieve fault isolation and improve the accuracy of isolation results. The advantages of the proposed scheme are illustrated by a numerical case and the Multimode Tennessee Eastman Process (MTEP) benchmark.


---

## Body

## 1. Introduction

Process monitoring plays a crucial role in the safe and efficient production of process industries such as the chemical industry and oil refining industry. With the increasing market demand, the complexity and the difficulty of process monitoring also raise. Massive industrial data are stored with the assistance of distributed control system (DCS). The data-driven process monitoring method is to discover potential regularities from colossal historical data for modeling and monitoring of intricate industrial processes. The most fundamental method of data-driven process monitoring is multivariate statistical process monitoring (MSPM) [1–5].In process industries, the seasonal variation and demand fluctuation will lead to the change of the equilibrium point during production, which causes the process to have the characteristic of a multimode. In the multimode process, once the operation mode switches, the model based on a single-mode process will not be able to achieve effective monitoring. One typical multimode process is the penicillin fermentation process, which is divided into two modes. The conventional MSPM method is only suitable for the process monitoring in the production process of a single mode. Mode one was the growth and propagation mode of Penicillium, accompanied by the rapid consumption of grape sugar base, which lasted about 43 hours. Mode two is the synthetic penicillin mode, which requires continuous glucose supplementation to the fermenter in order to maintain the high yield and high quality of penicillin [6]. The methods of multimode process monitoring mainly include the model integration method based on a rigid partition strategy and the global modeling method based on a flexible partition strategy [7, 8]. Due to the inertia of equipment operation and the hysteresis response of closed-loop systems, the output of the system has the characteristic of time delay, which is ubiquitous in process industry production [9]. Because of the autocorrelation between the successive measurements, the modeling method for static objects is no longer applicable to the dynamic system. Modeling methods for dynamic systems commonly include augmented matrix, time series analysis, subspace identification, and deep learning [10–14]. The augmented matrix method is to form the extended matrix by combining the current sample and the earlier measurements, which has been proven to be able to extract both the autocorrelation of samples and the mutual correlation of variables. Ku first proposed dynamic principal component analysis (DPCA) based on an augmented matrix to solve the modeling problem of linear dynamic systems [15]. The augmented matrix method is still an effective strategy for dynamic process modeling in the current research studies. Dynamic multimodal process monitoring has been studied in recent years.t-distribution stochastic neighbor embedding (t-SNE) is an improved manifold learning [16, 17] algorithm based on stochastic neighbor embedding (SNE). Both global and local data structures are preserved when t-SNE is used to accomplish feature extraction. t-SNE is early used in the field of data visualization [18, 19]. Due to the high computational complexity of the t-SNE algorithm, various methods have been proposed to accelerate the t-SNE algorithm [20, 21]. The t-SNE algorithm has been used successfully in various domains due to its superior performance. However, the application of t-SNE in process monitoring is still infrequent. Zhang et al. proposed a quality-related process monitoring method for the hot rolling process based on parameter t-SNE (P-t-SNE) and modified minimum error minimax probability machine (MMEMPM) method which achieved satisfying results [22]. Tang and Yan used the t-SNE algorithm to realize process monitoring through data visualization, which can more intuitively show the trend of process data when an anomaly occurs [23]. In recent years, dynamic multimodal processes have been studied. Qin et al. proposed a mixture slow feature analysis (MSFA) method to deal with the health degradation problem. MSFA integrates the local model to establish the global model, to achieve the equipment fault detection in the multimode process [24]. Zhang et al. introduced a novel sparse dynamic inner principal component analysis (SDiPCA) based monitoring for multimode dynamic processes, which can update the model of sequential modes by memorizing the significant features of existing modes and provide outstanding performance for successive mode [25]. Deep learning is a hot topic in process monitoring and it is feasible to use deep learning to realize multimodal process monitoring. A novel self-adaptive deep learning method based on local adaptive standardization and variational auto-encoder bidirectional long short-term memory (LAS-VB) was proposed by Wu and Zhao and showed excellent performance for multimode process monitoring [26].The t-SNE algorithm first calculates the Euclidean distance between high-dimensional data and then converts the distance into conditional probability. The dimensional reduction is realized by using the low-dimensional conditional probability to approximate the high-dimensional conditional probability. The original t-SNE method of feature extraction is only suitable for the modeling of static single-mode process. In this article, an improved t-SNE method calledt-distribution stochastic neighbor embedding algorithm based on augmented kernel Mahalanobis distance (AKMD-t-SNE) is proposed for dynamic multimode chemical process monitoring. The augmented matrix strategy is employed to expand the samples, which can ensure the samples contain information from the adjacent sampling instants. Moreover, the kernel Mahalanobis distance between the augmented samples is calculated to eliminate the scale differences among different modes [27]. Furthermore, samples are nonidentically distributed in different modes, the kernel Mahalanobis distance transforms the samples in different modes into independent identically distributed problems via the kernel covariance matrix. Then, the augmented kernel Mahalanobis distance is introduced into the t-SNE algorithm for dimensionality reduction. The model decomposes the high-dimensional data into low-dimensional feature subspace (FS) and residual subspace (RS). The FS reflects the main information of the process and the RS contains the redundant information of the process. Since the data is non-Gaussian, the SVDD method is adopted in this paper to construct statistics in the FS and RS, which further improve the accuracy of fault detection [28–30].Fault detection is the fundamental step of process monitoring. After the abnormal condition is detected effectively, it is necessary to trace the variable that causes the abnormal status, i.e., to conduct a fault isolation step. The most typical fault isolation methods based on the projection model are contribution plots and reconstruction-based contribution (RBC) [31]. The contribution plots indicate the contribution value of each variable to the detection statistics, but the disturbance of one variable may affect the contribution of other uncorrelated variables to the statistics; i.e., the smearing effect, even if the fault amplitude is large, there is still the probability of misdiagnosis. The fault isolation method based on reconstruction was proposed by Dunia et al. in 1996. This approach needs to give the fault direction to realize reconstruction. By comparing the reconstructed statistics, we can judge whether the given fault direction is correct and determine the isolated variable consequently [32]. The RBC method highly relies on prior knowledge, and the computational complexity of reconstruction is high [33, 34]. Verron et al. proposed a fault isolation method based on mutual information (MI), which is computationally efficient and does not need to rely on historical data, but the smearing effect still exists [35]. To overcome the above-mentioned limitations of fault isolation, this work proposes a hybrid correlation coefficient (HCC) method based on the MI to achieve fault isolation and exclude suspicious variables, which ensures the isolated variables are more explicit and consistent with the responsible variables. In short, this paper provides a global modeling method that integrates multiple modes and dynamic characteristics.


## 2. Materials and Methods

### 2.1. Process Feature Extraction Method

Due to the autocorrelation of process data, an augmented matrix is established to ensure the new samples contain the cross-correlation between variables and the time correlation between observations. The high-dimensional data is=x1,x2,…,xN∈RN×D, where D represents the number of variables and N represents the number of observations. Suppose the current sampling instant is k and the sample at k is xk. xk and past values with orders of data dynamics L forms an augmented matrix:(1)xkL=xk,xk−1,⋯,xk−L∈R1×L+1D.In equation (1), the sampling time k must be larger than the time delay L, and the form of the augmented matrix Xa formed is as follows:(2)Xb=x1+Lxl⋯x1x2+Lx1+L⋯x2⋮⋮⋱⋮xNxN−1⋯xN−L∈RN−L×L+1D,k>L,when the sampling time k is smaller than the time delay L, the form of the augmented matrix Xb is as follows:(3)Xa=x1x1⋯x1x2x1⋯x1⋮⋮⋱⋮xLxL−1⋯x1∈RL×L+1D,k≤L.Then, the final augmented matrix formed by the training dataX as XL:(4)XL=Xa;Xb∈RN×L+1D.The AKMD-t-SNE algorithm is used to map the augmented high-dimensional dataXL∈RN×l+1D to the low-dimensional feature subspace YL=y1,y2,...,yN∈Rd×N, where XL is a matrix with N observations and l+1D dimensions. YL is the matrix with N observations and d dimensions.The Mahalanobis distance could eliminate the discrepancies among different modes of training data.xi,xj are two samples in the high-dimensional data. The Mahalanobis distance between high-dimensional samples xi and xj is calculated as(5)DMxi,xj=xi−xjTΣ−1xi−xj,where Σ is the covariance matrix of XL.Since Mahalanobis distance is not suitable to deal with nonlinear problems, the kernel trick is introduced for the purpose of capturing the nonlinearity. By virtue of a nonlinear transformationϕk=Fxk, the augmented vector is mapped into a higher-dimensional feature space. Denote the mapped matrix as Φ=ϕ1,ϕ2,⋯,ϕN. The augmented kernel Mahalanobis distance is determined as(6)DKMϕi,ϕj=ϕi−ϕjTΣϕ−1ϕi−ϕj,where ϕi,ϕj∈Φ, Σϕ represents the covariance matrix of Φ:(7)Σϕ=1N−1∑k=1Nϕk−μϕϕk−μϕΤ=1N−1Φ¯Φ¯Τ,(8)μϕ=1N∑k=1Nϕk,where Φ¯=ϕ¯1,ϕ¯2,⋯,ϕ¯N satisfies Φ¯=ΦH, H is a symmetric and idempotent matrix:(9)H=IN−1NININΤ.The dot product of mapped vectors is calculated by the Gaussian kernel function:(10)Kxi,xj=ϕiΤϕj=exp−12σ2xi−xjΤxi−xj.Then, the centered kernel matrix is defined as(11)K¯=Φ¯ΤΦ¯=HΤΦΤΦH=HKH.A full rank factorization via singular value decomposition (SVD) method is analyzed:(12)Φ¯=USVΤ,where U and V are both orthogonal matrices.From (7), we have(13)Σϕ=1N−1USSΤUΤ.Then, its inverse matrix is calculated:(14)Σϕ−1=N−1USSΤ−1UΤ.The centered kernel matrix and its Moore-Penrose pseudoinverse matrix are given as(15)K¯=SVΤVSΤ,(16)K¯+=VSΤS+VΤ.Combining (15) with (17), we obtain(17)Σϕ−1=N−1Φ¯K¯+2Φ¯Τ.The augmented kernel Mahalanobis distance can be rewritten as(18)DKMϕi,ϕj=ϕi−ϕjΣϕ−1ϕi−ϕj=N−1ϕ¯i−ϕ¯jΦ¯K¯+2Φ¯Τϕ¯i−ϕ¯jΤ=N−1k⌢i,jΤHK¯+2Hk⌢i,j,where k⌢i,j=Φϕi−Φϕj.The conditional probability is calculated between samples in higher-dimensional space. The conditional probability of sampleϕj as the nearest neighbor of ϕi in higher-dimensional space is denoted as pji:(19)pji=exp−DKM2ϕi,ϕj/2σi2∑k≠iexp−DKM2ϕi,ϕk/2σi2,i≠j,0,i=j,where σi represents the variance of the Gaussian distribution with ϕi as the center, and ϕk represents one of the nearest neighbors of ϕi.The selection ofσi value is related to the degree of Perplexity, which is defined as(20)PrepPi=2HPi,where Hpi is the Shannon entropy of pi:(21)HPi=−∑jpj|ilog2pj|i.Assuming that the mapping samples of higher-dimensional observationsϕi and ϕj in the low-dimensional space are yi and yj, and the variance of the Gaussian distribution with yi as the center is set to 1/2, the corresponding conditional probability that yj is the nearest sample of yi denoted as qji:(22)qj|i=exp−yi−yj2∑k≠iexp−yi−yk2,i≠j,0,i=j,where yk represents one of the nearest neighbors of yi.Kullback–Leibler Divergence is used to measure the similarity between the two individual distributions. The similarity cost functionC as:(23)C=∑iKLPiQi=∑i∑jpj|ilogpj|iqj|i.The cost functionC is minimized by gradient descent, and the gradient of the cost function C is obtained:(24)δCδyi=2∑jpj|i−qj|i+pi|j−qi|jyi−yj.The conditional probability is replaced by the joint probability distribution, such that for anyi and j, qji=qij, pji=pij. The following symmetric method is adopted to simplify the objective function:(25)pij=pj|i+pi|j2N.The cost functionC is transformed into the following expression:(26)C=KLPQ=∑i∑jpijlogpijqij.In order to solve the crowding problem of low-dimensional data in classification,t-distribution is used to replace Gaussian distribution. Redefine qij with a t-distribution of 1 degree of freedom, it follows that:(27)qij=1+yi−yj2−1∑k≠l1+yk−yl2−1.The gradient formula is converted to the following equation:(28)δCδyi=4∑jpij−qijyi−yj1+yi−yj2−1.To avoid falling into the local optimum in the gradient descent process, a momentum term is added in the gradient update, as follows:(29)yt=yt−1+ηδCδy+αtyt−1−yt−2,where yt represents the solution of the t th iteration, αt represents the momentum of the t th iteration, and η represents the learning efficiency.The low-dimensional feature subspaceYL is obtained through the above steps. It is assumed that the higher-dimensional data are projected to the lower-dimensional features under the action of the mapping matrix A, which can be obtained from the Moore–Penrose pseudoinverse matrix XL+ of XL:(30)A=YLXL+.The feature subspaceY is recalculated via the mapping matrix and the residual subspace is denoted as E:(31)Y=AXL,E=XL−XL^=XL−A+Y,where XL^ is the high-dimensional projection of low-dimensional data and A+ is the Moore–Penrose pseudoinverse of A.For testing dataxon, the augmented matrix xonL is also constructed to ensure that the dimensions of online data are consistent with training data. The low-dimensional features of the augmented online data are calculated via the mapping matrix A:(32)yonL=AxonL.The residual ofxonL is given as(33)eonL=xonL−xonL^=xonL−A+yonL.


### 2.2. Fault Detection Method

After dimension-reduction by the AKMD-t-SNE scheme, the low-dimensional data will be expressed ast-distribution, seriously deviating from the normal distribution. Whereas the effective establishment of traditional statistics needs to be carried out under the condition that the data obey the Gaussian distribution. Consequently, the establishment of statistics using T2 and SPE will affect the estimation of the control limit, leading to the unsatisfactory effect of fault detection. Hence, this section uses the SVDD method to construct statistics. The SVDD algorithm is a single classification model proposed by Tax. The target of SVDD is to seek a hypersphere with center a and radius R, which can contain all the training data. SVDD statistics were established in the FS (AKMD-t-SNE-FS-SVDD) and RS (AKMD-t-SNE-RS-SVDD), respectively.To establish the hypersphere model for the FS of training data is to solve the following optimization problem:(34)minR,a,ξRy2+C1∑i=1nξi,s.t.ΦYi−ay2≤Ry2+ξi,ξi≥0,where Yi is the i th support vector of the FS; ξi is the slack variable to improve the fault tolerance of the SVDD model; C1 is the penalty factor, the size of the C1 reflects the extent of choosing outliers for a data set.To achieve fault detection on the feature of online data (yonL) is to calculate the distance between the sample and spherical center ay and compare it with Ry.(35)Distancey=ΦyonL−ay=1−2∑i=1nαiKyonL,Yo+∑i=1n∑j=1nαiαjKYi,Yj>Rys.t.0⩽αi⩽C1,∑i=1Nαi=1,where Yo is one of the support vectors in the FS of training data, Yi and Yj is the i th and j th support vector in the FS of training data.For the residual of online data (eonL), SVDD statistics (Distancee) and control limit (Re) are calculated in the same way to determine whether a fault occurs.


### 2.3. Fault Isolation Method

After the abnormal condition of the system is detected, it is necessary to diagnose the cause of the fault. Generally, the location of the disturbance is analyzed by isolating the fault variable. It has been proved that the mutual information (MI) method can effectively realize fault isolation, but it still has some defects. For the purpose of eliminating the influence of suspicious variables on the isolation results, an isolation method based on the hybrid correlation coefficient (HCC) is proposed to modify the MI method.MI is a measurement of the correlation between variables. Suppose there are two random variablesVx and Vy, and the MI of Vx and Vy is defined as IVx,Vy:(36)IVX,VY=∑vx∈VX∑vy∈VYpvx,vylogpvx,vypvxpvy,where pvx and pvy are the marginal distribution of Vx and Vy respectively, and pvx,vy are the joint distribution of Vx and Vy.Once the abnormal is detected in the augmented online data, the augmented online data containing normal data and abnormal data is denoted asZonL=Zll=1L+1∈RL+1D×N1, where N1 represents the number of samples and L+1D represents the number of variables. Since the augmented matrix contains the data of the previous L sampling instants, ZonL is composed of submatrices Zl, and the expression is given as(37)ZonL=z11,⋯,z1D⏟Z1∈RD×N1,z21,⋯,z2D⏟Z2∈RD×N1,⋯,zL+11,⋯,zL+1D⏟ZL+1∈RD×N1.The FS statistics ofZonL constitute the vector Zy. Compute the MI between Zy and each vector in each submatrix as(38)MIylq=IZlq,Zyq=1,2,⋯,D.Compute the MI between each vector in the online data andZy as(39)MIyq=∑l=1L+1MIylqL+1q=1,2,⋯,D.On this basis, a fault isolation method based on correlation coefficient (CC) is proposed. In order to understand the correlation level of variables expressed byMIy, the MI between Zy and itself is calculated as IZy,Zy. The correlation coefficient of FS (CCy) is established as(40)CCyq=MIy3qI2Zy,Zyq=1,2,⋯,D.The impact of each variable in the FS on the occurrence of disturbances was judged by the size ofCCyq.Similarly, the variable composed of RS statistics is denoted asZe, and the CC of RS (CCe) is calculated in the same way to achieve fault isolation.The hybrid mutual information (HMI) and the hybrid correlation coefficient (HCC) were defined as(41)HMIq=ωMIyq∑q=1DMIyq+1−ωMIeq∑q=1DMIeq×100%,HCCq=ωCCyq∑q=1DCCyq+1−ωCCeq∑q=1DCCeq×100%.Calculate the weight coefficientω:(42)ω1=ZyRy,ω2=ZeRe,ω=ω1ω1+ω2.The process monitoring flow diagram of the proposed method is shown in Figure1.Figure 1
Dynamic multimode process monitoring flow chart based on AKMD-t-SNE-SVDD and HCC.


## 2.1. Process Feature Extraction Method

Due to the autocorrelation of process data, an augmented matrix is established to ensure the new samples contain the cross-correlation between variables and the time correlation between observations. The high-dimensional data is=x1,x2,…,xN∈RN×D, where D represents the number of variables and N represents the number of observations. Suppose the current sampling instant is k and the sample at k is xk. xk and past values with orders of data dynamics L forms an augmented matrix:(1)xkL=xk,xk−1,⋯,xk−L∈R1×L+1D.In equation (1), the sampling time k must be larger than the time delay L, and the form of the augmented matrix Xa formed is as follows:(2)Xb=x1+Lxl⋯x1x2+Lx1+L⋯x2⋮⋮⋱⋮xNxN−1⋯xN−L∈RN−L×L+1D,k>L,when the sampling time k is smaller than the time delay L, the form of the augmented matrix Xb is as follows:(3)Xa=x1x1⋯x1x2x1⋯x1⋮⋮⋱⋮xLxL−1⋯x1∈RL×L+1D,k≤L.Then, the final augmented matrix formed by the training dataX as XL:(4)XL=Xa;Xb∈RN×L+1D.The AKMD-t-SNE algorithm is used to map the augmented high-dimensional dataXL∈RN×l+1D to the low-dimensional feature subspace YL=y1,y2,...,yN∈Rd×N, where XL is a matrix with N observations and l+1D dimensions. YL is the matrix with N observations and d dimensions.The Mahalanobis distance could eliminate the discrepancies among different modes of training data.xi,xj are two samples in the high-dimensional data. The Mahalanobis distance between high-dimensional samples xi and xj is calculated as(5)DMxi,xj=xi−xjTΣ−1xi−xj,where Σ is the covariance matrix of XL.Since Mahalanobis distance is not suitable to deal with nonlinear problems, the kernel trick is introduced for the purpose of capturing the nonlinearity. By virtue of a nonlinear transformationϕk=Fxk, the augmented vector is mapped into a higher-dimensional feature space. Denote the mapped matrix as Φ=ϕ1,ϕ2,⋯,ϕN. The augmented kernel Mahalanobis distance is determined as(6)DKMϕi,ϕj=ϕi−ϕjTΣϕ−1ϕi−ϕj,where ϕi,ϕj∈Φ, Σϕ represents the covariance matrix of Φ:(7)Σϕ=1N−1∑k=1Nϕk−μϕϕk−μϕΤ=1N−1Φ¯Φ¯Τ,(8)μϕ=1N∑k=1Nϕk,where Φ¯=ϕ¯1,ϕ¯2,⋯,ϕ¯N satisfies Φ¯=ΦH, H is a symmetric and idempotent matrix:(9)H=IN−1NININΤ.The dot product of mapped vectors is calculated by the Gaussian kernel function:(10)Kxi,xj=ϕiΤϕj=exp−12σ2xi−xjΤxi−xj.Then, the centered kernel matrix is defined as(11)K¯=Φ¯ΤΦ¯=HΤΦΤΦH=HKH.A full rank factorization via singular value decomposition (SVD) method is analyzed:(12)Φ¯=USVΤ,where U and V are both orthogonal matrices.From (7), we have(13)Σϕ=1N−1USSΤUΤ.Then, its inverse matrix is calculated:(14)Σϕ−1=N−1USSΤ−1UΤ.The centered kernel matrix and its Moore-Penrose pseudoinverse matrix are given as(15)K¯=SVΤVSΤ,(16)K¯+=VSΤS+VΤ.Combining (15) with (17), we obtain(17)Σϕ−1=N−1Φ¯K¯+2Φ¯Τ.The augmented kernel Mahalanobis distance can be rewritten as(18)DKMϕi,ϕj=ϕi−ϕjΣϕ−1ϕi−ϕj=N−1ϕ¯i−ϕ¯jΦ¯K¯+2Φ¯Τϕ¯i−ϕ¯jΤ=N−1k⌢i,jΤHK¯+2Hk⌢i,j,where k⌢i,j=Φϕi−Φϕj.The conditional probability is calculated between samples in higher-dimensional space. The conditional probability of sampleϕj as the nearest neighbor of ϕi in higher-dimensional space is denoted as pji:(19)pji=exp−DKM2ϕi,ϕj/2σi2∑k≠iexp−DKM2ϕi,ϕk/2σi2,i≠j,0,i=j,where σi represents the variance of the Gaussian distribution with ϕi as the center, and ϕk represents one of the nearest neighbors of ϕi.The selection ofσi value is related to the degree of Perplexity, which is defined as(20)PrepPi=2HPi,where Hpi is the Shannon entropy of pi:(21)HPi=−∑jpj|ilog2pj|i.Assuming that the mapping samples of higher-dimensional observationsϕi and ϕj in the low-dimensional space are yi and yj, and the variance of the Gaussian distribution with yi as the center is set to 1/2, the corresponding conditional probability that yj is the nearest sample of yi denoted as qji:(22)qj|i=exp−yi−yj2∑k≠iexp−yi−yk2,i≠j,0,i=j,where yk represents one of the nearest neighbors of yi.Kullback–Leibler Divergence is used to measure the similarity between the two individual distributions. The similarity cost functionC as:(23)C=∑iKLPiQi=∑i∑jpj|ilogpj|iqj|i.The cost functionC is minimized by gradient descent, and the gradient of the cost function C is obtained:(24)δCδyi=2∑jpj|i−qj|i+pi|j−qi|jyi−yj.The conditional probability is replaced by the joint probability distribution, such that for anyi and j, qji=qij, pji=pij. The following symmetric method is adopted to simplify the objective function:(25)pij=pj|i+pi|j2N.The cost functionC is transformed into the following expression:(26)C=KLPQ=∑i∑jpijlogpijqij.In order to solve the crowding problem of low-dimensional data in classification,t-distribution is used to replace Gaussian distribution. Redefine qij with a t-distribution of 1 degree of freedom, it follows that:(27)qij=1+yi−yj2−1∑k≠l1+yk−yl2−1.The gradient formula is converted to the following equation:(28)δCδyi=4∑jpij−qijyi−yj1+yi−yj2−1.To avoid falling into the local optimum in the gradient descent process, a momentum term is added in the gradient update, as follows:(29)yt=yt−1+ηδCδy+αtyt−1−yt−2,where yt represents the solution of the t th iteration, αt represents the momentum of the t th iteration, and η represents the learning efficiency.The low-dimensional feature subspaceYL is obtained through the above steps. It is assumed that the higher-dimensional data are projected to the lower-dimensional features under the action of the mapping matrix A, which can be obtained from the Moore–Penrose pseudoinverse matrix XL+ of XL:(30)A=YLXL+.The feature subspaceY is recalculated via the mapping matrix and the residual subspace is denoted as E:(31)Y=AXL,E=XL−XL^=XL−A+Y,where XL^ is the high-dimensional projection of low-dimensional data and A+ is the Moore–Penrose pseudoinverse of A.For testing dataxon, the augmented matrix xonL is also constructed to ensure that the dimensions of online data are consistent with training data. The low-dimensional features of the augmented online data are calculated via the mapping matrix A:(32)yonL=AxonL.The residual ofxonL is given as(33)eonL=xonL−xonL^=xonL−A+yonL.


## 2.2. Fault Detection Method

After dimension-reduction by the AKMD-t-SNE scheme, the low-dimensional data will be expressed ast-distribution, seriously deviating from the normal distribution. Whereas the effective establishment of traditional statistics needs to be carried out under the condition that the data obey the Gaussian distribution. Consequently, the establishment of statistics using T2 and SPE will affect the estimation of the control limit, leading to the unsatisfactory effect of fault detection. Hence, this section uses the SVDD method to construct statistics. The SVDD algorithm is a single classification model proposed by Tax. The target of SVDD is to seek a hypersphere with center a and radius R, which can contain all the training data. SVDD statistics were established in the FS (AKMD-t-SNE-FS-SVDD) and RS (AKMD-t-SNE-RS-SVDD), respectively.To establish the hypersphere model for the FS of training data is to solve the following optimization problem:(34)minR,a,ξRy2+C1∑i=1nξi,s.t.ΦYi−ay2≤Ry2+ξi,ξi≥0,where Yi is the i th support vector of the FS; ξi is the slack variable to improve the fault tolerance of the SVDD model; C1 is the penalty factor, the size of the C1 reflects the extent of choosing outliers for a data set.To achieve fault detection on the feature of online data (yonL) is to calculate the distance between the sample and spherical center ay and compare it with Ry.(35)Distancey=ΦyonL−ay=1−2∑i=1nαiKyonL,Yo+∑i=1n∑j=1nαiαjKYi,Yj>Rys.t.0⩽αi⩽C1,∑i=1Nαi=1,where Yo is one of the support vectors in the FS of training data, Yi and Yj is the i th and j th support vector in the FS of training data.For the residual of online data (eonL), SVDD statistics (Distancee) and control limit (Re) are calculated in the same way to determine whether a fault occurs.


## 2.3. Fault Isolation Method

After the abnormal condition of the system is detected, it is necessary to diagnose the cause of the fault. Generally, the location of the disturbance is analyzed by isolating the fault variable. It has been proved that the mutual information (MI) method can effectively realize fault isolation, but it still has some defects. For the purpose of eliminating the influence of suspicious variables on the isolation results, an isolation method based on the hybrid correlation coefficient (HCC) is proposed to modify the MI method.MI is a measurement of the correlation between variables. Suppose there are two random variablesVx and Vy, and the MI of Vx and Vy is defined as IVx,Vy:(36)IVX,VY=∑vx∈VX∑vy∈VYpvx,vylogpvx,vypvxpvy,where pvx and pvy are the marginal distribution of Vx and Vy respectively, and pvx,vy are the joint distribution of Vx and Vy.Once the abnormal is detected in the augmented online data, the augmented online data containing normal data and abnormal data is denoted asZonL=Zll=1L+1∈RL+1D×N1, where N1 represents the number of samples and L+1D represents the number of variables. Since the augmented matrix contains the data of the previous L sampling instants, ZonL is composed of submatrices Zl, and the expression is given as(37)ZonL=z11,⋯,z1D⏟Z1∈RD×N1,z21,⋯,z2D⏟Z2∈RD×N1,⋯,zL+11,⋯,zL+1D⏟ZL+1∈RD×N1.The FS statistics ofZonL constitute the vector Zy. Compute the MI between Zy and each vector in each submatrix as(38)MIylq=IZlq,Zyq=1,2,⋯,D.Compute the MI between each vector in the online data andZy as(39)MIyq=∑l=1L+1MIylqL+1q=1,2,⋯,D.On this basis, a fault isolation method based on correlation coefficient (CC) is proposed. In order to understand the correlation level of variables expressed byMIy, the MI between Zy and itself is calculated as IZy,Zy. The correlation coefficient of FS (CCy) is established as(40)CCyq=MIy3qI2Zy,Zyq=1,2,⋯,D.The impact of each variable in the FS on the occurrence of disturbances was judged by the size ofCCyq.Similarly, the variable composed of RS statistics is denoted asZe, and the CC of RS (CCe) is calculated in the same way to achieve fault isolation.The hybrid mutual information (HMI) and the hybrid correlation coefficient (HCC) were defined as(41)HMIq=ωMIyq∑q=1DMIyq+1−ωMIeq∑q=1DMIeq×100%,HCCq=ωCCyq∑q=1DCCyq+1−ωCCeq∑q=1DCCeq×100%.Calculate the weight coefficientω:(42)ω1=ZyRy,ω2=ZeRe,ω=ω1ω1+ω2.The process monitoring flow diagram of the proposed method is shown in Figure1.Figure 1
Dynamic multimode process monitoring flow chart based on AKMD-t-SNE-SVDD and HCC.


## 3. Results and Discussion

### 3.1. Numerical Case

To testify the effectiveness of the proposed method, a numerical dynamic multimode process is simulated. The dynamic multimode model is created as follows:(43)x1=0.6352s1+0.1254s2+e1,x2=0.6951s1+0.1557ln1+s22+e2,x3=0.5028s1s2+0.8312ln2+coss2+e3,x4=0.4152s1+0.2617s22+e4,x5=0.0569x3+0.7962sinx4+e5,where s1 and s2 are latent variables, e1,e2,e3,e4,e5T are independent zero-mean Gaussian white noises with and a standard deviation of 0.01. Variable x5 is the correlation variable of variable x3 and variable x4. Set three modes for the latent variable as(44)Mode1:s1τ:Uniform−5,−1s2τ:N−2,2,Mode2:s1τ=0.8cos2s1τ−1+5esinτ−1s2τ=1.5coses2τ−1+ln1+s22τ−1,Mode3:s1τ=0.1sin2s1τ−1+2.8ecos2s1τ−1s2τ=3.25e0.3sins2τ−1+0.1coss2τ−1,where τ represents the index of the current sampling instant. The initial state of mode 2 and mode 3 is 0, i.e., s11=s21=0.Five hundred samples were generated in each mode, that is, a total of 1500 samples were generated in the three modes to form the training data. In order to show the spatial distribution characteristics of the data of the three modes intuitively, Figure 2 shows the scatter plot of variables x2, x3, and x4 in three modes.Figure 2
Numerical case 3D scatter plot.Set two types of fault:Type 1: the system operates normally for 200 sampling instants in mode 2. At the 201st sampling instant (τ=201), a step signal with an amplitude of +1.35 is added to variable x2, and 800 sampling instants are continued in this state.Type 2: the system operates normally for 200 sampling instants in mode. At the 201st sampling instant (τ=201), the ramp signal of 0.03τ−200 is added to variable 3 and continues to operate for 800 sampling instants in this state.Each of the two fault types has 1000 test samples. The augment matrix-based dynamic locally linear embedding (DLLE), t-SNE, t-SNE based on Mahaobanobis distance (M-t-SNE), and AKMD-t-SNE are used to extract the features of the established dynamic multimode numerical model. SVDD statistics are established in all the four methods for fault detection. According to experience, the dimension reduction of the four methods is 2, the nearest neighbor of DLLE is 25, the penalty coefficient of SVDD is 0.01, and the kernel width is 300. From the perspective of accuracy and running time, the number of iterations of t-SNE, M-t-SNE, and AKMD-t-SNE were all selected as 300, and thePerplexity was set as 30. The time delay of DLLE and AKMD-t-SNE is 10; The confidence level of the fault detection control limit is 99%.Figures3 and 4 display the fault detection plots for fault type 1 and fault type 2, respectively. As can be seen from Figure 3, the FS of the four methods can detect the occurrence of faults to a certain extent. The AKMD-t-SNE-FS-SVDD method can distinguish between normal data and abnormal data more obviously and exhibits a better detection effect over those of the alternative models. After the addition of disturbance, the statistics of the four methods all changed to varying degrees in Figure 4. The proposed method achieves superior detection results in both FS and RS.Figure 3
Numerical case fault type 1 detection results of different methods: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)Figure 4
Numerical case fault type 2 detection results of different methods: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)Generally, it is considered that the detection of six or more continuous statistics exceeding the control limit indicates that the fault is intrinsically detected. The first sampling instant that fault is intrinsically detected is denoted asDetectionTime, and the time at which the disturbance was added is denoted as OccurTime. Delay is defined as(45)Delay=DetectionTime−OccurTime.Table1 shows the results of fault detection performance indicators for numerical case, including the miss alarm rate (MAR), detection delay, and average false alarm rate (FAR) of the two types of faults. The detection delay is defined by the earlier Delay between the FS and the RS. Under each fault type in Table 1, the position where the method proposed in this paper obtains the most excellent detection performance indexes are marked in bold. Compared with other methods, the AKMD-t-SNE-SVDD method can capture the information of anomaly more rapidly and obtains the lowest MAR.Table 1
Results of numerical case fault detection performance indicators.
Fault typeMAR (%) and delay in numerical caseDLLE-SVDDt-SNE-SVDDM-t-SNE-SVDDAKMD-t-SNE-SVDDFSRSDelayFSRSDelayFSRSDelayFSRSDelay140.6399.0011324.2596.25152.3896.6300.598.884282.132.38068.1343.3854056.6351.2549421.382.380Average FAR (%) in numerical case11.2541.56—0.2520.75—9.0019.75—4.7541.50—After the system is detected in an abnormal state, fault isolation is realized through HMI and HCC, respectively, based on the four methods, and the sampling interval isτ∈1,500. Figures 5 and 6 are fault isolation bar charts for type 1 and type 2, respectively. The reason for the occurrence of fault type 1 is that variable x2 is added with step disturbance. In Figure 5, the isolation results of the DLLE-SVDD method and the AKMD-t-SNE-SVDD method fit the fault description. Specifically, in the AKMD-t-SNE-SVDD-HCC method, the probability that x2 is judged to be the primary responsible variable exceeds 60%. Fault type 2 occurs because a ramp signal is added to variable x3. Since variable x5 is the related variable of x3, x3 is the main variable causing the fault and x5 is the secondary responsible variable. The method in Figure 6(c) cannot isolate the variables which caused the change. x3 and x5 can be isolated in Figures 6(a), 6(b), and 6(d), where Figure 6(d) is more consistent with the fault description. Compared with the HMI isolation method, the ©rrelevant variables of the HCC method have less impact on the isolation results, which verifies the effectiveness of the proposed method.Figure 5
Numerical case fault type 1 isolation bar charts of different methods based on HMI and HCC: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)Figure 6
Numerical case fault type 2 isolation bar charts of different methods based on HMI and HCC: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)


### 3.2. Multimode Tennessee Eastman Process

TE Process (TEP) is extensively researched in process monitoring owing to its complexity. The flow chart of TEP is shown in Figure7, which is composed of five major units, including the reactor, condenser, gas-liquid separator, stripper, and compressor [36]. A total of 6 working conditions are set in the multimode TE Process (MTEP), and 20 fault conditions are set for each mode. The detailed descriptions of 20 types of abnormal disturbances are shown in Table 2 [37]. There are 41 measurement variables and 11 manipulated variables in the TEP, among which the measurement variables are composed of 22 continuous measurement variables and 19 component variables. The stripper steam valve in the manipulated variables is constant at 1 under any mode, thus, 22 continuous measurement variables and 10 manipulated variables are selected to achieve modeling and monitoring. MTEP of the MATLAB simulation model can be downloaded from https://depts.washington.edu/control/LARRY/TE/download.html.Figure 7
TEP flowchart.Table 2
Disturbances in the TE process.
No.Variable descriptionTypeIDV (1)A/C feed ratio, B composition constantStepIDV (2)B composition, A/C feed ratio constantStepIDV (3)D feed temperatureStepIDV (4)Reactor cooling water inlet temperatureStepIDV (5)Condenser cooling water inlet temperatureStepIDV (6)A feed lossStepIDV (7)C Header pressure loss-reduced availabilityStepIDV (8)A, B, C feed compositionRandomIDV (9)D feed temperatureRandomIDV (10)C feed temperatureRandomIDV (11)Reactor cooling water inlet temperatureRandomIDV (12)Condenser cooling water inlet temperatureRandomIDV (13)Reaction kineticsSlow driftIDV (14)Reaction kineticsStickingIDV (15)Condenser cooling water valveStickingIDV (16)UnknownUnknownIDV (17)UnknownUnknownIDV (18)UnknownUnknownIDV (19)UnknownUnknownIDV (20)UnknownUnknownWhen the system is running normally, 300 samples are collected in each mode, and a total of 1800 samples in six modes constitute the training data set. The abnormal data of mode 1 and mode 2 were selected as testing data, with 1000 samples for each disturbance, and abnormal disturbance was added from the 201st sampling instant. As shown in Figure8, the variables v1, v2, and v3 of the training data were selected to draw a three-dimensional scatter diagram.Figure 8
MTEP 3D scatter diagram.The four models of the numerical case are also applied to fault detection of MTEP. According to historical experience, the dimensionality reduction of the four algorithms is set to 15, the iteration times of t-SNE, M-t-SNE, and AKMD-t-SNE are all selected as 500. Other parameters remain the same as in the above numerical case.Figures9 and 10 show the fault detection plots of IDV (3) in mode 1 and IDV (12) in mode 2, respectively. In Figure 9, DLLE and t-SNE methods cannot detect the occurrence of the disturbance, while the M-t-SNE method can briefly identify system anomalies, only AKMD-t-SNE-SVDD can effectively detect a persistent anomaly. As shown in Figure 10, the SVDD distances of the four methods have changed to different degrees when the disturbance occurs. Compared with the other three methods, AKMD-t-SNE-FS-SVDD has obtained a higher fault detection rate.Figure 9
MTEP mode 1 IDV (3) fault detection results of different methods: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)Figure 10
MTEP mode 2 IDV (12) fault detection results of different methods: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)Tables3 and 4, respectively, show the results of the fault detection performance indicators for MTEP mode 1 and mode 2, where △ indicates that the fault is not truly detected. Compared with other methods, the method proposed in this paper has obvious fault detection performance advantages in multiple fault types.Table 3
Results of fault detection performance indicators for MTEP mode 1.
IDVMAR (%) and delay in MTEP mode 1DLLE-SVDDt-SNE-SVDDM-t-SNE-SVDDAKMD-t-SNE-SVDDFSRSDelayFSRSDelayFSRSDelayFSRSDelay13.633.25261.382.88131.633.00142.383.2519229.6341.13326.7538.00272.8841.50263.8841.13323100.00100.00△96.2599.88△87.3899.882282.88100.00140.0098.63099.6397.7540.0098.1300.0098.6305100.00100.00△96.6397.25△97.6399.75△24.75100.0010060.000.7500.000.1300.000.5000.000.75070.0084.5000.1385.5000.0085.0000.0084.500815.1326.6311420.5028.3811612.0027.3810712.3826.631029100.00100.00△92.1396.007284.7599.388614.38100.00191025.63100.0017136.5099.2518723.0099.1316816.75100.001071111.8860.502468.6365.502511.0060.75233.2560.63211299.8874.757978.7562.756967.3877.258027.7575.00791316.1312.5010012.5012.139811.0012.009312.5012.501001417.0099.507171.2593.881177.6396.75657.6399.506615100.00100.00△100.0098.88△88.6399.2522832.75100.004116100.00100.00△100.0097.88△87.1399.6322818.13100.0022173.8878.252413.3878.00403.1380.13252.1378.13171871.8895.1346498.3893.3868560.1396.1322824.8895.251071912.0093.136455.2591.631514.1392.75432.7593.13212026.7557.6316530.8847.5023819.0057.7516710.7557.6342Average FAR (%) in MTEP mode 11.051.53—3.301.83—2.530.68—10.481.18—Table 4
Results of fault detection performance indicators for MTEP mode 2.
IDVMAR (%) and delay in MTEP mode 2DLLE-SVDDt-SNE-SVDDM-t-SNE-SVDDAKMD-t-SNE-SVDDFSRSDelayFSRSDelayFSRSDelayFSRSDelay14.002.13173.252.00231.251.50141.252.1310215.382.251815.501.88143.001.63144.382.25183100.0098.25123100.0096.50△83.7589.0019317.8898.2510240.1392.001100.0097.3860.0086.1300.0092.1305100.0098.6325899.5098.00△87.0090.6319427.2598.75560.000.1300.000.2500.000.0000.000.13070.0079.2500.1377.6300.0076.0000.0079.380818.7512.7510220.6324.3811813.5013.259912.8812.751029100.0098.38123100.0098.88△84.6388.7511832.0098.381031026.1396.6312350.0091.1321022.8886.0016214.1396.631051122.6333.882559.0071.882812.8842.13234.6333.88231299.6361.003087.0084.002871.3864.752715.6361.00181314.6310.008010.0010.638510.139.008012.7510.0080143.0043.502168.8888.88△2.5063.25200.8843.251215100.0096.88122100.0091.2511886.6389.2511817.8896.881131699.8896.1390100.00100.00△94.3887.6322820.5096.00891714.0074.632632.6392.13373.1377.50252.7574.50201871.5091.2512390.1396.6364545.7585.3819116.0091.381141920.8883.886754.8883.501188.3877.88553.0083.88242026.0022.7512227.8840.3810419.3827.1311812.5022.75100Average FAR (%) in MTEP mode 20.000.03—2.632.03—0.337.03—5.080.03—After the system anomaly is detected, fault isolation is realized through HMI and HCC respectively, and the data sampling interval isτ∈1,500. Figure 11 shows the fault isolation results of IDV (3) in mode 1. The fault of IDV (3) is described as the step disturbance added to the D feed temperature, which directly affects the change of the D feed flow valve (v23) under the regulation of the control system. The change of feed D causes the change of E feed flow valve (v24) and A feed flow valve (v25) at the same time, and the change of feed also affects the relevant variables of the upstream stripper and downstream separator (v13, v16, v29, v30). The main variables isolated by AKMD-t-SNE-SVDD-HCC conform to the fault description. Figure 12 shows the comparison of the trajectories of the three main isolation variables of IDV (3) in mode 1 under normal and abnormal conditions.Figure 11
MTEP mode 1 IDV (3) fault isolation bar charts of different methods based on HMI and HCC: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)Figure 12
Isolated variable trajectories of IDV (3) in mode 1: (a)v30, (b) v29, and (c) v23.
(a)(b)(c)Figure13 shows the isolation results of IDV (12) in mode 2. The reason for the occurrence of IDV (12) is that the condenser cooling water inlet temperature adds random disturbance. As the condenser acts directly on the output product of the reactor, it will directly affect the variable of reactor pressure (v7). The change of the output product of the reactor will also cause the alteration of the upstream stripper pressure (v16). The gas-liquid separator pressure (v13) in the downstream unit of the condenser is also affected, which causes the fluctuation of the separator temperature (v11), the purge valve (v28), and the purge rate (v10). Thus, the isolated responsible variables by the proposed scheme are consistent with the fault description of IDV (12). Figure 14 shows the comparison of the trajectories of the three major isolation variables of IDV (12) in mode 2 under the normal and abnormal conditions. The fault isolation results of the proposed method outperform the other methods in the above-given experiments.Figure 13
MTEP mode 2 IDV (12) fault isolation bar charts of different methods based on HMI and HCC: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)Figure 14
Isolated variable trajectories of IDV (12) in mode 2: (a)v16, (b) v13, and (c) v7.
(a)(b)(c)Table5 shows the fault isolation results of mode 1 and mode 2 based on the HCC method. The table presents the three most important isolation variables and HCC values of each type of fault disturbance. Due to the different working conditions, there are some differences in the main variables of isolation under the two modes, and the isolated responsible variables in the two modes basically correspond to the disturbance description.Table 5
Fault isolation results of MTEP mode 1 and mode 2 based on HCC (%).
IDVMode 1Mode 2VHCCVHCCVHCCVHCCVHCCVHCC1117.52187.22297.06268.18287.15296.8922817.311016.982113.97289.86109.86209.7933021.562915.732312.762820.782517.96117.0243177.80296.37306.043168.25288.37107.2752318.683013.112411.25118.622518.241814.1962517.28116.771011.25119.282518.772110.7772639.821821.062015.562636.001815.942015.598188.23267.98117.79268.28117.25187.2092323.932411.261310.542816.792516.36115.05101849.313010.53235.631829.912513.002811.64113162.42215.8295.503165.88214.8974.73121318.81718.781617.781618.271316.41714.55132812.791012.75219.953012.412910.19167.8014925.633124.002120.383123.872121.67921.37152317.102415.163011.262519.442818.59117.49162316.323013.782912.392825.202518.39115.12173136.072130.32911.452145.493137.4693.93183016.212314.392410.732521.102819.45118.49193052.742929.00244.663026.022822.422919.39201612.321812.321312.101610.671310.421810.31


## 3.1. Numerical Case

To testify the effectiveness of the proposed method, a numerical dynamic multimode process is simulated. The dynamic multimode model is created as follows:(43)x1=0.6352s1+0.1254s2+e1,x2=0.6951s1+0.1557ln1+s22+e2,x3=0.5028s1s2+0.8312ln2+coss2+e3,x4=0.4152s1+0.2617s22+e4,x5=0.0569x3+0.7962sinx4+e5,where s1 and s2 are latent variables, e1,e2,e3,e4,e5T are independent zero-mean Gaussian white noises with and a standard deviation of 0.01. Variable x5 is the correlation variable of variable x3 and variable x4. Set three modes for the latent variable as(44)Mode1:s1τ:Uniform−5,−1s2τ:N−2,2,Mode2:s1τ=0.8cos2s1τ−1+5esinτ−1s2τ=1.5coses2τ−1+ln1+s22τ−1,Mode3:s1τ=0.1sin2s1τ−1+2.8ecos2s1τ−1s2τ=3.25e0.3sins2τ−1+0.1coss2τ−1,where τ represents the index of the current sampling instant. The initial state of mode 2 and mode 3 is 0, i.e., s11=s21=0.Five hundred samples were generated in each mode, that is, a total of 1500 samples were generated in the three modes to form the training data. In order to show the spatial distribution characteristics of the data of the three modes intuitively, Figure 2 shows the scatter plot of variables x2, x3, and x4 in three modes.Figure 2
Numerical case 3D scatter plot.Set two types of fault:Type 1: the system operates normally for 200 sampling instants in mode 2. At the 201st sampling instant (τ=201), a step signal with an amplitude of +1.35 is added to variable x2, and 800 sampling instants are continued in this state.Type 2: the system operates normally for 200 sampling instants in mode. At the 201st sampling instant (τ=201), the ramp signal of 0.03τ−200 is added to variable 3 and continues to operate for 800 sampling instants in this state.Each of the two fault types has 1000 test samples. The augment matrix-based dynamic locally linear embedding (DLLE), t-SNE, t-SNE based on Mahaobanobis distance (M-t-SNE), and AKMD-t-SNE are used to extract the features of the established dynamic multimode numerical model. SVDD statistics are established in all the four methods for fault detection. According to experience, the dimension reduction of the four methods is 2, the nearest neighbor of DLLE is 25, the penalty coefficient of SVDD is 0.01, and the kernel width is 300. From the perspective of accuracy and running time, the number of iterations of t-SNE, M-t-SNE, and AKMD-t-SNE were all selected as 300, and thePerplexity was set as 30. The time delay of DLLE and AKMD-t-SNE is 10; The confidence level of the fault detection control limit is 99%.Figures3 and 4 display the fault detection plots for fault type 1 and fault type 2, respectively. As can be seen from Figure 3, the FS of the four methods can detect the occurrence of faults to a certain extent. The AKMD-t-SNE-FS-SVDD method can distinguish between normal data and abnormal data more obviously and exhibits a better detection effect over those of the alternative models. After the addition of disturbance, the statistics of the four methods all changed to varying degrees in Figure 4. The proposed method achieves superior detection results in both FS and RS.Figure 3
Numerical case fault type 1 detection results of different methods: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)Figure 4
Numerical case fault type 2 detection results of different methods: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)Generally, it is considered that the detection of six or more continuous statistics exceeding the control limit indicates that the fault is intrinsically detected. The first sampling instant that fault is intrinsically detected is denoted asDetectionTime, and the time at which the disturbance was added is denoted as OccurTime. Delay is defined as(45)Delay=DetectionTime−OccurTime.Table1 shows the results of fault detection performance indicators for numerical case, including the miss alarm rate (MAR), detection delay, and average false alarm rate (FAR) of the two types of faults. The detection delay is defined by the earlier Delay between the FS and the RS. Under each fault type in Table 1, the position where the method proposed in this paper obtains the most excellent detection performance indexes are marked in bold. Compared with other methods, the AKMD-t-SNE-SVDD method can capture the information of anomaly more rapidly and obtains the lowest MAR.Table 1
Results of numerical case fault detection performance indicators.
Fault typeMAR (%) and delay in numerical caseDLLE-SVDDt-SNE-SVDDM-t-SNE-SVDDAKMD-t-SNE-SVDDFSRSDelayFSRSDelayFSRSDelayFSRSDelay140.6399.0011324.2596.25152.3896.6300.598.884282.132.38068.1343.3854056.6351.2549421.382.380Average FAR (%) in numerical case11.2541.56—0.2520.75—9.0019.75—4.7541.50—After the system is detected in an abnormal state, fault isolation is realized through HMI and HCC, respectively, based on the four methods, and the sampling interval isτ∈1,500. Figures 5 and 6 are fault isolation bar charts for type 1 and type 2, respectively. The reason for the occurrence of fault type 1 is that variable x2 is added with step disturbance. In Figure 5, the isolation results of the DLLE-SVDD method and the AKMD-t-SNE-SVDD method fit the fault description. Specifically, in the AKMD-t-SNE-SVDD-HCC method, the probability that x2 is judged to be the primary responsible variable exceeds 60%. Fault type 2 occurs because a ramp signal is added to variable x3. Since variable x5 is the related variable of x3, x3 is the main variable causing the fault and x5 is the secondary responsible variable. The method in Figure 6(c) cannot isolate the variables which caused the change. x3 and x5 can be isolated in Figures 6(a), 6(b), and 6(d), where Figure 6(d) is more consistent with the fault description. Compared with the HMI isolation method, the ©rrelevant variables of the HCC method have less impact on the isolation results, which verifies the effectiveness of the proposed method.Figure 5
Numerical case fault type 1 isolation bar charts of different methods based on HMI and HCC: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)Figure 6
Numerical case fault type 2 isolation bar charts of different methods based on HMI and HCC: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)


## 3.2. Multimode Tennessee Eastman Process

TE Process (TEP) is extensively researched in process monitoring owing to its complexity. The flow chart of TEP is shown in Figure7, which is composed of five major units, including the reactor, condenser, gas-liquid separator, stripper, and compressor [36]. A total of 6 working conditions are set in the multimode TE Process (MTEP), and 20 fault conditions are set for each mode. The detailed descriptions of 20 types of abnormal disturbances are shown in Table 2 [37]. There are 41 measurement variables and 11 manipulated variables in the TEP, among which the measurement variables are composed of 22 continuous measurement variables and 19 component variables. The stripper steam valve in the manipulated variables is constant at 1 under any mode, thus, 22 continuous measurement variables and 10 manipulated variables are selected to achieve modeling and monitoring. MTEP of the MATLAB simulation model can be downloaded from https://depts.washington.edu/control/LARRY/TE/download.html.Figure 7
TEP flowchart.Table 2
Disturbances in the TE process.
No.Variable descriptionTypeIDV (1)A/C feed ratio, B composition constantStepIDV (2)B composition, A/C feed ratio constantStepIDV (3)D feed temperatureStepIDV (4)Reactor cooling water inlet temperatureStepIDV (5)Condenser cooling water inlet temperatureStepIDV (6)A feed lossStepIDV (7)C Header pressure loss-reduced availabilityStepIDV (8)A, B, C feed compositionRandomIDV (9)D feed temperatureRandomIDV (10)C feed temperatureRandomIDV (11)Reactor cooling water inlet temperatureRandomIDV (12)Condenser cooling water inlet temperatureRandomIDV (13)Reaction kineticsSlow driftIDV (14)Reaction kineticsStickingIDV (15)Condenser cooling water valveStickingIDV (16)UnknownUnknownIDV (17)UnknownUnknownIDV (18)UnknownUnknownIDV (19)UnknownUnknownIDV (20)UnknownUnknownWhen the system is running normally, 300 samples are collected in each mode, and a total of 1800 samples in six modes constitute the training data set. The abnormal data of mode 1 and mode 2 were selected as testing data, with 1000 samples for each disturbance, and abnormal disturbance was added from the 201st sampling instant. As shown in Figure8, the variables v1, v2, and v3 of the training data were selected to draw a three-dimensional scatter diagram.Figure 8
MTEP 3D scatter diagram.The four models of the numerical case are also applied to fault detection of MTEP. According to historical experience, the dimensionality reduction of the four algorithms is set to 15, the iteration times of t-SNE, M-t-SNE, and AKMD-t-SNE are all selected as 500. Other parameters remain the same as in the above numerical case.Figures9 and 10 show the fault detection plots of IDV (3) in mode 1 and IDV (12) in mode 2, respectively. In Figure 9, DLLE and t-SNE methods cannot detect the occurrence of the disturbance, while the M-t-SNE method can briefly identify system anomalies, only AKMD-t-SNE-SVDD can effectively detect a persistent anomaly. As shown in Figure 10, the SVDD distances of the four methods have changed to different degrees when the disturbance occurs. Compared with the other three methods, AKMD-t-SNE-FS-SVDD has obtained a higher fault detection rate.Figure 9
MTEP mode 1 IDV (3) fault detection results of different methods: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)Figure 10
MTEP mode 2 IDV (12) fault detection results of different methods: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)Tables3 and 4, respectively, show the results of the fault detection performance indicators for MTEP mode 1 and mode 2, where △ indicates that the fault is not truly detected. Compared with other methods, the method proposed in this paper has obvious fault detection performance advantages in multiple fault types.Table 3
Results of fault detection performance indicators for MTEP mode 1.
IDVMAR (%) and delay in MTEP mode 1DLLE-SVDDt-SNE-SVDDM-t-SNE-SVDDAKMD-t-SNE-SVDDFSRSDelayFSRSDelayFSRSDelayFSRSDelay13.633.25261.382.88131.633.00142.383.2519229.6341.13326.7538.00272.8841.50263.8841.13323100.00100.00△96.2599.88△87.3899.882282.88100.00140.0098.63099.6397.7540.0098.1300.0098.6305100.00100.00△96.6397.25△97.6399.75△24.75100.0010060.000.7500.000.1300.000.5000.000.75070.0084.5000.1385.5000.0085.0000.0084.500815.1326.6311420.5028.3811612.0027.3810712.3826.631029100.00100.00△92.1396.007284.7599.388614.38100.00191025.63100.0017136.5099.2518723.0099.1316816.75100.001071111.8860.502468.6365.502511.0060.75233.2560.63211299.8874.757978.7562.756967.3877.258027.7575.00791316.1312.5010012.5012.139811.0012.009312.5012.501001417.0099.507171.2593.881177.6396.75657.6399.506615100.00100.00△100.0098.88△88.6399.2522832.75100.004116100.00100.00△100.0097.88△87.1399.6322818.13100.0022173.8878.252413.3878.00403.1380.13252.1378.13171871.8895.1346498.3893.3868560.1396.1322824.8895.251071912.0093.136455.2591.631514.1392.75432.7593.13212026.7557.6316530.8847.5023819.0057.7516710.7557.6342Average FAR (%) in MTEP mode 11.051.53—3.301.83—2.530.68—10.481.18—Table 4
Results of fault detection performance indicators for MTEP mode 2.
IDVMAR (%) and delay in MTEP mode 2DLLE-SVDDt-SNE-SVDDM-t-SNE-SVDDAKMD-t-SNE-SVDDFSRSDelayFSRSDelayFSRSDelayFSRSDelay14.002.13173.252.00231.251.50141.252.1310215.382.251815.501.88143.001.63144.382.25183100.0098.25123100.0096.50△83.7589.0019317.8898.2510240.1392.001100.0097.3860.0086.1300.0092.1305100.0098.6325899.5098.00△87.0090.6319427.2598.75560.000.1300.000.2500.000.0000.000.13070.0079.2500.1377.6300.0076.0000.0079.380818.7512.7510220.6324.3811813.5013.259912.8812.751029100.0098.38123100.0098.88△84.6388.7511832.0098.381031026.1396.6312350.0091.1321022.8886.0016214.1396.631051122.6333.882559.0071.882812.8842.13234.6333.88231299.6361.003087.0084.002871.3864.752715.6361.00181314.6310.008010.0010.638510.139.008012.7510.0080143.0043.502168.8888.88△2.5063.25200.8843.251215100.0096.88122100.0091.2511886.6389.2511817.8896.881131699.8896.1390100.00100.00△94.3887.6322820.5096.00891714.0074.632632.6392.13373.1377.50252.7574.50201871.5091.2512390.1396.6364545.7585.3819116.0091.381141920.8883.886754.8883.501188.3877.88553.0083.88242026.0022.7512227.8840.3810419.3827.1311812.5022.75100Average FAR (%) in MTEP mode 20.000.03—2.632.03—0.337.03—5.080.03—After the system anomaly is detected, fault isolation is realized through HMI and HCC respectively, and the data sampling interval isτ∈1,500. Figure 11 shows the fault isolation results of IDV (3) in mode 1. The fault of IDV (3) is described as the step disturbance added to the D feed temperature, which directly affects the change of the D feed flow valve (v23) under the regulation of the control system. The change of feed D causes the change of E feed flow valve (v24) and A feed flow valve (v25) at the same time, and the change of feed also affects the relevant variables of the upstream stripper and downstream separator (v13, v16, v29, v30). The main variables isolated by AKMD-t-SNE-SVDD-HCC conform to the fault description. Figure 12 shows the comparison of the trajectories of the three main isolation variables of IDV (3) in mode 1 under normal and abnormal conditions.Figure 11
MTEP mode 1 IDV (3) fault isolation bar charts of different methods based on HMI and HCC: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)Figure 12
Isolated variable trajectories of IDV (3) in mode 1: (a)v30, (b) v29, and (c) v23.
(a)(b)(c)Figure13 shows the isolation results of IDV (12) in mode 2. The reason for the occurrence of IDV (12) is that the condenser cooling water inlet temperature adds random disturbance. As the condenser acts directly on the output product of the reactor, it will directly affect the variable of reactor pressure (v7). The change of the output product of the reactor will also cause the alteration of the upstream stripper pressure (v16). The gas-liquid separator pressure (v13) in the downstream unit of the condenser is also affected, which causes the fluctuation of the separator temperature (v11), the purge valve (v28), and the purge rate (v10). Thus, the isolated responsible variables by the proposed scheme are consistent with the fault description of IDV (12). Figure 14 shows the comparison of the trajectories of the three major isolation variables of IDV (12) in mode 2 under the normal and abnormal conditions. The fault isolation results of the proposed method outperform the other methods in the above-given experiments.Figure 13
MTEP mode 2 IDV (12) fault isolation bar charts of different methods based on HMI and HCC: (a) DLLE-SVDD, (b) t-SNE-SVDD, (c) M-t-SNE-SVDD, and (d) AKMD-t-SNE-SVDD.
(a)(b)(c)(d)Figure 14
Isolated variable trajectories of IDV (12) in mode 2: (a)v16, (b) v13, and (c) v7.
(a)(b)(c)Table5 shows the fault isolation results of mode 1 and mode 2 based on the HCC method. The table presents the three most important isolation variables and HCC values of each type of fault disturbance. Due to the different working conditions, there are some differences in the main variables of isolation under the two modes, and the isolated responsible variables in the two modes basically correspond to the disturbance description.Table 5
Fault isolation results of MTEP mode 1 and mode 2 based on HCC (%).
IDVMode 1Mode 2VHCCVHCCVHCCVHCCVHCCVHCC1117.52187.22297.06268.18287.15296.8922817.311016.982113.97289.86109.86209.7933021.562915.732312.762820.782517.96117.0243177.80296.37306.043168.25288.37107.2752318.683013.112411.25118.622518.241814.1962517.28116.771011.25119.282518.772110.7772639.821821.062015.562636.001815.942015.598188.23267.98117.79268.28117.25187.2092323.932411.261310.542816.792516.36115.05101849.313010.53235.631829.912513.002811.64113162.42215.8295.503165.88214.8974.73121318.81718.781617.781618.271316.41714.55132812.791012.75219.953012.412910.19167.8014925.633124.002120.383123.872121.67921.37152317.102415.163011.262519.442818.59117.49162316.323013.782912.392825.202518.39115.12173136.072130.32911.452145.493137.4693.93183016.212314.392410.732521.102819.45118.49193052.742929.00244.663026.022822.422919.39201612.321812.321312.101610.671310.421810.31


## 4. Conclusions

Process monitoring is of great significance to the safe and stable production of the chemical process industry. Most of the research studies on process monitoring mainly focus on the static single-mode process. In this paper, based on the t-SNE algorithm, an AKMD-t-SNE method is proposed to realize feature extraction for dynamic multimode chemical process, and the SVDD method is used to construct statistics for non-Gaussian data fault detection. After the system anomaly is detected, fault isolation is achieved by the HCC method. The results of the numerical case and multimode TE process demonstrate that the proposed method is effective for fault detection and fault isolation in a dynamic multimode chemical process.



---

*Source: 8460463-2022-12-29.xml*