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

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 deﬁned by the machining position and machining allowance information. Then, radial basis function artiﬁcial neural networks (RBFANNs) are used to predict the position-dependent modal parameters in X 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 ( a p lim ) 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 ( a p , n ) pairs to guarantee the milling stability, validating the feasibility of the proposed robust milling chatter stability prediction method.


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][2][3]. e regenerative chatter is caused by the forces generated during the dynamic cutting process and is not dependent on the external periodic forces. e 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]. erefore, 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][7][8].
e 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 timedelay term simultaneously and proposed an updated fulldiscretization 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. en, the position-dependent tool point FRFs result in uncertain prediction of milling stability. erefore, the positiondependent 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.
en, a state-space model was obtained and incorporated into a time domain simulation of the positiondependent dynamic behaviour. e 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. is 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.
en, at the optimal position, modal parameter variations within the machining space defined by threedirectional 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. erefore, 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. en, a robust position-dependent SLD for selecting robust chatterfree 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 positiondependent milling stability based on the edge theorem and zero exclusion condition method is first presented in Section 2.
e 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. is is followed by a summarized conclusion in Section 5.

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. en, 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.  Figure 1(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 ΔF x and ΔF y can be can be described as follows: a xx a xy a xy a yy where a p is the axial cutting depth, N is the tool tooth number, K tc 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: where K rt 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. In frequency domain, vibration displacements Δx and Δy can be described based on the transfer function matrix G(iω): where τ is the tool tooth passing period, G xx (iω) and G yy (iω) are the direct tool point FRFs in X and Y directions, and G yx (iω) and G yx (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: An eigenvalue analysis is performed on equation (4) to obtain its following characteristic equation: en, if the cross tool point FRFs in two directions are ignored, the characteristic equation can be rewritten as follows: where Λ is the eigenvalue and its real and imaginary parts are named Λ R and Λ I . e critical stability point occurs at the chatter frequency ω c , and it is substituted into equation (6) to obtain the corresponding limiting axial cutting depth a plim and spindle speed n: erefore, in the spindle speed range, a traditional stability lobe diagram shown in Figure 1(b) can be plotted according to equation (7). e lobe curve is formed by different combinations of the limiting axial cutting depth a plim and spindle speed n, which is used to divide the stable and unstable zooms. e (a p , n) pairs below the curve correspond to stable machining conditions, and the (a p , n) pairs above the curve result in unstable machining conditions.

Robust Milling Stability Prediction by Edge eorem and
Zero Exclusion Condition. Equations (1) to (7) show that the SLD mainly depends on the tool point FRFs G xx (iω) and G yy (iω). A general FRF G(iω) with N m modes can be described using the following equation [25]: where s � iω is the Laplace variable and ω nr , ξ er , and k er 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 Section 2.1. e 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 polynomial P 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][23][24]: 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: e 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. en, the left and right boundary values of the natural frequency, damping ratio, and stiffness variation internals in X and Y directions are defined as ω nxmin , ω nxmax , ω nymin , ω nymax , ξ xmin , ξ xmax , ξ ymin , ξ ymax , k xmin , k xmax , k ymin , and k ymax .
ese 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: where i is the total number of the combinations with extreme values and i will be 2 Nv if the number of varying parameters is N v . As one dominant mode has three modal parameters, there are 6 uncertain parameters considering X and Y directions. us, N v equals 6, and there are 64 combinations containing extreme values of modal parameters. en, values of the extreme combinations are calculated using equation (10) to produce vertices in the complex plane when a frequency is given. e 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. erefore, the zero exclusion condition is introduced to efficiently evaluate the stability of each edge. e 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. e 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. e two uncertain parameters yield four vertices, which form a quadrilateral in the complex plane as shown in Figure 2. e 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. e symbols + andshown in Figure 2 and Table 1 mean the maximum and minimum quadrant angles in one quadrant, respectively. erefore, 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 speed n, a smaller axial cutting depth a p is first defined. With the given (n, a p ) 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. en, 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 (a p , n) pairs.

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. erefore, 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.

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 Figure 3, which consists of one input layer, one hidden layer, and one output layer [26][27][28]. e input layer contains the input variable vector X � [X 1 , X 2 ,. . .,X m ] T , where its dimension m is equal to the number of neurons in this layer. e input vector is transferred into the nonlinear hidden layer by the radial basis functions, changing the original nonlinear problem into a separable linear problem. e 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 e 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 X k can be described as follows: 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.
where d max is the maximum distance between the defined centers and N umh is the neuron number of the nonlinear hidden layer. Neurons in the hidden layer are linked to the neurons in the output layer by the weights shown in Figure 3. en, the output layer is a linear weighting superposition of the output of each hidden layer neuron. Assuming that the output response vector Y ki � [y k1 , y k2 ,. . .,y kq ] T has q dimensions, the neuron number in the output layer is equal to q. erefore, the output of one neuron in the output layer can be described as follows: Mathematical Problems in Engineering where k is the kth input training sample, j is the jth output layer neuron, i is the ith hidden layer neuron, and w ij is the weight between the ith and jth neurons in the hidden and output layers, respectively. (12) to (14) show that the basic function center X c , 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. e self-organizing selection center method is used in this paper to determine these parameters with the following two procedures. First, the center X c and the standard deviation σ are obtained through the unsupervised learning process, and then the linking weights w are optimized by a supervised learning process.

Milling condition
Distributions of the vertices in the complex plane 1 Stable All vertices are located in the same quadrant, and the origin is outside the polygon. 2 Stable Vertices are located in two adjacent quadrants, such as I-II, II-III, III-IV, or IV-I.

Stable
Vertices 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.

Stable
Vertices 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: Unstable Vertices are located in four quadrants, and the origin is inside the polygon.

Input Output
Normalization Antinormalization  Figure 3: e topologic structure of a generalized RBFANN. e 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. ree stages are summarized to illustrate the application of the K-means clustering method as follows [29]: (1) Randomly select L different vectors from the input training sample information to be the initial cluster centers X c . (2) Randomly select one vector X k 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 X k and the ith cluster X ci (s) after s iterations, X k is classified to the ith cluster center, and then the ith cluster center is updated according to equation (16) as a new vector X k joins in: 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.
erefore, X ci (s) at the end of the iteration is the final cluster center. en, 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: where E j is the error of the jth neuron in the output layer, N umt is the total number of the training sample, and e jk is the error between the predicted and actual output values of the jth output layer neuron for the kth training sample.
where y rjk is the actual value of the jth output layer neuron for the kth training sample and N umh is the total neuron number of the hidden layer. en, the optimal linking weights are searched for minimizing the cross function with the gradient descent algorithm.
where η is the learning rate with the value between 0 and 1.

Application of the RBFANN in Robust Milling Stability
Prediction. According to the modal theory, traditional milling chatter stability model, and edge theorem described in Section 2, 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. erefore, 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. e 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. e tool point FRFs for each combination are obtained and further used to identify corresponding modal parameters. en, 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. e 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. e 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. e cluster centers and the standard deviation are determined based on the K-means clustering algorithm described in Section 3.2. us, 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. e 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. e 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. erefore, with the input and output sample information, the second RBFANN can be trained and validated according to the procedures (2) to (5) summarized above.

Case Study
e proposed robust milling chatter stability prediction method based on the RBFANN and edge theorem is applied to a real three-axis vertical machining center. e worktable, saddle, and spindle system move along the X, 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. en, 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.

Sample Information Based on Experiment Design Method.
Extreme values of the x, 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. e moving parts were driven to the 64 machining positions with the arranged displacements in sequence, and the tool point FRFs in X and Y directions at each position were obtained by impact testing shown in Figure 4(a) [31]. e 64 FRF curves for each direction are described in Figure 4(b), where the natural frequencies and amplitudes show obvious differences. ere are four obvious modes of each X-directional FRF and three obvious modes of each Ydirectional tool point FRF.
en, 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.

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 the X-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. erefore, RBFANNs are utilized to predict the natural frequencies, modal damping ratios, and modal stiffnesses, respectively. en, 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. e training process of a single RBFANN for predicting the modal stiffnesses in X direction is described to give an illustration.
e input layer has 3 neurons standing for the displacements X, 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. e remaining 6 combinations of the position information and the corresponding modal stiffnesses were defined as the input and output testing samples, respectively. en, 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 positiondependent modal stiffness. Moreover, the varying fourthorder modal stiffness values within the Y-Z plane with the Xdirectional displacement of 300 mm are described in Figure 5(c) to show the necessity of considering the position effects.

RBFANN for Extreme Value Prediction of the Dominant
Modal Parameters. Section 4.2 shows that the modal parameters vary within the machine tool work volume. en, 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. us, the edge theorem and exclusion condition described in Section 2 are used to perform the robust milling stability prediction. e application of the edge theorem needs to obtain the maximum and minimum values of each modal parameter within the machining space. e 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. erefore, the RBFANN is again used to predict the position and machining allowance-dependent extreme values of the modal parameters.
According to the edge theorem in Section 2.2, there will be 2 21 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]. erefore, 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. en, the defined 64 positions in Section 4.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 x s , y s , and z s , and limiting coordinates of three directions are defined as x lim , y lim , and z lim . en, sample coordinates [x s , x 1 , . . ., x o , . . ., x lim ], [y s , y 1 , . . . , y p , . . ., y lim ], and [z s , z 1 , . . ., z q , . . ., z lim ] of each direction are determined from the intervals [x p x lim ], [y p y lim ], and [z p z lim ], and the adjacent coordinates of each direction have a 40 mm difference. en, the distance between x s and x o , the distance between y s and y p , and the distance between z s and z q 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. en, 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 for x direction and the other one was established to predict 6 extreme values in y direction. en, 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  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.

Robust Milling Chatter Stability Prediction with the Edge eorem and RBFANNs.
Considering the machining within the machining space, six modal parameters of the 4th and 2nd dominant modes in X and Y directions, respectively, are taken as the uncertainties, including ω x4th , k x4th , ξ x4th , ω y2nd , k y2nd , and ξ y2nd . us, 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. e 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. e 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. diameter of 20 mm, a radial cutting width of 8 mm, and tangential and radial cutting force coefficients of 1769 MPa and 1219 MPa. e obtained RSLD is described as a red lobe curve in Figure 7(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. e 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. erefore, the (a p , n) pairs below the red curve is selected to guarantee the milling stability during the real machining process.
Down milling tests have been performed based on the three-axis vertical machining center shown in Figure 4(a) to illustrate the feasibility of the RSLD. Eighteen (a p , 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 (a p , 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 (a p , n) pair has caused a chatter vibration [32]. Comparing the distributions of the stable and unstable (a p , n) pairs shown in Figure 7(a), all (a p , 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 (a p , n) pairs above the blue curve are corresponding to unstable milling conditions, which should be avoided when selecting chatter-free machining parameters; (a p , 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. erefore, the main benefit of the robust milling chatter prediction is that it can provide conservative   a p values to guarantee reliable milling stabilities within the machining space.

Conclusions
is 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. en, 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 in X 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. erefore, 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. e 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 in X, Y, and Z directions. Impact testing was performed to obtain the tool point FRFs and corresponding modal parameters of the 64 positions. en, 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. e robust milling chatter stability prediction was developed based on the 4th sample position, and the RSLD was plotted and compared with the traditional SLD. e RSLD is below the SLD and provides more conservative limiting axial cutting depths. Several (a p , 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. e (a p , n) pairs below the RSLD all correspond to stable conditions, and the stable and unstable (a p , n) pairs are both observed between two lobe curves, validating the feasibility of the RSLD for predicting reliable chatter-free machining parameters. erefore, 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.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.