Bayesian stability and force modeling for uncertain machining processes

ABSTRACT Accurately simulating machining operations requires knowledge of the cutting force model and system frequency response. However, this data is collected using specialized instruments

in an ex-situ manner. Bayesian statistical methods instead learn the system parameters using cutting test data, but to date, these approaches have only considered milling stability. This

paper presents a physics-based Bayesian framework which incorporates both spindle power and milling stability. Initial probabilistic descriptions of the system parameters are propagated

through a set of physics functions to form probabilistic predictions about the milling process. The system parameters are then updated using automatically selected cutting tests to reduce

parameter uncertainty and identify more productive cutting conditions, where spindle power measurements are used to learn the cutting force model. The framework is demonstrated through both

numerical and experimental case studies. Results show that the approach accurately identifies both the system natural frequency and cutting force model. SIMILAR CONTENT BEING VIEWED BY

OTHERS ENERGY EFFICIENCY IDENTIFICATION AND SURFACE ROUGHNESS PREDICTION USING CUTTING FORCE SIGNAL FOR COMPUTER NUMERICAL CONTROLLED MACHINE SYSTEMS Article Open access 16 August 2024

NEWTON–SIMPSON-BASED PREDICTOR–CORRECTOR METHODS FOR MILLING CHATTER STABILITY PREDICTION Article Open access 02 January 2025 MONITORING OF MANUFACTURING PROCESS USING BAYESIAN EWMA CONTROL

CHART UNDER RANKED BASED SAMPLING DESIGNS Article Open access 25 October 2023 INTRODUCTION Milling is a common machining process used to produce high-precision components1. Producing parts

for maximized profit requires selecting milling parameters, such as depth of cut and spindle speed, that minimize cost. However, there are physical constraints on this optimal parameter

selection. Some sets of cutting parameters may cause chatter which damages the part and tool2,3. Others may be stable, but not provide the required part geometric accuracy or surface

finish4,5. Others could exceed the machine tool’s capabilities, such as the maximum spindle power output6. Therefore, selecting optimal milling parameters requires the ability to predict

whether a cut will be stable and produce acceptable parts under realizable conditions. To enable pre-process selection of milling parameters, milling simulation has been widely studied in

the literature. For example, chatter can be predicted by calculating a stability map based on the cutting force model and system dynamics7,8,9, while surface location error predictions

estimate part errors due to the relative vibration of the tool and workpiece10,11. The simulation performance predictions can be used to optimize the milling process by selecting the most

productive, realizable parameters. However, these predictions are subject to both epistemic and aleatoric uncertainty. Epistemic uncertainty is attributed to a lack of information about the

system and the model inputs, such as the tool tip frequency response function (FRF) and cutting force model required to calculate the stability map12. Measuring these parameters requires

specialized instruments, such as cutting force dynamometers and impact hammers/accelerometers, which are not widely available in the industry. Therefore, in practice, these parameters may

not be measured. Furthermore, even if the measurements are performed, there is always uncertainty associated with the measurement results. For example, accelerometer mass loading13 and

changes in spindle dynamics with spindle speed14 can introduce uncertainty into FRF measurements. Cutting force measurements are subject to dynamometer dynamics and can vary with milling

parameters or as the tool wears15. While epistemic uncertainty can be reduced by measurements, aleatoric uncertainty is unattributable stochastic uncertainty which must be modeled as

essentially random. This is traditionally assumed to be unpredictable variations that may change between tests, such electrical noise in a sensor signal or cutting force variations due to

tool wear12. However, it can also be used to describe errors in the machining simulations themselves. Stability algorithms are uncertain due to assumptions or factors not included in the

model, which are separate from uncertainty in the model inputs. As an example, the zero-order approximation (ZOA) frequency domain method developed by Altintas and Budak is widely used

because of its computational efficiency, but it is less accurate for low radial immersion cuts16. Alternate stability algorithms may be more accurate, but they are generally more

computationally expensive2,3,17,18. Due to these difficulties, data-driven and machine learning methods for identifying the stability map have become a major research focus. These approaches

extrapolate the stability map from a small set of cutting tests. For example, Kawai et al. used Gaussian processes19, Greis et al. compared support vector machines, k-nearest neighbors, and

neural networks20, Postel et al. proposed both deep neural networks pre-trained by physics models21 and a genetic algorithm approach for inverse parameter identification22, and Deng et al.

trained their model using limiting axial depths rather than binary stable/unstable test results to reduce the required number of data points23. However, these methods have some limitations.

First, they usually require many cutting tests to learn the stability map and, more critically, may not identify tests for optimal learning. As a result, these approaches are best applied to

correct errors in a calculated stability map, rather than as a replacement for system measurements. Second, while these approaches can learn the stability map for one specific case, they

are not typically generalizable. For example, while the approaches can predict stability under the same conditions that were used for training, they may not extend to a different work

material or radial depth of cut. Third, these methods have focused on stability and have not considered other inputs, including the available spindle power, required surface finish, and

dimensional tolerances. Bayesian modeling is a probabilistic approach to modeling and updating problems. The approach assigns a probabilistic prediction about the system parameters and then

updates it based on a series of observations. Authors have applied this approach to model cutting forces24, surface roughness25, and surface location error26,27. Bayesian stability models

assign probabilistic predictions about how likely a given set of cutting parameters is to be stable, rather than a binary stable/unstable classification predicted by a stability algorithm.

This approach was first proposed by Karandikar et al. with a non-physics-informed Bayesian stability model. They addressed the first issue by allowing the user to select cutting tests which

provide new information about the system28. However, they did not connect the final stability predictions to system parameters. Physics-based Bayesian machining modeling is a more recent

area of study that directly connects the stability modeling uncertainty to the uncertainty in the machining system parameters. These approaches establish an initial probabilistic input

(called the prior) for the system dynamics and cutting force model and then update those probabilistic inputs based on a series of cutting tests to reduce the system uncertainty. Li et al.

proposed the first method in 2020, updating a measured FRF to account for shifts in spindle dynamics at higher spindle speeds29. Cornelius et al. demonstrated that Bayesian stability

modeling can be used to select productive machining parameters and showed how assumptions in the ZOA stability algorithm can cause errors in the predicted stability map30. Chen et al. used a

surrogate machine learning model to reduce computational time by approximating the output of the stability algorithm31. Schmitz et al. established the prior for the dynamics using

receptance coupling substructure analysis (RCSA) and incorporated the measured chatter frequency into the Bayesian updating to learn more from each test cut32. Ahmadi et al. applied Bayesian

learning to turning and used the measured tool vibrations to help learn the system dynamics from stable cuts33. Akbari et al. extended the work in ref. 32 by demonstrating a new test

selection approach focused on reducing uncertainty in the stability map and showing how the measured parameters can be used to model cutting test results across multiple machine tools34.

Cornelius et al. scaled the Bayesian approach to utilize high-performance computing resources and learn process damping35. Compared with the previous data-driven approaches, physics-guided

Bayesian methods address two challenges. First, the current knowledge about the system can be used to select informative test cuts, either to improve productivity30,32 or to reduce stability

uncertainty34. This minimizes the number of required cutting tests and the amount of required machine time. Studies have shown that highly productive cutting parameters can be selected in

less than 10 cuts32. Second, because these methods learn the underlying system parameters, they can theoretically be used to calculate the stability map for different conditions, e.g., at a

different radial width than the original tests or across multiple machines34. However, multiple combinations of system parameters generate similar stability maps, so it may not be possible

to uniquely identify the true system based solely on stability results. To date, this has not been fully investigated. The Bayesian models have, so far, focused on stability and have not

considered other constraints on machining parameters. This paper describes a physics-based Bayesian approach for statistically simulating machining operations under a state of uncertainty

and establishing probabilistic predictions about the process. It extends prior research by considering other factors about the machining process in addition to stability. These predictions

are used to select a series of cutting tests that provide new information about the system and identify increasingly productive machining parameters. Compared with existing Bayesian

algorithms, this paper provides three new contributions: * 1. demonstrating how Bayesian machining modeling can be extended to include cutting force modeling and spindle power measurements

and evaluating how this affects convergence to the true system parameters; * 2. proposing a likelihood function based on aleatoric uncertainty that increases explainability and relates it

directly to the Bayesian algorithm’s predictive capabilities; and * 3. defining a new, computationally efficient test selection method designed to reduce uncertainty in the underlying system

parameters. The paper is organized as follows. “Results” Section describes the proposed Bayesian framework and describes the results of numerical and experimental case studies. “Discussion”

Section discusses the results and provides suggestions for future research. “Methods” Section provides detailed information on the case studies necessary for replication. Code for the case

studies can be found at the Github repository linked in the Code Availability Statement at the end of this paper. The supplementary information available along with this paper contains

additional figures showcasing the uncertainty reduction test selection method. RESULTS This section describes the Bayesian milling model and demonstrates it using two case studies. There are

five general steps, illustrated in Fig. 1. * 1. The user defines a prior distribution describing their initial uncertainties about the frequency response function (FRF) and cutting force

model. * 2. The prior samples are propagated through a stability map algorithm to create probabilistic predictions about the system and test cuts. * 3. The predictions are used to recommend

and perform an informative test cut based on the user’s goal, such as improving productivity or reducing system uncertainty. Testing halts if there is no recommendation. * 4. A new posterior

distribution is defined which incorporates the information gained from the previous test cut using the likelihood function. * 5. New samples are drawn from the posterior distribution and

returned to propagate through the physics algorithms and create updated predictions. This framework is modular and can be customized depending on the user’s goal and knowledge about the

system. For example, new physics algorithms can be added to simulate different aspects of the milling operation, the likelihood function can be customized depending on the type of

information gathered from each test cut, and different methods can be employed to sample the posterior distribution. These options are described in more detail in Section “Prior definition”

through “Markov Chain Monte Carlo sampling”. Section “Verification studies” shows two verification studies which demonstrate the framework, one numerical and one experimental. PRIOR

DEFINITION Let \(\theta\) be a variable that defines the machining system and provides enough information to calculate the physics models (described in Section “Physics models and

probabilistic simulation”). This variable can take many different forms, but for stability modeling \(\theta\) must at minimum define the cutting force model and FRF. This paper uses the

form shown in Eq. 1. The cutting force is modeled using a mechanistic cutting force model, where \({K}_{s}\) is the specific cutting force in N/mm2, \(\beta\) is the cutting force angle in

degrees describing the ratio between the normal and tangential cutting forces \({F}_{n}\) and \({F}_{t}\) respectively, and \({k}_{{te}}\) is the tangential edge coefficient in N/mm1. Eq. 2

shows how these are used to calculate \({F}_{t}\) and \({F}_{n}\). This omits the normal edge coefficient \({k}_{{ne}}\) since it is not relevant for calculating stability or the spindle

power. The FRF is defined as a mass-spring-damper system with natural frequency \({f}_{n}\), stiffness \(k\), and damping ratio \(\zeta\). This can also be extended by including additional

modes, as shown in the experimental case study. $$\theta =\left[{K}_{s},\beta ,{k}_{{te}},{f}_{n},k,\zeta \right]$$ (1) $${F}_{t}={K}_{s}\sin \left(\beta

\right){bh}+{k}_{{te}}b,{F}_{n}={K}_{s}\cos \left(\beta \right){bh}$$ (2) These parameters can be converted to different forms for convenient calculation. For example, the modal parameters

\({\theta }_{{f}_{n}}\), \({\theta }_{k}\), \({\theta }_{\zeta }\) can be converted into the mass-spring-damper values \({\theta }_{m}\), \({\theta }_{k}\), \({\theta }_{c}\) respectively

using Eq. 3. The force model parameters \({K}_{s}\) and \(\beta\) can be converted into the tangential and normal cutting force coefficients \({k}_{{tc}}\), \({k}_{{nc}}\) as shown in Eq. 4.

These coefficients model the tangential and normal forces shown in Fig. 2. The remainder of this paper uses these conversions to provide the most convenient form for any equation. $${\theta

}_{m}=\frac{{\theta }_{k}}{{\left(2\pi {\theta }_{{f}_{n}}\right)}^{2}},{\theta }_{c}={\theta }_{\zeta }\cdot 2\sqrt{{\theta }_{k}{\theta }_{m}}$$ (3) $${\theta }_{{k}_{{tc}}}={\theta

}_{{K}_{s}}\sin \left({\theta }_{\beta }\right),{\theta }_{{k}_{{nc}}}={\theta }_{{K}_{s}}\cos \left({\theta }_{\beta }\right)$$ (4) Once the general form for \(\theta\) has been determined,

the user defines a prior distribution \(P\left(\theta \right)\) which encodes the epistemic uncertainty for the system, e.g., by defining a mean and standard deviation for each variable.

There are many ways to set the prior. Li et al. measured the system dynamics by tap testing29, while Schmitz et al. used receptance coupling substructure analysis32 to predict the tool tip

FRF. This paper uses relatively uninformed priors to demonstrate how the Bayesian algorithm can identify the system despite having little initial information. Once the prior is established,

a set of samples can be drawn from that distribution. Each sample \({\theta }_{i}\) is a given set of parameters that define the system. Since the prior is generally a standard distribution

like a multivariate normal or uniform distribution, there are functions available to draw samples from the distribution. Figure 3 shows 4000 samples drawn from the prior distribution for the

numerical case study shown in Section 2.6.1. This specific prior is intentionally set to be very broad to help evaluate whether the Bayesian model can converge under worst-case conditions.

PHYSICS MODELS AND PROBABILISTIC SIMULATION Once the prior has been established, it is propagated through a series of physics models via Monte Carlo sampling. Each physics model takes one of

the sampled system parameters \({\theta }_{i}\) and uses it to predict the system or the outcome of a cut, such as the cut stability (stable or chatter) or the cut spindle power. In this

study, three physics algorithms are used to estimate the system frequency response function (FRF), cutting stability and chatter frequency, and spindle power. SYSTEM FREQUENCY RESPONSE The

complex-valued FRF is calculated first for a given frequency _f_ in Hz using Eq. 5, assuming the system has symmetric dynamics in the \(x\) and \(y\) directions (i.e., that

\({FR}{F}_{{xx}}={FR}{F}_{{yy}}\)). For multiple modes, the FRF is calculated for each mode and then summed. Performing this calculation for each sample drawn from the prior gives a

probabilistic prediction for the system dynamics. Figure 4 shows the prior distribution of natural frequencies for the numerical case study shown in Section “Numerical study”. In that

example, the prior is inaccurate and overestimates the natural frequency. Alternatively, the FRF could be calculated using RCSA34. In this case, the \(\theta\) variable would contain

information on the model parameters, such as the connection between the tool and holder, rather than a set of modal parameters. $${FRF}\left(\theta ,f\right)=\frac{1}{-{\left(2\pi

f\right)}^{2}{\theta }_{m}+i\left(2\pi f\right){\theta }_{c}+{\theta }_{k}}$$ (5) STABILITY PREDICTION The system parameters are then used to make predictions about the outcome of the

milling operation. A straight-line cut can be defined by its radial depth of cut \(a\), axial depth of cut \(b\), spindle speed \(n\), and feed per flute \({f}_{z}\). Let \(\omega\) be a

vector variable containing all cutting parameters for a given cut. The spindle speed is then \({\omega }_{n}\), for example. The stability of any given cut \(S\left(\theta ,\omega \right)\)

can be predicted using a variety of algorithms. This paper employs the well-established ZOA stability algorithm, which removes the time-dependent component of the stability problem by

truncating the frequency domain dynamic forces to only the mean component7. See refs. 1,7 for a full derivation of this algorithm. Compared with other stability algorithms, the ZOA is fast

to calculate but is less accurate at lower radial depths of cut16. For each spindle speed, radial depth, and set of system parameters \(\theta\), it calculates a limiting axial depth

\({b}_{lim}\left(\theta ,{\omega }_{n},{\omega }_{a}\right)\). Cuts at axial depths lower than \({b}_{lim}\) are assumed to be stable (\(S\left(\theta ,\omega \right)=1\)), while cuts above

\({b}_{lim}\) are predicted to be unstable/chatter (\(S\left(\theta ,\omega \right)=0\)). The most productive cutting parameters offer the highest material removal rate (MRR) without

exceeding the stability boundary, as shown in Fig. 5. Each sampled FRF is used to calculate a stability map and the combination of all those stability maps defines the probabilistic

stability map, where \({P}_{S}\left(\omega \right)\) is the probability of stability for a given set of cutting parameters. The probability of stability can be calculated in two ways. The

simpler method is to treat the stability predictions as deterministic and calculate what percentage of samples predict a given set of cutting parameters will be stable, as shown in Eq. 6.

$${P}_{S}\left(\omega \right)=\frac{1}{N}\mathop{\sum }\limits_{i=1}^{N}S\left({\theta }_{i},\omega \right)$$ (6) However, each stability prediction has aleatoric uncertainty, caused by

uncertainty in stability classification and assumptions in the stability algorithm. While the stability algorithm predicts that \({b}_{lim}\) is the true stability limit, \({b}_{lim}^{* }\)

may be either above or below that value, even if the system parameters \(\theta\) were known perfectly. Rather than assigning a binary stable/unstable prediction, the stability algorithm

instead assigns a likelihood of stability \({L}_{s}\left({S}_{{ex}},\omega |\theta \right)\) which states how likely a cut \(\omega\) is to experience a stable or unstable result

\({S}_{{ex}}\) for a fixed set of system parameters \(\theta\). Let \({\sigma }_{b}\) be the standard deviation associated with the predicted \({b}_{lim}\) value. There is a 68% probability

that \({b}_{lim}^{* }\) lies within the range \({b}_{lim}\pm {\sigma }_{b}\), a 95% probability that \({b}_{lim}^{* }\) is between \({b}_{lim}\pm 2{\sigma }_{b}\), and so on. The likelihood

that a given axial depth \(b\) is above \({b}_{lim}^{* }\) (and will chatter) is given by the cumulative density function (CDF) for a normal distribution with mean \(\mu\) and standard

deviation \(\sigma\) evaluated at point \(x\), as shown in Eq. (7). Figure 6a shows how this represents the true stability boundary as an uncertainty centered around the nominal predicted

stability limit. Figure 6b compares the proposed aleatoric likelihood to the sigmoid likelihood used by Akbari et al.34. and the Gaussian dropoff used by Karandikar, Schmitz, and Cornelius

et al.30,32,35. In practice, these functions behave similarly, but the CDF approach provides a more intuitive definition for the likelihood hyperparameter by associating it directly with the

aleatoric uncertainty for the stability prediction and measurement. This function can be customized depending on the specific form of the stability uncertainty, for example, to express

\({\sigma }_{b}\) as a percentage of \({b}_{lim}\) rather than as a fixed value. $$\begin{array}{c}{L}_{S}\left({S}_{{ex}},\omega ,|,\theta

\right)=\left\{\begin{array}{c}\begin{array}{cc}1-{CDF}\left(x={\omega }_{b},\mu ={b}_{lim}\left(\theta ,{\omega }_{n}\right),\sigma ={\sigma }_{b}\right) &

{\rm{Stable}}\left({S}_{{ex}}=1\right)\\ {CDF}\left(x={\omega }_{b},\mu ={b}_{lim}\left(\theta ,{\omega }_{n}\right),\sigma ={\sigma }_{b}\right) &

{\rm{Unstable}}\left({S}_{{ex}}=0\right)\end{array}\end{array}\right.\end{array}$$ (7) Equation 8 shows how the probability of stability can be calculated using the stability likelihood

function. Figure 7 compares probabilistic stability maps generated by these two approaches, showing how incorporating the aleatoric uncertainty decreases the confidence in the stability

boundary. To make the difference obvious, the aleatoric uncertainty has been set to \({\sigma }_{b}=\) 1 mm. This value is so high that the probability of stability at the stability boundary

local minima never reaches 100%, even at zero axial depth. In practice, \({\sigma }_{b}\) should be set based on the user’s expectations about how closely the stability algorithm will

predict the true behavior of the machine. $${P}_{S}\left(\omega \right)=\frac{1}{N}\mathop{\sum }\limits_{i=1}^{N}{L}_{S}\left(S=1,\omega |{\theta }_{i}\right)$$ (8) The ZOA stability

algorithm also provides a chatter frequency estimate for a given cut, \({f}_{c}\left(\theta ,\omega \right)\). Refer to ref. 32 for more details on how the dominant chatter frequency is

extracted for a given sample. Similar to the stability classification, there is aleatoric uncertainty for these chatter frequency predictions. Let \({f}_{c}\left(\theta ,\omega \right)\)

have a standard deviation \({\sigma }_{{f}_{c}}\). Based on this, the likelihood of encountering a specific chatter frequency \({{f}_{c}}_{{ex}}\) is given by Eq. 9, where

\({\mathscr{N}}\left(x,\mu ,\sigma \right)\) is the probability density of a normal distribution with mean \(\mu\) and standard deviation \(\sigma\) evaluated at a value \(x\). In practice,

the likelihood should be between 0 and 1, so the non-normalized version shown in Eq. (10) can be used instead. This is equivalent to the Gaussian dropoff approach presented by Schmitz et

al.32. $${L}_{{f}_{c}}\left({{f}_{c}}_{{ex}},\omega |\theta \right)\propto \left\{\begin{array}{c}\begin{array}{cc}{\mathscr{N}}\left({f}_{{c}_{{ex}}},{f}_{c}\left(\theta ,\omega

\right),{\sigma }_{{f}_{c}}\right) & {\rm{Unstable}}\\ 1 & {\rm{Stable}}\end{array}\end{array}\right.$$ (9) $${L}_{{f}_{c}}\left({{f}_{c}}_{{ex}},\omega |\theta

\right)=\left\{\begin{array}{c}\begin{array}{cc}\exp \left(-\displaystyle\frac{1}{2}{\left(\frac{{f}_{c}\left(\theta ,\omega \right)-{{f}_{c}}_{{ex}}}{{\sigma }_{{f}_{c}}}\right)}^{2}\right)

& {\rm{Unstable}}\\ 1 & {\rm{Stable}}\end{array}\end{array}\right.$$ (10) Equation 11 gives the probability of encountering a given chatter frequency considering both aleatoric and

epistemic uncertainty. Based on this, Fig. 8 shows a probabilistic prediction for the chatter frequency. $${P}_{{f}_{c}}\left(\omega ,{{f}_{c}}_{{ex}}\right)=\frac{1}{N}\mathop{\sum

}\limits_{i=1}^{N}{L}_{{f}_{c}}\left({{f}_{c}}_{{ex}},\omega |{\theta }_{i}\right)$$ (11) SPINDLE POWER The cutting forces on the tool govern the power that the spindle motor must output to

perform a given cut. If the motor is not able to output the required power, the spindle will stall. Due to the inertia of the rotating tool assembly, cuts are generally limited by the

average cutting power, which can be estimated based on the mechanistic force model1,6. The tangential cutting force applied to a single flute on the tool is given by Eq. (12), where

\({k}_{{tc}}\) is the tangential cutting force coefficient (see Eq. (3)) and \({k}_{{te}}\) is the tangential edge coefficent. The first term models the cutting force as proportional to the

chip area. Since the chip area is related to the volume of material removed from the part, the cutting power \({\bar{P}}_{c}\) associated with this term is given by Eq. (13), where

\({N}_{t}\) is the number of flutes the tool. $${F}_{t}={k}_{{tc}}{bh}+{k}_{{te}}b$$ (12) $${\bar{P}}_{c}\left(\omega \right)={k}_{{tc}}\cdot {MRR}={k}_{{tc}}\cdot {abn}{f}_{z}{N}_{t}$$ (13)

The second term in Eq. (11) models the edge rubbing force and is proportional to the total length of the cutting edge engaged with the part. Since the flute is moving at a surface speed of

\({v}_{c}=\frac{n\pi d}{60}\), the rubbing power while the flute is engaged with the part is \({V}_{c}\cdot {k}_{{te}}b\). In milling, each flute is only in contact with the part some

fraction of the time \({p}_{{engage}}\), given by Eq. (14), where \(d\) is the diameter of the cutting tool. Equation (15) gives the average rubbing power consumption for all the flutes on

the tool. Summing this power with the cutting power for a given set of system parameters \(\theta\) and cutting parameters \(\omega\) gives the total power consumption \(\bar{P}\left(\theta

,\omega \right)\), shown in Eq. (16). Note, however, that this equation is only applicable for stable cuts. Unstable cuts can have high amplitude vibration that causes the tool to disengage

from the workpiece, altering the force. The spindle power model in Eq. (16) also neglects tool deflection, which can cause the actual engagement of the tool to be slightly larger or smaller

than nominal. $${p}_{{engage}}=\frac{1}{2\pi }\cdot {\rm{acos}}\left(1-\frac{2a}{d}\right)$$ (14) $${\bar{P}}_{e}={k}_{{te}}b\cdot \frac{{v}_{c}}{2\pi }\cdot {N}_{t}\cdot

{\rm{acos}}\left(1-\frac{2a}{d}\right)=\frac{{k}_{{te}}{bnd}{N}_{t}}{120}{\rm{acos}}\left(1-\frac{2a}{d}\right)$$ (15) $$\bar{P}\left(\theta ,\omega

\right)={\bar{P}}_{c}+{\bar{P}}_{e}={\theta }_{{k}_{{tc}}}\cdot {MRR}\left(\omega \right)+\frac{{\theta }_{{k}_{{te}}}{\omega }_{b}{\omega

}_{n}d{N}_{t}}{120}{\rm{acos}}\left(1-\frac{2{\omega }_{a}}{d}\right)$$ (16) To incorporate the aleatoric uncertainty, let \({\sigma }_{\bar{P}}\) be the standard deviation for the spindle

power prediction. The likelihood function for experiencing a specific spindle power \({\bar{P}}_{{ex}}\) is given by the same approach as the chatter frequency prediction, shown in Eq. (17).

Alternatively, the power uncertainty can be expressed as a percentage of the predicted value \({\sigma }_{\bar{P} \% }\) using Eq. (18). Note that the spindle power likelihood can be

omitted when the test cut is unstable since the spindle power estimation is less accurate during chatter. See the numerical case study in Section “Numerical study” as an example.

$${L}_{\bar{P}}\left({\bar{P}}_{{ex}},\omega ,|,\theta \right)=\left\{\begin{array}{c}\begin{array}{cc}\exp \left(-\displaystyle\frac{1}{2}{\left(\frac{\bar{P}\left(\theta ,\omega

\right)-{\bar{P}}_{{ex}}}{{\sigma }_{\bar{P}}}\right)}^{2}\right) & {\rm{Stable}}\\ 1 & {\rm{Unstable}}\end{array}\end{array}\right.$$ (17)

$${L}_{\bar{P}}\left({\bar{P}}_{{ex}},\omega ,|,\theta \right)=\left\{\begin{array}{c}\begin{array}{cc}\exp \left(-\displaystyle\frac{1}{2}{\left(\frac{\bar{P}\left(\theta ,\omega

\right)-{\bar{P}}_{{ex}}}{\bar{P}\left(\theta ,\omega \right)\cdot {\sigma }_{\bar{P} \% }}\right)}^{2}\right) & {\rm{Stable}}\\ 1 & {\rm{Unstable}}\end{array}\end{array}\right.$$

(18) Propagating the prior samples through Eq. (16) gives a probabilistic estimate for how much spindle power any cut will require (due to uncertainty in \({k}_{{tc}}\) and \({k}_{{te}}\)).

Figure 9 gives an example, showing the distribution of possible spindle powers for one set of cutting parameters. TEST SELECTION METHODS Once the probabilistic predictions for milling

stability and spindle power have been established, they can be used to select a new cutting test based on the user’s goal. Here, two different test selection strategies are described. First,

the MRR maximizing approach presented in ref. 32 is used. This approach selects the most productive possible cutting test parameters which conform to the user’s selected risk tolerance for

encountering chatter. A risk tolerance of \({R}_{{stable}}=\) 0.75, for example, indicates that the user is willing to perform cutting tests only if there is at least a 75% probability of

the cut being stable. Based on this risk tolerance, the most productive cutting test is selected using Eq. (19), where \(\Omega\) is the set of all possible cutting parameters for the tool.

Figure 10 demonstrates this approach, showing how the test is selected using lines of constant MRR and the risk tolerance boundary. $${\omega }_{{test}}=\mathop{{\rm{arg}}\,

{\rm{max}}}\limits_{\omega \in \Omega }\left({MRR}\left(\omega \right)\right)\,{\rm{s}}.\,{\rm{t}}.\,{P}_{S}\left(\omega \right)\ge {R}_{{stable}}$$ (19) However, this MRR maximizing

approach requires a certain level of confidence to select new tests. For highly uncertain systems, there may not be any productive points with a high probability of stability. Instead, tests

can be selected to provide new information about the system. The following algorithm selects cutting tests based on the expected reduction in uncertainty for some variable, e.g., the

natural frequency for the tool tip FRF. This builds on a method described by Karandikar et al. for reducing process damping uncertainty36 by extending it to other variables, providing a

faster method for estimating the posterior uncertainty and incorporating the chatter frequency and spindle power into the uncertainty estimates. Let \(u\left(x\right)\) be the standard

deviation for some variable \(x\), calculated using the distribution of the samples taken from the prior. The conditional uncertainties \(u\left(x|S\left(\omega \right)\right)\) and

\(u\left(x|{\neg S}\left(\omega \right)\right)\) are the posterior uncertainties for \(x\) if a stable or unstable cut was observed for the cutting parameters \(\omega\), respectively. These

conditional uncertainties can be quickly approximated by weighting each prior sample according to the stability likelihood and calculating the weighted standard deviation. Equation 20 gives

the weight for each sample for the stable and unstable cases, where \({w}_{i}\) is the prior weight for the sample \(i\). The sample weights can be subsequently used to calculate the

weighted uncertainty \({\hat{\sigma }}_{w}\) using Eq. (21), where \({\mu }_{w}={\sum }_{i=1}^{N}{w}_{i}^{* }{x}_{i}/{\sum }_{i=1}^{N}{w}_{i}^{* }\) is the weighted mean for all the samples.

\(u\left(x|S\left(\omega \right)\right)\) is given by evaluating \({\hat{\sigma }}_{w}\) using the stable weights, while \(u\left(x|{\neg S}\left(\omega \right)\right)\) uses the unstable

weights. Once those conditional uncertainties have been estimated, the expected reduction in uncertainty \(E\left(\Delta u\left(x|\omega \right)\right)\) is calculated using Eq. (22).

$${w}_{i}^{* }={w}_{i}\cdot \left\{\begin{array}{c}\begin{array}{cc}{L}_{S}\left(\omega |{\theta }_{i}\right) & {\rm{Stable}}\; {\rm{case}}\\ 1-{L}_{S}\left(\omega |{\theta }_{i}\right)

& {\rm{Unstable}}\; {\rm{case}}\end{array}\end{array}\right.$$ (20) $$\begin{array}{c}{\hat{\sigma }}_{w}^{2}=\frac{\mathop{\sum}\nolimits_{i=1}^{N}{w}_{i}^{* }{\left({x}_{i}\,-\,{\mu

}_{w}\right)}^{2}}{\mathop{\sum}\nolimits_{i=1}^{N}{w}_{i}^{* }}\end{array}$$ (21) $$E\left(\varDelta u\left(x|\omega \right)\right)=u\left(x\right)-\left[u\left(x|S\left(\omega

\right)\right)\cdot {P}_{S}\left(\omega \right)+u\left(x|{{\neg }}S\left(\omega \right)\right)\cdot \left(1-{P}_{S}\left(\omega \right)\right)\right]$$ (22) This process is repeated for all

sets of cutting parameters \(\omega \in \Omega\) to determine which parameters will provide the most new information about \(x\). This method for estimating \(E\left(\varDelta

u\left(x|\omega \right)\right)\) is computationally inexpensive since it is performed using cached stability maps for each of the prior samples. With \(N=4000\) the expected reduction in

uncertainty for every set of cutting parameters can be calculated in less than a minute. This approach can be extended to consider additional observations like spindle power and chatter

frequency. First, the weight for each sample is updated to include the probability of seeing a specific chatter frequency and spindle power in addition to the stability prediction as shown

in Eq. (23). This specific formulation omits the power likelihood since the spindle power calculation described in Section “Physics models and probabilistic simulation” is less accurate for

unstable cuts. Those weights can be used to evaluate the posterior uncertainty conditioned on different test results, with \(u\left(x|S\left(\omega \right),\bar{P}\right)\) being the

uncertainty after observing a stable cut with a measured power of \(\bar{P}\) and \(u\left(x|{\neg S}\left(\omega \right),{f}_{c}\right)\) giving the uncertainty after an unstable cut with a

chatter frequency of \({f}_{c}\). $${w}_{i}^{* }\left({f}_{c},\bar{P}\right)={w}_{i}\cdot \left\{\begin{array}{c}\begin{array}{cc}{L}_{S}\left(\omega |{\theta }_{i}\right)\cdot

{L}_{\bar{P}}\left(\bar{P},\omega |\theta \right) & {\rm{Stable}}\; {\rm{case}}\\ \left(1-{L}_{S}\left(\omega |{\theta }_{i}\right)\right)\cdot {L}_{{f}_{c}}\left({f}_{c},\omega |\theta

\right) & {\rm{Unstable}}\; {\rm{case}}\end{array}\end{array}\right.$$ (23) Those uncertainties are evaluated for every chatter frequency and spindle power predicted by the prior

samples, giving \(u\left(x|\neg S\left(\omega \right),{f}_{c}\left({\theta }_{j},\omega \right)\right)\) and \(u\left(x|S\left(\omega \right),\bar{P}\left({\theta }_{j},\omega

\right)\right)\) for \(j=1:N\). The uncertainties are then averaged using Eqs. 24 and 25 to give the average expected stable and unstable uncertainties \(u\left(x,|,S\left(\omega

\right)\right)\) and \(u\left(x|{\neg S}\left(\omega \right)\right)\), respectively. The prior sample weights \({w}_{j}\) are included since a higher weight indicates a higher likelihood of

encountering that chatter frequency or spindle power. These stable and unstable uncertainties are then used to calculate the expected uncertainty reduction using Eq. (22).

$$u\left(x|S\left(\omega \right)\right)=\frac{1}{{\sum }_{j=1}^{N}{w}_{j}}\mathop{\sum }\limits_{j=1}^{N}{w}_{j}\cdot u\left(x|S\left(\omega \right),\bar{P}\left({\theta }_{j},\omega

\right)\right)$$ (24) $$u\left(x|{{\neg }}S\left(\omega \right)\right)=\frac{1}{{\sum }_{j=1}^{N}{w}_{j}}\mathop{\sum }\limits_{j=1}^{N}{w}_{j}\cdot u\left(x|{{\neg }}S\left(\omega

\right),{f}_{c}\left({\theta }_{j},\omega \right)\right)$$ (25) This calculation is significantly slower than the stability-only case since the uncertainty must be estimated for each

predicted chatter frequency and spindle power, causing the calculation time to scale with the square of the sample count. To reduce calculation time, this can be done with a subset of the

samples. Figure 11 shows how including the chatter frequency gives a higher uncertainty reduction for the natural frequency. To reduce calculation time, only 1000 samples were considered,

rather than the full 4000, and required 9 min. The supplemental material for this article provides the expected uncertainty reduction maps for different variables. Figure 11 shows an

example, demonstrating how a cutting test can be selected to reduce uncertainty for the system natural frequency. However, this method does not consider the safety of any proposed test and

may recommend cuts that could cause damage to the tool or spindle/machine. The example shown in Fig. 11 will select the test with the highest probability of instability since that gives the

highest probability of observing the chatter frequency. Instead, the user can specify a desired reduction in probability and the algorithm can select the least risky test expected to provide

at least that much information. For example, Eq. (26) selects the cutting test with the lowest axial depth that provides at least a 75% expected reduction in \(u\left({f}_{n}\right)\). This

is similar to the safe test selection method presented by Akbari et al.34, but it selects cutting tests based on the absolute amount of information provided, rather than the ratio of

informativeness compared to the most productive test. If no test provides the requested uncertainty reduction, the user can adjust their cutoff, change to a different test selection

strategy, or select a test manually. The supplementary materials for this article provide example uncertainty reduction maps to demonstrate how tests can be selected to provide information

about different variables and how obtaining more information in each test cut can help reduce uncertainty. $${\omega }_{{test}}=\mathop{{\rm{arg }}\,{\rm{min}} }\limits_{\omega \in \Omega

}\left({\omega }_{b}\right)\,{\rm{s}}.\,{\rm{t}}.\,E\left(\varDelta u\left({f}_{n}|\omega \right)\right)\ge 0.75\cdot u\left({f}_{n}\right)$$ (26) LIKELIHOOD FUNCTIONS After the cutting test

has been performed, it is used to construct a new posterior probability distribution which incorporates the new information obtained from the test. Let \(\alpha\) be a variable containing

all observed results for the cutting test, including the cut stability \({\alpha }_{S}\), chatter frequency \({\alpha }_{{f}_{c}}\), and measured spindle power \({\alpha }_{\bar{P}}\).

Equation 27 shows the posterior distribution defined using Bayes’ rule, where \(P\left(\theta |\alpha \right)\) is the conditional distribution for \(\theta\) after observing the test

result, \(P\left(\theta \right)\) is the prior probability, and \(P\left(\alpha |\theta \right)\) is the likelihood function. For a given prior probability, higher likelihoods increase the

probability of a given \(\theta\) value, while lower likelihoods decrease the probability. Note that this posterior distribution may not be normalized, i.e., \({\int }_{\theta }P\left(\theta

|\alpha \right)\,\ne\, 1\), which violates the definition of a probability distribution. This will be addressed later in Section 2.5. This equation can be chained to include an arbitrary

number of tests by using the posterior distribution for the first test as the prior for the next. Equation 28 shows the posterior after \({N}_{{test}}\) tests have been run. $$P\left(\theta

|\alpha \right)\propto P\left(\theta \right)\cdot P\left(\alpha |\theta \right)$$ (27) $$P\left(\theta |{\alpha }_{1:{N}_{{test}}}\right)=P\left({\alpha }_{{N}_{{test}}}|\theta \right)\cdot

P\left(\theta |{\alpha }_{1:{N}_{{test}}-1}\right)=P\left(\theta \right)\cdot \mathop{\prod }\limits_{i=1}^{{N}_{{test}}}P\left({\alpha }_{i}|\theta \right)$$ (28) The overall likelihood

function \(P\left(\alpha ,|,\theta \right)\) is the probability that a given test result \(\alpha\) will occur for a fixed set of system parameters \(\theta\). It can be calculated by

combining the likelihood functions derived in Section 2.2, as shown in Eq. 29. This can be extended to include any arbitrary measurements about the test cuts. The examples shown here

(stability, chatter frequency, and spindle power) utilize process measurements which are practical to measure in production with minimal instrumentation, but other measurements (e.g., tool

vibration amplitude) can be added, if available. $$P\left(\alpha |\theta \right)={L}_{S}\left({\alpha }_{S},{\alpha }_{\omega }|\theta \right)\cdot {L}_{{f}_{c}}\left({\alpha

}_{{f}_{c}},{\alpha }_{\omega }|\theta \right)\cdot {L}_{\bar{P}}\left({\alpha }_{\bar{P}},{\alpha }_{\omega }|\theta \right)$$ (29) MARKOV CHAIN MONTE CARLO SAMPLING Once the posterior

distribution is defined, new samples are drawn to establish new probabilistic predictions using Markov Chain Monte Carlo sampling (MCMC). MCMC is a common method for drawing samples from the

high-dimensional non-parametric distributions obtained from Bayesian updating37,38. Rather than sampling the target distribution directly, MCMC methods create a chain of correlated samples

by randomly “walking” around the probability space and accepting or rejecting steps based on whether they are moving to higher or lower probability regions. As the chain is executed, the

marginal distribution of samples converges toward the target distribution. Compared with alternative sampling methods like rejection sampling or Laplace approximation, MCMC is more efficient

for sampling high-dimensional distributions and can accurately sample complex non-parametric or multimodal distributions. There are a wide variety of MCMC samplers available. Previous

physics-based Bayesian stability modeling has selected different samplers. For example, Li et al. used a parallel ensemble algorithm that combines multiple chains29, while Ahmadi et al.

applied a transitional MCMC sampler which weights the output chain33. This paper uses a modified version of the parallel adaptive MCMC sampler used by Schmitz et al., which is specifically

designed for iterative Bayesian updating and gives good performance without manual tuning for different distributions32. This MCMC sampler has two main steps. First, rejection sampling is

performed using the test results on all the samples taken from the prior distribution to select the proposal distribution and the initial start point. Samples are accepted or rejected based

on the likelihood they are assigned to the actual test results, with samples that more accurately predicted the test result being more likely to be retained. After this step, the remaining

samples approximate the posterior distribution \(P\left(\theta |\alpha \right)\). However, the rejection sampling assumes that the posterior will be reasonably similar to the prior such that

the remaining samples approximate the posterior distribution. If the distribution shifts too rapidly due to the rejection of a large number of samples, the initial proposal distribution

would recommend small steps that will not efficiently explore the target distribution. This is a significant issue for this study because 1) the spindle power is also used to update the

likelihood function, causing the posterior to shift more quickly; and 2) uninformed priors are used, so only a relatively small number of samples will accurately predict the test results. A

slight modification to the rejection sampling algorithm is presented which requires the user select a minimum number of samples to retain (for example, \(m=\) 100). If many samples are

retained after rejection sampling, the samples approximate the posterior distribution. However, when many samples being rejected, it increases the size of the proposal distribution to avoid

mode collapse, at the cost of decreased acceptance rate. The following steps describe the modified rejection sampling step: * For each prior sample \({\theta }_{i}\), calculate the

likelihood for the most recent cutting test \(L\left(\omega |{\theta }_{i}\right)\). * Find the \(m\) th largest likelihood \({L}_{m}\), then calculate the relative likelihoods

\(\hat{L}\left(\omega |{\theta }_{i}\right)=\frac{L\left(\omega |{\theta }_{i}\right)}{{L}_{m}}\). * For each sample, draw a random number \({u}_{i}\) from the uniform distribution

\(U\left(\mathrm{0,1}\right)\). If \({u}_{i} \,>\, \hat{L}\left(\omega |{\theta }_{i}\right)\) then discard that sample. Since samples with likelihood ≥ 1 will always be accepted there

will always be a minimum of \(m\) samples remaining after this step. The remaining samples \({\theta }^{* }\) will approximate the posterior distribution. Next, the MCMC sampling step is

initiated, using an adaptive proposal sampler which automatically tunes the proposal distribution based on the covariance of the current chain39. To improve performance, this sampler is

parallelized using the general framework proposed by Calderhead, enabling it to propose and accept multiple samples at every step40. This process is identical to the sampler defined in ref.

32. First, the sampler is initialized: * Randomly select the MCMC starting point from the retained samples \({\theta }^{* }\) and add it to the output chain as the first point \({\theta

}_{0}\). * Calculate the covariance matrix for the retained samples \({\Sigma }_{{start}}={\rm{Cov}}\left({\theta }^{* }\right)\) and define the initial proposal distribution

\({P}_{{prop}}^{{start}}\,{=}\,{\mathscr{N}}\left({\theta }_{0},\frac{5.76}{d}\cdot {\Sigma }_{{start}}\right)\), where \(d\) is the number of elements in the \(\theta\) vector. The user

defines a desired number of unique samples \({N}_{{sample}}\) and how many samples should be proposed and accepted for every iteration of the sampler, denoted by \(j\) and \(s\)

respectively. This paper uses \({N}_{{sample}}=4000\), \(j=s=60\). While there are less than \({N}_{{sample}}\) unique samples in the output chain, the following steps are repeated: * Draw

\(j\) samples \({\theta }_{1:j}^{{cand}}\) from the proposal distribution and calculate their posterior probabilities \(P\left({\theta }_{i}^{{cand}}|\alpha \right)\) using Eq. (28). This

process can be performed in parallel to save time since it requires calculating the physics algorithms for each sample. * Create the stationary distribution \(I\left({\theta

}_{i}^{{cand}}\right)=\frac{P\left({\theta }_{i}^{{cand}}|\alpha \right)}{{\sum }_{l=0}^{j}P\left({\theta }_{i}^{{cand}}|\alpha \right)}\) for \(i=0:j\), where \({\theta

}_{0}^{{cand}}={\theta }_{{last}}\) is the last element in the output chain. This represents the relative probabilities for each sample, with higher probability samples having a higher

representation in \(I\). * Draw \(s\) samples from the stationary distribution \(I\) and append them to the output chain. If \({\theta }_{0}^{{cand}}\) is selected, then another copy of

\({\theta }_{{last}}\) is appended. Note that any given sample can be accepted multiple times, with the candidates that have the highest likelihoods being the most likely to be added. To

reduce memory usage, the samples can be weighted based on how many times they have been accepted. * Calculate the covariance matrix for the union of \({\theta }^{* }\) and the current output

chain \(\varSigma ={\rm{Cov}}\left(\theta \cup {\theta }^{* }\right)\). Define a new proposal distribution \({P}_{{prop}}{=}{\mathscr{N}}\left({\theta }_{{last}},\frac{5.76}{d}\cdot

\varSigma \right)\). It is generally preferable to perform these calculations using the logarithm of the probability distribution to provide higher accuracy41. The final set of samples are

then used to generate new probabilistic predictions. As an example, Fig. 12 shows how the process stability predictions update after observing the results of a test cut. These updated

predictions are then used to select a new cutting test, iterating until the user’s stopping criterion has been met. VERIFICATION STUDIES This section presents the numerical and experimental

studies used to demonstrate the proposed model. These studies compare the ability of the Bayesian model to converge to the true system parameters with and without including the spindle power

in the Bayesian updating NUMERICAL STUDY This section describes a numerical study which shows how the algorithm performs under optimal conditions. The simulation modeled a three-flute 12.7 

mm diameter endmill with symmetric single degree-of-freedom dynamics and an aluminum workpiece. Test cuts were selected using the MRR maximizing algorithm with a risk tolerance of 50%, and

then simulated using a time domain simulation to determine stability, chatter frequency, and spindle power. Since both the rubbing and cutting power are proportional to the spindle speed and

the axial depth of cut, performing cutting tests at only one feed rate is not sufficient to learn both the cutting and edge coefficients. Therefore, each test cut was simulated at two

different feed rates, 0.05 mm, and 0.1 mm, enabling the algorithm to learn both coefficients. See Section 4.1 for details on the time domain simulation and how stability and cutting power

were estimated from the simulation results. Table 1 describes the settings for the Bayesian model used for this study. The aleatoric uncertainties were generally set low for this case since

it was expected that the time domain simulation and model predictions would closely agree. The prior distribution was set manually, using normal distributions for each variable. These prior

uncertainties were intentionally offset from the nominal values used in the time domain simulations. Most notably, the \({f}_{n}\) mean was set to be two standard deviations higher than the

nominal value. This approach evaluates the ability of the Bayesian model to converge even when given inaccurate priors. See Table 3 for the full prior along with the posterior distributions.

The list of selected test cuts is shown in Table 2. Bayesian updating was performed using these test cuts, both with and without using spindle power learning. Figure 13 shows the prior and

posterior stability maps compared to the nominal stability boundary calculated using the ZOA stability algorithm. This example identified the stability peak in only three recommended test

cuts (six total cuts since each cut was run at two feed rates), and there is no significant difference between the two different posterior stability maps. While using spindle power did not

significantly affect the stability map, it did help the model learn the tangential cutting force model. Table 3 compares the prior and posterior uncertainties for the two cases. Including

the spindle power in the updating reduced uncertainty for both \({k}_{{te}}\) and \({k}_{{tc}}\), but did not significantly lower the \({k}_{{nc}}\) or \(\beta\) uncertainties. Additionally,

the lower cutting force uncertainty reduces the stiffness uncertainty, since the limiting axial depth is determined by both the cutting force and stiffness. The model also learns to

accurately predict the spindle cutting power, even for tests where the spindle power was not used for learning. Figure 14 compares the 1σ probabilistic predictions for the spindle power for

each cut before and after the cutting test series was run. Both the prior and no power learning cases significantly overestimate the spindle power for each test, largely due to the high

prior estimate for \({k}_{{te}}\). When spindle power was incorporated in the learning process, however, the Bayesian model accurately predicts the power requirements for each cut. Note that

the spindle power model presented in Section “Physics models and probabilistic simulation” does not perfectly match the time domain power simulations. The green dots show that model’s

nominal predictions for the spindle power for each cut. Severe chatter was encountered during test 1 (\(M=\) 155 μm at 0.1 mm feed rate) which caused the tool to disengage with the

workpiece, reducing the power requirements. Test 2 also experienced chatter but had less error since the chatter had lower amplitude (\(M=\) 29 μm at 0.1 mm feed rate). The stable cuts also

are not perfectly accurate since the spindle power model ignores tool deflection and vibration, but the predictions are generally close to the simulated values. EXPERIMENTAL STUDY An

experimental case study was performed to validate the model and demonstrate how it accommodates uncertainty in the physics algorithms. It builds on the results presented in ref. 41. Cutting

tests were performed on a 1080 mild steel workpiece on a vertical machining center. The audio and spindle power were recorded for each test cut. Specific details on the experiments and

measurements are found in Section “Experimental study”. Table 4 lists the settings for the Bayesian model. The stability aleatoric uncertainty was set higher than during the numerical study

since the tool had a 2 mm corner radius which can shift the stability boundary. The prior uncertainties were set manually and intended to reflect a worst-case scenario case where the user

has very little information about the machining system and has not performed any measurements for the cutting forces or dynamics. The cutting force prior is set based on author's

experience machining steel, and in this case, significantly underestimates the true cutting forces. The dynamics prior is set as a uniform distribution which assumes the system has two modes

with natural frequencies between 1000 and 2000 Hz and similar stiffness and damping, but otherwise provides no information about the modes. However, if the priors for both modes are

identical then the Bayesian algorithm will find two solutions for the model, one where \({f}_{n1} \,>\, {f}_{n2}\) and another where \({f}_{n1} \,<\, {f}_{n2}\). To avoid this, the

prior was multiplied by the penalty function shown in Eq. 30, which constrains the system such that \({f}_{n1}\ge {f}_{n2}\). See Table 6 for the full prior along with the posterior

distributions. $${P}_{{penalty}}\left(\theta \right)=1-\exp \left(-\frac{1}{2}{\left(\frac{\max \left({\theta }_{{f}_{n2}}-{\theta }_{{f}_{n1}},0\right)}{100}\right)}^{2}\right)$$ (30)

First, a series of test cuts were selected and performed without using the spindle power in the Bayesian updating. As with the numerical study, cutting tests were performed at two feed rates

to effectively learn both the cutting and edge force coefficients. Initially, tests were selected using the MRR maximizing algorithm. However, after the first two selected test points (four

total test cuts at different feed rates) the MRR maximizing algorithm was unable to suggest a new test due to high \({f}_{n1}\) uncertainty. While \({f}_{n2}\) was identified by observing

the chatter frequency, \({f}_{n1}\) was not observed and the algorithm must infer \({f}_{n1}\) from the fact that \({f}_{n2}\) has a lower limiting axial depth at 7000 rpm, resulting in a

complex multimodal distribution as shown in Fig. 15. Since there is no clear stability peak, the MRR maximizing algorithm cannot recommend a new test. Therefore, the next cutting test was

selected to reduce uncertainty in \({f}_{n1}\). After this, the algorithm was confident enough about the location of the stability peak to switch back to MRR maximization until the stopping

criteria was met. Table 5 shows the final test cut series. Next, updating was performed using the same series of test cuts including the measured spindle power in the likelihood function.

Since the same series of test cuts are used, this establishes how the spindle power measurement helps learn the cutting force model. Two different power uncertainties were used: a low

uncertainty case with \({\sigma }_{\bar{P}}=\) 5 W based on the noise in the power measurement, and a high uncertainty case with \({\sigma }_{\bar{P} \% }=\) 25% based on the uncertainty of

the power prediction. When combined with no spindle power measurement, these three cases give the posterior probabilistic stability maps shown in Fig. 16. Note that they all converged to

similar stability predictions. Table 6 shows the prior and posterior distributions for each case. The posterior distributions are calculated by calculating the mean and standard deviation of

the samples drawn using the MCMC method. The nominal values for the tool dynamics and force model were measured using modal tap testing and cutting force dynamometer tests, respectively.

The no-spindle power case underestimates \({k}_{{tc}}\) by 44%. This is because the true stability limit is significantly higher than predicted by the ZOA algorithm, causing the Bayesian

algorithm to increase the system stiffness and decrease the cutting force to compensate. This inaccuracy in the nominal prediction is likely due to approximations in the ZOA model. First,

the ZOA does not consider the tool’s 2 mm corner radius. Due to that corner radius, the actual engagement and forces between the tool and workpiece are lower than expected. Second, the ZOA

ignores process damping, where vibrational energy is dissipated due to interference between the tool and workpiece as they vibrate relative to each other42. Both these factors can cause the

ZOA to underestimate the true stability limit. Furthermore, the no-spindle power case entirely failed to learn the edge coefficient \({k}_{{te}}\) since it has no effect on stability or

chatter frequency. In contrast, the low uncertainty power case accurately learned \({k}_{{tc}}\) and \({k}_{{te}}\). However, the posterior for the system stiffness is less accurate than the

original. When the spindle power is measured, the cutting force cannot be artificially decreased to raise the stability limit, so the stiffness must be increased further. The high

uncertainty case was between the other cases. The high uncertainty for the power allowed the Bayesian algorithm to artificially decrease the power, albeit not as much as when power was

omitted entirely. None of the cases accurately identified \({f}_{n1}\) since no chatter frequency was observed for that mode. This could occur either because every test cut could have lined

up with a stability peak, or because mode 1 could be much stiffer than mode 2 such that it has a higher limiting axial depth even at its stability troughs. In this case, the second scenario

is correct, where mode 1 is 3.8 times stiffer than mode 2. However, the prior assumes that both modes have similar stiffnesses and therefore finds a natural frequency that aligns the

stability peak with the tests. The spindle speed for the \(j\) th stability peak for some given natural frequency can be calculated using Eq. 31. For a natural frequency of 1648 Hz the 5th

stability peak occurs at roughly 6600 rpm, allowing it to plausibly encompass all the test cuts. $${n}_{{peak}}={f}_{n}\frac{60}{j{N}_{t}}$$ (31) While the low power uncertainty case learned

the cutting force coefficients, it was not able to accurately predict the cutting test results due to high uncertainty in the power predictions themselves. Figure 17 compares the measured

power and nominal predictions based on the measured cutting force coefficients to the probabilistic power predictions for each case. On average, there is a 25% absolute deviation between the

measured spindle power and the nominal prediction for each test cut, with different cases both over and underestimating the power. This error is significantly larger than the numerical case

study due to the 2 mm corner radius on the tool. The cutting force model presented in Section “Physics models and probabilistic simulation” assumes that the tool has a square corner and

thus does not accurately predict the spindle power in all cases. The low uncertainty case is not able to accommodate this model uncertainty and confidently predicts a narrow range of

possible values for each test, often not including the true value. In contrast, the high uncertainty case captures the full range of test outcomes, with six of the 10 test cuts appearing

within the 1σ 68% confidence interval. DISCUSSION The case studies demonstrate how the Bayesian model can be extended to include other predictions and inputs. Incorporating the chatter

frequency and measured spindle power helps learn the true tangential cutting force model and system natural frequency from a small number of cutting tests, which is not possible when using

only binary stable/unstable test results. Cutting tests can also be selected to provide maximum information about the underlying system, allowing the user to reduce uncertainty as quickly as

possible. However, other parameters like the normal force model, stiffness, and damping are relatively difficult to learn through test cuts. The normal cutting forces have little effect on

stability or spindle power and thus cannot be directly obtained from the data collected in this study. The stiffness can be inferred from the limiting axial depth, but determining this

accurately requires many cutting tests and is susceptible to bias from physics model assumptions. In the future, this could be resolved by incorporating additional signals to the Bayesian

updating, for example, the axis drive currents or measured surface location error. In addition to identifying system parameters, the Bayesian approach can make probabilistic predictions

about test outcomes such as the spindle cutting power. Incorporating both epistemic and aleatoric uncertainty into the predictions accommodates uncertainty in measurements and physics models

and establishes reasonable uncertainty bounds, even when there is high model uncertainty as in the experimental case study. This is critical for safe test selection since it can allow users

to fully understand the risks associated with performing a given test cut, e.g., to allow test selection to be constrained based on the machine tool’s available spindle power. However, this

method of incorporating model assumptions into the aleatoric uncertainty has some limitations. First, it assumes that physics uncertainties are normally distributed and cannot account for

consistent bias in predictions. This can result in skewed predictions for the system parameters, e.g., in the experimental case study where the ZOA algorithm consistently underpredicted the

stability limit. Furthermore, because the Bayesian model is solely physics-based, the assumptions of the physics algorithm constrain its predictions. For example, in the experimental case

study, the stability boundary depended on feed rate, but since the ZOA algorithm assumes that stability is independent of feed it fits an approximate stability boundary to both results.

Similarly, unlike more general machine learning algorithms (e.g., neural networks), the model uncertainty cannot be reduced as more tests are added, as shown in the experimental case study

where there was significant uncertainty about the power prediction for a given test even after it has already been performed and observed. Future work can investigate new ways to model or

reduce bias, for example by incorporating more accurate machining simulations (e.g., high-accuracy time domain simulations43 with non-linear cutting force models44 and process damping42),

automatically learning the aleatoric uncertainty from the test cut results by incorporating \({\sigma }_{\bar{P}}\) into the \(\theta\) vector, incorporating learnable bias terms separate

from the likelihood function, and automatically learning the aleatoric uncertainties from test results. METHODS This section describes details of the experimental verification necessary to

replicate the verification studies. Both studies used the Bayesian updating model described in Section “Prior definition” through “Markov Chain Monte Carlo sampling” to iteratively select

test cuts and update the prior distribution of system parameters, repeating until the stopping criteria were reached. Full MATLAB code for both case studies is available on GitHub under an

MIT license: see the Code Availability Statement at the end of this article. A brief introduction to the codebase and detailed descriptions on implementation and optimization can be found in

ref. 41. The remainder of this section describes specific experimental details for each study. NUMERICAL STUDY The numerical study simulated a three flute 12.7 mm diameter endmill cutting

an aluminum workpiece using climb (down) milling. Each test was simulated using time domain simulation1. This approach estimates the instantaneous cutting forces on the tool at discrete

timesteps. The cutting force is calculated using the mechanistic force model based on the current engagement between the tool and workpiece, factoring in the current vibrational displacement

of the tool. Those forces are then used to estimate the vibration of the tool at the next timestep estimated using Euler integration. This process is repeated until the desired simulation

length has been achieved, at which point the simulated displacements and cutting forces can be analyzed. The simulation used a timestep of 1/30,000th of a second and was run for 100 tool

revolutions (e.g., if the tool was rotating at 100 rotations per second then the simulation would run for one second and have 30,000 total timesteps). These settings were selected to provide

sufficient temporal resolution to fully capture the tool vibration while having enough length to reach steady-state vibration. The time domain simulation function is included in the

published MATLAB code. Each time domain simulation was classified as stable or unstable using the periodic sampling \(M\)-metric shown in Eq. 32 with a cutoff of \(M=\) 1 μm45. In that

formula, \({N}_{s}\) is the total number of once-per-rev samples drawn from the signal and \({x}_{s}(i)\) is the \(i\) th displacement sample. This approach samples from the tool’s

displacement signal once-per revolution and calculates the average absolute difference between subsequent samples. In stable cuts, the tool experiences forced vibration which repeats with

every flute pass, resulting in a low \(M\) value. Unstable cuts show variation from flute to flute and have high \(M\) values, indicating that the tool is experiencing self-excited

vibration. Figure 18 compares the simulated deflections for stable and unstable milling conditions. In Fig. 18a the once-per-revolution samples in the stable cut repeat nearly perfectly

every revolution. This indicates that the tool is experiencing forced vibration, and the cut is therefore classified as stable. In contrast, Fig. 18b shows a case where self-excited

vibrations cause the tool to vibrate at a frequency which is not a harmonic of the flute passing frequency. As a result, the once-per-revolution samples do not repeat and the cut is

classified as unstable. $$M=\frac{1}{{N}_{s}-1}\mathop{\sum }\limits_{i=2}^{{N}_{s}}\left|{x}_{s}\left(i\right)-{x}_{s}\left(i-1\right)\right|$$ (32) EXPERIMENTAL STUDY In the experimental

case study, test cuts were performed on a vertical machining center to demonstrate the proposed framework. Table 7 details specific testing conditions. Each test was recorded using a

microphone located inside of the machine enclosure to obtain the chatter frequency and help classify the test stability. Spindle power was measured using a Rogowski coil setup to measure the

spindle motor current and voltage. Figure 19a, b shows the microphone and Rogowski coil setups inside of the machine tool. For each test, the spindle cutting power was calculated by

subtracting the average power before the tool entered the cut from the average steady-state power during the machining cut, as shown in Fig. 19c. DATA AVAILABILITY Data from the experimental

test cuts and code for the numerical and experimental studies are available in the GitHub repository found at: https://doi.org/10.5281/zenodo.13381781. CODE AVAILABILITY The underlying code

for this study is available in GitHub and can be accessed via this link: https://doi.org/10.5281/zenodo.13381781. REFERENCES * Schmitz, T. L. & Smith, K. S. Machining dynamics:

frequency response to improved productivity. (Springer International Publishing, Cham, 2019). https://doi.org/10.1007/978-3-319-93707-6. * Altintas, Y., Stepan, G., Budak, E., Schmitz, T.

& Kilic, Z. M. Chatter stability of machining operations. _J. Manuf. Sci. Eng._ 142, 110801, https://doi.org/10.1115/1.4047391 (2020). Article  Google Scholar  * Quintana, G. &

Ciurana, J. Chatter in machining processes: a review. _Int. J. Mach. Tools Manuf._ 51, 363–376 (2011). Article  Google Scholar  * Wu, G. et al. A state-of-art review on chatter and geometric

errors in thin-wall machining processes. _J. Manuf. Process._ 68, 454–480 (2021). Article  Google Scholar  * Yan, B., Hao, Y., Zhu, L. & Liu, C. Towards high milling accuracy of turbine

blades: a review. _Mech. Syst. Signal Process._ 170, 108727 (2022). Article  Google Scholar  * Rai, J. K., Brand, D., Slama, M. & Xirouchakis, P. Optimal selection of cutting parameters

in multi-tool milling operations using a genetic algorithm. _Int. J. Prod. Res._ 49, 3045–3068 (2011). Article  Google Scholar  * Altintaş, Y. & Budak, E. Analytical prediction of

stability lobes in milling. _CIRP Ann._ 44, 357–362 (1995). Article  Google Scholar  * Dong, X. & Qiu, Z. Stability analysis in milling process based on updated numerical integration

method. _Mech. Syst. Signal Process._ 137, 106435 (2020). Article  Google Scholar  * Hajdu, D., Borgioli, F., Michiels, W., Insperger, T. & Stepan, G. Robust stability of milling

operations based on pseudospectral approach. _Int. J. Mach. Tools Manuf._ 149, 103516 (2020). Article  Google Scholar  * Ozoegwu, C. & Eberhard, P. A literature review on prediction

methods for forced responses and associated surface form/location errors in milling. _J. Vib. Eng. Technol._ 12, 5905–5934 (2024). Article  Google Scholar  * Schmitz, T. L. & Mann, B. P.

Closed-form solutions for surface location error in milling. _Int. J. Mach. Tools Manuf._ 46, 1369–1377 (2006). Article  Google Scholar  * Packard, M. D. & Clark, B. B. Mitigating

versus managing epistemic and aleatory uncertainty. _AMR_ 45, 872–876 (2020). Article  Google Scholar  * Özşahin, O., Özgüven, H. N. & Budak, E. Analysis and compensation of mass loading

effect of accelerometers on tool point FRF measurements for chatter stability predictions. _Int. J. Mach. Tools Manuf._ 50, 585–589 (2010). Article  Google Scholar  * Cao, H., Li, B. &

He, Z. Chatter stability of milling with speed-varying dynamics of spindles. _Int. J. Mach. Tools Manuf._ 52, 50–58 (2012). Article  Google Scholar  * Sarıkaya, M. et al. A state-of-the-art

review on tool wear and surface integrity characteristics in machining of superalloys. _CIRP J. Manuf. Sci. Technol._ 35, 624–658 (2021). Article  Google Scholar  * Honeycutt, A. &

Schmitz, T. L. Milling bifurcations: a review of literature and experiment. _J. Manuf. Sci. Eng._ 140, 120801 (2018). Article  Google Scholar  * Yue, C., Gao, H., Liu, X., Liang, S. Y. &

Wang, L. A review of chatter vibration research in milling. _Chin. J. Aeronautics_ 32, 215–242 (2019). Article  Google Scholar  * Zhu, L. & Liu, C. Recent progress of chatter

prediction, detection and suppression in milling. _Mech. Syst. Signal Process._ 143, 106840 (2020). Article  Google Scholar  * Kawai, K. et al. A prediction model for high efficiency

machining conditions based on machine learning. _Procedia CIRP_ 101, 54–57 (2021). Article  Google Scholar  * Greis, N. P., Nogueira, M. L., Bhattacharya, S., Spooner, C. & Schmitz, T.

Stability modeling for chatter avoidance in self-aware machining: an application of physics-guided machine learning. _J. Intell. Manuf._ 34, 387–413 (2023). Article  Google Scholar  *

Postel, M., Bugdayci, B. & Wegener, K. Ensemble transfer learning for refining stability predictions in milling using experimental stability states. _Int. J. Adv. Manuf. Technol._ 107,

4123–4139 (2020). Article  Google Scholar  * Postel, M., Bugdayci, B., Kuster, F. & Wegener, K. Neural network supported inverse parameter identification for stability predictions in

milling. _CIRP J. Manuf. Sci. Technol._ 29, 71–87 (2020). Article  Google Scholar  * Deng, C. et al. Improved milling stability analysis for chatter-free machining parameters planning using

a multi-fidelity surrogate model and transfer learning with limited experimental data. _Int. J. Prod. Res._ 62, 1126–1143 (2024). Article  Google Scholar  * Karandikar, J. M., Schmitz, T. L.

& Abbas, A. E. Application of Bayesian inference to milling force modeling. _J. Manuf. Sci. Eng._ 136, 021017 (2014). Article  Google Scholar  * Kong, D., Zhu, J., Duan, C., Lu, L.

& Chen, D. Bayesian linear regression for surface roughness prediction. _Mech. Syst. Signal Process._ 142, 106770 (2020). Article  Google Scholar  * Cao, L., Zhang, X. M., Huang, T.

& Ding, H. Online monitoring machining errors of thin-walled workpiece: a knowledge embedded sparse bayesian regression approach. _IEEE/ASME Trans. Mechatron._ 24, 1259–1270 (2019).

Article  Google Scholar  * Sun, H. et al. In-situ prediction of machining errors of thin-walled parts: an engineering knowledge based sparse Bayesian learning approach. _J. Intell. Manuf._

35, 387–411 (2024). Article  Google Scholar  * Karandikar, J., Traverso, M., Abbas, A. & Schmitz, T. Bayesian inference for milling stability using a random walk approach. _J. Manuf.

Sci. Eng._ 136, 031015 (2014). Article  Google Scholar  * Li, K. et al. Bayesian uncertainty quantification and propagation for prediction of milling stability lobe. _Mech. Syst. Signal

Process._ 138, 106532 (2020). Article  Google Scholar  * Cornelius, A., Karandikar, J., Gomez, M. & Schmitz, T. A Bayesian framework for milling stability prediction and reverse

parameter identification. _Procedia Manuf._ 53, 760–772 (2021). Article  Google Scholar  * Chen, G., Li, Y., Liu, X. & Yang, B. Physics-informed Bayesian inference for milling stability

analysis. _Int. J. Mach. Tools Manuf._ 167, 103767 (2021). Article  Google Scholar  * Schmitz, T., Cornelius, A., Karandikar, J., Tyler, C. & Smith, S. Receptance coupling substructure

analysis and chatter frequency-informed machine learning for milling stability. _CIRP Ann._ 71, 321–324 (2022). Article  Google Scholar  * Ahmadi, K. Bayesian updating of modal parameters

for modeling chatter in turning. _CIRP J. Manuf. Sci. Technol._ 38, 724–736 (2022). Article  Google Scholar  * Akbari, V. O. A., Kuffa, M. & Wegener, K. Physics-informed Bayesian machine

learning for probabilistic inference and refinement of milling stability predictions. _CIRP J. Manuf. Sci. Technol._ 45, 225–239 (2023). Article  Google Scholar  * Cornelius, A.,

Karandikar, J., Tyler, C. & Schmitz, T. Process damping identification using Bayesian learning and time domain simulation. _J. Manuf. Sci. Eng._ 146, 081002 (2024). Article  Google

Scholar  * Karandikar, J., Tyler, C. & Schmitz, T. Process damping coefficient identification using Bayesian inference. _Trans. North Am. Manuf. Res. Inst. SME_ 41, 55–65 (2013). Google

Scholar  * Andrieu, C., De Freitas, N., Doucet, A. & Jordan, M. I. An introduction to MCMC for machine learning. _Mach. Learn._ 50, 5–43 (2003). Article  Google Scholar  * Handbook for

Markov Chain Monte Carlo. (Taylor & Francis, Boca Raton, 2011). https://doi.org/10.1201/b10905 * Haario, H., Saksman, E. & Tamminen, J. An adaptive metropolis algorithm. _Bernoulli_

7, 223 (2001). Article  MathSciNet  Google Scholar  * Calderhead, B. A general construction for parallelizing Metropolis−Hastings algorithms. _Proc. Natl Acad. Sci. USA_ 111, 17408–17413

(2014). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Cornelius, A. Milling stability map identification and machining parameter optimization using Bayesian inference.

(University of Tennessee Knoxville, 2024). https://trace.tennessee.edu/utk_graddiss/10108 * Tyler, C. T. & Schmitz, T. L. Analytical process damping stability prediction. _J. Manuf.

Process._ 15, 69–76 (2013). Article  Google Scholar  * Dombovari, Z., Sanz-Calle, M. & Zatarain, M. The basics of time-domain-based milling stability prediction using frequency response

function. _JMMP_ 4, 72 (2020). Article  CAS  Google Scholar  * Rubeo, M. A. & Schmitz, T. L. Mechanistic force model coefficients: a comparison of linear regression and nonlinear

optimization. _Precis. Eng._ 45, 311–321 (2016). Article  Google Scholar  * Honeycutt, A. & Schmitz, T. L. A new metric for automated stability identification in time domain milling

simulation. _J. Manuf. Sci. Eng._ 138, 074501 (2016). Article  Google Scholar  Download references ACKNOWLEDGEMENTS This study was funded by the Department of Energy under grant

DE-EE0009400. The funder played no role in study design, data collection, analysis, and interpretation of data, or the writing of this manuscript. AUTHOR INFORMATION AUTHORS AND AFFILIATIONS

* Department of Mechanical, Aerospace, and Biomedical Engineering, University of Tennessee, Knoxville, Knoxville, TN, USA Aaron Cornelius & Tony Schmitz * Oak Ridge National Laboratory,

Oak Ridge, TN, USA Jaydeep Karandikar & Tony Schmitz Authors * Aaron Cornelius View author publications You can also search for this author inPubMed Google Scholar * Jaydeep Karandikar

View author publications You can also search for this author inPubMed Google Scholar * Tony Schmitz View author publications You can also search for this author inPubMed Google Scholar

CONTRIBUTIONS A.C., T.S., and J.K. collaborated on designing the overall Bayesian learning model. A.C. wrote code, designed and performed experiments, and wrote the initial manuscript draft.

J.K. and T.S. edited the draft. All authors read and approved the final manuscript. CORRESPONDING AUTHOR Correspondence to Aaron Cornelius. ETHICS DECLARATIONS COMPETING INTERESTS The

authors declare no competing interests. ADDITIONAL INFORMATION PUBLISHER’S NOTE Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional

affiliations. SUPPLEMENTARY INFORMATION SUPPLEMENTARY MATERIAL RIGHTS AND PERMISSIONS OPEN ACCESS This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives

4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original

author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted

material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise

in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the

permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/. Reprints and

permissions ABOUT THIS ARTICLE CITE THIS ARTICLE Cornelius, A., Karandikar, J. & Schmitz, T. Bayesian stability and force modeling for uncertain machining processes. _npj Adv. Manuf._ 1,

11 (2024). https://doi.org/10.1038/s44334-024-00011-y Download citation * Received: 27 August 2024 * Accepted: 18 November 2024 * Published: 18 December 2024 * DOI:

https://doi.org/10.1038/s44334-024-00011-y SHARE THIS ARTICLE Anyone you share the following link with will be able to read this content: Get shareable link Sorry, a shareable link is not

currently available for this article. Copy to clipboard Provided by the Springer Nature SharedIt content-sharing initiative