Multi-model optimization algorithm based on signal feature parameters and its application in fault diagnosis

Abstract

To achieve accurate signal analysis, this paper designs a rule for constructing a multi-model fitness function by integrating fault feature ratio (FFR), index of orthogonality (IO), index of energy conservation (IEC), envelope entropy, and kurtosis. Based on this rule, a vibration signal analysis method named multi-objective optimization algorithm-based variational mode decomposition with Bayesian Information Criterion (MOA–VMD–BIC) is proposed. First, Model_A’s fitness function is constructed using envelope entropy and kurtosis, and an optimization algorithm is employed to obtain the optimal VMD parameters {k, α}. Simultaneously, the values of Model_B (built with FFR, IO, and IEC) are calculated. Second, the intersection of solutions from Model_A and Model_B is derived to determine the optimal parameters {k, α} for the multi-objective optimization framework. Next, the optimal parameters are applied to decompose the original signal X into intrinsic mode functions (IMFs) via VMD. A signal sequence is then constructed using X and its components, and the Bayesian Information Criterion (BIC) is utilized to estimate the number of vibration sources in X. Finally, the reconstructed vibration source signals are subjected to fast Fourier transform (FFT) and Hilbert envelope analysis to diagnose mechanical faults. Simulation signals and experimental datasets validate the effectiveness of the proposed method.

Introduction

Bearings are critical components in the transmission systems of rotating machinery such as wind turbines and railway systems, and their operational condition directly affects the safe operation of the entire equipment. During actual service, bearings are subjected to harsh operating conditions including low-speed heavy-load, alternating loads, poor lubrication conditions, and contact friction with other components, making them highly susceptible to failure. When a bearing develops defects or incurs damage due to friction with other elements, periodic impulses are generated in the vibration signals during repeated passes over the damaged area. ¹ These impulses serve as the primary basis for identifying rolling bearing faults.

Under actual operating conditions, the fault characteristics exhibited by vibration signals differ from those of laboratory-simulated periodic pulse faults, demonstrating non-stationary, nonlinear, and strongly coupled properties. ² The key to fault diagnosis lies in rapidly extracting fault features from such non-stationary, nonlinear, and strongly coupled vibration signals.

In response to this, scholars have investigated signal decomposition methods based on the nonlinear and non-stationary characteristics of vibration signals, and proposed an empirical mode decomposition (EMD) method capable of decomposing a signal into the sum of multiple intrinsic mode functions (IMFs) and a trend component. ³ However, scholars have identified that the application of EMD to signal decomposition may result in modal aliasing and end effects.^4,5 Therefore, Lei et al. ⁶ proposed an ensemble empirical mode decomposition (EEMD) method, which utilizes the properties of white noise to assist in decomposing vibration signals. The method combining EEMD with diagnostic models has been adopted by numerous scholars.^7–9 However, with the application of EEMD, scholars have found that although it can effectively suppress the aliasing phenomenon, it is essentially a recursive modal decomposition method and thus cannot avoid the end effect and modal aliasing phenomenon after all. To address this, Dragomiretskiy and Zosso ¹⁰ proposed a non-recursive adaptive variational mode decomposition (VMD) method. By constructing a constrained variational solution approach, VMD transforms the signal decomposition problem into an optimization problem, thereby decomposing the target signal into multiple component signals. Thus, some scholars have conducted reconstruction analysis based on VMD-decomposed signals,^11,12 or carried out relevant diagnostic work by combining with machine learning.^13,14

Although VMD is an effective method for analyzing nonlinear, non-stationary, and complex vibration signals and can effectively suppress mode mixing, it is necessary to determine the number of decomposition modes k, the quadratic penalty factor a, and the selection criteria for valid IMFs from the decomposed components prior to spectral analysis, in order to ensure accurate fault signal identification. Currently, to address this issue, researchers commonly adopt a hybrid approach combining population-based optimization algorithms with signal analysis metrics. For instance, Kong et al. ¹⁵ utilized the particle swarm optimization (PSO) algorithm to optimize the parameters of VMD, employing Wasserstein distance (WD) and kurtosis index to select IMFs and reconstruct vibration signals, thereby achieving effective fault feature extraction and identification. Similarly, Shi et al. ¹⁶ applied the differential search algorithm (DSA) to optimize VMD parameters, combined with the correlation kurtosis index for IMF screening and signal reconstruction, thereby achieving effective fault feature extraction and identification. Additionally, Chen and Zhao ¹⁷ implemented the whale optimization algorithm (WOA) for VMD parameter optimization, integrating threshold denoising techniques to accomplish fault feature extraction and identification. Numerous similar cases exist, yet these methods invariably utilize envelope entropy or isolated signal indicators (e.g. kurtosis) as the algorithm’s fitness function, neglecting interrelationships and mutual influences among multiple indicators. Furthermore, during the IMF screening process, analyses rely solely on signal indicators such as kurtosis, while ignoring the multi-source nature of signals.

To address the aforementioned issues, incorporates the time-frequency domain statistical indicators of vibration signals and constructs multiple models targeting the optimal solution, thereby reducing the error of single-indicator evaluation. In this paper, a comprehensive rule for constructing a multi-model fitness function integrating time-frequency domain statistical indicators was designed. An improved genetic algorithm (IGA) was utilized to optimize the parameters {k, a} of the VMD algorithm, ensuring global search capability and avoiding local optima. Using the Bayesian Information Criterion (BIC) algorithm to calculate the minimum cost of VMD, and reconstruct the signal based on the determined number of vibration sources of the vibration signal. Meanwhile, the envelope spectrum algorithm is used to analyze the reconstructed signal after fourth-order high-pass filtering, thereby realizing the fault feature extraction of the vibration signal. In this paper, the accuracy of the proposed method is evaluated by calculating the relative error between the extracted frequency and the target frequency. Compared with the conventional GA–VMD method, the algorithm proposed in this paper mainly has the following advantages:

1. It comprehensively considers the influence of statistical indicators and parameters of time-frequency domain signals on signal diagnosis, thereby avoiding the limitations imposed by a single indicator.

2. Through improvements to crossover and mutation operations, the IGA enhances the diversity of population iteration and mutation, overcoming the tendency of the traditional GA to fall into local optima.

3. It employs a multi-objective optimization model to define the fitness function of the optimization algorithm, and improves the reliability of the optimization objective through comprehensive indicators.

4. After optimizing VMD via the IGA, signal reconstruction and analysis are not performed directly. Instead, an additional BIC-based source identification layer is incorporated to filter out irrelevant signals, thereby improving the accuracy of the target signal.

Brief introduction of VMD, MOA, and BIC

Variational modal decomposition (VMD)

As a powerful tool in the field of signal processing, VMD is used to decompose a real-valued input signal X into a discrete number of sub-signals (modes) u _k , essentially based on the three concepts of Wiener Filtering, Hilbert Transform and Analytic Signal, Frequency Mixing and Heterodyne Demodulation. ¹⁰ Since VMD has carried out a detailed analysis and description in Dragomiretskiy and Zosso, ¹⁰ this paper will not repeat it, but only describes the implementation process of VMD. The process is described as follows:

• Step 1: Initialize modal components $u_k=u_1,\cdots ,u_k$ , the center frequency ${\omega }_k={\omega }_1,\cdots ,{\omega }_k$ , and the Lagrangian multipliers λ(t).

• Step 2: Update the modal components u _k , the center frequency ω _k , and the Lagrangian multipliers λ(t). The (n + 1)th iteration process is shown in equations (1) to (3). Where, α represents quadratic penalty factor, n represents nth iteration, τ represents relaxation coefficient, x represents vibration signal.

$$u_k^{n+1}(\omega )=\frac{x(\omega )-\sum _{i<k}{u}_k^{n+1}(\omega )-\sum _{i<k}{u}_k^n(\omega )+\frac{{\lambda }^n(\omega )}{2}}{1+2\alpha {(\omega -{\omega }_k^n)}^2}$$

(1)

$${\omega }_k^{n+1}=\frac{\underset{0}{\overset{\infty }{\int }}\omega {|u_k^{n+1}(\omega )|}^{2}d\omega }{\underset{0}{\overset{\infty }{\int }}{|u_k^{n+1}(\omega )|}^{2}d\omega }$$

(2)

$${\lambda }^{n+1}(\omega )={\lambda }^n(\omega )+\tau (x(\omega )-\sum _k{u}_k^{n+1}(\omega ))$$

(3)

• Step 3: If $\sum _k(\Vert u_k^{n+1}-u_k^n{\Vert }_2^2/\Vert u_k^n{\Vert }_2^2)<\epsilon$ , and the variational modal solution process ends, several modal components signal with finite bandwidth are obtained. Otherwise, Step 2 is repeated. Where, ε represents red limit value of error.

From the above derivations, the signal decompose can be converted to selection of the parameters {k, α} of VMD. Traditionally, they are optimized by evolutionary algorithms with a single fitness function.^15–18 To improve the robustness of optimal VMD, multi-method optimization is employed and presented in the next subsection.

Multi-model optimization algorithm (MOA)

In the most real-world optimization problems may involve more than one objective, these objectives are possibly conflict. In order to solve these objectives, numerical simulation methods are usually used. However, different from ordinary multi-objective optimization problems, it is necessary to comprehensively consider the time-frequency domain statistical indicators ¹⁹ in the process of solving, such as: fault feature ratio (FFR), the index of orthogonality (IO), and the index of energy conservation (IEC) in the process of solving.

• Definition 1: The FFR serves as a quantitative metric to select specific fault-related frequency components within a signal. A higher FFR value indicates that the selected frequency components contain richer periodic impulsive information. This information is critical for identifying mechanical faults. Such faults include bearing spalling and gear tooth breakage. The calculation of FFR is expressed as follows equation (4):

$$\mathrm{FF}{R}_f=\sum _{k=1}^{K}S(\mathrm{kf})/S$$

(4)

Where, in equation (4), f denotes the fault characteristic frequency. S represents the sum of amplitudes of the envelope spectrum of the time-domain signal X(t). And S(kf) indicates the envelope spectral amplitudes corresponding to each harmonic of the fault characteristic frequency.

• Definition 2: The IO and IEC are metrics used to evaluate the effectiveness of signal decomposition. Specifically, a smaller IO value indicates better orthogonality among signal components and higher decomposition accuracy, while an IEC value closer to 1 suggests fewer energy leakage during decomposition and more optimal decomposition outcomes. ²⁰ The computational methods for IO and IEC are defined in equations (5) and (6):

$$\mathrm{IO}=\sum _{t=0}^T\sum _{i=1}^{n+1}\sum _{j=1}^{n+1}{x}_i(t)x_j(t)/X^2(t)i\ne j$$

(5)

$$\mathrm{IEC}=\sum _{t=0}^T\left(\frac{\sum _{i=1}^n{x}_i^2(t)}{{(X(t)-r_n(t))}^2}\right)$$

(6)

Where, in equations (5) and (6), x _i (t) and x _j (t) represent the ith and jth components of signal X(t). X(t) denotes the original signal. And r _n (t) denotes the trend term.

• Definition 3: The envelope entropy value of a signal reflects the periodic fault characteristic information contained within. If the signal contains a higher amount of fault characteristic information, it indicates greater sparsity, leading to a smaller envelope entropy value. Conversely, if the signal contains fewer fault characteristics, it cannot exhibit significant sparsity, resulting in a larger envelope entropy value. The calculation method for envelope entropy is shown in equation (7):

$$E{E}_i=\sum _{i=1}^N(\mathrm{Hilbert}(x_i(t))\lg \mathrm{Hilbert}(x_i(t)))$$

(7)

Where, in equation (7), x _i (t) represent the ith component of signal X(t). Hilbert(·) denotes the envelope signal obtained through Hilbert transform of a signal.

• Definition 4: Kurtosis is a dimensionless parameter. This parameter exhibits high sensitivity to subtle variations in fault signals. Such sensitivity makes kurtosis particularly suitable for early fault diagnosis of bearings. Under normal bearing operation, vibration signals approximate a normal distribution, with kurtosis values approaching 3. Localized faults may develop in rolling bearings. When this happens, vibration signals deviate from normality. The deviation is caused by impact-induced distortions. Meanwhile, fault-related information is predominantly contained within amplitude-modulated signals. These amplitude-modulated signals are generated by the aforementioned periodic impacts. Consequently, elevated kurtosis values indicate enhanced impact energy concentration and richer fault signatures, thereby facilitating more effective feature extraction. The kurtosis of component is mathematically defined as:

$$q{d}_i=\frac{E{(x_i(t)-u_i)}^4}{{\sigma }_i^4}$$

(8)

Where, in equation (8), x _i (t) represent the ith component of signal X(t). E(·) represents the expected value of the discrete signal sequence x _i (t). u _i represent the mean amplitude of the discrete signal sequence x _i (t). σ _i represent the standard deviation of the discrete signal sequence x_i(t).

When solving the optimal solution for signal decomposition, we integrate the characteristic indicators of the aforementioned signals to establish a multi-objective optimization model. This model comprises two sub-models (Model_A and Model_B), where Model_A serves as the foundation for iterative optimization of signal decomposition, while Model_B operates based on the computational results from Model_A. Each iteration of Model_A and Model_B generates a solution set. We first normalize the results of both models. Meanwhile, we use the parameter values {k, α} derived from Model_A as the baseline reference. On this basis, we identify the intersection set of the two models. Through this process, we obtain the optimal solution for the multi-objective optimization model. The proposed multi-objective optimization framework is formulated as follows:

$$\begin{array}{c}\mathrm{subject}\mathrm{to}:\{k,a\}=\mathrm{Model}A\cap \mathrm{Model}B\hfill \\ \mathrm{Model}A:(\left\{k\right\},\left\{a\right\},\left\{\mathrm{fitness}\right\})=F(\mathrm{fit},f)\hfill \\ \mathrm{Model}B:(\left\{k\right\},\left\{a\right\},\left\{\mathrm{value}\right\})=2+\mathrm{IO}\left(f\right)-\mathrm{IEC}\left(f\right)-\mathrm{FFR}\left(f\right)\hfill \\ \mathrm{fit}=min\sum _{i=1}^{N}E{E}_i+\frac{1}{|q{d}_i-3|}0<i\le N\hfill \\ f=\{\begin{array}{c}G(X(t),i,\beta )\mathrm{if}\mathrm{Model}\_A0<\beta \le M\hfill \\ \multicolumn{1}{c}{G(X(t),k,a)\mathrm{if}\mathrm{Model}\_B}\end{array}\hfill \end{array}$$

(9)

Where, k, {k} denote the number of modes and set of modes in the VMD algorithm, respectively. The a, {a} represent the regularization parameter and set of regularization parameters in the VMD algorithm. The function Model_A ∩ Model_B refers to the intersection of Model_A and Model_B, which corresponds to the {k, a} parameters shared by both models. The symbol ∩ is standardized as ∩ in set theory. F(·) represents the optimization algorithm (e.g. genetic algorithm). fit is the fitness value calculation function used to evaluate the performance of candidate solutions. {fitness} is the set of fitness values generated during the optimization process. f is a general function, and the parameters incorporated in its calculation depend on whether Model_A or Model_B is used. G(·) denotes the signal analysis algorithm VMD. X(t) is the raw input signal. i and β are the instantaneous values of the number of modes and regularization parameter, respectively, with constraints: 0 < i ≤ N (where N is the maximum mode count), 0 < β < M (where M is the upper bound for regularization). FFR(·) denote the fault feature ratio defined in equation (4). IO(·) denote the orthogonality index defined in equation (5). IEC(·) denote energy conservation index defined in equation (6).

To further elaborate on the multi-objective optimization model, equation (9) is further explained herein. Firstly, the objective of the algorithm is to find the intersection corresponding to Model_A and Model_B. The intersection of Model_A and Model_B serves as the optimal solution of the VMD algorithm. Secondly, the solution of Model_A is obtained as follows: by using an optimization algorithm (the optimized genetic algorithm is adopted in this paper), within the set ranges of the number of modes (0 < i ≤ N) and the regularization parameter (0 < β < M), the solution is calculated in accordance with the fitness function fit. The fitness function fit consists of the envelope quotient and kurtosis. (It should be noted that what is obtained from the solution of Model_A is a solution set {{k}, {a}, {fitness}} rather than a specific solution. The reason why it is a solution set instead of a specific solution lies in that during each population iteration and optimization process of the genetic algorithm, the optimal {a} will be found on the basis of {k}.) Thirdly, the solution of Model_B is derived as follows: the solution set {{k}, {a}} obtained from Model_A is substituted into the VMD algorithm. According to the calculation formula of Model_B, the {value} corresponding to the solution set {{k}, {a}} is obtained. Finally, with {k} as the baseline, the intersection of Model_A and Model_B is found, so as to achieve the optimal effect under the multi-model scenario. The solution of the algorithm can be referred to Table 1.

Table 1.

Multi-model optimization algorithm.

Algorithm 1: Multi-model optimization algorithm

1. Input: signal X, and initialized parameters of Model_A including i, β and parameters of optimization algorithm 2. Calculation: base on VMD to calculate {xi,i,β} = G(X,i,β), 0<i≤N,0<β<1 3. Procedure: according Model_A to calculate _fitness to obtain {k=i,a=β,fitness=_fitness}      for i= 0:N        for j = 0:Maxstep          rand(β,[0,M])         {i,β,_fitness}=F(fit,G(X,i,β))        if fitness > _fitness           fitness=_fitness,k=i,a=β            {k,a,fitness}            Model_B: {k,a,value}=2+IO(G(X,k,a))-IEC(G(X,k,a))-FFR(G(X,k,a))         else          fitness=fitness,k=i,a=β            {k,a,fitness}           Model_B: {k,a,value}=2+IO(G(X,k,a))-IEC(G(X,k,a))-FFR(G(X,k,a))        end       {k,a,fitness}       {k,a,value}      end 4. End procedure 5. Procedure: calculation the common solutions between Model_A and Model_B common solutions as solution of multi-Model optimization algorithm: {k,a} 6. End procedure

Bayesian Information Criterion (BIC)

Vibration signals are typically composed of multiple sources. To accurately analyze fault information in vibration signals, a clear understanding of their composition is essential. Therefore, this study employs the Bayesian information criterion (BIC) for signal source analysis. First, the signal X(t) is decomposed into k mode components u _k using VMD optimized by a multi-model optimization algorithm. Second, X(t) and the mode components u _k are combined to form a multidimensional signal matrix M  = (X(t), u ₁ , u ₂ , …, u _k ). The covariance matrix R of M is then computed, and singular value decomposition (SVD) is performed on R . Finally, based on the number of non-zero eigenvalues of R , the n minimum-cost components are calculated to derive the n sources of X(t). ²¹ The computational steps are as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{obj}.\mathrm{BIC}(n)={\left(\underset{j=1}{\overset{n}{\Pi }}{\lambda }_j\right)}^{-N/2}{\tilde{\sigma }}_n^{-N(m-n)/2}{N}^{-(d_n+n)/2}}\\ \multicolumn{1}{c}{s.t.\{\begin{array}{c}{\tilde{\sigma }}_n=\sqrt{(\sum _{j=n+1}^m{\lambda }_j)/(m-n)}\hfill \\ \multicolumn{1}{c}{d_n=\mathrm{mn}-n(n+m)/2}\end{array}}\end{array}$$

(10)

Where, m denotes the number of non-zero eigenvalues of the covariance matrix R . N represents the length of the data associated with R . And λj represents the jth eigenvalue of R .

Proposed MOA–VMD–BIC for extracting fault impulses signal

Analysis of VMD parameters

To demonstrate why multi-model is introduced to solve the problem of vibration signal feature extraction, a simulated signal ²² which can be as the research object to analyze the influence of parameter selection of VMD. The simulated signal can be expressed as equation (11):

$$\{\begin{array}{c}\multicolumn{1}{c}{s(t)=p(t)+h(t)+g(t)+n(t)}\\ \multicolumn{1}{c}{p(t)=\sum _i(1+A_0\cos (2\pi f_r{t}_{\mathrm{i}}))·\exp (-\lambda ·(t_i-i{T}_d-r_i))·\cos (2\pi f_n·(t_i-i{T}_d-r_i))}\\ \multicolumn{1}{c}{h(t)=\sum _i{R}_i·\exp (-\lambda ·(t_i-\epsilon ))·\cos (2\pi f_n·(t_i-\epsilon ))}\\ \multicolumn{1}{c}{g(t)=\sum _l{A}_l\cos (2\pi f_r{t}_i+{\theta }_l)+\sum _k{A}_k\cos (2\pi z{f}_r{t}_i+{\theta }_k)}\end{array}$$

(11)

Where, these four parts in equation (11) stand for bearing fault impulses signal p(t), high-amplitude impulses signal h(t), harmonic vibration of the rotor and cyclostationary interference from gear meshing g(t), and Gaussian noise n(t), respectively. f _r denotes the rotating frequency, f _n denotes the impulsive frequency, z denotes the gear teeth number, T _d is the fault period, r _i denotes the roller random slips of (−0.02 T _d , 0.02 T _d ). A₀, A_l, A _k represents the amplitude of bearing fault, rotor vibration, gear meshing, respectively. The values of all the parameters in equation (11) are listed in Table 2.

Table 2.

Multi-model optimization algorithm.

p(t)					h(t)					g(t)
A 0	f r	λ	T d	r i	f n	R i	λ	ε	f n	A l	fr	θ l	A k	z	θ k
0.3	10	560	1/80	(−0.02 T d , 0.02 T d )	2400	1.5	480	0.12	4500	0.1	10	π/2	0.2	36	π/2

The simulated signal s(t) and its components are plotted in Figure 1, in which the sampling frequency Fs is 10,000 Hz. And its frequency spectrum, squared envelope spectrum (SES) are plotted in Figure 2, in which f _m is gear meshing frequency.

Figure 1.

The simulated signal s(t) and its components.

Figure 2.

Frequency spectrum, squared envelope spectrum of simulated signal.

It can be seen from Figure 1 that in the simulated signal (5), bearing fault impulses signal (1) has been submerged by Gaussian noise (4). The rotational frequency and gear meshing frequency can be accurately identified in Figure 2, but the bearing fault frequency cannot be identified from frequency spectrum and SES. The value of envelope entropy can represent signal sparsity. And the lower the signal entropy value, the higher the sparsity of the signal, and the more periodic components the signal contains. So, the signal is decomposed to study the characteristics of IMF of VMD. Taking the envelope entropy of the IMF as an indicator, the signal entropy values for different values of k and α are plotted in Figure 3 (note: k ∈ (1, 15), α ∈ (1000, 5000)). The minimum entropy values corresponding to each k are plotted in Figure 4.

Figure 3.

Envelope entropy trend diagram under different {k, α}.

Figure 4.

Minimum envelope entropy for different k.

As illustrated in Figure 4, the minimum envelope entropy value is 5.167, with the optimal parameter combination {k, α} = {15, 5000}. These parameters were then input into the VMD algorithm to decompose the signal s(t). The components and spectral characteristics of signal s(t) are illustrated in Figure 5. As shown in Figure 5, the signal s(t) is over-decomposed, and the component signals only reveal the presence of the shaft rotational frequency (f _r ) and gear meshing frequency (f _m ), while failing to detect bearing fault frequencies. Therefore, when selecting parameters for VMD, it is insufficient to rely solely on envelope entropy, other factors must be comprehensively considered.

Figure 5.

VMD decomposition and spectrum analysis of signal s(t).

Multi-model optimization algorithm (MOA) and Bayesian Information Criterion (BIC)

To address the aforementioned issues, this paper proposes a multi-model optimization algorithm based on signal characteristic parameters. The optimization algorithm employs a multi-model optimization approach (as detailed in Section 2.2) to obtain the VMD parameters {k, α}, combined with the Bayesian Information Criterion (BIC) to identify the vibration source of signal s(t). Based on the vibration source, s(t) is reconstructed, and spectral analysis is applied to extract fault characteristics from the reconstructed signal. In order to facilitate memory, the method proposed in this paper is referred to as MOA–VMD–BIC, and its algorithm flow chart is as follows (Figure 6):

Figure 6.

Flow chart of MOA–VMD–BIC model algorithm.

Note that the optimization algorithm F(·) used in the MOA–VMD–BIC model is the optimized genetic algorithm (IGA) proposed in Pi et al. ²³ Specifically, the function F(·) in equation (9) should be substituted with this optimized genetic algorithm (IGA). The computational steps of the MOA–VMD–BIC model are as follows:

• Step 1: The number of modes k VMD is set to range from 1 to 15, and the regularization parameter α is set to range from 1000 to 5000, with iterations starting from 1. The parameters of the improved genetic algorithm (IGA) are initialized (e.g. maximum number of iterations maxgen = 100, population size sizepop = 20, crossover probability p_cross = 0.6, mutation probability p_mutation = 0.05, custom weight w = 0.65), where the number of modes k and regularization parameter α of VMD serve as the optimization targets of the IGA.

• Step 2: Perform crossover operation according to equation (12), mutation operation according to equation (13), and generate new individuals.

$$\begin{array}{c}\mathrm{obj}.{\mathit{s}}_{\mathrm{best}}={\mathit{s}}_j\mathrm{when}\underset{j=1}{\overset{4}{min}}f\mathrm{it}({\mathbf{s}}_j)\hfill \\ s.t.\{\begin{array}{c}\multicolumn{1}{c}{{\mathit{s}}_1=\{({\mathit{f}}_i+{\mathit{m}}_i)/2\}}\\ \multicolumn{1}{c}{{\mathit{s}}_2=\{g_{max}(1-w)+w·max({\mathit{f}}_i+{\mathit{m}}_i)\}}\\ \multicolumn{1}{c}{{\mathit{s}}_3=\{g_{min}(1-w)+w·min({\mathit{f}}_i+{\mathit{m}}_i)\}}\\ \multicolumn{1}{c}{{\mathit{s}}_4=\left\{\frac{(g_{max}+g_{min})·(1-w)+w·({\mathit{f}}_i+{\mathit{m}}_i)}{2}\right\}}\end{array}\hfill \end{array}$$

(12)

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{obj}.\mathit{s}=p_i({\mathit{s}}_{\mathrm{best}})\mathrm{when}\underset{i=1}{\overset{3}{min}}f\mathrm{it}(p_i({\mathit{s}}_{\mathrm{best}}))}\\ \multicolumn{1}{c}{s.t.p_i({\mathit{s}}_{\mathrm{best}})=\{{\mathrm{ch}}_1{\mathrm{ch}}_2{\mathrm{ch}}_3\cdots {\mathrm{ch}}_{\mathrm{num}}\}}\\ \multicolumn{1}{c}{+\{{\beta }_1{\delta }_1{\beta }_2{\delta }_2{\beta }_3{\delta }_3\cdots {\beta }_{\mathrm{num}}{\delta }_{\mathrm{num}}\}}\end{array}$$

(13)

Where, f _i and mi denote the ith gene in the selected father and mother individuals’ chromosomes, respectively. g_max and g_min represent the maximum and minimum values of genes within the chromosome. s _j indicates the new individual chromosome generated through four distinct crossover strategies. s_best refers to the optimal individual after crossover operations, selected based on fitness evaluation. fit(·) is the fitness function, defined in equation (9). w is a user-defined weight within the range (0, 1). s denotes the optimal individual after mutation. p _i (·) represents three types of mutation operators, including local mutation, global mutation, and backward reasoning mutation, and detailed mutation rules are specified in Pi et al. ²³ β takes only binary values 0 or 1. δ is a random number uniformly distributed in (−0.1, 0.1). ch _l (l = 1, 2, …, num) denotes the gene value of the individual, where num is the total number of genes in the chromosome.

• Step 3: Record the optimal parameter {k, α} and fitness value, and obtain the solution vector A {k, α, fitness} of Model_A in a single iteration.

• Step 4: The optimal parameters {k, α} are substituted into VMD to decompose the signal s(t), and the target value of Model_B is calculated. Thus, the solution vector B {k, α, value} of Model_B in a single iteration is obtained.

• Step 5: In each iteration, the solution vectors of Model_A and Model_B, are respectively, constructed into solution matrix A and B, respectively. Check whether the termination conditions (maxstep = maxgen and k ≤ 15) are met. If no, return to Step 2. If so, the output matrix A and B.

• Step 6: The fitness and value of matrix A and B are normalized, respectively, and the corresponding {k_best, a_best} when fitness = value is solved.

Since the {k, α} calculated by IGA may not be arranged in order during the solution process, the matrix A and B are sorted according to k in the solution first. The solution sets of Model_A and Model_B, as well as the time consumed for solving the solution sets, are presented in Table 3. After applying smoothing processing to these solution sets, their graphical representations are illustrated in Figure 7. As shown in Figure 7, both Model_A and Model_B exhibit intersecting solutions at the modal components of VMD with k = 5.525 and k = 10.64. It can be observed from Table 3 that when k = 10.62, the fitness values of Model_A and Model_B are much larger than those when k = 5.525. Additionally, the computation time required for k = 10.62 is also longer than that for k = 5.525. According to the definitions of Model_A and Model_B in equation (9), the parameter selection should prioritize the minimum fitness value. Meanwhile, computational complexity, that is, computation time, should be taken as an indicator for parameter selection. Based on the above analysis, the value of k = 5. (Note: Since the number of modal components is an integer, the value is directly rounded during the solution process.)

Table 3.

The solutions of Model_A and Model_B.

No.	The solution of VMD {k, a}	The fitness of Model_A (normalization)	The value of Model_B (normalization)	Calculation time consumption (s)
1	(1, 2065)	1	0.510841780	97.57196572
2	(2, 2250)	0.882767596	0.489177418	99.25192599
3	(3, 2033)	0.775056896	0.488738003	112.2883463
4	(4, 3806)	0.484637282	0.198037309	113.3846554
5	(5, 4155)	0.314436722	0.184744127	114.1451714
6	(6, 2229)	0.259868986	0.471713847	115.2572631
7	(7, 2700)	0.138427264	0.681739120	115.4912426
8	(8, 1653)	0.423583611	0.704142585	120.1943283
9	(9, 2938)	0.641764017	0.725143243	120.2677429
10	(10, 1137)	0.727845764	0.739785704	127.0668625
11	(11, 4241)	0.767134956	0.769302630	135.2064879
12	(12, 4565)	0.723553154	0.796796970	154.4600914
13	(13, 4338)	0.668215826	1	155.9750478
14	(14, 2437)	0.386286050	0.833149013	161.4072389
15	(15, 2005)	0.259547282	0.794548424	211.2690515

Figure 7.

The solutions of Model_A and Model_B.

After applying the optimal decomposition parameters (k = 5, a = 4155) to the VMD, the signal s(t) is decomposed into five IMFs, denoted as IMF  = {imf ₁ , imf ₂ , imf ₃ , imf ₄ , imf ₅ }. The signal s(t) and the IMF are combined to form a multidimensional signal matrix M  = {s(t), imf ₁ , imf ₂ , imf ₃ , imf ₄ , imf ₅ }. Based on the BIC algorithm described in Section 2.3, the covariance matrix R of M is computed. Subsequently, the matrix R is decomposed by singular value to obtain the eigenvalues and target costs as shown in Table 4.

Table 4.

Eigenvalue of matrix R and calculation results of BIC.

No.	Eigenvalues, λ i	BIC (k)
1	0.004333360	2.11 × 1018
2	0.004159969	5.81 × 1014
3	0.004073234	7.37 × 1011
4	0.001812689	2.47 × 1011
5	0.018181824	4.38 × 1022
6	0.133037359	2.023 × 1029

In Section 3.1, the simulated signal s(t) is composed of four components: bearing fault impulses signal, high-amplitude impulses signal, harmonic vibration of the rotor, and cyclostationary interference from gear meshing, Gaussian noise, indicating that the actual number of vibration sources equals 4. As shown in Table 4, the Bayesian objective cost reaches its minimum value when the estimated number of vibration sources is also 4. And the parameter selection comparison table is shown in Table 5.

Table 5.

Parameter k selects the comparison table.

Models	Number of signal sources	MOA–VMD	BIC
Value	4	5	4

Consequently, the intrinsic mode functions {imf ₁ , imf ₂ , imf ₃ , imf ₄ } were reconstructed into the vibration signal s′(t), and envelope spectrum analysis was performed on the reconstructed signal s′(t). The results of the envelope spectrum analysis are illustrated in Figure 8. According to the simulation signal expression in equation (11), it can be seen that the bearing fault frequency is 80 Hz. And its harmonics can now be easily identified from the SES shown in Figure 8.

Figure 8.

Envelope spectrum analysis results based on MOA–VMD–BIC model.

Applications in fault diagnosis of bearings

Bearing experiment description

The dataset used in this study was sourced from Case Western Reserve University (CWRU). ²⁴ As shown in Figure 9, the experimental setup in the CWRU Bearing Test Rig contains 6205-2RS JEM SKF deep groove ball bearings installed at both the fan end and drive end of the motor. The bearing dimensions are illustrated in Figure 10. In the experimental, outer race fault diameter is 0.1778 mm (0.007 inches) and depth is 0.2794 mm (0.011 inches) at the drive-end bearing. In addition, the motor speed is 1797 rpm (fr = 1797/60 = 29.95 Hz), the sampling frequency is 12 kHz, and the sampling number is 10,240. According to equation (14), the outer ring bearing frequency is calculated as 107.36 Hz = (1797/60) × 3.5848 Hz, that is, 3.5848 times the rotation frequency.

Figure 9.

Schematic diagram of bearing test bench at Case Western Reserve University.

Figure 10.

SKF6205 JEM SKF.

6205-2RS JEM SKF deep groove ball bearing outer ring fault characteristic frequency calculation formula is as follows:

$$f_d=\frac{(n\times z)}{120}\times (1-\frac{l\times \cos \beta }{B})$$

(14)

Where, n represents the speed, and its unit is rpm. z is the number of rolling elements. l is the diameter of the rolling element, and its unit is mm. B represents the pitch circle diameter of the rolling bearing (B = (D₁ + d₁)/2), and its unit is mm. β is the contact angle. The collected time domain, frequency domain, and envelope spectrum signals of the faulty bearing are shown in Figure 11. As can be seen from Figure 11, the outer ring bearing frequency (107.36 Hz) cannot be directly obtained from the spectrum or the envelope spectrum. Therefore, further analysis of the signal is necessary in order to obtain the outer ring bearing frequency (107.36 Hz)

Figure 11.

The time domain, frequency domain, and envelope spectrum signals of 6205-2RS JEM SKF.

Vibration signal analysis of 6205-2RS JEM SKF

Introduction of comparison algorithm

To validate the effectiveness of the proposed algorithm, we adopted widely recognized methodologies in academia as benchmark comparisons. ²⁵ The implementation process is structured as follows, first the minimum envelope entropy of the signal was utilized as the fitness function for the optimization algorithm to determine the optimal decomposition parameters {k, a} of VMD. The kurtosis values of the component signals were computed based on equation (15), and components with kurtosis values exceeding 3 were selected for signal reconstruction to enhance feature extraction.

$$k=\frac{\sum _{i=1}^n{(x_i-\mu )}^4/n}{{\left({\sigma }^2\right)}^2}-3$$

(15)

Where, μ represents the average value of the component signals. σ ² represents the variance of the signal component data. x _i represents the ith point of the component signals data. (Note: When the kurtosis in the signal is significantly higher than 3 (the kurtosis of a standard normal distribution), this indicates a nonlinear effect or local fault in the signal.)

In the comparative algorithm, the optimization algorithm, data preprocessing method (DC component removal), and filter are identical to those employed in the proposed method. The distinction between them lies in the fitness function of the optimization algorithm and the signal reconstruction method. For brevity, the comparative algorithm is abbreviated as CA, and its flowchart is illustrated Figure 12.

Figure 12.

Flow chart of CA.

Vibration signal analysis results

Through computational, the MOA–VMD–BIC algorithm obtained the optimal VMD decomposition parameters as {k = 4, α = 4239}, while the CA algorithm derived the optimal VMD parameters as {k = 6, α = 3519}. The vibration signal components decomposed by the MOA–VMD–BIC and CA algorithms are sequentially illustrated in Figures 13 and 14.

Figure 13.

Vibration signal VMD decomposition results of MOA–VMD–BIC.

Figure 14.

Vibration signal VMD decomposition results of CA.

In the MOA–VMD–BIC algorithm, the covariance matrix R of the multidimensional signal matrix M is first calculated based on the BIC algorithm described in Section 2.3. Subsequently, the eigenvalues of R and the corresponding Bayesian objective cost are computed, with the results presented in Table 6. Analysis of Table 6 reveals that the minimum Bayesian objective cost value is 2. Therefore the IMF1 and IMF2 components of VMD are selected to reassemble the reconstructed signal rx ₁ (t). In the CA algorithm, the kurtosis values of each component are first calculated based on equation (15), as listed in Table 7. Subsequently, the IMF5 and IMF6 components with kurtosis values >3 are selected from Table 6 to reconstruct the signal rx ₂ (t).

Table 6.

Eigenvalue of matrix R and calculation results of BIC.

No.	Eigenvalues, λi	BIC (k)
1	0.003908296	5.48 × 108
2	0.001238243	1.15 × 108
3	0.009266714	4.49 × 1010
4	0.017691220	1.37 × 1014
5	0.103587328	9.34 × 1029

Table 7.

The kurtosis of the vibration signal vibration signal component.

Signal component	IMF1	IMF2	IMF3	IMF4	IMF5	IMF6
Kurtosis	2.175	2.522	2.881	2.571	3.691	3.764

Prior to envelope spectrum analysis of the vibration signals rx ₁ (t) and rx ₂ (t), a fourth-order high-pass filter was applied to eliminate the influence of DC components in the signal. The reconstructed signals rx ₁ (t) and rx ₂ (t) obtained through the MOA–VMD–BIC and CA algorithms, along with their envelope spectrum analyses, are shown in Figure 15.

Figure 15.

Envelope spectrum analysis by MOA–VMD–BIC, CA, and no algorithm.

As shown in Figure 15, the spectral analysis of the results from CA and no algorithm shows that both methods can extract low-order frequency components of rotor rotation and fault (1× and 2×). However, by using the MOA–VMD–BIC method, it is possible to obtain higher-order rotor rotation and fault frequency components compared to CA and no algorithm. The analysis results show that there is a deviations exist between the extracted frequency and their object frequency. The percentage deviation between the extracted and object frequency is calculated using equation (16), and the calculation results are summarized in the Table 8 below.

Table 8.

Analysis results of rotor frequency extraction.

	Rotor frequency (f r )
			No algorithm		MOA–VMD–BIC		CA
		Obj (Hz)	Extraction (Hz)	Error (%)	Extraction (Hz)	Error (%)	Extraction (Hz)	Error (%)
	1×	29.95	29	3.17	30.07	0.4007	31.45	5.0083
	2×	59.90	58	3.17	60.13	0.3840	62.91	5.0250
	3×	89.95	/	/	90.18	0.3673	94.36	5.0195

The term error in Table 8 denotes the relative error, which can be calculated using the following equation (16):

$$\mathrm{Error}=\left|\frac{f_{\mathrm{ex}}-f_{\mathrm{obj}}}{f_{\mathrm{obj}}}\right|\times 100\%$$

(16)

Where, f _ex represents the frequency value extracted from the reconstructed signal (rx ₁ (t) or rx ₂ (t)) using an algorithm (MOA–VMD–BIC or CA). f _obj represents the object frequency value of signal (rx ₁ (t) or rx ₂ (t)). Analysis of frequency extraction in vibration signals based on Table 7 calculation results: The proposed MOA–VMD–BIC method demonstrates superior applicability for spectrum analysis of vibration signals. Specifically, in terms of rotational frequency extraction for rotor systems, the MOA–VMD–BIC method maintains a relative error within 0.5%, whereas the CA method exhibits significantly larger deviations exceeding 5% in all scenarios.

As illustrated in the spectrum map shown in Figure 15, the fault frequencies are more readily identifiable in the results obtained from the MOA–VMD–BIC algorithm. In contrast, the spectrum generated by the CA algorithm exhibits challenges in fault frequency detection, as these frequencies may be obscured by neighboring frequency components. To systematically address this issue, we conducted a statistical analysis comparing the spectral amplitudes of fault characteristic frequency harmonics between the two algorithms. The statistical results of harmonic components and their corresponding amplitudes for MOA–VMD–BIC and CA algorithms are summarized in the Table 9 below.

Table 9.

Fault frequency and frequency amplitude analysis results.

	Fault frequency (f d )
			No algorithm		MOA–VMD–BIC		CA
		Obj (Hz)	Extraction (Hz)	Error (%)	Extraction (Hz)	Error (%)	Extraction (Hz)	Error (%)
	1×	107.36	99	7.78	107.4	0.04	107.4	0.04
	2×	214.72	216	0.6	214.8	0.04	214.8	0.04
	3×	322.08	/	/	322.2	0.04	/	/
	5×	536.80	/	/	538	0.22	/	/

As shown in Table 9, the relative errors between the harmonic frequencies identified by the two algorithms and the target harmonic frequencies range from 0.04% to 0.3%. The spectral analysis of the results from CA and no algorithm shows that both methods can extract low-order fault frequencies (1× and 2×). However, the fault frequencies obtained through CA are significantly more accurate than those obtained without using the algorithm. The spectral analysis of the results from MOA–VMD–BIC and CA reveals that the reconstructed signal rx ₂ (t) obtained by the CA algorithm still contains extraneous frequency components, which hinders effective identification of higher-order fault frequencies (e.g. 3× and 5× harmonics). In contrast, the MOA–VMD–BIC algorithm leverages the source number identification capability of the BIC algorithm to effectively filter out irrelevant frequency components, thereby facilitating the detection of higher harmonics such as 3× and 5× frequencies.

In conclusion, based on the analysis of reconstructed signals using MOA–VMD–BIC and CA algorithms, the identification of rotational frequencies and fault characteristic frequencies demonstrates that the MOA–VMD–BIC algorithm is more suitable for vibration signal-based fault feature extraction.

Conclusions

To address the challenges of multi-source vibration signals masking fault-induced components and low fault recognition rate in vibration signal analysis, a novel multi-objective optimization algorithm-based variational mode decomposition with Bayesian Information Criterion (MOA–VMD–BIC) method was proposed for bearing vibration signal analysis. The following conclusions were drawn:

1. In this paper, a multi-objective optimization model algorithm is proposed. This model comprises two sub-models (Model_A and Model_B), where Model_A serves as the foundation for iterative optimization of signal decomposition, while Model_B operates based on the computational results from Model_A. Each iteration of Model_A and Model_B generates a solution set. By normalizing the results of both models and using the parameter values {k, α} derived from Model_A as the baseline reference, we identify the intersection set of the two models, thereby obtaining the optimal solution for the multi-objective optimization model.

2. In the multi-objective optimization model algorithm, the characteristic parameter set of vibration signals (e.g. FFR, IO, IEC, envelope entropy, kurtosis) is utilized to construct the fitness function of the multi-objective optimization model.

3. Taking simulated signals as the research object, this study elaborates on the application of MOA–VMD–BIC in vibration signal fault identification. The simulation results demonstrate that MOA–VMD effectively decomposes signals and reconstructs vibration signals using the BIC-based source number identification capability. Subsequent envelope spectrum analysis of the reconstructed signals enables accurate identification of fault characteristics.

4. To validate the practical applicability of the proposed MOA–VMD–BIC algorithm, experimental verification was conducted using the Case Western Reserve University (CWRU) bearing dataset, with CA as the comparative algorithm. Throughout the analysis, identical optimization algorithms, data preprocessing methods (including DC component removal), and filters were employed. The key distinctions lie in the fitness function and signal reconstruction methodology of MOA–VMD–BIC compared to CA. From the identification results of rotor frequency, it can be seen that both MOA–VMD–BIC and CA algorithms can identify the 1×, 2×, and 3× frequency components of rotor rotation. However, deviations exist between the extracted frequency and their object frequency. Specifically, in terms of rotational frequency extraction for rotor systems, the MOA–VMD–BIC method maintains a relative error within 0.5%, whereas the CA method exhibits significantly larger deviations exceeding 5% in all scenarios. From the identification results of fault frequency, it can be seen that the relative errors between the harmonic frequencies identified by the two algorithms and the target harmonic frequencies range from 0.04% to 0.3%. However, spectral analysis of the results from MOA–VMD–BIC and CA reveals that the reconstructed signal rx₂(t) obtained by the CA algorithm still contains extraneous frequency components, which hinders effective identification of higher-order fault frequencies (e.g. 3× and 5× harmonics). In contrast, the MOA–VMD–BIC algorithm leverages the source number identification capability of the BIC algorithm to effectively filter out irrelevant frequency components, thereby facilitating the detection of higher harmonics such as 3× and 5× frequencies. In conclusion, based on the analysis of reconstructed signals using MOA–VMD–BIC and CA algorithms, the identification of rotational frequencies and fault characteristic frequencies demonstrates that the MOA–VMD–BIC algorithm is more suitable for vibration signal-based fault feature extraction.

In future research, we will further simplify the research results presented in this paper, apply them to industrial scenarios, and achieve experience accumulation and result transformation in this field. Meanwhile, we will strengthen the comparative study of algorithms related to VMD and EMD, and introduce the application research of emerging technologies such as machine learning in fault diagnosis.