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 named modified 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.
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 called t-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, where represents the number of variables and represents the number of observations. Suppose the current sampling instant is and the sample at is . and past values with orders of data dynamics forms an augmented matrix:
In equation (1), the sampling time must be larger than the time delay , and the form of the augmented matrix formed is as follows:when the sampling time is smaller than the time delay , the form of the augmented matrix is as follows:
Then, the final augmented matrix formed by the training data as :
The AKMD-t-SNE algorithm is used to map the augmented high-dimensional data to the low-dimensional feature subspace , where is a matrix with observations and dimensions. is the matrix with observations and dimensions.
The Mahalanobis distance could eliminate the discrepancies among different modes of training data. are two samples in the high-dimensional data. The Mahalanobis distance between high-dimensional samples and is calculated aswhere is the covariance matrix of .
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 , the augmented vector is mapped into a higher-dimensional feature space. Denote the mapped matrix as . The augmented kernel Mahalanobis distance is determined aswhere , represents the covariance matrix of :where satisfies , is a symmetric and idempotent matrix:
The dot product of mapped vectors is calculated by the Gaussian kernel function:
Then, the centered kernel matrix is defined as
A full rank factorization via singular value decomposition (SVD) method is analyzed:where and are both orthogonal matrices.
From (7), we have
Then, its inverse matrix is calculated:
The centered kernel matrix and its Moore-Penrose pseudoinverse matrix are given as
Combining (15) with (17), we obtain
The augmented kernel Mahalanobis distance can be rewritten aswhere .
The conditional probability is calculated between samples in higher-dimensional space. The conditional probability of sample as the nearest neighbor of in higher-dimensional space is denoted as :where represents the variance of the Gaussian distribution with as the center, and represents one of the nearest neighbors of .
The selection of value is related to the degree of Perplexity, which is defined aswhere is the Shannon entropy of :
Assuming that the mapping samples of higher-dimensional observations and in the low-dimensional space are and , and the variance of the Gaussian distribution with as the center is set to , the corresponding conditional probability that is the nearest sample of denoted as :where represents one of the nearest neighbors of .
Kullback–Leibler Divergence is used to measure the similarity between the two individual distributions. The similarity cost function as:
The cost function is minimized by gradient descent, and the gradient of the cost function is obtained:
The conditional probability is replaced by the joint probability distribution, such that for any and , , . The following symmetric method is adopted to simplify the objective function:
The cost function is transformed into the following expression:
In order to solve the crowding problem of low-dimensional data in classification, t-distribution is used to replace Gaussian distribution. Redefine with a t-distribution of 1 degree of freedom, it follows that:
The gradient formula is converted to the following equation:
To avoid falling into the local optimum in the gradient descent process, a momentum term is added in the gradient update, as follows:where represents the solution of the th iteration, represents the momentum of the th iteration, and represents the learning efficiency.
The low-dimensional feature subspace 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 , which can be obtained from the Moore–Penrose pseudoinverse matrix of :
The feature subspace is recalculated via the mapping matrix and the residual subspace is denoted as :where is the high-dimensional projection of low-dimensional data and is the Moore–Penrose pseudoinverse of .
For testing data , the augmented matrix 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 :
The residual of is given as
2.2. Fault Detection Method
After dimension-reduction by the AKMD-t-SNE scheme, the low-dimensional data will be expressed as t-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 and radius , 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:where is the th support vector of the FS; is the slack variable to improve the fault tolerance of the SVDD model; is the penalty factor, the size of the reflects the extent of choosing outliers for a data set.
To achieve fault detection on the feature of online data () is to calculate the distance between the sample and spherical center and compare it with .where is one of the support vectors in the FS of training data, and is the th and th support vector in the FS of training data.
For the residual of online data (), SVDD statistics () and control limit () 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 variables and , and the MI of and is defined as :where and are the marginal distribution of and respectively, and are the joint distribution of and .
Once the abnormal is detected in the augmented online data, the augmented online data containing normal data and abnormal data is denoted as , where represents the number of samples and represents the number of variables. Since the augmented matrix contains the data of the previous sampling instants, is composed of submatrices , and the expression is given as
The FS statistics of constitute the vector . Compute the MI between and each vector in each submatrix as
Compute the MI between each vector in the online data and as
On this basis, a fault isolation method based on correlation coefficient (CC) is proposed. In order to understand the correlation level of variables expressed by , the MI between and itself is calculated as . The correlation coefficient of FS () is established as
The impact of each variable in the FS on the occurrence of disturbances was judged by the size of .
Similarly, the variable composed of RS statistics is denoted as , and the CC of RS () is calculated in the same way to achieve fault isolation.
The hybrid mutual information (HMI) and the hybrid correlation coefficient (HCC) were defined as
Calculate the weight coefficient :
The process monitoring flow diagram of the proposed method is shown in Figure 1.
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:where and are latent variables, are independent zero-mean Gaussian white noises with and a standard deviation of 0.01. Variable is the correlation variable of variable and variable . Set three modes for the latent variable aswhere represents the index of the current sampling instant. The initial state of mode 2 and mode 3 is 0, i.e., .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 , , and in three modes.
Set two types of fault: Type 1: the system operates normally for 200 sampling instants in mode 2. At the 201st sampling instant (), a step signal with an amplitude of +1.35 is added to variable , 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 (), the ramp signal of 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 the Perplexity 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%.
Figures 3 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.
(a)
(b)
(c)
(d)
(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 as , and the time at which the disturbance was added is denoted as . is defined as
Table 1 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 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.
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 . 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 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 is judged to be the primary responsible variable exceeds 60%. Fault type 2 occurs because a ramp signal is added to variable . Since variable is the related variable of , is the main variable causing the fault and is the secondary responsible variable. The method in Figure 6(c) cannot isolate the variables which caused the change. and 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.
(a)
(b)
(c)
(d)
(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 Figure 7, 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.
When 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 Figure 8, the variables , , and of the training data were selected to draw a three-dimensional 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.
Figures 9 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Tables 3 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.
After the system anomaly is detected, fault isolation is realized through HMI and HCC respectively, and the data sampling interval is . 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 () under the regulation of the control system. The change of feed D causes the change of E feed flow valve () and A feed flow valve () at the same time, and the change of feed also affects the relevant variables of the upstream stripper and downstream separator (, , , ). 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
Figure 13 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 (). The change of the output product of the reactor will also cause the alteration of the upstream stripper pressure (). The gas-liquid separator pressure () in the downstream unit of the condenser is also affected, which causes the fluctuation of the separator temperature (), the purge valve (), and the purge rate (). 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
Table 5 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.
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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported in part by the National Nature Science Foundation of China (Grant no. NSFC; 61403256), and Nature Science Foundation of Shanghai (Grant no. 19ZR1455200).