A Fault Diagnosis Method for CNC Machine Feed System Based on DoubleEnsemble-LightGBM Model

1. Introduction

According to statistics, among the failures of CNC, the mechanical body failure accounts for about 57%, and the electrical system failure accounts for about 37.5%.CNC system failure accounts for only 5.5%, and most of the current CNC machines have the self-diagnosis function of electrical and CNC systems [1]. Fault of mechanical body is the key and difficult point of current research. With the continuous development of data acquisition technology, information technology and artificial intelligence technology, Fault diagnosis methods have also experienced the development process from artificial experience diagnosis to intelligent diagnosis, from single sensor diagnosis to multi-sensor fusion diagnosis.

Grether [2] took Siemens CNC machines as the research object. According to the expert knowledge in the field of fault diagnosis, an ontology-based knowledge representation structure is proposed, and then the SimRank algorithm is used to calculate the similarity between the fault phenomenon and the fault cause in the case base to realize the fault diagnosis of the CNC machine. Wang [3] established the fault tree model of CNC machines, and on this basis, a deep neural network model is constructed to classify and identify the features. The average recognition rate of BP network after feature reduction is 86%. Kemal [4] used Morlet Wavelet analysis to extract the features of vibration signals of CNC machines, and then proposed a deep LSTM model for fault classification, which effectively improves the classification accuracy. Shan et al. [5] proposed to arrange multiple sensors at different positions of the ball screw. The fault location of the ball screw is realized by carrying out weight distribution on the fault sensitivity indexes of different sensors and combining with a convolutional neural network, Finally, the effectiveness of the method is verified by the ball screw test bench. Zhang et al. [6] applied a new unsupervised learning method, generalized normalized sparse filtering, to rolling bearing intelligence under complex working conditions. The experiment proves that the method can obtain higher diagnosis accuracy with less training samples. Chen et al. [7] proposed a multi-scale feature alignment CNN for bearing fault diagnosis under different working conditions. The displacement invariance of CNN is improved, and the effectiveness and advancement of the method are verified by rolling bearing experiments. Moslem [8] proposed a domain adaptive method based on deep learning for cross-domain ball screw fault diagnosis. A deep convolutional neural network is used for feature extraction, and the maximum average difference metric is proposed to measure and optimize the data distribution under different working conditions. The effectiveness of the proposed method is proved by the experiment with the monitoring data of the ball screw under real working conditions. Pandhare et al. [9] Collected the vibration acceleration signals at five different positions on the ball screw test bench. A data domain adaptive fault diagnosis method based on convolutional neural network is proposed. The method minimize that maximum average difference of the high-level representation between the source domain data and the target domain data. The average diagnostic accuracy of the model is 98.25%. Jin et al. [10] proposed an end-to-end adaptive anti-noise neural network framework without manual feature selection and denoising processing. The convolutional feature extraction part of the network takes the exponential linear unit as the activation function, the extracted features are learned and classified by a gated recurrent neural network improved by an attention mechanism, The accuracy of bearing fault diagnosis under the conditions of noise and variable load is effectively improved. Patel T H et al. [11] modeled the mixed fault, analyzed its vibration signal, and then recognized the mixed fault pattern. Abbasion S et al. [12] applied the combination of Wavelet packet decomposition and support vector machine to the mixed fault diagnosis of bearings. Lei Y G et al. [13] proposed a classification method based on adaptive fuzzy neural inference to diagnose the composite faults of electric locomotives. Delgado. M et al. [14] extracted fault features from the motor current signal and vibration signal, and used partial least squares to reduce the dimensionality of the extracted features and construct feature vectors. Finally, they used a support vector machine (SVM) model to achieve the diagnosis of motor interturn short-circuit fault; Yang et al. [15] first derived a prediction formula for motor faults, and then used a Kalman filter to track and estimate the change trend of motor speed to determine whether the motor has failed. Wang et al. [16] used a multi-task shared classifier based on incremental learning to achieve better fault diagnosis of support bearings under various working conditions. Li et al. [17] proposed a method based on attention mechanism to solve the problem of low accuracy and poor stability of the model caused by unbalanced data sets. Experiments show that the method has good diagnosis effect under unbalanced data conditions. Xu et al. [18] used an improved method of combining multi-scale convolutional neural network with feature attention mechanism to improve the generalization ability of the model. Zhang et al. [19] proposed an instance-based transfer learning method to solve the problem of insufficient labeled samples in the application of ball screw fault diagnosis. Azamfar et al. [20] proposed a deep adaptive fault diagnosis method for cross-domain fault diagnosis of ball screw. Wu et al. [21] adopted a fault diagnosis method combining domain antagonistic neural network and attention mechanism. Experimental results show that this method has great potential in cross-domain diagnosis of rolling bearings. Huang Min et al. [22] proposed a method to solve the problem of data distribution deviation in the fault diagnosis of support bearing migration. The experimental results show that the method can support bearing migration fault diagnosis suitable for different working conditions.

Based on the comprehensive analysis of the research status of fault diagnosis of feed system, this paper’s primary contributions can be summarized as follows:

• In order to solve the problem of fault diagnosis of CNC machine feed system under variable speed conditions, a fault diagnosis method based on multi monitoring signals, multi domain feature extraction and DoubleEnsemble-LightGBM integrated learning model is proposed. The experimental results show that this method achieves better diagnosis effect than Xgboost and other advanced integrated learning models;
• A variety of monitoring signals including vibration signal, noise signal and current signal are collected. The monitoring signals are preprocessed by singularity elimination, trend item elimination, Wavelet threshold denoising, and then the time-domain, frequency-domain feature indexes and IMF information entropy of the monitroring signals are extracted. Finally, the multi-dimensional mixed domain feature set is constructed;
• Based on the LightGBM model, the DoubleEnsemble-LightGBM fault diagnosis model is constructed by introducing the sample reweighting mechanism based on learning trajectory and the feature selection mechanism based on shuffling technology, which realizes the intelligent fault diagnosis of CNC machine feed system.

2. Relevant Theories

2.1. CEEMDAN Decomposition and IMF Information Entropy

2.1.1. CEEMDAN Decomposition

CEEMDAN（Complete Ensemble Empirical Mode Decomposition with Adaptive Noise) overcomes the mode mixing problem of EMD by adding adaptive white noise. And can effectively reduce the residual white noise in the IMF components obtained after decomposition [23].

The specific process of CEEMDAN decomposition is as follows:

1. adding k times of random Gaussian white noise with mean value is 0 into the signal $x(t)$ to be decomposed, construct the sequence $x_i(t)$ of the k times experiment according to formula (1), where $i=1,2,\cdot \cdot \cdot ,K$.

$$x_i(t)=x(t)+{\epsilon }_0{\delta }_i(t)$$

(1)

Where: ${\delta }_i(t)$ is the random Gaussian white noise added in the ith experiment; ${\epsilon }_0$ is the weight coefficient of the Gaussian white noise.

2. carry out EMD decomposition on that sequence $x_i(t)$, take the average value of the first IMF component obtained from the k times experiment as the first IMF component obtained from CEEMDAN decomposition, and see formula (2) for calculation. See formula (3) for the calculation of the residual signal after the first decomposition.

$${\overline{IMF}}_1(t)=\frac{1}{K}{\displaystyle \sum _{i=1}^{K}IM{F}_1^i(t)}$$

(2)

$$r_1(t)=x(t)-{\overline{IMF}}_1(t)$$

(3)

3. a new sequence $r_1(t)+{\epsilon }_1{E}_1\left({\delta }_i(t)\right)$ is obtained by adding k times specific noise to $r_1(t)$, and the EMD decomposition is carried out, calculate the second IMF component obtained by CEEMDAN decomposition according to formula (4).

$${\overline{IMF}}_2(t)=\frac{1}{K}{\displaystyle \sum _{i=1}^K{E}_1\left(r_1(t)+{\epsilon }_1{E}_1\left({\delta }_i(t)\right)\right)}$$

(4)

Where: $E_1(\cdot )$ is the first IMF component obtained after EMD decomposition; ${\epsilon }_1$ is the weight coefficient for adding noise to $r_1(t)$.

4. calculating a margin signal $r_m(t)$$m=2,\cdot \cdot \cdot ,M$ according to the formula (5), and obtaining the m+1th IMF component of the CEEMDAN in the same way as the step (3), See formula (6) for calculation.

$$r_m(t)=r_{m-1}(t)-{\overline{IMF}}_m(t)$$

(5)

$${\overline{IMF}}_{m+1}(t)=\frac{1}{K}{\displaystyle \sum _{i=1}^K{E}_1\left(r_m(t)+{\epsilon }_m{E}_m\left({\delta }_i(t)\right)\right)}$$

(6)

In the formula, $E_m(\cdot )$ represents the mth IMF component obtained after EMD decomposition of a certain sequence; ${\epsilon }_m$is the weight coefficient for adding noise to $r_m(t)$.

5. repeating the step (4) to calculate other IMF components of the CEEMDAN until the number of extreme points of $r_m(t)$is less than 2. Eventually, the signal $x(t)$ is decomposed into m IMF components and a residual component $R(t)$.

$$R(t)=x(t)-{\displaystyle \sum _{m=1}^M{\overline{IMF}}_m(t)}$$

(7)

2.1.2. False Modal Component Rejection

The IMF components obtained by CEEMDAN decomposition may contain false modal components, and the spurious modal components need to be rejected. The correlation coefficient can be used to describe the degree of correlation between the IMF component and the original signal. The closer the correlation coefficient is to 1, the more useful information the component contains, the stronger the correlation with the original signal is. Therefore, the false modal components after CEEMDAN decomposition can be adaptively eliminated through the correlation coefficient.

The correlation coefficient $C_m$ between the mth IMF component and the original signal is calculated as follows:

$$C_m=\frac{{\displaystyle \sum _{i=1}^N(x_i-\overline{x})(y_i-\overline{y})}}{\sqrt{{\displaystyle \sum _{i=1}^N{(x_i-\overline{x})}^2}}\sqrt{{\displaystyle \sum _{i=1}^N{(y_i-\overline{y})}^2}}}$$

(8)

Where: $x_i$is the ith element value in the original signal sequence. $\overline{x}$ is the average value of the original signal sequence. $y_i$ is the value of the ith element in the mth IMF component. $\overline{y}$is the average value of the mth IMF component. $N$ is the signal sequence length.

Albert et al. [24] gave a formula for calculating the adaptive threshold of the correlation coefficient, see Equation (9). If $C_m<\mu$, then the mth IMF component will be rejected.

$$\mu =\frac{\max (C_m)}{10\times \max (C_m)-3},m=1,2,\cdot \cdot \cdot ,M$$

(9)

In the formula, M is the number of IMF components decomposed from the original signal and is $\max (C_m)$ the maximum correlation coefficient value.

2.1.3. Calculation of IMF Information Entropy

In the field of fault diagnosis, entropy can effectively reflect the complexity of the signal and describe its nonlinear characteristics. A single entropy value is often difficult to fully describe the signal characteristics, so extract multiple information entropy eigenvalues are extracted simultaneously. Assume that K effective IMF components are obtained after the signal $x(t)$ is decomposed by CEEMDAN, denoted as $u_i(t),i=1,2,\cdot \cdot \cdot ,k$.

1. Energy entropy of IMF

Energy entropy is an index that can characterize the energy complexity of a signal. The IMF energy entropy is calculated as follows:

First, the energy value of each effective IMF component is calculated by Equation (10):

$$E_i={\displaystyle {\int }_{-\infty }^{+\infty }{\left|u_i(t)\right|}^2}dt,i=1,2,\cdot \cdot \cdot ,k$$

(10)

Then, the total energy value is calculated by Equation (11):

$$E={\displaystyle \sum _{i=1}^k{E}_i}$$

(11)

Finally, the IMF energy entropy is calculated by Equation (12):

$$H_{\mathrm{E}}=-{\displaystyle \sum _{i=1}^k{h}_i\lg \begin{array}{cc}{h}_i& \end{array}}$$

(12)

Where: $h_i=E_i/E$represents the proportion of the energy value of the ith IMF component to the total energy value.

2. Power spectrum entropy of IMF

Power spectrum entropy can reflect the change of signal energy in frequency domain. The IMF power spectrum entropy is calculated as follows: First, each effective IMF component $u_i(t)$is Fourier transformed to obtain $u_i(\omega ),i=1,2,\cdot \cdot \cdot ,k$. Then, the power spectrum of each effective IMF component is calculated by Equation (13):

$$S_i=\frac{1}{2\pi N}{\left|u_i(\omega )\right|}^2$$

(13)

Finally, the IMF power spectrum entropy is calculated by Equation (14):

$$H_{\mathrm{F}}=-{\displaystyle \sum _{i=1}^k{p}_i\lg \begin{array}{cc}{p}_i& \end{array}}$$

(14)

Where: $p_i=S_i/{\displaystyle \sum _{i=1}^k{S}_i}$ represents the proportion of the power spectrum of the ith IMF component to the total power spectrum.

3. The singular spectral entropy of IMF

Singular spectral entropy can quantitatively describe the complex state characteristics of time series. The calculation of the IMF singular spectral entropy is as follows:

First, each IMF component is formed into a characteristic matrix A:

$$\mathit{A}={\left[u_1(t),\cdot \cdot \cdot ,u_k(t)\right]}^{\mathrm{T}}$$

(15)

Then, the singular values ${\lambda }_i,i=1,2,\cdot \cdot \cdot ,k$ of the characteristic matrix A are computed.

Finally, the IMF singular spectral entropy is calculated by Equation (16):

$$H_{\mathrm{S}}=-{\displaystyle \sum _{i=1}^k{q}_i\lg \begin{array}{cc}{q}_i& \end{array}}$$

(16)

Where: $q_i={\lambda }_i/{\displaystyle \sum _{i=1}^k{\lambda }_i}$ represents the proportion of the ith singular value to the sum of all singular values.

2.2. LightGBM Algorithm

Light GBM [25] (Light Gradient Boosting Machine) is a lightweight gradient lifting model. It is an optimized implementation based on the classical ensemble learning model GBDT [26]. The principle of GBDT is shown in Figure 1.

The basic idea is to use the decision tree as a weak classifier, a plurality of weak classifiers are iteratively trained through a gradient lifting strategy, and all the weak classifiers are combined in a linear addition mode to form a strong classifier with better classification effect.

On the basis of GBDT, LightGBM is optimized as follows:

(1) Gradient-based One-Side Sampling (GOSS) algorithm is used to compress the training data samples without loss of accuracy, and its basic idea is to discard some samples that are not helpful to the calculation of information gain. Then the data calculation amount can be reduced, and the operation cost is greatly reduced.

(2) The Exclusive Feature Bundling (EFB) algorithm is used to merge the mutually exclusive features in high-dimensional data into one feature, which can effectively reduce the feature dimension and reduce the computational load.

(3) Histogram algorithm is used to improve the node segmentation strategy of decision tree. The basic idea is to discretize the continuous floating-point eigenvalues into K integers and construct a histogram with width K. This can greatly reduce the computational time and memory consumption, and has little impact on the overall classification accuracy of the model under the framework of gradient boosting. At the same time, it has the effect of regularization, which can prevent the model from overfitting and enhance the stability and robustness of the model.

(4) LightGBM uses a leaf-by-leaf growth strategy, which splits each time by finding the one with the largest split gain from all current leaf nodes, and so on. In the case of the same number of splits, the strategy of growing by leaves can reduce more errors and get better accuracy. However, in this way, a deeper decision tree may be grown, resulting in overfitting of the model, so the maximum depth limit is added to the growth strategy by leaf.

To sum up, LightGBM not only inherits the advantages of GBDT, but also greatly improves the training efficiency and memory consumption. Compared with other integrated learning models, this model is easier to deal with large-scale data, and requires low computing power of the equipment. Therefore, LightGBM is used as the basic model for mechanical fault diagnosis of feed system.

2.3. DoubleEnsemble Algorithm

Double Ensemble is a new ensemble algorithm framework that can be used with a variety of machine learning models. It includes two key technologies, one of which is the sample reweighting technology based on learning trajectory, which can give different weights to different samples in the process of model training, thus reducing the interference of simple samples and noise samples and enhancing the training of key samples. The feature selection technology based on shuffling mechanism can help the model to automatically screen sensitive features in the training process, thus effectively improving the accuracy of the model and reducing the risk of overfitting.

The algorithm flow (pseudocode) of DoubleEnsemble is shown in Table 1.

The algorithm sequentially trains K machine learning sub-models, denoted as${\mathit{M}}^1$、···、${\mathit{M}}^K$, all sub-models are weighted and integrated according to formula (17), and the integrated model ${\overline{\mathit{M}}}^K（\cdot ）$ is taken as the final output of the algorithm.

$${\overline{\mathit{M}}}^K（\cdot ）=\frac{1}{K}{\displaystyle \sum _{i=1}^K{a}_i{\mathit{M}}^i（\cdot ）}$$

(17)

Where: $a_i$is the weight coefficient of the ith sub-model ${\mathit{M}}^i$.

The training data consists of a feature matrix X and a label vector y. $X={[x_1,\cdot \cdot \cdot ,x_N]}^{\mathrm{T}}\in {\mathrm{R}}^{N\times F}$,${\mathit{x}}_i$ represents the feature set of the ith sample, N is the total number of training samples, and F is the dimension of the feature set. $\mathit{y}=(y_1,\cdot \cdot \cdot ,y_N)$, $y_i$ represents the fault label of the ith sample. For the first sub-model ${\mathit{M}}^1$, the algorithm will use all the feature indexes in the feature set of the training data for training, i.e., ${\mathit{f}}^1$= [F], the initial sample weights are set to${\mathit{w}}^1$ = (1,···,1). The subsequent sub-models are trained based on the newly selected feature set ${\mathit{f}}^k\subseteq [F]$and the updated sample weights${\mathit{w}}^k=(w_1^k,\cdot \cdot \cdot ,w_N^k)$, Where ${\mathit{w}}^k$ and ${\mathit{f}}^k$ are obtained through sample re-weighting based on learning trajectory and feature selection based on shuffling mechanism algorithm [27], respectively.

3. Model: Multi-Domain Feature and DoubleEnsemble-LightGBM

The CNC machine feed system is a complex system with multi-mechanical components, and it is difficult to describe its fault state by the characteristics in a single domain. In order to reflect the health status of the feeding system more comprehensively, the time domain characteristic indexes, the frequency domain characteristic indexes and the time-frequency domain characteristic indexes of various monitoring signals including vibration signals, noise signals and current signals are firstly extracted, and a multi-dimensional mixed domain feature set as shown in Figure 2 is constructed.

In addition, considering that the total dimension of the multi-dimensional mixed domain feature set reaches hundreds of dimensions, it may contain invalid features, which will have a negative impact on model training.In addition, there may be simple samples and useless high noise samples in the collected training samples, which leads to poor training performance of the model and easy to overfit.Therefore, the fault diagnosis model is further optimized, and multiple LightGBM classification sub-models are trained and integrated through the DoubleEnsemble algorithm, so as to construct the DoubleEnsemble-LightGBM model as shown in Figure 3 for intelligent identification of the fault mode of the feed system.

4. Experimental Results

4.1. Data Set Description

4.1.1. University of Ottawa Variable Speed Bearing Failure Open Data Set

The vibration data of ER16K deep groove ball bearing under different speed conditions were collected from the variable speed bearing fault data set of the University of Ottawa in Canada, and the sampling frequency was 200 kHz. The fault types of bearings include: health, inner ring fault, outer ring fault, rolling element fault, and compound fault of inner and outer rings and rolling elements. Speed changes include: speed up (from 846 R/min to 1428 r/min), speed down (from 1734 r/min to 822 r/min), speed up first and then speed down (from 882 r/minute to 1518 r/minute and then to 1260 r/minutes), First decrease and then increase (from 1452 r/min to 888 r/min and then to 1236 r/min).

Firstly, five kinds of original data collected from the data set under four speed conditions (speed up, speed down, speed up and then speed down, speed down and then speed up) are divided into samples, and each sample contains 2000 data points. Then, the obtained samples are divided into the training set and the test set in a ratio of 8:2. The sample distribution of the data set and the corresponding relationship of the fault labels are shown in Table 2.

4.1.2. Dataset of Feed System Test Bench

Based on the transmission principle and mechanical structure of the X-direction feed system of the vertical machining center, a feed system test bench as shown in Figure 4 is built.The test bench is made of heavy steel. The model and specification of the key parts used in the test are the same as those of the vertical machining center. The model of the ball screw pair is Taiwan Shangyin R4010FSI, the model of the rolling bearing is Japan NSK angular contact ball bearing 30TAC62B, the guide rail pair is a roller type rail with good rigidity, and the driving motor is a three-phase AC servo motor.

The model and parameters of the data acquisition equipment used in the experiment are shown in Table 3. Among them, the data acquisition instrument uses a high-precision distributed acquisition instrument developed by Beijing Dongfang Vibration Research Institute. The device has Ethernet and WIFI interfaces, supports multiple synchronous cascades, and can perform data acquisition through DASP software. The sensors used are three-directional vibration acceleration sensors, noise sensors produced by Beijing Dongfang Vibration Research Institute, and open-loop Hall current sensors produced by Beijing Senshe Electronics Co., Ltd.

According to the historical fault statistics of feed system, the fault frequency of rolling bearing is the highest, accounting for 42% of all faults, and the fault frequency of ball screw pair is the second, accounting for 26% [28].Therefore, in order to collect the data of common fault types of rolling bearings and ball screw pairs, the tools such as file and electric grinding needle are used to make different degrees of wear or damage scars on the inner and outer rings of bearings and the raceways of screws, and the bearing balls are polished with sandpaper to make wear faults. Figure 5 shows the tools used and some of the manufactured fault parts.

In this experiment, the health data and fault data of three common feeding conditions were collected respectively. The feed rates of the cases 1 to 3 are set as 1000 mm/min, 2000 mm/min and 3000 mm/min, respectively. Fault types include: bearing inner ring fault, bearing outer ring fault, bearing ball fault, screw wear, screw bending, screw wear and bearing inner ring composite fault, screw wear and bearing outer ring composite fault, and screw wear and bearing ball composite fault. The collected signals include vibration signal, noise signal and current signal. The sampling frequency is 10 kHz, and the sampling time for each fault is 120 s. The fault data set divided by 2000 data points per sample is shown in Table 4.

4.2. Signal Preprocessing

(1) Elimination of singular point

By setting the upper and lower threshold limits for the signal, the abnormal values outside the threshold range are eliminated. The empirical formula for the upper and lower limits of the threshold is the signal mean ± 4 signal standard deviations. Taking the noise sensor signal shown in Figure 6 (a) as an example, the calculated upper and lower threshold values are 5 and−5, respectively. The signal after removing the singular points is shown in Figure 6 (b).

(2) Elimination of trend term

In order to ensure the accuracy of the original data as much as possible, the signal trend line is fitted by the least square method and subtracted. Figure 7 (a) and (b) show the comparison of the X-direction vibration signals before and after the removal of the trend item.

(3) Wavelet threshold denoising

Wavelet threshold denoising is a nonlinear denoising method based on Wavelet transform. It is very suitable for processing non-stationary fault signals of CNC machines. In industrial signals, the fault signal mostly exists in the low-frequency component of the signal, while the noise is usually a high-frequency signal with a small amplitude [29]. The process of Wavelet threshold denoising is shown in Figure 8.

Sym5 is selected as the Wavelet base for signal denoising, and the original signal is decomposed by three-layer Wavelet. Then the soft and hard threshold compromise method is used for noise reduction, and the expression of the threshold function is shown as formula 18.

$${\sigma }_{\lambda }(\omega )=\{\begin{array}{c}\mathrm{sgn}(\omega )·(\left|\omega \right|-\alpha \lambda ),\left|\omega \right|\ge \lambda \\ 0,\left|\omega \right|<\lambda \end{array}$$

(18)

Where: $\omega$ is the Wavelet coefficient; $\lambda$ is the threshold; $\alpha$ is the scaling factor, The value of $\alpha$ in this paper is 0.5.

Figure 9 shows the comparison between the original vibration signal and the signal after the above wavelet threshold denoising. It can be seen that this method effectively eliminates the high-frequency noise while retaining the main characteristic information of the original signal, and the denoising effect is good.

4.3. Signal Feature Extraction

(1) Time domain feature extraction

In order to reflect the overall situation of the signal and reflect the sudden change in the signal, 13 time domain characteristic indexes with dimension and non-dimension are extracted, as shown in Table 5. In the table, $\left\{x_i\right\}$ is the discrete signal, and $i=1,2,\cdot \cdot \cdot ,N$ ,N is the number of sampling points.

(2) feature extraction in frequency domain

Spectrum analysis can reflect the distribution and change of frequency components in the signal and provide effective fault information in the signal. The three extracted frequency domain characteristic indexes and their calculation are shown in Table 6.

(3) feature extraction in time-frequency domain

CEEMDAN is used to decompose the preprocessed signal to extract and select the effective IMF components, and then the information entropy of IMF components, such as energy, power spectrum and singular spectrum, is calculated. Taking the X-direction vibration signal of bearing ball wear fault as an example, the result of CEEMDAN decomposition is shown in Figure 10. The correlation coefficient between each IMF component and the original signal is shown in Table 7, and the correlation coefficient threshold can be calculated as 0.178 according to formula (9). Therefore, IMF 1, IMF 9 and IMF 10 will be removed, and IMF 2 ~ IMF8 will be retained as effective IMF components for subsequent signal feature extraction.

4.4. Experimental Environment, Hyper-Parameter Setting and Model Evaluation Index

(1) Experimental environment configuration

The experiment uses a self-configured server with Intel core i9 11900k CPU, 128GB running memory and 64-bit Windows 10 operating system. The development environment is LightGBM 3.2.1.99, Python 3.8.

(2) Hyperparameter setting

The training hyperparameters of the DoubleEnsemble-LightGBM fault diagnosis model are set as follows: the number of iterations (num_iterations) is 100,The learning_rate is 0.14, the maximum depth of the decision tree (Max _ depth) is 7, the number of leaf nodes (num _ leaves) is 21,The minimum sample number of leaf nodes (min_data _in_leaf) is 30, the number of sub-models is 5, and the weight of sub-models is (1,1,1,1,1).The number of sample subsets is 4, the feature sampling ratio is 80%, and the loss function is the classification cross entropy loss.

(3) Model evaluation index

Confusion matrix [30] is often used to judge the performance of multi-classification models. Table 8 shows the confusion matrix of the fault category prediction results. Where the number in the main diagonal position indicates the number of samples that the model correctly classifies for each fault, Larger indicates better model diagnostic performance. The numbers in the remaining positions represent the number of misclassified samples, and the smaller the number, the better the diagnostic performance of the model. Through the confusion matrix, it can be clearly distinguished which kind of fault is easy to be confused by the model.

The overall diagnosis accuracy and individual diagnosis accuracy are used as the evaluation indexes of the fault diagnosis model. The overall diagnostic accuracy can reflect the overall diagnostic performance of the model, as calculated in equation (19). The individual diagnostic accuracy can reflect the diagnostic performance of the model for a specific type of fault, as calculated in equation (20).

$$T=\frac{{\displaystyle \sum _{i=1}^N{a}_{ii}}}{{\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=1}^N{a}_{ji}}}}$$

(19)

$$I_i=\frac{a_{ii}}{{\displaystyle \sum _{j=1}^N{a}_{ji}}},i=1,2,3,\cdot \cdot \cdot N$$

(20)

Where: $T$is the overall diagnostic accuracy rate; $I_i$ is the individual diagnostic accuracy rate; $a_{ji}$is the element value of the ith column of the jth row in the confusion matrix.

4.5. Analysis of Experimental Results

4.5.1. Analysis of Experimental Results of Public Data Set

Considering the influence of random factors on model training and testing, 10 repeated experiments were carried out. Figure 11 shows the confusion matrix of the last experimental test result.

The overall diagnostic accuracy and individual diagnostic accuracy of the DoubleEnsemble-LightGBM model under each speed condition are calculated by analyzing the confusion matrix, and the calculation results are shown in Table 9.It can be seen from the table that the overall diagnostic accuracy of the model is 90.96% after averaging the calculation results under four speed conditions, which can achieve better overall diagnostic performance.The individual diagnostic accuracy of category 1 to category 5 is 96.46%, 91.88%, 88.54%, 87.92% and 90% respectively. It can be seen that the diagnostic accuracy of the model for category 1 (health) is the highest, and the diagnostic accuracy for category 3 (bearing ball failure) and category 4 (bearing outer ring failure) is lower.

In addition, the diagnosis performance of the constructed DoubleEnsemble-LightGBM model is compared with that of the original LightGBM model and three other ensemble learning models with excellent performance in the field of fault diagnosis. The RF model used in [31], the AdaBoost model used in [32], and the XGBoost model in [33].The average value of the overall fault diagnosis accuracy of 10 experiments is taken as the evaluation index, and the experimental comparison results are shown in Table 10.It can be seen from the table that the average overall diagnostic accuracy of DoubleEnsemble-LightGBM model is the highest, which is increased by 6.57%, 6.61%, 3.42% and 4.06% respectively compared with RF model, AdaBoost model, XGBoost model and LightGBM original model. Figure 12 shows the comparison of the overall diagnostic accuracy of the five models under different speed conditions. The diagnostic performance of the Double Ensemble-Light GBM model is significantly better than that of other models.

4.5.2. Analysis of Experimental Results of Feed System Test Bench Data Set

Divide the feed system fault data set established by the feed system test bench in Section 4.1.2 into the training set and the test set at a ratio of 8:2. The distribution of the divided samples and the corresponding relationship of the fault labels are shown in Table 11.

In order to ensure the reliability of the model, 10 repeated experiments were also carried out. Figure 13 shows the confusion matrix for the last experimental test result.

The overall diagnostic accuracy and individual diagnostic accuracy of the DoubleEnsemble-LightGBM model under each feed condition are calculated by analyzing the confusion matrix, and the calculation results are shown in Table 12.In the table, the feed speed corresponding to working condition 1, working condition 2 and working condition 3 is 1000 mm/min, 2000 mm/min and 3000 mm/min respectively. It can be seen from the table that after averaging the calculation results under the three feeding conditions, the overall diagnostic accuracy of the model is 98.06%, and the individual diagnostic accuracy of categories 1 to 9 is 100%, 97.78%, 98.06%, 95%, 99.45%, 95.55%, 98.61%, 99.17% and 98.89% respectively. The results show that the DoubleEnsemble-LightGBM model can achieve high-precision fault diagnosis, and the classification accuracy of health data (class 1) reaches 100%.

In addition, RF model, AdaBoost model, XGBoost model and LightGBM original model are also selected to compare the diagnostic performance with DoubleEnsemble-LightGBM model. The average value of the overall fault diagnosis accuracy of 10 experiments is taken as the evaluation index, and the experimental comparison results are shown in Table 13. It can be seen from the table that compared with the original LightGBM model, the average overall diagnostic accuracy of the constructed DoubleEnsemble-LightGBM model is improved by 2.91% under three feeding conditions, indicating that the introduction of sample reweighting and feature selection mechanism can effectively improve the overall diagnostic performance of the model. Compared with RF model, AdaBoost model and XGBoost model, the average overall diagnostic accuracy of DoubleEnsemble-LightGBM model is still the highest, which is improved by 4.48%, 3.87% and 2.66%, respectively. Figure 14 shows more intuitively the comparison of the overall diagnostic accuracy of the five models at different feed rates. The diagnostic performance of the Double Ensemble-Light GBM model is significantly better than that of the other models.

5. Conclusions and Future Work

In order to solve the problem of intelligent fault diagnosis of CNC machine feed system under variable speed conditions, a variety of signals such as current signal, vibration signal and noise signal are used as monitoring data. Firstly, the above signals are preprocessed by singularity elimination, trend item elimination, Wavelet threshold de-noising. Then, time-domain analysis and frequency-domain analysis are carried out for each signal, and 13 time-domain characteristic indexes and 3 frequency-domain characteristic indexes are extracted. The time-frequency analysis of the signal is carried out through CEEMDAN, and three IMF information entropy are calculated. The multi-dimensional mixed domain feature set is constructed by stitching the above multiple feature indexes into feature vectors. Finally, LightGBM is selected as the basic fault diagnosis model. In addition, in order to further improve the training performance of the model and improve the diagnosis accuracy, the sample reweighting mechanism based on learning trajectory and the feature selection mechanism based on shuffling technology are introduced to build a DoubleEnsemble-LightGBM fault diagnosis model. The experimental results show that the average diagnostic accuracy of the DoubleEnsemble-LightGBM model is 91.07% on the public variable speed bearing fault data set, and 98.06% on the self-built fault data set of the feed test bench. Compared with the RF, AdaBoost, Xgboost and other advanced ensemble learning models and the original LightGBM model, the proposed DoubleEnsemble-LightGBM model effectively improves the diagnostic accuracy on both data sets.

Due to the limitation of experimental conditions, the fault data of the key mechanical components of the feed system are mainly collected by building a feed test bench and artificially making simulated faults. The follow-up research will try to accumulate the real fault data of the actual work and production of the CNC machine feed system, and further expand the types of faults.