Bayesian Model Updating for Chatter in Milling

1. Introduction

Regenerative chatter is a common and detrimental form of unstable self-excited vibration in machining. If not prevented, it can cause several adverse effects such as machine tool damage, poor surface finish, and accelerated tool wear, ultimately constraining the material removal rate and limiting overall machining productivity [1,2,3]. To suppress regenerative chatter in milling, stability analysis is employed to predict the boundary between stable and unstable cutting parameters. The results are typically presented in the form of a Stability Lobe Diagram (SLD), which determines regions of stability and instability as functions of spindle speed and axial depth of cut. These lobes are constructed by modeling the dynamic interaction between the machine tool structure and the cutting process, requiring the structural dynamics which are usually represented by the Frequency Response Function (FRF) or modal parameters measured at the tooltip or another highly flexible component of the setup [4,5].

The conventional approach to predicting modal parameters is through Experimental Modal Analysis (EMA), such as an impulse hammer test. However, a critical limitation arises because these tests are conducted while the machine is idle, with the spindle stationary and no cutting forces present. This idle state fails to capture the significant changes in system dynamics that occur under actual operational conditions. Consequently, an SLD based on an FRF from an idle-state EMA may be inaccurate in predicting chatter-free cutting parameters. To overcome this discrepancy and improve predictive accuracy, it is essential to calibrate or update the model parameters based on experimental observations collected while the machine is under operational conditions.

An effective approach to improving the chatter model’s accuracy is inverse chatter analysis, where the model parameters are tuned based on experimentally observed stability borders [6,7,8,9,10]. Machine Learning (ML) models have also been employed to construct SLD based on in-process sensor data [11,12,13,14,15]. By relying on operational measurements rather than offline modeling, both these approaches offer greater adaptability and accuracy in capturing the time-varying dynamics of machining processes. However, a key limitation is their reliance on extensive training datasets, which often require a significant amount of data collected under unstable cutting conditions, raising safety and productivity concerns in industrial applications. Alternatively, EMA can be performed using machining forces instead of an impulse hammer [16,17]. With this approach, the FRFs are measured during operation. However, since machining forces are harmonic at the spindle rotation frequency, specialized workpieces or experiments must be designed to broaden the excitation frequency spectrum.

Operational Modal Analysis (OMA) [18] has been adopted as a powerful technique for identifying system dynamics directly from vibration signals measured during the cutting process. Unlike modal testing, which characterizes the system’s structural dynamics, OMA captures the poles of the closed-loop machining system, reflecting the combined effects of the structural dynamics and the regenerative cutting process. Kim and Ahmadi [19] applied OMA to in-process vibration measurements during turning to extract the dominant system pole approaching instability. Building on this, a clustering-based algorithm was later introduced by Ebrahimi-Tirtashi and Ahmadi [20] to automate the OMA process for potential online applications. Zorlu et al. [21] further advanced the use of OMA by adapting it to milling operations, where the dynamics are Linear Time-Periodic (LTP), unlike the Linear Time-Invariant (LTI) dynamics in turning. To address this challenge, they adopted the lifting method, which enables mapping the LTP system into an equivalent LTI representation.

In this study, the OMA approach of Zorlu et al. [21] is coupled with Bayesian Model Updating (BMU) [22,23,24] to calibrate the modal parameters of the machining system based on the vibrations measured during standard, stable milling operations. Ahmadi introduced this approach for the LTI chatter dynamics in turning, which is often formulated accurately in the frequency domain [25]. However, its feasibility and practicality for the LTP chatter dynamics in milling are unexplored. Discrete-time domain methods such as the Semi-Discretization Method (SDM) [26,27] are widely used to link the structural modal parameters of the machine tool to the poles of the closed-loop time-periodic system, enabling deterministic forward prediction of stability. BMU solves the inverse problem of adjusting the modal parameters (natural frequency, damping ratio, and stiffness) based on the closed-loop poles extracted by OMA from measured vibration signals. This results in a more accurate representation of the system’s in-process dynamics and improves the reliability of chatter stability predictions. We show that, unlike in turning, directly applying BMU in milling can be computationally prohibitive because it requires repeatedly estimating the eigenvalues of the state transition matrix in discrete-time domain models. To address this issue, we augment the BMU workflow with a surrogate function that approximates those eigenvalues accurately and efficiently enough to be used in the model updating algorithm.

In the following sections, we begin with a concise overview of milling chatter dynamics formulated in the discrete-time domain. We then describe how OMA is applied using the lifting method and TDPR framework developed by Zorlu et al. [21]. The poles identified in this process serve as inputs to the BMU framework presented in Section 2.2. To validate the effectiveness of the proposed method in predicting the stability boundary, both numerical simulations and an experimental case study are presented in Section 3.

2. Milling Dynamics

To model the dynamic interaction between the tool and workpiece during milling, a simplified representation of the structural system is first considered. Single-degree-of-freedom (SDOF) model is used here to simplify the explanation of the proposed approach. However, the method is general and can be extended to multi-degree-of-freedom (MDOF) systems to capture more complex structural dynamics.

Figure 1 illustrates an SDOF dynamic model representing the relative vibrations between the tool and the workpiece in milling. The structure is assumed to be flexible in the feed (x) direction, and the tool rotates around its axis at an angular speed of $\mathsf{\Omega }$. The regenerative effect introduces a delay $\tau ={\displaystyle \frac{2\pi }{\mathsf{\Omega }N}}$ in the system response, resulting in the following Delay Differential Equation (DDE) [28]:

$$\ddot{x}\left(t\right)+2{\zeta }_x{\omega }_n\dot{x}\left(t\right)+\left[{\omega }_n^2-{\displaystyle \frac{aH\left(t\right){\omega }_n^2}{k_x}}\right]x\left(t\right)={\displaystyle \frac{-aH\left(t\right){\omega }_n^2}{k_x}}x(t-\tau )$$

(1)

where ${\zeta }_x$, ${\omega }_n$, and $k_x$ are damping ratio, natural frequency, and stiffness, respectively. a is the nominal depth of cut and $H\left(t\right)$ is cutting force directional coefficient, which is periodic with respect to the tooth passing period $\tau$, i.e., $H\left(t\right)=H(t+\tau )$ and is defined as:

$$\begin{array}{cc}\hfill H\left(t\right)& =\sum _{j=1}^N{g}_j\left(t\right)\left[-K_{t}cos{\theta }_j\left(t\right)-K_{n}sin{\theta }_j\left(t\right)\right]sin{\theta }_j\left(t\right)\hfill \\ \hfill g_j\left(t\right)& ={\displaystyle \frac{1}{2}}[sgn({\theta }_j-{\theta }_{st})-sgn({\theta }_j-{\theta }_{ex})]\hfill \end{array}$$

(2)

where $K_t$ and $K_n$ are tangential and normal cutting coefficients, $g_j\left(t\right)$ is a sign function that determines if the tooth is in or out of cut. N is the number of teeth and ${\theta }_{st}$ and ${\theta }_{ex}$ are start and exit angles of a tooth to and from the cut, respectively.

One of the ways to determine the stability of vibrations governed by Eq. 1 is using SDM, in which the delay interval $\tau$ is divided into equal subintervals of length $\mathsf{\Delta }t=\tau /r$. Within each subinterval, the delayed state $x(t-\tau )$ is assumed to be constant. This approximation transforms the original infinite-dimensional DDE into a finite-dimensional discrete-time system:

$$\begin{array}{cc}\hfill {\mathbf{z}}_{i+1}& ={\mathbf{D}}_i{\mathbf{z}}_i\hfill \\ \hfill {\mathbf{z}}_i& ={[x_i,{\dot{x}}_i,x_{i-1},\dots ,x_{i-r}]}^{\top }\hfill \\ \hfill {\mathbf{D}}_i& =\left[\begin{array}{ccccc}exp\left({\mathbf{A}}_i\mathsf{\Delta }t\right)& 0& \dots & 0& \left(exp\left({\mathbf{A}}_i\mathsf{\Delta }t\right)-I\right){\mathbf{A}}_i^{-1}{\mathbf{B}}_i\\ I& 0& \dots & 0& 0\\ 0& I& \dots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \dots & I& 0\end{array}\right]\hfill \end{array}$$

(3)

where ${\mathbf{z}}_i$ is an augmented state vector that includes the state $x\left(t\right)$ and its delayed values at each subinterval, and $D_i$ is the corresponding state transition matrix constructed from system matrices ${\mathbf{A}}_i$, ${\mathbf{B}}_i$ as:

$${\mathbf{A}}_i=\left[\begin{array}{cc}0& 1\\ {\omega }_n^2(1-aH\left(t\right)/k_x)& -2{\zeta }_x{\omega }_n\end{array}\right];{\mathbf{B}}_i=\left[\begin{array}{c}0\\ -aH\left(t\right){\omega }_n^2/k_x\end{array}\right]$$

(4)

By recursively applying the transition matrices over one full delay period, the system’s monodromy matrix $\mathsf{\Phi }$ is obtained:

$$\mathsf{\Phi }={\mathbf{D}}_{r-1}{\mathbf{D}}_{r-2}\dots {\mathbf{D}}_0$$

(5)

The eigenvalues of matrix $\mathsf{\Phi }$ correspond to the Floquet multipliers of the delayed periodic system described by Eq. 1. Based on Floquet’s theorem for delayed LTP systems, the periodic motion loses stability when any eigenvalue of the monodromy matrix has a modulus greater than one.

Experimentally, the dominant multiplier of the milling dynamics can be obtained using the OMA approach presented by Zorlu et al. [21]. A summary of this approach is presented in the following subsection, and the resulting multipliers are used as observations in the proposed BMU framework to improve the chatter model.

2.1. Operational Modal Analysis

The lifting technique [29,30] is a method to transform LTP systems into an equivalent LTI system by sampling the system’s response at fixed intervals synchronized with its periodicity. This results in a “lifted” system that evolves in discrete time, with each step representing one full period of the original system.

Considering the tool’s vibration along the x-axis, if r samples are taken during each tooth-passing period $\tau$, the corresponding sampling time is given by $\mathsf{\Delta }t=\tau /r$. Assuming that the observation window spans $N_c$ full tooth-passing intervals, the vibration measurements collected during the $c^{th}$ cycle ($c=1,\dots ,N_c$) can be organized as follows:

$${\mathbf{M}}_c={\left[x_{c,1},x_{c,2},\dots ,x_{c,r}\right]}^{\top }$$

(6)

Here, $x_{c,j}$ refers to the $j^{th}$ sampled data point in the $c^{th}$ cycle, while $x_{c,r}$ denotes the final measurement in that cycle.

Although the measurement vector in Eq. 6 is formed using data from one physical accelerometer, each row of this vector can be interpreted as the response of a virtual sensor that captures the tool tip motion at a consistent angular position within each tooth-passing cycle. Consequently, when this vector serves as the output of the lifted LTI model, the system’s periodic characteristics are effectively removed from the output. Thus the the lifted response matrix, $\mathbf{X}$ can be organized as below, where each row shows the signal sampled once at period $\tau$. For more details, refer to [21].

$$\mathbf{X}=\left[{\mathbf{M}}_1\dots {\mathbf{M}}_{N_c}\right]={\left[\begin{array}{cccc}{x}_{1,1}& x_{1,2}& \dots & x_{1,r}\\ ⋮& ⋮& \ddots & ⋮\\ x_{N_c,1}& x_{N_c,2}& \dots & x_{N_c,r}\end{array}\right]}^{\top }$$

(7)

Various OMA methods can be employed to identify the dominant pole of the lifted time-independent dynamics, which also represents the dominant Floquet multiplier of milling dynamics. The Time Domain Polyreference (TDPR) method is an OMA technique used to identify the system poles of a vibratory system from its response under ambient or process-induced excitation [18]. In the context of machining, TDPR is particularly valuable for estimating the dominant poles of the closed-loop system dynamics during cutting operations without requiring external excitation sources [19,21].

TDPR constructs an autoregressive (AR) model of the system’s correlation functions (CF) derived from the measured response signals.

$${\mathbf{R}}_X\left(c\right)={\displaystyle \frac{1}{N_c-c}}\sum _{k=1}^{N_c-c}{\mathbf{X}}_k{\mathbf{X}}_{k+c}^{\top },c=0,1,\dots ,N_c-1$$

(8)

where ${\mathbf{X}}_k$ denotes the $k^{\mathrm{th}}$ column of the matrix $\mathbf{X}$, representing the output of the $k^{\mathrm{th}}$ lifted vibration signal. The companion matrix form of the AR model can be used to obtain the closed-loop pole of the system:

$$\begin{array}{cc}\hfill {\mathbf{U}}_d(c+1)& =\mathbf{G}{\mathbf{U}}_d\left(c\right)\hfill \\ \hfill {\mathbf{U}}_d\left(c\right)& ={\left[\begin{array}{ccc}{\mathbf{R}}_X{(c-n_a+1)}^{\top }& \dots & {\mathbf{R}}_X{\left(c\right)}^{\top }\end{array}\right]}^{\top }\hfill \end{array}$$

(9)

The eigenvalues of $\mathbf{G}$ represent the closed-loop poles of the system and quantify the stability of the process experimentally.

To distinguish physical poles from numerical artifacts introduced by overfitting, a stabilization diagram is constructed by incrementally increasing the model order. A pole is considered stabilized, and thus likely to represent a true system mode, if its frequency and damping ratio remain consistent across successive orders of the AR model. Given a discrete-time pole $\mu$ as an eigenvalue of $\mathbf{G}$, the corresponding continuous-time pole is obtained as $\lambda ={\displaystyle \frac{ln\left(\mu \right)}{\mathsf{\Delta }t}}$. The corresponding estimated natural frequency $\widehat{\omega }$ and the damping ratio $\widehat{\zeta }$ are then computed as:

$$\widehat{\omega }=\left|\lambda \right|,\widehat{\zeta }=-{\displaystyle \frac{Re\left(\lambda \right)}{\left|\lambda \right|}}$$

(10)

The dominant pole identified in this manner is used as a proxy for the system’s proximity to regenerative instability and can subsequently serve as the basis for the likelihood function in the BMU framework.

2.2. Model Updating

Let $\mathbf{D}$ represent the observations, i.e., dominant pole frequencies and damping ratios identified via OMA during the milling process. For a given model class $\mathcal{M}$, the objective is to use these data to update the relative probability of each candidate model characterized by the parameter vector through the application of Bayes’ Theorem [23,24]:

$$p\left(\mathit{\theta }\right|\mathbf{D},\mathcal{M})={\displaystyle \frac{p\left(\mathbf{D}\right|\mathit{\theta },\mathcal{M}\left)p\right(\mathit{\theta }\left|\mathcal{M}\right)}{p\left(\mathbf{D}\right|\mathcal{M})}}$$

(11)

where $p\left(\theta \right|\mathbf{D},\mathcal{M})$ is the posterior PDF given the observation and model class $\mathcal{M}$, $p\left(\theta \right|\mathcal{M})$ is the prior PDF which expresses the initial belief about the system and quantifies the prior plausibility of each model in the model class $\mathcal{M}$, $p(\mathbf{D}\mid \mathcal{M})$ is the evidence or marginal likelihood for the model class $\mathcal{M}$ provided by data $\mathbf{D}$ and serves as the normalizing constant in Bayes’ theorem to ensure the posterior distribution integrates to one. $p\left(\mathbf{D}\right|\theta ,\mathcal{M})$ is the likelihood function which quantifies the agreement between the experimental observations (i.e., the dominant closed-loop pole’s damping ratio and frequency identified by OMA ) and their theoretical counterpart predicted by the SDM for a given set of parameter variations $\theta$. Specifically, let ${\widehat{\omega }}_{p,j}$ and ${\widehat{\zeta }}_{p,j}$ denote the frequency and damping ratio of the dominant pole identified from the $j^{\mathrm{th}}$ OMA repetition at the $p^{\mathrm{th}}$ combination of cutting depth and spindle speed.

The likelihood is then given by:

$$p(\mathbf{D}\mid \mathit{\theta },\mathcal{M})=\prod _{j=1}^{N_s}\prod _{p=1}^{N_p}p({\widehat{\omega }}_{p,j}\mid \mathit{\theta },\mathcal{M})p({\widehat{\zeta }}_{p,j}\mid \mathit{\theta },\mathcal{M})$$

(12)

where $N_p$ is the number of tested cutting pairs, and $N_s$ is the number of OMA repetitions per pair.

The connection between the measured dominant pole through OMA and the predicted one is given in terms of the prediction error as:

$${\widehat{\omega }}_{p,j}={\omega }_p+e_{{\omega }_{p,j}},{\widehat{\zeta }}_{p,j}={\zeta }_p+e_{{\zeta }_{p,j}}$$

(13)

and the prediction errors are assumed to follow independent zero-mean Gaussian distributions with variances ${\sigma }_{\omega }^2$ and ${\sigma }_{\zeta }^2$, respectively. Thus, the likelihood function can be written as:

$$\begin{array}{cc}\hfill p\left(\mathbf{D}\right|\mathit{\theta },\mathcal{M})& ={\displaystyle \frac{1}{2\pi \sqrt{{\sigma }_{\omega }{\sigma }_{\zeta }}}}exp\left({\displaystyle \frac{-1}{2}}\sum _{p=1}^{N_p}{J}_p\left(\mathit{\theta }\right)\right),\hfill \\ \hfill J_p\left(\mathit{\theta }\right)& =\sum _{j=1}^{N_s}\left({\displaystyle \frac{{({\widehat{\omega }}_{p,j}-{\omega }_p\left(\mathit{\theta }\right))}^2}{{\sigma }_{\omega }^2}}+{\displaystyle \frac{{({\widehat{\zeta }}_{p,j}-{\zeta }_p\left(\mathit{\theta }\right))}^2}{{\sigma }_{\zeta }^2}}\right)\hfill \end{array}$$

(14)

When Bayesian updating is applied with the parameter vector $\theta$ including the prediction error variances for frequency and damping ratio ${\sigma }_{\omega }^2$ and ${\sigma }_{\zeta }^2$, it is important to select their prior distributions carefully. The chosen priors must ensure that the relative influence of frequency and damping components in the likelihood function is properly balanced, so that neither dominates the posterior update and both contribute meaningfully to the identification of modal parameters. To manage this balance, following the evidence-based approach proposed by Goller et. al. [24], multiple model classes $\mathcal{M}\left({\alpha }_k\right),k=1..N_k$ are defined, in which the ratio of prediction-error variances for damping and frequency is fixed at ${\alpha }_k={\displaystyle \frac{{\sigma }_{\zeta }}{{\sigma }_{\omega }}}$. The likelihood function for the model class ${\mathcal{M}}_k$, where ${\mathcal{M}}_k$ is shorthand for $\mathcal{M}\left({\alpha }_k\right)$, is obtained by substituting ${\alpha }_k$ into Eq. 14 as:

$$J_p\left(\mathit{\theta }\right)=\sum _{j=1}^{N_s}\left({\displaystyle \frac{{({\widehat{\omega }}_{p,j}-{\omega }_p\left(\mathit{\theta }\right))}^2}{{\sigma }_{\omega }^2}}+{\displaystyle \frac{{({\widehat{\zeta }}_{p,j}-{\zeta }_p\left(\mathit{\theta }\right))}^2}{{\left({\alpha }_k{\sigma }_{\omega }\right)}^2}}\right)$$

(15)

Here, $\mathit{\theta }=[{\theta }_{\omega },{\theta }_{\zeta },{\theta }_k,{\sigma }_{\omega }]$, where the random variables ${\theta }_{\omega }$, ${\theta }_{\zeta }$, and ${\theta }_k$ represent the relative variations of the system’s true modal parameters from their respective nominal values, defined as:

$$\begin{array}{c}\hfill {\omega }_n={\overline{\omega }}_n(1+{\theta }_{\omega })\\ \hfill {\zeta }_x={\overline{\zeta }}_x(1+{\theta }_{\zeta })\\ \hfill k_x={\overline{k}}_x(1+{\theta }_k)\end{array}$$

(16)

Beyond updating the structural modal parameters within each individual model class ${\mathcal{M}}_k$, a second layer of Bayesian inference is employed to evaluate and compare multiple model classes corresponding to different values of the weighting parameter $\alpha$. This approach allows us to determine which model class is more likely to describe the observed data. In doing so, it identifies the most appropriate value of $\alpha$, which defines the optimal balance between the contributions of closed-loop frequency and damping ratio prediction errors in the likelihood function. The appropriate Bayes’ theorem at the model class level is [24]:

$$p({\mathcal{M}}_k\mid \mathbf{D},\mathcal{M})={\displaystyle \frac{p\left(\mathbf{D}\right|{\mathcal{M}}_k)p\left({\mathcal{M}}_k\right|\mathcal{M})}{p\left(\mathbf{D}\right|\mathcal{M})}}$$

(17)

where $\mathcal{M}$ is a set of $N_k$ model classes ${\mathcal{M}}_k$.

Calculating the model evidence (marginal likelihood) analytically in Bayesian inference is often impractical, due to the high-dimensional integration required over the parameter space. To address this challenge, sampling-based methods are commonly employed to approximate the posterior distribution without explicitly computing the evidence [31,32]. Among these, Transitional Markov Chain Monte Carlo (TMCMC) [33] is particularly suitable for our application, as it not only facilitates efficient sampling through a sequence of intermediate distributions bridging the prior and posterior but also provides a natural and accurate estimate of the model evidence. This capability makes TMCMC especially advantageous when performing Bayesian model selection, where comparing the marginal likelihoods of candidate models $p\left(\mathbf{D}\right|{\mathcal{M}}_k)$ is essential for determining the most plausible explanation of observed data $p\left(\mathbf{D}\right|\theta ,\mathcal{M})$. [34]

To evaluate the likelihood function in this study, we need to predict the dominant closed-loop pole, or in other words, its frequency and damping ratio, for each candidate set of modal parameters. This is achieved by building the discrete-time transition matrix of the system and calculating its eigenvalues. However, doing this many times during the TMCMC sampling becomes computationally prohibitive due to the size and complexity of the matrices involved. To reduce the computational cost, a surrogate model is trained to estimate the relationship between the modal parameters and the dominant pole characteristics. In this study, we use a Gaussian Process Regression (GPR) model [35] with a Matérn kernel to predict the closed-loop frequency and damping ratio based on the modal parameters. The surrogate model is trained using a set of parameter samples selected from a uniform grid covering the range of the prior distribution. For each sample, the full-order model is run once to calculate the dominant eigenvalue of the system’s transition matrix. This eigenvalue is then converted into its corresponding continuous-time frequency and damping ratio, which are used as the target outputs for training the GPR model. Once the surrogate model is trained, it can quickly estimate the system’s dynamic response for any set of modal parameters. This makes it possible to evaluate the likelihood function much faster than computing eigenvalues directly, greatly reducing the overall cost of the TMCMC-based Bayesian updating process.

3. Results

3.1. Numerical Simulation

To verify the effectiveness of the Bayesian Model Updating (BMU) approach in calibrating the modal parameters and quantifying the uncertainties in the milling process, a numerical simulation was performed on a single degree of freedom (SDOF) system.

Figure 2 shows the block diagram used to simulate the single degree-of-freedom (SDOF) milling dynamics numerically, considering $K_t=650$ MPa and $K_n=110$ MPa for tangential and normal cutting force coefficients, and the dynamics of tool-workpiece assembly are considered as ${\omega }_n=2\pi \times 194$ rad/s, ${\zeta }_x=0.005$, and $k_x=2.35$$N/\mu m$ for Natural frequency, damping ratio, and stiffness, respectively.

The stability lobe diagrams of the setup, obtained by the SDM, are shown in Figure 3. The input force consists of periodic machining forces superimposed by broadband white noise to account for the random forces during the process. Simulations were conducted in MATLAB Simulink 2023b. Process vibrations were simulated for milling with 2000 and 2100 rev/min spindle speed and various stable DOC values marked with circles on the SLD in Figure 3.

The power spectral density (PSD) of simulated acceleration at 2100 rev/min and $0.4$ mm depth is displayed in Figure 4. The PSD of the simulated vibrations shows peaks at the tooth-passing frequency harmonics and a non-harmonic peak at around $195.8$ Hz.

The measured signals are lifted at the tooth passing period and organized in the measurement matrix, $\mathbf{X}$, described in Eq. 7. Subsequently, treating each row of the $\mathbf{X}$ matrix as a separate lifted sensor, we computed the PSD matrix of the lifted sensors; the top singular value of the resulting PSD is shown in Figure 5. The Nyquist frequency of the signals after lifting reduces to half of the tooth-passing frequency, $f_{tp}/2=35$ Hz.

Since milling dynamics is periodic, according to Floquet theory, the spectrum exhibits frequency components at $f_c\pm k\times f_{tp}$ and $-f_c\pm k\times f_{tp}$, where k is an integer, $f_c$ represents the oncoming chatter frequency, and $f_{tp}$ denotes the toothpassing frequency [36]. Considering that the response before lifting shows a chatter frequency of $195.8$ Hz, after lifting, this frequency is aliased to one of its harmonics at around $14.2$ Hz, which corresponds to $-f_c+3\times f_{tp}$.

The stabilization diagram resulting from TDPR identification of the lifted response at point P with various model orders $n_a$ is shown in Figure 5. While the diagram includes a mix of physical and spurious poles, the column of circles with similar shading indicates a physical pole at the corresponding frequency.

The dominant pole was identified at $N_p=7$ combinations of spindle speed and cutting depth, and identification was repeated $N_s=10$ times at each combination. The mean value and standard deviation of the identified poles are shown in Table 1.

To reduce the high computational cost associated with repeatedly evaluating the likelihood function in TMCMC, particularly the eigenvalue analysis of the state transition matrix $\mathbf{G}$, a GPR surrogate model with a Matérn kernel was employed. This surrogate was trained to approximate the dominant closed-loop frequency and damping ratio as functions of the model parameters $[{\theta }_{\omega },{\theta }_{\zeta },{\theta }_k]$, which represent expected variations in the natural frequency, damping ratio, and mode shape scaling, respectively. These parameter ranges were uniformly sampled to reflect the prior distribution based on initial uncertainty. A training dataset consisting of 1000 combinations of model parameters was used to construct the GPR model. Once trained, the surrogate provides rapid predictions of the closed-loop poles for any arbitrary parameter combination required by the TMCMC algorithm, thereby enabling efficient and accurate sampling of the posterior distribution without recomputing expensive matrix operations at every iteration.

The prior distributions of the model parameters, $\mathit{\theta }$, are assumed to be uniform over the range $[-0.3,0.3]$, and the prior distribution of the prediction error variance, ${\sigma }_{\omega }^2$, is taken as uniform over the interval $[0,1]$.

The TMCMC algorithm was employed to generate 5000 samples from the posterior distribution of the model parameters across model classes with ${\alpha }_k$ over the range of $0.001$ to 1. The model evidence (log marginal likelihood) for each class was also computed using TMCMC and is presented in Figure 6. The evidence reaches its maximum at $\alpha =0.01$, indicating that this model class has the highest plausibility among the candidates. Consequently, the model class corresponding to this amount is selected for parameter updating.

For the selected model class, samples were generated from the posterior distribution of the three model parameters. The diagonal plots in Figure 7 show the histograms of each parameter, while the off-diagonal plots display scatter views of the samples from the prior (gray points) and the posterior (green points).

Since the analysis is conducted entirely through numerical simulation, the dominant closed-loop poles computed mathematically from the discretized state-space model and those identified through OMA yield identical results. This equivalence is evident in the plotted figure, where all estimated parameter variations converge to zero. The convergence indicates that the posterior distributions of the model parameters remain centered at their nominal values, confirming that the closed-loop poles derived from the numerical discretization perfectly match those extracted via OMA. This outcome confirms the consistency of the surrogate-based Bayesian updating framework when applied to simulated data, where no modeling error or measurement noise is present.

Figure 8 presents the histogram of the frequency prediction error, illustrating the distribution of discrepancies between the predicted and observed dominant frequencies.

3.2. Experimental Case Study

The experimental tests presented in this section are based on the milling setup and data reported in Zorlu et al. [21]. A brief description of the experimental configuration is provided below, and the closed-loop poles identified via OMA in that study are used as inputs to the Bayesian inference framework developed in this work.

A set of milling tests is conducted using a single-degree-of-freedom milling setup shown in Figure 9a. The setup consists of a 6061-T6 Aluminum workpiece mounted on a uniaxial compliant fixture, machined by a two-fluted carbide cutting tool with a diameter of $12.7$ mm. As shown in Figure 9a, the setup is flexible in the x (feed) direction, and the workpiece is machined using a feedrate of 0.2 mm/tooth in a slotting (full immersion) operation. The tangential and normal linear mechanistic cutting force coefficients were identified by the orthogonal to oblique transformation as $K_t=650$ MPa and $K_n=110$ MPa, respectively [37]. The piezoelectric accelerometer (PCB 352C22) is mounted on the flexure in the x-direction. A Hall-effect sensor is used to detect spindle rotation instants and provides the trigger signal for synchronous sampling in dSpace MicroLabBox. Figure 9b illustrates the schematic of the data acquisition setup. This signal provides a precise reference for identifying spindle rotation events and synchronizing the vibration signal with each revolution. Therefore, the accelerometer’s analog signal is sampled at the same tool rotational angles in each period. The number of data points after each trigger is determined by the tooth-passing period and the sampling period. For instance, at 2100 rev/min, one full tooth-passing period contains 285 data points when sampled at a rate of 20000 samples per second.

Before removing any material from the workpiece, the modal parameters of the setup in the feed direction are determined by impulse hammer tests as ${\omega }_n=2\pi \times 195.3$ rad/s, $k_x=2.34$$N/\mu m$, and ${\zeta }_x=0.004$. The stability lobe diagram obtained using the SDM is shown in Figure 10. Also shown in the figure are the experimentally determined stable and unstable points. The proposed OMA analysis was performed for seven of the tested stable points, which are marked as P1 to P7 in the figure.

Four sets of full-immersion milling experiments were conducted at spindle speeds between 2000 to 2400 rev/min. Figure 11 presents the OMA stabilization diagrams for two representative measurements at 2200 rev/min with $0.6$ mm depth (point P4) and $0.7$ mm depth (point P5), identified using the TDPR method. More results are reported in ref. [21]. Figure 11a and Figure 11c are the PSD of the measured machining sound at point P4 and P5, and Figure 11b and Figure 11d are the corresponding PSD of the lifted vibration signals. In both cases, the vibrations are stable, and the sound spectrum is dominated by harmonics of the tooth-passing frequency. Additionally, the power spectral density (PSD) of the lifted vibration signals exhibits a distinct peak near 20 Hz, which the TDPR method identifies as a stable pole. A similar analysis was performed for all test conditions ($N_p=7$) indicated in Figure 10, and the extracted dominant poles are provided in Table 2.

Since each test was repeated multiple times ($N_s=10$) to ensure consistent results, a substantial amount of material was removed from the workpiece after completing the tests at each cutting pair. This material removal led to a notable change in the system’s modal parameters. As a result, the dominant poles identified by OMA at 2200 and 2400 rev/min deviated significantly from the chatter stability predictions based on the SLD constructed using modal parameters obtained from the initial impact hammer test. This discrepancy indicates that the system’s modal parameters evolved during the course of testing, and therefore, they need to be updated to reflect the in-process measurements.

The prior distribution of the modal parameters variations ${\mathit{\theta }}_m=[{\theta }_{\omega },{\theta }_{\zeta },{\theta }_k]$ are assumed to be uniform across $[-0.1,0.1]$ for ${\theta }_{\omega }$ and ${\theta }_k$, and $[-0.5-0.5]$ for ${\theta }_{\zeta }$. The prior distribution of the prediction error variance ${\sigma }_{\omega }^2$ is assumed to be uniform across 0 and 1. A training dataset consisting of 3375 different combinations of model parameters was used to construct the GPR model, enabling rapid estimation of the eigenvalues of the monodromy matrix for any given parameter set required by the TMCMC algorithm. The TMCMC generated 5000 samples from the posterior distribution of the parameters based on the assumed priors and the likelihood function. A set of model class with ${\alpha }_k$ in the range of $0.001$ and 1 are considered and the log marginal likelihood for each $\alpha$ is shown in Figure 12. The model class with $\alpha =0.006$ is selected accordingly, and the sampling process is repeated for 10000 samples according to the selected $\alpha$.

For the selected model class, samples were generated from the posterior distribution of the three model parameters. Figure 13 depicts the joint and marginal distributions of the model parameters ${\theta }_{\omega }$, ${\theta }_{\zeta }$, and ${\theta }_k$ used in the Bayesian model updating process. The diagonal plots show the marginal posterior distributions of each parameter, where the histograms indicate high concentration and reduced uncertainty relative to the prior. The off-diagonal plots display the pairwise joint distributions: the gray background points represent the prior samples drawn from uniform distributions, while the green points illustrate the updated posterior samples. In all off-diagonal plots, the posterior samples cluster tightly, indicating strong parameter identifiability and convergence toward the updated values below. The histogram of the frequency prediction error and its joint distribution with ${\theta }_{\omega }$ are shown in Figure 14, further confirming the convergence of the TMCMC algorithm.

$$\mathbb{E}[{\theta }_{\omega },{\theta }_{\sigma },{\theta }_k]=[3.6\times {10}^{-4},0.35,0.003]$$

(18)

Figure 15 shows samples of the posterior distributions for the same example, except that the prior is normal distribution for $\mathit{\theta }$ with mean value $[0,0.35,0]$ and $[0.005,0.01,0.01]$ for variances and uniform distribution for prediction error variance, ${\sigma }_{\omega }^2$ across 0 and 1.

Using the mean values of the updated parameters from Eq. 18, the revised SLD is presented in Figure 16. This figure also includes the experimentally determined stable (circles) and unstable (crosses) combinations of spindle speed and axial depth of cut, all of which show good agreement with the stability boundaries predicted by the updated model. Notably, the updated SLD is consistent with the stability boundary obtained from the impact hammer test conducted after machining, as reported by Zorlu et al. [21]. Figure 17 compares the frequency response function computed from the updated parameters with that obtained through post-process modal testing.

4. Conclusions

This study applies a Bayesian inference framework to refine the modal parameters of a milling system using in-process vibration data. It addresses a well-known limitation of conventional modal testing, which is typically conducted on idle machines and therefore fails to capture the true dynamics present under operational conditions. This discrepancy often results in inaccurate chatter stability predictions. To overcome this issue, OMA was employed to extract the closed-loop poles of the system from vibration signals measured during stable milling. The experimentally identified poles were then used to update the prior distributions of the system’s modal parameters. An experimental case study demonstrated that the updated parameters improve the accuracy of the stability lobe diagram in comparison with those generated from offline modal testing data.

An important advantage of this approach is its potential to eliminate the need for impact hammer testing. The experimental case study showed that, even for a uniform prior distribution, which reflects minimal initial knowledge of the system’s dynamics, the Bayesian updating approach can still converge to the true modal parameters using only stable operational data. This makes the method particularly attractive for model calibration in industrial environments, where minimizing manual intervention and reducing experimental effort are critical.

Despite its demonstrated effectiveness, the approach has certain limitations. First, it is not currently applicable for online monitoring or real-time parameter updating, as the process relies on computationally intensive inference. Second, when extended to multi-mode systems, the method results in a high-dimensional parameter space. This is because each additional mode introduces multiple parameters, which significantly increase the computational burden and the complexity of convergence. These factors limit the scalability of the approach for systems with multiple dominant modes.