NGO-MRE combined with DELM for bearing fault diagnosis

2.1. Multi scale range entropy

RE is a nonlinear method proposed by Amir in 2018, which can more accurately and effectively measure complex nonlinear time series. However, due to the fact that rolling bearing fault signals are generally distributed on multiple scales, MRE is used to obtain fault information at multiple scales [24]. The specific formula is as follows:

Reconstruct the phase space and give a time series of length $L$, is $x_i=x_1,x_2,x_3,\cdots ,x_N$. Select $m$ data points and reconstruct the sequence in the original order to obtain $Z_i^{m,d}$, as shown in the calculation:

where, $m$ and $d$ represent the embedding dimension and delay time, respectively. When $d=1$, the original data is split into a new sequence $Z_i^m=\{x_i,x_{i+1},x_{i+2},\text{...},x_{i+m-1}\}$ to obtain the reconstructed sequence.

Using rescaled range analysis, calculate the distance $D_{range}$ between $Z_i^m$ and other subsequences $Z_j^m$. The formula is:

Compare the difference between each $D_{range}\left(Z_i^m,Z_j^m\right)$ and similar capacity $r$, and calculate $B_i^{m,r}$ that satisfies the sum of $D_{range}\left(Z_i^m,Z_j^m\right)>r$, as follows:

where $j\ne i$, $\mathrm{\Psi }\left(\cdot \right)$ is the Heaviside function, calculated as follows:

Add all A and take the mean to obtain the expression for B:

Repeat the above steps to obtain $B_{m+1}^r$ from $m=m+1$, and then calculate RE using the following formula:

Coarse grained $x_i=\{x_1,x_2,x_3,\text{...},x_N\}$ using the following formula:

where, $\tau$ is the scaling factor, $j$ is the coarse-grained time series, and $\left[\right]$ is the rounding function. When $\tau =1$, $y_j^{\left(\tau \right)}$ is the original time series.

Calculate the coarse-grained time series for each RE, and the MRE is obtained from the following formula:

2.2. NGO optimization algorithm

NGO is a new optimization algorithm proposed by Mohammad et al. [16] in 2021, and its mathematical model mainly consists of the following three parts:

1) Initialize the population. Randomly distribute the population in the search space, and the population matrix is:

where, $X_i$ represents the position of the $i$ northern eagle, $N$ represents the population size, $m$ represents the dimension for solving the problem, and $x_{i,j}$ represents the position of the $i$ dimension of the $j$ northern eagle.

The objective function value is expressed as a vector:

where, $F$ is the objective function vector of the population, and $F_i$ is the objective function value of the $i$ northern eagle.

2) Phase 1: Search and recognition phase. The northern eagle randomly selects a prey in the search space and conducts a global search to obtain the optimal region. Its mathematical model is as follows:

$where$, $P_i$ represents the prey position of the $i$ northern eagle, $p_{i,j}$ represents the prey position of the $i$ northern eagle in the $j$ dimension, and $k$ is a random integer of $\left[1，N\right]$; $X_i^{new,P1}$ is the new position of the $i$ northern eagle, $x_{i,j}^{new.P1}$ is the $j$ new position of the $i$ northern eagle, and $F_i^{new,P1}$ is its corresponding fitness value; The parameters $r$ and $I$ are new random numbers updated during the search and iteration process, where the $r$ interval is a random number of [0, 1] and the value of $I$ is 1 or 2;

3) Phase 2: Chasing and Escaping. During the pursuit of prey by the Northern Goshawk, the prey attempts to escape, and the Northern Goshawk will chase at an extremely fast speed and ultimately capture the prey. Assuming that the hunting range of the Northern Goshawk is radius $R$. This simulation method increases the local search ability of NGOs in the search space, and its mathematical model is:

where, $x_{i,j}^{new,P2}$ represents the $j$ updated position of the $i$ northern eagle after $P2$ update; $t$ is the current number of iterations, and $T$ is the maximum number of iterations; $x_i^{new,P2}$ refers to the new position of the $i$ northern eagle in $P2$, and $F_i^{new,P2}$ is the objective function value of the $i$ northern eagle updated by $P2$.

2.3. Deep extreme learning machine

Extreme Learning Machine (ELM) is a single hidden layer feedforward neural network applied to supervised learning problems. Given $N$ samples $\left(x_i,y_i\right)\in R^n\times R^m$, $n$ and $m$ represent the number of features and categories, respectively [25]. Its model is as follows:

where, ${\beta }_i=[{\beta }_{i1},{\beta }_{i2},\cdots ,{\beta }_{im}{]}^T$ is the output weight of the $i$ hidden layer node, $a_i=[a_{i1},a_{i2},\cdots ,a_{in}{]}^T$ is the input weight of the $i$ hidden layer node, $b_i$ is the threshold of the $i$ hidden layer node, and $G\left(•\right)$ is the activation function. This article uses the sigmoid function.

The matrix expression of ELM:

By using the least squares solution $\widehat{\beta }$, make:

Among them, $H$ is the output matrix of the hidden layer, $T$ is the expected output matrix of the sample, and the formula for the output weight is:

Among them, $H^{+}$ is the generalized inverse matrix of $H$.

Applying the autoencoder AE to ELM, its weight $\beta$ is converted to:

Stacking multiple extreme learning machine autoencoders ELM-AE to form a DELM, the weights and biases of each ELM-AE random input will affect the final training results. In subsequent experiments, the Northern Eagle Optimization Algorithm will be introduced to optimize the parameters of the DELM.

3.1. Experimental data

The experimental data utilizes the publicly available Southeast University dataset, with bearing operating conditions categorized as 20-0 (20 Hz, 1200 rpm, unloaded, 0 V, 0 N/m) and 30-2 (30 Hz, 1800 rpm, loaded, 2 V, 7.32 N/m). The bearing conditions are divided into five types, namely: ball failure, inner ring failure, outer ring failure, inner and outer ring combination failure, and health status. The experimental setup is shown in the figure.

The data in its public dataset file has a total of 8 columns of vibration signals, and the 4th column of data in each type of dataset is selected for experiments in this paper. The summarized data types are 10 in total, and the selected data are labeled as follows in Table 1.

Fig. 1Southeast university bearing test stand

In Table 1, the bearing conditions selected are categorized into ball failure, inner race failure, outer race failure, outer race combined failure, and healthy state. The working environments are classified into two conditions: 1) rotational speed of 20 Hz (1200 rpm) with no load (0 V, 0 N/m), and 2) rotational speed of 30 Hz (1800 rpm) with a load of 2 V (7.32 N/m). These categories are presented across five rows. The reasons for these selections are as follows:

Representativeness: The chosen fault types cover common bearing failures, each with distinct vibration characteristics, allowing for a comprehensive evaluation of the diagnosis method.

Combination Failure: The inclusion of inner and outer ring combination failure is crucial as it presents more complex vibration patterns, making it harder to identify compared to single faults. This adds practical challenge to the experiment, as combination failures are common in real-world scenarios.

Health Status as Baseline: The health status serves as a reference for normal conditions, enabling comparison against fault states.

Data Diversity: These 5 rows offer sufficient fault variety while maintaining manageable data complexity for efficient analysis.

3.2. Experimental steps

1) Using VMD to decompose the fault signal and obtain all IMF.

2) Select a portion of the intrinsic mode components and use the NGO optimization algorithm to optimize the $m$ and $r$ values of multi-scale range entropy based on the minimum margin factor.

3) Using the optimal sum value, calculate the range entropy of all data at different scales to obtain feature signals.

4) Using ISOMAP to reduce the dimensionality of feature signals and remove excess fault signals;

5) Using the DELM optimization algorithm of NGO to identify and classify the fault features after dimensionality reduction.

Fig. 2Experimental flow chart

3.3. Description of NGO algorithm parameter settings

Population Size: Set to 50 individuals. This setting aims to balance optimization stability and computational efficiency.

Maximum Number of Iterations: Set to 150 iterations, as the algorithm fully converges at this iteration count.

Search Space Boundaries: The search range for the optimization variable (DELM network weights) is limited to [–5, 5], based on the characteristics of the Sigmoid activation function.

Optimization Variable Dimension: The dimension is automatically determined by the DELM network structure.

Objective Function: The optimization objective is to minimize the classification error rate of the training set.

3.4. Parameter settings of Harris eagle optimization (HHO) algorithm

3.5. Fault diagnosis analysis

Select the bearing status data with a rotational speed of 20 Hz- no load (20-0) (including health status, ball faults, combined faults of inner and outer rings, inner ring faults and outer ring faults), process the original signal using VMD decomposition, select the default values of parameters, the number of modes $k=$ 4 and the penalty factor $\alpha =$ 2000 [27], and obtain its time-domain graph and frequency domain graph as follows.

Fig. 3Time domain diagrams of VMD decomposition for different faults

Fig. 4Spectrum diagram of VMD decomposition for different faults

Fig. 3 reveals that the initial input signals overlap significantly prior to VMD decomposition, making it difficult to distinguish between their fault types. After performing VMD decomposition, partial separation of fault signal information becomes evident in IMF4. The accompanying time-domain plot and spectrogram further confirm the extraction of some time-frequency domain features. However, the figure also clearly shows that a substantial number of signals remain overlapping and difficult to identify, even after decomposition. This suggests that relying solely on VMD decomposition is insufficient for completely distinguishing between these five fault types under identical operating conditions.

Moreover, analyzing the same type of fault (outer race fault) across different operating conditions further highlights this limitation. For example, comparing outer race fault signals under Condition 1 (20 Hz (1200 RPM) speed with 0 V (0 N/m) load) versus Condition 2 (30 Hz (1800 RPM) speed with 2 V (7.32 N/m) load) demonstrates that using VMD decomposition alone does not yield more effective fault information.

Specific results are shown in Figs. 5 to 6.

Therefore, after decomposing the VMD signal, further use the fault information obtained from MRE. However, due to the fact that the values of the embedding dimension $m$ and tolerance parameter $r$ in the MRE index will affect the final effectiveness of obtaining fault information, optimization algorithms are used to optimize the two parameters $m$, $r$ of MRE before extracting fault information. Among them, there are 20 populations and 15 optimization iterations.

Fig. 5Time domain diagram after VMD decomposition of the same fault

Fig. 6Spectrum diagram of VMD decomposition for the same fault

Select a single fault signal with a length of 1024, and use two optimization algorithms, NGO and HHO, to conduct experiments with the objective function of minimum skewness squared and minimum margin factor, where $m\in$ [1, 5] is an integer and $r\in$ [0.01, 1] is a random number. The iteration curve of the experiment is shown below.

From the comparison of Fig. 7(a, c), and (b, d), it can be seen that the convergence speed of the NGO algorithm under the same objective function is significantly faster than that of the HHO algorithm.

Select 100 fault signals with a length of 1024 as experimental data, and use MRE to extract fault characteristic signals based on the parameters obtained through optimization. The final fault characteristic signal dataset size is 100×128. Divide the feature signal dataset into training and testing sets, with sizes of 50×128, and conduct experiments.

Take the average entropy value of the original data and select the data from the first 50 sample points of one type of fault for analysis. Fig. 8 reveals two key findings regarding the performance of the HHO and NGO optimization algorithms. First, their performance on the profit rate indicators is strikingly similar, with the curves completely overlapping. This convergence likely arises from shared characteristics: both algorithms are based on swarm intelligence principles, optimize similar objectives (minimizing classification error, maximizing feature accuracy), and were tested under identical experimental conditions (population size, iterations) using the same dataset. The limited complexity of the dataset may have also restricted the diversity of potential solutions, causing both algorithms to converge on the same optimal point.

Fig. 7Iteration curve

a) HHO minimum skewness square

b) HHO minimum margin factor

c) NGO minimum skewness square

d) NGO minimum margin factor

Fig. 8Entropy values of different parameters

Second, the analysis of average entropy values for fault data underscores the critical role of the chosen objective function. When minimizing the squared skewness, entropy values were excessively high, indicating poor discrimination of fault information, which resulted in suboptimal parameter performance and significant entropy differences. In contrast, optimizing the minimum margin factor produced a much higher overlap in entropy values, with values that were more appropriate for fault feature extraction. This suggests superior parameter accuracy, leading to more effective dimensionality reduction and improved classification accuracy.

Therefore, although the similar experimental conditions led both HHO and NGO to achieve identical profit rate outcomes in this specific test, the entropy analysis clearly shows that the effectiveness of parameter optimization is fundamentally dependent on the choice of objective function. The minimum margin factor was found to be significantly more effective than the squared skewness for generating robust and informative parameters from the data.

By using a two-dimensional contour map, further examine the distribution of mean entropy values for 100 sample points of 10 types of faults.

Fig. 9Minimum skew square contour map

a) HHO minimum skewness square

b) NGO minimum skewness square

Fig. 10Minimum margin factor contour map

a) HHO minimum margin factor

b) NGO minimum margin factor

It can be clearly seen from the graph that with the minimum skewness square as the objective function, the entropy values obtained by different optimization algorithms vary greatly, indicating that their performance as the objective function for optimization is not ideal; On the contrary, after extracting entropy information from the MRE parameters obtained with the minimum margin factor as the objective function, it can be seen that the fault characteristic signals of different optimization algorithms have minimal changes and almost no significant changes; At the same time, it can be found that the contour lines of the same type of horizontal faults obtained with the minimum skewness as the objective function are more encrypted, indicating that the entropy distribution obtained is wide, the data is discrete, and good fault information is not obtained; The entropy value obtained by the minimum margin factor is significantly more concentrated in the same fault category, indicating that the parameter effect obtained is better.

Further analyze the obtained 10 fault characteristic signals pairwise. Taking the results obtained by the NGO optimization algorithm as an example, compare the difference between the minimum skewness square and the minimum margin factor as objective functions through the distribution of Pearson correlation coefficients. Perform pairwise correlation tests on 10 types of fault characteristic signals, take the absolute value of the obtained correlation data, and classify it according to 0-0.05, 0.05-0.1, 0.1-0.15, 0.15-0.2, ≥ 0.2 to obtain Figs. 11, 12.

The dark blue area in the figure represents the proportion of the minimum correlation coefficient. It is evident from the figure that the MRE parameter obtained with the minimum margin factor as the objective function has a significantly larger proportion of the minimum correlation coefficient between the extracted fault feature signals than the minimum skewness square. The proportion of the dark blue area is as high as 44 %, indicating a high degree of discrimination in fault information.

There were ten types of original fault data, and the entire dataset size was 1000×128. The data was too large and the fault information was not clear enough, which not only was not conducive to subsequent classification, but also consumed too much time. Therefore, all fault category feature data with a size of 100×128 were subjected to dimensionality reduction using the ISOMAP manifold dimensionality reduction method [26], resulting in a data size of 100×32 after dimensionality reduction.

Fig. 11NGO skewness

Fig. 12NGO margin factor

Taking the NGO optimization algorithm as an example, further compare the differences in fault signals extracted from MRE parameters obtained by different objective functions after dimensionality reduction. Select the feature data of 1000×32 of the fault feature signal set after dimensionality reduction for visualization processing, and obtain Figs. 13, 14.

Fig. 13NGO minimum margin factor scatter plot

Fig. 14NGO minimum skewness squared scatter plot

Through the visualization scatter plot of the data, it can be clearly seen that the discrimination of all feature data with the minimum margin factor is significantly better than the minimum skewness square. The fault feature signals of the same category are more concentrated, while the discrimination between different categories is more obvious and the overlap is low, which is beneficial for subsequent fault identification.