Natural Frequency Prediction Method for 6R Machining Industrial Robot

: The industrial robot machining performance is highly dependent on dynamic behavior of the robot, especially the natural frequency. This paper aims at introducing a method to predict the natural frequency of a 6R industrial robot at random conﬁguration, for improving dynamic performance during robot machining. A prediction model of natural frequency which expresses the mathematical relation between natural frequency and conﬁguration is constructed for a 6R robot. Joint angles are used as input variables to represent the conﬁgurations in the model. The quantity and range of variables are limited for e ﬃ ciency and practicability. Then sample conﬁgurations are selected by central composite design method due to its capacity of disposing nonlinear e ﬀ ects, and natural frequency data is acquired through experimental modal test. The model, which is in form of regression equation, is ﬁtted and optimized with sample data through partial least square (PLS) method. The proposed model is veriﬁed with random conﬁgurations and compared with the original model and a model ﬁtted by least square method. Prediction results indicate that the model ﬁtted and optimized by PLS method has the best prediction ability. The universality of the proposed method is validated through implementation onto a similar 6R robot.


Introduction
Due to the rapid development of industrial automation, the global amount of industrial robot keeps increasing significantly in recent years [1], and new applications attract more attention [2][3][4][5]. Because of the multiple degree of freedom, large workspace, and relative low cost, industrial robot becomes a potential alternative of traditional machine tool in hard material drilling [6], boring [7], and milling [8,9] processes, which are not maturely and widely applied.
For machine tool, in the whole workspace, the stiffness and dynamic behavior nearly remain steady, so machining stability and chatter suppression, which are concerned with multiple factors, are research emphasis [10][11][12]. However, for industrial robot, current researches focus on the characteristic of robot itself. Common industrial robot is of open-loop articulated serial structure with six revolute joints (6R); therefore, the stiffness and dynamic performance are much inferior to machine tool, which results in vibration issues and poor surface quality in machining process [13]. Static stiffness issues are studied to avoid the robot deformation in machining, and several stiffness indexes are proposed to evaluate the static stiffness of robot in the workspace [14][15][16][17]. On the other hand, dynamic behavior of robot is studied to evaluate the capacity of dealing with time-variant load in machining process, which may help to suppress vibration [18,19]. As 6R robot realizes movement through adjusting joint angles, dynamic behavior keeps varying with configurations during motion process. For that reason, research object is gradually expanded from one configuration to the whole workspace. Bisu et al. [20] selected three configurations from a linear movement path of a 6R robot with interval of 50 cm for dynamic behavior study and vibration analysis. It was found that different natural frequencies led to dynamic performance variation. Mousavi et al. [21] performed a similar four-point-study and analyzed the change of stability region through modal information. Several configurations may present a linear movement, but the usage of robot should not be such limited. Karim et al. [22] executed modal test with 63 measured points for 23 configurations, through which mode shapes were obtained and approximate distribution of first two step natural frequencies in a plane was depicted. Natural frequency is the vital dynamic parameter in robot machining, but massive tests are inefficient although the accuracy is well ensured. A method was proposed by Glogowski et al. [23] to predict natural frequency in the whole workspace efficiently, in which sample configurations were selected by design of experiment (DoE) method and prediction model was fitted by least square (LS) method with the natural frequency data of sample configurations, but the prediction deviation reached up to 8.9%, making the method inadequate for practical application.
The low order natural frequencies of 6R robot are relatively small and the first order is merely about 10 Hz, which is covered in frequency range of most machining processes, so 6R robot tends to forced vibration during machining processes [24]. To obtain basis for vibration suppression, structure optimization and path planning, low order natural frequencies are particularly concerned among dynamic parameters. It is proposed that low order natural frequencies of robot are decided by the configurations, while the high order natural frequencies are mainly related to the machining system, which do not vary with robot configurations [25], this phenomenon is demonstrated in Figure 1.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 2 of 14 process. For that reason, research object is gradually expanded from one configuration to the whole workspace. Bisu et al. [20] selected three configurations from a linear movement path of a 6R robot with interval of 50 cm for dynamic behavior study and vibration analysis. It was found that different natural frequencies led to dynamic performance variation. Mousavi et al. [21] performed a similar four-point-study and analyzed the change of stability region through modal information. Several configurations may present a linear movement, but the usage of robot should not be such limited. Karim et al. [22] executed modal test with 63 measured points for 23 configurations, through which mode shapes were obtained and approximate distribution of first two step natural frequencies in a plane was depicted. Natural frequency is the vital dynamic parameter in robot machining, but massive tests are inefficient although the accuracy is well ensured. A method was proposed by Glogowski et al. [23] to predict natural frequency in the whole workspace efficiently, in which sample configurations were selected by design of experiment (DoE) method and prediction model was fitted by least square (LS) method with the natural frequency data of sample configurations, but the prediction deviation reached up to 8.9%, making the method inadequate for practical application. The low order natural frequencies of 6R robot are relatively small and the first order is merely about 10 Hz, which is covered in frequency range of most machining processes, so 6R robot tends to forced vibration during machining processes [24]. To obtain basis for vibration suppression, structure optimization and path planning, low order natural frequencies are particularly concerned among dynamic parameters. It is proposed that low order natural frequencies of robot are decided by the configurations, while the high order natural frequencies are mainly related to the machining system, which do not vary with robot configurations [25], this phenomenon is demonstrated in Figure 1. Thus, a method to predict natural frequency of 6R industrial robot efficiently and precisely is proposed in this paper. As shown in Figure 2, sample configurations are selected to present the practical workspace (variable range) for fast prediction, and a regression equation fitted and optimized by partial least square (PLS) method is used as prediction model for higher accuracy. A 6R robot, ABB IRB4600, is used to illustrate and validate the method. Another ABB IRB6400 is used verify the universality of proposed method. Thus, a method to predict natural frequency of 6R industrial robot efficiently and precisely is proposed in this paper. As shown in Figure 2, sample configurations are selected to present the practical workspace (variable range) for fast prediction, and a regression equation fitted and optimized by partial least square (PLS) method is used as prediction model for higher accuracy. A 6R robot, ABB IRB4600, is used to illustrate and validate the method. Another ABB IRB6400 is used verify the universality of proposed method.

Sample Preparation
The axes of 6R robot are noted as a1, a2, …, a6, corresponding joint angles are q1, q2, …, q6. The mapping relation between natural frequency (fF) and configuration can be expressed as With samples, a regression equation can be fitted to represent the mapping relation. However, it is unnecessary and inefficient to include the entire workspace of the robot as the configurations near the boundary could not be adopted in machining process. Therefore, range of variables is limited before selecting samples through DoE method. Then natural frequency information of each sample configuration is acquired through modal test.

Variable Range Determination
a1, which is the basic rotation axis of the robot structure, and a6, which is a 360° rotation axis controlling the flange, are proved to be nonsignificant in changing natural frequency by preliminary researches on 6R robots including KUKA KR60, ABB IRB6400, and IRB 4600. As low order natural frequencies is mainly decided by robot configuration, none end-effector is fixed on the flange in above cases for general research. Hereby, q1 and q6 are excluded from model (q1 = q6 = 0° in practice) for the simplification, and q2-q5 are applied as variables. Prediction model can be adapted to To ensure the efficiency of model, partial workspace that may cause movement interference in machining process is rejected by simulation software. The variable range is finally determined and shown in Figure 3, where the solid area indicates the new efficient workspace, and the dashed area represents the original workspace.

Sample Preparation
The axes of 6R robot are noted as a 1 , a 2 , . . . , a 6 , corresponding joint angles are q 1 , q 2 , . . . , q 6 . The mapping relation between natural frequency (f F ) and configuration can be expressed as With samples, a regression equation can be fitted to represent the mapping relation. However, it is unnecessary and inefficient to include the entire workspace of the robot as the configurations near the boundary could not be adopted in machining process. Therefore, range of variables is limited before selecting samples through DoE method. Then natural frequency information of each sample configuration is acquired through modal test.
2.1.1. Variable Range Determination a 1 , which is the basic rotation axis of the robot structure, and a 6 , which is a 360 • rotation axis controlling the flange, are proved to be nonsignificant in changing natural frequency by preliminary researches on 6R robots including KUKA KR60, ABB IRB6400, and IRB 4600. As low order natural frequencies is mainly decided by robot configuration, none end-effector is fixed on the flange in above cases for general research. Hereby, q1 and q6 are excluded from model (q 1 = q 6 = 0 • in practice) for the simplification, and q 2 -q 5 are applied as variables. Prediction model can be adapted to To ensure the efficiency of model, partial workspace that may cause movement interference in machining process is rejected by simulation software. The variable range is finally determined and shown in Figure 3, where the solid area indicates the new efficient workspace, and the dashed area represents the original workspace. Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 14

Sample Configuration Selection
DoE method has advantages in study influence of different factors, especially for basic research of robot machining, e.g., Simoes et al. [26] used DoE to study influence of robot milling parameters. Different types of selection methods are compared in pilot study including three-level full factorial design (3kD), orthogonal experimental design, Box-Behnken design, and central composite design (CCD). The prediction performances of 3kD and CCD are nearly matched and superior to the other two, while CCD demands a much smaller sample amount [23]. Meanwhile CCD owns conspicuous ability in disposing the nonlinear influence that joint angles may lead to natural frequency. Therefore, CCD method is adopted to select sample configurations in this paper, corresponding factors and levels are shown in Table 1. Experimental modal test (EMT) is applied so that the natural frequency data of sample configurations can be identified accurately. According to principle of modal test, any element of frequency response function (FRF) can be written as independent of excitation and response combinations. Theoretically, one combination of excitation and response is enough to identify all the natural frequencies, provided the excitation and response points avoid the nodes. The sample configurations selected by CCD design vary greatly, for the convenience of implementing all the trials in the CCD design, hammer excitation is chosen, and a three-direction accelerometer is fixed on the flange of robot to acquire response signal. In this way, the extra mass brought by EMT device can be negligible as the accelerometer is small and light compared to the

Sample Configuration Selection
DoE method has advantages in study influence of different factors, especially for basic research of robot machining, e.g., Simoes et al. [26] used DoE to study influence of robot milling parameters. Different types of selection methods are compared in pilot study including three-level full factorial design (3kD), orthogonal experimental design, Box-Behnken design, and central composite design (CCD). The prediction performances of 3kD and CCD are nearly matched and superior to the other two, while CCD demands a much smaller sample amount [23]. Meanwhile CCD owns conspicuous ability in disposing the nonlinear influence that joint angles may lead to natural frequency. Therefore, CCD method is adopted to select sample configurations in this paper, corresponding factors and levels are shown in Table 1.

Sample Data Acquisition
Experimental modal test (EMT) is applied so that the natural frequency data of sample configurations can be identified accurately. According to principle of modal test, any element of frequency response function (FRF) can be written as which presents the response of i point caused by the excitation at j point. a ijk and a * ijk are related to mode shapes, their values are depending on relevant combination of excitation and response. p k and p * k contain the information of natural frequency and damping, which are global characteristics independent of excitation and response combinations. Theoretically, one combination of excitation and response is enough to identify all the natural frequencies, provided the excitation and response points avoid the nodes.
The sample configurations selected by CCD design vary greatly, for the convenience of implementing all the trials in the CCD design, hammer excitation is chosen, and a three-direction accelerometer is fixed on the flange of robot to acquire response signal. In this way, the extra mass brought by EMT device can be negligible as the accelerometer is small and light compared to the robot. Excitations are generated from three directions for each trial, which are corresponding to the three directions of the accelerometer, as shown in Figure 4. Then three FRFs will be obtained for natural frequency identification, avoiding incomplete excitation or excitation on the nodes. Moreover, for each direction, five groups of excitation and response signals are acquired to decrease the signal-to-noise ratio through average algorithm. With FRFs the natural frequencies are preliminarily identified, then complex modal indicating function is used to confirm the accurate values of natural frequencies.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 14 robot. Excitations are generated from three directions for each trial, which are corresponding to the three directions of the accelerometer, as shown in Figure 4. Then three FRFs will be obtained for natural frequency identification, avoiding incomplete excitation or excitation on the nodes. Moreover, for each direction, five groups of excitation and response signals are acquired to decrease the signalto-noise ratio through average algorithm. With FRFs the natural frequencies are preliminarily identified, then complex modal indicating function is used to confirm the accurate values of natural frequencies. With the ABB IRB4600 industrial robot being the object, a EMT system for sample data acquisition is set up, including a LIXIE hammer, a PCB three-direction accelerometer, an ECON device for vibration signal processing and a PC (as shown in Figure 5). The configuration of robot is adjusted according to the CCD design in turn. All the identified natural frequency information will be used as the sample data to construct a prediction model.

Prediction Model Construction and Optimization
Because of the serial structure, the relation between configuration and joint angles could not be linear. Thus, Equation (2) is fitted to a general second order equation in this case according to normal engineering application. PLS method is adopted to solve the issues of nonlinear influence and multicorrelation between joint angles in regression equation fitting. Meanwhile, the fitting error caused by the limitation of sample quantity can be remarkably reduced, due to the advantage of PLS method in coping with small sample quantity. In addition, PLS has good interpretability for the terms of regression equation, which can be used as optimization criterion. With the ABB IRB4600 industrial robot being the object, a EMT system for sample data acquisition is set up, including a LIXIE hammer, a PCB three-direction accelerometer, an ECON device for vibration signal processing and a PC (as shown in Figure 5). The configuration of robot is adjusted according to the CCD design in turn. All the identified natural frequency information will be used as the sample data to construct a prediction model. robot. Excitations are generated from three directions for each trial, which are corresponding to the three directions of the accelerometer, as shown in Figure 4. Then three FRFs will be obtained for natural frequency identification, avoiding incomplete excitation or excitation on the nodes. Moreover, for each direction, five groups of excitation and response signals are acquired to decrease the signalto-noise ratio through average algorithm. With FRFs the natural frequencies are preliminarily identified, then complex modal indicating function is used to confirm the accurate values of natural frequencies. With the ABB IRB4600 industrial robot being the object, a EMT system for sample data acquisition is set up, including a LIXIE hammer, a PCB three-direction accelerometer, an ECON device for vibration signal processing and a PC (as shown in Figure 5). The configuration of robot is adjusted according to the CCD design in turn. All the identified natural frequency information will be used as the sample data to construct a prediction model.

Prediction Model Construction and Optimization
Because of the serial structure, the relation between configuration and joint angles could not be linear. Thus, Equation (2) is fitted to a general second order equation in this case according to normal engineering application. PLS method is adopted to solve the issues of nonlinear influence and multicorrelation between joint angles in regression equation fitting. Meanwhile, the fitting error caused by the limitation of sample quantity can be remarkably reduced, due to the advantage of PLS method in coping with small sample quantity. In addition, PLS has good interpretability for the terms of regression equation, which can be used as optimization criterion.

Prediction Model Construction and Optimization
Because of the serial structure, the relation between configuration and joint angles could not be linear. Thus, Equation (2) is fitted to a general second order equation in this case according to normal engineering application. PLS method is adopted to solve the issues of nonlinear influence and multi-correlation between joint angles in regression equation fitting. Meanwhile, the fitting error caused by the limitation of sample quantity can be remarkably reduced, due to the advantage of PLS method in coping with small sample quantity. In addition, PLS has good interpretability for the terms of regression equation, which can be used as optimization criterion.

PLS Method Fitting Process
The principle of PLS method in multivariate regression fitting is to extract principal components t and u respectively from independent variables X and dependent variables Y while retaining the most correlation, which can be mathematically described as The process of regression equation fitting by PLS method [27] is briefly introduced as follows.
Current principal component t h is validated by following indexes: where S SPE,h is the square sum of prediction error of Y h−1 , S SE,h is square sum of deviation of Y h−1 . R 2 Y (cum) indicates the ability of principal component to explain the variability of Y h−1 , and its upper limit value is 1. Q 2 h represents the contribution of t h to the model. When it is more than (1-0.95) 2 = 0.0975, t h has significant effect on prediction model. Meanwhile, Q 2 h (cum) increases, meaning that the comprehensive significance in explaining Y is enhanced.

4.
Subsequent principal components calculation. The residual matrixes from Step c. are continued to compute new principal components, until t Q+1 deduces Q 2 h (cum). The principal components before t Q+1 are validated ones.

5.
Regression equation 6. Data reduction. The final regression equation for prediction can be obtained through the reverse process of data standardization.

Interpretability of PLS Method for Regression Equation
Variable importance in projection (VIP), which indicates the importance of each item in regression equation, is defined as where w jh is the jth component of w h . The larger VIP j is, the more important x j is in constructing Y. The effects of different items can be compared quantitatively through VIP values, so that some redundant items can be rejected to realize regression equation optimization.

Prediction Model Construction
The sample data from Section 2.1.3 is utilized to construct prediction model, taking the first order natural frequency of ABB IRB4600 robot as example.
As mentioned, the general regression equation concludes q 2 , q 3 , q 4 , q 5 , and their cross items and square items. The general equation is marked as M. After the second principal component is extracted, validation indicates are shown in Table 2. The second principal component deduces Q 2 h (cum), meaning that comprehensive ability is weakened. Though R 2 Y (cum) increases, the second principal component is still abandoned.   Figure 6. All the one-degree items should be retained as basic items. All the square items are important enough to stay in the model because of their high VIP values. As VIP values of cross items are relatively low, different combination groups should be tested. A new model M1 is constructed by cutting off cross items except q2q3, because q2q3 has the largest VIP value among cross items; meanwhile, q2 and q3 are most two important variables according to VIP values. Based on M1, the other five cross items are added respectively in descending order of their VIP values. The new five models are marked as M2-M6. Validation results of above models are shown in Table 3. All the six models are constructed with two principal components. 2    All the one-degree items should be retained as basic items. All the square items are important enough to stay in the model because of their high VIP values. As VIP values of cross items are relatively low, different combination groups should be tested. A new model M 1 is constructed by cutting off cross items except q 2 q 3 , because q 2 q 3 has the largest VIP value among cross items; meanwhile, q 2 and q 3 are most two important variables according to VIP values. Based on M 1 , the other five cross items are added respectively in descending order of their VIP values. The new five models are marked as M 2 -M 6 . Validation results of above models are shown in Table 3. All the six models are constructed with two principal components. R 2 Y (cum) and Q 2 h (cum) values of the six models are all improved obviously compared to M. M 2 , formed by adding q 3 q 5 to M 1 , has the best R 2 Y (cum) and Q 2 h (cum) values. M 3 -M 6 , which get lower R 2 Y (cum) and Q 2 h (cum) values compared to M 1 , are abandoned. Based on M 2 , four new models are introduced by adding q 2 q 4 , q 3 q 4 , q 2 q 5 , and q 4 q 5 in turn. However, the largest Q 2 h (cum) only reaches up to 0.764. All the other combination groups of above four cross items are add to M 2 for test, but none of them can acquire better validation result. Thus, the optimized model is finally confirmed to be M 2 . The coefficients of items are listed in Table 4 together with M.

Model Verification with Random Configurations
Sixteen groups of q 2 -q 5 values are picked randomly from the variable range defined in Section 2.1.1, through which 16 corresponding robot configurations are defined for model verification. The simulation images of the 16 configurations are reranked by joint angles q 2 and q 3 as shown in Figure 7, and the measured natural frequency is listed in the bottom left corner of each image. The configurations are nearly uniformly distributed in the variable range, which ensures the comprehensiveness of verification objects and rationality of verification result.  1.099 × 10 −5 -

Model Verification with Random Configurations
Sixteen groups of q2-q5 values are picked randomly from the variable range defined in Section 2.1.1, through which 16 corresponding robot configurations are defined for model verification. The simulation images of the 16 configurations are reranked by joint angles q2 and q3 as shown in Figure  7, and the measured natural frequency is listed in the bottom left corner of each image. The configurations are nearly uniformly distributed in the variable range, which ensures the comprehensiveness of verification objects and rationality of verification result. The deviation between prediction value and measured value is defined as prediction error, the details are listed in Table 5; Table 6.  The deviation between prediction value and measured value is defined as prediction error, the details are listed in Table 5; Table 6. The average error and average relative error of optimized model M 2 are both much lower than those of original model M. The average relative error of M 2 is less than 5% which brings more reliability. On the other hand, seven high relative errors (more than 5%) which are underlined show up in the prediction process of M, and the maximum error reaches up to 1.78 Hz. While M 2 has no predictions with over 5% relative error, and the maximum error is 0.51 Hz which is acceptable. Obviously, optimized model M 2 has better prediction ability.
A model (M LS ) fitted by least square (LS) method with the same sample data is taken as the contrast. Moreover, the prediction errors are displayed in Table 7. Average error and average relative error of M LS are both higher than those of M and M 2 , that is, PLS method have better performance in construct prediction model of natural frequency than LS method. In conclusion, model fitted and optimized by PLS method has the best prediction ability.

Model Construction Quality Analysis
In PLS method, when model is constructed by principal components, information in the last residual matrix is ignored, causing the error between fitted model and original data. Standardized distance between sample point and model can be used to evaluate the construction quality. (DModX,N) i or (DModY,N) i is the ratio of construct quality of the ith sample point and the average construct quality. When the ratio is less than two, the construction for the ith sample point is reasonable. Construction quality of M 2 is shown in Figure 8. (DModX,N) i and (DModY,N) i are all less than two, that means the construction of M 2 is rational.

Method Universality Verification
To testify the universality of the natural frequency prediction method proposed in this paper, it is completely applied to an ABB IRB6400 industrial robot (as shown in Figure 9). Sample data are obtained according to Chapter 2.1 as the structure is similar to IRB4600. PLS is used to fit and optimize the prediction model M6400. In this case, q2q4, q2q5, and q4q5 are abandoned according to VIP values. Validation results can be seen in Table 8, values of 2 (cum) Y R and 2 h (cum) Q indicate that the fitting is rational. Construction quality is eligible as shown in Figure 10. The result of model verification through random configurations is shown in Table 9, average relative error is merely 2.982%, indicating the good prediction performance.

Method Universality Verification
To testify the universality of the natural frequency prediction method proposed in this paper, it is completely applied to an ABB IRB6400 industrial robot (as shown in Figure 9). Sample data are obtained according to Chapter 2.1 as the structure is similar to IRB4600. PLS is used to fit and optimize the prediction model M 6400 . In this case, q 2 q 4 , q 2 q 5 , and q 4 q 5 are abandoned according to VIP values. Validation results can be seen in Table 8, values of R 2 Y (cum) and Q 2 h (cum) indicate that the fitting is rational. Construction quality is eligible as shown in Figure 10. The result of model verification through random configurations is shown in Table 9, average relative error is merely 2.982%, indicating the good prediction performance.

Method Universality Verification
To testify the universality of the natural frequency prediction method proposed in this paper, it is completely applied to an ABB IRB6400 industrial robot (as shown in Figure 9). Sample data are obtained according to Chapter 2.1 as the structure is similar to IRB4600. PLS is used to fit and optimize the prediction model M6400. In this case, q2q4, q2q5, and q4q5 are abandoned according to VIP values. Validation results can be seen in Table 8, values of 2 (cum) Y R and 2 h (cum) Q indicate that the fitting is rational. Construction quality is eligible as shown in Figure 10. The result of model verification through random configurations is shown in Table 9, average relative error is merely 2.982%, indicating the good prediction performance.        The two models constructed for two different robots through the method proposed in this paper both own good prediction ability. That is, the universality of the proposed method is testified.

Discussion
Once the prediction model of natural frequency is constructed and verified, the prediction result can be used as optimization parameter to improve the machining performance. An application case of the prediction model is explained as follows. a 2 and a 3 control the two longest links, by which the basic configuration of robot is determined. Through VIP values, items involving q 2 and q 3 are demonstrated to act dominated role in prediction model. a 4 and a 5 decide the direction of robot flange, which affect the configuration less than a 2 and a 3 . Taking M 2 for an example, the influence of q 2 and q 3 on natural frequency is specifically studied and illustrated in Figure 11. When q 2 = −55 • and q 3 = 40 • (P max ), the maximum first order natural frequency is obtained, and the minimum appears when q 2 = −2 • and q 3 = −56 • (P min ).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 14 The two models constructed for two different robots through the method proposed in this paper both own good prediction ability. That is, the universality of the proposed method is testified.

Discussion
Once the prediction model of natural frequency is constructed and verified, the prediction result can be used as optimization parameter to improve the machining performance. An application case of the prediction model is explained as follows.
a2 and a3 control the two longest links, by which the basic configuration of robot is determined. Through VIP values, items involving q2 and q3 are demonstrated to act dominated role in prediction model. a4 and a5 decide the direction of robot flange, which affect the configuration less than a2 and a3. Taking M2 for an example, the influence of q2 and q3 on natural frequency is specifically studied and illustrated in Figure 11. When q2 = −55° and q3 = 40° (Pmax), the maximum first order natural frequency is obtained, and the minimum appears when q2 = −2° and q3 = −56° (Pmin).  Figure 12. The first order natural frequencies decrease in turn from P1 to P3 according to model M2. The setup of milling test is shown in Figure 13. High speed milling with short and straight path are executed, so that robot configuration can be ignored in one milling path. Milling parameters are as follows: the spindle rotation frequency is 800 Hz, f is 2.4 mm/s, ae is 1 mm, and ap is 4 mm. Acceleration signals and milled surface are analyzed. Three configurations near P max , P1 (q 2 = −48 • , q 3 = 50 • ), P2 (q 2 = −38 • , q 3 = 45 • ), and P3 (q 2 = −28 • , q 3 = 40 • ) are chosen for milling test (as shown in Figure 12. The first order natural frequencies decrease in turn from P1 to P3 according to model M 2 . The setup of milling test is shown in Figure 13. High speed milling with short and straight path are executed, so that robot configuration can be ignored in one milling path. Milling parameters are as follows: the spindle rotation frequency is 800 Hz, f is 2.4 mm/s, a e is 1 mm, and a p is 4 mm. Acceleration signals and milled surface are analyzed.  Milling acceleration signals are treated with short-time fourier transform (STFT), spectrograms of the range 5-100 Hz are displayed in Figure 14. For all three configurations, several peaks are conspicuous around 50-90 Hz during the whole milling process. In low frequency stage, a peak about 16 Hz can be found in P2 and a 15 Hz peak for P3, which are corresponding to their predicted natural frequencies. That no conspicuous can be found in P1, may because the certain frequency is not significantly impacted by the milling parameters. As for the milled surface, an obvious tool recession appears at P3, as shown in the red frame in Figure 15, while the situations are better in at P1 and P2. In addition, the quality of milled surfaces at P3 is the poorest compared with P1 and P2, and the vibration phenomenon on milled surface is getting severer from P1 to P3, indicating that configuration with higher first order frequency may lead to better milling performance, as mentioned in [20]. Whether there exits practical relevance between the value of natural frequency and machining performance can be a further research topic to develop more application of natural frequency prediction.  Milling acceleration signals are treated with short-time fourier transform (STFT), spectrograms of the range 5-100 Hz are displayed in Figure 14. For all three configurations, several peaks are conspicuous around 50-90 Hz during the whole milling process. In low frequency stage, a peak about 16 Hz can be found in P2 and a 15 Hz peak for P3, which are corresponding to their predicted natural frequencies. That no conspicuous can be found in P1, may because the certain frequency is not significantly impacted by the milling parameters. As for the milled surface, an obvious tool recession appears at P3, as shown in the red frame in Figure 15, while the situations are better in at P1 and P2. In addition, the quality of milled surfaces at P3 is the poorest compared with P1 and P2, and the vibration phenomenon on milled surface is getting severer from P1 to P3, indicating that configuration with higher first order frequency may lead to better milling performance, as mentioned in [20]. Whether there exits practical relevance between the value of natural frequency and machining performance can be a further research topic to develop more application of natural frequency prediction. Milling acceleration signals are treated with short-time fourier transform (STFT), spectrograms of the range 5-100 Hz are displayed in Figure 14. For all three configurations, several peaks are conspicuous around 50-90 Hz during the whole milling process. In low frequency stage, a peak about 16 Hz can be found in P2 and a 15 Hz peak for P3, which are corresponding to their predicted natural frequencies. That no conspicuous can be found in P1, may because the certain frequency is not significantly impacted by the milling parameters.  Milling acceleration signals are treated with short-time fourier transform (STFT), spectrograms of the range 5-100 Hz are displayed in Figure 14. For all three configurations, several peaks are conspicuous around 50-90 Hz during the whole milling process. In low frequency stage, a peak about 16 Hz can be found in P2 and a 15 Hz peak for P3, which are corresponding to their predicted natural frequencies. That no conspicuous can be found in P1, may because the certain frequency is not significantly impacted by the milling parameters. As for the milled surface, an obvious tool recession appears at P3, as shown in the red frame in Figure 15, while the situations are better in at P1 and P2. In addition, the quality of milled surfaces at P3 is the poorest compared with P1 and P2, and the vibration phenomenon on milled surface is getting severer from P1 to P3, indicating that configuration with higher first order frequency may lead to better milling performance, as mentioned in [20]. Whether there exits practical relevance between the value of natural frequency and machining performance can be a further research topic to develop more application of natural frequency prediction. As for the milled surface, an obvious tool recession appears at P3, as shown in the red frame in Figure 15, while the situations are better in at P1 and P2. In addition, the quality of milled surfaces at P3 is the poorest compared with P1 and P2, and the vibration phenomenon on milled surface is getting severer from P1 to P3, indicating that configuration with higher first order frequency may lead to better milling performance, as mentioned in [20]. Whether there exits practical relevance between the value of natural frequency and machining performance can be a further research topic to develop more application of natural frequency prediction. Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 14 Figure 15. Milled surfaces of P1-P3.

Conclusions
It is impossible to measure all the natural frequency information of the whole robot workspace, but the information is vital in improving robot milling performance. In view of that, a method to predict the natural frequency is proposed in this paper. The core content of the method proposed is to construct a prediction model with sample configurations, which takes joint angles as input variables. Two insignificant variables q1 and q6 are abandoned and workspace is limited for practicability. Considering the nonlinear influence of joint angles on natural frequency, CCD method is used to select sample configurations from limited variable range, and EMT is applied to acquire the sample data of the example robot. Several models are fitted through the sample data. The model validation procedure proves that the model fitted and optimized by PLS method has the best prediction ability. Then the method proposed is applied completely onto another robot for universality verification, and the prediction performance turns out to be outstanding. Thus, this method can be applied on any 6R machining industrial robot to evaluate its natural frequency distribution, which can be used for in machining performance improvement.

Conclusions
It is impossible to measure all the natural frequency information of the whole robot workspace, but the information is vital in improving robot milling performance. In view of that, a method to predict the natural frequency is proposed in this paper. The core content of the method proposed is to construct a prediction model with sample configurations, which takes joint angles as input variables. Two insignificant variables q 1 and q 6 are abandoned and workspace is limited for practicability. Considering the nonlinear influence of joint angles on natural frequency, CCD method is used to select sample configurations from limited variable range, and EMT is applied to acquire the sample data of the example robot. Several models are fitted through the sample data. The model validation procedure proves that the model fitted and optimized by PLS method has the best prediction ability. Then the method proposed is applied completely onto another robot for universality verification, and the prediction performance turns out to be outstanding. Thus, this method can be applied on any 6R machining industrial robot to evaluate its natural frequency distribution, which can be used for in machining performance improvement.