Research on non-parametric prediction method of reciprocating compressor time series based on prediction credibility scale

Abstract

Aiming at the long-term unpredictability of the reciprocating compressor vibration signal, a non-parametric prediction method of reciprocating compressor time series based on the prediction credibility scale is proposed in this paper. The method is to take the multifractal singular spectrum as the prediction parameter and use the Smoothness Priors Approach (SPA) method to obtain the singular spectrum parameters of different components, and construct the phase space reconstruction dynamic modeling domains. It enables the prediction model to reflect the real-time characteristics of the dynamics evolution of complex systems and highlights the independent influence of each component on the prediction. Meanwhile, the information entropy saturation principle is introduced into the K-Nearest Neighbor (KNN) model to establish the improved K neighborhood dynamic non-parametric prediction model based on the maximum prediction credibility scale, which improves the credibility of the prediction results. Finally, a complete SPA&PSR_KNN prediction algorithm is proposed. Through example validation and error analysis, compared with KNN, BP, and SVM, it can be seen that the prediction results of spectral characteristic parameters obtained by this algorithm have smaller error and higher reliability, and faster operation speed. Thus, the prediction of vibration signal time series of reciprocating compressor is realized.

Introduction

Reciprocating compressor is a typical nonlinear system, and its vibration signal exhibits the characteristics of periodic oscillation and chaotic fractal, that is, long-term unpredictability. ¹ For this kind of signal prediction, it is obviously undesirable to pursue the error between the vibration predicted amplitude and the measured value at a certain time, while ignoring the maximum prediction credibility scale. Meanwhile, the prediction credibility scale is time-varying restricted by the initial sensitivity of the system. Therefore, the prediction results based on fixed-parameter models, which are limited to a few finite discrete point values at a certain moment in time, have no practical significance because they cannot show the refined characteristics of the complete signal cycle.

Currently, the parametric-based prediction methods include time series modeling, Kalman filtering, regression analysis, and neural network analysis.^2–9 The prediction based on time series can be divided into regression analysis models such as auto regressive (AR), moving average (MA), and auto regressive integrate moving average (ARMA). This type of prediction model usually needs to perform steps such as sequence model identification, model order determination, parameter estimation, and model verification. Although the model technology is relatively mature, it is more suitable for stationary time series prediction because its parameters cannot be transplanted and its initialization is complicated. Kalman filter method has a large amount of parameter estimation and high precision for processing stationary data, which is not suitable for state prediction of nonlinear and stochastic systems. The regression analysis model takes the factors determining the predicted objects as independent variables and the predicted objects as dependent variables to establish a mathematical regression equation of the relationship between independent variables and dependent variables. According to the different parameters of the regression equation, it can be divided into univariate or multivariate linear prediction, nonlinear regression prediction model and so on. Its prediction model is simple structure and high calculating speed, but for complex nonlinear system and the insufficient external accidental factors to consider, there is a big limitation. ¹⁰ The methods based on Back Propagation (BP), fuzzy or high-order neural network can optimize the model through a large amount of training, but the complexity of the training time and function mapping make it impossible to show the advantages of real-time predictions in the complex nonlinear system model construction mechanism. ¹¹ Comprehensive above the analysis of the prediction method based on the model parameters, it can be clarified that the initial sensitivity of the complex nonlinear mechanical system with chaotic characteristics and the sudden change of the characteristic state caused by the “singular point” in the operating state, make the Eager type algorithm based on the model and parameters is difficult to play its advantages of training and learning.

For the prediction of dynamic nonlinear system, the main methods include artificial neural network prediction, fuzzy prediction, support vector machine, and their fusion algorithms. Most of the above prediction model take the single-step prediction mode and the rolling prediction mode with fixed constants as the research object. Due to the limitation of the existing model order or parameters, the adjustment and training of the model are time-consuming and laborious, and as a result, the adaptability is constrained by the uncertainty factors within the system state, and the model cannot be distorted according to the change of system state characteristics.

In summary, the non-parametric regression-based prediction method is a modeling method suitable for uncertain, nonlinear dynamic systems that has emerged in recent years. ⁸ The core of the algorithm is not to establish an approximate model for historical data and train, but to obtain information directly from historical data, and to deduce future results according to the relationship of intrinsic factors contained in historical data. The principle is internally consistent with the prediction system based on multivariate metric and wave theory proposed by Weng. ¹² The former does not require prior knowledge, and the core of the principle is to use historical data to find similar nearest neighbors to the current point to predict the value of the next moment, which effectively avoids the model error brought by smoothing. The latter proposes “fidelity algorithm” and “frequency division technique,” that is, for the complex periodic data, the technique of frequency separation and classification of data of different depths is used to improve the accuracy of prediction by comparison and auxiliary judgment.

Parzen window estimation and K-Nearest Neighbor (KNN) estimation are mainly used for non-parametric prediction pattern recognition, which is more applicable than Bayesian estimation that requires prior determination of prior probability and probability density assumption. The Parzen window estimates the probability density from the sample, while the KNN algorithm directly omits the estimation of the unknown probability density function. Since the selection and trade-off of the window function is an obstacle to the application of the Parzen algorithm, while the KNN algorithm is based on the posterior probability of the sequence element distance or similarity measure, which is more generalizable.^13–15

The principle of the prediction algorithm based on KNN estimation is to construct a certain segment vector sequence, find K neighbors by the method of distance or similarity measurement, and obtain a new prediction value through a certain average or weighting process. Due to the simple and intuitive algorithm of non-parametric prediction methods, when the sample value is large enough, its error is low, and the prediction effect under special (fault or singular point) conditions is more accurate than parametric modeling. This type of method has been widely used in the fields of transportation, molecular biology, finance, and meteorology for prediction. Scholars have carried out related research based on this method. Zhu Hong combined the KNN algorithm with the Bayesian algorithm to improve the time efficiency and accuracy of sample classification. ¹⁶ Feng and Meng ¹⁷ predicted flight departure time based on KNN and Support Vector Regression (SVR). Compared with the Support Vector Regression model alone, the relative absolute error was at a lower level, and the prediction accuracy was significantly improved. In view of the low prediction accuracy of non-parametric models for complex systems, Lin ¹⁸ improved the typical KNN method and combined the Empirical Mode Decomposition (EMD) method to predict the daily closing prices of different stocks, which effectively improved the prediction accuracy.

At present, nonparametric regression-based prediction algorithms have achieved certain results in many fields, but there are many complications in the field of fault prediction that restrict their further development. For example, the modeling process completely ignores the internal evolutionary state of the system, the requirements of the predicted parameters and prediction time for the model, the selection of the number of neighboring values and the mechanism of model construction, and so on. How to improve the model interpretability and prediction accuracy of the algorithm, and make it suitable for the prediction and evaluation of reciprocating machinery fault diagnosis is an urgent problem to be solved.

The KNN algorithm is applicable to predict the feasibility study of complex systems, but also requires a certain optimization means, to enhance the adaptability and interpretability of the algorithm, namely, the need for nonlinear time-varying systems based on the analysis of the dynamic evolution trend, with reciprocating machinery vibration signal to build effective prediction of time series sequence method. Based on the G_P algorithm proposed by Grassberger and Procaccia, fractal scale relationship can be found from the perspective of phase space reconstruction, and the maximum prediction credibility scale can be calculated by combining the principle of information entropy saturation. It can be seen that the information entropy theory can be used to analyze the evolution process of nonlinear systems and improve the prediction reliability of nonlinear complex systems. Combined with the theory of in-phase space reconstruction, the prediction model of time series complex systems can be constructed.

Aiming at the adaptability of the reciprocating compressor vibration signal prediction model and the initial sensitivity of the nonlinear system, a non-parametric prediction method of reciprocating compressor time series based on the prediction credibility scale is proposed in this paper. Using the multifractal singular spectrum as the prediction parameter, an improved K neighborhood dynamic prediction model based on the maximum prediction credibility scale is proposed. The principle of information entropy saturation is introduced to the maximum prediction credibility scale, which improves the credibility of the prediction results. The phase space reconstruction type dynamic modeling domain is constructed with the trend and de-trend component spectrum parameters decomposed by the SPA, and a complete SPA&PSR_KNN prediction algorithm flow is proposed. The prediction model reflects the real-time characteristics of the dynamic evolution of the complex system, and highlights the independent influence of each modal component on the prediction. The effectiveness of the prediction model proposed in this paper is verified by fitting regression and error analysis. Thus, the prediction of vibration signal time series of reciprocating compressor is realized.

Prediction credibility time scale model

Multifractal and Renyi entropy

According to the knowledge of fractal theory, the correlation dimension can be expressed as a second-order information dimension, namely

$$D_2=\underset{r\to 0}{lim}\frac{\ln \sum _{i=1}^N{p}_i^2\left(r\right)}{\ln \left(r\right)}$$

(1)

In the G_P algorithm for calculating the correlation dimension, it is necessary to reconstruct the phase space of the time series first. Take the phase distribution of dimension m as an example,

$$\left[\begin{array}{cccc}x(t_1)\hfill & x(t_2)\hfill & ...\hfill & x(t_n-(m-1)\tau )\hfill \\ \multicolumn{1}{c}{x(t_1+\tau )}& x(t_2+\tau )\hfill & ...\hfill & x(t_n-(m-2)\tau )\hfill \\ \multicolumn{1}{c}{...}& ...\hfill & ...\hfill & ...\hfill \\ \multicolumn{1}{c}{x(t_1+(m-1)\tau )}& x(t_2+(m-1)\tau )\hfill & ...\hfill & x(t_n)\hfill \end{array}\right]$$

(2)

In equation (1), $t_1,t_2,L,t_n$ are the sampling time, τ is the time delay, and m is the embedding dimension, and the matrix under this phase distribution is the reconstructed phase space matrix.

The m points connected in equation (2) are used as a vector point of the continuation phase space, $k_i$ is set as the number of the point sequence (N points) in the area of the i-th unit, that is, it exists $P{\left(r\right)}_i=\underset{N\to \infty }{lim}\frac{k_i}{N}$ , considering the unit scale r (in the G_P algorithm, r is also called the similarity tolerance) and the number of points in the unit $S_i$ , the second-order information dimension can be approximated as:

$$\begin{array}{c}{D}_2=\underset{r\to 0}{lim}\frac{\ln C_m\left(r\right)}{\ln \left(r\right)}=\frac{\underset{N\to 0}{lim}\frac{1}{N^2}\sum _{i=1}^k(k_i^2-k_i)}{\ln \left(r\right)}\hfill \\ =\frac{\sum _{i=1}^k(\underset{N\to 0}{lim}\frac{k_i^2}{N^2}-\underset{N\to 0}{lim}\frac{k_i}{N^2})}{\ln \left(r\right)}\approx \frac{\sum _{i=1}^k{P}_i^2\left(r\right)}{\ln \left(r\right)}\hfill \end{array}$$

(3)

Among them, $C_m\left(r\right)$ is the proportion of “point pairs” smaller than r. According to the definition and form of the multifractal generalized dimension, we can give the definition of q-order Renyi entropy:

$$K_q=-\underset{r\to 0}{lim}\underset{k\to \infty }{lim}\left(\frac{1}{q-1}\frac{1}{k\tau }\ln \sum _{i=1}^k{P}_i^q\left(r\right)\right)$$

(4)

Obviously, q = 2, the above formula is:

$$K_2=-\underset{r\to 0}{lim}\underset{k\to \infty }{lim}\left(\frac{1}{k\tau }\sum _{i=1}^k{P}_i^2\left(r\right)\right)$$

(5)

Information entropy and prediction time

In Shannon’s information theory, the calculation formula of entropy is expressed as:

$$I_k=-\sum _i^k{P}_i\ln P_i$$

(6)

For the phase space point pair of the m-dimensional time series in equation (2), when the idea of G_P algorithm is introduced, that is, it is divided into k units of $r″\prime$ size, the information entropy of the generation probability of the unit sequence $S_i$ is $I_k$ , the entropy value must be proportional to the information required for the system unit area sequence $S_i$ (or track) determined with accuracy.

In other words, when it is known that the area sequence that the system was $s_1^{*},L,s_i^{*}L,s_n^{*}$ , and $s_{n+1}^{*}$ is the information to be predicted, this means that the difference at $r\to 0$ , $I_{k+1}-I_k$ measures the information loss of the system from time $(k-1)\tau$ to time $k\tau$ .

Then, the information loss value is Kolmogorov entropy (referred to as K-entropy), and the entropy value can be expressed as:

$$\begin{array}{c}K=\underset{\tau \to 0}{lim}\underset{k\to \infty }{lim}\left(\frac{1}{k\tau }\sum _{i=1}^k(I_{i+1}-I_i)\right)=\hfill \\ -\underset{\tau \to 0}{lim}\underset{k\to \infty }{lim}\left(\frac{1}{k\tau }\sum _{i=1}^k{P}_i\ln P_i\right)\hfill \end{array}$$

(7)

K-entropy can also be obtained by limiting the first-order Renyi entropy of equation (5). The greater the loss of characterizing information, the greater the uncertainty of the system. K-entropy is an effective tool to describe the complexity of the system.

If the change value of K-entropy is used to calculate the maximum predictable time of the system, the following methods and steps can be used:

1. Suppose the amount of information at time t is $I\left(t\right)$ , and the amount of information after passing $\mathrm{\Delta }t$ is $I(t+\mathrm{\Delta }t)$ ;

2. Then there is, $I(t+\mathrm{\Delta }t)=0$ , $I(t+\mathrm{\Delta }t)=I\left(t\right)-K\mathrm{\Delta }t$ ;

3. If the order $I\left(t\right)=1$ , and at that time, the maximum predictable time $T=1/K$ of the system is obtained.

The Literature ¹⁹ believes that the prediction within T is accurate. If the prediction time is greater than this value, it can only be said that the prediction is statistical.

Combining the above ideas based on K-entropy loss and G_P algorithm, to calculate the maximum predictive credibility scale T, because of the nature of Renyi entropy is $K_q\le K_{q+1}$ . In practical applications, due to the difficulty in calculating Kolmogorov entropy, the lower bound value of Renyi entropy K₂ is generally used as an approximate estimate of K-entropy.

Suppose the time series interval is ${\tau }_0$ , when $\tau ={\tau }_0$ , from equation (7) is that

$$\begin{array}{c}\multicolumn{1}{c}{K_2=-\underset{r\to 0}{lim}\underset{k\to \infty }{lim}\left(\frac{1}{k\tau }\sum _{i=1}^k{P}_i^2\left(r\right)\right)}\\ =-\underset{r\to 0}{lim}\underset{k\to \infty }{lim}\left(\frac{1}{k{\tau }_0}\sum _{i=1}^k{C}_m\left(r\right)\right)\hfill \end{array}$$

(8)

That exists: $C_m\left(r\right)\propto \exp (-k_m{\tau }_0{K}_2)$ . Therefore,

$$K_2=\underset{m\to \infty }{lim}\underset{r\to 0}{lim}\frac{1}{{\tau }_0}\ln \frac{C_m\left(r\right)}{C_{m+1}\left(r\right)}$$

(9)

In summary, the calculation flow chart for calculating the maximum predictable time scale of complex system time series is as follows in Figure 1.

Figure 1.

Calculation flow chart of maximum predicted confidence time of time series in systems.

According to the above calculation process, the maximum predictable time $T\approx 1/K$ of the system is obtained.

Phase space reconstruction KNN prediction

The principle of the traditional KNN prediction algorithm is to construct m sequences of the original n data and find the nearest k values, and then use the weighted value of the k points as the predicted value of the n+1th data. The algorithm expresses that the inherent complex mapping relationships of the system are all contained in historical data. It neither builds a training model for historical data nor processes historical data, and is particularly suitable for rapid classification of nonlinear problems. However, this method has poor interpretability and cannot show the evolution trend and process of the system from a theoretical perspective. The phase space reconstruction theory is introduced into the KNN prediction algorithm, and the chaotic time series prediction principle can be described as follows:

For discrete time series: $\{x\left(t_1\right),x\left(t_1\right),L,x\left(t_n\right);y\}$ , where the former n term represents the characteristics of this discrete sampling sequence at n different moments, and the last term y represents a new classification target based on the predicted data of the sequence space. After combining with the chaotic time series reconstruction algorithm, the phase spatial distribution of the data points at the sampling time can be reconstructed into the following phase distribution:

$$\left(\begin{array}{cccc}x\left(t_1\right)\hfill & x\left(t_2\right)\hfill & L\hfill & x(t_n-(m-1){\tau }_0)\hfill \\ \multicolumn{1}{c}{x(t_1+{\tau }_0)}& x(t_2+{\tau }_0)\hfill & L\hfill & x(t_n-(m-2){\tau }_0)\hfill \\ \multicolumn{1}{c}{M}& M\hfill & M\hfill & M\hfill \\ \multicolumn{1}{c}{x(t_1+(m-1){\tau }_0)}& x(t_2+(m-1){\tau }_0)\hfill & L\hfill & x\left(t_n\right)\hfill \end{array}\right)$$

Take each behavioral prediction vector, that is, construct m data vector $a=\left\{a_1{a}_{2}L{a}_{m-1}{a}_m\right\}$ as the modeling domain, where each vector is used as a sub-sequence of the vector $a_{m+1}$ to be predicted, and find its nearest k neighbors. When calculating the neighboring distance of the sub-sequence, in order to avoid the phenomenon of spatial high dimensionality Euclidean distance “big eats small,” the method of calculating the cosine distance between vectors is adopted. Through calculation, the nearest neighbor $a_m$ found from the vector $a_1,a_2,L,a_{m-1}$ is recorded as ${\beta }_1,{\beta }_2,\dots ,{\beta }_{k-1},{\beta }_k$ , because $\{x(t_1+(m-1){\tau }_0),x(t_2+(m-1){\tau }_0),L,x_n\}$ is used for prediction $x_{n+1}$ , so under this time delay, The next value of the last component of these k vectors is considered the nearest neighbor $x_{n+1}$ .

For example, if $\beta =\{x\left(t_1\right),x\left(t_2\right),\dots ,x(t_n-(m-1){\tau }_0)\}$ is a nearest neighbor of the data $a_m$ , the next sampling point value $x(t_{n+1}-(m-1){\tau }_0)$ based on the sampling interval is taken as the next nearest neighbor $m_1=x(t_{n+1}-(m-1){\tau }_0)$ of $x_{n+1}$ , so that the k nearest neighbors $m_1,m_2,\dots m_{k-1},m_k$ of $x_{n+1}$ are obtained, and then the average of the k numbers is the predicted value can be obtained. which is:

$$x_{n+1}=\frac{\sum _{i=1}^k{m}_i}{k}$$

(10)

MSS dynamic SPA_KNN prediction model based on prediction credibility scale

SPA

Dr. Karjalainen ²⁰ from Kuopio University in Finland proposes a nonlinear de-trend method—Smoothness Priors Approach (SPA). It assumes that the original data signal, namely the time series Z,

$$Z=Z_{\mathrm{s}}+Z_{\mathrm{t}}$$

(11)

In equation (11), Z_s is a stationary term, Z_t is a nonlinear low-frequency trend component,

$$Z_{\mathrm{t}}=H\theta +\mathrm{v}$$

(12)

In equation (12), $H\in R_{M\times N}$ is the observation matrix, $\theta \in R_M$ is the regression parameter, and v is the observation error, thus the task is transformed into an optimization method to estimate the parameter θ, so that ${\hat{Z}}_{\mathrm{t}}=H\hat{\theta }$ can be used to estimate the trend term in the original signal. The usual method to estimate the parameter is the least square method. The SPA method adds the differential item $\Vert D_{\mathrm{d}}(H\theta )\Vert$ to the optimization process and minimizes it to ensure that $H\theta$ filters out the trend part of the signal:

$${\hat{\theta }}_{\lambda }=\arg \underset{\theta }{min}\{{\Vert H\theta -Z\Vert }^2+{\lambda }^2{\Vert D_d(H\theta )\Vert }^2\}$$

(13)

Where $\lambda$ is the positive regularization parameter and D_d is the discrete d order differential operation matrix. The solution method of the matrix D_d is as follows:

Set sequence Z contain N local extremun points, which can be represented by column quantities $R=[R_1,R_2,\cdots ,R_N{]}^T\in R^N$ , the first-order trend of R is expressed in a discrete way as follows:

$$R_1=[R_2-R_1,R_3-R_2,\cdots ,R_N-R_{N-1}{]}^T\in R^{(N-1)}$$

The second-order trend of R is expressed in a discrete way as follows:

$$\begin{array}{c}{R}_2=[R_3-R_2-(R_2-R_1),R_4-R_3-(R_3-R_2),\cdots ,\hfill \\ R_N-R_{N-1}-(R_{N-1}-R_{N-2}){]}^T\in R^{(N-2)}\hfill \end{array}$$

(14)

By that analogy, the discrete representation of any order trend of R can be obtained, that is, the d-order derivative of R can be represented by D_d:

$$D_{\mathrm{d}}=\left[\begin{array}{c}d(R_d{)}_1/d{R}_1\cdots d(R_d{)}_1/d{R}_N\hfill \\ \multicolumn{1}{c}{⋮⋮}\\ \multicolumn{1}{c}{d(R_d{)}_{N-d}/d{R}_1\cdots d(R_d{)}_{N-d}/d{R}_N}\end{array}\right]$$

(15)

The solution in Formula (5) is

$$\begin{array}{c}{\hat{\theta }}_{\lambda }=(H^{T}H+{\lambda }^2{H}^T{D}_{\mathrm{d}}^T{D}_{d}H{)}^{-1}{H}^{T}Z\hfill \\ {\hat{Z}}_{\mathrm{t}}=H{\hat{\theta }}_{\lambda }\hfill \end{array}$$

(16)

Where t is the estimate of the trend item that needs to be removed. The matrix H can be obtained by analyzing the characteristics of the original signal Z. In order to facilitate the analysis, H adopts the identity matrix $I\in R^{N\times N}$ . For the matrix D_d, the aperiodic trend item in the signal can be well estimated when the order is 2, so the order of D_d is 2, $D_2\in R^{(N-2)\times N}$ can be expressed as:

$$D_2=\left[\begin{array}{c}1-210\cdots 0\hfill \\ \multicolumn{1}{c}{01-21\cdots ⋮}\\ \multicolumn{1}{c}{⋮⋮⋮⋮0}\\ \multicolumn{1}{c}{0\cdots 01-21}\end{array}\right]$$

(17)

Then, the stationary part of the original signal after removing the trend item can be expressed as

$${\hat{Z}}_{\mathrm{s}}=Z-H{\hat{\theta }}_{\lambda }=[I-{(I+{\lambda }^2{D}_2^T{D}_2)}^{-1}]Z=\mathrm{LZ}$$

(18)

Where $L=[I-{(I+{\lambda }^2{D}_2^T{D}_2)}^{-1}]$ , then there is ${\hat{Z}}_{\mathrm{s}}=\mathrm{LZ}$ .

MSS

The DFA (Detrend Fluctuation Analysis) analysis method was first proposed by Peng et al. ²¹ based on the DNA mechanism. This method uses the scale index to filter out the various trend components in the molecular chain sequence of the DNA, and successfully detects the degree of relevance and possibility in the sequence. Kantelhardt et al. proposed a multifractal de-trend fluctuation analysis method based on DFA, abbreviated as MFDFA. This method can easily calculate the scale index and multifractal spectrum, and perform effective analysis. Based on the MFDFA correlation analysis, it has been widely used in life sciences, environmental sciences, economics, and mechanical fault diagnosis.^22–24 Among them, Among them, Hurst index correlation analysis has been proven to be one of the most reliable and important tools for studying the long-term correlation of non-stationary time series. The specific algorithm steps of MFDFA are as follows:

Remove the mean of the sequence $\{x_k|k=1,2,...,N\}$ and construct a summation sequence

$$Y\left(i\right)=\sum _{k=1}^i(x_k-〈x〉),i=1,2,\cdots ,N$$

(19)

$$〈x〉=\frac{1}{N}\sum _{k=1}^N{x}_k$$

(20)

$Y\left(i\right)$ is divided into $N_s=\mathrm{int}(N/s)$ non-overlapping sub-intervals with a length of s, and obtain $2N_s$ sub-intervals in order to ensure the integrity of the information. Fitting the k-order polynomial by the least square method, we can get:

Then $m=1,2,...N_s$ ,

$$F^2(s,m)=\frac{1}{s}\sum _{i=1}^s{\{Y[(m-1)s+i]-y_m(i)\}}^2$$

(21)

At that time, $m=N_s+1,N_s+2,...2N_s$

$$F^2(s,m)=\frac{1}{s}\sum _{i=1}^s{\{Y[N-(m-1)s+i]-y_m(i)\}}^2$$

(22)

Calculate the mean value of $F^2(s,m)$ and express it as a q-order volatility function $F_q\left(s\right)$

$$F_q(s)={\left\{\frac{1}{2N_s}\sum _{m=1}^{2N_s}{[F^2(s,m)]}^{\frac{q}{2}}\right\}}^{\frac{1}{q}}$$

(23)

When q > 0, $F_q\left(s\right)$ reflects the large volatility trend of the time series; and q = 2, MFDFA is degenerated into DFA; when q < 0, $F_q\left(s\right)$ reflects the small volatility trend of the time series; and q = 0

$$F_0\left(s\right)=\exp \left\{\frac{1}{4N_s}\sum _{m=1}^{2N_s}\ln \left[F^2(s,m)\right]\right\}\propto s^{h(0)}$$

(24)

If the time series are self-similar, then the q-order fluctuation trend of $F_q\left(s\right)$ and time scale s satisfy the following power law relationship:

$$F_q(s)\propto s^{h(q)}$$

(25)

In equation (25), $h\left(q\right)$ is the generalized Hurst exponent. By analyzing the double logarithmic function of $F_q\left(s\right)$ and s, the scale of the fluctuation function can be obtained. When analyzing the time series with multifractal characteristics, $h\left(q\right)$ appears as a function that depends on q.

The Hausdorff dimension of the multifractal set of the time series is:

$$f(\alpha (q))=\underset{\delta \to 0}{lim}\frac{\sum _i{\mu }_i(q,\delta )\ln [{\mu }_i(q,\delta )]}{\ln \delta }$$

(26)

The overall singularity average of the sequence is:

$$\alpha (q)=\underset{\delta \to 0}{lim}\frac{\sum _i{\mu }_i(q,\delta )\ln [p_i(\delta )]}{\ln \delta }$$

(27)

The relationship of above two and $h\left(q\right)$ can be expressed as:

$$\alpha =\frac{d\tau \left(q\right)}{\mathrm{dq}}=h(q)+\mathrm{qh}\prime (q)$$

(28)

$$f\left(\alpha \right)=q\alpha -\tau \left(q\right)=q[\alpha -h\left(q\right)]+1$$

(29)

Multifractal Singulity Spectrum (MSS) can be obtained by the MFDFA method. The spectral parameters can describe the dynamic behavior of the multifractal time series. At the same time, the singularity index reflects the growth probability of the fractal in a small area, and characterizes the degree of unevenness of the time series distribution in the local probability measure.

Figure 2 shows the morphological characteristics and physical meaning of the “inverted bell-shaped” singularity spectrum. The slope ${\alpha }_{-\infty }$ of the left end of the singularity spectrum corresponds to the maximum volatility of the singular index at that point $q\to +\infty$ ; The slope ${\alpha }_{+\infty }$ of the right end of the singularity spectrum corresponds to the minimum volatility of the singular index at that point $q\to -\infty$ . The spectral width $\mathrm{\Delta }\alpha ={\alpha }_{-\infty }-{\alpha }_{+\infty }$ reflects the uneven degree of the probability measure distribution of the time series in the entire fractal structure; the slope of the extreme point $q=0$ makes the extreme value of the multifractal singular spectrum equal to the Hausdorff dimension. $\alpha \left(0\right)$ reflects the uneven distribution of the local probability measure in the time series when q = 0.

Figure 2.

Characteristic parameters of MSS.

Improved algorithm prediction process

Starting from the nature of prediction and the state characteristics of complex mechanical systems, this paper introduces the concept of predictable time scale and SPA algorithm into the KNN prediction algorithm, and selects feature parameters based on the singular spectrum distribution to construct the sequence to be predicted, a dynamic SPA_KNN prediction method based on the MSS eigenvalues of the reciprocating compressor vibration time series is proposed. The basic steps are as follows:

1. Extract the “short-term” periodic vibration signals of the measuring points at equal intervals, and construct the characteristic time series X₀ to be predicted, $X_0=\{x_1,x_2,\dots ,x_p;x_i,x_{i+1},\dots ,x_{i+p-1}\cdots x_{n-p-1},x_{n-p},\dots ,x_n\}$ ;

2. Perform SPA decomposition on X₀, and calculate the predicted time scale of the trend and de-trend component of each sampling segment according to the sampling period [T₁, T₂…T_n];

3. Extract the “short-term” multifractal singular spectrum eigenvalues of the trend and de-trend component, and form a feature sequence according to the extraction interval;

4. Create the initial modeling domain of the component spectrum sequence of different trend and de-trend items, and calculate the embedding dimension m and time delay τ according to the phase reconstruction algorithm;

5. Calculate the k neighbors of the phase distribution based on the prediction model, and obtain the predicted spectrum value with the mean deviation;

6. Repeat steps (4)–(5), the modeling domain rolls forward, and then predict;

7. Taking min[T₁, T₂…T_n] of the trend and de-trend item component as the cut-off time to form the predicted characteristic spectrum sequence.

The MSS parameter prediction flowchart of vibration time series based on the maximum credible scale is shown in Figure 3.

Figure 3.

Prediction process of vibration time series based on maximum prediction scales for MSS.

Reciprocating compressor prediction example

According to the above process, select 2D12 type reciprocating compressor measuring point data (first and second bearing failures), combined with the development and evolution trend of its vibration time series, analyze the influence of the sequence length, sampling frequency and other parameters on the predicted credible time, Based on parameter optimization analysis and taking the gap fault of the two-stage connecting rod large-end bearing as an example, the effectiveness of the prediction method is analyzed.

A 2D12-70/0.1–13 double-acting reciprocating compressor was studied in Figure 4. Its shaft power is 500 kW, the exhaust volume is 70 m³/min, the exhaust pressure is 0.2746–0.2942 MPa, the piston stroke is 240 mm, the motor speed is 496 r/min, and the corresponding motor rotation frequency is 8.27 Hz. In order to facilitate the selection and analysis of sensitive measurement points for bearing vibration signals, eight measurement points were arranged at the cylinder end cover (inside and outside), crankcase, and crosshead to install the integrated circuit piezoelectric (ICP) acceleration sensors (CT1010LC), respectively, as shown in Figure 5. In this paper, we mainly study the reciprocating compressor secondary connecting rod sliding bearing clearance fault, and the time domain and spectrum plots obtained from the normal bearing clearance and excessive large head bearing clearance were collected at measurement points 6, 7, and 8, respectively, as shown in Figures 6 and 7.

Figure 4.

Two-stage double-acting reciprocating compressor of type 2D12: (a) photos of the test bench; (b) and (c) sensor on the top of the crankcase and the crosshead guide surface.

Figure 5.

Distribution of measurement points for reciprocating compressor.

Figure 6.

Signal time-frequency diagram for 6, 7, and 8 measurement points at normal condition.

Figure 7.

Signal time-frequency diagram for 6, 7, and 8 measurement points at fault clearance condition.

Since the bearing clearance fault signal is obtained through a complex transmission path on the surface of the body, the sample signal sampling point should be set to make the shortest transmission path from the source of the excitation force to the measurement point. That is, in the same failure mode, compared with the normal state, the energy increment of the sensitive measurement point will increase more than other measurement points, and this principle can be used to determine the sensitive measurement points in various bearing clearance fault modes. Table 1 shows the energy gain of each measurement point before and after the compressor secondary connecting rod bearing clearance failure, and the locations of the different fault sensitive measurement points are indicated in bold. It can be seen from Table 1 that for the secondary connecting rod large head bearing clearance fault sensitive measurement point location for the measurement point No. 6.

Table 1.

Vibration signal energy difference for each measurement point under different fault states.

Fault location	Energy difference of vibration signal at each sampling point
	1	2	3	4	5	6	7	8
Primary connecting rod big head	128.1	145.8	324.5	55.9	98.9	65.2	102.3	57.6
Primary connecting rod small head	224.4	103.3	67.5	45.4	95.4	45.2	159.5	54.1
Secondary connecting rod big head	78.3	97.7	114.2	67.3	64.9	315.9	143.5	132.4
Secondary connecting rod small head	86.9	94.9	57.2	76.8	56.2	207.8	132.4	298.8

In the simulation experiments carried out, the clearance states of each component of the secondary connecting rod of the 2D12 reciprocating compressor were studied, and a total of four states of normal clearance, bearing and crosshead failure were set. According to the technical manual for the use of piston compressors and the API-618 standard, a typical moderate wear failure was simulated for the secondary bearing assembly in the failure test, that is, the clearance of the secondary connecting rod small head bearing and crosshead were set to 0.15 mm, and the clearance of the secondary connecting rod large head bearing was 0.25 mm. The aim is to find an effective classification and identification method for the difficult to identify moderate wear faults, and to provide the underlying diagnostic algorithm support for the subsequent life prediction.

Prediction time and parameter analysis

Influence of parameters on the prediction time scale

According to the calculation process based on the chaotic time series model, taking the sampling sequence of the vibration signal of the moderate wear-sensitive measuring point of the secondary large-head bearing as an example, the influence parameters of the maximum predictive credibility scale are analyzed.

Figure 8 shows the saturation of the second-order Reyin entropy K₂ with the embedding dimension m under different similarity tolerances r for the sampling sequence, that is, there is a downward jump near the embedding dimension m = 13. Search for the K₂ value that does not change in the vicinity of m, as shown in Figure 9 (where r is the sequence index, not the actual value), and the calculated prediction time is $T\approx 1/K_2=32$ , that is, the prediction credible period is about 32 sampling intervals, If the sampling interval is increased, the mid-to-long-term prediction and life analysis of the vibration state can be realized.

Figure 8.

The relation between entropy of K₂ and embedded dimension m.

Figure 9.

K ₂ value of the saturation similar tolerance r.

Considering the periodicity of the compressor vibration signal, combining the key phase signal and sampling frequency to divide the signal into different periods, considering the sampling frequency and sequence length, and analyzing the influence of each parameter on the embedding dimension m and the maximum prediction time T of different sampling sequences, The results are summarized in Table 2, which V_t represents the sampling interval.

Table 2.

Influence parameters analysis of prediction dependable scale for vibration signal.

Sequence length	Sampling frequency (KHz)	Saturation dimension $m_0$	K 2 estimated value	Maximum prediction scale	Sampling period
12050	25	14	0.030	32Vt	4
6024	25	14	0.029	35Vt	2
6024	50	13	0.032	33Vt	1
1205	10	13	0.028	34Vt	1
1205	5	10	0.020	50Vt	0.5
603	5	8	0.016	62Vt	1

It can be seen from Table 2 that for bearing moderate wear faults, when the sampling sequence meets the fluctuation trend and resolution, ensure that one vibration period is sampled, that is, the sequence length is more than 1000 points, and the maximum predictable scale is basically saturated at about 30 sampling interval which reflects the period-like characteristics of the vibration sequence, and the length of the sequence is consistent with the minimum suitable time sequence length required in the phase space reconstruction theory of the time sequence.

Fix the sequence length and sampling frequency, select representative bearing measurement points, analyze and compare the normal state and the fault state (fault type and working condition simulation are the same as described in the previous section), and observe that the various parameters are paired under different fault conditions. The results of the maximum prediction scale are summarized in Table 3.

Table 3.

Influence parameters analysis of prediction dependable scale for different working condition.

Data category	Data length	Saturated dimension $m_0$	K 2 estimated value	Maximum prediction scale	Sampling period
Secondary normal bearing	1205	8	0.016	61Vt	1
Slightly worn secondary bearing	1205	11	0.027	39Vt	1
Moderate wear of secondary bearing	1205	14	0.031	32Vt	1
Severe wear of secondary bearing	1205	–	–	–	1

Table 3 shows that for different types of failures, the predicted credible time values differ significantly. For the bearing measuring point signal, the prediction credibility scale under normal conditions is relatively large. In the case of mild to moderate failure, due to the weakening of the quasi-period feature and the increase of the scale-free area of the fractal feature, the entropy value increases, and the prediction credibility time is shortened. In a severely worn state, the entropy value of the system does not converge and tends to a chaotic state, which is unpredictable.

It can be seen that the prediction model for predicting the credible time scale reflects the dynamic evolution trend of the system, and the model for predicting the credible time scale is a prerequisite guarantee for improving the credibility and accuracy of the prediction algorithm.

SPA_KNN algorithm based on predictive credibility scale

Through the prediction parameters analysis of the above time series, for the vibration signal prediction, limited by the prediction credibility scale, even if the vibration acceleration value is obtained at a finite time with a small error, the finite discrete prediction data value cannot be used for feature extraction and analyse.

Therefore, we consider that when constructing the sampling sequence, the “short-term” periodic sampling sequence at each sampling moment is used to replace the instantaneous value to construct the sequence to be modeled, so as to ensure that each “short-term” sampling moment contains at least one vibration period while obtaining a representative sampling sequence for a certain period of time. Two components of trend item and de-trend item are obtained through SPA decomposition, and characteristic parameters such as the spectrum width $\mathrm{\Delta }\alpha$ , peak difference $\mathrm{\Delta }f$ , center value ${\alpha }_0$ , and symmetry parameter B of the multifractal singular spectrum of the two components are calculated. And the robust and separable spectral eigenvalues are extracted, consider the predicted credibility scale of each component to construct the phase space sequence based on the time series of the spectral eigenvalues at each “short-term” sampling time, and construct the variable parameter KNN modeling domain to achieve prediction and spectral feature extraction with a high degree of confidence.

Prediction examples and method evaluation

Following the SPA_KNN prediction process proposed in this section, the error analysis of the prediction results is carried out through examples.

SPA_KNN algorithm prediction example

Take the 2D12 type reciprocating compressor with a moderate wear failure of the secondary big head bearing as an example. To ensure the effective length of the sampling, while taking into account the evolution trend of the nonlinear chaotic sequence and the time limit factor of the sampling interval for the effectiveness of the prediction, each interval of the long-term sampling sequence 2 min intercepts about 2410 sampling points in two cycles (sampling frequency 10 kHz), and obtain 2-h interval sampling data. The number of SPA decompositions in each sampling period is two, and the maximum predicted credible time of each component sequence is calculated as $T=\left[43V_{t}30V_t\right]$ , we can see that at least the predicted value within 60 min after the calculation point is credible.

Since almost all the computations in the KNN algorithm occur in the classification phase, and the magnitude of the K value is directly related to the final classification accuracy. In general, the larger the K value, the more accurate the classification result. If the K value is chosen too large and the sample set is small, it will increase the computation and lead to the decrease of classification accuracy. The value of K is too small, but the number of neighbors obtained is too small, which will also reduce the classification accuracy. There is no definite formula for the K value, and generally by manually taking multiple K values for multiple training, the K value with the smallest classification result error is taken. In this paper, K values are taken as shown in the following table, and K = 3 is selected by the accuracy magnitude of K values in Table 4.

Table 4.

Diagnostic accuracy at different K values.

K	1	3	5	7	9
Accuracy	95%	97.5%	92.5%	92.5%	92.5%

The phase space of the sequence is updated rolling, and a new modeling domain is constructed to re-predict the sequence to reproduce the dynamic evolution trend of the system in real time.

In the multifractal theory, it can be known that the singular spectrum index has an equivalent transformation relationship with the generalized dimension, and it is an effective tool for analyzing the partial shape characteristics of the signal. The effective extraction and analysis of the characteristic parameters based on the spectrum shape has been widely used in many fields.

Singular spectrum characteristic parameters: ${\alpha }_0$ , $\mathrm{\Delta }\alpha$ , B, $\mathrm{\Delta }f=f({\alpha }_{+\infty })-f({\alpha }_{-\infty })$ , their definition and physical meaning are: ${\alpha }_0$ is defined as the central value, representing the long-range correlation of the signal, the larger the value, the stronger the correlation; $\mathrm{\Delta }\alpha$ is the spectral width, describing the fluctuation degree of the vibration signal, the value The larger the value, the more severe the fluctuation; B is called the degree of symmetry, which is defined as the fitting value of the nearby spectral curve. When it is greater than 0, the curve shape tilts to the left, which corresponds to a small singularity index. On the contrary, the curve has strong singularity; $\mathrm{\Delta }f=f({\alpha }_{+\infty })-f({\alpha }_{-\infty })$ is the peak difference. It is a true reflection of the proportion of the change rate of the peak value of the vibration signal. If the value is less than 0, the number of the largest probability subset is greater than the number of the smallest probability subset.

The sensitivity of feature parameters to state is an important indicator for evaluating state feature extraction methods. How to improve the separability between different state features is the prerequisite for the selection and modeling of evaluation indicators. First, extract the 30-period multifractal singular spectrum feature quantity of the sensitive measuring point of the secondary connecting rod bearing moderate wear fault, and compare and analyze the stability and separability of each feature vector. The result is shown in Figure 10.

Figure 10.

Eigenvalue of multifractal singular spectrum for vibration signals.

Before prediction, according to the analysis of the stability and separability of the MSS eigenvalues in Figure 10, the spectral center value and the spectral width are stable, and each eigenvalue has obvious nonlinear fractal characteristics, and the phase space reconstruction can be used. The improved KNN prediction model is constructed to predict, and the obtained prediction results can be further analyzed to achieve fault classification, performance evaluation, and life prediction.

The predicted values obtained are shown by the dotted lines in Figures 11 and 12. From the predicted trend of the feature parameters of the two components, it can be seen that the spectral parameters obtained by the prediction algorithm show the fluctuation and fractal characteristics of the vibration signal. Spectral width and central value of the two spectral parameters, for different components, the predicted value and actual value show good follow ability. In particular, the predictions within the first 20–25 prediction credible time intervals, that is, the predicted values within 50 min, have a high degree of agreement, which proves that the predicted values obtained under this sampling interval have a high degree of credibility.

Figure 11.

Comparison between the predicted values and the measured values of $a_0$ for trend and de-trend items.

Figure 12.

Comparison between the predicted values and the measured values of $\mathrm{\Delta }a$ for trend and de-trend items.

Method comparison

The following is a comparative analysis of the three prediction methods of KNN, SPA_KNN and SPA&PSR_KNN, respectively, using the same measured data and regression analysis and error analysis for evaluation. Both methods select two spectral characteristic parameters (spectral width $\mathrm{\Delta }a$ and spectral center value $a_0$ )with better separability and robustness as the evaluation objects. Among them, the regression analysis adopts the method of linear fitting regression to construct the fitting equation $Y_{\mathrm{predict}}=R\times X_{\mathrm{actual}}+\mathrm{ERR}$ , where $Y_{\mathrm{predict}}$ and $X_{\mathrm{actual}}$ respectively represent the predicted value and the measured value, and ERR is the error term; R is the slope of the curve. The closer its value is to 1, the better the curve prediction effect. The comparison results of regression analysis are summarized in Figures 13 and 14.

Figure 13.

Regression analysis of $a_0$ for different forecasting methods: (a) trend item and (b) de-trend item.

Figure 14.

Regression analysis of $\mathrm{\Delta }a$ for different forecasting methods: (a) trend item and (b) de-trend item.

It can be seen from Figures 13 and 14 that the prediction regression curve slopes of the SPA&PSR_KNN method proposed in the four modes are all higher than those of the first two methods. Further analysis can also notice that the method based on SPA decomposition is better than the traditional KNN method in predicting the center value ${\alpha }_0$ of the singular spectrum with better stability, especially the spatial reconstruction type SPA&PSR_KNN prediction method within the effective prediction confidence time, the prediction effect of the spectrum center value ${\alpha }_0$ is the best.

In the error analysis, two indicators are selected: Standard Deviation (STD) and Root Mean Square Error (RMSE). The errors of each modal parameter are shown in Tables 5 and 6.

Table 5.

Comparison between prediction and measurement of $a_0$ error for vibration signal.

		Multifractal singular spectrum center value $a_0$
		Trend item	De-trend item
KNN	STD	0.1482	0.1242
	RMSE	0.1787	0.1499
	Time(s)	1017	1013
SPA—KNN	STD	0.0827	0.1004
	RMSE	0.0759	0.0954
	Time(s)	579	572
SPA&PSR—KNN	STD	0.0460	0.0477
	RMSE	0.0596	0.0527
	Time(s)	396	398
SPA—BP	STD	0.0571	0.0779
	RMSE	0.0691	0.0807
	Time(s)	487	492
SPA—SVM	STD	0.0463	0.0485
	RMSE	0.0601	0.0636
	Time(s)	447	451

Table6.

Comparison between prediction and measurement of $\mathrm{\Delta }a$ error for vibration signal.

		Multifractal singular spectrum width $\mathrm{\Delta }a$
		Trend item	De-trend item
KNN	STD	0.4637	0.5719
	RMSE	0.4659	0.2467
	Time(s)	1023	1019
SPA—KNN	STD	0.1617	0.1613
	RMSE	0.1756	0.2520
	Time(s)	637	631
SPA&PSR—KNN	STD	0.0671	0.0967
	RMSE	0.0352	0.1802
	Time(s)	425	431
SPA—BP	STD	0.0952	0.1192
	RMSE	0.1197	0.2070
	Time(s)	516	521
SPA—SVM	STD	0.0719	0.0998
	RMSE	0.0810	0.1886
	Time(s)	482	489

It can be seen from Table 5 and Table 6 that the error between the SPA&PSR_KNN predicted sequence and the measured sequence is the smallest. After singular spectrum feature extraction and analysis, Considering the error of the predicted value within 20–25 sampling intervals, the prediction of vibration characteristics in the future 40–50 min can be achieved. The predicted spectrum value $\mathrm{\Delta }\alpha -{\alpha }_0$ is close to the singular spectrum clustering center in a severely worn state, and its performance degradation has been severely degraded, which is consistent with the actual state of the fault. At this time, it should be shut down for maintenance. In order to better verify the superiority of the proposed method in this paper, SVM and BP methods were used to analyze the results according to the research process above, and the results are shown in Tables 5 and 6. Based on the data values in the table, it can be seen that the error between the SPA&PSR_KNN predicted sequence and the measured sequence is smaller than that of both SPA-BP and SPA-SVM predictions, and the computation time of the proposed method is obviously better than other method, further proving that the results predicted using the SPA&PSR_KNN method will be more accurate and faster when performing fault prediction for the fault state feature vector dataset obtained by the method in this paper.

Conclusions

Starting from the factors affecting the accuracy of complex nonlinear predictions, this paper analyzes the applicability of non-parametric prediction methods in nonlinear classification and prediction, studies the evolution of system chaotic time series, and proposes a combination of SPA and KNN based on the maximum predictive credibility scale, the multifractal singular spectrum feature value is used as the prediction parameter to realize the prediction of the time series of the vibration signal of the reciprocating compressor.

1. Based on the theory of phase space reconstruction, a maximum predictive credible time algorithm is proposed from the perspective of information entropy minimization, and an improved time series phase space reconstruction KNN prediction model is constructed.

2. On the basis of singular spectrum distribution analysis, the spectrum parameter prediction sequence is constructed by the improved phase space reconstruction KNN algorithm, and the central value of the singular spectrum spectrum and the spectrum width of each component are selected as the prediction parameters, and the K neighbor dynamic extrapolation is performed within a credible scale, which improves the interpretability of prediction from the perspective of system evolution

3. Construct the variable parameter modeling domain with the predictive credibility scale as the precondition, and propose a complete SPA&PSR_KNN prediction algorithm flow. The example verification and error analysis show that compared with the traditional KNN algorithm, BP and SVM, the prediction results of the spectral characteristic parameters obtained by the algorithm have smaller errors and higher reliability, and the computation time of the proposed method is obviously better than other method. Moreover, if the spectral characteristic parameters can be combined with the fault classification threshold, which can further realize the degradation state assessment and life cycle analysis.