Chatter Detection in Robotic Milling Using Entropy Features

Abstract

Chatter detection in robotic milling is a difficult issue due to the complex dynamic behavior of robots. In this paper, a novel approach to detecting chatter in the robotic milling process is proposed. The method of improved complete ensemble empirical mode decomposition with adaptive noise is introduced for decomposing the milling vibration signals into a series of intrinsic mode functions (IMFs). The effective IMFs are chosen according to the correlation between the original signals and each IMF. Signal reconstruction is conducted using the selected IMFs. The weighted refined composite multiscale dispersion entropy is extracted from the reconstructed signals in order to characterize the chatter states. Then, a classification model is established for chatter detection. Experimental results prove that the proposed method is feasible for chatter detection in the robotic milling process under different robot configurations and machining parameters.

1. Introduction

Industrial robots are essential for manufacturing. Robotic machining is a promising technique and attracts increasing interest [1,2,3,4]. In particular, robot milling is a suitable way to deal with large parts thanks to the high dexterity and wider working space of robots. However, the rigidity of industrial serial robots is much lower than conventional machine tools [5]. Chatter is more likely to occur in the robotic milling process.

When chatter occurs, the machined surface becomes rougher and the machining accuracy cannot be guaranteed [6]. Chatter has been regarded as a bottleneck that significantly limits the spread of robotic milling techniques. In the past few decades, the existing body of research has concentrated on revealing the mechanism of chatter generation during robotic milling. Pan et al. [7] suggested that mode coupling chatter is an important form of chatter during robotic milling. Cordes et al. [8] reported that chatter is caused by the regeneration effect.

However, the reasons for chatter in robotic milling can be regarded as a complicated interaction behavior containing regeneration effects and mode coupling effects [9]. Investigators have proposed several methods aimed at offline chatter suppression, such as optimizing processing parameters [10,11] and optimizing configuration [12,13]. Nevertheless, it is almost impossible to quantify the dynamic characteristics of a robot system precisely due to the existence of multiple flexible bodies. The results of offline chatter suppression may not be effective. Therefore, it is necessary to conduct online detection of chatter in robotic milling, and subsequently, a vibration control method can be carried out [14,15,16].

The issue of chatter detection in the conventional machining process has attracted a lot of interest [17,18]. However, significant differences can be found in dynamic properties between the milling process with robots and conventional machine tools. In robotic milling, the configurations of robots can apparently influence the generation of chatter, so chatter has time-varying characteristics along a machining trajectory. Existing studies in the field of chatter detection mainly focus on milling with conventional machine tools, which may be not applicable for chatter detection in robotic milling.

Until now, far too little attention has been paid to chatter detection in the robotic machining process. Cen et al. [19] developed a novel approach for real-time chatter monitoring using discrete wavelet transform. The chatter fault threshold was determined based on the standard deviation of the wavelet signals. Tao et al. [20] proposed a new method for chatter detection in robotic drilling combining the notch filter and local maximum synchrosqueezing transform (LMSST). Chen et al. [21] used entropy features extracted from adaptive variational mode decomposition (VMD) for chatter detection in robotic milling. Wang et al. [22] proposed a chatter detection model based on VMD, information entropy, and SVM.

However, there are some deficiencies in the existing works. Firstly, fast and reliable detection of chatter still remains a challenge. Secondly, fixed thresholds of chatter indicators are commonly set for detecting chatter [19,23,24,25,26]. In practice, cutting parameters and robot configurations may be different with the change in machining conditions. The applicability of using a manual setting threshold for variable machining conditions requires further validation.

Taking into consideration the time-varying characteristics of chatter, time-frequency analysis methods show advantages in online chatter detection. Typical time-frequency domain analysis methods, such as VMD methods [21,27], suffer the complicated process of parameter optimization. Ensemble empirical mode decomposition (EEMD) [25,28,29] has the shortcoming of residual noise and difficulty in averaging operations with different numbers of IMFs. To overcome the shortcomings of EEMD and improve the decomposition performance, several improved methods based on EEMD have been proposed [30,31,32]. Among them, the improved complete ensemble empirical mode decomposition with adaptive noise (ICEEMDAN) method has a number of attractive features, including less residual noise, less spurious modes, and more physical meanings. Recently, the ICEEMDAN method has proven to be feasible for fault diagnosis [33,34].

A set of IMFs can be obtained from the ICEEMDAN method. Signal reconstruction can be conducted using IMFs. However, the reconstructed signals cannot be directly used for chatter detection. Appropriate features should be extracted and used as chatter indicators. Fortunately, entropy, which is a measure of the uncertainty and irregularity of the data series, has been widely used as a chatter indicator [35,36,37,38]. Specifically, dispersion entropy (DE) is another algorithm for evaluating and quantifying the randomness and complexity of one-dimensional time series, which is highly sensitive to signal changes [39,40]. The above entropy features are commonly extracted under single-scale entropy. However, due to the strong nonstationary characteristics in vibration signals, the vibration modes are distributed among multi-scale. Therefore, useful information about chatter will be lost when only single-scale entropy is extracted. To address this issue, the refined composite multiscale dispersion (RCMDE) method has been proposed for estimating the complexity of time series under multi-scale [41].

In this work, the ICEEMDAN method is introduced for decomposing the milling vibration signals into a series of IMFs. Signal reconstruction is performed using IMFs, which are highly correlated with original signals. Weighted refined composite multiscale dispersion entropy (wRCMDE) is extracted from obtained reconstructed signals and acts as the chatter indicator. Then, a classification model is established for chatter detection. The remainder of this paper is organized as follows: In Section 2, the theoretical basis and framework for the chatter detection system are given. In Section 3, the experimental setup for robotic milling is presented, and an analysis of the vibration signals is performed. In Section 4, the proposed method is verified, and in Section 5, conclusions are drawn.

2. Theoretical Basis and Framework for the Chatter Detection System

2.1. Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (ICEEMDAN)

The ICEEMDAN method, which is based on the EMD decomposition principle, can further reduce the residual noise in the intrinsic mode functions (IMFs) and improve the decomposition effect. Considering the original signal, x, the main idea of the ICEEMDAN method can be expressed as follows [31]:

Step 1: White Gaussian noise is added to the original signal, x, and X₁⁽ⁱ⁾ = x + β₀E₁(w⁽ⁱ⁾) (i = 1, 2, 3, …, N) indicates the first set of the realization of the signal-plus-noise, where w⁽ⁱ⁾ denotes the noise, E_k(·) denotes the kth mode obtained from EMD, and β₀ = ε₀ σ(x)/σ(E₁(w⁽ⁱ⁾)) means a desired signal–to–noise ratio (SNR) at first decomposition. Then, the first residue, s₁, can be calculated by:

$$s_1=\langle M\left(X_1{}^{\left(i\right)}\right)\rangle =\langle X_1{}^{\left(i\right)}-E_1\left(X_1{}^{\left(i\right)}\right)\rangle \left(i=1,2,3,\dots ,N\right)$$

(1)

where M(·) denotes the local mean of the signal and 〈·〉 denotes the action of averaging throughout the realizations.

Step 2: The first mode can be calculated by:

I₁ = x − s₁.

(2)

Step 3: For k ≥ 2, constructing the kth set of signal-plus-noise as X_k⁽ⁱ⁾ = s_k−1 + β_k−1E_k(w⁽ⁱ⁾), where β_k−1 = ε₀ σ(x)/σ(s_k−1). Therefore, the kth residue, s_k, can be estimated as:

$$s_k=\langle M\left(X_k{}^{\left(i\right)}\right)\rangle =\langle X_k{}^{\left(i\right)}-E_k\left(X^{\left(i\right)}\right)\rangle \left(i=1,2,3,\dots ,N\right)$$

(3)

Step 4: The kth mode can be calculated by:

I_k = s_k−1 − s_k

(4)

Step 5: Repeat the process for the next k, beginning at Step 3.

2.2. Weighted Refined Composite Multiscale Dispersion Entropy (wRCMDE)

2.2.1. Dispersion Entropy (DisEn)

Considering the given time series {u(i), i = 1, 2, 3, …, L} with length, L, DisEn can be calculated as follows:

Step 1: Each element, u(i), is mapped into c classes. To be more specific, u(i) is first normalized to {y(i), i = 1, 2, 3, …, L}, y(i) ∊ [0, 1] using the normal cumulative distribution function, and y(i) is transformed into an integer from 1 to c using:

$$z_i^c=round\left(c\times y\left(i\right)+0.5\right)$$

(5)

Step 2: Defining embedding vectors as:

$$z_i^{m,c}=\left\{z_i^c,z_{i+\tau }^c,\dots ,z_{i+(m-1)\tau }^c\right\}i=1,2,\dots ,L-\left(m-1\right)\tau$$

(6)

where m is the embedding dimension and τ is the time delay. Each embedding vector corresponds to a dispersion pattern, ${\pi }_{v_0{v}_1\dots v_{m-1}}$, where $v_0=z_i^c,v_1=z_{i+\tau }^c,\dots ,v_{m-1}=z_{i+(m-1)\tau }^c$. The total number of possible dispersion patterns is c^m.

Step 3: The relative probability of ${\pi }_{v_0{v}_1\dots v_{m-1}}$ is calculated by:

$$p\left({\pi }_{v_0{v}_1\dots v_{m-1}}\right)=\frac{\left\{i|i\le L-\left(m-1\right)\tau ,z_i^{m,c}\mathrm{has}\mathrm{type}{\pi }_{v_0{v}_1\dots v_{m-1}}\right\}}{L-\left(m-1\right)\tau }$$

(7)

Step 4: DisEn can be calculated by:

$$DisEn\left(u,m,c,\tau \right)=-{\displaystyle \sum _{\pi =1}^{c^m}p({\pi }_{v_0{v}_1\dots v_{m-1}})}\ln \left(p({\pi }_{v_0{v}_1\dots v_{m-1}})\right)$$

(8)

DisEn is a statistical measurement that reflects the complexity of the time series. In general, the more complexity a time series has, the larger the DisEn value.

2.2.2. Multiscale Dispersion Entropy (MDE)

To overcome the shortcoming of unreasonable results when using DisEn, MDE is proposed for further characterizing the randomness of the time series.

MDE is calculated as follows:

Step 1: The “coarse-graining” process is applied to the time series {u(i), i = 1, 2, 3, …, L}, which is given by:

$$y_j^{\left(s\right)}=\frac{1}{s}{\displaystyle \sum _{i=\left(j-1\right)s+1}^{js}u\left(i\right)},1\le j\le \left[\frac{L}{s}\right]$$

(9)

where s denotes the scale factor and [·] indicates the rounding operation.

Step 2: The DisEn is calculated for each ‘‘coarse-grained” time series under the same parameters with different scale factors, s.

$$MDE\left(u,m,c,\tau ,s\right)=DisEn\left(y^s,m,c,\tau \right)$$

(10)

2.2.3. Weighted Refined Composite Multiscale Dispersion Entropy (wRCMDE)

To overcome the shortcoming of inaccurate results when using MDE for short signals, RCMDE is proposed, which is calculated as follows [41]:

Step 1: The “coarse-graining” process of RCMDE can be expressed by:

$$y_{k,j}^{\left(s\right)}=\frac{1}{s}{\displaystyle \sum _{b=k+\left(j-1\right)s}^{k+js-1}u\left(b\right)},1\le j\le \left[\frac{L}{s}\right],1\le k\le s$$

(11)

Step 2: For each scale factor, s, the DisEn of all of the obtained time series $y_k^{\left(s\right)},1\le k\le s$ is calculated, and the RCMDE is given by:

$$RCMDE(u,m,c,\tau ,s)=\frac{1}{s}{\displaystyle \sum _{k=1}^{s}DisEn\left(y_k^{\left(s\right)},m,c,\tau \right)}$$

(12)

Furthermore, kurtosis of the time series, KT(u), is given by:

$$KT\left(u\right)=\frac{\frac{1}{L}{\displaystyle \sum _{i=1}^L{\left(u(i)-\overline{u}\right)}^4}}{{\left(\frac{1}{L}{\displaystyle \sum _{i=1}^L{\left(u(i)-\overline{u}\right)}^2}\right)}^2}$$

(13)

is used for weighting the RCMDE. Therefore, the defined wRCMDE is calculated by:

$$wRCMDE\left(u,m,c,\tau ,s\right)=KT\left(u\right)\times RCMDE\left(u,m,c,\tau ,s\right)$$

(14)

Step 3: s = s + 1; repeating steps (2) and (3) under the same parameters until s = s_m, where s_m is the largest scale factor.

2.3. Framework for the Chatter Detection System for Robotic Milling

The framework for the proposed chatter detection system is shown in Figure 1. It consists of three steps:

Step 1: Acquiring online data

Extensive experiments is conducted to collect the vibration signals in robotic milling under variable robot configurations and cutting parameters.

Step 2: Extracting features

The ICEEMDAN method is performed on the vibration signals, and a set of IMFs is obtained. The effective IMFs are selected according to the correlation between the original signals and each IMF. Subsequently, the wRCMDE is extracted from the reconstructed signals using the selected IMFs.

Step 3: Constructing classification model

Based on the extracted wRCMDE features, the classification model can be constructed and trained for chatter detection in robotic milling.

3. Experimental Verification and Analysis

3.1. Robotic Milling System and Cutting Parameters

The experimental setup is illustrated in Figure 2. In the present study, milling tests, under variable robot configurations and cutting parameters, were performed on a robotic milling system. The system consisted of a Comau industrial robot. An electric spindle (Jäger, F120-H824) acted as the spindle unit. The vibration signals were measured using a PCB 356A15 acceleration sensor and recorded by a computer. The sampling rate was 5120 Hz. Φ 10 mm milling tools and aluminum alloy workpieces were used for tests. Two configurations were chosen for the experimental verification, as illustrated in Figure 3. The milling parameters are given in

Table 1. Milling parameters and the chatter state.

No.	Configuration	n (r/min)	ap (mm)	ae (mm)	Chatter State
1	I	3000	1	8	Chatter
2	I	4000	1	8	Chatter-free
3	I	6000	1	8	Chatter-free
4	I	8000	1	8	Chatter-free
5	I	3000	5	1	Chatter
6	I	3000	5	2	Chatter
7	I	3000	5	3	Chatter
8	I	4000	5	2	Chatter
9	I	6000	5	1	Chatter-free
10	I	6000	5	2	Chatter
11	I	8000	5	1	Chatter-free
12	II	3000	1	8	Chatter
13	II	6000	1	8	Chatter-free
14	II	3000	5	1	Chatter-free
15	II	3000	5	2	Chatter
16	II	3000	5	3	Chatter
17	II	4000	5	3	Chatter-free
18	II	6000	5	3	Chatter

. The feed speed was fixed at 2 mm/s. It is noted that the different tests are defined as configuration−n−a_e−a_p.

3.2. Vibration Signal Analysis

The frequency domain features can be clearly revealed by exploiting fast Fourier transform (FFT) [42]. In this work, vibration signals in the machining direction are used for the subsequent analysis. Generally speaking, in the stable cutting process, the spectrum peaks of FFT results mainly appear at the spindle rotation frequency (SRF), which is defined as n/60, and its integer multiples. When chatter occurs, the spectrum peaks will appear at chatter frequencies. The representative FFT results are shown in Figure 4. Test No. 7 (I−3000−5−3) is unstable because several peaks of the spectrum can be seen in the low-frequency range (16.25/140.9/165.3/177.2 Hz), and the peaks of the spectrum are not integer multiples of SRF. Test No. 10 (I−6000−5−2) is unstable because obvious peaks of the spectrum appear at 2076 and 2110 Hz. Test No. 12 (II−3000−1−8) is unstable because the peaks of the spectrum clearly appear at 2062 Hz. Test No. 14 (II−3000−5−1) is stable because only SRF and its integer multiples can be seen in the frequency domain plots.

4. Results and Discussion

4.1. Features Extraction

Several parameters should be manually set before performing ICEEMDAN. The parameter ε₀ affects the desired SNR. According to [31], ε₀ should be set as 0.2. Then, the parameter N determines the averaging number. If N is large, too many averaging operations need to be conducted, which is a time cost. If N is small, reasonable results cannot be obtained. For fast detection of chatter, we used N = 5. A set of IMFs can be obtained from the ICEEMDAN. Further, only IMFs which are highly correlated with the original signals contain useful information about chatter. Therefore, correlation analyses between original signals and each IMF needed to be conducted. The IMFs with correlation coefficients greater than 0.1 were taken into consideration and used to reconstruct the original vibration signals.

Next, the wRCMDE could be extracted from the reconstructed signals. Further, it was critical that some parameters were predetermined. The number of classes, c, was the first parameter of focus. When c is too small, two data points with large differences in amplitude may be classified into the same class, which leads to underestimating the irregularity of the signal. Conversely, when c is too large, the estimation results are easily disturbed by noise and the robustness is poor. To obtain an ideal estimation result, c was set to 5 in this paper based on the existing research [40]. The embedding dimension, m, was the second parameter to be considered. With the decrease in the embedding dimension, m, the embedding vectors become shorter, and more embedding vectors can be obtained, which cannot accurately reflect the overall dynamic changes in the signals. However, fewer embedding vectors would be obtained with an increase in m, which may result in the loss of detailed information for the time series. According to [43], m was set to 5 in this paper. Commonly, τ is set as 1 in practice. The largest scale factor, s_m, is usually set at s_m ≥ 10 [43]. However, unstable estimation results would be observed if s was too large [44]. Consequently, we chose s_m = 15.

The wRCMDE was extracted from the reconstructed signals corresponding to all of the robotic milling tests. The typical wRCMDE values under different robot configurations and cutting parameters are shown in Figure 5. It can be seen from Figure 5 that the wRCMDE values showed a downward trend with the increase in the scale factor. It is worth noting that the wRCMDE values obtained from the chatter state were larger than those obtained from the chatter-free state when the scale factor was s ≤ 4. Therefore, wRCMDE values obtained at s ≤ 4 were chosen as the features for constructing the chatter detection model.

4.2. Classification Models Construction

In this paper, classification models were developed using the “classification learner” toolkit provided by MATLAB. Several classical classification learner models were trained, including support vector machine (SVM), k-nearest neighbor (KNN), and decision tree. A total of 540 groups of wRCMDE samples for the two chatter states (chatter and chatter-free) were selected as training samples. The wRCMDE samples were be normalized first, and 5-fold cross-validation was adopted to modify the classification accuracy.

During the training of each classification model, the kernel parameters were optimized based on the training accuracy. The training accuracy of the selected models is listed in

Table 2. Training accuracy of classification models.

Classification Model	Parameters		Training Accuracy
SVM	Kernel function	Gaussian	95.2%
	Kernel scale	1.45
KNN	Number of neighbors	5	95.0%
	Distance metric	Chebyshev
Decision Tree	Max number of splits	7	93.9%
	Split criterion	Gini’s diversity index

. The SVM model with the Gaussian kernel function was chosen as the final classification model.

Furthermore, the SVM models based on ICEEMDAN and RCMDE (ICEEMDAN-RCMDE), ICEEMDAN and MDE (ICEEMDAN-MDE), EEMD and wRCMDE (EEMD-wRCMDE), EEMD and RCMDE (EEMD-RCMDE), and EEMD and MDE (EEMD-MDE) were also constructed for comparison. It can be seen from

Table 3. Training accuracy of different methods.

Method	Training Accuracy
ICEEMDAN-wRCMDE (proposed method)	95.2%
ICEEMDAN-RCMDE	91.9%
ICEEMDAN-MDE	90.0%
EEMD-wRCMDE	90.6%
EEMD-RCMDE	87.8%
EEMD-MDE	86.9%

that the proposed model (named ICEEMDAN-wRCMDE) showed superiority over the other methods.

Further, to show the performance details of the constructed SVM model based on ICEEMDAN-wRCMDE, the parallel coordinates plot, the receiver operating characteristic (ROC) plot, and the confusion matrix are presented in Figure 6, Figure 7 and Figure 8, respectively.

4.3. Classification Performance

During the chatter detection process, the signals used for processing each time contain 1024 sample points. To verify the effectiveness of the proposed method, chatter detection was conducted by checking the classification results over time. Meanwhile, the vibration signals were filtered to eliminate the SRFs and processed using the short-time Fourier transform (STFT) to illustrate the time–frequency spectrogram. The original vibration signals, chatter classification results, and time–frequency spectrograms of each experiment are shown in Figure 9 and Figure 10. It is noted that label 1 denotes a chatter state and label 0 denotes a chatter-free state. The test groups in Figure 9 come from the experiments listed in Table 1. Furthermore, the test groups in Figure 10 are additional experiments that differ from those listed in Table 1.

It can be seen from Figure 9 and Figure 10 that the proposed method can realize fast and reliable chatter detection with different robot configurations and cutting parameters. Only intermittent false results can be observed, which does not affect the final decision of the chatter state and is acceptable for practice.

Concerning the influence of the acceleration sensor position and location on the raw data acquisition, additional experiments were conducted. Configuration III denotes another robot configuration (γ₁ = −10°, γ₂ = 22°, γ₃ = −164°, γ₄ = −10°, γ₅ = −74°, γ₆ = 1.5°). The position of two acceleration sensors is shown in Figure 11a, where acceleration sensor 1 is mounted on the spindle and acceleration sensor 2 is mounted on the robot end. The comparison of time-domain signals obtained from the two sensors is shown in Figure 11b. As can be seen, the vibration signals obtained from acceleration sensor 1 show larger amplitudes than acceleration sensor 2, due to acceleration sensor 1 being located closer to the cutting area. Furthermore, as can be seen from Figure 12, the results of the wRCMDE values are similar to those in Figure 5; this indicates that the acceleration sensor position and location have little effect on the proposed chatter detection method.

5. Conclusions

A novel approach to detecting chatter in the robotic milling process is proposed. The ICEEMDAN method was used to decompose the milling vibration signals into a series of IMFs. Based on the correlations between the original signals and each IMF, the effective IMFs were obtained. Signal reconstruction was conducted using the obtained effective IMFs. Weighted refined composite multiscale dispersion entropy was extracted from the obtained reconstructed signals to characterize the chatter states. Then, a classification model based on SVM was established for chatter detection. Extensive case studies prove that the proposed method is feasible for chatter detection in the robotic milling process under different robot configurations and cutting parameters.