Prediction of Tool Point Frequency Response Functions Within Machine Tool Work Volume Considering the Position and Feed Direction Dependence

: A chatter vibration in milling process results in poor surface finish and machining efficiency. To avoid the chatter vibration, the stability lobe diagram (SLD) which is the function of tool point frequency response functions (FRFs) is adopted to predict the chatter-free machining parameters. However, the tool point FRF varies with the changes of machining positions and feed directions within machine tool work volume. Considering this situation, this paper presents a method to predict the position and feed direction-dependent tool point FRF. First, modal parameters of the tool point FRFs obtained at some typical positions and feed directions are identified by the modal theory and matrix transformation method. With the sample information, a back propagation (BP) neural network whose inputs are the position coordinates and feed angle and outputs are the modal parameters can be trained with the aid of the particle swarm optimization (PSO) algorithm. Then, modal parameters corresponding to any position and feed direction can be predicted by the trained BP neural network and used to reorganize the tool point FRFs with the modal fitting technique. A case study was performed on a real vertical machining center to demonstrate the accurate prediction of position and feed direction-dependent tool point FRFs. Furthermore, the position and feed direction-dependent milling stability was researched and origin-symmetric distributions of the limiting axial cutting depths at each machining position were observed.


Introduction
Self-excited vibration occurring in the milling process of machine tool is the key factor that results in worsening the machined surface quality, limiting the productivity, deteriorating tool wear and shortening machine tool service life span [1][2][3]. Generally, the stability lobe diagram (SLD) has been proposed to select appropriate chatter-free machining parameters. To plot the SLD, frequency response functions (FRFs) at the tool point should be obtained in advance [4][5][6]. However, since the machine tool spatial structure and the joints contact positions will vary with the movements of the machine tool moving parts, such as the spindle system and worktable, the tool point FRFs are not constant in the machine tool work volume and result in uncertain SLD predictions. Therefore, further research has been developed to predict the position-dependent tool point FRFs, and study variations of the milling stability in the machine tool work volume [7][8][9].
The tool point FRFs can usually be obtained through two approaches. One is the impact testing performed on real machine tool structure, and the other is the modal analysis performed on the mathematic model or finite element model (FEM) in virtual environment [10,11]. Luo et al. [12] considered the influence of machine tool structure change on the tool dynamic properties and proposed a mathematic model to predict the natural frequency by modifying the mass matrix. However, this research has not mentioned how to obtain the tool point FRFs when the machine tool structure varies. Law et al [13] proposed a substructurally synthesized reduced order machine model to realize the theoretical rapid prediction of the position-dependent tool point FRFs by reconstructing the machine tool dynamic model at each position, and then the milling stability for different combinations of machining positions and feed directions were analyzed. Tunc et al. [14] measured the tool tip FRFs at different machining positions and along different feed directions, and these FRFs were further used to plot the SLDs to determine the optimal feed direction to increase the chatterfree material removal rate.
Since the aforementioned methods are still time-consuming and it is difficult to reconstruct the machine tool real structure or FEM to measure or simulate the machine tool dynamics at each machining position, some researches combining the approximate model and experimental design method have been proposed. Baumann et al. [15] determined 23 positions in the machine tool workplace to measure the dynamics of the spindle, and then a triangle-interpolating model was established to predict the dynamics of the spindle at other positions. Deng et al. [16] obtained FRFs of the machine tool frame-holder base at 27 sample positions, and then a Kriging model whose inputs were the moving parts displacements in x, y and z directions was established. And in addition, the predicted FRFs at the base points were further used to predict the FRFs for different tools based on the receptance coupling substructure analysis (RCSA). Similarly, considering the different toolholder assemblies, Liu and Chen et al. [17,18] used transfer learning to predict the pose-dependent tool tip dynamics by integrating domain adaptation and adaptive weighting. Once the posedependent tool tip dynamics were obtained by sufficient impact tests, only few impact tests to measure the new tool tip dynamics were required.
Majority of researches discussed the effects of the varying machining positions on the tool point FRFs and the machining stability. However, the feed direction which also affects the tool point dynamics has been addressed little. Besides the aforementioned research developed by Law and Tunc, Yang et al. [19] investigated the SLDs with different spindle positions and feed directions based on the modal tests. And in addition, the performed cutting tests also validated that the machining stability differences caused by the feed directions. However, these methods are still time-consuming to investigate the position and feed direction-dependent milling stability within whole machine tool work volume. Therefore, this paper proposes a method to predict the tool point FRFs at any machining position and feed direction without complicated impact modal experiments. First, typical machining position and feed direction combinations are determined to obtain the related tool point FRFs with the modal tests and matrix transformation method; then the sample position coordinates and feed angles are taken as the input information and related modal parameters are taken as the output information to train a BP neural network with the aid of the PSO algorithm; thus, the trained and validated BP neural network can be used to predict modal parameters at other combinations of machining position and feed direction, which are further used to reorganize the corresponding tool point FRFs with the modal fitting technique. Then reorganized tool point FRFs are utilized to plot stability lobe diagrams to investigate the position and feed direction-dependent machining stability.
Henceforth, this paper is organized as follows. The principles of the proposed method for predicting the position and feed direction-dependent tool point dynamics are provided in Section 2. In Section 3, a case study describing the application of the proposed PSO-BP method is performed on a real vertical machining center. The position and feed direction-dependent machining stability is investigated through some case studies in Section 4. Finally, conclusions from the current research are summarized in Section 5.

Theoretical Background
Three steps are defined to efficiently predict the position and feed direction-dependent tool point FRFs. First, the modal theory and matrix transformation are combined to obtain the FRFs considering feed directions and identify the corresponding modal parameters. Then, with the information of machining position, feed direction and modal parameters, the particle swarm optimization algorithm is used to benefit the training of the BP neural network. Finally, the trained network is used to predict the position and feed direction-dependent modal parameters, and the related tool point FRFs are reorganized by using the modal fitting technique. A flowchart for describing the aforementioned whole research is provided in Figure 1.

Modal theory and matrix transformation
According to the modal theory, the frequency response function can be expressed as a superposition of multiple single modes. When the excitation force is applied to the system at point p, the vibration response at point l is described as follows [20]: where ωr, ξr, Mr and Ker are the rth modal natural frequency, damping ratio, mass and stiffness respectively.
If the system vibrates near the rth mode, the dynamic response can be obtained by Equation (2): Then the real and imaginary parts of the frequency response function in Equation (2) can be written as follows. λ ξ , and λ rth min can be equal to (1+ξr) if ξr 2 is ignored. Therefore, the vibration frequencies corresponding to the maximum and minimum real part values are ω1 r =ωr(1-ξr) and ω2 r =ωr(1+ξr) respectively, and then the ξr can be obtained by ω1 r and ω2 r using the following Equation (4).
Accordingly, once a frequency response function curve is obtained through impact testing or virtual simulation, the modal parameters for each mode can be identified according equations (4) and (5). The real and imaginary parts of a measured FRF curve are taken as an instance to describe the identification process. The ω1 r , ω2 r and A r for each mode are labeled in Figure 2(a) and used to identify related modal parameters. The identified modal parameters are used to reorganize the measured FRF and well consistence is observed in Figure 2(a).
Generally, the horizontal (x) and normal (y) axes are used to define a machine tool coordinate system (XMT and YMT) shown in Figure 2(b), and the x and y directional tool point FRFs in the machine tool coordinate system are measured to analyze vibrations in related directions. The direct and cross FRFs at the tool point in x and y directions compose the frequency response function matrix Gxy described in Equation (6). However, in real machining process, the feed direction is not always consistent with the x and y axes defined in machine tool coordinate system (XMT and YMT). Therefore, based on the machine tool coordinate system, a rotation angle θ and the feeding directions (u and v) are defined to establish a feed direction coordinate system (Ufeed and Vfeed) as shown in Figure 2 According to the relationship between the machine tool and feed direction coordinate systems, the frequency response function matrix Gxy in machine tool coordinate system is modified by a transformation matrix R described in Equation (6) to obtain the frequency response function matrix Guv in feed direction coordinate system [13,21]. Therefore, if the position-dependent matrix Gxy is measured or simulated, it can be further used to calculate the position and feed direction-dependent matrix Guv with Equation (6).

BP neural network and PSO algorithm
A two degree of freedom (DOF) dynamic model shown in Figure 3(a) is usually established to describe the milling process. In addition, the milling stability is analyzed based on a modal model of the machine and the following characteristic equation [22,23].
where ΛR and ΛI are the real and imaginary parts of the complex eigenvalue Λ, Nt is the tool tooth number, Kt is the tangential cutting force coefficient, ap is the axial cutting depth, ωc is the chatter frequency and T is the tooth passing period determined by the spindle speed n. The G0(iωc) is the directional matrix can be written as Equation (8).
Then, with equations (7) and (8), the limiting axial cutting depth considering the feed direction aplim-feed and the related spindle speed n can be analytically calculated using Equation (9).
According to equations (7) to (9), the prerequisite to determine the milling stability within the machine tool work volume is obtaining the position and feed direction-dependent Guv. Since the matrix Guv is composed by the tool point FRFs, its variations can be represented by the variations of the corresponding FRFs. Considering that the FRFs are organized by the modal parameters as described in the aforementioned section 2.1, variations of the tool point FRFs can be represented by the variations of modal parameters, and further ascribed to the position and feed direction changes. Therefore, a mathematic model needs to be established to predict the modal parameters at any combination of machining position and feed direction.

BP neural network
Since the BP neural network can theoretically approximate any nonlinear continuous function with reasonable topological structure and appropriate structure parameters, it is selected to establish the mathematic model for predicting the position and feed direction-dependent modal parameters [24]. Figure 3(b) shows a three layer BP neural network, and it makes use of error gradient descent algorithm to minimize the mean square error between the actual output value and the output value predicted by network.
As shown in Figure 3(b), each layer is composed by a single or multiple neurons, and the neurons in adjacent layers are connected by the linking weights. The input layer contains the input variable X=[x1, x2,…,xm], and its neurons number equals the number of input variables. The output layer contains the output variable vector Y=[y1, y2,…,yt], and its neuron number equals the number of output variables. Since this paper aims to obtain the modal parameter at any position and feed direction, the displacements of the moving parts in x, y and z directions and the feed angle θ described in Figure 2(b) are taken as the inputs of the BP neural network, and the modal parameters are taken as the outputs of the BP neural network. The neurons number in the hidden layer usually determined according to the empirical equation: l m t a = + + (10) where m and t are the number of neurons in the input and the output layer respectively, and a is an arbitrary integer number between 1 and 10.
For the neuron in hidden or output layer, its output can be calculated by the following Equation (11).
where Oik is the output of the kth neuron in ith layer, f(·) is the excitation function, wkj is the linking weight between the jth neuron in i-1th layer and the kth neuron in ith layer, Ψik is the threshold of the kth neuron in ith layer. The linking weight is modified based on the error function descried by the actual and predicted output values in Equation (12).
where Δw is the linking weight adjustment, η is the constant between 0 and 1 reflecting the learning rate of the BP neural network, and yaq and ypq are the qth actual and predicted output value respectively.
Therefore, if the basic parameters for establishing a BP neural network are determined, such as the neuron number for each layer, the linking weight between each pair of neurons, the threshold value of each neuron, learning rate for the back training algorithm and so on, the BP neural network can be trained to approximate the mathematic relationship between the input and output values.

Particle swarm optimization algorithm
BP neural network has the advantage in approximating the nonlinear relationship between the inputs and outputs, but it shows slow training speed and is easy to fall into the local minimum in the training process. To overcome these shortcomings, the particle swarm optimization algorithm and the genetic algorithm (GA) have been first selected to determine the structure parameters of the BP neural network, and the BP neural network trained with the aid of PSO algorithm has shown higher prediction accuracy. Thus, the particle swarm optimization algorithm which has the advantages of fast convergence, strong robustness and global search ability is finally selected in this paper to help training the BP neural network [25,26].
In the D-dimensional search space, m particles form a population. For each particle, it has two attributes including the position and velocity. For the ith particle, its position vector can be expressed as Xi=(xi1, xi2,…,xiD). During the searching process, each current particle position is a candidate solution to the corresponding optimization problem, and the fitness function is determined to judge the historical optimal position of each particle and the historical optimal position of the population. Each current individual optimal position and the global optimal position are further used to determine the corresponding velocity vector of each particle. For the ith particle, the determined velocity vector Vi=(vi1, vi2,…,viD) is used to modify the current particle position to close to the historical optimal one. Then, the new position for each particle is obtained to continue the next iteration. During the iteration, the position and velocity of each particle are repeatedly updated using the following Equation (13). The iteration terminates until the best global fitness value or the total iteration number meets the termination condition.
where k is the current iteration number, ω is the inertia weight, c1 and c2 are the acceleration factors, r1 and r2 are the random numbers between 0 and 1. In order to establish the BP neural network with the aid of PSO algorithm, the weights and thresholds are taken as the particles. Then positions and velocities of the particles are updated iteratively to realize the optimal positions. The application flow chart of the combined BP neural network and PSO algorithm for predicting the position and feed direction-dependent modal parameters is shown in Figure 4, and the related procedures are summarized as follows.
Step 1: Initialize the BP neural network. Define the neuron numbers of the input layer, hidden layer and output layer according to section 2.2.1, and normalize the input and output sample data between -1 and 1. Other basic parameters should also be determined, such as the learning rate, excitation function, training goal and so on.
Step 2: Initialize required parameters of the PSO algorithm. First, the particle dimension D is determined by the total number of the weights and thresholds using Equation (14).
where Ninput and Noutput are the neuron numbers in the input and output layers, Nht is the total number of hidden layers, and Nhi is the neuron number of the ith hidden layer. Then, the population size, iteration number, inertia weight, initial positions and velocities, acceleration factors, and rand numbers are initialized. Moreover, the lower and upper boundaries of these positions and velocities are also defined. The error function described in Equation (12) is taken as the fitness function.
Step 3: Obtain the optimal initial weights and thresholds of the BP neural network. During each iteration of the PSO algorithm, the individual and global optimal positions are updated and recoded by comparing the so far calculated particle fitness values, and the corresponding individual and global optimal fitness values are also updated and recoded. When the iteration terminates, the global optimal position is taken as the optimal initial weights and thresholds to train the BP neural network according to section 2.2.1. In addition, the accuracy of the trained BP neural network is verified by the testing sample including the information of the input and output variables.

Case study
In this section, the proposed BP-PSO method for predicting the position and feed directiondependent tool point FRFs was applied to a real three-axis vertical machining center shown in Figure  5, the worktable, saddle and spindle system of which move along the x, y, and z directions respectively. The movement ranges of the worktable, saddle, and spindle system are 0-550mm, 0-400mm and 0-400mm respectively, and the corresponding displacements of the three axes are combined to represent the machining position.

Impact testing based on orthogonal experiment design
To establish the BP neural network for predicting the tool point FRFs within the machine tool work volume, the commonly used orthogonal experimental design is first introduced to determine some representative sample machining positions and feed directions. The orthogonal experimental design has the characteristics of uniform distribution and good comparability, and it can benefit obtaining the typical position and feed direction-dependent tool point FRFs with fewer modal experiments [27]. Then, the displacements in x, y and z directions and the feed direction angle θ are taken as the factors, and eight levels for each factor are determined within its variation range. The corresponding factors and levels are listed in Table 1, and an orthogonal table L64 (8 4 ) is used to determine 64 experiment schemes shown in Table 2.
The first four columns of Table 2 represent the typical combinations of machining position and feed direction. At each combination, the worktable, saddle, and spindle system are driven to achieve the machining position, and the impact testing was performed at the tool point of each machining position. During the experiment, the impact hammer with the type PCB 086D05, sensibility 0.23

BP-PSO method for modal parameters prediction
With the obtained machining position and feed direction-dependent tool point FRFs in section 3.1, corresponding modal parameters can be identified according to equations (4) and (5). Analyzing these tool point FRF curves, three obvious modes for each curve can be observed, and then the natural frequencies, modal stiffnesses and modal damping ratios for the three modes were identified. One FRF in v direction is taken as an example in Figure 7(a), and the identified modal parameters are listed in Table 3. Considering the u and v directions, the BP neural network should be established to predict the modal parameters of six modes. Predicting 18 modal parameters simultaneously will result in over-complicated structure and slow convergence speed. Therefore, six BP neural networks were established, and each neural network is used to predict the natural frequency, modal stiffness, and modal damping ratio of each mode. Then, the modeling process of a single BP neural network with the aid of PSO algorithm is described to give an illustration. Taking the prediction of the second modal parameters of the tool point FRF in u direction as an instance, the input layer has 4 neurons, which represent the displacements of x, y and z directions and the feed direction angle θ; the output layer has 3 neurons, including the natural frequency ω, modal stiffness K and modal damping ratio ζ; one hidden layer is defined, and 3 neurons are determined according to Equation (10). Then, in the PSO algorithm, 100 particles are defined to compose one population, and each particle has 27 dimensions calculated according to Equation (14) and represents the initial linking weights and thresholds of the BP neural network; the fitness function is the mean squared errors between the real and predicted modal parameters; the minimum and maximum velocity are -0.5 and 0.5; the inertia weight variation range is defined as [0.4, 1.2], and the inertia weight is first defined as 1.2 and decreases with the iteration number linearly; the initial velocity vector of each particle was defined randomly, and the acceleration constants were 1.45.
According to section 3.1, 64 combinations of three directional displacements and feed direction angle (x, y, z, θ) in Table 2 were defined as the input sample, and the corresponding 64 combinations of the modal parameters were taken as the output sample. Moreover, 58 combinations were selected randomly to be the training sample, and the other 6 combinations were determined as the testing sample. Then, the BP neural network was trained based on the PSO algorithm and the training sample. Using the testing sample to verify the prediction accuracy of the trained BP neural network, comparisons of the original and predicted modal parameters values for the testing sample are described in Table 4. Relatively small errors are observed, validating the accuracy of the trained BP neural network, which can be used to predict the position and feed direction-dependent modal parameters.
After the six BP neural networks were established and trained, the predicted 18 modal parameters can be further used to reorganize related position and feed direction-dependent tool point FRFs. Figure 7(b) is an example of the measured and predicted FRFs along u and v directions at one position, and good consistencies are observed. Therefore, the predicted position and feed directiondependent tool point FRFs can be used to analyze the milling stability within the whole machine tool work volume. Table 4. Comparisons between the measured and predicted modal parameters for the second mode in u direction for the six testing samples.

Sample
No.

Chatter stability analysis considering uncertain position and feed direction
Main purpose of the milling stability prediction is selecting appropriate machining parameters to avoid the chatter vibration. Since the milling stability depends on the tool point FRFs which is a function of the machining positions and feed direction, the milling stability will vary in the whole machine tool work volume. Therefore, this section describes the position and feed directiondependent milling stability simulation for downing milling the workpiece whose material is the commonly used ASTM 1045 steel with tool diameter 20mm and tooth number Nt=4, tangential cutting force coefficient Kt=1799MPa and radial cutting force coefficient Kr=760 MPa.
First, the left-most, middle and right-most positions of each direction are defined, and the feed direction-dependent milling stabilities related to the three positions were researched based on the equations (7) to (9) and the established BP neural networks in section 3. When the position of one moving part varies along one direction, no movements of the other two directions were defined. Taking the spindle speed 8400 rpm as an instance, the calculated axial limiting cutting depth aplim_feed values are plotted in Figure 8.   Figure 8(a), when the feed direction angle varies from 0º to 360º at the three positions, the maximum and minimum aplim_feed values are 7.20 mm and 6.42 mm respectively, and the variation rate is 12.15%. In Figure 8(b), the maximum and minimum aplim_feed values are 5.66 mm and 5.24 mm respectively, and the variation rate is 8.02%. In Figure 8(c), the maximum and minimum aplim_feed values are 9.08 mm and 6.02 mm respectively, and the variation rate is 50.83%. Comparing three figures, Figure 8 Furthermore, since the limiting axial cutting depth is affected by the tool passing period T which depends on the spindle speed n shown in Equation (7), the spindle speed was taken as a factor and its eight levels were determined and listed in Table 1; then, the first five columns of Table 2 were used to perform the milling stability analysis, and the calculated axial limiting cutting depth aplim_feed values are shown in Table 2. Thus, the range analysis and variance analysis described as follows were adopted to comprehensively study the effects of the machining position (x, y, z), feed directions θ and spindle speed n on the milling stability.
In the range analysis, the range R was calculated according to Equation (15), which were used to investigate the effect degree of each factor on the axial limiting cutting depth aplim_feed and determine the optimal level of each factor. A factor with higher range value indicates that it shows a greater effect on the aplim_feed [29].
where Rj is the range value of the jth factor, kij is the sum of the aplim_feed values calculated using the ith level of the jth factor, i=1, 2, …, m, j=1, 2, …, n, and m and n are the level and factor numbers of the orthogonal table respectively.
In the variance analysis, the F-test was used to determine the effect degree of each factor on the axial limiting cutting depth aplim_feed, and the parameters needed to perform the F-test were calculated according to Equation (16) [30].
where xp is the aplim_feed value calculated using the pth scheme designed in the orthogonal table, Nump is the total number of designed schemes, SSj is the sum of squared deviation of the jth factor, dfT is the total degree of freedom, dfj is the degree of freedom for the jth factor, Sj is the variance of the jth factor, and Fj is the F value of the jth factor. The significance level of each factor is determined by comparing the obtained F value with the standard F value (Fα(dfj, dfe)). When the significant level α is given, the standard value Fα(dfj, dfe) can be attained from the F-table in accordance with the degree of freedom [30]. If the F value is higher than the standard Fα value, the factor is regarded as significant. And in addition, the factor with a higher F value means that it has more evident impact on the aplim_feed. Therefore, the information in Table 2 were used to perform the range analysis and variance analysis with the equations (15) and (16). Comparing the kij value of each factor listed in Table 5, the level with a higher kij value is regarded as the optimal level of the related factor. Thus, the 4th, 4th, 1st, 6th, 8th levels of the factor x, y, z, θ and n are determined as the optimal combination to have a bigger aplim_feed, and specific values of the subscripts are listed in Table 1. For the variance analysis, the calculated SST is 1029.10, SSe is 83.31, dfT is 63, dfe is 28, and other corresponding parameters are listed in Table 5. In this current research, the significant level α is defined as 0.05, and the obtained standard F values Fα are listed in Table 5. Comparing the F and Fα values of each factor listed in Table 5, the F values of the spindle speed n, z directional displacement, feed direction angle θ and x directional displacement are bigger than the F0.05 (7,28)=2.359, which indicates that these four factors have significant effects on aplim_feed with a confidence level of 95%. Further comparing the F values of the five factors, the influence degree of the machining position (x, y, z), feed directions θ and spindle speed n on axial limiting cutting depth aplim_feed is n> z >θ> x > y. Analyzing the R and F values in Table  5, the spindle speed shows the dominant effects on the aplim_feed, the displacement in z direction and the feed direction angle θ have similarly significant effects on the aplim_feed, the displacement in x direction shows less evident influence on the aplim_feed, and the displacement in y direction has no evident effect on the aplim_feed. Accordingly, when designing the tool path, an optimal combination of z directional position and feed direction needs to be determined, and the movements in x and y directions should be first taken into consideration.

Conclusions
Considering the machining position and feed direction have significant effects on the tool point FRFs and thus on the milling stability, this paper proposes a method for predicting the machining position and feed direction-dependent tool point FRFs based on the modal theory, matrix transformation, BP neural network and PSO algorithm. The feasibility of the proposed method was verified by case studies on a real vertical machining center.
(1) The displacements of the worktable, saddle and spindle system along the x, y and z directions represents the machining position, and the feed direction is represented by the rotation angle θ around the machine tool coordinate system. 8 levels of each factor were determined within its variation range, and then an orthogonal table with 64 schemes was designed. For each scheme, the machine tool was driven to the related spatial structure to measure the x and y directional tool point FRFs, and the matrix transformation method was used to calculate the related FRFs in u and v directions. (2) Three obvious modes for each tool point FRF were observed, and the related modal parameters were identified according to modal theory. Thus, considering the u and v directions, 18 modal parameters should be predicted for one machining position and feed direction combination. To avoid the complex training process and improve the convergence speed, 6 BP neural networks were established for predicting the modal parameters of six modes respectively. Then these BP neural networks were trained with the information of the randomly selected 58 schemes from the Table 2 and the optimized initial linking weights and thresholds using the PSO algorithm. The accuracies of established BP neural networks were verified by the measured modal parameters related to the combinations of machining position and feed direction included in the testing sample. (3) Moreover, these BP neural networks were used to research the milling stability within the machine tool work volume. Distributions of the axial limiting cutting depth within the feed direction angle range [0°, 360°] at 9 selected machining positions show that the z directional displacement and feed direction θ have significant effects on the limiting axial cutting depth. The same phenomenon was also observed from results of the range analysis and variance analysis performed according to the information in orthogonal Table 2. The optimal level for each factor was determined to obtain an optimal combination of machining position and feed direction with higher limiting axial cutting depth. Therefore, the proposed method for predicting the machining position and feed directiondependent FRFs can benefit an optimal process planning according to real milling conditions at the design stage. In further research, the proposed method can be extended to take more factors into consideration, such as the spindle speed and tool change, to make the FRFs prediction more fit in with the actual milling process.