*Article* **Natural Frequency Prediction Method for 6R Machining Industrial Robot**

#### **Jiabin Sun 1,\*, Weimin Zhang <sup>1</sup> and Xinfeng Dong <sup>2</sup>**


Received: 10 October 2020; Accepted: 13 November 2020; Published: 17 November 2020

**Abstract:** 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 configuration, for improving dynamic performance during robot machining. A prediction model of natural frequency which expresses the mathematical relation between natural frequency and configuration is constructed for a 6R robot. Joint angles are used as input variables to represent the configurations in the model. The quantity and range of variables are limited for efficiency and practicability. Then sample configurations are selected by central composite design method due to its capacity of disposing nonlinear effects, and natural frequency data is acquired through experimental modal test. The model, which is in form of regression equation, is fitted and optimized with sample data through partial least square (PLS) method. The proposed model is verified with random configurations and compared with the original model and a model fitted by least square method. Prediction results indicate that the model fitted 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.

**Keywords:** machining robot; natural frequency prediction; model optimization; dynamic performance

#### **1. 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–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–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–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.

**Figure 1.** Frequency Response Functions of different robot configurations.

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.

**Figure 2.** Prediction method of robot natural frequency.

#### **2. Materials and Methods**

#### *2.1. 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

$$f\_{\mathcal{F}} = \mathbf{f}(q\_1, q\_2, \dots, q\_6),\tag{1}$$

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*<sup>1</sup> = *q*<sup>6</sup> = 0◦ in practice) for the simplification, and *q*2–*q*<sup>5</sup> are applied as variables. Prediction model can be adapted to

$$f\_{\mathcal{F}} = \mathbf{f}(q\_2, q\_3, q\_4, q\_5),\tag{2}$$

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.

**Figure 3.** Partial variable range.

#### 2.1.2. 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.


**Table 1.** Factors and Levels of central composite design (CCD).

#### 2.1.3. 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

$$h\_{\vec{i}\vec{j}}(\omega) \;=\sum\_{k=1}^{n} \left( \frac{a\_{\vec{i}\vec{j}k}}{(j\omega - p\_k)} + \frac{a\_{\vec{i}\vec{j}k}^\*}{(j\omega - p\_k^\*)} \right) \tag{3}$$

which presents the response of *i* point caused by the excitation at *j* point. *aijk* and *a*<sup>∗</sup> *ijk* are related to mode shapes, their values are depending on relevant combination of excitation and response. *pk* and *p*∗ *<sup>k</sup>* 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.

**Figure 4.** Directions of accelerometer and excitations.

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.

**Figure 5.** The setup of experimental modal test (EMT) system.

#### *2.2. 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.

#### 2.2.1. 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

$$\text{Cov}(\mathbf{t}, \mathfrak{u}) = \sqrt{\text{Var}(\mathbf{t}) \text{Var}(\mathfrak{u})} \mathbf{r}(\mathfrak{t}, \mathfrak{u}) \to \max,\tag{4}$$

The process of regression equation fitting by PLS method [27] is briefly introduced as follows.


Current principal component *t*<sup>h</sup> is validated by following indexes:

$$R\_Y^2(\text{cum}) \underset{\mathbf{h} = \mathbf{h}}{\text{ }}^\text{m} \mathbf{r}^2(\mathbf{y}, \mathbf{t}\_\mathbf{h})\text{.}\tag{5}$$

$$Q\_{\rm h}^2 = 1 - \frac{S\_{SPE,\rm h}}{S\_{SE,\rm h}},\tag{6}$$

$$Q\_{\rm h}^2(\text{cum}) = 1 - \prod\_{\mathbf{k}=1}^{\text{h}} \frac{S\_{\rm SE,k}}{S\_{\rm SE,k}} \, ^\prime \tag{7}$$

where *SSPE*,h is the square sum of prediction error of *<sup>Y</sup>*h−1, *SSE*,h is square sum of deviation of *<sup>Y</sup>*h−1.

*R*2 *<sup>Y</sup>*(cum) indicates the ability of principal component to explain the variability of *<sup>Y</sup>*h−1, and its upper limit value is 1. *Q*<sup>2</sup> <sup>h</sup> represents the contribution of *<sup>t</sup>*<sup>h</sup> to the model. When it is more than (1–0.95)<sup>2</sup> = 0.0975, *t*<sup>h</sup> has significant effect on prediction model. Meanwhile, *Q*<sup>2</sup> <sup>h</sup>(cum) increases, meaning that the comprehensive significance in explaining *Y* is enhanced.


$$\mathbf{Y}\_0 = \mathbf{t}\_1 \mathbf{r}\_1^\mathrm{T} + \mathbf{t}\_2 \mathbf{r}\_2^\mathrm{T} + \dots + \mathbf{t}\_\mathrm{Q} \mathbf{r}\_\mathrm{Q}^\mathrm{T} + \mathbf{Y}\_\mathrm{Q} \tag{8}$$

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

#### 2.2.2. 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

$$\text{VIP}\_{\text{j}} = \sqrt{\frac{p}{R\_{\text{Y}}^2 \text{(cum)}\_{\text{h}}}} \sum\_{\text{h}=1}^{\text{m}} \text{r}^2(\text{Y}, \text{t}\_{\text{h}}) w\_{\text{h}\text{'}}^2 \tag{9}$$

where *w*jh is the jth component of *w*h.

The larger VIPj is, the more important *x*<sup>j</sup> 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.

#### 2.2.3. 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*<sup>2</sup> <sup>h</sup>(cum), meaning that comprehensive ability is weakened. Though *R*<sup>2</sup> *<sup>Y</sup>*(cum) increases, the second principal component is still abandoned.

**Table 2.** Principal component validation of M.


#### 2.2.4. Prediction Model Optimization

*Q*2 <sup>h</sup>(cum) of M is merely 0.639, which means M is inadequate to represent *f* F. To optimize the model, VIP values of items in M is calculated and illustrated in Figure 6.

**Figure 6.** VIP values of items in M.

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 *q*2*q*3, because *q*2*q*<sup>3</sup> has the largest VIP value among cross items; meanwhile, *q*<sup>2</sup> and *q*<sup>3</sup> 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. *R*<sup>2</sup> *<sup>Y</sup>*(cum) and *<sup>Q</sup>*<sup>2</sup> <sup>h</sup>(cum) values of the six models are all improved obviously compared to M. M2, formed by adding *q*3*q*<sup>5</sup> to M1, has the best *R*<sup>2</sup> *<sup>Y</sup>*(cum) and *<sup>Q</sup>*<sup>2</sup> <sup>h</sup>(cum) values. M3–M6, which get lower *R*<sup>2</sup> *<sup>Y</sup>*(cum) and *<sup>Q</sup>*<sup>2</sup> <sup>h</sup>(cum) values compared to M1, are abandoned.

Based on M2, four new models are introduced by adding *q*2*q*4, *q*3*q*4, *q*2*q*5, and *q*4*q*<sup>5</sup> in turn. However, the largest *Q*<sup>2</sup> <sup>h</sup>(cum) only reaches up to 0.764. All the other combination groups of above four cross items are add to M2 for test, but none of them can acquire better validation result. Thus, the optimized model is finally confirmed to be M2. The coefficients of items are listed in Table 4 together with M.


**Table 3.** Validation results of M1–M6.


**Table 4.** Coefficients of items in M and M2.

#### **3. Results**

#### *3.1. Prediction Model Verification and Analysis*

#### 3.1.1. Model Verification with Random Configurations

Sixteen groups of *q*2–*q*<sup>5</sup> 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*<sup>2</sup> and *q*<sup>3</sup> 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.

**Figure 7.** Random configurations for model verification.

The first order natural frequencies of the 16 configurations are obtained by EMT. Then, original model M and optimized model M2 are utilized to predict the first order natural frequencies. The deviation between prediction value and measured value is defined as prediction error, the details are listed in Table 5; Table 6.


**Table 5.** Prediction error of M.


**Table 6.** Prediction error of M2.

The average error and average relative error of optimized model M2 are both much lower than those of original model M. The average relative error of M2 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 M2 has no predictions with over 5% relative error, and the maximum error is 0.51 Hz which is acceptable. Obviously, optimized model M2 has better prediction ability.

A model (MLS) 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 MLS are both higher than those of M and M2, 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.


**Table 7.** Prediction error of MLS.

#### 3.1.2. 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. (DMod*X*,N)*<sup>i</sup>* or (DMod*Y*,N)*<sup>i</sup>* is the ratio of construct quality of the *<sup>i</sup>*th sample point and the average construct quality. When the ratio is less than two, the construction for the *i*th sample point is reasonable. Construction quality of M2 is shown in Figure 8. (DMod*X*,N)*<sup>i</sup>* and (DMod*Y*,N)*<sup>i</sup>* are all less than two, that means the construction of M2 is rational.

#### *3.2. 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, *q*2*q*4, *q*2*q*5, and *q*4*q*<sup>5</sup> are abandoned according to VIP values. Validation results can be seen in Table 8, values of *R*<sup>2</sup> *<sup>Y</sup>*(cum) and *<sup>Q</sup>*<sup>2</sup> <sup>h</sup>(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.

**Figure 9.** ABB IRB6400 industrial robot.

**Table 8.** Validation results of M6400.


**Figure 10.** Construction quality of M6400.


**Table 9.** Prediction error of M6400.

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.

#### **4. 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*<sup>2</sup> and *a*<sup>3</sup> control the two longest links, by which the basic configuration of robot is determined. Through VIP values, items involving *q*<sup>2</sup> and *q*<sup>3</sup> are demonstrated to act dominated role in prediction model. *a*<sup>4</sup> and *a*<sup>5</sup> decide the direction of robot flange, which affect the configuration less than *a*<sup>2</sup> and *a*3. Taking M2 for an example, the influence of *q*<sup>2</sup> and *q*<sup>3</sup> on natural frequency is specifically studied and illustrated in Figure 11. When *q*<sup>2</sup> = −55◦ and *q*<sup>3</sup> = 40◦ (Pmax), the maximum first order natural frequency is obtained, and the minimum appears when *q*<sup>2</sup> = −2◦ and *q*<sup>3</sup> = −56◦ (Pmin).

**Figure 11.** Variation of natural frequency with *q*<sup>2</sup> and *q*3.

Three configurations near Pmax, P1 (*q*<sup>2</sup> = −48◦, *q*<sup>3</sup> = 50◦), P2 (*q*<sup>2</sup> = −38◦, *q*<sup>3</sup> = 45◦), and P3 (*q*<sup>2</sup> = −28◦, *q*<sup>3</sup> = 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 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.

**Figure 13.** Milling test system setup.

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.

**Figure 14.** Short-time fourier transform (STFT) of vibration signals of P1–P3.

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.

**Figure 15.** Milled surfaces of P1–P3.

#### **5. 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*<sup>1</sup> and *q*<sup>6</sup> 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.

**Author Contributions:** Conceptualization, J.S., W.Z. and X.D.; methodology, J.S. and X.D.; validation, J.S.; formal analysis, J.S.; investigation, J.S.; resources, W.Z.; data curation, J.S.; writing—original draft preparation, J.S.; writing—review and editing, J.S. and W.Z.; visualization, J.S.; supervision, D.X and W.Z.; project administration, W.Z.; funding acquisition, W.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is supported by the National key R&D project—The Construction, Reference Implementation, and Verification Platform of Reconfigurable Intelligent Production System, China (Grant No. 2017YFE0101400).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
