SDM for Milling Dynamics & Chatter Control

• The paper presents SDM, a numerical tool that rigorously computes stability lobes for milling processes by discretizing delay differential equations.
• It outlines an algorithmic framework that integrates matrix averaging, delayed state interpolation, and eigenvalue analysis to delineate stability boundaries.
• The method underpins modern process control strategies by enabling online SLD estimation and machine learning for surface roughness prediction.

The semi-discretization method (SDM) is a numerical technique for stability analysis of milling processes governed by delay differential equations. It enables rigorous computation of the Stability Lobe Diagram (SLD), identifying spindle-speed–depth-of-cut combinations that avoid regenerative chatter—a self-excited vibration that degrades surface quality and accelerates tool wear. SDM is foundational for process control approaches, such as adaptive controllers leveraging online SLD estimation and surface roughness prediction via machine learning (Huang et al., 22 Nov 2025).

1. Continuous Time-Delay Model of Milling Dynamics

The milling system is modeled as a two-degree-of-freedom spindle–tool assembly in the $XY$–plane. The dynamical state consists of the vector $\bm q(t)=[q_x(t)\;q_y(t)]^T$, representing tool tip displacement, along with its derivatives. Motion is governed by the coupled delay differential equation: $$\bm M\,\ddot{\bm q}(t) + \bm C\,\dot{\bm q}(t) + \bm K\,\bm q(t) = a_p\,\bm H_d(t)\,[\bm q(t)-\bm q(t-\tau)] + \bm u(t)$$ where $\bm M$, $\bm C$, and $\bm K$ are the mass, damping, and stiffness matrices respectively, $N$ the number of teeth, $\omega_{sp}$ the spindle speed, $\tau=60/(N\omega_{sp})$ the tooth-pass delay, $a_p$ the axial depth of cut, and $\bm H_d(t)$ the time-varying cutting-force coefficient matrix. External control input $\bm u(t)$ is typically zero for pure stability analysis.

Transforming to first-order form yields

$$\dot{\bm x}(t) = \bm A(t)\,\bm x(t) + \bm B(t)\,\bm x(t-\tau) + \bm D(t)\,\bm u(t)$$

where $$\bm x(t) = [\bm q(t); \dot{\bm q}(t)] \in \mathbb R^4$$, and the system matrices incorporate the mechanical and cutting-force parameters.

2. Semi-Discretization Procedure

SDM approximates the continuous-delay model by discretizing the delay $\tau$ into $m$ uniform sub-intervals of length $\epsilon$, such that $\tau = m\epsilon$. For each interval $[i\epsilon, (i+1)\epsilon]$, the approach entails:

1. Averaging system matrices: $$\bm A_i = \frac{1}{\epsilon} \int_{i\epsilon}^{(i+1)\epsilon} \bm A(t)\,dt, \quad \bm B_i = \frac{1}{\epsilon} \int_{i\epsilon}^{(i+1)\epsilon} \bm B(t)\,dt, \quad \bm D_i = \frac{1}{\epsilon} \int_{i\epsilon}^{(i+1)\epsilon} \bm D(t)\,dt$$
2. Approximating the delayed state using linear interpolation at midpoints: $$\bm x(t_i-\tau) \approx \tfrac12(\bm x_{i-m} + \bm x_{i-m+1}), \quad \bm x_i \equiv \bm x(i\epsilon)$$
3. Integrating over each sub-interval with constant delayed state, resulting in the discrete map: $$\bm x_{i+1} = \bm P_i\,\bm x_i + \tfrac12\,\bm R_i(\bm x_{i-m} + \bm x_{i-m+1}) + \bm Q_i\,\bm u_i$$ with

$$\bm P_i = \exp(\bm A_i\epsilon),\quad \bm R_i=(\bm P_i-\bm I)\bm A_i^{-1}\bm B_i,\quad \bm Q_i=(\bm P_i-\bm I)\bm A_i^{-1}\bm D_i$$

1. Forming a finite-dimensional state by augmenting past $m$ delayed states: $$\bm y_i = \operatorname{col}(\bm x_i, \bm x_{i-1}, \dots, \bm x_{i-m}) \in \mathbb R^{4(m+1)}$$ which evolves as

$\bm y_{i+1} = \bm E_i\,\bm y_i + \bm G_i\,\bm u_i$

where $\bm E_i$ is block-upper-Hessenberg incorporating $\bm P_i$, $\tfrac12\bm R_i$, and identity shifts.

2. Global transition mapping over one delay period: $$\bm y_{i+m} = \bm\Phi_i\,\bm y_i + \bm\Gamma_i\,\bm u_i, \quad \bm\Phi_i = \bm E_{i+m-1}\cdots\bm E_i$$ $\bm\Phi_i$ is the central object for stability analysis.

3. Characteristic Equation and Stability Criterion

For autonomous ($\bm u=0$), periodic-coefficient systems, one drops the interval index and analyzes $\bm\Phi$. Asymptotic stability—absence of chatter—holds iff

$$\sigma(\bm\Phi) \subset \{z\in\mathbb C : |z| < 1\}$$

i.e., all eigenvalues are inside the unit circle. The characteristic equation for SDM reads

$\det(z\,\bm I - \bm\Phi) = 0$

or equivalently, with delay incorporated,

$$\det\left[z\,\bm I_{4(m+1)} - \bm P - \tfrac12\,\bm R (z^{-m} + z^{-m-1})\right]=0$$

where $\bm P$ propagates instantaneous dynamics and $\bm R$ encodes delay feedback from previous intervals, with boundary roots $|z|=1$ demarcating the stability threshold.

4. Stability Lobe Boundary and Critical Depth of Cut

On the boundary of stability (chatter onset), the dominant eigenvalue $z=e^{j\theta}$ of $\bm\Phi$ satisfies $|z|=1$. Setting $z=e^{j\theta}$, the boundary condition becomes: $$\det\Bigl[e^{j\theta}\bm I-\bm P-\tfrac12\,\bm R(e^{-j m\theta}+e^{-j(m+1)\theta})\Bigr]=0$$ This yields, for the single-degree-of-freedom case,

$$a_p^{\rm cr}(\omega_{sp}) =\frac{2\,\omega_n\,M\,\zeta\,\sin\theta }{|\bm h_d(e^{j\theta})|}, \quad \theta = \frac{\omega_n\,\tau}{2} = \frac{\pi\,\omega_n}{N\,\omega_{sp}}$$

where $\bm h_d(e^{j\theta})$ denotes the discrete-Fourier symbol of cutting stiffness. For multi-DOF systems, $a_p^{\rm cr}$ is determined numerically via root-finding or contour interpolation of the spectral radius criterion.

5. Numerical and Algorithmic Implementation

A direct computational approach for milling stability via SDM proceeds as follows:

• Parameter sweep: Discretize $\omega_{sp}$ across $[\omega_{\min},\omega_{\max}]$ (e.g., 200 points) and $a_p$ over $[0, a_{\max}]$ (e.g., 100 points), or solve $a_p$ for fixed $\omega_{sp}$.
• Stepwise matrix formation: For each $(\omega_{sp}, a_p)$, compute $\tau$, sub-interval $\epsilon = \tau/m$ ($m \sim 20$), assemble $\bm A_i$, $\bm B_i$ via mid-point/trapezoidal rule, calculate $\bm P_i$, $\bm R_i$, and $\bm E_i$, then multiply to obtain $\bm\Phi$.
• Eigenvalue analysis: Compute maximum absolute eigenvalue $\lambda_{\max} = \max|\sigma(\bm\Phi)|$ using standard eigensolvers (MATLAB eig, Python SciPy linalg.eig).
• Stability determination: Mark parameter combination as stable if $\lambda_{\max}<1$, unstable if $\lambda_{\max}>1`; interpolate the contour$\lambda_{\max}=1$$to construct the SLD.</li> </ul> <p>Recommended tools include MATLAB (using <code>expm</code>, <code>eig</code>, <code>contour</code>), Python (NumPy/SciPy <code>expm</code>, <code>eig</code>), and matplotlib for visualization. For efficient root-finding in$$a_p$, use methods such as bisection or secant on$\lambda_{\max}(a_p)=1$.

  6. Integration with Process Control and Machine Learning

  The SDM is integral to advanced process controllers capable of real-time chatter suppression. In this context, machine learning frameworks are developed to estimate the SLD and surface roughness from sensor data online. These estimates feed into adaptive controllers that adjust spindle speed to maintain operation within the stable region defined by the SLD, maximizing surface finish and minimizing tool wear. The efficacy of such controllers is supported by simulations and experimental data demonstrating improvements over prior approaches (Huang et al., 22 Nov 2025).

  7. Research Context and Further Applications

  The semi-discretization method bridges rigorous time-delay systems analysis with pragmatic milling process stability assessment. It provides a reproducible computational foundation for researchers and practitioners aiming to delineate stability boundaries, design chatter-resistant operation schedules, and integrate online estimation schemes for industrial control. Application of SDM is central to contemporary approaches in cyberphysical manufacturing systems, especially where online adaptation is enabled by sensor-driven ML frameworks. Its formalism is extensible to multi-degree-of-freedom systems and complex tool geometries, supporting ongoing innovation in milling dynamics and chatter suppression.