## 2 Mathematical Model  
Without loss of generality, the milling dynamic of the process with multiple and distributed delays, taking into account the regeneration effects, can be expressed by an n-dimensional linear DDE [7; 9; 32; 50] as  
\[\dot{\mathbf{y}}(t)\ =\ \mathbf{A}\mathbf{y}(t)+\mathbf{B}(t)\mathbf{y}(t)-\sum_{j\ =\ 1}^{N}\int_{a_{1}}^{a_{2}}\mathbf{C}_{j}(t,\phi)\mathbf{y}(t-\theta_{j}(t,\phi))d\phi \tag{1}\]  
where \(A\) is a constant matrix, \(N\) is the number of the distributed delays, and \(\mathbf{B}(t)\), \(\mathbf{C}_{j}(t,\phi)\) and \(\theta_{j}(t,\phi)\) are periodic coefficient matrices that respectively satisfy \(\mathbf{B}(t)\ =\ \mathbf{B}(t+T)\), \(\mathbf{C}_{j}(t)\ =\ \mathbf{C}_{j}(t+T)\) and \(\theta_{j}(t,\phi)\ =\ \theta_{j}(t+T,\phi)\) with the system period \(T\) and \(j\) = 1,2,..., \(N\). Starting with the idea in Ref. [7], Equation (1) can be written in  
\[\dot{\mathbf{y}}(t)\ =\ \mathbf{A}\mathbf{y}(t)+\mathbf{B}(t)\mathbf{y}(t)-\sum_{j\ =\ 1}^{N}\sum_{k\ =\ 1}^{f}\widetilde{\mathbf{C}}_{j,k}(t)\mathbf{y}(t-\widetilde{\theta}_{j,k}(t)) \tag{2}\]  
where  
\[\widetilde{\mathbf{C}}_{j,k}(t)\ =\ \left\{ \begin{array}{cc}\int_{a_{1}}^{a_{1}+ b\mathbf{t}}\mathbf{C}_{j}(t,\phi)d\phi&k\ =\ 1,2,\ldots,f-1\\ \int_{a_{1}+(f-1)h}^{a_{2}}\mathbf{C}_{j}(t,\phi)d\phi&k\ =\ f\\ \widetilde{\theta}_{j,k}(t)\ =\ \left\{\begin{array}{cc}\theta_{j}(t,(k-\frac{1}{2})h) &k\ =\ 1,2,\ldots,f-1\\ \theta_{j}(t,(a_{2}+fh-h)/2)&k\ =\ f\\ \end{array}\right.\end{array} \right. \tag{3}\]  
Firstly, one can divide the system period \(T\) into \(l\) equal parts, and obtain \(l\) number of time intervals with length \(\Delta t\ =\ T/l\). Then, for the \(i-\)th interval \([t_{i},t_{i+1}]\), the ODE can be represented as [7]  
\[\dot{\mathbf{y}}(t)\ =\ \mathbf{A}\mathbf{y}(t)+\widetilde{\mathbf{B}}(t)\widetilde{\mathbf{y}}(t)- \sum_{j\ =\ 1}^{N}\sum_{k\ =\ 1}^{f}\widetilde{\mathbf{C}}_{j,k}(t)\widetilde{\mathbf{y}}(t- \widetilde{\theta}_{j,k,i}). \tag{4}\]  
where \(\widetilde{\mathbf{B}}(t)\) and \(\widetilde{\mathbf{C}}_{j,k}(t)\) respectively are the linear interpolation approximations of \(\mathbf{B}(t)\) and \(\widetilde{\mathbf{C}}_{j,k}(t)\) on boundary \(t\ =\ t_{i}+\xi\ (\xi\in[0,\ \Delta t])\), and they can be defined by  
\[ \begin{array}{cc}\widetilde{\mathbf{B}}(t)\ =\ \widetilde{\mathbf{B}}(\xi+t_{i})\ =\ \mathbf{B}_{i}+\frac{\mathbf{B}_{i+1}-\mathbf{B}_{i}}{ \Delta t}\underline{\xi},\\ \widetilde{\widetilde{C}}_{j,k}(t)\ =\ \widetilde{\widetilde{C}}_{j,k}(\xi+t_{i})= \widetilde{\mathbf{C}}_{j,k,i}+\frac{\widetilde{\mathbf{C}}_{j,k,i+1}-\widetilde{ \mathbf{C}}_{j,k,i}}{\Delta t}\underline{\xi}.\end{array}  \tag{5}\]  
and the delay item \(\widetilde{\theta}_{j,k,i}\) is equal to \(\widetilde{\theta}_{j,k,i}\ =\ \int_{t_{i}}^{t_{i+1}}\widetilde{\theta}_{j,k}(t)dt/\Delta t\). When the number of intervals corresponding to \(\widetilde{\theta}_{j,k,i}\) is defined by \(m_{j,k,i}\ =\ \text{ceil}(\widetilde{\theta}_{j,k,i}/\Delta t)\), the time-delayed item \(\widetilde{\mathbf{y}}(t-\widetilde{\theta}_{j,k,i})\) in Equation (2) can be obtained by Lagrange quadratic interpolation [15; 32] at node points \(\mathbf{y}_{i-m_{j,k,i}}\), \(\mathbf{y}_{i-m_{j,k,i}+1}\) and \(\mathbf{y}_{i-m_{j,k,i}+2}\), i.e.,  
\[\widetilde{\mathbf{y}}(t-\widetilde{\theta}_{j,k,i})\ =\ \widetilde{\mathbf{y}}(\xi+t_{i}- \widetilde{\theta}_{j,k,i})\ =\ \mathbf{L}_{j,k,i}(\xi)\mathbf{y}_{i-m_{j,k,i}+M_{j,k,i}(\xi)}\mathbf{y}_{i-m_{j,k,i}+1}+N_{j,k,i}(\xi)\mathbf{y}_{i-m_{j,k,i}+2}, \tag{6}\]where  
\[ \begin{array}{l}\mathbf{L}_{j,k,j}(\xi)\ =\ \frac{(\xi-v_{j,k,j}\Delta t)*(\xi-v_{j,k,j} \Delta t-\Delta t)\mathbf{1}}{2\Delta t^{2}},\\ \mathbf{M}_{j,k,i}(\xi)\ =\ \frac{(\xi+w_{j,k,j}\Delta t)*(\xi-v_{j,k,j}\Delta t- \Delta t)\mathbf{1}}{\Delta t^{2}},\\ \mathbf{N}_{j,k,i}(\xi)\ =\ \frac{(\xi-\eta_{j,k,i}\Delta t)*(\xi+w_{j,k,j}\Delta t )\mathbf{1}}{2\Delta t^{2}},\\ w_{j,k,i}\ =\ m_{j,k,i}-\widetilde{\tau}_{j,k,i}/\Delta t,\ v_{j,k,i}\ =\ 1-w_{j,k,i}. \end{array}