Milling chatter detection based on an optimized iTransformer-BiGRU-random forest ensemble model

Chatter represents the most detrimental manifestation of self-excited vibrations in milling operations, precipitating severe degradation of workpiece surface integrity, accelerated tool wear, and potentially catastrophic equipment failure, thereby substantially compromising machining efficiency and product quality. To address this critical challenge, researchers have predominantly pursued two complementary avenues of investigation: chatter prediction [1] and chatter detection [2]. The predictive approach leverages stability models to facilitate the judicious selection of machining parameters. However, owing to the inherent complexity and variability of practical machining conditions, chatter phenomena may still materialize even with theoretically optimized parameters. In contrast, chatter detection methodologies enable real-time monitoring of the machining state with instantaneous warning capabilities, garnering considerable attention within the research community due to their superior temporal responsiveness and diagnostic accuracy.

Signal acquisition constitutes the foundational step in chatter detection frameworks. Commonly employed sensing modalities encompass vibration signals [3], cutting force signals [4], and acoustic emissions [5], [6]. Empirical evidence substantiates that both cutting force and acceleration signals exhibit remarkable efficacy for chatter identification [7]. Nevertheless, subsequent investigations have revealed that acceleration signals demonstrate superior suitability for chatter detection applications relative to cutting force measurements, primarily attributable to their reduced computational complexity and implementation requirements [8]. Presently, acceleration-based sensing has emerged as the predominant measurement paradigm in chatter detection methodologies.

Chatter feature extraction constitutes the subsequent critical step in chatter detection. Prevalent methodologies encompass time-domain methods, frequency-domain methods, and time-frequency domain approaches. Time-domain features predominantly comprise dimensionless characteristics, such as root mean square [2], [9], [10], mean values [2], [9]. Gao et al. [9] formulated a composite time-domain feature by integrating mean value, standard deviation, and kurtosis. Frequency-domain features predominantly include spectral mean [2], [9] and spectral centroid frequency [2]. Time-frequency domain methodologies have garnered extensive adoption owing to their inherent capability of furnishing simultaneous temporal and spectral information, encompassing techniques such as wavelet energy entropy [2], continuous wavelet transform [11], and synchrosqueezing transform [12], [13].

Chatter detection classification algorithms predominantly encompass traditional machine learning paradigms and deep learning architectures. Among conventional machine learning approaches, support vector machines [14], [15] and random forests [5], [10] have achieved the most widespread deployment. Tran et al. [5] employed random forest algorithms to train multi-sensor fusion data for machining state identification. Yan et al. [10] leveraged random forests in conjunction with optimally selected feature subsets to accomplish intelligent diagnosis of machining stability. Nevertheless, these methodologies exhibit fundamental dependence on manual feature engineering, thereby constraining their capacity to autonomously extract deep nonlinear patterns and temporal dependencies.

In recent years, deep learning has attracted considerable scholarly attention by virtue of its formidable automatic feature learning and classification capabilities, progressively permeating the domain of chatter detection. Gao et al. [16] developed a hybrid neural network architecture integrating dual-scale parallel convolutional neural networks, BiGRU, and multi-head attention mechanisms for chatter detection. Unver et al. [17] proposed a transfer learning framework that synergistically combines analytical solutions with convolutional neural networks for chatter identification. Jauhari et al. [18] introduced a novel network architecture that amalgamates Stem modules, Inception modules, and embedded residual attention mechanisms for chatter detection. Despite demonstrating exceptional performance in feature representation, deep learning approaches remain susceptible to overfitting under small-sample scenarios and exhibit insufficient model interpretability.

Existing chatter detection methodologies are encumbered by three critical limitations. Firstly, prevailing research predominantly relies upon monolithic models, which struggle to maintain robust performance across complex and variable operational conditions. Secondly, although ensemble models have demonstrated substantial efficacy in other fault diagnosis domains [19], [20], their application in chatter detection remains notably limited. Thirdly, contemporary hyperparameter optimization techniques [21] are primarily tailored for individual models, lacking systematic mechanisms for the synergistic optimization of hyperparameters across multiple constituent models within ensemble architectures, as well as their corresponding weighting coefficients.

To address the aforementioned challenges, this study presents an optimized iTransformer-BiGRU-Random Forest (iTBU-RF) ensemble model for milling chatter detection. The principal contributions of this methodology encompass: (1) The simultaneous optimization of hyperparameters across both the iTBU model and the RF model through an integrated optimization algorithm, whereby a global search strategy facilitates the acquisition of optimal model configurations while circumventing the suboptimal solutions inherent in conventional stepwise optimization procedures. (2) The development of a dynamic weighting ensemble framework that synergistically integrates the iTBU model and the RF model, wherein the contribution weights of individual constituent models are adaptively modulated in accordance with varying operational conditions and data characteristics, thereby substantially enhancing the generalization capability and detection accuracy of the ensemble architecture.

The Ivy algorithm was proposed by Ghasemi et al. in 2024 [22]. This algorithm primarily emulates the distinct life stages of Ivy plants, including growth, climbing, and propagation within the Ivy population.

Upon algorithm initialization, the initial positions of the Ivy population within the search space are stochastically determined using Eq. (1).

where, $I_{\max }$ and $I_{\min }$ denote the upper and lower bounds of the search space, respectively. $\odot$ represents the Hadamard product of two vectors. ${rand}$(1, $D$) denotes a D-dimensional vector comprising uniformly distributed random numbers within the interval [0, 1]. $N_{pop}$ signifies the total population size of Ivy plants.

Drawing upon data-intensive experimental and simulation procedures, the difference equation governing the growth velocity $G v_i(t)$ of member $I_i$ is formulated as Eq. (2).

Vectors $\Delta G v_i(t + 1)$ and $\Delta G v_i(t)$ represent the growth rates of the discrete-time system at time instants $t$ and $t+1$, respectively. $rand^2$ denotes a random number from a random variable whose probability density function is equivalent to $1/(2 \sqrt{x}$). $N$(1, $D$) signifies a random vector of dimensionality $D$.

The equation governing the climbing and logical movement of member $I_i$ toward the light source direction via leveraging member $I_{i i}$ is expressed in Eq. (3).

where, $|N(1\text{,}\ D)|$ denotes the absolute value of $N(1\text{,}\ D)$. $\oslash$ represents the Hadamard division of two vectors.

Member $I_i$ explores the vicinity of member $I_{\text {Best }}$ in search of superior optimal solutions. The computational procedure is presented in Eq. (4).

When the objective function value $f\left(I_i\right)$ of member $I_i$ is less than a multiple of $f\left(I_{\text {Best }}\right)$, the parameter $\beta$ = (2 + $rand$)$/$2. Subsequently, the Ivy commences expansion of its branch and leaf width (as governed by Eq. (3)). Conversely, the Ivy undergoes upward growth and climbing (as prescribed by Eq. (4)).

The iTransformer model constitutes a refined neural network architecture predicated upon the Transformer framework [23]. In contrast to conventional Transformer models, iTransformer treats individual time series as tokens and leverages the self-attention mechanism to capture multivariate correlations. Furthermore, it employs layer normalization and feedforward neural network modules to more effectively acquire the global representation of sequences for time series forecasting applications. These architectural enhancements enable iTransformer to exhibit superior performance in processing long-sequence tasks with enhanced versatility. The specific architecture of the iTransformer model is illustrated in Figure 1.

The architecture of iTransformer encompasses local layer normalization, feedforward neural networks, and self-attention modules. The normalization operation is applied to the time series representations of variables, normalizing all tokens to a normal distribution, thereby facilitating the reduction of variance introduced by measurements and addressing the non-stationarity problem.

The multivariate attention mechanism ensures that highly correlated variables obtain elevated weighting values in subsequent representation interactions. This computational process involves calculating attention scores, scaling, applying the softmax function, and performing weighted summation. The formulation is expressed in Eq. (5).

where, $Q$, $K$, and $V$ denote the query vector, key vector, and value vector, respectively. $d_k$ represents the projection dimensionality.

The feedforward neural network constitutes another essential component within the iTransformer model, comprising linear transformations and nonlinear activation functions. Its computational procedure is formulated in Eq. (6).

where, $x$ denotes the input vector, $W_1$, $b_1$, $W$, and $b_2$ represent the learnable weight matrices and bias terms, $\text{max} (0\text{,}\ \cdots)$ signifies the ReLU activation function, and $FFN(x)$ designates the output of the feedforward neural network.

BiGRU addresses the propagation directionality issue by employing two GRU layers [24]. The forward GRU processes the input sequence from the initial to the terminal time step, while the backward GRU processes it in the reverse direction. Each GRU layer maintains its respective hidden state, updating it based on the current input and the hidden state from the preceding time step.

The computational formulation for the hidden layer units within BiGRU is expressed in Eq. (7).

where, $\overrightarrow{h}_t$ and $\overleftarrow{h}_t$ denote the forward and backward states of the hidden layer at time instant $t$, respectively. $W_{\overleftarrow{h}}$ and $W_{\overrightarrow{h}}$ represent the forward and backward weights of the hidden layer at time instant $t$, while $b_t$ signifies the bias term of the hidden layer at time instant $t$.

RF constitutes an ensemble learning methodology whose fundamental principle involves employing bootstrap aggregating (Bagging) to aggregate multiple weak classifiers into a robust ensemble classifier [25]. The ultimate prediction is determined through majority voting or averaging across independently distributed constituent weak classifiers.

The core mechanism of RF entails constructing binary decision trees through recursive binary partitioning of the data, wherein each internal node represents a feature attribute, and each terminal node corresponds to either a class label or a numerical value.

ITBU architecture demonstrates superior capability in capturing patterns across diverse temporal scales, thereby enhancing predictive accuracy and robustness. The Random Forest model, characterized by its expeditious training efficiency and formidable generalization capability, furnishes stable predictive performance on time-frequency domain features. Consequently, this study adopts an ensemble framework integrating the iTBU model and RF model for machining state identification.

Initially, the entire dataset is partitioned into training, validation, and test subsets. The iTBU model and RF model are independently trained on the training set, yielding machining state detection models $M_{\mathrm{iTBU}}$ and $M_{\mathrm{RF}}$, respectively. The validation set data are subsequently fed into models $f_{\mathrm{iTBU}}$ and $f_{\mathrm{RF}}$ to obtain accuracies $A_V^{\mathrm{iTBU}}$ and $A_V^{\mathrm{RF}}$. An optimization algorithm is employed to select optimal model hyperparameters, thereby acquiring the optimal predictive models $M_V^{\mathrm{iTBU}}$ and $M_V^{\mathrm{RF}}$ with corresponding maximum accuracies $A_{\max }^{\mathrm{iTBU}}$ and $A_{\max }^{\mathrm{RF}}$.

The ensemble weights are determined based on the validation set errors.

When a substantial accuracy discrepancy exists between the two models on the validation set, ensemble integration would compromise overall predictive accuracy. Consequently, an error threshold is established. If the accuracy differential between the two models exceeds the predefined threshold, the superior-performing model is designated as the final detection model. Conversely, if the accuracy disparity falls below the threshold, the models are integrated through weighted combination.

The detailed procedure is illustrated in Figure 2.

The validation set accuracy is employed as the fitness function for hyperparameter optimization of the model (Eq. (9)).

Performance evaluation metrics including accuracy, precision, and recall are adopted to assess the predictive capabilities of the model, as formulated in Eq. (10).

The machining operations were conducted on a VDL-1000E three-axis CNC machine tool manufactured by Dalian Machine Tool Group. A flat-end milling cutter with four flutes and a diameter of 10 mm was employed for the cutting trials. The workpiece material comprised TC4 titanium alloy with dimensions of 200 × 200 × 5 mm. Acceleration signals during the machining process were acquired using a PCB accelerometer featuring a sensitivity of 10.42 mV/g, coupled with a Donghua DH5922 data acquisition system operating at a sampling frequency of 5000 Hz. The milling configuration consisted of climb milling under dry cutting conditions. The experimental setup is depicted in Figure 3. The simulation experiments were executed on a Windows 10 (64-bit) operating system with a hardware platform comprising an Intel® Core™ i9-12900K CPU, NVIDIA GeForce RTX 3080 graphics processing unit, operating at 3.2 GHz clock frequency with 32 GB RAM. The computational environment was established using MATLAB R2023b.

The radial depth of cut was maintained at 0.5 mm, with a feed per tooth of 0.1 mm/tooth. The remaining cutting parameters, machining conditions, and chatter characteristics are comprehensively documented in Reference [21].

The Pearson correlation coefficient constitutes a widely adopted quantitative metric in feature sensitivity analysis, providing intuitive numerical characterization of the correlation strength between features and target variables [26]. However, this metric exclusively reflects the magnitude of association while lacking the capability to assess the statistical reliability of the observed correlations. The significance level, ascertained through rigorous statistical hypothesis testing, enables the discrimination of genuine correlations from those potentially arising from stochastic errors, thereby possessing well-defined statistical interpretation and effectively compensating for this inherent limitation of the Pearson correlation coefficient. Predicated upon these considerations, the present investigation synergistically employs both the Pearson correlation coefficient and significance level to conduct feature sensitivity analysis, thereby simultaneously achieving the dual objectives of correlation strength screening and correlation reliability validation. The resultant correlation coefficients and significance levels are illustrated in Figure 4.

As depicted in Figure 4, the Gravity frequency of frequency (GFFS) and Mean square frequency of frequency spectrum (MSFS) features exhibit negative correlations with the machining state, whereas the remaining time-frequency domain features demonstrate positive correlations. From the perspective of overall association strength across feature dimensions, temporal-domain features manifest the most pronounced correlation intensity, followed by frequency-domain features. In accordance with the correlation strength classification criteria established in the literature [27], correlation coefficients residing within the interval of 0.6 to 0.8 are designated as strong correlations. Should the Pearson correlation coefficient threshold be established above 0.7, approximately 41.7% of the strongly correlated features identified in this investigation would be eliminated, thereby incurring substantial loss of efficacious information. Considering both the definitional criteria for strong correlation and the imperative for information preservation, the present study establishes the Pearson correlation coefficient screening threshold at 0.6, concurrently mandating that correlations attain an extremely significant level. Only features simultaneously satisfying these dual criteria are retained; otherwise, they are discarded. Following the aforementioned correlation screening and significance validation procedures, the following features are ultimately identified as sensitive indicators of machining state: Standard deviation (STD), Root mean square (RMS), Absolute mean (ABS), Peak to peak value (PTP), Gravity frequency of frequency (GFFS), Wavelet packet energy entropy 3 (WPE 3), Mean square frequency of frequency spectrum (MSFS), Mean of frequency spectrum (MAF), Wavelet packet energy entropy 2 (WPE 2), Wavelet packet energy entropy 7 (WPE 7), Wavelet packet energy entropy 4 (WPE 4), and Wavelet packet energy entropy 5 (WPE 5).

Figure 5 illustrates the mean classification accuracy achieved via ten-fold cross-validation on the test set across varying threshold configurations.

As illustrated in Figure 5, the mean classification accuracy exhibits a progressive ascending trend as the threshold increases from 0.02 to 0.05, culminating in a peak value at the threshold of 0.05, whereupon it subsequently diminishes with further threshold augmentation. Concurrently, the standard deviation demonstrates a gradual reduction from thresholds 0.02 to 0.04, attaining its minimum at 0.04; a marginal increase manifests at the threshold of 0.05, though the magnitude remains at a comparatively modest level. Upon comprehensive evaluation, the threshold of 0.05 corresponds to the maximal mean accuracy among all investigated threshold values, whilst simultaneously maintaining the root mean square error at a consistently low magnitude, thereby ensuring both optimal overall classification accuracy and robust result stability. Consequently, the threshold value of 0.05 is identified as the optimal selection for the ensemble model.

The detailed specifications of iTBU-RF model parameters optimized by the Ivy algorithm, encompassing the parameter value ranges and optimization outcomes, are presented in Table 1.

As delineated in Table 1, the iTBU model encompasses 1, 128 trainable parameters, of which 400 constitute task-specific parameters, whilst the RF model comprises a total of 552 parameters. The dataset employed in the present investigation incorporates 2,880 training samples, 360 validation samples, and 360 testing samples. The ratio of sample size to the parameter count of the iTBU model approximates 7.2:1, whereas the corresponding ratio for the RF model stands at 5.2:1, both of which conform to the generally established requisites for deep learning model training [28].

Under the optimized hyperparameter configuration of the iTBU-RF model, the classification results for the training, validation, and test sets are illustrated in Figure 6.

As illustrated in Figure 6 (a), the classification accuracy achieved on the training set reaches 98.75%. Specifically, 6 instances of stable cutting were erroneously classified as slight chatter, 12 instances of slight chatter were misidentified as severe chatter, 12 instances of slight chatter were incorrectly recognized as stable cutting, and 6 instances of severe chatter were mistakenly categorized as slight chatter. Figure 6 (b) demonstrates that the validation set attains a classification accuracy of 96.67%, wherein 6 stable cutting cases were misclassified as slight chatter, and 6 slight chatter instances were erroneously identified as stable cutting. As depicted in Figure 6 (c), the test set exhibits a classification accuracy of 95%. The misclassification patterns encompass 6 stable cutting samples incorrectly recognized as slight chatter, 6 slight chatter instances misidentified as stable cutting, and 6 severe chatter cases erroneously categorized as slight chatter.

The ten-fold cross-validation results for machining state detection prior to and subsequent to hyperparameter optimization are presented in Table 2.

As evident from Table 2, the model detection performance exhibits comprehensive enhancement following hyperparameter optimization. The accuracies of the training, validation, and test sets are improved by 1.78%, 2.40%, and 1.78%, respectively. The optimized model demonstrates more balanced performance across different datasets, with the accuracy disparity between the training and test sets reduced from 2.94 percentage points prior to optimization to 2.40% points thereafter.

The ten-fold cross-validation results for machining state detection before and after feature selection are presented in Table 3.

As demonstrated in Table 3, the comparative analysis before and after feature optimization reveals that when the feature dimensionality is reduced from 20 to 12, the model exhibits marginal performance enhancement across all datasets. Specifically, the training set accuracy increases slightly from 98.69% to 98.76%, the validation set accuracy improves from 96.61% to 96.63%, while the test set accuracy demonstrates a modest increase from 95.15% to 95.16%. Notably, following feature selection, the model training time is substantially reduced from 153.83 seconds to 116.74 seconds, representing an approximate reduction of 24.11%. These findings indicate that the 12 critical features identified through the feature selection algorithm effectively capture the intrinsic data structure and informational content, thereby achieving substantial computational efficiency improvements while maintaining classification performance.

To validate the reliability and superiority of the proposed model, comparative evaluations are conducted against RF, iTransformer, iTransformer-BiLSTM (iTBM), and iTBU. To mitigate potential errors, the mean results from ten-fold cross-validation are adopted as the final evaluation metrics. The comparative results for the training set are illustrated in Figure 7, while the test set comparison outcomes are depicted in Figure 8.

As illustrated in Figure 7, the proposed methodology demonstrates superior performance across all evaluated metrics. Specifically, the accuracy attains 0.984, precision reaches 0.983, recall achieves 0.985, F1-score stands at 0.983, and specificity measures 0.989. Relative to the benchmark models, the proposed approach realizes improvements of up to 2.9% in accuracy, 2.9% in precision, 3.0% in recall, 2.7% in F1-score, and 1.5% in specificity.

Figure 8 reveals that the proposed methodology sustains optimal performance on the test set, wherein the accuracy attains 0.951, precision reaches 0.952, recall achieves 0.955, F1-score stands at 0.970, and specificity measures 0.986. Compared with the baseline models, enhancements of up to 3.9% in accuracy, 3.8% in precision, 4.0% in recall, 4.2% in F1-score, and 1.2% in specificity are achieved.

The proposed methodology exhibits minimal performance degradation between the training and test sets, with an average decline of approximately 2.19% across the five performance indicators, substantially lower than that observed for alternative comparative methods. This phenomenon substantiates the exceptional generalization capability of the proposed approach and its efficacy in mitigating overfitting phenomena. Notably, the performance decrements in F1-score and specificity are particularly modest (1.32% and 0.30%, respectively), thereby corroborating the methodology's capacity to maintain consistently high-level performance across disparate datasets.

(1) The proposed milling chatter detection methodology based on the optimized iTBU-RF ensemble model demonstrates efficacy in identifying three distinct machining states during the milling process: stable cutting, slight chatter, and severe chatter, thereby furnishing reliable technical support for chatter prevention in practical machining operations.

(2) Analysis employing Pearson correlation coefficients and significance levels reveals that the features STD, RMS, ABS, PTP, GFFS, WPE 3, MSFS, MAF, WPE 2, WPE 7, WPE 4, and WPE 5 exhibit robust correlations with machining states.

(3) Following hyperparameter optimization, the model detection performance demonstrates comprehensive enhancement. The accuracies of the training, validation, and test sets improve by 1.78%, 2.40%, and 1.78%, respectively. Concurrently, the optimized model exhibits more balanced performance across diverse datasets.

(4) Upon feature selection, the feature dimensionality is reduced from 20 to 12, achieving a reduction in training time from 153.83 seconds to 116.74 seconds (approximately 24.11% decrease) while maintaining classification performance.

(5) In comparison with baseline methodologies, the proposed model achieves optimal performance across all evaluation metrics. Notably, while maintaining high accuracy (95.1%), it attains superior F1-score (0.970) and specificity (0.986), demonstrating significant advantages in balanced identification of diverse chatter states.