State monitoring and fault prediction of centrifugal compressors based on long short–term memory and principal component analysis (LSTM-PCA)

Construction of monitoring indicators and threshold setting

In this paper, the training sample is 1,500 samples, and the monitoring sample window is a large sample of 131 data, and the two samples are independent of each other and conform to the normal distribution. We use a PCA-based process monitoring method to make judgments using the Hotelling-statistic and the square prediction error (SPE) statistic (Bakdi & Kouadri, 2017). The process is as follows:

Build a sampling matrix Y_t with future and past data at each sampling point.

Y_t by sliding window for 100 past and current moment data fragments {x(t_i)... x(t_i−100)} and 30 predicted snippets of future data {y(t_i)... y(t_i+30)} Compose. ${\displaystyle Y_t=\left(\begin{array}{cccc}\hfill x_1\left(t_{i-100}\right)\hfill & \hfill x_2\left(t_{i-100}\right)\hfill & \hfill \cdots \hfill & \hfill x_{12}\left(t_{i-100}\right)\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill x_1\left(t_{i-1}\right)\hfill & \hfill x_2\left(t_{i-1}\right)\hfill & \hfill \cdots \hfill & \hfill x_{12}\left(t_{i-1}\right)\hfill \\ \hfill x_1\left(t_i\right)\hfill & \hfill x_2\left(t_i\right)\hfill & \hfill \cdots \hfill & \hfill x_{12}\left(t_i\right)\hfill \\ \hfill y_1\left(t_{i+1}\right)\hfill & \hfill y_2\left(t_{i+1}\right)\hfill & \hfill \cdots \hfill & \hfill y_{12}\left(t_{i+1}\right)\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill y_1\left(t_{i+30}\right)\hfill & \hfill y_2\left(t_{i+30}\right)\hfill & \hfill \cdots \hfill & \hfill y_{12}\left(t_{i+30}\right)\hfill \end{array}\right)}$

where x(t_i)…x(t_i−100) are the 100 historical data before the current moment and y(t_i+1)…y(t_i+30) are the 30 future data predicted at the current moment.

Then, the data is filtered by fault-tolerant filtering (Shaolin, Na & Wenming, 2016) to obtain the matrix $Y_t^{\prime }$. Secondly, the data is normalized. The mean variance uses the mean variance of the training data, and the data in each column of the matrix $Y_t^{\prime }$ is normalized to data with a mean value of 0 and a root mean square of 1, and the matrix ${\hat{Y}}_t$ is obtained.

Firstly, according to the offline data, the normal data H is selected, and the H is {H(t₁)…H(t₁₅₀₀) }, and the dimension is the same as that of the monitoring matrix. Then, the matrix is filtered for fault tolerance (Shaolin, Na & Wenming, 2016), and normalized to obtain the matrix $\tilde{H}$ The matrix $\hat{X}$ is obtained by normalizing the data in each column of the matrix $\tilde{H}$ to data with a mean value of 0 and a root mean square of 1.

Calculate the covariance matrices of the filtered and normalized training data and monitoring data respectively: (1)${\displaystyle R_{ta}=\frac{1}{n-1}{\hat{X}}^T\hat{X}}$ (2)${\displaystyle R_{ts}=\frac{1}{h-1}{\hat{Y}}_{\mathrm{t}}^T{\hat{Y}}_t}$

where n is the number of training samples, h is the number of monitoring samples, and the resulting covariance matrix is a feature 12*12-dimensional matrix.

The eigenvalues were sorted from large to small, and the top k features with a sum of more than 85% of the eigenvalues were selected for PCA dimensionality reduction, which was the eigenvalue. (3)${\displaystyle \frac{\sum _{i=1}^k{\lambda }_i}{\sum _{i=1}^m{\lambda }_i}\ge 0.85.}$

Let the first k eigenvalues from large to small form a diagonal matrix, and k corresponding eigenvectors form a matrix. Namely:

The diagonal matrix and the corresponding eigenvectors of the training set form a matrix: (4)${\displaystyle S_{k_1\ast k_1}=diag\left({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_{k_1}\right)P_{m\ast k_1}=\left[p_1,p_2,\dots ,p_{k_1}\right].}$

The diagonal matrix of the monitoring set and the corresponding eigenvectors form the matrix: (5)${\displaystyle {\hat{S}}_{k_2\ast k_2}=diag\left({\hat{\lambda }}_1,{\hat{\lambda }}_2,\dots ,{\hat{\lambda }}_{k_2}\right){\hat{P}}_{m\ast k_2}=\left[{\hat{p}}_1,{\hat{p}}_2,\dots ,{\hat{p}}_{k_2}\right].}$

For $\hat{X},{\hat{Y}}_t$, the number of samples is still n samples, but the number of features becomes K, and after the dimensionality reduction, it is: (6)${\displaystyle {\tilde{X}}_{n\ast k1}={\hat{X}}_{n\ast m}\ast P_{m\ast k_1}}$ (7)${\displaystyle {\tilde{Y}}_{h\ast k_2}={\hat{Y}}_t\ast {\hat{P}}_{m\ast k_2}.}$ The formula for calculating the matrix of $\breve{X},\breve{Y}$ obtained by reconstructing $\hat{X}$ and ${\hat{Y}}_t$ is as follows: (8)${\displaystyle \breve{X}={\tilde{x}}_{n\ast k_1}{P}_{m\ast k_1}^T={\hat{X}}_{n\ast m}{P}_{m\ast k_1}{P}_{m\ast k_1}^T}$ (9)${\displaystyle \breve{Y}={\tilde{Y}}_{h\ast k_2}{\hat{P}}_{m\ast k_2}^T={\hat{Y}}_t{\hat{P}}_{m\ast k_2}{\hat{P}}_{m\ast k_2}^T}$

where n is the number of rows in the training set, h is the number of rows in the monitoring set, m is the number of corresponding key variables, and k is the number of corresponding variables after dimension reduction.

Utilize PCA-based process monitoring methods for fault monitoring and fault trend prediction. There are two traditional statistical indices: the Hotelling-T² index and the square prediction error (SPE) statistical index Q. Establish improved T² statistical monitoring indices for $\breve{Y}$, specifically constructing T² statistics for ${\breve{Y}}_{i,:}$ (i = 1,2...h), with the calculation formula as follows: (10)${\displaystyle T^2={\breve{Y}}_{i,:}\ast {\hat{P}}_{m\ast k_2}\ast {\hat{S}}_{k_2\ast k_2}^{-1}\ast {\hat{P}}_{m\ast k_2}^T\ast {\breve{Y}}_{i,:}^T}$

where I_mm is the unit matrix and m is the dimension of the monitoring data. ${\breve{Y}}_{i,:}$ represents the ith row vector of $\breve{Y}$.

Since $\breve{X,}\breve{Y}$ obeys a normal distribution with a mean of zero, and the sample sampling points are independent of each other, $F_{\alpha }\left(k_1,n-k_1\right)$ obeys the F distribution with the first degree of freedom k1 and the second degree of freedom n-k1, and at a specific confidence degree of 1 − α, usually α is 0.01, and the $F_{\alpha }\left(k_1,n-k_1\right)$ value can be obtained by looking up the table. Then the statistical control limit is determined (Ahmad & Ahmed, 2021), and the calculation formula is as follows: (11)${\displaystyle T_{\alpha }=\frac{k_1\left(n^2-1\right)}{n\left(n-k_1\right)}{F}_{\alpha }\left(k_1,n-k_1\right)}$

where 1 − α is the confidence level, n is the number of samples in the training dataset, and k1 is the number of features selected after PCA.

The SPE can be used to measure the difference between the observed value and the predicted value based on the PCA model, thus identifying data points that may be anomalous. It measures the sum of squares of the differences between observed values and model predictions. Specifically, for each observation, it can be mapped to the principal component space using the PCA model, and the reconstructed predicted value can be obtained from the inverse transformation of the principal component space. Then, the difference between the observed value and the reconstructed predicted value is calculated, and the squares of the difference are added up to obtain the SPE statistic. A higher SPE value indicates a large difference between the observed value and the model prediction and may indicate an anomaly or change. The control limits for the SPE statistic were proposed by Jackson and Mud Holkar in 1979. Establish an improved SPE statistic (Zhou, Parkj & Liu, 2016) for monitoring $\breve{Y}$, and construct the statistical index Q for ${\breve{Y}}_{i,:}$ (i = 1,2...h): (12)${\displaystyle Q={\breve{Y}}_{i,:}\ast \left(I_{m\ast m}-{\hat{P}}_{m\ast k_2}{\hat{P}}_{m\ast k_2}^T\right)\ast {\breve{Y}}_{i,:}^T}$

where ${\breve{Y}}_{i,:}$ is the row vector of row i of $\breve{Y}$, I_m∗m is the identity matrix, and ${\hat{P}}_{m\ast k_2}$ is the matrix ${\hat{P}}_{m\ast k_2}$ of k eigenvectors corresponding to k eigenvalues from large to small from Eqs. (12) and (13).

Since $\breve{X},\breve{Y}$ obeys a normal distribution with zero mean, and the sample sampling points are independent of each other, it has been proved in the literature (Jacksonj & Mudholkarg, 1979) that Eq. (18) is approximately in line with the normal distribution. (13)${\displaystyle {\left(Q/{\theta }_1\right)}^{h_0}\sim N\left[1+\frac{{\theta }_2{h}_0\left(h_0-1\right)}{{\theta }_1^2},\frac{2{\theta }_2{h}_0^2}{{\theta }_1^2}\right]}$ (14)${\displaystyle {\theta }_r=\sum _{j=k1+1}^m{\lambda }_j^r,r=1,2,3}$ (15)${\displaystyle h_0=1-\frac{2{\theta }_1{\theta }_3}{3{\theta }_2^2}.}$

Therefore, the approximate mean of c is a normal distribution of 0, which is calculated as follows: (16)${\displaystyle c=\frac{{\theta }_1\left[{\left(Q/{\theta }_1\right)}^{h_0}-1-{\theta }_2{h}_0\left(h_0-1\right)/{\theta }_1^2\right]}{\sqrt{2{\theta }_2{h}_0^2}}.}$

Then the control limit of Q is: (17)${\displaystyle Q_{\alpha }={\theta }_1{\left[\frac{c_{\alpha }{h}_0\sqrt{2{\theta }_2}}{{\theta }_1}+1+\frac{{\theta }_2{h}_0\left(h_0-1\right)}{{\theta }_1^2}\right]}^{\frac{1}{h_0}}.}$

c_α is the confidence limit of the standard normal distribution, c_α is 0.99 confidence, λ is the eigenvalue of the corresponding covariance matrix R_ta, and h₀ is the adjustment coefficient.

Condition monitoring and fault prediction

(1) According to Eqs. (5) and (6), the T²-statistic and T²-statistic limit T_α of each sampling point were calculated.

(2) According to Eqs. (7) and (12), the SPE statistic Q and the statistical limit Q_α of each sampling point of the monitoring data were calculated, respectively.

(3) Fault determination

If the system is operating normally, the T²-value of the sample should meet T² <T_α, otherwise it is considered to be faulty.

If the system is functioning normally, the Q-value of the sample should meet Q<Q_α, otherwise it is considered to be faulty.

LSTM predictive analytics

First, a training data matrix X_t is constructed at each sampling point, which is {x(t_i)…x(t_i−400)} represents the first 400 steps of the current moment in the sliding window input vector, where t represents the current moment.

Then, the data X_t is filtered by fault-tolerant filtering (Shaolin, Na & Wenming, 2016) to obtain the matrix $X_t^{\prime }$.

We set the telemetry fragment {x_tk,1,…, x_tk,15},(t_k = t₀ + kh, k =1 , 2, 3… , t₀ is the sampling start time, h is the sampling interval) in the sliding time window [t_k − 7, t_k + 7] with a window radius of 7, and sorts it from the lowest to the largest value as {${\overline{x}}_{t_{k,1}}$,…, ${\overline{x}}_{t_{k,15}}$}, and constructs a quartile mean operator. (18)${\displaystyle \begin{array}{c}\hfill Q\left(x_{t_{k,1}},\dots ,x_{t_{k,15}},15\right)\hfill \\ \hfill =\frac{1}{q_3-q_1+1}\sum _{i=q_1}^{q_3}{\overline{x}}_{t_{k,i}}\hfill \end{array}}$

where: q₁ and q₃ denote the lower and upper quartile ordinal numbers of { ${\overline{x}}_{t_{k,1}},\dots ,{\overline{x}}_{t_{k,15}}$ }, respectively. Based on the good fault-tolerant ability of the quartile operator, a fault-tolerant smoothing filtering algorithm for sliding window quartile combinations can be constructed:

Select the appropriate window width radius H for the data slice: (19)${\displaystyle S_k=\left\{x_{t_i},t_i\in \left[\right.t_k-H,t_k+H\left.\right)\right\}}$

Perform quartile filtering:

(20)${\displaystyle x_{t_i}^{\prime }=Q\left(S_k,15\right).}$

The specific forecast calculation steps are as follows:

(1) The combination of input $X_t^{\prime }$ and hidden states h_t−1 is handled by a hyperbolic tangent function tanh: (21)${\displaystyle tanh\left(\mathrm{x}\right)=\left(\mathrm{e}\mathrm{x}-\mathrm{e}-\mathrm{x}\right)/\left(\mathrm{e}\mathrm{x}+\mathrm{e}-\mathrm{x}\right)}$ (22)${\displaystyle g_t=tanh\left(W_g{X}_t^{\prime }+U_g{h}_{t-1}+b_g\right).}$

(2) Forget the door f_t will determine what information should be removed from the storage unit C_t via an element-based sigmoid function (23)${\displaystyle f_t=\sigma \left(w_f{X}_t^{\prime }+u_f{h}_{t-1}+b_f\right).}$

(3) The input gate will determine what information will be stored in the storage cell C_t via the sigmoid function: (24)${\displaystyle i_t=\sigma \left(w_i{X}_t^{\prime }+u_i{h}_{t-1}+b_i\right).}$

(4) Update the information in the memory unit C_t−1 by partially forgetting the information in the previous memory unit, the method is as follows: (25)${\displaystyle C_t=f_t\ast C_{t-1}+i_t\ast g_t}$

where ∗ represents the element-by-element product function of the two vectors.

(5) Update the output hidden state h_t based on the calculated cell state C_t as follows: (26)${\displaystyle O_{\mathrm{t}}=\sigma \left(w_o{X}_t^{\prime }+u_o{h}_{t-1}+b_o\right)}$ (27)${\displaystyle h_t=O_t\ast tanh\left(C_t\right).}$

(6) The output y_t of the current moment is shown below: (28)${\displaystyle y_t=softmax\left(W_y{h}_t+b_y\right).}$

The predicted data y_tcan be mathematically abstracted as {y(t_i), t_i = t₀ + kh₁, k =1 , 2, 3… }, where t₀ and h₁ is the starting time of the measurement process and the sampling interval, respectively.

Through the training process, the network input weight w_g,i,f,o,y, cyclic weight u_g,i,f,o and bias b_g,i,f,o,y are calculated, and the trained network is used to predict the time series data. In this paper, we use the LSTM model for multi-step prediction of eigenvariables. Specifically, we used a 401-step sliding window as the input to the LSTM, and an output with a 30-step sliding window, meaning that we were predicting 30 data points in the future. To achieve this, we employ a recursive single-step prediction approach and use the predicted values to update the network state. This approach allows us to predict the next 0.5 h of data at once, thus meeting the need for real-time monitoring of future samples. In practical application, we make forecasts every 10 min to ensure the accuracy and timeliness of forecasts. This method not only improves the reliability of the monitoring judgment of each sampling point of the monitoring index, but also provides more response time for the site engineer.

LSTM predictive simulation analysis

Model evaluation

Before the fault, when the fault is detected and after the fault is detected, the predicted root mean square error, mean square error, mean absolute error, and mean absolute percentage error are as shown in Tables 1, 2, 3 and 4 below, 1–12 represents the predicted 12-dimensional variables.

Fault monitoring verification based on combining LSTM and PCA algorithm

Analysis of data causes

The motor fault and its related indicators are selected for verification. After the motor fault, the faulty motor cannot maintain a relatively constant speed and power as before the fault occurred, and the motor stator temperature also gradually changes with the motor fault. They can serve as key indicators of motor failure. The motor was in a shutdown state before 10:00 on 2022/10/14, and ran normally after 10:50. An abnormality occurred around 19:30, and returned to normal around 5:30 on 10/15. The temperature of the motor stator remains unchanged when the motor is stopped, rises after the motor operates normally, drops after an abnormality occurs in the motor, and recovers after the motor recovers. The motor speed and power graphs with the motor stator temperature are shown in Figs. 4 and 5, and the motor temperature anomaly graph is shown in Fig. 6. The colors gradually transition from warm (red and yellow) to cool (blue and green).

Comparison of improved methods

The monitoring effect before the algorithm improvement is shown in Fig. 8. Prior to the improvement, each point only contained the monitoring effect of the current point, leading to false detections, false alarms, and insufficient information content.

After the improvement, each sampling point now includes 100 historical data points prior to the sample and 30 future data points following it, greatly enhancing the reliability and accuracy of the monitoring alarms. If the 10 data points before the monitoring moment are all abnormal, an anomaly is considered to have occurred, and the forecast results for the time after the current moment are provided. The normal monitoring results for each sampling point are shown in Fig. 9, thereby allowing the previously mentioned false detections to be monitored correctly.

Table 5:

Motor faults and solutions.

Cause of failure	Solutions
There is something wrong with the pump	The main oil pump does not start: there is no working medium, notify the responsible department. The auxiliary oil pump does not start when the oil pressure drops: electrical failure, electrical failure of the pump automatic equipment (notify the responsible department)
Oil pipeline leak	Leakage at flange connection: seal must be replaced. Oil line rupture hazard: Fire hazard due to contact with hot parts
Cooler, filter or strainer is dirty	Conversion cooler, filter cleaning coarse filter
Defective oil pressure balancing valve or pressure reducing valve	Check valve and replace if necessary
Oil pressure too low	The pump is defective, the oil pipeline is leaking, the cooler, filter or strainer is dirty, the oil pressure balance valve or pressure reducing valve is defective
The oil level in the high tank is too low	Add enough oil
Seal oil pressure difference is too low	Adjust control valve
Oil temperature is too low	turn on heater
Fault	Motor failure

DOI: 10.7717/peerjcs.2433/table-5