Nonstationary Process Monitoring Based on Cointegration Theory and Multiple Order Moments

: In industrial processes, process data often exhibit complex characteristics, such as nonstationarity and nonlinearity, which brings challenges to process monitoring. In this study, a monitoring strategy for nonstationary processes is proposed based on cointegration theory and multiple order moments. Considering the nonstationarity presented in some variables, cointegration analysis (CA) is applied to obtain long-term equilibrium relationships among these nonstationary variables, which are then combined with stationary variables to form a new stationary dataset. For the purpose of process monitoring, a new monitoring index that contains multiple order moments is proposed to capture different statistical features of a previously obtained stationary data set. Moving windows are applied to capture changes of local statistical characteristics to implement online monitoring. Case studies on simulation data and an industrial dataset are presented to illustrate the effectiveness of the proposed method for nonstationary process monitoring. Comparing with the PCA and common CA-based monitoring methods, the proposed method has better performance with a lower false alarm rate and earlier alarm time.


Introduction
With increasing attention to industrial safety, process monitoring and control have attracted broad public attention. In modern industrial processes, huge amounts of process data have been accumulated due to the development of distributed control system and computer technology. However, such massive data may present complex characteristics such as nonstationarity and nonlinearity. Thus, many monitoring methods have been proposed to deal with this issue.
Due to frequent load adjustments and external disturbances, process data may present clear nonstationary characteristic in industrial processes [1]. Thus, fault information can be easily covered by nonstationary changes in variables, which makes it difficult to properly extract fault features and accurately identify current operating condition. A common solution is to update the model recursively to adapt to normal nonstationary behavior [2][3][4][5]. However, the model will be wrongly adapted if a slow-varying fault is regarded as normal variation.
In addition, differential data preprocessing is used to convert nonstationary data to stationary data, which will lose the dynamic information of process. The cointegration analysis (CA) theory proposed by Granger is an alternative to analyze long-term equilibrium relationships among original multivariate nonstationary time series, and this theory exhibits superiority since no differential data pretreatment is needed [6]. According to the cointegration theory, time-invariant features can be extracted from nonstationary variables as a monitoring index if nonstationary variables show common trends, which indicates that there is at least one stable long-term dynamic equilibrium relationship among these variables.
In recent years, cointegration theory has been introduced to industrial nonstationary process monitoring and has obtained better performance. Li et al. [7] applied CA to separate multivariate nonstationary series into stochastic common trends and equilibrium error, which are used to construct monitoring indices to realize the monitoring of process operation status. However, the inherent serial correlation of process variables is ignored in the above methods. Lin et al. [8] introduced a nonstationary monitoring strategy that removed serial correlations by applying autoregressive filters and a compensation scheme after extracting the common stationary and nonstationary factors by CA.
Most methods are only focused on nonstationary series in process; however, even in a nonstationary process, there are variables that exhibit stationary characteristics that should also be considered in modeling. Zhao et al. [9] proposed a full-condition monitoring method based on cointegration and slow feature analysis where stationary variables and nonstationary variables were separated first. CA was applied to nonstationary variables to extract the long-term equilibrium relations among variables. Then, slow feature analysis was applied to the combination of residuals after CA and stationary variables. Six monitoring indicators were established to monitor static variation and dynamic variation, respectively, by which the changes of operation conditions and real fault can be distinguished.
However, the cointegration relationship may be broken under certain circumstances. Yu et al. [10] proposed a recursive cointegration analysis (RCA) for nonstationary industrial processes, in which the monitoring model was updated when the long-term equilibrium relationship of process variables extracted by cointegration analysis changed. Thus, the method can avoid frequent model updating compared with common adaptive methods. On this basis, Zhang et al. [11] applied recursive principal component analysis (RPCA) to capture remaining short-term dynamic information after extracting the long-term equilibrium relationship by RCA. Thus, a comprehensive monitoring framework was built.
The above studies show that cointegration analysis is an effective tool in monitoring nonstationary industrial processes. However, only the linear relationship among variables is concerned in cointegration analysis, while the nonlinear relationship among variables is neglected. In industrial processes, the relationship among variables is more complicated than a linear relationship. Thus, dynamic relationship among system variables cannot be described comprehensively by linear cointegration analysis. More precisely, higherorder dynamic relationship of variables will be remained in model residual. Thus, the residual of long-term equilibrium relationship of multivariate nonstationary time series will contain dynamic components and exhibit non-Gaussian traits, which will lead to incorrect monitoring results.
In practical applications, it is impossible to obtain all data, which also leads to a huge workload. In this way, the statistics of part of the samples are widely applied to estimate the overall features. Thus, statistics are often used as monitoring indices to estimate the condition of industrial processes. In process monitoring, the T 2 statistic and SPE statistic are commonly used statistics for principal component analysis [12].
Statistics are mainly divided into two types: the measure of central tendency and measure of variation. The former is an index used to express the central tendency among variables, such as the mean; while the latter is an index to express the degree of dispersion within variables, such as the standard deviation and variance. The statistics are usually first order and second-order statistics, which are assumed to represent the characteristics of data that obey normal distribution. However, when the data disobey a normal distribution, higher-order statistics are required to better describe the characteristics of the data.
For non-Gaussian or non-linear signals, higher-order statistic techniques are effective tools to capture statistical characteristics of the signals, which have been widely used in the field of signals analysis [13]. Choudhury et al. [14] applied higher order statistics to analyze the control loop performance. Two new higher order statistics indices were proposed to detect and quantify the potential non-Gaussian and nonlinearity of signals in process.
Kara et al. [15] introduced a statistic moment approach to analyze structures with variable parameters. Higher order statistics measure the extent of nonlinearity and quantify the non-Gaussian of the probability distribution of the process variables. However, the signal processing methods only focus on one dimensional data, which are not appropriate to be applied in multivariate process monitoring.
Thus, Wang et al. [16] proposed a novel approach based on higher order cumulants to monitor a multivariate process. The method effectively extracts both amplitude and phase information by applying higher order cumulants, especially for non-Gaussian distribution. Wang et al. [17] proposed a method that replaced original variables with the statistics of each variable when applying PCA to realize process monitoring. By applying different order statistics to characterize process behavior, the results showed better performance in nonlinear process.
Considering the complex characteristics presented in industrial process data, a monitoring strategy based on cointegration analysis and comprehensive statistics is proposed. For nonstationary characteristics of some variables, cointegration analysis is first applied to these nonstationary series to extract time-invariant characteristics, which are then combined with stationary variables as a new data set. A new monitoring index containing multiple order moments, such as the mean, variance and skewness, is employed to capture different features of the generated data set, which makes it sensitive to fault information in process monitoring.
The rest of this paper is organized as follows. Section 2 introduces the basic theory of stationary test, cointegration analysis and multiple order moments. Section 3 presents the monitoring strategy based on proposed method. In Section 4, the proposed method is applied to the numerical data and industrial data. The final section gives our conclusions.

Stationary Test
The stationarity of time series can be determined by the method of statistical testing. A common stationarity test is the Augmented Dickey-Fuller (ADF) unit root test. For a p-order autoregressive AR (p) process with no deterministic trend, y t = a 1 y t−1 + a 2 y t−2 + · · · + a p y t−p + u t (1) where p is the lag order and u t is the stochastic part of the series, u t ∼ N 0, σ 2 . The corresponding vector error correction model is obtained by parameter transformation and adjustment to Equation (1), which can be written as follows: The stability of the sequence is determined by the hypothesis test of the coefficient ρ, which is estimated by ordinary least square method. The AR (p) process has a unit root if ρ = 1, which leads to nonstationarity of the process.
If a nonstationary sequence can be transformed into a stationary variable by multiple first-order differences, the sequence is called an integral sequence, and the time of differences is called an integral order. The stationary time series is called "zero-order integration" and denoted as I (0). If the first-order difference of time series is a stationary process, it is called "first-order integration", denoting I (1) and is known as the "unit root process".

Cointegration Theory
The study on the cointegration relationship is mainly focused on the I (1) process, in which the cointegration method is used to find the common random trend among variables. The basic idea is that, although some variables exhibit nonstationary characteristics individually, they may have a common stochastic trend, which means that there is a stable long-term equilibrium relationship among them. Thus, one or more stationary sequences are obtained by specific linear combinations. For example, it is assumed that sequences in Y t = (y 1t , y 2t , · · · , y mt ) T are all integrated series of first order where m is the number of variables.
For a p-order vector auto-regression model (VAR): where p is the lag order and ε t is a p-dimensional joint normal process, ε t ∼ I N(0, Ω ε ). The corresponding vector error correction model is obtained by parameter transformation and adjustment on Equation (3), which can be written as follows: The maximum Likelihood Estimation (MLE) is applied to estimate the parameters [18] in Equation (4), in which matrix Φ contains the main characteristics and information of the cointegration system. The rank of matrix Φ is assumed to be r. If r = 0, Equation (4) represents a common VAR model of first order difference, which means no cointegration relationships existing among variables. Thus, the rank of matrix Φ equal to the number of independent cointegration relationships, which can be obtained by testing the significance of characteristic roots of matrix Φ [6]. A commonly used testing method is the eigenvalue trace test and the test statistic λ trace is calculated as follows, The null hypothesis of the eigenvalue trace test is that the number of independent cointegration vectors is less than or equal to r. Clearly, λ trace (r) = 0 if the null hypothesis is true. The critical value of λ trace can be obtained by Monte Carlo simulation. Thus, the number of independent cointegration vectors is obtained by such a statistical test. If the number of cointegration vectors is r, the first r lines ofβ estimated from the maximum likelihood estimation is the cointegration coefficient matrix. If 0 < r < m, matrix Φ can be decomposed as: where α, β R m×r , and both α and β are full rank matrices. Thus, the "stationary part" ξ t can be calculated by: in which ξ t is a stationary series of I(0), and β r is called the cointegration coefficient matrix. The cointegration relationships among variables exist when r > 0, which indicates that there are r long-term equilibrium relationships among these nonstationary sequences.

Multiple Order Moments
Moments are statistics that can be employed to describe the probability distribution of variables.
The first order moment is the mean of a set of samples: The second order central moment is the variance: Processes 2022, 10, 169 5 of 12 The first-order moment corresponds to the mean and the theorem of large numbers, and the second-order central moment corresponds to the variance and the central limit theorem. The first two moments correspond to the two most important theorems of statistics, which determine the most important properties of the probability distribution. The weak stationarity in time series analysis is to see whether the first and second moments are stable.
The third order standard moment is called skewness: The fourth order standard moment is kurtosis: The third-order moments and above are called standard moments, and, just as the influence of variance should be removed from the mean, the influence of variance should also be removed from skewness and kurtosis. Thus, these moments are assumed to be "independent" or "orthogonal" to each other. Moments commonly used in general statistics are up to the fourth order.
As the order of moments increases, each moment provides more detailed probability distribution information. Together with lower-order moments, a more complete description of the probability distribution (from the mean, variance, skewness, kurtosis . . . ) can be obtained with the similar idea for Taylor series and Fourier series. However, the higher the order of the moment, the more data is needed for estimation and the more difficult it is to describe and understand the moment.

The Procedure of the Proposed Monitoring Framework
In the nonstationary process, some variables exhibit obvious nonstationary characteristics, while other variables remain stationary. Thus, for nonstationary process monitoring, the variables are first separated into stationary and nonstationary variables by the ADF stationary test. Cointegration analysis is then applied to nonstationary variables to extract stationary features among variables, which are combined with stationary variables as a new dataset.
The new dataset is divided into several windows to capture changes in the local characteristics for the purpose of online detection, and the window width is selected based on trials with different window widths. If the window width is too small, the statistic characteristics of variables cannot be captured properly, which will be greatly influenced by the fluctuation of the data, while, if the window width is too large, the monitoring model is insensitive, and thus the fault cannot be detected in a timely fashion.
The proper window width l is selected within a certain range with the balance of model sensitivity and well feature extraction. Then mean matrix M, variance matrix V and skewness matrix S are obtained, which are combined into A = [M, V, S] as a set of multiple order moments. Based on the calculation of A, monitoring index m is calculated by where a is one of the row vectors in A. Thus, different features reflected from each order of statistics can be expressed by the combined monitoring index, which is sensitive to the abnormal process. The control limit of the monitoring index is obtained by kernel density estimation [19] of the training data. The Gaussian kernel function is selected as the kernel function. The diagram of nonstationary process monitoring based on cointegration theory and multiple order moments is shown in Figure 1.
where is one of the row vectors in . Thus, different features reflected from each order of statistics can be expressed by the combined monitoring index, which is sensitive to the abnormal process. The control limit of the monitoring index is obtained by kernel density estimation [19] of the training data. The Gaussian kernel function is selected as the kernel function.
The diagram of nonstationary process monitoring based on cointegration theory and multiple order moments is shown in Figure 1.

Offline Modelling
Data under normal conditions is selected as training data with samples and variables, which is normalized and denoted as × . The ADF unit root test is applied as a stationary test method to separate into stationary data set × 1 and nonstationary dataset × 2 , in which 1 and 2 are the number of nonstationary variables and the number of stationary variables, respectively.
The Johansen cointegration test is applied to nonstationary dataset , in which the lag order of the cointegration analysis model is determined by the AIC criterion. The number of cointegration vectors is obtained through the trace eigenvalue test, and the cointegration coefficient matrix = ( 1 , 1 , ⋯ , ) is obtained from the maximum likelihood estimation, the stationary part = ( 1 , 2 , ⋯ , ) extracted from nonstationary dataset is obtained with Equation (7).
The stationary part of the nonstationary variables are combined with the stationary variables = ( 1 , 2 , ⋯ , 2 ) to construct a new dataset denoted as = ( 1 , 2 , ⋯ , , 1 , 2 , ⋯ , 2 ), ×( + 2 ) . A window with width of is applied to , and each window is a × ( + 2 ) matrix denoted as , = 0,1, ⋯ , − . Calculate multiple order moments of each window matrix . Take, for example, the mean of each window is stacked vertically to a matrix denoted by ( − )×( + 2 ) with each row representing the mean of , = 0,1, ⋯ , − . The variance matrix and skewness matrix are obtained the same as the mean matrix . A set of multiple order moments is constructed by stacking the statistical moment matrix , , horizontally, which is denoted as = Figure 1. The diagram of the proposed monitoring method.

Offline Modelling
Data under normal conditions is selected as training data with n samples and m variables, which is normalized and denoted as X t R n×m . The ADF unit root test is applied as a stationary test method to separate X t into stationary data set X s t R n×m 1 and nonstationary dataset X n t R n×m 2 , in which m 1 and m 2 are the number of nonstationary variables and the number of stationary variables, respectively.
The Johansen cointegration test is applied to nonstationary dataset X n t , in which the lag order p of the cointegration analysis model is determined by the AIC criterion. The number of cointegration vectors r is obtained through the trace eigenvalue test, and the cointegration coefficient matrix B = (β 1 , β 1 , · · · , β r ) is obtained from the maximum likelihood estimation, the stationary part ξ t = (ξ 1t , ξ 2t , · · · , ξ rt ) extracted from nonstationary dataset is obtained with Equation (7).
The stationary part ξ t of the nonstationary variables are combined with the stationary variables X s t = (x 1 , x 2 , · · · , x m 2 ) to construct a new dataset denoted as Y t = (ξ 1t , ξ 2t , · · · , ξ rt , x 1 , x 2 , · · · , x m 2 ), Y t R n×(r+m 2 ) . A window with width of l is applied to Y t , and each window is a l × (r + m 2 ) matrix denoted as y it , i = 0, 1, · · · , n − l. Calculate multiple order moments of each window matrix y it . Take, for example, the mean of each window is stacked vertically to a matrix denoted by M R (n−l)×(r+m 2 ) with each row representing the mean of y it , i = 0, 1, · · · , n − l.
The variance matrix V and skewness matrix S are obtained the same as the mean matrix M. A set of multiple order moments is constructed by stacking the statistical moment matrix M, V, S horizontally, which is denoted as A = (a 1 , a 2 , · · · , a n−l ) T , A R (n−l)×3(r+m 2 ) . The monitoring index is calculated by m i = a T i a i , i = 1, 2, · · · , n − l, and the kernel density estimation is applied to estimate the control limit C α with a significance level of α.

Online Monitoring
We scaled the testing data with the mean and variance obtained from the training data and separated the scaled testing data into a stationary dataset and nonstationary dataset denoted as X s new and X n new . We projected the nonstationary dataset X n new to the cointegration coefficient matrix B to obtain the stationary part ξ new = B T X n new of the new window in each Processes 2022, 10, 169 7 of 12 cointegration direction, which was then combined with the stationary variables X s new to construct a new dataset denoted as Y new . We calculated the monitoring index m new and compared it with the control limit C α to determine the current operating conditions.

Numerical Example
To verify the effectiveness of the proposed method, a multivariate nonstationary process Y t = (y 1t , y 2t , y 3t , y 4t ) containing four variables was constructed as follows. Each sequence contains the same random walk process x t as a simulation of common stochastic trend. The common stochastic trend and other four sequences were constructed as follows: In Equations (13)- (17), e 1t ∼ N(0, 1), e 2t , e 3t , e 4t ∼ N(0, 0.5), the number of sample points is 2000. A soft fault is introduced to {y 2t } at 1000 sample points: To compare with the proposed method, a principal component analysis (PCA)-based monitoring method and common CA-based method [7] were also adopted.
The first two subplots in Figure 2 depict the monitoring results based on PCA, in which the blue solid line represents the statistics, and the red solid line represents the control limit. The T 2 statistic and Q statistics are the monitoring indices, which are computed by Equations (19) and (20), respectively.
where t represents the score vector of PCA and Λ k is the covariance of t. The third subplot in Figure 2 shows the monitoring result based on the common CAbased method, and the monitoring index is the T 2 statistic that is computed by Equation (21) in which B is the cointegration matrix, x t is the original multiple series, and Λ ξ is the covariance of B T x t . The fourth subplot in Figure 2 depicts the result based on the proposed method, in which the blue solid line represents the monitoring index that is calculated by Equation (12), and the red horizontal line is the control limit. A validation set is used to determine the proper window width, which is set to be 15, 30, 60, 90 and 120, respectively. The monitoring results of the alarm time and false alarm rate are shown in Table 1, and 60 is considered to be the suitable window width with weighing the alarm time and false alarm rate. (12), and the red horizontal line is the control limit. A validation set is used to determine the proper window width, which is set to be 15, 30, 60, 90 and 120, respectively. The monitoring results of the alarm time and false alarm rate are shown in Table 1, and 60 is considered to be the suitable window width with weighing the alarm time and false alarm rate.  As can be seen in Figure 2, for the multivariate process with nonlinearity and nonstationarity, the monitoring results based on the proposed method exhibit better performance compared with the monitoring methods based on PCA and the common CA method. The monitoring result based on the PCA method shows that the fault occurs after the 1500th sample both in the principal part and residual part, which is at least 500 samples too late for an alarm.
The method based on the common CA method only explores the linear long-term equilibrium among nonstationary series, thus, when a nonlinear relationship exists among variables, the false acceptance rate is high as is shown in the monitoring results in the third subplot in Figure 2. However, the monitoring methods based on the proposed monitoring strategy exhibit great properties to identify the normal conditions and real fault with alarms in time and no false alarms.  As can be seen in Figure 2, for the multivariate process with nonlinearity and nonstationarity, the monitoring results based on the proposed method exhibit better performance compared with the monitoring methods based on PCA and the common CA method. The monitoring result based on the PCA method shows that the fault occurs after the 1500th sample both in the principal part and residual part, which is at least 500 samples too late for an alarm.
The method based on the common CA method only explores the linear long-term equilibrium among nonstationary series, thus, when a nonlinear relationship exists among variables, the false acceptance rate is high as is shown in the monitoring results in the third subplot in Figure 2. However, the monitoring methods based on the proposed monitoring strategy exhibit great properties to identify the normal conditions and real fault with alarms in time and no false alarms.

Industrial Case
In this section, the proposed method is applied to monitor the industrial data collected from a continuous catalytic reformer (CCR) unit, and the monitoring results are compared with the monitoring methods based on PCA and the common CA method respectively.
The main monitoring purpose is to identify the abnormal increase of hot end pressure drop of the reforming feed heat exchanger in a catalytic reforming unit. In this work, 3240 samples were used to validate the effectiveness of the proposed mothed. The sampling frequency was 1 min, and each sample was composed of 26 process variables. The first 1800 samples were used as a training dataset to establish the monitoring model, and the remaining samples were used as test datasets that contained abnormal working conditions. The change of the pressure drop mentioned above is mainly affected by the changes in the feed volume of naphtha and circulating hydrogen. As shown in Figure 3, the subplots from top to bottom are the variation of the pressure drop, the feed volume of naphtha and the feed volume of circulating hydrogen over time in the sampling interval. The pressure drop of the heat exchanger shows an upward trend after 2700 sample points, in which the feed volume has not increased while the pressure drop has a clear increasing trend; therefore, this is considered as an abnormal rise in the pressure drop.

spectively.
The main monitoring purpose is to identify the abnormal increase of hot end pressure drop of the reforming feed heat exchanger in a catalytic reforming unit. In this work, 3240 samples were used to validate the effectiveness of the proposed mothed. The sampling frequency was 1 min, and each sample was composed of 26 process variables. The first 1800 samples were used as a training dataset to establish the monitoring model, and the remaining samples were used as test datasets that contained abnormal working conditions.
The change of the pressure drop mentioned above is mainly affected by the changes in the feed volume of naphtha and circulating hydrogen. As shown in Figure 3, the subplots from top to bottom are the variation of the pressure drop, the feed volume of naphtha and the feed volume of circulating hydrogen over time in the sampling interval. The pressure drop of the heat exchanger shows an upward trend after 2700 sample points, in which the feed volume has not increased while the pressure drop has a clear increasing trend; therefore, this is considered as an abnormal rise in the pressure drop. The ADF unit root test is first applied to the 26 process variables in the training dataset to separate the original dataset into nonstationary variables and stationary variables for modeling. The test results are shown in Table 2.  The ADF unit root test is first applied to the 26 process variables in the training dataset to separate the original dataset into nonstationary variables and stationary variables for modeling. The test results are shown in Table 2. If the p-value is less than 0.05, the variable is stable at the 5% significance level, and the test results show that 4 variables in bold are nonstationary as shown in Table 2. Then the first-order difference of these variables is tested for stationarity. The test results show that the first-order differences of these variables are stable, and the results are not shown here.
The four first-order single integer variables were used as nonstationary modeling variables. The Johansen cointegration test was performed on these variables, and the AIC criterion was used to determine the lag order of the VAR model. The lag order was determined to be 2 after calculation.
The number of cointegration relations was tested by trace statistics, and the test results are shown in Table 3. The results in Table 3, at the 5% significance level, indicate rejecting the hypothesis of r ≤ 1 and accepting the hypothesis of r ≤ 2. Thus, the number of cointegration relations can be described as r = 2-that is, there are two forms of linear combinations to eliminate the common stochastic trend among these nonstationary variables.
The testing data and cointegration coefficient matrix obtained from the training data were substituted into Equation (7). A labeled dataset containing abnormal conditions was applied as the validation set to obtain the optimal window width in the same way as described in Section 4.1. The selected window width was 80 in this case. The monitoring index and control limits were calculated according to the method mentioned in Section 3. The monitoring results based on PCA, the common CA method and the proposed method are shown in Figure 4, in which the blue line represents the statistics, and the red horizontal line represents the control limits. If the blue line exceeds the red line, it is considered that the current condition is abnormal.  As shown in Figure 4, after about 2800 samples, the statistics of all three monitoring methods exceed the control limits, which indicates that an abnormal rise in the pressure drop after 2800 samples was identified by all these methods. However, in this process, several variables exhibited obvious nonstationary characteristics, while the monitoring method based on PCA assumed that the variables in the process were all stable. Thus, the nonstationary characteristics of certain variables led to a high false alarm rate when PCA was applied.
As shown in the first subplot in Figure 4, the alarm time of the 2 statistic was later than the cointegration analysis-based methods. While the Q statistic exceeded the control limit at a normal condition as shown in the second subplot in Figure 4, it was ineffective to monitor nonstationary process with a high false alarm rate.
Cointegration analysis is an effective tool to deal with multivariate nonstationary se- As shown in Figure 4, after about 2800 samples, the statistics of all three monitoring methods exceed the control limits, which indicates that an abnormal rise in the pressure drop after 2800 samples was identified by all these methods. However, in this process, several variables exhibited obvious nonstationary characteristics, while the monitoring method based on PCA assumed that the variables in the process were all stable. Thus, the nonstationary characteristics of certain variables led to a high false alarm rate when PCA was applied.
As shown in the first subplot in Figure 4, the alarm time of the T 2 statistic was later than the cointegration analysis-based methods. While the Q statistic exceeded the control limit at a normal condition as shown in the second subplot in Figure 4, it was ineffective to monitor nonstationary process with a high false alarm rate.
Cointegration analysis is an effective tool to deal with multivariate nonstationary series by eliminating the common trend among variables to extract time-invariant features, which are then used as an index to monitor the operation conditions. The commonly used statistic is T 2 , which is used to monitor the principal part of the extracted time-invariant features among nonstationary variables. However, the methods only focus on the nonstationary variables and their linear long-term relationship among them, which is not sufficient for complex industrial processes. The monitoring results in the third subplot in Figure 4 show that the statistics of several samples exceeded the control limit before fault occurred, which led to false alarms.
The monitoring results of the proposed method as shown in the fourth subplot in Figure 4 exhibited superiority with no false alarms, and the statistics exceeded the control limit in time when the pressure drop rose abnormally. With stationary variables combined in the model, the information included in the monitoring index is more comprehensive. In addition, multiple order moments were applied to represent different features of the process, which is appropriate for a monitoring index. Thus, the proposed method exhibits its effectiveness in nonstationary process monitoring.

Conclusions
For nonstationary processes that cannot be monitored by traditional multivariate statistical methods, a monitoring strategy based on cointegration theory and multiple order moments is proposed. Unlike most cointegration-based methods that only focus on nonstationary variables, stationary variables are also included in the monitoring model in this work, which ensures that the model contains more comprehensive information of the process. In addition, multiple order moments of the stationary part extracted by cointegration analysis, such as the mean, variance and skewness, were used to construct a new monitoring index that contains different features represented by each order moment.
Moving windows were applied to capture changes in the local statistical characteristics for the purpose of online monitoring. Thus, the proposed monitoring method is sensitive to the abnormal conditions and can realize early identification of faults with slow features. The case study on simulation data indicates that the proposed method is effective to monitor nonstationary processes with low false alarm rates even when nonlinear relationships exist among variables.
For industrial data, compared with PCA and the common cointegration analysis monitoring method with the T 2 statistic, the monitoring results of the proposed method demonstrated the superiority in the fault alarm rate and alarm time. Therefore, for the monitoring of nonstationary processes, the use of multiple order moments as a monitoring index based on cointegration analysis can provide early alarms for abnormal conditions and can effectively identify normal changes and abnormalities.