# Robust Chatter Stability Prediction of the Milling Process considering Uncertain Machining Positions

Authors: Congying Deng; Wei Zhou; Kai Yang; Zhiyu Huang; Qian Tang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Milling stability is a function of the tool point frequency response functions (FRFs), which vary with the movements of the moving parts within the whole machine tool work volume. The position-dependent tool point FRFs result in uncertain prediction of the stability lobe diagram (SLD) for chatter-free machining parameter selection. Taking the variations of modal parameters to represent the variations of tool point FRFs, this paper introduces the edge theorem to predict the robust milling chatter stability. The application of the edge theorem requires the minimum and maximum modal parameters within the machining space defined by the machining position and machining allowance information. Then, radial basis function artificial neural networks (RBFANNs) are used to predict the position-dependent modal parameters inX and Y directions based on the sample information of machining positions and related modal parameters at the tool point. Moreover, sample machining spaces are determined based on the aforementioned sample positions, and the trained RBFANNs are used to obtain corresponding sample extreme modal parameters. On this basis, RBFANNs for predicting the position and machining allowance-dependent extreme modal parameters can also be trained, and they are combined with the edge theorem and zero exclusion condition to calculate robust pairs of the spindle speed (n) and limiting axial cutting depth (aplim) and then plot the robust SLD (RSLD). A case study was performed on a real three-axial vertical machining center, and the plotted RSLD considering position variations was compared with the traditional SLD. Results of the chatter tests show that the RSLD can provide more reliable (ap, n) pairs to guarantee the milling stability, validating the feasibility of the proposed robust milling chatter stability prediction method.


---

## Body

## 1. Introduction

Milling is one of the major machining processes that can obtain finished products with desired geometry, dimensions, and surface roughness by material removal. During milling operations, the most common self-induced vibration defined as regenerative chatter is a significant obstacle limiting improvements of the production efficiency and processing quality [1–3]. The regenerative chatter is caused by the forces generated during the dynamic cutting process and is not dependent on the external periodic forces. The disastrous nature of chatter vibration brings numerous negative effects including poor surface finish, aggravated tool wear, excessive noise, and damages to machine tool components [4, 5]. Therefore, extensive research studies have been conducted to prevent the chatter occurrence by analytical stability prediction, real-time detection, chatter control technique, and so on [6–8].The well-known stability lobe diagram (SLD), which provides stable and unstable cutting zones divided by pairs of limiting axial cutting depth and spindle speed, is a significant approach benefiting selections of chatter-free machining parameters, such as spindle speed, axial cutting depth, radial cutting width, and feed rate per tooth [9]. Generally, the SLD can be predicted by theoretical analyses in time and frequency domains or based on the results of chatter tests [10]. Guo et al. [11] have taken the effect of the cutter’s helix angle on the chatter into consideration and proposed a new integral interpolation method to predict the stability lobes with multidelays. Tang et al. [12] considered the high-order interpolation of the state item and the time-delay term simultaneously and proposed an updated full-discretization method for predicting the milling stability more accurately. Altintas and Budak [13, 14] developed a two-degree-of-freedom (2DOF) dynamic model for the milling vibration system and proposed an analytical chatter prediction theory to obtain the SLD by introducing the Fourier approximation method. As this proposed method was proved to be efficient in obtaining the SLD, it lays a foundation in following milling stability predictions considering various affecting factors.To analyze the milling stability using these aforementioned methods, the frequency response functions (FRFs) at the tool points are required. However, as the machine tool structure is assembled by key parts and their related joints, the tool point FRFs are a function of the moving part positions. Then, the position-dependent tool point FRFs result in uncertain prediction of milling stability. Therefore, the position-dependent machine tool dynamic behaviour and milling stability analysis have attained a lot of attention in recent years [15, 16]. Law et al. [17] developed an efficient position-dependent multibody dynamic model of a machine tool and proposed related methodology to compute the tool point FRFs based on reduced model substructural synthesis, and the rapidly calculated tool point FRFs were used to research the position-dependent milling stability. Semm et al. [18] also divided the machine tool into reduced substructures and assembled at the desired axis positions and used the local damping models and substructuring approaches to efficiently model the changing damping distribution and eigenfrequencies. Then, a state-space model was obtained and incorporated into a time domain simulation of the position-dependent dynamic behaviour. The approximate models were also used in predicting position-dependent milling stability. Zhang et al. [19] considered the moving part displacements in three axes and established a response surface model to predict the minimum position-dependent critical axial cutting depth.Majority of research studies validate that the SLD reflecting the milling stability varies within the machine tool work volume and provides methods to obtain the optimal position with chatter-free machining parameters corresponding to higher machining efficiency [20]. However, in real milling operations, deviations between the predicted and measured chatter occurring conditions are still observed. This problem may be ascribed to fact that the predicted chatter-free machining parameters are only useful at one specific position and its adjacent space. When the relative moving displacements between the tool and workpiece are large, the change of machine tool structure causes the variations of tool point FRFs, making the predicted chatter-free machining parameters at one position invalid and bringing in the milling chatter occurrence. Robust milling stability predictions considering the tool point FRF uncertainties have been proposed, but few research studies focusing on the uncertainties caused by machining position variations were addressed [21]. Graham et al. and Park and Qin [22, 23] have defined specific changes to modal parameters of the tool point FRFs considering effects of cutting force coefficients, high spindle speed, and so on and combined the edge theorem and zero exclusion condition method to obtain the robust SLD for machining parameter selection. However, only a simple ±5% or ±10% variation was taken to each modal parameter. Deng et al. [24] predicted position-dependent modal parameters by the Kriging model and obtained an optimal position with a higher axial cutting depth through particle swarm optimization. Then, at the optimal position, modal parameter variations within the machining space defined by three-directional machining allowances were determined, and the edge theorem was used to predict the robust chatter-free machining parameters. However, this method is time consuming to repeatedly calculate these position-dependent modal parameters within the machining space when the machining allowances or optimal positions change.Therefore, this paper presents a theoretical approach to fascinate rapid prediction of robust chatter-free machining parameters within the machining space determined by the machining allowances. Typical positions within whole machine tool work volume are determined by the experimental design method, and tool point FRFs at these positions are measured to identify the modal parameters. Radial basis function artificial neural networks (RBFANNs) are used to predict position-dependent modal parameters based on the information of sample positions and modal parameters. For each typical position, the trained RBFANNs are used to predict the maximum and minimum values of each dominant modal parameter within different machining spaces determined by sample machining allowances. On this basis, RBFANNs whose inputs are the information of sample positions and their related machining allowances are also trained to predict the extreme modal parameters. Then, a robust position-dependent SLD for selecting robust chatter-free machining parameters is predicted using the edge theorem and zero exclusion condition method.Henceforth, this paper is organized as follows. An analytical methodology for predicting the robust position-dependent milling stability based on the edge theorem and zero exclusion condition method is first presented in Section2. The radial basis function artificial neural network to predict the position-dependent modal parameters and their extreme values within the machining space is described in Section 3. In Section 4, a case study is performed on a real vertical machining center for illustrating application of the proposed method. This is followed by a summarized conclusion in Section 5.


## 2. The Edge Theorem Applied in Traditional Milling Chatter Stability Model

Research studies have been conducted to accurately establish and solve the chatter stability model, and the chatter stability has been shown to be more accurately analyzed by solving eigenvalue problem proposed by Budak and Altintas. Then, uncertainties in the dynamic milling process system can be implemented in this accurate chatter stability model by applying the edge theorem and zero exclusion condition.


### 2.1. DOF Milling Stability Model

A two-degree-of-freedom (DOF) dynamic model shown in Figure1(a) is generally used to describe the milling process system, in which the tool is assumed to have two orthogonal degrees of freedom in X and Y directions, and m, ξ, and k represent corresponding mass, damping, and stiffness coefficients. According to the analytical milling stability theory presented by Budak and Altintas [13, 14], dynamic variations of the cutting forces in two directions ΔFx and ΔFy can be can be described as follows:(1)ΔFxtΔFyt=apKtcN4πaxxaxyaxyayyΔxΔy=apKtcN4πAΔxΔy,where ap is the axial cutting depth, N is the tool tooth number, Ktc is the tangential cutting force coefficient, Δx and Δy are the X- and Y-directional vibration displacements, and [A] is the average directional factor which is expressed as follows:(2)axx=−cos2φjl+2Krtφjl−Krtsin2φjlφstφex,axy=2φjl+sin2φjl−Krtcos2φjlφstφex,ayx=−2φjl+sin2φjl−Krtcos2φjlφstφex,ayy=cos2φjl+2Krtφjl+Krtsin2φjlφstφex,where Krt is the ratio of radical cutting force coefficient to the tangential one and φst and φex are the start and exit angles of the cutting tooth.Figure 1
(a) The 2DOF dynamic milling system. (b) The traditional stability lobe diagram.
(a)(b)In frequency domain, vibration displacements Δx and Δy can be described based on the transfer function matrix G(iω):(3)ΔxΔy=1−e−iωτGiωFeiωt⟵Giω=GxxiωGxyiωGyxiωGyyiω,where τ is the tool tooth passing period, Gxx(iω) and Gyy(iω) are the direct tool point FRFs in X and Y directions, and Gyx(iω) and Gyx(iω) are the cross tool point FRFs in X and Y directions. Appling equation (3) into equation (1), the dynamic cutting forces are expressed using the following equation:(4)Feiωt=apKtcN4πA1−e−iωτGiωFeiωt.An eigenvalue analysis is performed on equation (4) to obtain its following characteristic equation:(5)detI−apKtcN4π1−e−iωτAGiω=0.Then, if the cross tool point FRFs in two directions are ignored, the characteristic equation can be rewritten as follows:(6)a0Λ2+a1Λ+1=0Λ=−apKtcN4π1−e−iωτ=ΛI+ΛR⟵a0=GxxiωGyyiωaxxayy−axyayxa1=axxGxxiω+ayyGyyiω,where Λ is the eigenvalue and its real and imaginary parts are named ΛR and ΛI. The critical stability point occurs at the chatter frequency ωc, and it is substituted into equation (6) to obtain the corresponding limiting axial cutting depth aplim and spindle speed n:(7)aplim=−4πΛR+iΛIKtcNcosωcτ−isinωcτ,n=60ωcN2k+1π−2arctanΛR/ΛI.Therefore, in the spindle speed range, a traditional stability lobe diagram shown in Figure1(b) can be plotted according to equation (7). The lobe curve is formed by different combinations of the limiting axial cutting depth aplim and spindle speed n, which is used to divide the stable and unstable zooms. The (ap, n) pairs below the curve correspond to stable machining conditions, and the (ap, n) pairs above the curve result in unstable machining conditions.


### 2.2. Robust Milling Stability Prediction by Edge Theorem and Zero Exclusion Condition

Equations (1) to (7) show that the SLD mainly depends on the tool point FRFs Gxx(iω) and Gyy(iω). A general FRF G(iω) with Nm modes can be described using the following equation [25]:(8)Giω=∑r=1Nmωnr/kers2+2ξrωnrs+ωnr2,where s = iω is the Laplace variable and ωnr, ξer, and ker are the rth modal natural frequency, damping ratio, and stiffness, respectively. During a milling operation, the machine tool structure and its joint contact positions vary with movements of the moving parts, causing variations in tool point FRFs and related modal parameters and then leading to uncertain SLD prediction.Considering that the modal parameters will vary within the machining space determined by the workpiece machining allowances, the edge theorem and zero exclusion condition are applied to develop a robust milling stability prediction by extending traditional milling chatter stability theory described in Section2.1.The edge theorem can be used to determine the stability of the time-delay system with uncertainties varying within a specific range defined by the minimum and maximum values. According to the edge theorem, for a polynomialP with uncertain parameters in the Laplace domain, the minimum and maximum values of each uncertain parameter can be selected to form different combinations, and a corresponding family of polynomials can be obtained by substituting these combinations into the polynomial. Calculating each extreme polynomial at a given frequency will obtain different vertices on the complex plane, and the adjacent vertices can be connected to form a polygon. If the polygon edges are all stable, robust stability of the system within the boundary of edges can be guaranteed.For applying the edge theorem to the 2DOF milling system with uncertain modal parameters, the characteristic equation of the chatter stability model expressed in equation (6) is first rewritten in the Laplace domain as follows [22–24]:(9)0=Gxxs,ωnGyys,ωnaxxayy−axyayx−NapKtc4π1−e−sτ2+axxGxxs,ωn+ayyGyys,ωn×−NapKtc4π1−e−sτ+1.Since the application of edge theorem in performing the robust stability evaluation requires the system characteristic equation to be in a polynomial form, the denominator of one dominant mode shown in equation (8) is multiplied to the left and right sides of equation (9) to yield equation (10) in a polynomial form:(10)0=Gxxs,ωnGyys,ωnaxxayy−axyayx−NapKtc4π1−e−sτ2+axxGxxs,ωn+ayyGyys,ωn−NapKtc4π1−e−sτ+1×s2+2ξiωn,is+ωn,i2.The variations of dominant modal parameters are taken to represent the machining position changes, and they are assumed to vary within its minimum and maximum values. Then, the left and right boundary values of the natural frequency, damping ratio, and stiffness variation internals inX and Y directions are defined as ωnxmin, ωnxmax, ωnymin, ωnymax, ξxmin, ξxmax, ξymin, ξymax, kxmin, kxmax, kymin, and kymax. These extreme values of modal parameters are selected to form different combinations, which are substituted into equation (10) to determine a family of polynomials P as expressed in the following equation:(11)P=p1s,v,p2s,v,…,pis,vv=ωnx,ξx,kx,ωny,ξy,ky⟸ωn∈ωnmin,ωnmax,ξ∈ξmin,ξmax,k∈kmin,kmax,where i is the total number of the combinations with extreme values and i will be 2Nv if the number of varying parameters is Nv. As one dominant mode has three modal parameters, there are 6 uncertain parameters considering X and Y directions. Thus, Nv equals 6, and there are 64 combinations containing extreme values of modal parameters. Then, values of the extreme combinations are calculated using equation (10) to produce vertices in the complex plane when a frequency is given. The adjacent vertices are connected to form the edges of a polygon, and the robust stability of the milling system is dependent on the edges according to the edge theorem. Therefore, the zero exclusion condition is introduced to efficiently evaluate the stability of each edge.The zero exclusion condition is a graphic approach to investigate the relationships among the edge locations and the origin at a given frequency in the complex plane. The system is stable if the origin is located inside the polygon formed by the edges; otherwise, the system is unstable. A system with two uncertain parameters is assumed to have a better illustration of the zero exclusion condition. The two uncertain parameters yield four vertices, which form a quadrilateral in the complex plane as shown in Figure2. The quadrant angle α of each vertex is defined, and distributions of the quadrant angles are used to discuss whether the quadrilateral encircles the origin. Five cases are summarized in Table 1 according to distributions of the quadrant angles. The symbols + and ‒ shown in Figure 2 and Table 1 mean the maximum and minimum quadrant angles in one quadrant, respectively.Figure 2
The quadrant angles in the complex plane.Table 1
Stable or unstable conditions for different distributions of the vertices in the complex plane.
Case no.Milling conditionDistributions of the vertices in the complex plane1StableAll vertices are located in the same quadrant, and the origin is outside the polygon.2StableVertices are located in two adjacent quadrants, such as I-II, II-III, III-IV, or IV-I.3StableVertices are located in diagonal quadrants such as I–III or II–IV, and the system is stable only whenαm+2+-αm_ < π or αm+2+-αm+ > π. m is the mth quadrant, and αm+2+ and αm_ are the minimum and maximum quadrant angles.4StableVertices are located in three quadrants as I-II-III, II-III-IV, III-IV-I, and IV-I-II where only the following four conditions correspond to the stable system:α3+-α1_ < π in I-II-III; α4+-α2_ < π in II-III-IV; α3_-α1+ < π in III-IV-I; α4_-α2+ < π in IV-I-II.5UnstableVertices are located in four quadrants, and the origin is inside the polygon.Therefore, the robust milling stability prediction is transformed into a graphic problem. Modal parameters corresponding to the positions within the machining space are predicted to determine the minimum and maximum values of each dominant modal parameter. At each spindle speedn, a smaller axial cutting depth ap is first defined. With the given (n, ap) pair, the vertices and stabilities of the edges at each frequency within the defined variation range are determined according to equations (9) to (11) and the zero exclusion condition described in Table 1. If the formed polygons are stable for all frequencies, the axial cutting depth is increased by a small increment to repeat the aforementioned calculation until an unstable condition occurs. Then, for this spindle speed value, the corresponding axial cutting depth is the limiting one and the frequency is the chatter frequency. After the repeated calculations sweep all the defined spindle speed values, a robust SLD can be plotted by the detected limiting (ap, n) pairs.


## 2.1. DOF Milling Stability Model

A two-degree-of-freedom (DOF) dynamic model shown in Figure1(a) is generally used to describe the milling process system, in which the tool is assumed to have two orthogonal degrees of freedom in X and Y directions, and m, ξ, and k represent corresponding mass, damping, and stiffness coefficients. According to the analytical milling stability theory presented by Budak and Altintas [13, 14], dynamic variations of the cutting forces in two directions ΔFx and ΔFy can be can be described as follows:(1)ΔFxtΔFyt=apKtcN4πaxxaxyaxyayyΔxΔy=apKtcN4πAΔxΔy,where ap is the axial cutting depth, N is the tool tooth number, Ktc is the tangential cutting force coefficient, Δx and Δy are the X- and Y-directional vibration displacements, and [A] is the average directional factor which is expressed as follows:(2)axx=−cos2φjl+2Krtφjl−Krtsin2φjlφstφex,axy=2φjl+sin2φjl−Krtcos2φjlφstφex,ayx=−2φjl+sin2φjl−Krtcos2φjlφstφex,ayy=cos2φjl+2Krtφjl+Krtsin2φjlφstφex,where Krt is the ratio of radical cutting force coefficient to the tangential one and φst and φex are the start and exit angles of the cutting tooth.Figure 1
(a) The 2DOF dynamic milling system. (b) The traditional stability lobe diagram.
(a)(b)In frequency domain, vibration displacements Δx and Δy can be described based on the transfer function matrix G(iω):(3)ΔxΔy=1−e−iωτGiωFeiωt⟵Giω=GxxiωGxyiωGyxiωGyyiω,where τ is the tool tooth passing period, Gxx(iω) and Gyy(iω) are the direct tool point FRFs in X and Y directions, and Gyx(iω) and Gyx(iω) are the cross tool point FRFs in X and Y directions. Appling equation (3) into equation (1), the dynamic cutting forces are expressed using the following equation:(4)Feiωt=apKtcN4πA1−e−iωτGiωFeiωt.An eigenvalue analysis is performed on equation (4) to obtain its following characteristic equation:(5)detI−apKtcN4π1−e−iωτAGiω=0.Then, if the cross tool point FRFs in two directions are ignored, the characteristic equation can be rewritten as follows:(6)a0Λ2+a1Λ+1=0Λ=−apKtcN4π1−e−iωτ=ΛI+ΛR⟵a0=GxxiωGyyiωaxxayy−axyayxa1=axxGxxiω+ayyGyyiω,where Λ is the eigenvalue and its real and imaginary parts are named ΛR and ΛI. The critical stability point occurs at the chatter frequency ωc, and it is substituted into equation (6) to obtain the corresponding limiting axial cutting depth aplim and spindle speed n:(7)aplim=−4πΛR+iΛIKtcNcosωcτ−isinωcτ,n=60ωcN2k+1π−2arctanΛR/ΛI.Therefore, in the spindle speed range, a traditional stability lobe diagram shown in Figure1(b) can be plotted according to equation (7). The lobe curve is formed by different combinations of the limiting axial cutting depth aplim and spindle speed n, which is used to divide the stable and unstable zooms. The (ap, n) pairs below the curve correspond to stable machining conditions, and the (ap, n) pairs above the curve result in unstable machining conditions.


## 2.2. Robust Milling Stability Prediction by Edge Theorem and Zero Exclusion Condition

Equations (1) to (7) show that the SLD mainly depends on the tool point FRFs Gxx(iω) and Gyy(iω). A general FRF G(iω) with Nm modes can be described using the following equation [25]:(8)Giω=∑r=1Nmωnr/kers2+2ξrωnrs+ωnr2,where s = iω is the Laplace variable and ωnr, ξer, and ker are the rth modal natural frequency, damping ratio, and stiffness, respectively. During a milling operation, the machine tool structure and its joint contact positions vary with movements of the moving parts, causing variations in tool point FRFs and related modal parameters and then leading to uncertain SLD prediction.Considering that the modal parameters will vary within the machining space determined by the workpiece machining allowances, the edge theorem and zero exclusion condition are applied to develop a robust milling stability prediction by extending traditional milling chatter stability theory described in Section2.1.The edge theorem can be used to determine the stability of the time-delay system with uncertainties varying within a specific range defined by the minimum and maximum values. According to the edge theorem, for a polynomialP with uncertain parameters in the Laplace domain, the minimum and maximum values of each uncertain parameter can be selected to form different combinations, and a corresponding family of polynomials can be obtained by substituting these combinations into the polynomial. Calculating each extreme polynomial at a given frequency will obtain different vertices on the complex plane, and the adjacent vertices can be connected to form a polygon. If the polygon edges are all stable, robust stability of the system within the boundary of edges can be guaranteed.For applying the edge theorem to the 2DOF milling system with uncertain modal parameters, the characteristic equation of the chatter stability model expressed in equation (6) is first rewritten in the Laplace domain as follows [22–24]:(9)0=Gxxs,ωnGyys,ωnaxxayy−axyayx−NapKtc4π1−e−sτ2+axxGxxs,ωn+ayyGyys,ωn×−NapKtc4π1−e−sτ+1.Since the application of edge theorem in performing the robust stability evaluation requires the system characteristic equation to be in a polynomial form, the denominator of one dominant mode shown in equation (8) is multiplied to the left and right sides of equation (9) to yield equation (10) in a polynomial form:(10)0=Gxxs,ωnGyys,ωnaxxayy−axyayx−NapKtc4π1−e−sτ2+axxGxxs,ωn+ayyGyys,ωn−NapKtc4π1−e−sτ+1×s2+2ξiωn,is+ωn,i2.The variations of dominant modal parameters are taken to represent the machining position changes, and they are assumed to vary within its minimum and maximum values. Then, the left and right boundary values of the natural frequency, damping ratio, and stiffness variation internals inX and Y directions are defined as ωnxmin, ωnxmax, ωnymin, ωnymax, ξxmin, ξxmax, ξymin, ξymax, kxmin, kxmax, kymin, and kymax. These extreme values of modal parameters are selected to form different combinations, which are substituted into equation (10) to determine a family of polynomials P as expressed in the following equation:(11)P=p1s,v,p2s,v,…,pis,vv=ωnx,ξx,kx,ωny,ξy,ky⟸ωn∈ωnmin,ωnmax,ξ∈ξmin,ξmax,k∈kmin,kmax,where i is the total number of the combinations with extreme values and i will be 2Nv if the number of varying parameters is Nv. As one dominant mode has three modal parameters, there are 6 uncertain parameters considering X and Y directions. Thus, Nv equals 6, and there are 64 combinations containing extreme values of modal parameters. Then, values of the extreme combinations are calculated using equation (10) to produce vertices in the complex plane when a frequency is given. The adjacent vertices are connected to form the edges of a polygon, and the robust stability of the milling system is dependent on the edges according to the edge theorem. Therefore, the zero exclusion condition is introduced to efficiently evaluate the stability of each edge.The zero exclusion condition is a graphic approach to investigate the relationships among the edge locations and the origin at a given frequency in the complex plane. The system is stable if the origin is located inside the polygon formed by the edges; otherwise, the system is unstable. A system with two uncertain parameters is assumed to have a better illustration of the zero exclusion condition. The two uncertain parameters yield four vertices, which form a quadrilateral in the complex plane as shown in Figure2. The quadrant angle α of each vertex is defined, and distributions of the quadrant angles are used to discuss whether the quadrilateral encircles the origin. Five cases are summarized in Table 1 according to distributions of the quadrant angles. The symbols + and ‒ shown in Figure 2 and Table 1 mean the maximum and minimum quadrant angles in one quadrant, respectively.Figure 2
The quadrant angles in the complex plane.Table 1
Stable or unstable conditions for different distributions of the vertices in the complex plane.
Case no.Milling conditionDistributions of the vertices in the complex plane1StableAll vertices are located in the same quadrant, and the origin is outside the polygon.2StableVertices are located in two adjacent quadrants, such as I-II, II-III, III-IV, or IV-I.3StableVertices are located in diagonal quadrants such as I–III or II–IV, and the system is stable only whenαm+2+-αm_ < π or αm+2+-αm+ > π. m is the mth quadrant, and αm+2+ and αm_ are the minimum and maximum quadrant angles.4StableVertices are located in three quadrants as I-II-III, II-III-IV, III-IV-I, and IV-I-II where only the following four conditions correspond to the stable system:α3+-α1_ < π in I-II-III; α4+-α2_ < π in II-III-IV; α3_-α1+ < π in III-IV-I; α4_-α2+ < π in IV-I-II.5UnstableVertices are located in four quadrants, and the origin is inside the polygon.Therefore, the robust milling stability prediction is transformed into a graphic problem. Modal parameters corresponding to the positions within the machining space are predicted to determine the minimum and maximum values of each dominant modal parameter. At each spindle speedn, a smaller axial cutting depth ap is first defined. With the given (n, ap) pair, the vertices and stabilities of the edges at each frequency within the defined variation range are determined according to equations (9) to (11) and the zero exclusion condition described in Table 1. If the formed polygons are stable for all frequencies, the axial cutting depth is increased by a small increment to repeat the aforementioned calculation until an unstable condition occurs. Then, for this spindle speed value, the corresponding axial cutting depth is the limiting one and the frequency is the chatter frequency. After the repeated calculations sweep all the defined spindle speed values, a robust SLD can be plotted by the detected limiting (ap, n) pairs.


## 3. The RBFANN-Based Robust Position-Dependent Milling Chatter Stability

As the tool point FRF is a function of the displacements of the machine tool moving parts, the position dependency should be considered to predict milling stability more accurately. It is time consuming to repeatedly perform the impact testing or virtual simulation to obtain the tool point FRFs at any position within the whole machine tool work volume. Moreover, the machining space is also uncertain as the workpiece machining allowances vary in real milling operations, bringing many difficulties in determining the variation ranges of the modal parameters for applying the edge theorem in predicting the robust milling chatter stability. Therefore, the RBFANN is introduced in this paper to predict position-dependent modal parameters of tool point FRFs and their variation ranges corresponding to machining allowances.


### 3.1. A Generalized RBFANN Topologic Structure

RBFANN is an efficient three-layer forward neural network with the activation function of radial basis function, which can approximate to a nonlinear function accurately with simple topological structure and fast convergence speed. A generalized RBFANN structure is shown in Figure3, which consists of one input layer, one hidden layer, and one output layer [26–28]. The input layer contains the input variable vector X = [X1, X2,…,Xm]T, where its dimension m is equal to the number of neurons in this layer. The input vector is transferred into the nonlinear hidden layer by the radial basis functions, changing the original nonlinear problem into a separable linear problem. The neural number in the hidden layer is dependent on the system complexity, and it is no more than dimension of the input variable vector. For the ith neuron in the hidden layer, its output can be written as Φ(‖X−Xic‖), where Xic = [Xic1, Xic2,…,Xicm]T is the center of the basis function and ‖·‖ means the 2-nom. The Gauss function is generally the basis function of RBFANN, and then the output of the ith hidden layer neuron corresponding the kth input sample vector Xk can be described as follows:(12)ΦXk,Xic=exp−12σ2Xk−Xic2,where σ is the standard deviation of the Gauss function and can be determined according to equation (13) to avoid the Gauss curve to be too sharp or flat.(13)σ=dmax2Numh,where dmax is the maximum distance between the defined centers and Numh is the neuron number of the nonlinear hidden layer.Figure 3
The topologic structure of a generalized RBFANN.Neurons in the hidden layer are linked to the neurons in the output layer by the weights shown in Figure3. Then, the output layer is a linear weighting superposition of the output of each hidden layer neuron. Assuming that the output response vector Yki = [yk1, yk2,…,ykq]T has q dimensions, the neuron number in the output layer is equal to q. Therefore, the output of one neuron in the output layer can be described as follows:(14)ykj=w0j+∑i=1NumhwijΦXk,Xic=∑i=0NumhwijΦXk,Xic,j=1,2,…,q,where k is the kth input training sample, j is the jth output layer neuron, i is the ith hidden layer neuron, and wij is the weight between the ith and jth neurons in the hidden and output layers, respectively.


### 3.2. Learning Algorithm of a RBFANN

Equations (12) to (14) show that the basic function center Xc, the standard deviation of the Gauss function σ, and the weights linking the hidden and output layers are three important parameters to establish an accurate RBFANN. The self-organizing selection center method is used in this paper to determine these parameters with the following two procedures. First, the center Xc and the standard deviation σ are obtained through the unsupervised learning process, and then the linking weights w are optimized by a supervised learning process.The K-means clustering algorithm is one of the most common used clustering methods, which divides the data into different categories with similar internal properties to get representative centers. Three stages are summarized to illustrate the application of the K-means clustering method as follows [29]:(1)
Randomly selectL different vectors from the input training sample information to be the initial cluster centers Xc.(2)
Randomly select one vectorXk from the input training sample information and calculate the distance between it and each cluster center using equation (15). If the minimum distance is between Xk and the ith cluster Xci(s) after s iterations, Xk is classified to the ith cluster center, and then the ith cluster center is updated according to equation (16) as a new vector Xk joins in:(15)iXk=argminXk−Xcis,(16)Xcis+1=Xcis+λXk−Xcis,i=iXk,Xcis,i≠iXk,whereλ is the learning rate between 0 and 1.(3)
If the variations of the cluster centers are less than the defined small threshold, the K-means algorithm is regarded to have a convergence. Otherwise, the calculation will jump into the second stage to continue the iteration.Therefore,Xci(s) at the end of the iteration is the final cluster center. Then, after the cluster centers are determined, the standard deviation σ can be calculated using equation (13). On this basis, the linking weights between the hidden and output layer can be optimized by the following cost function:(17)Ej=12∑k=1Numtejk2,where Ej is the error of the jth neuron in the output layer, Numt is the total number of the training sample, and ejk is the error between the predicted and actual output values of the jth output layer neuron for the kth training sample.(18)ejk=yrjk−∑i=0Numh−1wijΦXk,Xic,where yrjk is the actual value of the jth output layer neuron for the kth training sample and Numh is the total neuron number of the hidden layer.Then, the optimal linking weights are searched for minimizing the cross function with the gradient descent algorithm.(19)∂Ev∂wijv=∑k=1NumteΦXk−Xic⟶wijv+1=wijv−η∂Ev∂wijv,where η is the learning rate with the value between 0 and 1.


### 3.3. Application of the RBFANN in Robust Milling Stability Prediction

According to the modal theory, traditional milling chatter stability model, and edge theorem described in Section2, if the position-dependent modal parameters and their variation ranges within the given machining space are obtained, the robust chatter milling stability can be predicted by the edge theorem. Therefore, RBFANNs for predicting two different types of outputs are defined in this paper to benefit the efficient robust milling stability prediction. One type is to predict the position-dependent modal parameters, and the other type is to predict the workpiece machining allowance-dependent variation ranges of modal parameters.The RBFANN for predicting the position-dependent modal parameter can be established according to the five procedures summarized as follows:(1)
Sample Data Processing. Different combinations of the moving parts’ displacements within the whole machine tool work volume are determined with the aid of the experimental design method. The tool point FRFs for each combination are obtained and further used to identify corresponding modal parameters. Then, different displacement combinations are defined as the input sample information, and the corresponding modal parameters are defined as the output sample information.(2)
Design of the Input and Output Layers. The obtained sample information in procedure (1) is divided into the training and testing sample, respectively, and the sample information is normalized between the interval (−1, 1) to have a fast convergence rate. The normalized input training sample is defined as the input data of the RBFANN, and the normalized output training sample is defined as the output data of the RBFANN. Neuron numbers of the input and output layers are determined by dimensions of the input and output variables, respectively.(3)
Design of the Hidden Layer. The cluster centers and the standard deviation are determined based on the K-means clustering algorithm described in Section 3.2. Thus, the number of neurons in the hidden layer is equal to the number of cluster centers.(4)
Training of the RBFANN. With the obtained input and output training sample information, cluster centers, and standard deviation, the linking weights between the hidden and output layers are optimized according to equations (17) to (19).(5)
Validation of the RBFANN. The trained RBFANN in procedure (4) is used to predict the modal parameters corresponding to the input testing sample, which are compared with those of the original output testing sample to validate the accuracy of the RBFANN.Procedures to establish the RBFANN for predicting the workpiece machining allowance-dependent variation ranges of the modal parameters are similar to the aforementioned five procedures. The differences are mainly lying in the first procedure, where different combinations of positions and machining allowances are defined as the input sample information and the obtained minimum and maximum values of the modal parameters within the corresponding machining spaces are defined as the output sample information. Therefore, with the input and output sample information, the second RBFANN can be trained and validated according to the procedures (2) to (5) summarized above.


## 3.1. A Generalized RBFANN Topologic Structure

RBFANN is an efficient three-layer forward neural network with the activation function of radial basis function, which can approximate to a nonlinear function accurately with simple topological structure and fast convergence speed. A generalized RBFANN structure is shown in Figure3, which consists of one input layer, one hidden layer, and one output layer [26–28]. The input layer contains the input variable vector X = [X1, X2,…,Xm]T, where its dimension m is equal to the number of neurons in this layer. The input vector is transferred into the nonlinear hidden layer by the radial basis functions, changing the original nonlinear problem into a separable linear problem. The neural number in the hidden layer is dependent on the system complexity, and it is no more than dimension of the input variable vector. For the ith neuron in the hidden layer, its output can be written as Φ(‖X−Xic‖), where Xic = [Xic1, Xic2,…,Xicm]T is the center of the basis function and ‖·‖ means the 2-nom. The Gauss function is generally the basis function of RBFANN, and then the output of the ith hidden layer neuron corresponding the kth input sample vector Xk can be described as follows:(12)ΦXk,Xic=exp−12σ2Xk−Xic2,where σ is the standard deviation of the Gauss function and can be determined according to equation (13) to avoid the Gauss curve to be too sharp or flat.(13)σ=dmax2Numh,where dmax is the maximum distance between the defined centers and Numh is the neuron number of the nonlinear hidden layer.Figure 3
The topologic structure of a generalized RBFANN.Neurons in the hidden layer are linked to the neurons in the output layer by the weights shown in Figure3. Then, the output layer is a linear weighting superposition of the output of each hidden layer neuron. Assuming that the output response vector Yki = [yk1, yk2,…,ykq]T has q dimensions, the neuron number in the output layer is equal to q. Therefore, the output of one neuron in the output layer can be described as follows:(14)ykj=w0j+∑i=1NumhwijΦXk,Xic=∑i=0NumhwijΦXk,Xic,j=1,2,…,q,where k is the kth input training sample, j is the jth output layer neuron, i is the ith hidden layer neuron, and wij is the weight between the ith and jth neurons in the hidden and output layers, respectively.


## 3.2. Learning Algorithm of a RBFANN

Equations (12) to (14) show that the basic function center Xc, the standard deviation of the Gauss function σ, and the weights linking the hidden and output layers are three important parameters to establish an accurate RBFANN. The self-organizing selection center method is used in this paper to determine these parameters with the following two procedures. First, the center Xc and the standard deviation σ are obtained through the unsupervised learning process, and then the linking weights w are optimized by a supervised learning process.The K-means clustering algorithm is one of the most common used clustering methods, which divides the data into different categories with similar internal properties to get representative centers. Three stages are summarized to illustrate the application of the K-means clustering method as follows [29]:(1)
Randomly selectL different vectors from the input training sample information to be the initial cluster centers Xc.(2)
Randomly select one vectorXk from the input training sample information and calculate the distance between it and each cluster center using equation (15). If the minimum distance is between Xk and the ith cluster Xci(s) after s iterations, Xk is classified to the ith cluster center, and then the ith cluster center is updated according to equation (16) as a new vector Xk joins in:(15)iXk=argminXk−Xcis,(16)Xcis+1=Xcis+λXk−Xcis,i=iXk,Xcis,i≠iXk,whereλ is the learning rate between 0 and 1.(3)
If the variations of the cluster centers are less than the defined small threshold, the K-means algorithm is regarded to have a convergence. Otherwise, the calculation will jump into the second stage to continue the iteration.Therefore,Xci(s) at the end of the iteration is the final cluster center. Then, after the cluster centers are determined, the standard deviation σ can be calculated using equation (13). On this basis, the linking weights between the hidden and output layer can be optimized by the following cost function:(17)Ej=12∑k=1Numtejk2,where Ej is the error of the jth neuron in the output layer, Numt is the total number of the training sample, and ejk is the error between the predicted and actual output values of the jth output layer neuron for the kth training sample.(18)ejk=yrjk−∑i=0Numh−1wijΦXk,Xic,where yrjk is the actual value of the jth output layer neuron for the kth training sample and Numh is the total neuron number of the hidden layer.Then, the optimal linking weights are searched for minimizing the cross function with the gradient descent algorithm.(19)∂Ev∂wijv=∑k=1NumteΦXk−Xic⟶wijv+1=wijv−η∂Ev∂wijv,where η is the learning rate with the value between 0 and 1.


## 3.3. Application of the RBFANN in Robust Milling Stability Prediction

According to the modal theory, traditional milling chatter stability model, and edge theorem described in Section2, if the position-dependent modal parameters and their variation ranges within the given machining space are obtained, the robust chatter milling stability can be predicted by the edge theorem. Therefore, RBFANNs for predicting two different types of outputs are defined in this paper to benefit the efficient robust milling stability prediction. One type is to predict the position-dependent modal parameters, and the other type is to predict the workpiece machining allowance-dependent variation ranges of modal parameters.The RBFANN for predicting the position-dependent modal parameter can be established according to the five procedures summarized as follows:(1)
Sample Data Processing. Different combinations of the moving parts’ displacements within the whole machine tool work volume are determined with the aid of the experimental design method. The tool point FRFs for each combination are obtained and further used to identify corresponding modal parameters. Then, different displacement combinations are defined as the input sample information, and the corresponding modal parameters are defined as the output sample information.(2)
Design of the Input and Output Layers. The obtained sample information in procedure (1) is divided into the training and testing sample, respectively, and the sample information is normalized between the interval (−1, 1) to have a fast convergence rate. The normalized input training sample is defined as the input data of the RBFANN, and the normalized output training sample is defined as the output data of the RBFANN. Neuron numbers of the input and output layers are determined by dimensions of the input and output variables, respectively.(3)
Design of the Hidden Layer. The cluster centers and the standard deviation are determined based on the K-means clustering algorithm described in Section 3.2. Thus, the number of neurons in the hidden layer is equal to the number of cluster centers.(4)
Training of the RBFANN. With the obtained input and output training sample information, cluster centers, and standard deviation, the linking weights between the hidden and output layers are optimized according to equations (17) to (19).(5)
Validation of the RBFANN. The trained RBFANN in procedure (4) is used to predict the modal parameters corresponding to the input testing sample, which are compared with those of the original output testing sample to validate the accuracy of the RBFANN.Procedures to establish the RBFANN for predicting the workpiece machining allowance-dependent variation ranges of the modal parameters are similar to the aforementioned five procedures. The differences are mainly lying in the first procedure, where different combinations of positions and machining allowances are defined as the input sample information and the obtained minimum and maximum values of the modal parameters within the corresponding machining spaces are defined as the output sample information. Therefore, with the input and output sample information, the second RBFANN can be trained and validated according to the procedures (2) to (5) summarized above.


## 4. Case Study

The proposed robust milling chatter stability prediction method based on the RBFANN and edge theorem is applied to a real three-axis vertical machining center. The worktable, saddle, and spindle system move along theX, Y, and Z directions, respectively, as shown in Figure 4(a), and the corresponding moving intervals are 0–550 mm, 0–400 mm, and 0–400 mm. Then, the X-directional displacement of the worktable, the Y-directional displacement of the saddle, and the Z-directional displacement of the spindle system are combined to represent the machining position.Figure 4
(a) The studied three-axial vertical machining center. (b) Tool point FRFs related to 64 positions.
(a)(b)


### 4.1. Sample Information Based on Experiment Design Method

Extreme values of thex, y, and z directions enclose the machine tool work volume. To efficiently obtain the tool point FRFs within the machine tool work volume, the orthogonal experimental design method is introduced to determine some representative machining positions. According to variation ranges of the x, y, and z directions, 8 specific values were selected uniformly for each directional displacement [30]. Factors and their levels are listed in Table 2; 64 combinations of the displacements in three directions are arranged in Table 3.Table 2
Factors and levels for the orthogonal experiment design method.
FactorLevel12345678x (mm)70140210280350420490550y (mm)50100150200250300350400z (mm)50100150200250300350400Table 3
Schemes arranged by the orthogonal experimental table.
Orderx (mm)y (mm)z (mm)Orderx (mm)y (mm)z (mm)Orderx (mm)y (mm)z (mm)17462354845368238624863461443717255264776446522677148845533127458495156441282655047675732988151256853730625527359432315625332410278321115463811467331225567412221345515672813414353425715514661368275858415872376165918816313382876021217854394856168318423401336237519753411666317720357422346478221647438362224344818The moving parts were driven to the 64 machining positions with the arranged displacements in sequence, and the tool point FRFs inX and Y directions at each position were obtained by impact testing shown in Figure 4(a) [31]. The 64 FRF curves for each direction are described in Figure 4(b), where the natural frequencies and amplitudes show obvious differences. There are four obvious modes of each X-directional FRF and three obvious modes of each Y-directional tool point FRF. Then, corresponding modal parameters of each tool point FRF were identified according to the modal theory. Two different positions are randomly selected from the 64 combinations, and their corresponding modal parameters are identified in Table 4. Difference of the natural frequency, damping ratio, and stiffness values at the two positions is observed, further indicating that the tool point FRFs are dependent on the machining position.Table 4
Modal parameters of two different positions.
DirectionModesPosition 1Position 2Natural frequencyω (Hz)Damping ratioζ (%)StiffnessK (×109 N/m)Natural frequencyω (Hz)Damping ratioζ (%)StiffnessK (×109 N/m)X15041.980.884845.370.6629121.976.589163.284.75312281.791.1312523.995.47413923.000.3914004.600.36Y13842.081.493804.213.7426361.890.136245.060.18314321.681.1414084.831.72


### 4.2. RBFANN for Position-Dependent Modal Parameter Prediction

Since each mode has three parameters including the natural frequency, modal damping ratio, and modal stiffness, there are total 21 modal parameters for the four and three obvious modes of theX- and Y-directional tool point FRFs. Predicting the 21 modal parameters at one time will introduce difficulties in training the RBFANN with higher accuracy and faster convergence speed. Therefore, RBFANNs are utilized to predict the natural frequencies, modal damping ratios, and modal stiffnesses, respectively. Then, 3 RBFANNs are trained to predict four natural frequencies, four modal damping ratios, and four modal stiffnesses in X direction, and another 3 RBFANNs are trained to predict three natural frequencies, three modal damping ratios, and three modal stiffnesses in Y direction. The training process of a single RBFANN for predicting the modal stiffnesses in X direction is described to give an illustration.The input layer has 3 neurons standing for the displacementsX, Y, and Z in three directions, and the output layer has 4 neurons standing for the modal stiffness of each mode. According to Section 4.1, 58 combinations of the position information were randomly selected from Table 2 to be the input training sample, and their corresponding modal stiffnesses were taken as the output training samples. The remaining 6 combinations of the position information and the corresponding modal stiffnesses were defined as the input and output testing samples, respectively. Then, the RBFANN for predicting the position-dependent four modal stiffnesses were trained based on the 58 input and output training samples, where the neuron number in the hidden layer is 55 and the standard deviation σ is 0.2. For the six testing samples, errors between the predicted and actual modal stiffness values are described in Figure 5(a). According to the absolute error distributions in Figure 5(b), the absolute errors all below 1.4% show that the trained RBFANN has a higher accuracy of predicting the position-dependent modal stiffness. Moreover, the varying fourth-order modal stiffness values within the Y-Z plane with the X-directional displacement of 300 mm are described in Figure 5(c) to show the necessity of considering the position effects.Figure 5
(a) Errors of the six testing samples . (b) Absolute error distributions of the six testing samples. (c) Variations of the modal stiffness within theY-Z plane.
(a)(b)(c)


### 4.3. RBFANN for Extreme Value Prediction of the Dominant Modal Parameters

Section4.2 shows that the modal parameters vary within the machine tool work volume. Then, the predicted chatter-free machining parameters with the tool point FRFs at one position cannot guarantee the stability during real milling process as the relative position of the tool and workpiece changes. Thus, the edge theorem and exclusion condition described in Section 2 are used to perform the robust milling stability prediction. The application of the edge theorem needs to obtain the maximum and minimum values of each modal parameter within the machining space. The workpiece and its machining allowances in X, Y, and Z directions are described in Figure 6(a), where the three-directional machining allowances form the machining space. At one position, the machining space will vary with the three-directional machining allowances, affecting the corresponding extreme values of the modal parameters. Therefore, the RBFANN is again used to predict the position and machining allowance-dependent extreme values of the modal parameters.Figure 6
(a) Machining space defined by the machining allowances. (b) An illustration of the sample position and its related machining space. (c) Absolute error distributions of the extreme modal parameters.
(a)(b)(c)According to the edge theorem in Section2.2, there will be 221 polynomials determined to check the system stability if all 21 uncertain modal parameters are considered, seriously increasing difficulties of the robust milling chatter stability calculation. Research studies have shown that the dominant modes of the tool point FRFs in X and Y directions exert more effects on the milling stability than other modes with smaller amplitudes [5]. Therefore, the 4th and 2nd modes with higher amplitudes in X and Y directions, respectively, are determined as the dominant modes, and their corresponding 6 modal parameters are defined as the uncertain parameters of the robust prediction.Then, the defined 64 positions in Section4.1 were used to determine some related machining spaces for obtaining the sample information of the RBFANN. A position described in Figure 6(b) is taken as an instance to illustrate how to obtain its related sample information. Coordinates of the example position are defined as xs, ys, and zs, and limiting coordinates of three directions are defined as xlim, ylim, and zlim. Then, sample coordinates [xs, x1, …, xo, …, xlim], [ys, y1, …,yp, …, ylim], and [zs, z1, …, zq, …, zlim] of each direction are determined from the intervals [xpxlim], [ypylim], and [zpzlim], and the adjacent coordinates of each direction have a 40 mm difference. Then, the distance between xs and xo, the distance between ys and yp, and the distance between zs and zq are taken as the length, width, and height of one machining space. Within each machining space, its corresponding machining allowance of each direction is divided uniformly to define some subpositions, and adjacent coordinates of each direction have a 5 mm difference. Then, dominant modal parameters of all subpositions are predicted by the trained RBFANNs described in Section 4.1 to determine the minimum and maximum values of the 6 dominant modal parameters. After the aforementioned procedures have been performed on the 64 sample positions in Table 3 in sequence, 8816 machining spaces are been obtained and defined as the sample information for training the RBFANNs.As each dominant modal parameter has two extreme values, 12 extreme values should be predicted for 6 dominant modal parameters in x and y directions. Using one RBFANN to predict 12 extreme values simultaneously will slow the convergence, so one RBFANN was established to predict 6 extreme values forx direction and the other one was established to predict 6 extreme values in y direction. Then, for one RBFANN, the input layer has 6 neurons including the position (X, Y, and Z) and the machining allowances (L, W, and H), and the output layer has 6 neurons including the minimum and maximum values of each dominant modal parameter. Taking the RBFANN for predicting the X-directional dominant modal parameters as an instance, 8700 combinations of the position and machining allowance sample information were defined as the input training sample, and their corresponding extreme values of the modal parameters were taken as the output training sample. Then, the RBFANN was trained, and its accuracy was validated by the remaining 116 combinations of the sample information. Distributions of the 116 × 6 = 696 errors between the predicted and actual extreme values of modal parameters are described in Figure 6(c). Most of the errors are less than 1%, and the maximum error is 1.52%, showing the accuracy of the trained RBFANN.


### 4.4. Robust Milling Chatter Stability Prediction with the Edge Theorem and RBFANNs

Considering the machining within the machining space, six modal parameters of the 4th and 2nd dominant modes inX and Y directions, respectively, are taken as the uncertainties, including ωx4th, kx4th, ξx4th, ωy2nd, ky2nd, and ξy2nd. Thus, when the position (X, Y, and Z) and machining allowance (L, W, and H) are defined, the extreme values of the 6 parameters can be calculated based on the two trained RBFANNs described in Section 4.3. The 4th position (420 mm, 250 mm, and 100 mm) in Table 3 and the machining allowances (100 mm, 40 mm, and 20 mm) are taken for instance. The modal parameters at the 4th position and the extreme values of the 6 dominant modal parameters within the machining space are listed in Table 5. Considering the error introduced by the RBFANNs, a ±2.9% variation calculated according to the maximum errors in Sections 4.2 and 4.3 was taken to the basic value of each modal parameter to obtain their extreme values. Then, with the final obtained minimum and maximum values of each modal parameter, the algorithm for the edge theorem and zero exclusion condition mentioned in Section 3 was developed in MATLAB software to predict and plot the robust stability lobe diagram (RSLD) on the basis of a four-teeth tool with a diameter of 20 mm, a radial cutting width of 8 mm, and tangential and radial cutting force coefficients of 1769 MPa and 1219 MPa.Table 5
The minimum and maximum values of the example position and machining space.
DirectionXYNatural frequencyω (Hz)Damping ratioζ (%)StiffnessK (×109 N/m)Natural frequencyω (Hz)Damping ratioζ (%)StiffnessK (×109 N/m)Basic value13922.590.396361.890.13Minimum value1336 × 1.0292.37 × 1.0290.24 × 1.029664 × 1.0291.71 × 1.0290.11 × 1.029Maximum value1468 × 0.9712.81 × 0.9710.43 × 0.971612 × 0.9712.05 × 0.9710.15 × 0.971The obtained RSLD is described as a red lobe curve in Figure7(a), which is formed by different robust pairs of spindle speed and its related maximum axial cutting depth stable for all frequencies within the focus frequency range from 0 Hz to 2000 Hz. The red lobe curve is below the blue lobe curve which is a traditional SLD plotted based on the basic modal parameters listed in Table 5 since it considers the variations of the modal parameters within the machining space and can provide more conservative limiting axial cutting depths. Therefore, the (ap, n) pairs below the red curve is selected to guarantee the milling stability during the real machining process.Figure 7
(a) Comparisons between the RSLD and SLD. (b) An example milling test. (c) The force frequency spectrum for a chatter condition.
(a)(b)(c)Down milling tests have been performed based on the three-axis vertical machining center shown in Figure4(a) to illustrate the feasibility of the RSLD. Eighteen (ap, n) pairs were defined according to the red and blue lobe curves in Figure 7(a), which were used to perform the chatter tests by end milling the 1045 steel with the four-teeth tool and in accordance with the 4th position. During the milling operations, the cutting force signals for each (ap, n) pair were measured as shown in Figure 7(b) and analyzed in the frequency domain to detect whether chatter occurred. A force spectrum of the milling condition with the spindle speed of 3000 rpm and axial cutting depth of 8.2 mm is described in Figure 7(c), and the dominant vibration frequency (718 Hz) is not the tool passing frequency (200 Hz) or its harmonics, indicating that the related (ap, n) pair has caused a chatter vibration [32]. Comparing the distributions of the stable and unstable (ap, n) pairs shown in Figure 7(a), all (ap, n) pairs below the red lobe curve correspond to stable milling conditions, validating that the proposed robust milling chatter prediction method can provide reliable chatter-free machining parameters; all the (ap, n) pairs above the blue curve are corresponding to unstable milling conditions, which should be avoided when selecting chatter-free machining parameters; (ap, n) pairs for stable and unstable milling conditions are both observed within the region between the red and blue curves since each uncertain parameter may fall anywhere within its minimum and maximum values and form a combination of uncertain modal parameters. Therefore, the main benefit of the robust milling chatter prediction is that it can provide conservative ap values to guarantee reliable milling stabilities within the machining space.


## 4.1. Sample Information Based on Experiment Design Method

Extreme values of thex, y, and z directions enclose the machine tool work volume. To efficiently obtain the tool point FRFs within the machine tool work volume, the orthogonal experimental design method is introduced to determine some representative machining positions. According to variation ranges of the x, y, and z directions, 8 specific values were selected uniformly for each directional displacement [30]. Factors and their levels are listed in Table 2; 64 combinations of the displacements in three directions are arranged in Table 3.Table 2
Factors and levels for the orthogonal experiment design method.
FactorLevel12345678x (mm)70140210280350420490550y (mm)50100150200250300350400z (mm)50100150200250300350400Table 3
Schemes arranged by the orthogonal experimental table.
Orderx (mm)y (mm)z (mm)Orderx (mm)y (mm)z (mm)Orderx (mm)y (mm)z (mm)17462354845368238624863461443717255264776446522677148845533127458495156441282655047675732988151256853730625527359432315625332410278321115463811467331225567412221345515672813414353425715514661368275858415872376165918816313382876021217854394856168318423401336237519753411666317720357422346478221647438362224344818The moving parts were driven to the 64 machining positions with the arranged displacements in sequence, and the tool point FRFs inX and Y directions at each position were obtained by impact testing shown in Figure 4(a) [31]. The 64 FRF curves for each direction are described in Figure 4(b), where the natural frequencies and amplitudes show obvious differences. There are four obvious modes of each X-directional FRF and three obvious modes of each Y-directional tool point FRF. Then, corresponding modal parameters of each tool point FRF were identified according to the modal theory. Two different positions are randomly selected from the 64 combinations, and their corresponding modal parameters are identified in Table 4. Difference of the natural frequency, damping ratio, and stiffness values at the two positions is observed, further indicating that the tool point FRFs are dependent on the machining position.Table 4
Modal parameters of two different positions.
DirectionModesPosition 1Position 2Natural frequencyω (Hz)Damping ratioζ (%)StiffnessK (×109 N/m)Natural frequencyω (Hz)Damping ratioζ (%)StiffnessK (×109 N/m)X15041.980.884845.370.6629121.976.589163.284.75312281.791.1312523.995.47413923.000.3914004.600.36Y13842.081.493804.213.7426361.890.136245.060.18314321.681.1414084.831.72


## 4.2. RBFANN for Position-Dependent Modal Parameter Prediction

Since each mode has three parameters including the natural frequency, modal damping ratio, and modal stiffness, there are total 21 modal parameters for the four and three obvious modes of theX- and Y-directional tool point FRFs. Predicting the 21 modal parameters at one time will introduce difficulties in training the RBFANN with higher accuracy and faster convergence speed. Therefore, RBFANNs are utilized to predict the natural frequencies, modal damping ratios, and modal stiffnesses, respectively. Then, 3 RBFANNs are trained to predict four natural frequencies, four modal damping ratios, and four modal stiffnesses in X direction, and another 3 RBFANNs are trained to predict three natural frequencies, three modal damping ratios, and three modal stiffnesses in Y direction. The training process of a single RBFANN for predicting the modal stiffnesses in X direction is described to give an illustration.The input layer has 3 neurons standing for the displacementsX, Y, and Z in three directions, and the output layer has 4 neurons standing for the modal stiffness of each mode. According to Section 4.1, 58 combinations of the position information were randomly selected from Table 2 to be the input training sample, and their corresponding modal stiffnesses were taken as the output training samples. The remaining 6 combinations of the position information and the corresponding modal stiffnesses were defined as the input and output testing samples, respectively. Then, the RBFANN for predicting the position-dependent four modal stiffnesses were trained based on the 58 input and output training samples, where the neuron number in the hidden layer is 55 and the standard deviation σ is 0.2. For the six testing samples, errors between the predicted and actual modal stiffness values are described in Figure 5(a). According to the absolute error distributions in Figure 5(b), the absolute errors all below 1.4% show that the trained RBFANN has a higher accuracy of predicting the position-dependent modal stiffness. Moreover, the varying fourth-order modal stiffness values within the Y-Z plane with the X-directional displacement of 300 mm are described in Figure 5(c) to show the necessity of considering the position effects.Figure 5
(a) Errors of the six testing samples . (b) Absolute error distributions of the six testing samples. (c) Variations of the modal stiffness within theY-Z plane.
(a)(b)(c)


## 4.3. RBFANN for Extreme Value Prediction of the Dominant Modal Parameters

Section4.2 shows that the modal parameters vary within the machine tool work volume. Then, the predicted chatter-free machining parameters with the tool point FRFs at one position cannot guarantee the stability during real milling process as the relative position of the tool and workpiece changes. Thus, the edge theorem and exclusion condition described in Section 2 are used to perform the robust milling stability prediction. The application of the edge theorem needs to obtain the maximum and minimum values of each modal parameter within the machining space. The workpiece and its machining allowances in X, Y, and Z directions are described in Figure 6(a), where the three-directional machining allowances form the machining space. At one position, the machining space will vary with the three-directional machining allowances, affecting the corresponding extreme values of the modal parameters. Therefore, the RBFANN is again used to predict the position and machining allowance-dependent extreme values of the modal parameters.Figure 6
(a) Machining space defined by the machining allowances. (b) An illustration of the sample position and its related machining space. (c) Absolute error distributions of the extreme modal parameters.
(a)(b)(c)According to the edge theorem in Section2.2, there will be 221 polynomials determined to check the system stability if all 21 uncertain modal parameters are considered, seriously increasing difficulties of the robust milling chatter stability calculation. Research studies have shown that the dominant modes of the tool point FRFs in X and Y directions exert more effects on the milling stability than other modes with smaller amplitudes [5]. Therefore, the 4th and 2nd modes with higher amplitudes in X and Y directions, respectively, are determined as the dominant modes, and their corresponding 6 modal parameters are defined as the uncertain parameters of the robust prediction.Then, the defined 64 positions in Section4.1 were used to determine some related machining spaces for obtaining the sample information of the RBFANN. A position described in Figure 6(b) is taken as an instance to illustrate how to obtain its related sample information. Coordinates of the example position are defined as xs, ys, and zs, and limiting coordinates of three directions are defined as xlim, ylim, and zlim. Then, sample coordinates [xs, x1, …, xo, …, xlim], [ys, y1, …,yp, …, ylim], and [zs, z1, …, zq, …, zlim] of each direction are determined from the intervals [xpxlim], [ypylim], and [zpzlim], and the adjacent coordinates of each direction have a 40 mm difference. Then, the distance between xs and xo, the distance between ys and yp, and the distance between zs and zq are taken as the length, width, and height of one machining space. Within each machining space, its corresponding machining allowance of each direction is divided uniformly to define some subpositions, and adjacent coordinates of each direction have a 5 mm difference. Then, dominant modal parameters of all subpositions are predicted by the trained RBFANNs described in Section 4.1 to determine the minimum and maximum values of the 6 dominant modal parameters. After the aforementioned procedures have been performed on the 64 sample positions in Table 3 in sequence, 8816 machining spaces are been obtained and defined as the sample information for training the RBFANNs.As each dominant modal parameter has two extreme values, 12 extreme values should be predicted for 6 dominant modal parameters in x and y directions. Using one RBFANN to predict 12 extreme values simultaneously will slow the convergence, so one RBFANN was established to predict 6 extreme values forx direction and the other one was established to predict 6 extreme values in y direction. Then, for one RBFANN, the input layer has 6 neurons including the position (X, Y, and Z) and the machining allowances (L, W, and H), and the output layer has 6 neurons including the minimum and maximum values of each dominant modal parameter. Taking the RBFANN for predicting the X-directional dominant modal parameters as an instance, 8700 combinations of the position and machining allowance sample information were defined as the input training sample, and their corresponding extreme values of the modal parameters were taken as the output training sample. Then, the RBFANN was trained, and its accuracy was validated by the remaining 116 combinations of the sample information. Distributions of the 116 × 6 = 696 errors between the predicted and actual extreme values of modal parameters are described in Figure 6(c). Most of the errors are less than 1%, and the maximum error is 1.52%, showing the accuracy of the trained RBFANN.


## 4.4. Robust Milling Chatter Stability Prediction with the Edge Theorem and RBFANNs

Considering the machining within the machining space, six modal parameters of the 4th and 2nd dominant modes inX and Y directions, respectively, are taken as the uncertainties, including ωx4th, kx4th, ξx4th, ωy2nd, ky2nd, and ξy2nd. Thus, when the position (X, Y, and Z) and machining allowance (L, W, and H) are defined, the extreme values of the 6 parameters can be calculated based on the two trained RBFANNs described in Section 4.3. The 4th position (420 mm, 250 mm, and 100 mm) in Table 3 and the machining allowances (100 mm, 40 mm, and 20 mm) are taken for instance. The modal parameters at the 4th position and the extreme values of the 6 dominant modal parameters within the machining space are listed in Table 5. Considering the error introduced by the RBFANNs, a ±2.9% variation calculated according to the maximum errors in Sections 4.2 and 4.3 was taken to the basic value of each modal parameter to obtain their extreme values. Then, with the final obtained minimum and maximum values of each modal parameter, the algorithm for the edge theorem and zero exclusion condition mentioned in Section 3 was developed in MATLAB software to predict and plot the robust stability lobe diagram (RSLD) on the basis of a four-teeth tool with a diameter of 20 mm, a radial cutting width of 8 mm, and tangential and radial cutting force coefficients of 1769 MPa and 1219 MPa.Table 5
The minimum and maximum values of the example position and machining space.
DirectionXYNatural frequencyω (Hz)Damping ratioζ (%)StiffnessK (×109 N/m)Natural frequencyω (Hz)Damping ratioζ (%)StiffnessK (×109 N/m)Basic value13922.590.396361.890.13Minimum value1336 × 1.0292.37 × 1.0290.24 × 1.029664 × 1.0291.71 × 1.0290.11 × 1.029Maximum value1468 × 0.9712.81 × 0.9710.43 × 0.971612 × 0.9712.05 × 0.9710.15 × 0.971The obtained RSLD is described as a red lobe curve in Figure7(a), which is formed by different robust pairs of spindle speed and its related maximum axial cutting depth stable for all frequencies within the focus frequency range from 0 Hz to 2000 Hz. The red lobe curve is below the blue lobe curve which is a traditional SLD plotted based on the basic modal parameters listed in Table 5 since it considers the variations of the modal parameters within the machining space and can provide more conservative limiting axial cutting depths. Therefore, the (ap, n) pairs below the red curve is selected to guarantee the milling stability during the real machining process.Figure 7
(a) Comparisons between the RSLD and SLD. (b) An example milling test. (c) The force frequency spectrum for a chatter condition.
(a)(b)(c)Down milling tests have been performed based on the three-axis vertical machining center shown in Figure4(a) to illustrate the feasibility of the RSLD. Eighteen (ap, n) pairs were defined according to the red and blue lobe curves in Figure 7(a), which were used to perform the chatter tests by end milling the 1045 steel with the four-teeth tool and in accordance with the 4th position. During the milling operations, the cutting force signals for each (ap, n) pair were measured as shown in Figure 7(b) and analyzed in the frequency domain to detect whether chatter occurred. A force spectrum of the milling condition with the spindle speed of 3000 rpm and axial cutting depth of 8.2 mm is described in Figure 7(c), and the dominant vibration frequency (718 Hz) is not the tool passing frequency (200 Hz) or its harmonics, indicating that the related (ap, n) pair has caused a chatter vibration [32]. Comparing the distributions of the stable and unstable (ap, n) pairs shown in Figure 7(a), all (ap, n) pairs below the red lobe curve correspond to stable milling conditions, validating that the proposed robust milling chatter prediction method can provide reliable chatter-free machining parameters; all the (ap, n) pairs above the blue curve are corresponding to unstable milling conditions, which should be avoided when selecting chatter-free machining parameters; (ap, n) pairs for stable and unstable milling conditions are both observed within the region between the red and blue curves since each uncertain parameter may fall anywhere within its minimum and maximum values and form a combination of uncertain modal parameters. Therefore, the main benefit of the robust milling chatter prediction is that it can provide conservative ap values to guarantee reliable milling stabilities within the machining space.


## 5. Conclusions

This paper considers the variations of tool point FRFs within the machining space and then provides a method to benefit the rapid prediction of the robust milling chatter stability with the edge theorem and radial basis function artificial neural networks. According to the modal theory, variations of the tool point FRFs can be represented by variations of the modal parameters. Then, the modal parameters are taken as uncertain parameters, and the edge theorem is introduced to predict robust pairs of spindle speed and limiting axial cutting depth in a graphic approach based on the traditional milling chatter stability model and zero exclusion condition. As the edge theorem needs the extreme values of each modal parameter in advance, RBFANNs for predicting the position-dependent modal parameters inX and Y directions are first presented and further used to obtain the minimum and maximum values of the dominant modal parameters within the machining space defined by the machining allowances in X, Y, and Z directions. Therefore, with the obtained extreme values of each modal parameter, the robust pairs of spindle speed and limiting axial cutting depth can be predicted and utilized to plot the robust stability lobe diagram.A case study was performed on a real vertical three-axial vertical machining center to illustrate the feasibility of the proposed method. The orthogonal experiment design method is used to determine 64 sample positions within the whole machine tool work volume defined by extreme moving distances of the saddle, worktable, and spindle system inX, Y, and Z directions. Impact testing was performed to obtain the tool point FRFs and corresponding modal parameters of the 64 positions. Then, 6 RBFANNs were trained to predict the natural frequencies, modal stiffnesses, and modal damping ratios in X and Y directions, respectively. Moreover, 8816 machining spaces were determined according to the 64 sample positions and the defined machining allowances, and then the minimum and maximum values of the dominant modal parameters within each machining space were obtained with the aid of the 6 RBFANNs established above. On this basis, two RBFANNs for predicting the position and machining allowance-dependent extreme values of each dominant modal parameter in X and Y directions were established, respectively. The robust milling chatter stability prediction was developed based on the 4th sample position, and the RSLD was plotted and compared with the traditional SLD. The RSLD is below the SLD and provides more conservative limiting axial cutting depths. Several (ap, n) pairs were determined according to the two lobe curves to perform the milling tests on the machine tool, and the milling cutting forces were measured to detect the chatter occurrences. The (ap, n) pairs below the RSLD all correspond to stable conditions, and the stable and unstable (ap, n) pairs are both observed between two lobe curves, validating the feasibility of the RSLD for predicting reliable chatter-free machining parameters. Therefore, the proposed method for predicting the robust milling chatter stability can benefit the reliable selection of chatter-free machining parameters considering uncertain machining position information at the design stage. In further research, the proposed method can be extended to take more uncertain parameters into consideration such as the temperature and cutting forces, and much more accurate models for predicting the position-dependent modal parameters and their extreme values can be studied to reduce the effects caused by the prediction errors.



---

*Source: 1398957-2020-10-12.xml*