Next Article in Journal
Spatial–Temporal Dynamic Graph Differential Equation Network for Traffic Flow Forecasting
Previous Article in Journal
Partial Eigenstructure Assignment for Linear Time-Invariant Systems via Dynamic Compensator
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Research on SVM-Based Bearing Fault Diagnosis Modeling and Multiple Swarm Genetic Algorithm Parameter Identification Method

1
Department of Automation, Guangdong University of Petrochemical Technology, Maoming 525000, China
2
Department of Computing, The Hong Kong Polytechnic University, Hong Kong 999077, China
3
School of Mechanical Engineering, Dongguan University of Technology, Dongguan 523000, China
4
The Hunan Provincial Key Laboratory of Health Maintenance for Mechanical Equipment, Hunan University of Science and Technology, Xiangtan 411201, China
*
Authors to whom correspondence should be addressed.
Mathematics 2023, 11(13), 2864; https://doi.org/10.3390/math11132864
Submission received: 31 May 2023 / Revised: 21 June 2023 / Accepted: 25 June 2023 / Published: 26 June 2023

Abstract

:
The bearing fault diagnosis of petrochemical rotating machinery faces the problems of large data volume, weak fault feature signal strength and susceptibility to noise interference. To solve these problems, current research presents a combined ICEEMDAN-wavelet threshold joint noise reduction, mutual dimensionless metrics and MPGA-SVM approach for rotating machinery bearing fault diagnosis. Firstly, we propose an improved joint noise-reduction method of an Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (ICEEMDAN) and wavelet thresholding. Moreover, the noise-reduced data are processed by mutual dimensionless processing to construct a mutual dimensionless index sensitive to bearing fault features and complete the fault feature extraction of the bearing signals. Furthermore, we design experiments on faulty bearings of multistage centrifugal fans in petrochemical rotating machinery and processed the input data set according to ICEEMDAN-wavelet threshold joint noise reduction and mutual dimensionless indexes for later validation of the model and algorithm. Finally, a support vector machine model used to effectively identify the bearing failures, and a multi-population genetic algorithm, is studied to optimize the relevant parameters of the support vector machine. The powerful global parallel search capability of the multigroup genetic algorithm is used to search for the penalty factor c and kernel parameter r that affect the classification performance of the support vector machine. The global optimal solutions of c and r are found in a short time to construct a multigroup genetic algorithm-support vector machine bearing fault diagnosis and identification model. The proposed model is verified to have 95.3% accuracy for the bearing fault diagnosis, and the training time is 11.1608 s, while the traditional GA-SVM has only 89.875% accuracy and the training time is 17.4612 s. Meanwhile, to exclude the influence of experimental data on the specificity of our method, the experimental validation of the Western Reserve University bearing failure open-source dataset was added, and the results showed that the accuracy could reach 97.1% with a training time of 14.2735 s, thus proving that the method proposed in our paper can achieve good results in practical applications.

1. Introduction

Petrochemical machinery and equipment are complex, with high speed, high power and high explosion-proof requirements. In this research, we investigated rolling bearings, which are very prone to failure in petrochemical rotating machinery, aiming to improve the fault diagnosis accuracy and series performance indices and study the noise reduction pre-processing of bearing vibration signals, bearing fault feature extraction and bearing fault identification and classification. At present, there are a series of bottlenecks in the research of fault diagnosis for petrochemical rotating machinery bearings. First, the petrochemical environment is complex and changeable, the collected bearing original data are easily disturbed by serious noise and the real bearing fault signal will be annihilated by the high-intensity noise. Second, the intensity performance of bearing fault characteristic signals is weak, making it difficult to extract the information reflecting the bearing fault characteristics from them. Finally, based on the industry characteristics, the number of bearing fault data collection requirements is higher than other industries; from a big data perspective, a model trained on a large number of data samples will be more accurate, but high data volume will affect computer processing efficiency. Therefore, the establishment of an efficient and reliable fault diagnosis model is one of the bottleneck problems.

1.1. Current Status of Domestic and International Research

1.1.1. Status of Research on Signal Noise Reduction Methods

Due to the complexity of its working environment, rolling bearings are inevitably affected by a variety of factors such as friction, impact and damping during operation, which leads to a large amount of noise in the obtained vibration signal and a strong non-linearity and non-smoothness. For this reason, it is necessary to denoise the signal to improve the extraction of its features. He et al. [1] proposed an improved wavelet total variation (IWATV) denoising method, which uses the wavelet coefficients obtained from wavelet decomposition of a noisy signal for local variance analysis to estimate the noise variance and provide scientific evaluation metrics. Luo et al. [2] proposed a probabilistic principal component analysis (PPCA) method based on GOA method and a comprehensive rolling bearing fault diagnosis method: Adaptive Probabilistic Principal Component Analysis-Enhanced morphological difference filter (APPCA-EMDF). APPCA can make up for the lack of EMDF in denoising ability and solve the interference problem caused by strong background noise. Yuan et al. [3] proposed a method based on multiple fractal fractional trend fluctuation analysis (MF-DFA) and smooth intrinsic time scale decomposition (SITD). This SITD method decomposes the vibration signal into several appropriate rotational components, overcomes the noise effect, maintains the effective signal and improves the signal-to-noise ratio. Huang et al. [4] proposed a new wavelet shrinkage operator to eliminate noise in nonlinear time series, which introduces a generic threshold to classify the details of signal and noise, and the proposed operator was used to solve the problems such as discontinuity of hard operators and large deviation of soft operators, which successfully reduced the noise in the nonlinear time series and also improved the classification accuracy of bearing faults in a noisy environment. REN et al. [5] proposed an in-band noise-reduction method based on iterative SVD (ISVD). The optimal frequency band of the vibration signal was determined using the enhanced wavelet packet transform kurtosis map, the kurtosis of each node was calculated using the envelope spectrum of the wavelet packet coefficient reconstruction signal, the envelope of the reconstruction signal was calculated using the Hilbert transform and the ISVD method was applied to reduce the in-band noise. Nguyen et al. [6] proposed a robust method for monitoring rolling bearing conditions, which uses a new empirical modal decomposition (EMD)-based method to eliminate high levels of noise from acoustic emission (AE) signals and a discrete wavelet packet transform (DWPT)-based envelope analysis technique to effectively identify bearing fault states. Ma et al. [7] proposed a new rolling bearing fault diagnosis method based on the improved VMD adaptive wavelet threshold, decomposing the VMD to obtain the modal components and constructing a dual determination criterion of sample entropy and correlation coefficient to filter the components. Then the adaptive wavelet thresholding function is proposed and quadratic noise reduction is applied to the hybrid imf, which in turn reconstructs each component to achieve joint noise reduction. Hu et al. [8] proposed a new method of rolling bearing fault diagnosis based on variational modal decomposition (VMD) and wavelet threshold denoising optimized by a genetic algorithm. Through theoretical calculation, numerical simulation and application research, the effectiveness and superiority of the method were verified and the method has some application prospects in rotating machinery fault diagnosis. Lee et al. [9] proposed an effective, noise-resistant fault diagnosis model for rolling elements. Signal noise processing uses product function (PF) selection and wavelet packet decomposition (WPD) for noise reduction, and the signal noise-reduction step can effectively remove high-frequency noise and extract fault information hidden under the noise. Since wavelets have the properties of low entropy, de-correlation and multi-resolution, it makes them very superior in analyzing and de-noising non-stationary signals. The wavelet thresholding method is simple in principle and effective in denoising, which is a common method for data denoising nowadays and is widely used for practical fault diagnosis.

1.1.2. Status of Research on Fault Feature Extraction Methods

Feature extraction is the basis of rolling bearing fault diagnosis. YU et al. [10] proposed a new feature extraction method combining the K-means method and standard deviation to select the most sensitive features, and also proposed an improved feature dimensionality reduction method to achieve a low-dimensional representation of the high-dimensional feature space and improve the accuracy of rotating machinery bearing fault diagnosis. Ju et al. [11] proposed a feature extraction method for improving multiscale entropy (IMSE) diagnosis, which improves the recognition accuracy of IMSE-extracted features and overcomes the information leakage problem when calculating the similarity of mechanical systems. Xue et al. [12] proposed a feature extraction method based on hierarchical analysis: Hierarchical Dispersive Entropy (HDE), which extracts features with better feature clustering effect and higher recognition rate than other methods, and can more accurately characterize the operating condition of bearings. Kuncan et al. [13] used a new feature extraction method for bearing faults, i.e., the one-dimensional ternary pattern (1D-TP) feature extraction method, which is an efficient and robust feature extraction method with high robustness for fast and real-time feature extraction of bearing faults. Li et al. [14] proposed a feature extraction and evaluation method for rotating machinery fault diagnosis, based on the central limit theory, and gave an extraction method to extract statistical features using existing signal processing tools, which can significantly improve the performance of fault classification and ensure the accuracy of fault classification. He et al. [15] proposed a new repetitive transient extraction method with group sparse structure for detecting rolling bearing faults. The proposed repetitive transient extraction algorithm (RTEA) can successfully extract the repetitive transients caused by bearing faults, and the method can effectively extract the features of outer ring defects and inner ring defects. Zhu et al. [16] constructed a multi-weight fusion-based singular value feature extraction method based on the contribution of singular value and entropy weights to solve the problem of it being difficult to effectively extract feature vectors in rolling bearing fault diagnosis, and combined it with SVM classifier to propose a feature fusion-based bearing fault diagnosis method, which can effectively extract sensitive features of bearings and reduce the interference of non-sensitive features to the diagnosis model. Xia et al. [17] proposed a spectral regression (SR) based method to extract fault features from the original features such as time, frequency, and time–frequency domain features of bearing accelerometer sensor signals, which can reduce the computational cost and preserve more structural information about different bearing faults and severity. There are many methods of feature extraction and various evaluation indexes based on different methods. Kumar et al. [18] proposed a migration learning based method for bearing fault detection, using a pre-trained ResNetV2 model as the base model to establish an effective bearing fault detection strategy, which reduces the need for manual feature extraction and selection of bearing faults. Combined with the actual working conditions of petrochemical production, exploring the fault feature extraction method of rotating mechanical equipment bearings is the focus of current research.

1.1.3. Status of Research on Fault Classification and Identification

The traditional bearing condition identification mainly relies on human experience to judge, which makes it difficult to diagnose the operation condition of the bearing quickly and accurately, and now the emergence of automated and intelligent algorithms has brought a new way of thinking for fault diagnosis and condition identification. Toumi et al. [19] designed an embedded intelligent system to identify the type and severity of bearing faults using envelope analysis to process the raw vibration signals obtained from the CWRU dataset to obtain the peak amplitudes corresponding to the bearing fault frequencies used as input to the MLP network classifier. Yan et al. [20] proposed a fault classification algorithm based on multi-domain feature optimization SVM, which consists of three main stages: multi-domain feature extraction, feature selection and fault identification. The method solves the problem that the current fault classification framework generally focuses on the pattern of single-domain feature classifier, which easily leads to insufficient feature extraction and low recognition accuracy. Liu et al. [21] introduced the higher-order cumulants (HOC) that can quantitatively describe the nonlinear feature signals with close relationships between mechanical faults, achieve noise reduction of the original bearing vibration signals to obtain a bispectral estimation image, and proposed a rolling bearing fault identification method that combines the basic higher-order spectrum (HOS) theory and the fuzzy clustering method in the field of data mining. Wan et al. [22] proposed a method based on group decomposition (SWD), morphological envelope dispersion entropy (MEDE) and random forest (RF) based on the ideas of signal denoising, feature extraction and pattern classification to achieve effective detection and intelligent identification of weak faults in rolling bearings. Fatima et al. [23] proposed a general method for fault classification of bearing imbalance in rotor bearing systems, which is based on a compensated distance-evaluation technique for feature selection, support vector machine and grid search methods for optimization of parameters in the fault classification process, which is a robust SVM technique that uses only time domain data and does not require any additional preprocessing for the balanced faults for classification. JIANG et al. [24] proposed a new convolutional neural networks (CNNs)-based bearing fault classification method in order to eliminate noise interference and take into account the possible connection between signal frames. Using the sensitivity to spectral kurtosis (SK) pulses, a SK-based filtering method was proposed to suppress the noise, which has strong classification ability under the interference of plant noise and Gaussian noise. Tajeddini et al. [25] used vibration and current sensors as a multi-sensor framework for high-precision bearing fault detection, and used wavelet packet transform for each sensor as an effective pre-processing method for signal denoising, and established a generic threshold based on the Steyn unbiased risk estimation method in the denoising algorithm, and trained a support vector machine for fault classification after extracting appropriate time-domain features from the denoised signal. Gong et al. [26] proposed a bubble entropy-based rolling bearing fault diagnosis method without parameter discussion by extracting fault features based on the refined composite multiscale bubble entropy (VMD-RCMBE) of variational pattern decomposition, which improves the fault identification accuracy of bearings while reducing or even eliminating the influence of parameter discussion. Attoui et al. [27] proposed a new time-frequency method for feature extraction based on classical wavelet packet decomposition features by using WPD and STFT as a time-frequency domain fault feature extraction scheme, LDA, LSDA or feedback envelope method as a feature reduction technique, and ANFIS as an automatic classification system to develop the proposed procedure, which achieves accurate detection and classification of bearing faults with different bearing conditions, different onset times, fault locations, fault types and fault severity. Tan et al. [28] proposed an intelligent bearing fault identification method based on SMCMDE and equilibrium optimizer-based support vector machine (EO-SVM), which not only extracts sensitive features of bearing faults, but also has higher accuracy and stability than other methods. Lu et al. [29] investigated an effective and reliable deep learning method based on convolutional neural networks (CNN) using cognitive theory to introduce the advantages of image recognition and visual perception into bearing fault diagnosis by simulating the cognitive processes in the cerebral cortex, and the method improved the accuracy of failure mode classification of rolling bearings with respect to environmental noise and fluctuating operating conditions. Tong et al. [30] proposed a multi-view learning (DAML)-based bearing fault diagnosis method which obtains a multi-view data set involving vibration and acoustic signals by performing a fast Fourier transform (FFT), and then implements a multi-view feature (MVF) representation including view-invariant information and class discrimination information in a common subspace based on typical correlation analysis (CCA) and uncorrelated linear discriminant analysis (ULDA). The multi-view feature (MVF) representation including view invariant information and class discriminant information is then implemented in a common subspace based on typical correlation analysis (CCA) and uncorrelated linear discriminant analysis (ULDA), and finally a bearing fault is identified using a K-nearest neighbor (KNN) classifier based on the multi-view feature. The petrochemical production process has complex working conditions and various types of faults, and intelligent algorithms are uniquely effective in fault diagnosis, which has always been the focus of scholars’ attention. Yang et al. [31] proposed a sample entropy based on mutual information to select the effective intrinsic mode function (IMF), and with the help of minimum entropy deconvolution, decomposed the denoised signal into several IMF by ICEEMDAN method as a pre-filter to increase the kurtosis value by a factor of about 3.2. The results show that compared with ICEEDAN and MED, the proposed hybrid scheme has better performance in weak fault feature extraction under strong noise background. Zou et al. [32] based on improved adaptive noise fully integrated Empirical Mode decomposition (ICEEMDAN) and inverse residual convolutional neural networks integrated with Ghost modules. An intelligent bearing fault diagnosis method is proposed. In this method, ICEEMDAN is used to decompose the signal, and the kurtosis, correlation coefficient and margin factor index are proposed to filter and reconstruct the important modal components, and then multi-scale sample entropy is used to construct the feature vector, and finally Ghost module, convolution block attention module and inverse residual bottleneck are combined to construct Ghost-IRCNN and realize bearing fault classification and diagnosis. Zhao et al. [33] proposed multi-point optimal minimum entropy deconvolution adjustment (MOMEDA). Bearing simulators were used to collect vibration signals of bearing inner and outer rings, enhanced by MOMEDA, decomposed by ICIEDAN into several intrinsic mode functions (IMFs), and analyzed by envelope demodulation. Finally, the shaft speed frequency, BPFI (ball-pass frequency, inner ring) and harmonic frequency, sideband interval, BPFO (ball-pass frequency, outer ring) and harmonic are obtained. This method can be used to accurately extract different frequency components of bearing fault vibration signals and diagnose different bearing fault locations. Li et al. [34] proposed a combined noise-reduction method for cut acoustic signals based on improved adaptive noise fully integrated Empirical Mode decomposition (ICEEMDAN) and singular value decomposition (SVD). The proposed method can effectively eliminate the influence of background noise in the signal and retain the characteristic frequency of the original cutting sound signal. Compared with traditional noise-reduction methods, ICEEMDAN-SVD combined noise-reduction method performs better in the noise-reduction evaluation criteria of signal-to-noise ratio and RMS error. Qin et al. [35] proposed a rolling bearing troubleshooting method based on improved adaptive white noise fully integrated Empirical Mode decomposition (ICEEMDAN) and Sparrow Search algorithm (SSA) to optimize correlation vector machine (RVM). The ICEEMDAN method is used to decompose the original vibration signal, and the ratio of power entropy and intrinsic mode component (IMF) power to the total power of the signal is taken as the fault characteristic. The width coefficient and hyperparameters of RVM were optimized by the SSA algorithm, and the RVM model was established to realize fault identification of rolling bearings.
Based on the petrochemical site environment and the above-mentioned literature, this paper proposes a rotating machinery bearing fault diagnosis method based on the combination of ICEEMDAN-wavelet threshold joint noise reduction, mutual dimensionless index and MPGA-SVM and designs the relevant experiments to obtain the raw data with the faulty bearings of multi-stage centrifugal fans of petrochemical rotating machinery. The raw data are processed by ICEEMDAN-wavelet threshold joint noise reduction and mutual dimensionless index feature extraction to obtain the standard input data set (including training set and test set), which is the key of the subsequent model and model parameter identification algorithm. Finally, based on the standard data set, different models and algorithms are applied to carry out the fault classification of bearings, so as to verify the superior performance of the MPGA-SVM model proposed in this paper.

1.1.4. Arrangement of the Content of This Paper

The rest of this paper is organized as follows. The second part introduces the research framework of this paper, which focuses on three aspects of wavelet-based threshold denoising, mutual dimensionless index extraction features and the MPGA-SVM-based bearing fault diagnosis method, and discusses the research content, technical routes and innovation points. The third part conducts experiments on a petrochemical multistage centrifugal fan fault diagnosis platform to obtain bearing fault data, verify the effectiveness of the designed denoising, feature extraction and fault diagnosis methods and compare them with traditional algorithms to prove the superiority of the research method. The fourth part discusses the main work and contributions of this paper.

2. Materials and Methods

2.1. Proposed Fault Diagnostic Framework

In this paper, we propose a fault diagnosis method for bearings based on ICEEMDAN-Wavelet threshold denoising and mutual dimensionless indicators and MPGA-SVM.
  • Part 1: A joint ICEEMDAN-wavelet thresholding based noise-reduction method. In this paper, a signal noise-reduction method based on the joint ICEEMDAN decomposition algorithm and wavelet thresholding is proposed to establish a bearing vibration signal noise-reduction model. After the bearing signal collected in the field has been decomposed by ICEEMDAN, the IMF components arranged from high frequency to low frequency are obtained; the method only denoises the high-frequency IMF components that contain more noise by wavelet threshold and retains the low-frequency IMF components that contain less noise; finally, the high- and low-frequency IMF components are reconstructed.
  • Part 2: Extraction of fault features using a mutual dimensionless index. The vibration signal of the rolling bearing of rotating machinery implies a large amount of bearing fault state information. After the noise reduction pre-processing, the bearing signal still cannot intuitively determine the specific operating state of the bearing, and the signal must be extracted by fault features. Mutual dimensionless indicators are widely used in this field because they are sensitive to bearing failures and can accurately reflect information about the characteristics of shaft failures.
  • Part 3: Construction of a fault identification model for rotating machine bearings based on the MPGA-SVM algorithm. The parallelism of multiple swarm genetic algorithms and the powerful global search capability are used to solve the optimal values of key parameters c and r, which affect the classification performance of the SVM, to build the SVM optimal rotating machine bearing fault classification model.
The overall research framework of the paper is shown in Figure 1.

2.2. Related Theories and Technology Roadmap

2.2.1. Wavelet Threshold Denoising and ICEEMDAN Decomposition

Bearing vibration signals are characterized by nonlinearity, multi-source heterogeneity and randomness, and there is serious noise interference. The common methods for noise reduction of the original bearing signal are time-domain filtering noise reduction, wavelet threshold noise reduction, empirical modal decomposition noise reduction, independent component analysis noise reduction and neural network method noise reduction. Each method has certain advantages and disadvantages, and there are many related studies, among which wavelet analysis is more suitable for the processing of local non-stationary signals. The advantages are, firstly, the wavelet transform achieves the improvement of local resolution in the frequency domain, and it is relatively easy to analyze the transients in the asymptotic signal and noise. Secondly, the wavelet transform can compress the signal into a local region in time and frequency, so that each wavelet term obtained after wavelet decomposition of a non-stationary signal can correspond to a local part of the analyzed signal, and the temporal changes of the signal can be well captured, which is very useful for analyzing the local characteristics of the signal. Thirdly, the wavelet decomposition algorithm outperforms other classical analysis methods in terms of efficiency and speed of time and frequency analysis.
Wavelet threshold denoising is a very widely used signal noise-reduction processing method. First, the original signal is decomposed into a series of wavelet coefficients. Second, a threshold is applied to each wavelet obtained after the decomposition; when the wavelet coefficients are below the threshold, they are considered noise to be eliminated. Finally, the wavelet coefficients after wavelet threshold processing is completed are reconstructed to obtain the noise-reduced signal. The principle of the method is shown in Figure 2.
Research shows that, although the wavelet threshold denoising method is effective in some aspects of noise processing, the method ignores that the decomposed wavelet coefficients with smaller thresholds also have valid information. Although the signal is weak, it cannot be ignored. If these weak signals are ignored, the denoising pre-processing will cause the loss of useful information in the weak energy.
Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN): Adaptive white noise is added several times in the decomposition stage. The decomposition process is averaged by the ensemble to have better completeness. Then, the residual components are used to reconstruct a more accurate intrinsic modal component (intrinsic mode function, IMF). The previous literature proposes an Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (ICEEMDAN), which solves the intrinsic modal component in a different manner from CEEMDAN. ICEEMDAN is based on the residual signal obtained from the previous iteration minus the residual signal of the current iteration and subsequently averaged together. The final decomposition of the eigenmode component shows a high frequency to low frequency ordering. In addition, the appearance of spurious modal components is effectively suppressed.
Based on the above theoretical basis, this paper proposes an ICEEMDAN decomposition algorithm combined with wavelet thresholding as a signal noise-reduction method to establish a bearing vibration signal noise-reduction model. After ICEEMDAN has decomposed the bearing signal collected in the field, the IMF components arranged from high frequency to low frequency are obtained, and the method only denoises the high-frequency IMF components, which contain more noise by wavelet thresholding, while retaining the low-frequency IMF components, which contain less noise. Finally, the high- and low-frequency IMF components are reconstructed. Figure 3 shows the process of the ICEEMDAN-wavelet threshold joint noise reduction.
The original signal with nonlinearity, which contains more noise, is decomposed by ICEEMDAN to obtain a series of IMF components from high frequency to low frequency.
The correlation coefficient method is used to calculate the correlation coefficients of the IMF components obtained from the ICEEMDAN decomposition. The high-frequency IMF components, which contain more noise (correlation coefficient value greater than 0.1), are selected.
The correlation coefficient enables one to measure the correlation between the IMF components after the ICEEMDAN decomposition and the original signal of the bearing after adding noise. A larger correlation coefficient value proves that the two are more similar and the IMF component contains more noise. On the contrary, a smaller correlation coefficient value proves that the two are less similar and the IMF component contains less noise. The correlation coefficient value is used as the basis to select the high-frequency components with more noise in the IMF components obtained by the ICEEMDAN decomposition. The Pearson correlation coefficient is used, and the correlation is defined as follows:
R ( I M F i , f ) = ( c o v ( I M F i , f ) σ I M F i σ f )
In Equation (1),
c o v ( I M F i , f ) = E [ ( I M F i E ( I M F i ) ) ( f E ( f ) ) ] σ I M F i = 1 N i = 1 N ( I M F i I M F i ¯ ) 2 σ f = 1 N i = 1 N ( f f ¯ ) 2
where c o v ( ) represents the covariance, f is the signal after white noise has been added to the bearing signal, f ¯ is the mean value of the signal, and σ I M F i and σ f are the standard deviations of I M F i and f respectively.
The IMF components with high frequency and more noise are wavelet-decomposed, and the appropriate decomposition scale and wavelet function should be selected before the decomposition.
A suitable threshold function is selected; then, a threshold is applied for the wavelet coefficients obtained after the wavelet decomposition for quantization.
The soft threshold value is shown in Equation (3):
I M F i ˜ = s g n ( I M F i ) ( I M F i λ ) I M F i λ 0 I M F i < λ
The hard threshold is shown in Equation (4):
I M F i ˜ = I M F i ( I M F i λ ) I M F i λ 0 I M F i < λ
where s g n ( ) is the symbolic function; | I M F i | is the unprocessed wavelet coefficients; λ is the threshold; I M F i ˜ is the wavelet coefficients after the threshold quantization process. Threshold λ is defined as:
λ = σ i 2 ln N = τ i 0.6745 2 ln N
In Equation (5), σ i is the standard deviation of the noise component layer i; 0.6745 is the empirical value; τ i is the absolute median of the wavelet coefficients in layer i. The noise-reduction effect is more satisfactory when the hard-threshold signal is used.
The wavelet coefficients after wavelet thresholding are analyzed. The coefficients in the wavelet coefficients that are larger than the wavelet threshold are selected. The wavelet reconstruction function waverec is used to reconstruct the signal. The high-frequency IMF component after noise reduction is obtained.
The noise-reduced high- and low-frequency IMF components are summed and reconstructed. Then, the noise-reduced signal x ( t ) ˙ is obtained.

2.2.2. Mutual Dimensionless Indicators

Let the collected rotating machinery bearing fault observable vibration signal be z ( t ) , which is composed of the fault-free vibration signal s ( t ) , bearing fault characteristic signal x ( t ) and Gaussian white noise v ( t ) . The fault-free vibration signal s ( t ) is the fault-free vibration signal collected by the sensor at the early stage of operation of the new rotating machinery after break-in. Figure 4 shows the combined model of the observable vibration signal z ( t ) .
Figure 4, z ( t ) [36] is expressed as
z ( t ) = c [ s ( t τ ) ] + P = 1 P x p ( t ) + v ( t )
In Equation (6), τ is the delay time constant, c is the correlation coefficient and P = 1 P x p ( t ) represents a series of fault characteristic signals in the bearing, such that x ( t ) = P = 1 P x p ( t ) is transformed to:
z ( t ) = c [ s ( t τ ) ] + x ( t ) + v ( t )
To construct the mutual dimensionless index, the fault features are first extracted from the bearing vibration signal. Figure 5 shows the specific process.
Fault feature extraction takes the observable vibration signal z ( t ) and fault-free vibration signal s ( t ) as the starting point and performs fast Fourier transform (FFT) on them to obtain Z ( k ) and S ( k ) in the frequency domain. Then, the method takes the complex conjugate S ( k ) of S ( k ) * . The correlation function M ( k ) in the frequency domain is the product of Z ( k ) and S ( k ) * . Then, the method performs a fast Fourier inverse transform (IFFT) on M ( k ) , and their correlation function m ( t ) in the time domain can be obtained. Afterward, the correlation coefficient c is calculated by correlation analysis of the delay time τ . Finally, y ( t ) is obtained by linearly subtracting z ( t ) and s ( t ) to complete the fault feature extraction. Where p ( y ) is the amplitude probability density of fault characteristic signal y ( t ) and p ( s ) is the amplitude probability density of non-fault vibration signal s ( t ) . The definition of the mutual dimensionless index is shown in Equation (8).
ϑ y s = [ + y l p ( y ) d y ] 1 l [ + s n p ( s ) d s ] 1 n = E ( y l ) l E ( s n ) n
Mutual waveform indicator: let l = 2 and n = 1
S y s = [ + y 2 p ( y ) d y ] 1 2 [ + s p ( s ) d s ] = E ( y 2 ) E ( s )
Interpulse indicators: let l = and n = 1
T y s = lim l [ + y l p ( y ) d y ] 1 l [ + s p ( s ) d s ] = lim l E ( y l ) l E ( s )
Mutual margin index: let l = and n = 1 2
M y s = lim l [ + y l p ( y ) d y ] 1 l [ + s 1 2 p ( s ) d s ] 2 = lim l E ( y l ) l [ E ( s ) ] 2
Mutual peak metrics: let l = and n = 2
C y s = lim l [ + y l p ( y ) d y ] 1 l [ + s 2 p ( s ) d s ] 1 2 = lim l E ( y l ) l E ( s 2 )
Mutual peak metrics: let l = 4 and n = 2
ϑ y s = + y 4 p ( y ) d y [ + s 2 p ( s ) d s ] 2 = E ( y 4 ) [ E ( s 2 ) ] 2
After the ICEEMDAN-wavelet threshold noise reduction and mutual dimensionless index processing, the mutual dimensionless indices of the bearing signals under different fault states are obtained.

2.2.3. MPGA-SVM-Based Rotating Machinery Bearing Fault Model Parameter Identification

In solving the classification problem of linear data samples, SVMs establish classification surfaces at relatively distant locations from the two types of samples. In solving the classification problem of nonlinear data samples, the SVM transforms the low-dimensional nonlinear classification problem into a high-dimensional spatially divisible one by performing a high-dimensional spatial transformation through the kernel function and effectively classifies linearly indivisible samples. In other words, it aims to find a hyperplane to partition the data samples with a partitioning principle of interval maximization. Figure 6 shows the SVM classification surface.
The classification performance of SVMs is mainly determined by the penalty factor c and parameters r (where r = 1 / 2 σ 2 ) of the kernel function. Currently, there is no specific approach to select c and r. Grid search and cross-test are the traditional approaches, but they are inefficient. If the range and step size of the parameters to be optimized are not properly set, the selection of optimal parameters will be affected. Thus, this paper will introduce the multiple swarm genetic algorithm to solve the optimal values of key parameters c and r, which affect the classification performance of SVMs, by using the parallelism of the multiple swarm genetic algorithm and powerful global search ability to establish the optimal classification model of SVMs.
The multiple population genetic algorithm performs the simultaneous optimization search for multiple populations, and the parameters are set differently for different subpopulations to realize the diversity of optimization results.
Meanwhile, each subpopulation can be linked using a migration operator to realize the co-evolution of multiple subpopulations, which yields the optimal solution. The optimal solution in the evolutionary generation of each subpopulation is preserved by the manual selection operator, which is also the criterion for determining whether the algorithm converges. In this paper, the minimum number of generations of optimal individuals is maintained as the basis for the end of algorithm operation. The principle of the multi-population genetic algorithm is shown in Figure 7.
The steps of the MPGA-optimized SVM are as follows:
  • Initialization: Multiple population genetic-algorithm-related parameters are initialized to determine the number of subpopulations and number of individuals in the subpopulation, while the SVM parameters c and r are binary encoded.
  • Determine the fitness function: The fitness function depends on the specific problem to be solved.
  • Selection: The selection operation is performed using the proportional selection method, and the individuals with larger fitness values in the current population will be replicated to the next generation.
  • Crossover and variation: The crossover probability PC is usually set in the range of 0.7–0.9, and single-point crossover is used. The variation is used to change a code position of a chromosome individual with a relatively small variation probability PM. The variation probability PM is usually set in the range of 0.001–0.1, and discrete variation is usually used.
  • Migration and manual selection: The migration operator is set to control the linkage between subpopulations, introduce the optimal individuals into the next subpopulation, make the linkage between subpopulations, and find the optimal result by co-evolution. The elite population is formed by manually selecting the best individuals in each generation to protect the optimal solution. The minimum number of generations of elite individuals is used as the termination condition of the algorithm. Figure 8 shows the flow chart of the MPGA-SVM algorithm.
Multiple swarm genetic algorithms can effectively search for the optimal values of key parameters c and r, which affect the classification performance of SVMs to obtain the optimization model of SVM. Hence, the methods mentioned in Section 2 are fused. Then, a bearing fault diagnosis method based on the combination of the ICEEMDAN-wavelet threshold, mutual dimensionless index and MPGA-SVM is proposed and experimentally validated using a dataset of petrochemical rotating machinery bearings under different fault states.

3. Experiments and Result Analysis

3.1. Experimental Environment

The petrochemical multi-stage centrifugal fan fault diagnosis platform is mainly composed of four parts: the inverter motor, load controller, gearbox and multi-stage centrifugal fan.
This paper focuses on the bearings of multi-stage centrifugal fan parts. The power supply adopts a three-phase alternating current and is connected to the inverter through a three-phase gate switch. When the three-phase gate switch is closed, the inverter motor is energized, runs and is connected to the fan through the coupling, load controller, gearbox and fan; thus, it drives the entire multi-stage centrifugal fan fault diagnosis unit to run. The experimental platform can simulate various multi-stage centrifugal fan rolling bearing fault states such as bearing inner ring wear, bearing outer ring wear and bearing’s lack of the rolling body. Figure 9 shows the petrochemical multi-stage centrifugal fan fault diagnosis platform.
The multi-stage centrifugal fan of the petrochemical unit uses the bearing type HRB-1310ATN, which is a class P5 grade. The inner diameter length is 50 mm, the outer diameter length is 110 mm, the bearing height is 27 mm, the bearing has 13 rolling elements, and the rolling element diameter is 14 mm.
In this paper, for the multi-stage centrifugal fan rolling bearing, three common types of failure are conducted for research and analysis: bearing inner ring wear, bearing outer ring wear and bearing’s lack of the rolling body. Figure 10 shows the petrochemical rotating machinery multi-stage centrifugal fan bearing failure parts.
The detailed experimental procedure is as follows:
  • In the first step, the sensor probe with a strong magnetic suction tip is vertically adsorbed above the bearing housing of the multistage centrifugal fan.
  • The second step is to set the sampling frequency of the EMT390 data collector to 1024 Hz and collect a set of data every 8 s, and then set the data to be saved in automatic mode.
  • The third step is to control the inverter motor to run at a constant speed of 1200 r/min and maintain a constant sampling frequency and number of sampling points for sampling.
  • The fourth step in the fan bearing seat position is to replace different fault parts respectively for the experiment.
  • Multi-stage centrifugal fan rolling bearing data collection is completed, shut down and the data are automatically stored.
In order to study the collected multi-stage centrifugal fan bearing vibration signal characteristics, respectively, from the time domain analysis and envelope spectrum analysis of two perspectives, we considered the normal bearing, inner ring wear, outer ring wear and bearing’s lack of rolling body—these being the types of bearing vibration data to analyze and study. The time domain waveforms and envelope spectra of the inner ring wear vibration signals show that the time domain waveforms are more complex when the inner ring wear fault occurs, and the envelope spectra can clearly observe that the fault characteristic frequencies are 71.77 Hz, 144.04 Hz at the second frequency, 215.8 Hz at the third frequency and 288.08 Hz at the fourth frequency. When the outer ring wear fault occurs, the appearance of the shock signal in the time domain waveform shows a certain periodicity, and the characteristic frequencies of the outer ring wear fault are 81.54 Hz, two-fold frequency 163.08 Hz and three-fold frequency 244.62 Hz. When the bearing is missing its rolling element, the appearance of the shock signal in the time domain waveform shows a greater randomness, and the characteristic frequencies of the bearing-missing-rolling-element fault are 73.24 Hz, two-fold frequency 146.48 Hz and three-fold frequency 146.48 Hz. The characteristic frequencies of the bearing-missing-rolling-element fault were 73.24 Hz, 146.48 Hz, 219.23 Hz and 292.48 Hz.

3.2. Experimental Verification of the Rolling Bearing Fault Diagnosis Method

The parameters related to the ICEEMDAN method are set as follows: the standard deviation of the added Gaussian white noise Nstd = 0.2, Number of times noise is added NR = 100, Maximum number of iterations MaxIter = 200. The bearing signals acquired in the experiment were decomposed by the ICEEMDAN method, and the results are shown in Figure 11, Figure 12, Figure 13 and Figure 14.
The signals of different bearing fault types are decomposed by ICEEMDAN to obtain a series of IMF components from high frequency to low frequency, as shown in Table 1. The correlation coefficient of each IMF component is calculated by the Pearson correlation coefficient formula, where the IMF components with correlation coefficient values greater than 0.1 are identified as high-frequency components with more noise and are subjected to wavelet threshold noise reduction, while IMF components with correlation coefficient values less than 0.1 contain less noise and are selected for retention.
Among the IMF components of the normal bearing signal obtained after ICEEMDAN decomposition, the phase relation values corresponding to the modal components IMF1–IMF9 are all greater than 0.1, indicating that these components contain a lot of noise and need to be denoised by the wavelet threshold. The phase relation values corresponding to IMF10–IMF14 are all less than 0.1, which indicates that these components contain less noise and are reserved.
In the next step, the wavelet threshold denoising method is used to denoise the selected high-frequency noise-containing components. The wavelet threshold denoising selects the db6 wavelet basis function, hard threshold function, heursure threshold and the decomposition layer number is 5. Firstly, the wavelet decomposition of high-frequency IMF components with noise is carried out, and then the threshold quantization of the wavelet coefficients obtained after decomposition is carried out. The wavelet coefficients after the wavelet threshold processing are analyzed and the wavelet coefficients greater than the wavelet threshold are selected. The wavelet coefficients of each high-frequency IMF component after the wavelet threshold quantization are reconstructed by using the wavelet reconstruction function instruction waverec. At this time, the noise reduction of high-frequency IMF weight has been completed. Then, all the high frequency IMF components and low frequency IMF components after noise reduction are added and reconstructed to complete the combined noise reduction of bearing vibration signals by ICEEMDAN-wavelet threshold.
In order to more fully reflect the difference between the bearing signal before and after noise reduction, two indexes, continuous mean square error and signal-to-noise ratio, are introduced to evaluate the signal after noise reduction.
1.
Root Mean Squared Error (RMSE): In the signal noise-reduction metric, the root mean square error is defined as the expected value of the squared difference between the un-denoised signal and the denoised signal recalculated as shown in Equation (14).
R M S E = 1 N i = 1 N ( x i x i ) 2 ( i = 1 , 2 , . . . , n )
2.
Signal-to-noise ratio (SNR): The signal-to-noise ratio is defined as shown in Equation (15).
S N R = 10 lg i = 1 N x i 2 i = 1 N ( x i x i ) 2 ( dB )
At the same time, in order to fully demonstrate the superiority of the ICEEMDAN-wavelet threshold noise-reduction method, the conventional wavelet threshold noise-reduction method is used to compare with it, and the results of the two noise-reduction methods are evaluated by the root mean square error and the signal-to-noise ratio. Among them, the wavelet thresholding noise reduction still uses the db6 wavelet basis function, heursure thresholding, hard thresholding function, and the number of decomposition layers is set to 5.
Based on the comparison results in Table 2 and the evaluation criteria of relevant noise-reduction indicators, it is clear that the ICEEMDAN-wavelet threshold combined denoising method proposed in this paper has a better effect on the noise reduction of bearing signals than that of the wavelet threshold combined denoising method, and the root-mean-square error of the ICEEMDAN-wavelet threshold combined denoising method is generally lower than that of the wavelet threshold reduction method. The signal-to-noise ratio is obviously higher in the former than in the latter.
The collected bearing data were subjected to ICEEMDAN-wavelet threshold noise reduction and subsequent mutual dimensionless processing to obtain a dataset consisting of five mutual dimensionless indicators, which are sensitive to bearing fault characteristics, as shown in Table 3. Then, this dataset is used as the input of the MPGA-SVM model. Then, the fault identification experiments are performed to verify the effectiveness of the fault diagnosis model.
To fully demonstrate the superiority of the MPGA-SVM algorithm proposed in this paper for bearing fault diagnosis in petrochemical rotating machinery, the SVM, GA-SVM and MPGA-SVM methods are simultaneously compared and experimented with. The experimental data are selected from the total data set of normal bearing, bearing outer ring wear, bearing inner ring wear and bearing missing ball. Eighty groups of samples are randomly selected for normal bearing, bearing outer ring wear, bearing inner ring wear and bearing missing ball, so there are 320 groups of samples in total. Table 4 shows the correspondence of the fault types, numbers and categories.
Conventional support vector machine: Among the data samples, 50% (160 groups) are selected as the training set for training the SVM model, and the remaining 50% (160 groups) are used as the test set for the SVM model. Parameter penalty factors and parameters of the kernel function in the SVM are randomly generated to build the SVM bearing fault identification classification model. The dataset was used as the input to the SVM model for several experimental validations, where one bearing fault identification result is shown in Table 5.
From Table 6, the results of six different experiments show that the traditional SVM model is not ideal for bearing fault identification, and the overall fault diagnosis accuracy is low. The randomly generated SVM parameters penalty factor c and kernel function parameters r are different each time, which leads to different fault diagnosis results in each experiment and longer algorithm running time. The analysis leads to the conclusion that the appropriateness of the values of c and r greatly impacts the classification performance of the SVM model. This result also shows the importance of incorporating an optimization algorithm to optimize the two SVM parameters c and r.
Genetic algorithm optimization support vector machine (GA-SVM): The SVM parameter penalty factor c and r kernel parameter are randomly generated, which leads to the poor classification effect of the SVM model for bearing fault identification. Hence, the standard genetic algorithm (SGA) is introduced to optimize the SVM parameter penalty factor c and kernel function parameter r to establish the GA-SVM bearing fault diagnosis model. The relevant parameters of the genetic algorithm are set as shown in Table 7.
Among the data samples, 50% (160 groups) are selected as the training set to train the SVM model. The optimized GA-SVM model is obtained using GA for the SVM parameter penalty factor c and kernel parameter r. The remaining 50% (160 groups) are used as the test set for the GA-SVM model. The GA-SVM model performance is validated. Several experimental validations are conducted, and one of the bearing fault identification classification results is shown in Table 8.
The optimal fitness value is 89.875%, at which time the optimal solutions of c and r are 25.2113 and 0.0425, respectively, although GA-SVM is more effective than SVM Although GA-SVM has improved effect compared to SVM, the overall fault identification accuracy remains low mainly because the genetic algorithm falls into “premature convergence” in the process of finding the best. The prediction results of each fault are shown in Table 9.
Multiple population genetic algorithm optimized support vector machine (MPGA-SVM): The next step is to validate the MPGA-SVM model proposed in this paper, where The parameter settings associated with the multi-population genetic algorithm are shown in Table 10.
  • In the first step, the parameters of the multi-population genetic algorithm are set: the maximum number of genetic generations M = 20, the number of individuals NING = 160, the crossover probability PC = 0.7 + (0.9−0.7) × rand(MP,1), set in the range [0.7, 0.9], the variation probability PM = 0.001 + (0.05−0.001) × rand(MP,1), set in the range [0.001, 0.05] setting the termination condition of the algorithm to the minimum number of maintained generations MAXGEN = 20.
  • The second step binary encodes the penalty factor c of the support vector machine and the kernel function-related parameter r with the number of bits set to 20 as the initial population and calculates the individual fitness in each population separately, where c,r ϵ [0, 100].
  • The third step trains the SVM on the randomly generated training set, with the objective function defined as the average accuracy obtained by cross-validation.
  • The fourth step performs individual subpopulation-independent crossover, mutation, selection and other related genetic operations on different subpopulations.
  • The fifth step calculates the individual fitness of each new population and generates linkage by migration between subpopulations through migration operator for the better individuals.
  • The sixth step establish elite populations by manually selecting the best individuals in each population.
  • The seventh step determines whether the multiple population genetic algorithm termination condition is satisfied, and if it is satisfied, the optimal result is output, and if it is not satisfied, it returns to the fourth step and continues the relevant genetic algorithm operation process.
  • Finally, the optimal combination of the penalty factor c and the kernel function parameter r obtained by the multi-population genetic algorithm is imported into the support vector machine to obtain a superior bearing fault identification classification model, which is then validated with a test sample set.
Among the data samples, 50% (160 groups) are selected as the training set to train the SVM model, and the remaining 50% (160 groups) are used as the test set of the MPGA-SVM model. The optimized MPGA-SVM model is obtained by seeking the SVM parameter penalty factor c and kernel parameter r through the parallelism of MPGA and the powerful global search capability. The performance of the MPGA-SVM model is verified. Figure 15 shows the MPGA-SVM evolution process. Figure 16 shows the experimental results.
Table 10 shows the initial value settings of MPGA-related parameters. Figure 16 shows the corresponding experimental results. Through its powerful parallelism and global search capability, MPGA searches for the optimal combination of SVM penalty factor c and kernel parameters r in a relatively short time. MPGA obtains the optimal c and r as 40.3468 and 6.3841, respectively, to establish the optimal MPGA-SVM bearing fault diagnosis model and is verified by the test set. It can well identify four types of faults in bearing data: normal bearing, inner ring wear, outer ring wear and bearing missing ball. The fault diagnosis accuracy is as high as 95.325%. In the MPGA process, the average fitness values of population 1 to population 10 are increasing, and the average fitness values of individual subpopulations approach the optimal fitness from the sixth generation onwards. The average fitness values of some subpopulations are equal to the optimal fitness values at the later stage of the iteration, which fully reflects the unique advantages of multiple population. The unique advantage of the genetic algorithm is that the subpopulations are linked by the migration operator, which can introduce the optimal individuals into the neighboring subpopulations; thus, it achieves the purpose of co-evolution. In addition, using the artificial selection method, the optimal individuals generated in each generation are artificially selected, and an elite population is established. The elite population is independent of other populations, and the preserved individuals are no longer subject to genetic operations. Thus, the optimal solutions generated are protected, which effectively overcome the shortcomings of the standard genetic algorithm of “premature convergence”. To reduce the influence of chance factors on the experimental results of bearing fault diagnosis and effectively prove the feasibility of the proposed algorithm in the rotating machinery bearing fault diagnosis engineering, the experimental data were randomly generated by the random function of MATLAB software as the input of MPGA-SVM classification model. Experimental were conducted several times, and the experimental comparison results are shown in Table 11 and Table 12.
The results of six random experiments with the MPGA-SVM fault diagnosis model show that the MPGA-SVM model can still effectively identify and classify common bearing faults with a high accuracy rate in six random sample experiments with randomly selected bearing data samples. Through random samples and multiple experiments, the proposed fault diagnosis method based on the ICEEMDAN-wavelet threshold, mutual dimensionless index and MPGA-SVM in this paper can be effectively applied to the engineering problem of bearing fault identification and classification in petrochemical rotating machinery. Next, the experimental results of random samples of these three SVM models are summarized, compared and analyzed using the experimental results. Table 13 shows the comparative experimental results.
First, in terms of the accuracy of bearing fault identification, the experimental results in Table 13 show that the MPGA-SVM model has obviously better accuracy than both GA-SVM and SVM. Thus, MPGA is more accurate in searching for the optimal solutions of the parameter penalty factor and kernel function parameter that affect the classification performance of SVM and MPGA constructs a more optimized bearing fault diagnosis model.
Second, in terms of the model training time, the MPGA-SVM model is more efficient than the GA-SVM and the other two models of SVM; the MPGA-SVM model has the shortest training time. With its powerful global parallel search capability, the MPGA can search for the parameter penalty factor and kernel function parameter that affect the SVM classification performance in a relatively short time by co-evolving multiple populations.
Finally, GA-SVM has shorter training time than the traditional SVM mainly because the traditional SVM uses the grid search method or cross-test method to find c and r. Moreover, these methods have certain blindness in the range of values of parameters to be found and the setting of the step size. When the step size is too small, the algorithm will take a lot of time to run; when the step size is too large, the search will miss the optimal solution. Therefore, in practical engineering applications, it is not desirable to use the traditional grid search method or cross-test method to find the optimal SVM-related parameters when encountering training samples with large amounts of data. As observed, the MPGA-SVM model is superior.
In order to further exclude the specificity of the dataset for this study, the Western Reserve University bearing failure open-source dataset was selected for experimental comparison, which is a common dataset for industrial bearing failure detection and diagnosis The dataset is from the Swedish Prognostics and Health Management (PHM) challenge and was collected, processed and labeled by professionals. It contains 4 different fault categories, 10 operating states and various data types such as acceleration signals, temperature signals and vibration signals. The drive end bearing is SKF6205, with a 0.007 inch radial clearance bearing as the object. Consider the following working conditions: sampling frequency is 12 KHz, load is 0, the fault occurred in the inner ring, rolling body and outer ring in three positions, speed is 1797 r/min, vibration signal to build the data set, with the application of ICEEMDAN-wavelet min value after noise reduction of the bearing signal components of each order. The corresponding correlation coefficient values are shown in Table 14. The correlation coefficients of IMF7 and IMF8 are lower than 0.1, so IMF1–IMF6 components are selected to reconstruct the data set and construct mutual dimensionless indicators to extract features.
Some of the mutual dimensionless index data are shown in Table 15, and this data set is used as the input of the MPGA-SVM model, and then the fault identification experiments are carried out to verify the effectiveness of the fault diagnosis model.
The MPGA-SVM model also searched the optimal combination of SVM penalty factor c and kernel parameter r in a relatively short time, and obtained the optimal solutions of c and r as 29.4533 and 5.8659, respectively. The optimal MPGA-SVM bearing fault diagnosis model was established and verified by the test set. It can identify four kinds of faults in the bearing data: normal bearing, inner ring wear, outer ring wear and lack of rolling body, and the fault diagnosis accuracy is as high as 97.1%.The prediction results of each fault are shown in Table 16.
The results prove that the method proposed in this paper is accurate and effective in the Western Reserve University bearing failure open-source dataset.

4. Conclusions

This paper addressed the problem of serious noise interference in rolling bearing vibration data, proposed an improved adaptive noise-complete empirical modal decomposition with wavelet threshold joint noise-reduction (ICEEMDAN-wavelet threshold) method to optimize the data set and completed the fault feature extraction of bearing signals through mutual dimensionless processing of data. Finally, a multi-population genetic algorithm was proposed to optimize the classification model of the SVM for petrochemical rotating bearing fault identification. The efficient fault diagnosis model for the rolling bearing of petrochemical rotating machinery is 95.3%, which is 5% higher than the fault diagnosis accuracy of the traditional GA-SVM, by using the powerful global parallel search capability of MPGA itself to search for key parameters and optimize the performance of the SVM.

Author Contributions

Conceptualization, M.L.; methodology, C.M. and H.H.; software, C.M.; validation, C.M. and T.Y.; formal analysis, C.M.; investigation, C.M.; resources, M.L. and Q.Z.; data curation, C.M.; writing—original draft preparation, F.Z. and H.H.; writing—review and editing, F.Z., H.H. and M.L.; visualization, H.H.; supervision, M.L.; funding acquisition, M.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China, grant number 62073091, the Key fields (new generation information technology) special projects of universities in Guangdong Province, grant number 2020ZDZX3042, the National Natural Science Foundation of China, grant number 61933013, and the Open Fund of Hunan Provincial Key Laboratory of Mechanical Equipment Health Maintenance, grant number 21903.

Data Availability Statement

The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the privacy policy of the organization.

Acknowledgments

The authors would like to thank the editors and anonymous reviewers for a lot of useful suggestions.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. He, W.; Yi, C.; Li, Y.; Xiao, H. Research on Mechanical Fault Diagnosis Scheme Based on Improved Wavelet Total Variation Denoising. Shock Vib. 2016, 2016, 3151802. [Google Scholar] [CrossRef] [Green Version]
  2. Luo, Y.; Chen, C.; Zhao, S.; Kong, X.; Wang, Z. Rolling Bearing Diagnosis Based on Adaptive Probabilistic PCA and the Enhanced Morphological Filter. Shock Vib. 2020, 2020, 8828517. [Google Scholar] [CrossRef]
  3. Yuan, Z.; Peng, T.; An, D.; Cristea, D.; Pop, M.A. Rolling bearing fault diagnosis based on adaptive smooth ITD and MF-DFA method. J. Low Freq. Noise Vib. Act. Control 2020, 39, 968–986. [Google Scholar] [CrossRef] [Green Version]
  4. Huang, Y.; Jin, W.; Li, L. A New Wavelet Shrinkage Approach for Denoising Nonlinear Time Series and Improving Bearing Fault Diagnosis. IEEE Sensors J. 2022, 22, 5952–5961. [Google Scholar] [CrossRef]
  5. Ren, Y.; Li, W.; Zhu, Z.; Jiang, F. ISVD-Based In-Band Noise Reduction Approach Combined With Envelope Order Analysis for Rolling Bearing Vibration Monitoring Under Varying Speed Conditions. IEEE Access 2019, 7, 32072–32084. [Google Scholar] [CrossRef]
  6. Nguyen, P.; Kang, M.; Kim, J.M.; Ahn, B.H.; Ha, J.M.; Choi, B.K. Robust condition monitoring of rolling element bearings using de-noising and envelope analysis with signal decomposition techniques. Expert Syst. Appl. 2015, 42, 9024–9032. [Google Scholar] [CrossRef]
  7. Ma, J.; Li, H.; Tang, B.; Wang, J.; Zou, Z.; Zhang, M. Rolling bearing fault diagnosis based on improved VMD-adaptive wavelet threshold joint noise reduction. Adv. Mech. Eng. 2022, 14, 168781322211283. [Google Scholar] [CrossRef]
  8. Hu, C.; Xing, F.; Pan, S.; Yuan, R.; Lv, Y. Fault Diagnosis of Rolling Bearings Based on Variational Mode Decomposition and Genetic Algorithm-Optimized Wavelet Threshold Denoising. Machines 2022, 10, 649. [Google Scholar] [CrossRef]
  9. Lee, C.Y.; Zhuo, G.L. Localization of Rolling Element Faults Using Improved Binary Particle Swarm Optimization Algorithm for Feature Selection Task. Mathematics 2021, 9, 2302. [Google Scholar] [CrossRef]
  10. Yu, X.; Dong, F.; Ding, E.; Wu, S.; Fan, C. Rolling Bearing Fault Diagnosis Using Modified LFDA and EMD With Sensitive Feature Selection. IEEE Access 2018, 6, 3715–3730. [Google Scholar] [CrossRef]
  11. Ju, B.; Zhang, H.; Liu, Y.; Liu, F.; Lu, S.; Dai, Z. A Feature Extraction Method Using Improved Multi-Scale Entropy for Rolling Bearing Fault Diagnosis. Entropy 2018, 20, 212. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Xue, Q.; Xu, B.; He, C.; Liu, F.; Ju, B.; Lu, S.; Liu, Y. Feature Extraction Using Hierarchical Dispersion Entropy for Rolling Bearing Fault Diagnosis. IEEE Trans. Instrum. Meas. 2021, 70, 1–11. [Google Scholar] [CrossRef]
  13. Kuncan, M.; Kaplan, K.; Minaz, M.R.; Kaya, Y.; Ertunç, H.M. A novel feature extraction method for bearing fault classification with one dimensional ternary patterns. ISA Trans. 2020, 100, 346–357. [Google Scholar] [CrossRef] [PubMed]
  14. Li, W.; Zhu, Z.; Jiang, F.; Zhou, G.; Chen, G. Fault diagnosis of rotating machinery with a novel statistical feature extraction and evaluation method. Mech. Syst. Signal Process. 2015, 50–51, 414–426. [Google Scholar] [CrossRef]
  15. He, W.; Ding, Y.; Zi, Y.; Selesnick, I.W. Repetitive transients extraction algorithm for detecting bearing faults. Mech. Syst. Signal Process. 2017, 84, 227–244. [Google Scholar] [CrossRef] [Green Version]
  16. Zhu, H.; He, Z.; Wei, J.; Wang, J.; Zhou, H. Bearing Fault Feature Extraction and Fault Diagnosis Method Based on Feature Fusion. Sensors 2021, 21, 2524. [Google Scholar] [CrossRef] [PubMed]
  17. Xia, Z.; Xia, S.; Wan, L.; Cai, S. Spectral Regression Based Fault Feature Extraction for Bearing Accelerometer Sensor Signals. Sensors 2012, 12, 13694–13719. [Google Scholar] [CrossRef] [Green Version]
  18. Kumar, P.; Kumar, P.; Hati, A.S.; Kim, H.S. Deep Transfer Learning Framework for Bearing Fault Detection in Motors. Mathematics 2022, 10, 4683. [Google Scholar] [CrossRef]
  19. Toumi, Y. FPGA Implementation of a Bearing Fault Classification System Based on an Envelope Analysis and Artificial Neural Network. Arab. J. Sci. Eng. 2022, 47, 13955–13977. [Google Scholar] [CrossRef]
  20. Yan, X.; Jia, M. A novel optimized SVM classification algorithm with multi-domain feature and its application to fault diagnosis of rolling bearing. Neurocomputing 2018, 313, 47–64. [Google Scholar] [CrossRef]
  21. Liu, W.; Han, J. Rolling Element Bearing Fault Recognition Approach Based on Fuzzy Clustering Bispectrum Estimation. Shock Vib. 2013, 20, 213–225. [Google Scholar] [CrossRef]
  22. Wan, S.; Peng, B. An Integrated Approach Based on Swarm Decomposition, Morphology Envelope Dispersion Entropy, and Random Forest for Multi-Fault Recognition of Rolling Bearing. Entropy 2019, 21, 354. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  23. Fatima, S.; Guduri, B.; Mohanty, A.; Naikan, V. Transducer invariant multi-class fault classification in a rotor-bearing system using support vector machines. Measurement 2014, 58, 363–374. [Google Scholar] [CrossRef]
  24. Jiang, Q.; Chang, F.; Sheng, B. Bearing Fault Classification Based on Convolutional Neural Network in Noise Environment. IEEE Access 2019, 7, 69795–69807. [Google Scholar] [CrossRef]
  25. Tajeddini, M.A.; Aalipour, A.; Safarinejadian, B. Decision Fusion Method for Bearing Faults Classification Based on Wavelet Denoising and Dempster–Shafer Theory. Iran. J. Sci. Technol. Trans. Electr. Eng. 2019, 43, 295–305. [Google Scholar] [CrossRef]
  26. Gong, J.; Yang, X.; Wang, H.; Shen, J.; Liu, W.; Zhou, F. Coordinated method fusing improved bubble entropy and artificial Gorilla Troops Optimizer optimized KELM for rolling bearing fault diagnosis. Appl. Acoust. 2022, 195, 108844. [Google Scholar] [CrossRef]
  27. Attoui, I.; Fergani, N.; Boutasseta, N.; Oudjani, B.; Deliou, A. A new time–frequency method for identification and classification of ball bearing faults. J. Sound Vib. 2017, 397, 241–265. [Google Scholar] [CrossRef]
  28. Tan, H.; Xie, S.; Liu, R.; Ma, W. Bearing fault identification based on stacking modified composite multiscale dispersion entropy and optimised support vector machine. Measurement 2021, 186, 110180. [Google Scholar] [CrossRef]
  29. Lu, C.; Wang, Z.; Zhou, B. Intelligent fault diagnosis of rolling bearing using hierarchical convolutional network based health state classification. Adv. Eng. Inform. 2017, 32, 139–151. [Google Scholar] [CrossRef]
  30. Tong, Z.; Li, W.; Zhang, B.; Gao, H.; Zhu, X.; Zio, E. Bearing Fault Diagnosis Based on Discriminant Analysis Using Multi-View Learning. Mathematics 2022, 10, 3889. [Google Scholar] [CrossRef]
  31. Yang, F.; Kou, Z.; Wu, J.; Li, T. Application of Mutual Information-Sample Entropy Based MED-ICEEMDAN De-Noising Scheme for Weak Fault Diagnosis of Hoist Bearing. Entropy 2018, 20, 667. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  32. Zou, X.; He, D.; Jin, Z.; Wei, Z.; Miao, J. Intelligent diagnosis method of bearing fault based on ICEEMDAN and Ghost-IRCNN. Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci. 2022, 237. [Google Scholar] [CrossRef]
  33. Zhao, L.; Zhang, Y.; Zhu, D. Feature extraction for rolling element bearing weak fault based on MOMEDA and ICEEMDAN. J. Vibroeng. 2018, 20, 2352–2362. [Google Scholar] [CrossRef] [Green Version]
  34. Li, C.; Peng, T.; Zhu, Y.; Lu, S. Noise reduction method of shearer’s cutting sound signal under strong background noise. Meas. Control 2022, 55, 783–794. [Google Scholar] [CrossRef]
  35. Qin, J.; Fei, T.; Qu, Z.; Zhang, Y.; Yu, X. Fault Diagnosis of Rolling Bearing Based on ICEEMDAN and SSA-RVM. J. Phys. 2023, 2419, 012077. [Google Scholar] [CrossRef]
  36. Xiong, J.; Zhang, Q.; Peng, Z.; Sun, G.; Xu, W.; Wang, Q. A Diagnosis Method for Rotation Machinery Faults Based on Dimensionless Indexes Combined with K -Nearest Neighbor Algorithm. Math. Probl. Eng. 2015, 2015, 563954. [Google Scholar] [CrossRef] [Green Version]
Figure 1. General technical route of fault diagnosis.
Figure 1. General technical route of fault diagnosis.
Mathematics 11 02864 g001
Figure 2. Wavelet threshold denoising process.
Figure 2. Wavelet threshold denoising process.
Mathematics 11 02864 g002
Figure 3. ICEEMDAN-wavelet threshold joint noise-reduction flowchart.
Figure 3. ICEEMDAN-wavelet threshold joint noise-reduction flowchart.
Mathematics 11 02864 g003
Figure 4. Observable vibration signal model.
Figure 4. Observable vibration signal model.
Mathematics 11 02864 g004
Figure 5. Fault feature extraction process.
Figure 5. Fault feature extraction process.
Mathematics 11 02864 g005
Figure 6. Schematic diagram of the optimal classification surface of a support vector machine.
Figure 6. Schematic diagram of the optimal classification surface of a support vector machine.
Mathematics 11 02864 g006
Figure 7. Principal block diagram of the multi-population genetic algorithm.
Figure 7. Principal block diagram of the multi-population genetic algorithm.
Mathematics 11 02864 g007
Figure 8. Multi-population genetic algorithm optimization support vector machine flow chart.
Figure 8. Multi-population genetic algorithm optimization support vector machine flow chart.
Mathematics 11 02864 g008
Figure 9. Multistage centrifugal fan fault diagnosis platform.
Figure 9. Multistage centrifugal fan fault diagnosis platform.
Mathematics 11 02864 g009
Figure 10. Normal and faulty bearings.
Figure 10. Normal and faulty bearings.
Mathematics 11 02864 g010
Figure 11. ICEEMDAN decomposition result of normal bearing data.
Figure 11. ICEEMDAN decomposition result of normal bearing data.
Mathematics 11 02864 g011
Figure 12. ICEEMDAN decomposition results of bearing inner ring wear.
Figure 12. ICEEMDAN decomposition results of bearing inner ring wear.
Mathematics 11 02864 g012
Figure 13. ICEEMDAN decomposition results of bearing outer ring wear.
Figure 13. ICEEMDAN decomposition results of bearing outer ring wear.
Mathematics 11 02864 g013
Figure 14. ICEEMDAN decomposition result of bearing lacking a rolling body.
Figure 14. ICEEMDAN decomposition result of bearing lacking a rolling body.
Mathematics 11 02864 g014
Figure 15. MPGA-SVM evolution process.
Figure 15. MPGA-SVM evolution process.
Mathematics 11 02864 g015
Figure 16. MPGA-SVM model fault diagnosis results.
Figure 16. MPGA-SVM model fault diagnosis results.
Mathematics 11 02864 g016
Table 1. Correlation coefficient value corresponding to each order component of bearing signal.
Table 1. Correlation coefficient value corresponding to each order component of bearing signal.
Normal Bearing Correlation CoefficientInner Ring Wear Correlation CoefficientOuter Ring Wear Correlation CoefficientLack of Scrolling Body Correlation Coefficient
IMF10.28010.25780.26110.2454
IMF20.40120.29310.34960.2788
IMF30.82580.80160.79800.8213
IMF40.78620.90040.79520.9008
IMF50.41780.11930.44730.1319
IMF60.42220.14440.44400.1529
IMF70.12640.14060.11460.1463
IMF80.11920.12380.11650.1457
IMF90.11900.11330.12010.1348
IMF100.08750.07830.07520.1030
IMF110.07940.05470.05160.0787
IMF120.05400.03790.02990.0467
IMF130.01940.02860.01530.0367
IMF140.01790.02270.02120.0033
Table 2. Performance indicators of different noise-reduction methods.
Table 2. Performance indicators of different noise-reduction methods.
Noise Reduction MethodWavelet Threshold DenoisingICEEMDAN-Wavelet Threshold Joint Denoising
Fault TypeRMSESNR/dBRMSESNR/dB
Normal bearing0.321716.32260.23519.3625
Inner ring wear0.461617.23150.312521.3424
Outer ring wear0.541815.23170.216718.2313
Lack of rolling body0.468919.21540.221920.1524
Table 3. Mutual dimensionless index values of the bearings in different fault states (part).
Table 3. Mutual dimensionless index values of the bearings in different fault states (part).
Bearing StatusSampleWaveform IndicatorsPulse IndicatorYield IndicatorPeak IndicatorsCliffness
Normal bearing10.65410.15350.13240.23860.1598
20.67680.15470.11500.22310.1788
30.68310.16190.12560.23710.568
40.69880.18370.15620.26240.6587
50.69750.14570.12270.25540.6587
Outer ring wear10.75710.21360.16700.27850.2847
20.76520.25820.21310.36290.2335
30.77540.29820.25030.38170.2584
40.76520.26340.20790.32710.2964
50.78650.25510.22450.33710.2567
Inner ring wear10.53210.27420.22470.35580.2585
20.77850.26850.22650.34020.3254
30.75680.24580.20600.31300.3287
40.87420.32840.29970.40970.3261
50.77660.23700.19190.30730.2562
lack of rolling body10.74250.21210.17560.27730.2564
20.76520.21570.17950.27890.2738
30.74210.21040.17200.27890.2365
40.71350.11330.09060.15880.2451
50.76230.15970.13240.20970.2134
Table 4. Correspondence of the fault types, numbers and categories.
Table 4. Correspondence of the fault types, numbers and categories.
CategoryCategory1Category2Category3Category4
Fault number1–8081–160161–240241–320
Table 5. SVM model fault diagnosis results.
Table 5. SVM model fault diagnosis results.
Number of ExperimentsFirst Time
Failure CategoryTest NumberNumber of Confirmed CasesAccuracy Rate
Normal bearing
Inner ring wear
Outer ring wear
Lack of rolling body
Total
38
42
40
40
160
33
35
34
30
152
86.8%
83.3%
85%
75%
82.5%
Table 6. Experimental results of the SVM multiple fault identification.
Table 6. Experimental results of the SVM multiple fault identification.
Number of ExperimentsPenalty Factor cKernel Function Parameters rAlgorithm Running TimeDiagnostic Accuracy
177.29620.369367.1106 s80.625%
28.21292.276956.1203 s78.33%
316.61726.506558.6515 s82.08%
459.73920.660070.8597 s77.69%
559.87154.086159.5072 s79.52%
648.12455.436863.1245 s76.82%
Table 7. Genetic algorithm related parameter settings.
Table 7. Genetic algorithm related parameter settings.
Number of IndividualsMaximum Number of Genetic GenerationsCrossover ProbabilityProbability of VariationGeneration Gap Factor
40400.70.050.9
Table 8. GA-SVM model fault diagnosis results.
Table 8. GA-SVM model fault diagnosis results.
Number of ExperimentsFirst Time
Failure CategoryTest NumberNumber of Confirmed CasesAccuracy Rate
Normal bearing
Inner ring wear
Outer ring wear
Lack of rolling body
Total
40
35
41
44
160
36
32
37
39
144
90%
91.4%
90.2%
88.6%
90%
Table 9. Prediction results of each fault of GA-SVM model.
Table 9. Prediction results of each fault of GA-SVM model.
Number of ExperimentsPenalty Factor cKernel Function Parameters rAlgorithm Running TimeDiagnostic Accuracy
125.21130.042517.4612 s89.875%
28.21292.276916.1203 s88.333%
316.61726.506518.6515 s82.508%
459.73920.660016.8597 s87.269%
559.87154.086119.5072 s89.532%
645.24579.843118.2541 s85.457%
Table 10. MPGA related parameter settings.
Table 10. MPGA related parameter settings.
Number of PopulationsNumber of IndividualsNumber of Binary Bits of VariablesCrossover ProbabilityProbability of VariationGeneration Gap Factor
1020400.70.050.9
Table 11. Comparison of three random sample fault diagnosis experiments (First experiment).
Table 11. Comparison of three random sample fault diagnosis experiments (First experiment).
Number of ExperimentsFirst TimeSecond TimeThird Time
Failure CategoryNumber of TestsNumber of Confirmed DiagnosisAccuracy RateNumber of TestsNumber of Confirmed DiagnosisAccuracy RateNumber of TestsNumber of Confirmed DiagnosisAccuracy Rate
Normal bearing4040100%3838100%4040100%
Inner ring wear454395.6%424095.2%434195.3%
Outer ring wear373594.5%363495.1%383694.7%
Lack of rolling body383489.5%444090.9%393590.4%
Total16015294.91%16015395.30%16015295.12%
Table 12. Comparison of three random sample fault diagnosis experiments (Second experiment).
Table 12. Comparison of three random sample fault diagnosis experiments (Second experiment).
Number of ExperimentsFirst TimeSecond TimeThird Time
Failure CategoryNumber of TestsNumber of Confirmed DiagnosisAccuracy RateNumber of TestsNumber of Confirmed DiagnosisAccuracy RateNumber of TestsNumber of Confirmed DiagnosisAccuracy Rate
Normal bearing4040100%393897.4%4141100%
Inner ring wear444397.7%413995.1%424095.2%
Outer ring wear383592.1%383592.1%373697.2%
Lack of rolling body383592.1%423992.8%403587.5%
Total16015395.47%16015194.35%16015394.97%
Table 13. Comparison of the results of three random sample experiments of different models.
Table 13. Comparison of the results of three random sample experiments of different models.
Model CategoryNumber of ExperimentsAccuracy RateTraining Timecr
SVM182.50%67.1106 s77.29620.3693
278.33%56.1203 s8.21292.2769
382.08%58.6515 s16.61726.5065
GA-SVM189.87%17.4612 s25.21130.0425
288.33%16.1203 s8.21292.2769
382.08%18.6515 s16.61726.5065
MPGA-SVM194.91%14.2346 s52.41267.8412
295.30%11.1608 s40.34686.3841
395.12%16.2457 s23.21491.2458
Table 14. Western Reserve University open source dataset correlation coefficient value corresponding to each order component of bearing signal.
Table 14. Western Reserve University open source dataset correlation coefficient value corresponding to each order component of bearing signal.
Inner Ring Wear Correlation CoefficientOuter Ring Wear Correlation CoefficientLack of Scrolling Body Correlation Coefficient
IMF10.31280.32010.3166
IMF20.67230.68310.7041
IMF30.82540.83810.8415
IMF40.38920.40110.4129
IMF50.19350.20980.1983
IMF60.13020.14270.1319
IMF70.05670.05820.0479
IMF80.04150.04640.0337
Table 15. Mutual dimensionless index values of bearings under different fault states (part).
Table 15. Mutual dimensionless index values of bearings under different fault states (part).
Bearing StatusSampleWaveform IndexPulse IndexMargin IndexPeak IndexKurtosis Index
Inner ring wear10.72410.22180.28190.34120.2614
20.75670.23710.30580.33080.3196
30.80750.26490.27690.30540.3055
Outer ring wear10.81510.24720.18820.29540.2151
20.79230.24100.20340.37690.2213
30.80630.21050.27550.39180.2428
Lack of rolling body10.81220.23470.15230.25490.2608
20.84690.28590.18870.27630.2577
30.75480.29560.17390.28010.2246
Table 16. Prediction results of each fault.
Table 16. Prediction results of each fault.
Number of ExperimentsPenalty Factor cKernel Function Parameters rAlgorithm Run TimeAlgorithm Run Time
129.45335.865914.2735 s97.1%
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Mo, C.; Han, H.; Liu, M.; Zhang, Q.; Yang, T.; Zhang, F. Research on SVM-Based Bearing Fault Diagnosis Modeling and Multiple Swarm Genetic Algorithm Parameter Identification Method. Mathematics 2023, 11, 2864. https://doi.org/10.3390/math11132864

AMA Style

Mo C, Han H, Liu M, Zhang Q, Yang T, Zhang F. Research on SVM-Based Bearing Fault Diagnosis Modeling and Multiple Swarm Genetic Algorithm Parameter Identification Method. Mathematics. 2023; 11(13):2864. https://doi.org/10.3390/math11132864

Chicago/Turabian Style

Mo, Changchun, Huizi Han, Mei Liu, Qinghua Zhang, Tao Yang, and Fei Zhang. 2023. "Research on SVM-Based Bearing Fault Diagnosis Modeling and Multiple Swarm Genetic Algorithm Parameter Identification Method" Mathematics 11, no. 13: 2864. https://doi.org/10.3390/math11132864

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop