Misalignment Fault Diagnosis for Wind Turbines Based on Information Fusion

2.1. D–S Evidence Theory

The D–S evidence theory is a method of uncertainty reasoning, proposed by Dempster in 1967 and later improved and developed by Shafer [28]. The D–S evidence method can produce a probability interval to an uncertain event by fusing multiple evidences with known probability distribution. As an indeterminate reasoning method, D–S evidence theory uses weaker conditions than Bayesian, and has the ability to quantify unknown and uncertainty [29]. The evidence theory contains three important functions: basic probability assignment function, belief function, and plausibility function. The basic probability assignment function is the probability distribution of all possible faults in each state, the belief function is the lower bound of the probability of the fault event, and the plausibility function is the upper bound of the probability of the fault event. The belief function and the plausibility function can be obtained by calculating the sum of the basic probability assignment function, and the final decision is made after combining multiple evidences from different sources.

The D–S evidence theory consists of the following parts [30].

• Frame of discernment:

A variety of possible mutually exclusive hypothesis $X_i\left(i=1,2,\cdots ,s\right)$ of a question constitute a finite and non-empty set, which is called the frame of discernment, denoted as $\Omega =\left\{X_1,X_2,\cdots ,X_s\right\}$.

• Basic probability assignment (BPA) function:

BPA function is also known as the mass function. Suppose $H$ is a subset of $\Omega$, if function $m\left(H\right)$ satisfying

then, function $m\left(H\right)$ is called the basic probability assignment of $H$ on $\Omega$.

• Belief function:

In the frame of discernment, the belief function represents the sum of the basic probability assignment functions of all subsets of $H$. The expression of the belief function is as follows:

• Plausibility function:

In the frame of discernment, the plausibility function represents the degrees of belief for not denying $H$, which is the sum of the basic probability assignments of all the subsets intersecting $H$. The expression of the plausibility function is as follows:

• Dempster’s rule of combination:

Dempster’s rule of combination is used to combine the BPA functions of multiple evidences. Although this rule is controversial at present, the authors of [31] have showed that it behaves perfectly when evidences do not conflict reciprocally. Only if we integrate conflicting evidences do we need to improve it. In this paper, there is no serious and complete conflict among the outputs from vibration signal, temperature signal, and stator current signal as evidences in this study. Therefore, Dempster’s rule is still used here.

Suppose there are $n$ independent evidences (sensors or expert opinions), $H_1,H_2,\cdots ,H_n$ (are subsets of $\Omega$), the BPA of them are $m_1,m_2,\cdots ,m_n$. Then, Dempster’s rule for the BPA functions on $\Omega$ is as follows:

Specifically, it can be expressed as follows:

where the expression of $K$ is as follows:

where K is the degree of conflict between evidences. When K = 1, the evidences are completely conflicted and cannot be synthesized by this formula; when K tends to 1, the evidences are highly conflicted, and synthesizing by this formula may lead to results contrary to fact [32].

• Decision rules:

The decision rule is to draw a diagnosis based on the uncertain interval $\left[bel\left(H\right),pl\left(H\right)\right]$ of the evidence. In the interval of $\left[0,1\right]$, the uncertainty of a proposition is shown in Figure 1.

In Figure 1, $\left[0,bel\left(H\right)\right]$ belongs to the support interval, $\left[0,pl\left(H\right)\right]$ is the accept interval, $\left[pl\left(H\right),1\right]$ is the rejection interval, and $\left[bel\left(H\right),pl\left(H\right)\right]$ is the uncertain interval.

When making a decision, choose a value in the uncertain interval as the final trustworthiness of the proposition. If this value has the highest trustworthiness among the possible hypothesis, this assumption is the final decision result.

2.2. Posterior Probability Least Squares Support Vector Machine

In the study, the fault samples collected are limited, while support vector machine (SVM) and LSSVM can obtain high diagnosis accuracy based on small sample data. Moreover, the speed of LSSVM is faster than that of SVM, so LSSVM is selected to be the initial classifier to judge the state. As the input parameters of the D–S evidence fusion are basic probability assignments in all classification spaces, the hard output (whether or not) of the traditional classifier has to be converted to a soft one (probability) [33]; that is, the output of the classifier must be changed to the posterior probability output. For the two-class problem, the posterior probability can be calculated using the sigmoid function to map the output $f\left(x\right) \left(+1,-1\right)$ of the LSSVM to the $\left[0,1\right]$ interval. Assuming that the probability is consistent with the sigmoid distribution, the posterior probability can be calculated [34]:

where $f$ is the classification result of the standard LSSVM, $p\left(y=1/f\right)$ is the probability when the classification is correct under the condition that the output value is $f, p\left(y=-1/f\right)$ is the probability when the classification is wrong under the condition that the output value is $f$, and $A$ and $B$ are parameters. So, the key to calculating the posterior probability is to obtain parameters $A$ and $B$. The posterior probability least squares support vector machine model is usually established by first establishing a standard LSSVM model, and then obtaining $A$ and $B$ on the training set $\left(f_i,t_i\right)$, where $t_i$ is the target probability output of the standard LSSVM:

where $N_{+}$ is the number of positive samples; $N_{-}$ is the number of negative samples; and the problem of obtaining parameters $A$ and $B$ is to solve the minimum likelihood optimization problem of the following, i.e.,

where

The Hessian matrix for solving (9) is as follows:

In order to get the minimum value of (9), the Hessian matrix must be positively determined. So, $A$ and $B$ are finally obtained by solving all the eigenvalues of the matrix that are greater than zero. The posterior probability can be obtained.

It is proved that the posterior probability least squares support vector machine with sigmoid function works well in practical applications [35], but this method can only be used for the two-class problem. The main methods for extending LSSVM from two-class to multi-class are the “one-versus-one” and “one-versus-all” methods. The Platt algorithm calculates the probability formula for each classifier as follows, where $p_m$ is the probability that sample $x$ belongs to the i-th class [35]:

2.3. The Improved Artificial Bee Colony

There are three kinds of kernel functions commonly used in LSSVM: linear kernel function, polynomial kernel function, and radial basis function (RBF) ($K\left(x_i,x_i\right)=exp\left(-\parallel x_i-x_i{\parallel }^2/{\sigma }^2\right)$, where $\sigma$ is the kernel width). Many studies and experiments [36] show that, compared with other kernel functions, RBF can map the original space into an infinite dimensional space and find the hyperplane better. It is a better choice as the kernel function. Therefore, it is necessary to select the regularization parameter $\mathsf{\gamma }$ (necessary for LSSVM, determining the trade-off between the training error minimization and smoothness) and the kernel squared bandwidth ${\sigma }^2$.

Choosing a better parameter value can greatly improve the performance of the LSSVM classifier and the accuracy of diagnosis. At present, the commonly used methods include trial and error, cross validation, grid search, and intelligent optimization algorithm [37]. Among them, the trial and error method not only consumes time and energy, but also the choice of parameters is greatly affected by subjective factors; the cross validation method divides the data set into training, validation, and testing, and different proportions will lead to different optimal models and optimal parameters; and the grid search method optimizes the model according to the set step size in the upper and lower limits of parameters, and then determines the optimal parameters, so the search speed is too slow and the precision is not high. Therefore, the advantages of the intelligent optimization algorithm are highlighted. It realizes the optimal distribution of food by simulating the behavior of animals in the population (interact information and cooperation among individuals). A swarm intelligence optimization algorithm is easy to implement and has high efficiency, so it is applied to the parameter optimization process of LSSVM.

Swarm intelligence optimization algorithms include genetic algorithm, particle swarm optimization, artificial fish swarm algorithm, artificial bee colony algorithm, and so on. Among them, artificial bee colony algorithm (ABC) is an optimization algorithm proposed in recent years, which not only has good optimization ability, but also controls less parameters in the process. Furthermore, it is simple, flexible, and easier to implement. The research [38] shows that the optimization performance of ABC is better than that of genetic algorithm and particle swarm algorithm, and the classification diagnosis accuracy of LSSVM optimized by ABC is higher than that of LSSVM optimized by genetic algorithm and particle swarm algorithm.

However, ABC has some shortcomings, such as slow convergence speed in the later stage of operation and the fact that it is easy to fall into local optimum. Therefore, in this paper, on the one hand, chaotic initialization is introduced in the artificial bee colony algorithm, which is used to initialize the population position to improve the diversity of the population and the ergodicity of the population search process. On the other hand, in the collecting bees stage of the artificial bee colony algorithm, the bees are divided into two parts: one part collects the optimal information of the region according to the original algorithm, and the other does Lévy flight around the global optimal solution to improve their global search capabilities. At the same time, in the observing bees stage, a search strategy based on the current local optimal solution (called pbest) is adopted to improve the local search ability of the algorithm.

(1) The logistic chaotic map is proposed to initialize the population. The equation for the logistic chaotic map is as follows:

In the formula,$y_t\in \left(0,1\right)$, $t$ is the number of iterations of the chaotic sequence, $\mu$ is the control parameter of the chaotic sequence, and the value range is $\left[3.75,4\right]$ [39].

(2) Lévy flight was introduced in the evolution strategy to improve the performance of the algorithm and achieve good results [39]. The calculation method is based on

where $\alpha$ is the characteristic index, which usually satisfies $0<\alpha <2$. $\Gamma (·)$ is the Gamma function defined as

Its update equation is as follows:

where $\alpha$ is the step length, which usually meets the standard normal distribution, and $L(·)$ is the random search path for Lévy flight.

(3) In the observing bees stage, for any current solution in each generation, the top $p\%$ solutions are randomly selected among all current solutions, and the best one (called pbest) can be used to balance global search capabilities and local development capabilities. The neighborhood search formula is as follows:

where $k\in \left\{1,2,\cdots ,S_N\right\}$, $S_N$ is the number of solutions for the bee colony, $j\in \left\{1,2,\cdots D\right\}$, $D$ is the dimension of the optimization problem, $k\ne i$, ${\phi }_{ij}\in \left[-1,1\right]$, and ${\varnothing }_{ij}\in \left[0,1.5\right]$.

4.1. Data Processing

After the vibration signal, temperature signal, and stator current signal under four working conditions are collected, in order to make better use of them and get good diagnosis results, the feature indexes in the time, frequency, and time–frequency domain are extracted. Table 2 shows a 21-dimension mixed feature library of the vibration signal.

Suppose signal x ($x_0, x_1, x_2,\cdots , x_{N-1}$) is a discrete time series with a finite length, the calculation formulas of time domain characteristic indexes are shown in Table 3, where $\overline{x}$ is the mean value of the signal, $x^{\prime }$ is the average amplitude, and $x_p$ is the peak value of the signal.

In signal analysis, power spectrum analysis is usually used to extract the frequency domain index. Center of gravity frequency, mean square frequency, root mean square frequency, and frequency variance are commonly used. The sampling frequency is set as $f_s$, and the calculation formula of each index is shown in Table 4, where $S\left(\omega \right)$ is the power spectrum of discrete time series, $S\left(\omega \right)=X\left(\omega \right)·\overline{X\left(\omega \right)}$, $X\left(\omega \right)=\sum _{i=0}^{N-1}x\left(i\right)e^{-j\pi \omega }, \omega$ is the angular frequency.

Time–frequency analysis is a fault diagnosis method that combines the law and reason of frequency changing with time. In this paper, image extended empirical mode decomposition (IEMD) is used to process the vibration signal, and dual tree complex wavelet transform (DTCWT) is used to process the stator current signal (see the literature [43] for details).

The gearbox tooth temperature $T_1$ and the generator rotor shaft temperature $T_2$ are selected as the characteristic values of the temperature signal. Construct a two-dimensional vector of the temperature signal: $X=\left[T_1,T_2\right]$.

Table 5 is a mixed feature library with a total of 29 dimensions in the time domain, frequency domain, and time–frequency domain of the stator current signal (see the literature [41] for details).

In order to eliminate the influence of different input dataset dimensions and large numerical differences, the original dataset is normalized, i.e.,

where $x$ is the value to be normalized, $y_{min}$ is the lower bound of the normalized interval, and $y_{max}$ is the upper bound of the normalized interval. In this paper, $y_{min}=0$, $y_{max}=1$, and the vector is normalized by column.

Because of the high dimensionality of the constructed vectors of the vibration signal and the stator current signal, not only does the amount of calculation increase, but also some difficulties are brought to fault diagnosis [44]. In order to make better use of various information and obtain good diagnostic results, the feature vectors are subjected to dimensionality reduction using t-SNE.

t-SNE based on conditional probability retains the similarity between high-dimensional and low dimensional space data and adopts symmetric objective function, and t distribution in low-dimensional space replaces Gaussian distribution, which solves the problem of crowding and clear visualization in low-dimensional space [45]. Its implementation steps are as follows:

(1) Define a high-dimensional data set: $x=\left\{x_1,x_2,\cdots ,x_n\right\}$.

(2) Compute the complexity parameter of the value equation $c$:

where $p_i$ is the conditional probability of data points (other than $x_i$) with respect to $x_i$, $p_{j∕i}$ is the conditional probability of high-dimensional data, $p_{ij}$ is the joint probability density in the high-dimensional space, and $q_{ij}$ is the joint probability density in the low-dimensional mapping space.

(3) Define the optimization parameters: the number of iterations $\mathrm{T}$, the learning rate $\mathsf{\eta }$, and the momentum factor at the $\mathrm{t}$th ($\mathrm{t}\le \mathrm{T}$) iteration $\mathsf{\alpha }\left(\mathrm{t}\right)$ ($0<\mathsf{\alpha }\left(\mathrm{t}\right)<1$). The value equation $\mathrm{c}$ is learned by the gradient descent method, and the low-dimensional mapping of the high-dimensional data is finally obtained:

where $y_i$ and $y_j$ are the mapping of the high-dimensional data $x_i$ and $x_j$ in the low-dimensional space.

In order to speed up the optimization process and prevent trapping into local minima, a relatively large momentum condition is imposed on the descent process. The current gradient value is summed to the previous gradient value for each iteration and then decays exponentially to determine the coordinates of the low-dimensional data. The momentum formula is as follows:

where $y$ is the data in the low-dimensional space.

4.2. The Fault Diagnosis Results

In this paper, “one-versus-all” is used to extend LSSVM from two classifications to multiple classifications. That is, each time, one fault is selected as one type, and the rest of the states are selected as another type. In order to produce the posterior probabilities of the four classifications in the vibration feature space, four two-class LSSVM are constructed, and each LSSVM calculates a set of $A$ and $B$, and then the corresponding posterior probability is calculated according to (5) and (6). In the same way, the probability vectors of the temperature and stator current signal classifiers for the four states can be obtained as the BPA of D–S evidence fusion.

The five-dimensional feature vectors of the vibration signal after t-SNE dimensionality reduction are used as the inputs, and the four working conditions of the transmission system are used as outputs to train the LSSVM, which is optimized by the improved artificial bee colony algorithm. The parameters of the four two-classification LSSVM in the vibration feature space are shown in Table 6. Four samples are selected, such as samples 5, 44, 82, and 130, and the corresponding BPA1 calculated is shown in Table 7.

The two-dimensional feature vectors of the temperature signal are used as the inputs, and the four operating states of the transmission system are used as the outputs to train the optimized LSSVM. The parameters of the four binary LSSVM in the temperature feature space are shown in Table 8. The BPA2 calculated from the same four samples is shown in Table 9.

The four-dimensional vectors after the dimensionality reduction of the stator current signal are used as inputs, and the four operating states of the transmission system are as outputs to train the optimized LSSVM. The parameters of the four two-class LSSVM in the stator current feature space are shown in Table 10, and the BPA3 calculated by the same four samples is shown in Table 11.

Then, the probability assignments are calculated after the fusion of the three BPAs. The category with the highest degree of belief is selected as belonging to the class of the fusion model. Table 12 shows the basic and the fusion probability of the three LSSVM outputs for the selected test samples. Table 13 shows the fusion and classification results of the four test samples. Figure 3 shows the test samples’ diagnosis results, in which “0” indicates normal operation, “1” indicates parallel misalignment, “2” indicates angular misalignment, and “3” indicates integrated misalignment.

In order to better evaluate the performance of the fault diagnosis method, three indexes are adopted: the training set classification accuracy, the testing set classification accuracy, and the fault false alarm rate. The fault false alarm rate means that the fault does not actually occur, but the fault detection alarm is given by the detection system. The false alarm rate equals the number of false alarm samples divided by the total number of actual fault-free samples. Table 14 compares the results of the sample sets diagnosed by the indexes of a single signal (vibration, temperature, or current signal) with the D–S evidence fusion.

From Table 14, it can be seen that the accuracy of D–S fusion is higher than that of any single signal, and the failure false alarm rate is equal to zero, lower than others, which proves the advantage of information fusion in the diagnosis of wind turbine misalignment fault.