Multimode Process Monitoring Method Based on Multiblock Projection Nonnegative Matrix Factorization

A multimode process monitoring method based on multiblock projection nonnegative matrix factorization (MPNMF) is proposed for traditional process monitoring methods which often adopt global model of data and ignore local information of data. Firstly, the training data set of each mode is partitioned by the complete link algorithm and the multivariate data space is divided into several subblocks. Then, the projection nonnegative matrix factorization (PNMF) algorithm is used to model each subspace of each mode separately. A joint probabilistic statistic index is defined to identify the running modes of the process data. Finally, the Bayesian information criterion (BIC) is used to synthesize the statistics of each subblock and construct a new statistic for process monitoring. The proposed process monitoring method is applied to the TE process to verify its effectiveness.


Introduction
With the rapid development of computer technology, the chemical process has become more automatic and intelligent. In recent years, accidents have occurred frequently in the chemical process, safety in the production process has become crucial, and process monitoring has received more and more attention [1]. Multivariate statistical process monitoring (MSPM) methods, such as principal component analysis (PCA) and partial least squares (PLS), are widely used in industrial processes due to their ability to extract effective characteristic information from process data and conduct process monitoring [2][3][4].
Traditional MSPM methods need to assume that the measurement data come from a single stable mode, while modern chemical enterprises have a variety of operating modes due to changes in product characteristics, set values, or raw material components. Therefore, it is assumed that the MSPM method under a single working condition cannot meet the requirements of multimode monitoring. This not only weakens the statistical characteristics of the process in different modes but also leads to inaccurate process performance analysis. Xu et al. proposed a PCA mixture model based on scale-incremental EM algorithm to realize multimode process monitoring [5]. Feital et al. proposed a multimodal modeling and monitoring method based on maximum likelihood principal component analysis and carried out modal identification [6]. Jiang et al. defined a joint probability statistical indictor to identify the running mode of the process [7]. Zhang et al. proposed a fully automatic offline pattern recognition method based on sliding window and k-means clustering [8]. Xiong et al. proposed a sliding window and differential algorithm to realize multicondition identification of industrial processes [9].
In addition to multimode processes, the modern chemical industry also contains a lot of non-Gaussian data. Therefore, the MSPM method assuming Gaussian conditions cannot meet monitoring requirements. Nonnegative matrix decomposition (NMF) is a new dimension reduction method. The advantage of this method is that the requirements for measurement data do not exceed nonnegative requirements so it is wider [10]. Yuan et al. proposed a projection nonnegative matrix factorization (PNMF) algorithm. The PNMF algorithm adds constraints to the NMF algorithm, reduces an unknown variable, and increases the sparsity of the matrix [11]. Ahmed et al. proposed an approximately optimal nonnegative matrix factor two-dimensional deconvolution (NMF2D) algorithm for separating underdetermined convolution integrals [12]. After that, they proposed a mono audio separation method based on multicomponent nonnegative matrix factor two-dimensional deconvolution (NMF2D) [13]. Woo et al. proposed an optimized complex nonnegative tensor factor two-dimensional deconvolution (CNTF2D) to separate sound sources mixed in an uncertain reverberation environment [14]. In recent decades, NMFbased methods have been used more and more widely. Gao et al. proposed a variational Bayesian subgroup adaptive sparse component extraction algorithm for thermal NDT & E imaging diagnosis [15]. They conducted experimental tests on artificial and natural defects of the robot and verified the effectiveness of the proposed method. Ahmed et al. proposed a thermal imaging CFRP defect detection method based on wavelet integral alternating sparse dictionary matrix factorization [16]. The above papers mainly study the thermal imaging diagnostic system for detecting product defects in the production process. Currently, there are few articles on NMF methods for industrial process fault monitoring.
Yang et al. proposed a method of fault diagnosis based on NMF and SVM [17]. They built an online monitoring model through NMF and then used SVM to build a fault classifier for further fault identification. Considering the imperfection of process data, Li et al. proposed a method based on robust nonnegative matrix projection (RNMP) to detect and diagnose faults in industrial processes [18]. Li et al. proposed an improved NMF (MNMF) method to extract potential variables in the process and combine it with process monitoring technology for fault detection [19]. However, the above process monitoring methods based on NMF only consider the global model of process data and ignore local information of the data, resulting in inaccurate monitoring results. To solve these problems, process data are usually divided into several subblocks, and then a monitoring model for each subblock is established separately. Jiang and Yan developed a multiblock principal component analysis method based on mutual information for fault detection and isolation [20]. Tong et al. introduced an improved multiblock principal component analysis (MBPCA) algorithm to extract the specificity of each block and the block fraction of correlation between different blocks [21]. Wang and Deng considered the local variable behavior of process data and proposed a multiblock PCA process monitoring method based on variable weight information [22]. Wang et al. proposed an adaptive partitioned nonnegative matrix factorization algorithm based on nonfixed subblock NMF model for fault monitoring in chemical processes [23]. Although they considered the local information of the data, they did not consider the mode identification of multimode processes. These multiblock monitoring methods have proved to be generally effective. However, multimode and non-Gaussian characteristics of the data should also be considered in process monitoring.
In this study, a multiblock projection nonnegative matrix factorization (MPNMF) algorithm for complex chemical processes in a multimode state is proposed. Firstly, the training data of each mode is divided into several subblocks by a complete link algorithm. Secondly, a PNMF model for each subblock is established for each mode, and a joint probability index is constructed to judge the current running mode. Finally, the statistics of each block in the recognition mode is combined by Bayesian information criterion to monitor whether abnormal events occur. The proposed method takes into account the multimode and non-Gaussian characteristics of process data, combines the advantages of partitioning idea and PNMF method, makes full use of local and global information of the data, enhances the sparsity of the data, and has good process monitoring performance.
The main contributions are as follows: (1) a multimode identification method based on joint probability is proposed, which can realize the identification of process operation mode. (2) A fault monitoring model based on multiblock PNMF algorithm is proposed. This paper uses the PNMF method to extract the main characteristics of the data and establish a monitoring model, which expands the application of PNMF in the field of process monitoring.

Preliminaries
This section briefly reviews the NMF and the PNMF algorithms, the process monitoring methods based on the PNMF method.
2.1. Nonnegative Matrix Factorization (NMF). As a new matrix factorization algorithm, in addition to nonnegative, the NMF method has no other requirements on the original data [10]. Since the decomposed matrix elements are only added and calculated, and the decomposed base matrix is sparse, NMF can learn the local features of the data and has better explanatory performance.
Supposing X = ½x 1 , x 2 , ⋯, x m ∈ R m×n is a nonnegative matrix, the core of NMF is to find W ∈ R m×k and H ∈ R k×n (belonging to a suitable set of nonnegative matrices) to satisfy the following formula: where m is the number of variables, n is the number of samples, k is the dimension of the low-dimensional space after dimensionality reduction, W is the base matrix, and H is the coefficient matrix. In addition, Eq. (1) should satisfy ðm + nÞk < mn.The solution of NMF can be summed up as a nonlinear optimization process. That is, an objective function is defined to describe the approximation of X and WH, and an appropriate set of iterative rules is found to solve the optimization problem. Lee and Seung proposed two kinds of objective functions to solve the optimization problem and gave iterative rules [10]. The first objective function "Euclidean distance" was selected in this paper.

Advances in Mathematical Physics
The objective function can be expressed as follows: The iterative rules of the above objective function are as follows: Under the above iterative rule, W retains the behavior information of the original matrix X, H is the lowdimensional approximation matrix of X. [11] proposed the PNMF algorithm. The basic idea of the algorithm is to replace the coefficient matrix H with W T X, which increases the sparsity of the matrix. PNMF improves the framework of NPF and reduces an unknown variable, which makes the solution complexity lower.

Projection Nonnegative Matrix Factorization (PNMF). In 2005, Yuan and Oja in
Let W T X = H, Eq. (1) can be expressed as follows: where the value of WW T should be infinitely close to the identity matrix. Similarly, Euclidean distance in Eq. (2) is used to measure the proximity of both sides of Eq. (5). The objective function of the PNMF algorithm is as follows: The iterative formula for the objective function is as follows [11]: At present, the initialization schemes often used are random value, singular value decomposition (SVD), and principal component analysis (PCA) [13,14,18]. Wang et al. pointed out that using the PCA's load matrix as the initial value of W will make NMF converge faster and minimize errors [18]. Therefore, this paper uses the method proposed by Wang et al. to obtain the initial value of NMF. The fault monitoring model of the PNMF is as follows: whereX is the estimation matrix of X, which contains the main feature information of matrix X.X is a residual matrix. H = W T X is a low-order approximation matrix. In general, the monitoring statistical indicators of PNMF are N 2 and SPE as follows [19]: withĥ, x,x, andx as the column vectors of the matricesĤ, X, X, andX, respectively. The N 2 statistic is the monitoring indicator of theX space, and the SPE statistic is the monitoring indicator of theX space. Some fault information are reflected in theX space, and some fault information are reflected in theX space. In order to monitor the fault more accurately, two monitoring statistics are usually monitored at the same time in many literatures [18][19][20][21][22][23]. N 2 lim and SPE lim are confidence limits of statistics N 2 and SPE. If statistic N 2 or SPE of the measurement sample exceeds the confidence limit, an abnormality is considered to have occurred and the system issues an alarm. Because the PNMF method does not assume that the data follow a certain fixed distribution, the kernel density estimation (KDE) can be used to calculate the confidence limits N 2 lim and SPE lim of statistics N 2 and SPE. The specific KDE algorithm can be referred to in [19].

Methodology
In this section, the proposed multimode process monitoring method based on MPNMF model will be introduced in detail. In this method, a complete link algorithm is used first to block the training data, then multiple PNMF models are established, and mode identification is performed using joint probability. Finally, two monitoring indicators are constructed using BIC to achieve fault detection.
3.1. Complete Link Algorithm. The complete link algorithm is a very popular method for calculating hierarchical clustering. Given a set X = ½x 1 , x 2 , ⋯, x m ∈ R m×n of variables, its core idea is to treat each set of variable data x i ði = 1, 2, ⋯, mÞ as a separate cluster. Then, merge the two closest clusters based on some distance metric to get the next cluster and repeat the merging process until only one cluster remains. In addition, the complete linking algorithm uses the distance between 3 Advances in Mathematical Physics the most distant sample pairs in two clusters as the distance between the two clusters [24].
The specific process of dividing subblocks based on the complete link algorithm is as follows: (1) Taking each variable data as a cluster, the matrix X is divided into m clusters (2) Calculating the distance dðx i , x j Þ between two different clusters x i ði = 1, 2, ⋯, mÞ and x j ðj = 1, 2, ⋯, mÞ.
is a common distance function, which is selected as the distance function of the complete link algorithm (3) Conducting distance matrix as following: Merging the two clusters with the smallest distance into the new cluster fx i , x j g and calculating the distance between fx i , x j g and the remaining clusters as follows: (4) Updating the distance matrix D and repeat Step (3) until the desired number of clusters is obtained 3.2. MPNMF Models. Supposing a set of historical normal data matrices X ∈ R m×n , where m and n are the number of variables and sample, and X can be divided into several subblocks using complete link algorithm. Therefore, the matrix X can be rewritten as: with B being the number of subblocks, X b ∈ R m b ×n , and m b the number of variables in the bth subblock. Then, the PNMF model of each subblock is established as follows: with W b ∈ R m b ×k b being the basis matrix of the bth subblock matrix decomposition, k b the dimension of the lowdimensional space after the dimensionality reduction of the b subblock, and the value of W b W T b is infinitely close to the identity matrix.
According to Eqs. (9) and (10), the monitoring statistical indicators N 2 b and SPE b of each PNMF model can be obtained as follows: with x b ∈ R m b ×1 as an measurement vector, N 2 b,lim and SPE b,lim are confidence limits of statistics N 2 b and SPE b . 3.3. Joint Probability for Mode Identification. Jiang and Yan used monitoring statistical indicators T 2 to identify the current sample mode and achieved good mode identification results [7]. Inspired by them, N 2 statistics with a large number of data features are used for mode identification in this paper.
For a current observation sample x new ∈ R m×1 , it is divided into B subblocks according to the partitioning results of the previous historical data, and get x new = ½x 1,new , x 2,new , ⋯,x B,new . Then, N 2 b,new is calculated for each subblock b in each corresponding mode. In order to identify the mode of the current subblock sample x b,new , the mode probability p is defined as follows [7]: with M i as the ith mode, and Þ, the probability of event x ∈ M i is calculated by the joint probability of each block as follows: If the observation sample x new ∈ R m×1 is assigned to a mismatch mode, the joint probability p is close to zero according to Jiang and Yan's study [7]. Therefore, it can be determined that the mode with the highest joint probability is the running mode of the current process.

Fault Monitoring Index Based on Bayesian Information
Criterion. After determining the mode of the current measurement sample x new , it is necessary to monitor it. Different subblocks have different monitoring statistics. To provide intuitive monitoring results, Bayesian information criterion (BIC) is often used to combine results with probabilistic methods [25]. The fault probability of x b,new corresponding to N 2 statistics is defined as follows: with F as the fault condition and N as the normal condition. The confidence level of p N 2 ðNÞ is αð0 ≤ α ≤ 1Þ and the 4 Advances in Mathematical Physics confidence level of p N 2 ðFÞ is 1 − α. Two conditional probabilities p N 2 ðx b,new jNÞ and p N 2 ðx b,new jFÞ are defined as: Using BIC to combine N 2 statistics in all subblocks to get the final statistic, as shown below: Similarly, using BIC to combine SPE statistic to get the final statistic, as shown below: The final statistical value within the confidence level α can be regarded as normal operation data.
3.5. The Proposed Method Implementation. In this paper, a multimode process monitoring method based on MPNMF is proposed. The method mainly includes offline modeling and online monitoring, as shown in Figures 1 and 2. The implementation of these two parts is as follows:

Simulation Study of TE Process
The Tennessee Eastman (TE) process simulation platform is developed by Eastman Chemical Company based on actual chemical processes. It can not only simulate the chemical production process and generate a large amount of normal data but also can obtain multiple fault data by setting multiple fault modes. Therefore, the TE process is often used to evaluate the effectiveness of various fault detection methods. As shown in Figure 3, the TE process is mainly composed of five operating units: a reactor, a condenser, a gas-liquid separator, a stripper, and a circulating compressor. A variety of control strategies have been proposed for the research of TE process. Ricker [26] proposed a decentralized control structure and gave six different operation modes according to the G/H ratio, as shown in Table 1, and the simulation data can be downloaded from http://depts.washington.edu/control/ LARRY/TE/download.html. Lyman and Georgakis [27] proposed a plant-wide control scheme, and simulation data can be downloaded from http://web.mit.edu/braatzgroup/ links.html.
Select mode 1 and mode 2 proposed by Richer and mode 3 proposed by Lyman and Georgakis to simulate multimodal processes. 33 variables (11 control variables, 22 measured variables) are selected for fault detection and diagnosis. The specific details of the variables can be found in [26].
Since the NMF-based method is very sensitive to the initial values of W ij and H ij , the selection of the initial values of W ij and H ij is important. At present, the initialization schemes often used are random value, singular value decomposition (SVD), and principal component analysis (PCA) [18,19,23]. Wang et al. pointed out that using the PCA's load matrix as the initial value of W ij will make NMF converge faster and minimize errors [23]. Therefore, this paper uses the method proposed by Wang et al. to obtain the initial value of NMF. 500 normal history samples were collected in each mode to build a monitoring model for each mode. In each mode, 21 faults are simulated, and the fault names are referenced in [27]. Each set of measurement data contains 960 samples, of which the first 160 samples are normal and the fault starts from the 161st sample and continues to the end of the process. The article uses 85% of the contribution of PCA to select the number of principal elements. The number of principal elements is used as the number of dimensionality reductions of all monitoring algorithms. The number of iterations of the NMF-based algorithm is 1000.

Advances in Mathematical Physics
All monitoring algorithms use kernel density estimation to calculate the control limits of statistical indicators, with a confidence level of 99%. In order to verify the validity of mode identification, two data sets were collected for testing, each of which contains 960 groups of samples, as shown below.
Case 1. Fault 0 in TE process represents a normal process. Normal data in different modes are used to test the performance of mode identification. In the system, samples range from 1 to 400 runs in mode 1, range from 401 to 560 runs in mode 2, and range from 561 to 960 runs in mode 3.
Case 2. Fault 1 in TE process is a step change in the A/C feed ratio, and the fault continues from the 161st point to the end of the process. The system runs from the 1st sample point to the 960th point in mode 2.
The joint probability of the three modes in Case 1. is shown in Figure 4(a). In Figure 4(a), it can be seen that as the operation mode changes, the joint probability will change. From point 1 to point 100, the joint probability of mode 1 is the highest, indicating that the system is operating in mode 1, followed by modes 2 and 3. The joint probability can be used to derive the operating mode of the system. The joint probability of the three modes in s2. is shown in Figure 4(b). Therefore, Figure 4(b) shows that the system is operating in mode 2.
Once the mode of the measurement sample is identified, it is necessary to monitor it. PCA method, NMF method, and PNMF method are selected as reference objects to verify the effectiveness of MPNMF method in fault monitoring. Fault     Table 2. The NMF-based algorithm used in Table 2 refers to the NMF software package published online by Professor Chih-Jen Lin of Tsinghua University in Taiwan (website: https://www.csie.ntu.edu.tw/ cjlin/nmf).
Compared with the results of PCA method, NMF method, and PNMF method, MPNMF method has the highest fault detection rate, especially in failure modes 3, 5, 10, 15, 20, and 21. This paper takes faults 5 and 10 as examples to illustrate the monitoring results of different methods.
Fault 5 is the step change of cooling water inlet temperature of condenser. Figure 5 shows the detection results of NMF, PNMF, and MPNMF methods. In Figure 5(a), it is easy to see that the fault detection rate of global NMF is lower. The   10 Advances in Mathematical Physics detection rates of N 2 statistics and SPE statistics are 33% and 36.1%, respectively. Figure 5(b) shows the detection results of PNMF with a detection rate of N 2 statistics is 35.9% and SPE statistics is 34.3%. In Figure 5(c), the detection results of MPNMF with the two detection rate of N 2 statistics and SPE statistics are 36.8% and 98.9%, respectively. Therefore, the detection performance of the proposed MPNMF method is superior to NMF and PNMF, and the SPE detection effect is more obvious, which proves that the method is effective. Fault 10: the fault is caused by the temperature variation of C feed. The faults are introduced from the 161th to the 960th sample points, and the fault type is random fault. The detection results of the NMF, PNMF, and MPNMF methods are shown in Figure 6. Figure 6(a) shows the detection results of NMF with the detection rates of N 2 statistic and SPE statistic are 55.4% and 57.4%, respectively. Figure 6(b) is the PNMF detection result that the detection rates of N 2 statistic and SPE statistic are 57% and 57.8%, respectively, and it is obvious that the detection effect is better than that of NMF. Figure 6(c) is the detection result of MPNMF, and the detection rates of N 2 statistic and SPE statistic are 63.6% and 59.5%, respectively. The overall detection effect of MPNMF is better than NMF and PNMF. In summary, in the monitoring of the TE process, the MPNMF method proposed in this paper has greater advantages than the traditional NMF method and PNMF method.
To compare the computational complexity of the various methods, the dataset was run on a PC with a 2.2 GHz dual-core processor (i5 Intel processor), 4GB RAM, and 1 TB HDD. Table 3 shows the computational time of each fault monitoring program. Computational time is in seconds. Since the fault monitoring process includes offline modeling and online monitoring, Table 3 gives the computational time of offline modeling and online monitoring. Although the proposed method has the longest offline modeling time, the online monitoring time is shorter. Industrial process monitoring is expected to have higher accuracy and lower monitoring time. As can be seen from Tables 2 and 3, the proposed method has better monitoring results.

Conclusion
Traditional process monitoring methods usually adopt a global model of data and ignore local information about the data. Compared with the NMF algorithm, the PNMF algorithm converges faster and has more advantages in data reduction and decomposition. Combining the advantages of block modeling and PNMF algorithm, a multimode process monitoring method based on MPNMF is proposed. In order to make full use of the local information of the process data, the entire process is divided into several subblocks by a complete link algorithm, a PNMF model of each subblock is established, and then a joint probability indicator is constructed to identify the current running mode of the process. Finally, BIC is adapted to synthesize the statistics of each subblock to realize process monitoring. The proposed process monitoring method is applied to the TE process to verify its effectiveness, and compared with the traditional PCA algorithm, NMF algorithm, and PNMF algorithm, the results show that the method can monitor the fault timely and effectively.
In future work, a more flexible model can be used instead of the PNMF model. With the help of historical fault data information, building a more optimized multiblock distributed monitoring model will become the focus of future research.

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 competing interests.