# Blind Parameter Identification of MAR Model and Mutation Hybrid GWO-SCA Optimized SVM for Fault Diagnosis of Rotating Machinery

Authors: Wenlong Fu; Jiawen Tan; Xiaoyuan Zhang; Tie Chen; Kai Wang
Journal: Complexity
 (2019)

---

## Abstract

As a crucial and widely used component in industrial fields with great complexity, the health condition of rotating machinery is directly related to production efficiency and safety. Consequently, recognizing and diagnosing rotating machine faults remain to be one of the main concerns in preventing failures of mechanical systems, which can enhance the reliability and efficiency of mechanical systems. In this paper, a novel approach based on blind parameter identification of MAR model and mutation hybrid GWO-SCA optimization is proposed to diagnose faults for rotating machinery. Signals collected from different types of faults were firstly split into sets of intrinsic mode functions (IMFs) by variational mode decomposition (VMD), the decomposing mode numberK of which was preset with central frequency observation method. Then the multivariate autoregressive (MAR) model of all IMFs was established, whose order was determined by Schwartz Bayes Criterion (SBC), and all parameters of the model were identified blindly through QR decomposition, where key features were subsequently extracted via principal component analysis (PCA) to construct feature vectors of different fault types. Afterwards, a hybrid optimization algorithm combining mutation operator, grey wolf optimizer (GWO), and sine cosine algorithm (SCA), termed mutation hybrid GWO-SCA (MHGWOSCA), was proposed for parameter selection of support vector machine (SVM). The optimal SVM model was later employed to classify different fault samples. The engineering application and contrastive analysis indicate the availability and superiority of the proposed method.


---

## Body

## 1. Introduction

Rotating machinery is a crucial part in modern industrial manufacture, and its failure may result in serious safety accidents and economic losses. Therefore, one of the major topics to be investigated in preventing failures of mechanical systems is recognizing and diagnosing rotating machine faults [1, 2]. In various rotating machineries, rolling bearings are the most widely used components, which are greatly related to the remaining lifetime of machineries [3]. Unfortunately, rolling bearings are prone to failure to some extent because of complex operating conditions and structural characteristics [4]. Given this fact, effective and feasible fault diagnosis methods need to be developed for rolling bearings, to promote the reliability of rotating machinery [5].Owing to rich information carried by vibration signals, most of fault diagnosis methods for rolling bearings rely on analyzing vibration signals [6]. However, vibration signal is commonly nonstationary in actual operation, which restricts the extraction effect for fault features, and thus affects accuracy of fault diagnosis. Addressing this issue, many time-frequency signal analysis methods are proposed to analyze nonstationary signal, including empirical mode decomposition (EMD) [7], ensemble empirical mode decomposition (EEMD) [8], empirical wavelet transform (EWT) [9] and variational mode decomposition (VMD) [10]. EMD has advantages in dealing with nonstationary signal and does not need to preset any basis function, but its performance is affected by end effects and mode mixing. To overcome these defects, EEMD as an improved version of EMD is proposed by introducing noise-assisted analysis method, which is also increasing the computational cost and cannot completely neutralize the added noise. EWT is a self-adaptive signal decomposition method that inherits advantages of EMD and wavelet transform, but the selection of an appropriate wavelet filter bank remains to be a problem. Unlike the methods mentioned above, VMD is a quasi-orthogonal signal decomposition method and its main structure is to construct and solve a constrained variational problem. It can effectively realize the separation of each frequency component in signal, avoiding the problem of mode mixing or filter selection. Its advancement and effectiveness have been demonstrated in many previous studies [11, 12]. As the nonstationarity of signal weakened by above signal processing methods, fault features can be smoothly extracted. Due to the high efficiency and consideration of variable correlation, multivariate autoregressive (MAR) model [13] is introduced to extract features from vibration signal, which has been widely adopted for feature extraction and shows superior performance [14, 15]. Further, key components among the features extracted from MAR model are screened out via principal component analysis (PCA) to finally construct feature vectors [16].The essence of fault diagnosis for rolling bearings is a pattern recognition problem; that is, the fault type is determined by fault features of samples. In engineering application, many alternative methods are available for addressing this issue. For instance, Bayesian decision has strong recognition ability by considering prior probability as well as class conditional probability, and its good accuracy requires assumption of a suitable prior model [17]. K-nearest neighbor (KNN) based on Euclidean or Manhattan distances is easy to realize and susceptible to distribution of samples [18]. Artificial neural network (ANN) has good ability in processing recognition problem with the support of a large number of samples [19]. Based on statistical learning and structural risk minimization theory, support vector machine (SVM) is an excellent machine learning method. SVM finds an optimal hyperplane that meets the classification requirements, where the hyperplane maximizes the interval between two classes, thus promoting classification performance. Currently, SVM combined with feature extraction method has been successfully employed in the field of pattern recognition [20, 21].Although SVM has superior ability in pattern recognition, its performance is affected by parameters. In view of this, various optimization algorithms are developed and applied to search the best parameters, such as genetic algorithm (GA) [22], particle swarm optimization (PSO) [23], sine cosine algorithm (SCA) [24], and grey wolf optimizer (GWO) [25]. As novel optimization algorithms, SCA and GWO proposed by Mirjalili et al. have proved effective in many previous works [26–28]. To achieve better performance in convergence precision, a mutation hybrid GWO-SCA (MHGWOSCA) optimization method is proposed in this paper. To sum up, based on blind parameter identification of MAR model and mutation hybrid GWO-SCA optimized SVM, a fault diagnosis method for rotating machinery is presented in the research. Firstly, VMD was applied to decompose the vibration signal into a set of intrinsic mode functions (IMFs), and the decomposing mode numberK was preset with central frequency observation method. Then the MAR model of all IMFs was established, whose order was ascertained by Schwartz Bayes Criterion (SBC). And all parameters within the model were identified blindly through QR decomposition, where the key features were subsequently extracted via PCA to construct feature vectors of different fault samples. Afterwards, a hybrid optimization algorithm combining mutation operator, GWO, and SCA theory, called mutation hybrid GWO-SCA (MHGWOSCA), was proposed for parameter selection of SVM. Finally, the optimal SVM model was employed to classify different fault samples. The engineering application and contrastive analysis reveal the availability and superiority of the proposed method.The paper is organized as follows: Section2 is dedicated to the basic knowledge of VMD, MAR, and SVM. Section 3 introduces the proposed fault diagnosis method based on blind parameter identification of MAR model and mutation hybrid GWO-SCA optimized SVM. Section 4 presents the engineering application, the result analysis of which demonstrates the availability of the proposed method. The conclusion is summarized in Section 5.


## 2. Fundamental Theories

### 2.1. Variational Mode Decomposition

Variational Mode Decomposition (VMD) is a new adaptive signal preprocessing method [29]. By setting a mode numberK in advance, the given signal can be decomposed intoK band-limited intrinsic mode functions (IMFs). The cyclic iterative method is used to obtain the optimal solution of constrained variational problem, through which the frequency center and bandwidth of each mode are determined. The constrained variation problem can be described as follows:(1)minmk,wk∑k∂tδt+jπt∗mkte-jwkt22s.t.∑k=1Kmkt=ft,k=1,2,…,Kwhere mk={m1,m2,…,mk} represents the set of K mode functions and wk={w1,w2,…,wk} represents the set of central frequencies, while ∂t and δ(t) are the partial derivative of timet for the function and unit pulse function, respectively. f(t) is the given real valued input signal.For the above variational problem, quadratic penalty function term and Lagrange multiplier are employed to transform it into an unconstrained problem. Then problem (1) can be specified as follows:(2)Lmk,wk,β=α∑k∂tδt+jπt∗mkte-jwkt22+ft-∑kmkt22+βt,ft-∑kmktwhere α and β(t) represent penalty factor and Lagrange multiplier, respectively, andf(t) is the given signal.Then alternate direction method is utilized to obtain the saddle point of Lagrange multiplier, that is, solving (2) by optimizingmk,wk, and β(t) alternately [30]. The optimization problems ofmk andwk are formulated as (3) and (4), respectively.(3)mkn+1t=min⁡α∂tδt+jπt∗mkte-jwkt22+ft-∑imit+βt222(4)wkn+1=min⁡∂tδt+jπt∗mkte-jwkt22Solving problems (3) and (4), the iterative equations are deduced as follows.(5)m^kn+1w=f^w-∑i≠km^iw+β^w/21+2αw-wk2(6)wkn+1=∫0∞wm^kw2dw∫0∞m^kw2dwThe Lagrange multipliers can be iterated with(7)β^n+1w=β^nw+τfw-∑km^kn+1wwhereτ is an updating parameter.The main steps of VMD can be summarized as follows:Step 1.
Initializemk1,wk1,β1,n=1.Step 2 :.
Execute loop,n=n+1.Step 3.
Updatemk andwk based on (5) and (6).Step 4.
Updateβ based on (7).Step 5.
If∑km^kn+1-m^kn22/m^kn22<ε, loop end, else turn to Step 2 for next iteration.


### 2.2. Multivariate Autoregressive Model

Autoregressive (AR) model deduces the regressive variables with itself, namely, using a linear combination of white noise at timen and values of previousp times to describe the output of system at timen. The equation can be described as follows [31]:(8)xn=α1xn-1+α2xn-2+⋯+αpxn-p+enwhere αii=1,2,…,p is the autoregressive coefficient, en is the white noise with mean 0 and varianceσ2, andp is the order of model AR(p).When multivariate modeling is required, AR model must be analyzed and established multiple times. For this reason, the model order determined by different variables may be inconsistent. Compared with AR model, multivariate autoregressive (MAR) model only models once and can reflect the internal relationship among different variables. Letxnk(k=1,2,…,m) be the value of thek-th variable at timen; obviously it relates to not only the values of previousp times and the white noise at timen, but also the values of previousp times of other variables [13]. The equation of MAR model can be described as(9)xn1=∑i=1pα1i1xn-i1+∑i=1pα1i2xn-i2+⋯∑i=1pα1imxn-im+en1xn2=∑i=1pα2i1xn-i1+∑i=1pα2i2xn-i2+⋯∑i=1pα2imxn-im+en2⋮xnm=∑i=1pαmi1xn-i1+∑i=1pαmi2xn-i2+⋯∑i=1pαmimxn-im+enmwherep is the order of MAR model, αkim is the autoregressive coefficient for thek-th variable when them-th variable is delayed with stepi, and enm is the random error of them-th variable.


### 2.3. Support Vector Machine

Support vector machine (SVM) is a machine learning model developed by Vapnik [32]. It can deduce the optimal solution between model complexity and learning ability with limited information, so as to obtain the best classification accuracy. SVM has unique advantages in solving problems of limited samples, nonlinear and high-dimensional data. For a given sample set xi,yi∣i=1,2,…,n, the main idea of SVM is to map samples into a high-dimensional feature space and find a hyperplane to solve the corresponding nonlinear problem in low-dimensional space. The hyperplane function can be constructed as(10)w·x+b=0wherew is weight vector andb is bias parameter, while w∙x represents the inner product ofw andx.Taking binary classification issue as an example, to classify samples correctly with sufficient classification interval, all samples are required to meet the following constraint.(11)w·xi+b>1foryi=1<-1foryi=-1In the sequel, classification interval can be calculated as2/w2, and maximizing the interval is equivalent to minimizing w2. To solve the linear indivisibility problem, slack termξ and penalty factorC are brought into (11), and hence construction of the optimal hyperplane is transformed into a constrained optimization problem.(12)minf=12w2+C∑i=1nξis.t.yiw·xi+b≥1-ξi,i=1,2,…,nWhen mapping samples to high-dimensional space, different inner product kernel functions will form different algorithms. Among them, radial basis function is widely employed in application of pattern recognition, which can be described as(13)Kxi,xj=ϕxi·ϕxj=exp⁡-gxi-xj2where g is kernel parameter.The solution of constrained optimization problem (12) is determined by saddle point of Lagrange function, and this quadratic programming problem is transformed into the corresponding dual problem; that is,(14)maxL=∑i=1nμi-12∑i,j=1nμiμjyiyjKxi,xjs.t.∑i=1nμiyi=0,μi≥0,i=1,2,…,nwhere μi is Lagrange multiplier.With the optimal value ofμi obtained from above dual problem, the optimal classification discriminant function can be ascertained as follows.(15)fx=sgn∑i=1nμiKxi,x+b


## 2.1. Variational Mode Decomposition

Variational Mode Decomposition (VMD) is a new adaptive signal preprocessing method [29]. By setting a mode numberK in advance, the given signal can be decomposed intoK band-limited intrinsic mode functions (IMFs). The cyclic iterative method is used to obtain the optimal solution of constrained variational problem, through which the frequency center and bandwidth of each mode are determined. The constrained variation problem can be described as follows:(1)minmk,wk∑k∂tδt+jπt∗mkte-jwkt22s.t.∑k=1Kmkt=ft,k=1,2,…,Kwhere mk={m1,m2,…,mk} represents the set of K mode functions and wk={w1,w2,…,wk} represents the set of central frequencies, while ∂t and δ(t) are the partial derivative of timet for the function and unit pulse function, respectively. f(t) is the given real valued input signal.For the above variational problem, quadratic penalty function term and Lagrange multiplier are employed to transform it into an unconstrained problem. Then problem (1) can be specified as follows:(2)Lmk,wk,β=α∑k∂tδt+jπt∗mkte-jwkt22+ft-∑kmkt22+βt,ft-∑kmktwhere α and β(t) represent penalty factor and Lagrange multiplier, respectively, andf(t) is the given signal.Then alternate direction method is utilized to obtain the saddle point of Lagrange multiplier, that is, solving (2) by optimizingmk,wk, and β(t) alternately [30]. The optimization problems ofmk andwk are formulated as (3) and (4), respectively.(3)mkn+1t=min⁡α∂tδt+jπt∗mkte-jwkt22+ft-∑imit+βt222(4)wkn+1=min⁡∂tδt+jπt∗mkte-jwkt22Solving problems (3) and (4), the iterative equations are deduced as follows.(5)m^kn+1w=f^w-∑i≠km^iw+β^w/21+2αw-wk2(6)wkn+1=∫0∞wm^kw2dw∫0∞m^kw2dwThe Lagrange multipliers can be iterated with(7)β^n+1w=β^nw+τfw-∑km^kn+1wwhereτ is an updating parameter.The main steps of VMD can be summarized as follows:Step 1.
Initializemk1,wk1,β1,n=1.Step 2 :.
Execute loop,n=n+1.Step 3.
Updatemk andwk based on (5) and (6).Step 4.
Updateβ based on (7).Step 5.
If∑km^kn+1-m^kn22/m^kn22<ε, loop end, else turn to Step 2 for next iteration.


## 2.2. Multivariate Autoregressive Model

Autoregressive (AR) model deduces the regressive variables with itself, namely, using a linear combination of white noise at timen and values of previousp times to describe the output of system at timen. The equation can be described as follows [31]:(8)xn=α1xn-1+α2xn-2+⋯+αpxn-p+enwhere αii=1,2,…,p is the autoregressive coefficient, en is the white noise with mean 0 and varianceσ2, andp is the order of model AR(p).When multivariate modeling is required, AR model must be analyzed and established multiple times. For this reason, the model order determined by different variables may be inconsistent. Compared with AR model, multivariate autoregressive (MAR) model only models once and can reflect the internal relationship among different variables. Letxnk(k=1,2,…,m) be the value of thek-th variable at timen; obviously it relates to not only the values of previousp times and the white noise at timen, but also the values of previousp times of other variables [13]. The equation of MAR model can be described as(9)xn1=∑i=1pα1i1xn-i1+∑i=1pα1i2xn-i2+⋯∑i=1pα1imxn-im+en1xn2=∑i=1pα2i1xn-i1+∑i=1pα2i2xn-i2+⋯∑i=1pα2imxn-im+en2⋮xnm=∑i=1pαmi1xn-i1+∑i=1pαmi2xn-i2+⋯∑i=1pαmimxn-im+enmwherep is the order of MAR model, αkim is the autoregressive coefficient for thek-th variable when them-th variable is delayed with stepi, and enm is the random error of them-th variable.


## 2.3. Support Vector Machine

Support vector machine (SVM) is a machine learning model developed by Vapnik [32]. It can deduce the optimal solution between model complexity and learning ability with limited information, so as to obtain the best classification accuracy. SVM has unique advantages in solving problems of limited samples, nonlinear and high-dimensional data. For a given sample set xi,yi∣i=1,2,…,n, the main idea of SVM is to map samples into a high-dimensional feature space and find a hyperplane to solve the corresponding nonlinear problem in low-dimensional space. The hyperplane function can be constructed as(10)w·x+b=0wherew is weight vector andb is bias parameter, while w∙x represents the inner product ofw andx.Taking binary classification issue as an example, to classify samples correctly with sufficient classification interval, all samples are required to meet the following constraint.(11)w·xi+b>1foryi=1<-1foryi=-1In the sequel, classification interval can be calculated as2/w2, and maximizing the interval is equivalent to minimizing w2. To solve the linear indivisibility problem, slack termξ and penalty factorC are brought into (11), and hence construction of the optimal hyperplane is transformed into a constrained optimization problem.(12)minf=12w2+C∑i=1nξis.t.yiw·xi+b≥1-ξi,i=1,2,…,nWhen mapping samples to high-dimensional space, different inner product kernel functions will form different algorithms. Among them, radial basis function is widely employed in application of pattern recognition, which can be described as(13)Kxi,xj=ϕxi·ϕxj=exp⁡-gxi-xj2where g is kernel parameter.The solution of constrained optimization problem (12) is determined by saddle point of Lagrange function, and this quadratic programming problem is transformed into the corresponding dual problem; that is,(14)maxL=∑i=1nμi-12∑i,j=1nμiμjyiyjKxi,xjs.t.∑i=1nμiyi=0,μi≥0,i=1,2,…,nwhere μi is Lagrange multiplier.With the optimal value ofμi obtained from above dual problem, the optimal classification discriminant function can be ascertained as follows.(15)fx=sgn∑i=1nμiKxi,x+b


## 3. Fault Diagnosis by Blind Parameter Identification of MAR Model and Mutation Hybrid GWO-SCA Optimized SVM

### 3.1. Blind Parameter Identification of MAR Model

Multivariate autoregressive (MAR) model is the multivariate version of autoregressive (AR) model. By synchronous autoregressive modeling and parameter identification of multiple variables, it possesses advantages of high efficiency and consideration of variable correlation. Parameter identification is a critical part of MAR analysis and appropriate parameters will contribute to a better fitting effect, during which process the order selection of MAR model not only relates to the extraction of essential features implicated in the sequence, but also has a great impact on recognition efficiency. In this research, Schwartz Bayes Criterion (SBC) is adopted to estimate the order of MAR model. The SBC function is as follows [33]:(16)Sp=lpm-1-npNln⁡Nwherenp=mp+1, lp=ln⁡det⁡dp(R22TR22), andR22 is the fourth part of the upper triangular matrix obtained by QR decomposition of K matrix of the model (see [33] for more details).The specific method of order selection is that, within the given range of model order,S(p) will gradually decrease with the order increasing. IfS(p) takes the minimum value or there is no significant decreasing inS(p) as the order increases, the correspondingp value is the model order.Considering the inner relation among different variable sequences, there will be a mass of parameters in MAR model. Generally, minimum variance criterion is used to estimate model parameters directly [34], but it would be computationally expensive for MAR model. To achieve better computational efficiency, a fast calculation method based on QR decomposition is applied to identify the parameters (see [33] for more details). Then all the identified parameters are combined to construct the feature vector of the multivariate sequences.


### 3.2. Mutation GWO-SCA Optimization

#### 3.2.1. Sine Cosine Algorithm

Sine cosine algorithm (SCA) is a newly swarm optimization algorithm, which is simple in structure and easy to realize. Using only sine and cosine functions, SCA achieves good optimization effect. Exploration and exploitation are the two key phases when SCA processes optimization problem [24]. With an elaborate mechanism of exploration and exploitation, SCA is able to search for feasible solutions quickly in the searching space. Let Zi=(Zi1,Zi2,…,Zid)T be the position ofi-th (i=1,2,…,m) individual, whereupon each solution of the optimization problem corresponds to a position of corresponding individual in the searching space, whered is the dimension of individuals. The core updating equation of SCA is as follows [24]:(17)Zil+1=Zil+r1×sin⁡r2×r3Pil-Zil,r4<0.5Zil+r1×cos⁡r2×r3Pil-Zil,r4≥0.5where Pi=(Pi1,Pi2,…,Pid)T is the best position of individuali. The parameter r1=a-l(a/lmax) determines the region of individualiat next iteration, where lmax,l, anda are the maximum number of iterations, the current number of iterations, and a constant, respectively. The parameterr2 is a random number in the scope of [0,2π], which defines the distance that next movement should be towards. To stochastically emphasize (r3>1) or deemphasize (r3<1) the effect of the best position of individual, a random weightr3 within [0,2] is brought in the equation. Finally, the parameterr4 is a random number in the range of [0,1] to switch fairly between sine and cosine components.


#### 3.2.2. Grey Wolf Optimizer

Grey wolf optimizer is a swarm intelligence algorithm based on strict organization system of grey wolves and its sophisticated cooperative hunting method. GWO algorithm simulates the hierarchy and hunting behavior of grey wolves in the wild. Wolves are divided into four subgroups, i.e.,α, β,δ, andω. The first three subgroups are in turn the three wolves with the best fitness values, which guide other wolves (ω) to search for the target. During the optimization process, the wolves update the positions of α, β,δ, andω. The updating equations of grey wolves can be described as [25](18)Distance  updating:Dα=C1·Xαl-XlDβ=C2·Xβl-XlDδ=C3·Xδl-Xl(19)Position  updating:X1l+1=Xαl-A1·DαX2l+1=Xβl-A2·DβX3l+1=Xδl-A3·Dδ(20)Xl+1=X1+X2+X33where Dα, Dβ, and Dδ are distances between the position of individual α, β, andδand an individual, respectively, whileC1,C2, andC3 are random numbers in [0,2].A1=2a∙r-a,A2=2a∙r-a, andA3=2a∙r-a, wherea=2-2(l/lmax) is an adaptive parameter which will decrease linearly to 0 with iteration andr is a real random number in [0,1]. Xα, Xβ, and Xδ represent the current position of individual α, β, andδ, respectively.X indicates the final position of an individual which would be in a random place within a circle defined by individual α, β, andδ. Finally,l is the current number of iterations.


#### 3.2.3. Mutation Hybrid GWO-SCA Optimization

In the updating strategy of GWO, the position ofα is of the most importance, which guides the search behavior of other individuals. Besides in SCA, sine and cosine functions in (17) possess outstanding ability of balancing exploration and exploitation process. Hence the combination of capabilities between GWO and SCA is expected to improve the convergence performance, which makes the hybrid algorithm, grey wolf optimizer-sine cosine algorithm (GWO-SCA), powerful in both exploring different regions and exploiting promising regions. The part of individual α in (18) is improved with the updating strategy of SCA [35]:(21)Dαl+1=r1×sin⁡r2×C1Xαl-Xl,r4<0.5r1×cos⁡r2×C1Xαl-Xl,r4≥0.5where values of parameterr1,r2, andr4 are the same as in SCA, and the meaning of Dα,C1, Xα, andX are the same as in GWO, whilel is the current time of iterations. Except for the subsequent mutation operation, the rest of the processing procedure of GWO remains unchanged.Analysis of GWO shows that diversity of individuals is gradually losing with iteration and it is likely to result in the local optimum in the later stage of iteration, which will affect the optimization effect. To solve this problem, mutation operators are introduced to enrich diversity of individuals [36]. The new equation containing the mutation operators is applied to enhance (19):(22)X1l+1=Xαl-A1Dαl,l≠nTX1l1+G0.5-unif,l=nTX2l+1=Xβl-A2Dβl,l≠nTX2l1+G0.5-unif,l=nTX3l+1=Xδl-A3Dδl,l≠nTX3l1+G0.5-unif,l=nTn=1,2,…whereG andT (T<lmax) are the given mutation amplitude and mutation period, respectively, whileunif represents a random number that conforms to the uniform distributionU(0,1). It can be observed that, during the updating process, the position of individual α, β, andδwill mutate periodically, which gives them a certain ability to jump out of local optimum even in the later stage of iterations.The main steps of proposed mutation hybrid GWO-SCA (MHGWOSCA) are shown below:Step 1.
Initialize the population randomly in the given range of variables and set mutation parameters.Step 2.
Calculate the fitness value of each individual.Step 3.
Choose and save the best three individuals asα, β, andδ.Step 4.
Update other individuals (ω) according to (21) and (22).Step 5.
Update parametersa,A, andC.Step 6.
If the maximum number of iterations is not met, go to Step 2.Step 7.
Output the position of individualα as the optimal solution.


### 3.3. Diagnosis Method Based on Blind Parameter Identification of MAR Model and Mutation Hybrid GWO-SCA Optimization

In this section, to promote the fault diagnosis accuracy, a fault diagnosis method based on blind parameter identification of MAR model and mutation hybrid GWO-SCA optimized SVM is proposed. The specific steps are detailed as follows:Step 1.
Collect the vibration signals.Step 2.
Determine the mode numberK of VMD by central frequency observation method.Step 3.
Decompose the samples into sets of IMFs with VMD.Step 4.
Establish the MAR model of all IMFs and apply SBC to determine the model order; then MAR parameter vectors of different fault samples can be identified by QR decomposition.Step 5.
Extract principal features from parameter vectors based on PCA; thus the feature vectors of different fault types are obtained.Step 6.
Optimize the parametersC and g for SVM with the proposed MHGWOSCA optimization method.Step 7.
Train the SVM model with optimal parametersC and g; then the optimal SVM model is obtained.Step 8.
Apply the optimal SVM model to classify different fault samples and evaluate the performance of the model.The flowchart of the proposed fault diagnosis model is shown in Figure1.Figure 1
The flowchart of the proposed fault diagnosis method.


## 3.1. Blind Parameter Identification of MAR Model

Multivariate autoregressive (MAR) model is the multivariate version of autoregressive (AR) model. By synchronous autoregressive modeling and parameter identification of multiple variables, it possesses advantages of high efficiency and consideration of variable correlation. Parameter identification is a critical part of MAR analysis and appropriate parameters will contribute to a better fitting effect, during which process the order selection of MAR model not only relates to the extraction of essential features implicated in the sequence, but also has a great impact on recognition efficiency. In this research, Schwartz Bayes Criterion (SBC) is adopted to estimate the order of MAR model. The SBC function is as follows [33]:(16)Sp=lpm-1-npNln⁡Nwherenp=mp+1, lp=ln⁡det⁡dp(R22TR22), andR22 is the fourth part of the upper triangular matrix obtained by QR decomposition of K matrix of the model (see [33] for more details).The specific method of order selection is that, within the given range of model order,S(p) will gradually decrease with the order increasing. IfS(p) takes the minimum value or there is no significant decreasing inS(p) as the order increases, the correspondingp value is the model order.Considering the inner relation among different variable sequences, there will be a mass of parameters in MAR model. Generally, minimum variance criterion is used to estimate model parameters directly [34], but it would be computationally expensive for MAR model. To achieve better computational efficiency, a fast calculation method based on QR decomposition is applied to identify the parameters (see [33] for more details). Then all the identified parameters are combined to construct the feature vector of the multivariate sequences.


## 3.2. Mutation GWO-SCA Optimization

### 3.2.1. Sine Cosine Algorithm

Sine cosine algorithm (SCA) is a newly swarm optimization algorithm, which is simple in structure and easy to realize. Using only sine and cosine functions, SCA achieves good optimization effect. Exploration and exploitation are the two key phases when SCA processes optimization problem [24]. With an elaborate mechanism of exploration and exploitation, SCA is able to search for feasible solutions quickly in the searching space. Let Zi=(Zi1,Zi2,…,Zid)T be the position ofi-th (i=1,2,…,m) individual, whereupon each solution of the optimization problem corresponds to a position of corresponding individual in the searching space, whered is the dimension of individuals. The core updating equation of SCA is as follows [24]:(17)Zil+1=Zil+r1×sin⁡r2×r3Pil-Zil,r4<0.5Zil+r1×cos⁡r2×r3Pil-Zil,r4≥0.5where Pi=(Pi1,Pi2,…,Pid)T is the best position of individuali. The parameter r1=a-l(a/lmax) determines the region of individualiat next iteration, where lmax,l, anda are the maximum number of iterations, the current number of iterations, and a constant, respectively. The parameterr2 is a random number in the scope of [0,2π], which defines the distance that next movement should be towards. To stochastically emphasize (r3>1) or deemphasize (r3<1) the effect of the best position of individual, a random weightr3 within [0,2] is brought in the equation. Finally, the parameterr4 is a random number in the range of [0,1] to switch fairly between sine and cosine components.


### 3.2.2. Grey Wolf Optimizer

Grey wolf optimizer is a swarm intelligence algorithm based on strict organization system of grey wolves and its sophisticated cooperative hunting method. GWO algorithm simulates the hierarchy and hunting behavior of grey wolves in the wild. Wolves are divided into four subgroups, i.e.,α, β,δ, andω. The first three subgroups are in turn the three wolves with the best fitness values, which guide other wolves (ω) to search for the target. During the optimization process, the wolves update the positions of α, β,δ, andω. The updating equations of grey wolves can be described as [25](18)Distance  updating:Dα=C1·Xαl-XlDβ=C2·Xβl-XlDδ=C3·Xδl-Xl(19)Position  updating:X1l+1=Xαl-A1·DαX2l+1=Xβl-A2·DβX3l+1=Xδl-A3·Dδ(20)Xl+1=X1+X2+X33where Dα, Dβ, and Dδ are distances between the position of individual α, β, andδand an individual, respectively, whileC1,C2, andC3 are random numbers in [0,2].A1=2a∙r-a,A2=2a∙r-a, andA3=2a∙r-a, wherea=2-2(l/lmax) is an adaptive parameter which will decrease linearly to 0 with iteration andr is a real random number in [0,1]. Xα, Xβ, and Xδ represent the current position of individual α, β, andδ, respectively.X indicates the final position of an individual which would be in a random place within a circle defined by individual α, β, andδ. Finally,l is the current number of iterations.


### 3.2.3. Mutation Hybrid GWO-SCA Optimization

In the updating strategy of GWO, the position ofα is of the most importance, which guides the search behavior of other individuals. Besides in SCA, sine and cosine functions in (17) possess outstanding ability of balancing exploration and exploitation process. Hence the combination of capabilities between GWO and SCA is expected to improve the convergence performance, which makes the hybrid algorithm, grey wolf optimizer-sine cosine algorithm (GWO-SCA), powerful in both exploring different regions and exploiting promising regions. The part of individual α in (18) is improved with the updating strategy of SCA [35]:(21)Dαl+1=r1×sin⁡r2×C1Xαl-Xl,r4<0.5r1×cos⁡r2×C1Xαl-Xl,r4≥0.5where values of parameterr1,r2, andr4 are the same as in SCA, and the meaning of Dα,C1, Xα, andX are the same as in GWO, whilel is the current time of iterations. Except for the subsequent mutation operation, the rest of the processing procedure of GWO remains unchanged.Analysis of GWO shows that diversity of individuals is gradually losing with iteration and it is likely to result in the local optimum in the later stage of iteration, which will affect the optimization effect. To solve this problem, mutation operators are introduced to enrich diversity of individuals [36]. The new equation containing the mutation operators is applied to enhance (19):(22)X1l+1=Xαl-A1Dαl,l≠nTX1l1+G0.5-unif,l=nTX2l+1=Xβl-A2Dβl,l≠nTX2l1+G0.5-unif,l=nTX3l+1=Xδl-A3Dδl,l≠nTX3l1+G0.5-unif,l=nTn=1,2,…whereG andT (T<lmax) are the given mutation amplitude and mutation period, respectively, whileunif represents a random number that conforms to the uniform distributionU(0,1). It can be observed that, during the updating process, the position of individual α, β, andδwill mutate periodically, which gives them a certain ability to jump out of local optimum even in the later stage of iterations.The main steps of proposed mutation hybrid GWO-SCA (MHGWOSCA) are shown below:Step 1.
Initialize the population randomly in the given range of variables and set mutation parameters.Step 2.
Calculate the fitness value of each individual.Step 3.
Choose and save the best three individuals asα, β, andδ.Step 4.
Update other individuals (ω) according to (21) and (22).Step 5.
Update parametersa,A, andC.Step 6.
If the maximum number of iterations is not met, go to Step 2.Step 7.
Output the position of individualα as the optimal solution.


## 3.2.1. Sine Cosine Algorithm

Sine cosine algorithm (SCA) is a newly swarm optimization algorithm, which is simple in structure and easy to realize. Using only sine and cosine functions, SCA achieves good optimization effect. Exploration and exploitation are the two key phases when SCA processes optimization problem [24]. With an elaborate mechanism of exploration and exploitation, SCA is able to search for feasible solutions quickly in the searching space. Let Zi=(Zi1,Zi2,…,Zid)T be the position ofi-th (i=1,2,…,m) individual, whereupon each solution of the optimization problem corresponds to a position of corresponding individual in the searching space, whered is the dimension of individuals. The core updating equation of SCA is as follows [24]:(17)Zil+1=Zil+r1×sin⁡r2×r3Pil-Zil,r4<0.5Zil+r1×cos⁡r2×r3Pil-Zil,r4≥0.5where Pi=(Pi1,Pi2,…,Pid)T is the best position of individuali. The parameter r1=a-l(a/lmax) determines the region of individualiat next iteration, where lmax,l, anda are the maximum number of iterations, the current number of iterations, and a constant, respectively. The parameterr2 is a random number in the scope of [0,2π], which defines the distance that next movement should be towards. To stochastically emphasize (r3>1) or deemphasize (r3<1) the effect of the best position of individual, a random weightr3 within [0,2] is brought in the equation. Finally, the parameterr4 is a random number in the range of [0,1] to switch fairly between sine and cosine components.


## 3.2.2. Grey Wolf Optimizer

Grey wolf optimizer is a swarm intelligence algorithm based on strict organization system of grey wolves and its sophisticated cooperative hunting method. GWO algorithm simulates the hierarchy and hunting behavior of grey wolves in the wild. Wolves are divided into four subgroups, i.e.,α, β,δ, andω. The first three subgroups are in turn the three wolves with the best fitness values, which guide other wolves (ω) to search for the target. During the optimization process, the wolves update the positions of α, β,δ, andω. The updating equations of grey wolves can be described as [25](18)Distance  updating:Dα=C1·Xαl-XlDβ=C2·Xβl-XlDδ=C3·Xδl-Xl(19)Position  updating:X1l+1=Xαl-A1·DαX2l+1=Xβl-A2·DβX3l+1=Xδl-A3·Dδ(20)Xl+1=X1+X2+X33where Dα, Dβ, and Dδ are distances between the position of individual α, β, andδand an individual, respectively, whileC1,C2, andC3 are random numbers in [0,2].A1=2a∙r-a,A2=2a∙r-a, andA3=2a∙r-a, wherea=2-2(l/lmax) is an adaptive parameter which will decrease linearly to 0 with iteration andr is a real random number in [0,1]. Xα, Xβ, and Xδ represent the current position of individual α, β, andδ, respectively.X indicates the final position of an individual which would be in a random place within a circle defined by individual α, β, andδ. Finally,l is the current number of iterations.


## 3.2.3. Mutation Hybrid GWO-SCA Optimization

In the updating strategy of GWO, the position ofα is of the most importance, which guides the search behavior of other individuals. Besides in SCA, sine and cosine functions in (17) possess outstanding ability of balancing exploration and exploitation process. Hence the combination of capabilities between GWO and SCA is expected to improve the convergence performance, which makes the hybrid algorithm, grey wolf optimizer-sine cosine algorithm (GWO-SCA), powerful in both exploring different regions and exploiting promising regions. The part of individual α in (18) is improved with the updating strategy of SCA [35]:(21)Dαl+1=r1×sin⁡r2×C1Xαl-Xl,r4<0.5r1×cos⁡r2×C1Xαl-Xl,r4≥0.5where values of parameterr1,r2, andr4 are the same as in SCA, and the meaning of Dα,C1, Xα, andX are the same as in GWO, whilel is the current time of iterations. Except for the subsequent mutation operation, the rest of the processing procedure of GWO remains unchanged.Analysis of GWO shows that diversity of individuals is gradually losing with iteration and it is likely to result in the local optimum in the later stage of iteration, which will affect the optimization effect. To solve this problem, mutation operators are introduced to enrich diversity of individuals [36]. The new equation containing the mutation operators is applied to enhance (19):(22)X1l+1=Xαl-A1Dαl,l≠nTX1l1+G0.5-unif,l=nTX2l+1=Xβl-A2Dβl,l≠nTX2l1+G0.5-unif,l=nTX3l+1=Xδl-A3Dδl,l≠nTX3l1+G0.5-unif,l=nTn=1,2,…whereG andT (T<lmax) are the given mutation amplitude and mutation period, respectively, whileunif represents a random number that conforms to the uniform distributionU(0,1). It can be observed that, during the updating process, the position of individual α, β, andδwill mutate periodically, which gives them a certain ability to jump out of local optimum even in the later stage of iterations.The main steps of proposed mutation hybrid GWO-SCA (MHGWOSCA) are shown below:Step 1.
Initialize the population randomly in the given range of variables and set mutation parameters.Step 2.
Calculate the fitness value of each individual.Step 3.
Choose and save the best three individuals asα, β, andδ.Step 4.
Update other individuals (ω) according to (21) and (22).Step 5.
Update parametersa,A, andC.Step 6.
If the maximum number of iterations is not met, go to Step 2.Step 7.
Output the position of individualα as the optimal solution.


## 3.3. Diagnosis Method Based on Blind Parameter Identification of MAR Model and Mutation Hybrid GWO-SCA Optimization

In this section, to promote the fault diagnosis accuracy, a fault diagnosis method based on blind parameter identification of MAR model and mutation hybrid GWO-SCA optimized SVM is proposed. The specific steps are detailed as follows:Step 1.
Collect the vibration signals.Step 2.
Determine the mode numberK of VMD by central frequency observation method.Step 3.
Decompose the samples into sets of IMFs with VMD.Step 4.
Establish the MAR model of all IMFs and apply SBC to determine the model order; then MAR parameter vectors of different fault samples can be identified by QR decomposition.Step 5.
Extract principal features from parameter vectors based on PCA; thus the feature vectors of different fault types are obtained.Step 6.
Optimize the parametersC and g for SVM with the proposed MHGWOSCA optimization method.Step 7.
Train the SVM model with optimal parametersC and g; then the optimal SVM model is obtained.Step 8.
Apply the optimal SVM model to classify different fault samples and evaluate the performance of the model.The flowchart of the proposed fault diagnosis model is shown in Figure1.Figure 1
The flowchart of the proposed fault diagnosis method.


## 4. Engineering Application

### 4.1. Data Collection

Vibration signals with different fault locations and sizes gathered from Bearings Data Center of Case Western Reserve University [37] were employed as the experiment data in this research. The experiment stand was mainly composed of a motor, an accelerometer, a torque sensor/encoder, and a dynamometer. The bearing was a deep groove ball bearing with model SKF6205-2RS. Accelerometers were applied to collect vibration signal; meanwhile, the torque sensor/encoder was for the measurement of power and speed. By using electro-discharge machining, the experiment device simulated three fault states of the rolling bearing: inner race fault, ball element fault, and outer race fault, and the depth of faults was 0.011 inches. In the experiment, vibration signals collected from the drive end (DE) were taken as the research objects, where the sample frequency was 12kHz and the rotation speed was 1750 rpm under the rated load of 2 hp. To fully verify the validity of the proposed fault diagnosis method, 9 types of samples were used in this paper, namely, inner race fault, ball fault, and outer race fault with diameters of 0.007, 0.014, and 0.021 inches (i.e., each of the three types of faults has three defect sizes). Further, vibration signal of each fault type contained 59 samples and each sample contained 2048 sample points. Image of the experiment device is shown in Figure 2, and the experimental data [5] are listed in Table 1.Table 1
Description of the experimental data.
Position of fault Defect size (inches) Label of classes Number of samples Inner race 0.007 L1 59 Inner race 0.014 L2 59 Inner race 0.021 L3 59 Ball 0.007 L4 59 Ball 0.014 L5 59 Ball 0.021 L6 59 Outer race 0.007 L7 59 Outer race 0.014 L8 59 Outer race 0.021 L9 59Figure 2
Experiment stand and its measurement sensors in bearing data center.


### 4.2. Evaluation Metrics

In this application, adjusted Rand index (ARI), normalized mutual information (NMI), F-measure (F), and accuracy (ACC) as four widely used evaluation measures are employed to evaluate the diagnosis results [38], which reflect the matching degree between classification results and real sample distribution information.LetΦ be the given actual label set and Ω be the classified result set. ARI considers the number of samples with the same labels and different labels. The scope of ARI is [-1,1] and the higher value indicates that more samples are labeled correctly. It can be expressed as(23)ARI=n11+n00Cn2wheren11 is the number of sample pairs with the same label in bothΦ andΩ, whilen00 is the number of sample pairs with different labels in Φ andΩ. Cn2 is all possible combinations of sample pairs.NMI is one of the common external evaluation metrics that can measure the degree of agreement between two data distributions [39]. The value of NMI lies between 0 and 1, and the higher value indicates a better classification performance. The calculation equation is as follows:(24)NMI=∑φ∈Φ∑ω∈Ωpφ,ωlog⁡pφ,ω/pφpω∑φ∈Φpφlog⁡pφ∑ω∈Ωpωlog⁡pωwherep(φ,ω) is the joint probability function ofΦ andΩ, andp(φ) andp(ω) are the probability distribution functions ofΦ andΩ, respectively.F and ACC are derivations from the confusion matrix whose structure is reflected in Table2 [40, 41]. F comprehensively considers the proportion of true positive result in positive actual label and in positive classified result. ACC represents the proportion of all correctly classified results in total results. The scope of these two metrics is both [0,1], and the higher value indicates a better classification performance. F and ACC can be defined as (25) and (26), respectively.(25)F=2TP/TP+FPTP/TP+FNTP/TP+FP+TP/TP+FN(26)ACC=TP+TNTP+FN+FP+FNTable 2
The confusion matrix.
Actual label Positive Negative Classification result Positive True positive (TP) False positive (FP) Negative False negative (FN) True negative (TN)


### 4.3. Application to Fault Diagnosis of Rolling Bearing

To demonstrate the superiority of the proposed VMD-MAR-MHGWOSCA-SVM method, the experiment was carried out in a comparative form at each stage, as follows: EMD was applied for comparison in the signal processing stage; fuzzy entropy (FE) [42] and permutation entropy (PE) [5] were used in the feature extraction stage for comparison; GWO, SCA, and HGWOSCA were employed for comparison in the parameter optimization stage. The parameters of the contrastive experiments were set in the same way as the proposed method.In the proposed method, the first step is to decompose signal of each fault type into a set of IMFs. The mode numberK of VMD needs to be determined in advance with central frequency observation method. IfK is too large, mode mixing would occur because central frequencies of adjacent IMFs are too close. In contrast, it is difficult to weaken the nonstationarity effectively ifK is too small. Central frequencies under differentK value are listed in Table 3, whereK was ascertained by the sample of inner race fault with diameter of 0.007 inches (L1). As illustrated in Table 3, whenKwas set to 5, the first two central frequencies were relatively close, which meant excessive decomposition had occurred; hence parameterKwas set to 4. The original signal waveforms of different fault types are shown in Figure 3 and decomposition results of 0.007 inches faults (L1, L4, L7) are illustrated in Figure 4. It can be observed that waveforms from 9 types of faults are quite different and IMFs of one fault type are also discriminative to a certain extent.Table 3
Normalized central frequencies with differentK values.
Number of modes Normalized central frequencies 2 0.2256 0.0907 3 0.2978 0.2275 0.0999 4 0.2471 0.1738 0.0860 0.0127 5 0.3066 0.2834 0.2271 0.1148 0.0525 6 0.3072 0.2854 0.2364 0.2142 0.1139 0.0519 7 0.3170 0.2981 0.2757 0.2270 0.1242 0.0982 0.0493 8 0.3178 0.2990 0.2777 0.2363 0.2147 0.1231 0.0975 0.0492 9 0.4022 0.3133 0.2970 0.2768 0.2363 0.2146 0.1230 0.0975 0.0492Figure 3
Time domain waveforms of signals with different faults: (a) inner race 0.007 inches, (b) inner race 0.014 inches, (c) inner race 0.021 inches, (d) ball 0.007 inches, (e) ball 0.014 inches, (f) ball 0.021 inches, (g) outer race 0.007 inches, (h) outer race 0.014 inches, (i) outer race 0.021 inches.
(a) (b) (c) (d) (e) (f) (g) (h) (i)Figure 4
The VMD decomposition results of 0.007 inches fault signals: (a) inner race, (b) ball, (c) outer race.
(a) (b) (c)When decomposed components of a sample were obtained, a multivariate autoregressive (MAR) model was established for all IMFs. The values of SBC valueS(p) are shown in Figure 5, where the order of MAR model was determined according to the sample of fault type L1. The figure shows that as the model orderp increases,S(p) decreases gradually, but the value ofS(p) is no longer significantly decreased whenp increases to 10. Hence, the MAR model orderp is selected as 10 in the experiment. After establishing the MAR model, all parameters were identified based on the QR decomposition and thus the fault feature vectors were formed. Further, principal features were extracted from all features by using principal component analysis (PCA), thus achieving a balance between the computational cost and the integrity of feature information. The principal component (PC) distribution of the feature vectors obtained by blind parameter identification of MAR model is illustrated in Figure 6 and Table 4, while the 3D projection of feature vectors for different faults is shown in Figure 7.Table 4
Principal components of feature vectors.
Fault label Principal components PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13 L1 -0.465 0.498 -0.055 -0.516 0.095 0.008 -0.101 0.334 0.151 0.100 -0.040 0.078 -0.097 -0.459 0.510 -0.076 -0.515 0.116 0.003 -0.082 0.297 0.137 0.083 -0.032 0.076 -0.095  L2 -0.362 -0.252 0.257 -0.034 0.268 -0.861 -0.699 -0.181 -0.131 -0.702 0.053 0.071 -0.185 0.396 -0.547 -0.194 -0.347 0.161 -0.432 -0.415 0.160 0.059 -0.438 0.510 -0.070 -0.393  L3 0.883 0.374 1.127 0.456 0.168 0.145 -0.120 0.229 -0.036 0.061 0.054 -0.031 0.020 0.903 0.393 1.113 0.462 0.140 0.162 -0.158 0.221 -0.054 0.065 0.081 -0.054 0.030  L4 -0.724 -1.063 0.770 0.037 -0.593 0.167 0.025 -0.271 -0.099 0.045 0.053 0.124 -0.244 0.277 -1.649 0.483 -0.168 0.614 1.125 0.160 -0.182 0.285 -0.320 0.017 0.495 0.519  L5 -0.387 0.075 -0.232 -0.161 0.232 0.566 0.057 -0.033 -0.077 -0.069 0.124 0.013 0.072 -0.334 -0.002 -0.292 -0.107 0.242 0.600 0.082 0.003 -0.120 -0.102 0.127 0.039 0.082  L6 0.350 -0.174 0.128 -0.874 -0.099 -0.441 0.413 0.533 -0.149 0.306 -0.065 0.204 0.002 -1.335 0.256 -0.407 0.647 -0.262 0.114 0.339 0.128 -0.068 0.374 -0.047 -0.265 -0.284  L7 1.847 0.174 0.321 0.410 0.250 -0.255 -0.026 -0.182 -0.008 0.198 -0.283 0.083 -0.011 1.995 0.217 0.045 0.276 0.170 -0.364 -0.070 -0.297 -0.048 0.130 -0.086 0.061 0.005  L8 -1.110 0.912 0.024 1.937 -1.912 -0.797 0.348 0.720 -0.036 -0.714 0.657 0.399 0.280 -1.337 0.116 -0.700 0.553 0.894 -0.020 0.247 -0.459 -0.077 0.289 -0.036 0.308 0.335  L9 -0.845 0.547 -0.570 0.382 -0.046 0.345 0.295 0.237 0.015 -0.127 -0.100 -0.218 0.005 0.753 0.129 -0.638 0.636 0.267 0.380 0.208 -0.340 -0.274 -0.473 -0.328 -0.078 -0.162Figure 5
Selection of MAR model orderp.Figure 6
Principal component distribution of feature vectors.Figure 7
3D projection of feature vectors for different faults.During optimization stage, the feature vectors of different fault types were divided randomly into two parts with 30 and 29 vectors, respectively, and the 30 ones were used for training while the 29 ones for testing. The proposed MHGWOSCA approach was applied to optimize the penalty factorC and the kernel parameter g, where individual number and iteration time were set to 30 and 100, respectively, while searching ranges ofC and g were both [2-10,210]. Considering the maximum number of iterations and the number of optimization parameters, mutation amplitude was set to 0.5 and mutation period was set to 5. During optimization process, fitness values were calculated by fivefold cross validation for the training samples. Hence, the optimization effect of the givenC and g can be measured by average accuracy of cross validation, and then the optimalCand g were obtained. The convergence procedure of MHGWOSCA is depicted in Figure 8, from which it can be observed that, with the progress of iteration, the average fitness value of individuals rises rapidly at first and then tends to be stable; that is, the individuals are close to the global optimal solution. In addition, due to the effect of the mutation operators, the best three individuals α, β, andδchange greatly and periodically, trying to jump out of the local optimum without affecting the overall convergence trend. The shaded part shows the distribution of the convergence curve in 10 runs. The convergence curves of different optimization methods are compared in Figure 9, while the average fitness values for different optimization methods in the last 20 iterations are listed in Table 5. The convergence data are the average of 10 runs. The figure shows that the convergence rate of SCA is significantly lower than other methods; GWO and HGWOSCA both have good convergence trend; MHGWOSCA first converges to a high level in the initial stage and shows the best convergence effect at later stage of iterations among all methods.Table 5
The average fitness value of individuals during iterations.
Iteration number Optimization methods MHGWOSCA HGWOSCA GWO SCA 81 88.12 97.25 97.15 95.95 82 97.32 97.25 97.17 96.06 83 97.32 97.25 97.17 95.87 84 97.32 97.25 97.17 96.03 85 97.33 97.26 97.17 96.13 86 89.70 97.25 97.17 96.16 87 97.32 97.25 97.17 96.04 88 97.33 97.25 97.18 96.02 89 97.32 97.25 97.18 96.08 90 97.32 97.25 97.17 96.09 91 88.87 97.26 97.19 96.25 92 97.32 97.25 97.18 96.33 93 97.32 97.25 97.18 96.34 94 97.32 97.25 97.19 96.29 95 97.32 97.26 97.19 96.43 96 88.99 97.26 97.19 96.37 97 97.32 97.26 97.19 96.40 98 97.32 97.26 97.18 96.43 99 97.32 97.26 97.19 96.46 100 97.32 97.26 97.19 96.46Figure 8
Convergence procedure of proposed optimization approach: (a) overall trend, (b) partial enlarged detail.
(a) (b)Figure 9
Comparison of different optimization methods: (a) overall trend, (b) partial enlarged detail.
(a) (b)With the optimal parametersCand g, the SVM model optimized by MHGWOSCA was trained and applied to recognize testing samples. For an in-depth verification of the proposed method, diagnosis results were averaged over 10 runs and training samples were chosen at random in every run. Furthermore, deviation of each result was calculated as a reference for result evaluation. The optimalC and g are presented corresponding to the best accuracy among 10 runs.Fault diagnosis results for different methods are shown in Table6 and Figure 10. It can be seen from Table 6 that the evaluation values of the proposed VMD-MAR-MHGWOSCA-SVM method are the highest among all four metrics, i.e., 09573, 0.9554, 0.9808, and 0.9808, and the deviations of the evaluation values are also at a small level. Specifically, the comparison between EMD-MAR-MHGWOSCA-SVM and the proposed method indicates the effectiveness of VMD; the comparison of VMD-FE-HGWOSCA-SVM, VMD-PE-HGWOSCA-SVM, and the proposed method shows MAR model possess a better feature extraction ability; similarly, it can be concluded from the contrastive experiment among VMD-MAR-GWO-SVM, VMD-MAR-HGWOSCA-SVM, VMD-MAR-SCA-SVM, and the proposed method that MHGWOSCA improves the diagnosis performance of SVM. Figure 11 is the boxplots of evaluation values, illustrating the performance of different diagnosis methods. As shown in the figure, the proposed method achieves the highest evaluation values and possesses outstanding stability among all methods.Table 6
Fault diagnosis results with different methods.
Processing methods C g Result evaluation Adjusted Rand index (ARI) Normalized mutual information (NMI) F-measure (F) Accuracy (ACC) VMD-FE-HGWOSCA-SVM 184.0764 2.2262 0.9001 0.9046 0.9519 0.9525 [-0.034, 0.034] [-0.032, 0.034] [-0.019, 0.017] [-0.018, 0.017] VMD-PE-HGWOSCA-SVM 2.5155 414.0703 0.9216 0.9221 0.9635 0.9636 [-0.037, 0.053] [-0.043, 0.048] [-0.017, 0.025] [-0.017, 0.025] EMD-MAR-GWO-SVM 172.6222 0.0010 0.9284 0.9244 0.9662 0.9667 [-0.031, 0.038] [-0.033, 0.039] [-0.017, 0.018] [-0.016, 0.018] EMD-MAR-MHGWOSCA-SVM 254.3134 0.0339 0.9296 0.9236 0.9661 0.9667 [-0.072, 0.053] [-0.074, 0.055] [-0.041, 0.026] [-0.039, 0.026] VMD-MAR-GWO-SVM 515.8658 0.1028 0.9325 0.9311 0.9691 0.9693 [-0.048, 0.033] [-0.054, 0.039] [-0.024, 0.016] [-0.023, 0.015] VMD-MAR-HGWOSCA-SVM 324.7595 0.0925 0.9380 0.9366 0.9719 0.9720 [-0.039, 0.037] [-0.053, 0.037] [-0.018, 0.017] [-0.018, 0.016] VMD-MAR-SCA-SVM 212.6435 0.3115 0.9431 0.9401 0.9741 0.9743 [-0.028, 0.032] [-0.041, 0.032] [-0.012, 0.014] [-0.013, 0.014] VMD-MAR-MHGWOSCA-SVM 261.4197 0.1909 0.9573 0.9554 0.9808 0.9808 [-0.034, 0.025] [-0.037, 0.023] [-0.015, 0.012] [-0.015, 0.011]Figure 10
Comparison of evaluation values for different methods.Figure 11
Boxplots of evaluation values for different methods; the x-axis tick labels correspond to 1: VMD-FE-HGWOSCA-SVM, 2: VMD-PE-HGWOSCA-SVM, 3: EMD-MAR-GWO-SVM, 4: EMD-MAR-MHGWOSCA-SVM, 5: VMD-MAR-GWO-SVM, 6 VMD-MAR-HGWOSCA-SVM, 7: VMD-MAR-SCA-SVM, and 8: VMD-MAR-MHGWOSCA-SVM.


## 4.1. Data Collection

Vibration signals with different fault locations and sizes gathered from Bearings Data Center of Case Western Reserve University [37] were employed as the experiment data in this research. The experiment stand was mainly composed of a motor, an accelerometer, a torque sensor/encoder, and a dynamometer. The bearing was a deep groove ball bearing with model SKF6205-2RS. Accelerometers were applied to collect vibration signal; meanwhile, the torque sensor/encoder was for the measurement of power and speed. By using electro-discharge machining, the experiment device simulated three fault states of the rolling bearing: inner race fault, ball element fault, and outer race fault, and the depth of faults was 0.011 inches. In the experiment, vibration signals collected from the drive end (DE) were taken as the research objects, where the sample frequency was 12kHz and the rotation speed was 1750 rpm under the rated load of 2 hp. To fully verify the validity of the proposed fault diagnosis method, 9 types of samples were used in this paper, namely, inner race fault, ball fault, and outer race fault with diameters of 0.007, 0.014, and 0.021 inches (i.e., each of the three types of faults has three defect sizes). Further, vibration signal of each fault type contained 59 samples and each sample contained 2048 sample points. Image of the experiment device is shown in Figure 2, and the experimental data [5] are listed in Table 1.Table 1
Description of the experimental data.
Position of fault Defect size (inches) Label of classes Number of samples Inner race 0.007 L1 59 Inner race 0.014 L2 59 Inner race 0.021 L3 59 Ball 0.007 L4 59 Ball 0.014 L5 59 Ball 0.021 L6 59 Outer race 0.007 L7 59 Outer race 0.014 L8 59 Outer race 0.021 L9 59Figure 2
Experiment stand and its measurement sensors in bearing data center.


## 4.2. Evaluation Metrics

In this application, adjusted Rand index (ARI), normalized mutual information (NMI), F-measure (F), and accuracy (ACC) as four widely used evaluation measures are employed to evaluate the diagnosis results [38], which reflect the matching degree between classification results and real sample distribution information.LetΦ be the given actual label set and Ω be the classified result set. ARI considers the number of samples with the same labels and different labels. The scope of ARI is [-1,1] and the higher value indicates that more samples are labeled correctly. It can be expressed as(23)ARI=n11+n00Cn2wheren11 is the number of sample pairs with the same label in bothΦ andΩ, whilen00 is the number of sample pairs with different labels in Φ andΩ. Cn2 is all possible combinations of sample pairs.NMI is one of the common external evaluation metrics that can measure the degree of agreement between two data distributions [39]. The value of NMI lies between 0 and 1, and the higher value indicates a better classification performance. The calculation equation is as follows:(24)NMI=∑φ∈Φ∑ω∈Ωpφ,ωlog⁡pφ,ω/pφpω∑φ∈Φpφlog⁡pφ∑ω∈Ωpωlog⁡pωwherep(φ,ω) is the joint probability function ofΦ andΩ, andp(φ) andp(ω) are the probability distribution functions ofΦ andΩ, respectively.F and ACC are derivations from the confusion matrix whose structure is reflected in Table2 [40, 41]. F comprehensively considers the proportion of true positive result in positive actual label and in positive classified result. ACC represents the proportion of all correctly classified results in total results. The scope of these two metrics is both [0,1], and the higher value indicates a better classification performance. F and ACC can be defined as (25) and (26), respectively.(25)F=2TP/TP+FPTP/TP+FNTP/TP+FP+TP/TP+FN(26)ACC=TP+TNTP+FN+FP+FNTable 2
The confusion matrix.
Actual label Positive Negative Classification result Positive True positive (TP) False positive (FP) Negative False negative (FN) True negative (TN)


## 4.3. Application to Fault Diagnosis of Rolling Bearing

To demonstrate the superiority of the proposed VMD-MAR-MHGWOSCA-SVM method, the experiment was carried out in a comparative form at each stage, as follows: EMD was applied for comparison in the signal processing stage; fuzzy entropy (FE) [42] and permutation entropy (PE) [5] were used in the feature extraction stage for comparison; GWO, SCA, and HGWOSCA were employed for comparison in the parameter optimization stage. The parameters of the contrastive experiments were set in the same way as the proposed method.In the proposed method, the first step is to decompose signal of each fault type into a set of IMFs. The mode numberK of VMD needs to be determined in advance with central frequency observation method. IfK is too large, mode mixing would occur because central frequencies of adjacent IMFs are too close. In contrast, it is difficult to weaken the nonstationarity effectively ifK is too small. Central frequencies under differentK value are listed in Table 3, whereK was ascertained by the sample of inner race fault with diameter of 0.007 inches (L1). As illustrated in Table 3, whenKwas set to 5, the first two central frequencies were relatively close, which meant excessive decomposition had occurred; hence parameterKwas set to 4. The original signal waveforms of different fault types are shown in Figure 3 and decomposition results of 0.007 inches faults (L1, L4, L7) are illustrated in Figure 4. It can be observed that waveforms from 9 types of faults are quite different and IMFs of one fault type are also discriminative to a certain extent.Table 3
Normalized central frequencies with differentK values.
Number of modes Normalized central frequencies 2 0.2256 0.0907 3 0.2978 0.2275 0.0999 4 0.2471 0.1738 0.0860 0.0127 5 0.3066 0.2834 0.2271 0.1148 0.0525 6 0.3072 0.2854 0.2364 0.2142 0.1139 0.0519 7 0.3170 0.2981 0.2757 0.2270 0.1242 0.0982 0.0493 8 0.3178 0.2990 0.2777 0.2363 0.2147 0.1231 0.0975 0.0492 9 0.4022 0.3133 0.2970 0.2768 0.2363 0.2146 0.1230 0.0975 0.0492Figure 3
Time domain waveforms of signals with different faults: (a) inner race 0.007 inches, (b) inner race 0.014 inches, (c) inner race 0.021 inches, (d) ball 0.007 inches, (e) ball 0.014 inches, (f) ball 0.021 inches, (g) outer race 0.007 inches, (h) outer race 0.014 inches, (i) outer race 0.021 inches.
(a) (b) (c) (d) (e) (f) (g) (h) (i)Figure 4
The VMD decomposition results of 0.007 inches fault signals: (a) inner race, (b) ball, (c) outer race.
(a) (b) (c)When decomposed components of a sample were obtained, a multivariate autoregressive (MAR) model was established for all IMFs. The values of SBC valueS(p) are shown in Figure 5, where the order of MAR model was determined according to the sample of fault type L1. The figure shows that as the model orderp increases,S(p) decreases gradually, but the value ofS(p) is no longer significantly decreased whenp increases to 10. Hence, the MAR model orderp is selected as 10 in the experiment. After establishing the MAR model, all parameters were identified based on the QR decomposition and thus the fault feature vectors were formed. Further, principal features were extracted from all features by using principal component analysis (PCA), thus achieving a balance between the computational cost and the integrity of feature information. The principal component (PC) distribution of the feature vectors obtained by blind parameter identification of MAR model is illustrated in Figure 6 and Table 4, while the 3D projection of feature vectors for different faults is shown in Figure 7.Table 4
Principal components of feature vectors.
Fault label Principal components PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13 L1 -0.465 0.498 -0.055 -0.516 0.095 0.008 -0.101 0.334 0.151 0.100 -0.040 0.078 -0.097 -0.459 0.510 -0.076 -0.515 0.116 0.003 -0.082 0.297 0.137 0.083 -0.032 0.076 -0.095  L2 -0.362 -0.252 0.257 -0.034 0.268 -0.861 -0.699 -0.181 -0.131 -0.702 0.053 0.071 -0.185 0.396 -0.547 -0.194 -0.347 0.161 -0.432 -0.415 0.160 0.059 -0.438 0.510 -0.070 -0.393  L3 0.883 0.374 1.127 0.456 0.168 0.145 -0.120 0.229 -0.036 0.061 0.054 -0.031 0.020 0.903 0.393 1.113 0.462 0.140 0.162 -0.158 0.221 -0.054 0.065 0.081 -0.054 0.030  L4 -0.724 -1.063 0.770 0.037 -0.593 0.167 0.025 -0.271 -0.099 0.045 0.053 0.124 -0.244 0.277 -1.649 0.483 -0.168 0.614 1.125 0.160 -0.182 0.285 -0.320 0.017 0.495 0.519  L5 -0.387 0.075 -0.232 -0.161 0.232 0.566 0.057 -0.033 -0.077 -0.069 0.124 0.013 0.072 -0.334 -0.002 -0.292 -0.107 0.242 0.600 0.082 0.003 -0.120 -0.102 0.127 0.039 0.082  L6 0.350 -0.174 0.128 -0.874 -0.099 -0.441 0.413 0.533 -0.149 0.306 -0.065 0.204 0.002 -1.335 0.256 -0.407 0.647 -0.262 0.114 0.339 0.128 -0.068 0.374 -0.047 -0.265 -0.284  L7 1.847 0.174 0.321 0.410 0.250 -0.255 -0.026 -0.182 -0.008 0.198 -0.283 0.083 -0.011 1.995 0.217 0.045 0.276 0.170 -0.364 -0.070 -0.297 -0.048 0.130 -0.086 0.061 0.005  L8 -1.110 0.912 0.024 1.937 -1.912 -0.797 0.348 0.720 -0.036 -0.714 0.657 0.399 0.280 -1.337 0.116 -0.700 0.553 0.894 -0.020 0.247 -0.459 -0.077 0.289 -0.036 0.308 0.335  L9 -0.845 0.547 -0.570 0.382 -0.046 0.345 0.295 0.237 0.015 -0.127 -0.100 -0.218 0.005 0.753 0.129 -0.638 0.636 0.267 0.380 0.208 -0.340 -0.274 -0.473 -0.328 -0.078 -0.162Figure 5
Selection of MAR model orderp.Figure 6
Principal component distribution of feature vectors.Figure 7
3D projection of feature vectors for different faults.During optimization stage, the feature vectors of different fault types were divided randomly into two parts with 30 and 29 vectors, respectively, and the 30 ones were used for training while the 29 ones for testing. The proposed MHGWOSCA approach was applied to optimize the penalty factorC and the kernel parameter g, where individual number and iteration time were set to 30 and 100, respectively, while searching ranges ofC and g were both [2-10,210]. Considering the maximum number of iterations and the number of optimization parameters, mutation amplitude was set to 0.5 and mutation period was set to 5. During optimization process, fitness values were calculated by fivefold cross validation for the training samples. Hence, the optimization effect of the givenC and g can be measured by average accuracy of cross validation, and then the optimalCand g were obtained. The convergence procedure of MHGWOSCA is depicted in Figure 8, from which it can be observed that, with the progress of iteration, the average fitness value of individuals rises rapidly at first and then tends to be stable; that is, the individuals are close to the global optimal solution. In addition, due to the effect of the mutation operators, the best three individuals α, β, andδchange greatly and periodically, trying to jump out of the local optimum without affecting the overall convergence trend. The shaded part shows the distribution of the convergence curve in 10 runs. The convergence curves of different optimization methods are compared in Figure 9, while the average fitness values for different optimization methods in the last 20 iterations are listed in Table 5. The convergence data are the average of 10 runs. The figure shows that the convergence rate of SCA is significantly lower than other methods; GWO and HGWOSCA both have good convergence trend; MHGWOSCA first converges to a high level in the initial stage and shows the best convergence effect at later stage of iterations among all methods.Table 5
The average fitness value of individuals during iterations.
Iteration number Optimization methods MHGWOSCA HGWOSCA GWO SCA 81 88.12 97.25 97.15 95.95 82 97.32 97.25 97.17 96.06 83 97.32 97.25 97.17 95.87 84 97.32 97.25 97.17 96.03 85 97.33 97.26 97.17 96.13 86 89.70 97.25 97.17 96.16 87 97.32 97.25 97.17 96.04 88 97.33 97.25 97.18 96.02 89 97.32 97.25 97.18 96.08 90 97.32 97.25 97.17 96.09 91 88.87 97.26 97.19 96.25 92 97.32 97.25 97.18 96.33 93 97.32 97.25 97.18 96.34 94 97.32 97.25 97.19 96.29 95 97.32 97.26 97.19 96.43 96 88.99 97.26 97.19 96.37 97 97.32 97.26 97.19 96.40 98 97.32 97.26 97.18 96.43 99 97.32 97.26 97.19 96.46 100 97.32 97.26 97.19 96.46Figure 8
Convergence procedure of proposed optimization approach: (a) overall trend, (b) partial enlarged detail.
(a) (b)Figure 9
Comparison of different optimization methods: (a) overall trend, (b) partial enlarged detail.
(a) (b)With the optimal parametersCand g, the SVM model optimized by MHGWOSCA was trained and applied to recognize testing samples. For an in-depth verification of the proposed method, diagnosis results were averaged over 10 runs and training samples were chosen at random in every run. Furthermore, deviation of each result was calculated as a reference for result evaluation. The optimalC and g are presented corresponding to the best accuracy among 10 runs.Fault diagnosis results for different methods are shown in Table6 and Figure 10. It can be seen from Table 6 that the evaluation values of the proposed VMD-MAR-MHGWOSCA-SVM method are the highest among all four metrics, i.e., 09573, 0.9554, 0.9808, and 0.9808, and the deviations of the evaluation values are also at a small level. Specifically, the comparison between EMD-MAR-MHGWOSCA-SVM and the proposed method indicates the effectiveness of VMD; the comparison of VMD-FE-HGWOSCA-SVM, VMD-PE-HGWOSCA-SVM, and the proposed method shows MAR model possess a better feature extraction ability; similarly, it can be concluded from the contrastive experiment among VMD-MAR-GWO-SVM, VMD-MAR-HGWOSCA-SVM, VMD-MAR-SCA-SVM, and the proposed method that MHGWOSCA improves the diagnosis performance of SVM. Figure 11 is the boxplots of evaluation values, illustrating the performance of different diagnosis methods. As shown in the figure, the proposed method achieves the highest evaluation values and possesses outstanding stability among all methods.Table 6
Fault diagnosis results with different methods.
Processing methods C g Result evaluation Adjusted Rand index (ARI) Normalized mutual information (NMI) F-measure (F) Accuracy (ACC) VMD-FE-HGWOSCA-SVM 184.0764 2.2262 0.9001 0.9046 0.9519 0.9525 [-0.034, 0.034] [-0.032, 0.034] [-0.019, 0.017] [-0.018, 0.017] VMD-PE-HGWOSCA-SVM 2.5155 414.0703 0.9216 0.9221 0.9635 0.9636 [-0.037, 0.053] [-0.043, 0.048] [-0.017, 0.025] [-0.017, 0.025] EMD-MAR-GWO-SVM 172.6222 0.0010 0.9284 0.9244 0.9662 0.9667 [-0.031, 0.038] [-0.033, 0.039] [-0.017, 0.018] [-0.016, 0.018] EMD-MAR-MHGWOSCA-SVM 254.3134 0.0339 0.9296 0.9236 0.9661 0.9667 [-0.072, 0.053] [-0.074, 0.055] [-0.041, 0.026] [-0.039, 0.026] VMD-MAR-GWO-SVM 515.8658 0.1028 0.9325 0.9311 0.9691 0.9693 [-0.048, 0.033] [-0.054, 0.039] [-0.024, 0.016] [-0.023, 0.015] VMD-MAR-HGWOSCA-SVM 324.7595 0.0925 0.9380 0.9366 0.9719 0.9720 [-0.039, 0.037] [-0.053, 0.037] [-0.018, 0.017] [-0.018, 0.016] VMD-MAR-SCA-SVM 212.6435 0.3115 0.9431 0.9401 0.9741 0.9743 [-0.028, 0.032] [-0.041, 0.032] [-0.012, 0.014] [-0.013, 0.014] VMD-MAR-MHGWOSCA-SVM 261.4197 0.1909 0.9573 0.9554 0.9808 0.9808 [-0.034, 0.025] [-0.037, 0.023] [-0.015, 0.012] [-0.015, 0.011]Figure 10
Comparison of evaluation values for different methods.Figure 11
Boxplots of evaluation values for different methods; the x-axis tick labels correspond to 1: VMD-FE-HGWOSCA-SVM, 2: VMD-PE-HGWOSCA-SVM, 3: EMD-MAR-GWO-SVM, 4: EMD-MAR-MHGWOSCA-SVM, 5: VMD-MAR-GWO-SVM, 6 VMD-MAR-HGWOSCA-SVM, 7: VMD-MAR-SCA-SVM, and 8: VMD-MAR-MHGWOSCA-SVM.


## 5. Conclusions

For rolling bearing fault diagnosis, a new method based on blind parameter identification of MAR model and mutation hybrid GWO-SCA optimized SVM is proposed in this study. Due to the nonstationarity of original signal, signals collected from different types of faults were firstly split into a set of IMFs by VMD, the decomposing mode numberKof which was determined with central frequency observation method. Then all the IMFs were employed to establish a MAR model, whose order was determined by SBC. Later, all parameters of the MAR model were identified blindly through QR decomposition and key features of parameter vectors were extracted via PCA, thus forming the fault feature vectors. Finally, a hybrid optimization method combining mutation operator, GWO, and SCA theory, termed mutation hybrid GWO-SCA (MHGWOSCA), was proposed to optimize SVM model which was subsequently applied to realize the pattern recognition for different fault samples. The contrastive experiments in Section 4 were implemented among the proposed VMD-MAR-MHGWOSCA-SVM method and other seven relevant methods, where nine types of faults with different locations and sizes were successfully diagnosed. The analysis of results indicates that the proposed method exhibits the best performance in all classification metrics, showing its outstanding precision and stability among all methods.



---

*Source: 3264969-2019-04-17.xml*