Modeling and Identiﬁcation of an Industrial Robot with a Selective Modal Approach

: The sti ﬀ ness properties of industrial robots are very important for many industrial applications, such as automatic robotic assembly and material removal processes (e.g., machining and deburring). On the one hand, in robotic assembly, joint compliance can be useful for compensating dimensional errors in the parts to be assembled; on the other hand, in material removal processes, a high Cartesian sti ﬀ ness of the end-e ﬀ ector is required. Moreover, low frequency chatter vibrations can be induced when low-sti ﬀ ness robots are used, with an impairment in the quality of the machined surface. In this paper, a compliant joint dynamic model of an industrial robot has been developed, in which joint sti ﬀ ness has been experimentally identiﬁed using a modal approach. First, a novel method to select the test conﬁgurations has been developed, so that in each conﬁguration the mode of vibration that chieﬂy involves only one joint is excited. Then, experimental tests are carried out in the selected conﬁgurations in order to identify joint sti ﬀ ness. Finally, the developed dynamic model of the robot is used to predict the variation of the natural frequencies in the workspace.


Introduction
In robotic processes, the compliance of the robot arm plays a very important role. In some conditions, for example, in robotic assembly, robot arm compliance can compensate for small position and orientation errors of the end-effector. In other processes, like machining, robot compliance may generate chatter vibrations, with an impairment in the quality of the machined surface. Indeed, still very few robots have been applied in this economically important application area [1], mainly due to their low stiffness [2]. In industrial robots, the compliance of the end-effector is mainly due to joint compliances [3][4][5], even if there are some examples of robots [6] having structural modes (dominated by link or bearing compliance) in the band of frequency that contains the modes dominated by joint compliance. The compliance properties of robots are very important in emerging fields as well, such as flexible assembly systems [7,8] and collaborative robotics [9]. Once robot joint stiffness has been identified, the kinematic redundancy of the arm with respect to the robotic task can be exploited to minimize (or maximize) the Cartesian compliance of the end-effector through optimization methods [10][11][12][13], as already proposed for the minimization of base reactions and kinetic energy.
In material removal processes, low stiffness causes imprecise products, due to the robot deflections during the robotic task. From a dynamic point of view, low frequency chatter vibrations [3,14,15] can be induced when low-stiffness robots are used, with an impairment in the quality of the machined surface. Moreover, vibrations cause a reduction of tool life and can lead to breakages in robot joint transmission. In particular, a low joint stiffness of the robot is one of the main issues in using robotic machining instead of computer numerical control (CNC) machining.
In robotic processes, the directions along which the robot arm shows large compliance are important. The stiffness matrix in the joint space does not directly give this piece of information; moreover, for a serial robot, the stiffness matrix in the Cartesian space is not diagonal and it is configuration-dependent. This means that the force and deformation in the Cartesian space are coupled and this can generate some counterintuitive results. A geometric and intuitive interpretation of the end-effector displacements in different robot operative conditions has been presented in [16].
Static tests are widely used to obtain the joint stiffness of industrial robots [17]. In the static tests, a set of forces is applied to the robot end-effector in different robot configurations, while the displacement sensors (laser sensors, vision systems, coordinate measuring machines) measure the static deformation of the end-effector. Therefore, the Cartesian stiffness of the robot can be calculated. Then, the stiffness of the joints is obtained through the analytical relation between the joint and Cartesian stiffness based on the kinematic model of robot. Various identification methods (least squares or genetic algorithms) are employed to numerically solve the problem.
Most of the researches (e.g., see [4,[18][19][20][21]) neglect link flexibility and identify the joints' rotational stiffness (modeled with linear rotational springs) using the above-mentioned method. In [22], the robot arm is modeled considering rigid links and three lumped rotational springs for each joint to take into account joint compliance, bearings compliance, and link deflections. In [23], an analysis method in which both joint and link stiffness are considered is presented. In [24,25], two types of robot dynamic models are presented by considering either only joint compliance or its combination with link compliance, in order to predict the robot dynamic behavior in different postures. In [6] a hybrid approach is adopted, in which a critical structural compliance in a direction perpendicular to a joint axis is taken into account, together with joint compliance.
In [26], experimental modal analysis is used to identify the joint and base stiffness of an industrial robot represented with a four-degrees-of-freedom (four-DOF) planar model. Experimental modal analysis was also used in [6,16], to identify the stiffness of the first three joints of an industrial robot.
The static tests for the identification of joint stiffness have some drawbacks. The equipment for the measurement of end-effector displacements along three directions may be complex and expensive, especially if high accuracy is required and the workspace of the robot is large. The application of vertical forces on the end-effector is rather simple, since weights can be used, but the application of horizontal forces may require specific actuators that have to be connected with the end-effector without introducing additional constraints. Since Cartesian stiffness is configuration dependent, many test configurations are needed to map the robot workspace. The identification of joint stiffness from the measured Cartesian stiffness requires, in most cases, an optimization method. Finally, the measurements are affected by the structural compliance of the links, and it is difficult to separate the contribution of end-effector displacement caused by joint compliance from the one due to link compliance.
For the above-mentioned reasons, this research aims to use the methods of impulsive modal analysis to identify the joint stiffness of industrial robot arms. In future, the same methodology could be extended to other kinds of robots e.g., walking robots [27]. Some factors foster the development of testing methods based on impulsive excitation. The basic equipment (which includes an instrumented hammer, an accelerometer, a data acquisition board, and a computer) is nowadays affordable for most laboratories. Hammer excitation can be exerted in three directions, even if the robot workspace is wide. Modal analysis identifies all the modes of vibration, both the ones dominated by joint compliance, and the ones dominated by link compliance (if present), hence, the modal method allows a better understanding of the importance of link compliance than the static method.
The modes of vibration of a robot typically involve the rotations of some joints and the natural frequencies are functions of the stiffness of these joints. Therefore, the main issue of the modal methods for the identification of joint stiffness is the possibility of finding robot configurations with modes of vibration dominated by the stiffness of only one joint. In this specific case, modal stiffness coincides with joint stiffness, and the latter can be identified from the measured natural frequency and the calculated value of moment of inertia. Other important issues are the choice of proper excitation directions for the selected configurations and the minimization of test configurations. This paper is organized as follows. In Section 2, a serial six-DOF robot is considered (Omron Adept Viper s650), since it is representative of industrial robots used for a wide range of operations, and a dynamic model of robot vibrations is developed. Section 3 deals with selective modal testing, and some criteria for finding robot configurations with modes of vibration dominated by only one joint, or by a couple of joints, are presented and discussed. Experimental tests carried out in the selected configurations are presented in Section 4, which significantly extends a previous work of the authors [16], in which only the first three joints of the robot were considered. Section 5 deals with the identification of joint stiffness, and an optimization technique is presented to cope with the cases in which the identified mode is dominated by a couple of joints. The use of the robot dynamic model with the identified stiffness values is presented in Section 6, and the variation of the natural frequencies in the robot workspace is analyzed. Finally, in Section 7, conclusions are drawn.

Dynamic Model
In the framework of this research an Adept Viper s650 robot manufactured by Omron Adept Inc. (Pleasanton, CA, USA) is considered. It is a medium size six-DOF industrial robot, its workable space is roughly a hemisphere with a radius of about 600 mm. Figure 1 shows that this robot has a typical anthropomorphic structure. There are three main joints that define the position of the wrist and three minor joints that define the orientation of the end-effector. All joints are driven by alternate current (AC) servomotors. The modes of vibration of a robot typically involve the rotations of some joints and the natural frequencies are functions of the stiffness of these joints. Therefore, the main issue of the modal methods for the identification of joint stiffness is the possibility of finding robot configurations with modes of vibration dominated by the stiffness of only one joint. In this specific case, modal stiffness coincides with joint stiffness, and the latter can be identified from the measured natural frequency and the calculated value of moment of inertia. Other important issues are the choice of proper excitation directions for the selected configurations and the minimization of test configurations. This paper is organized as follows. In Section 2, a serial six-DOF robot is considered (Omron Adept Viper s650), since it is representative of industrial robots used for a wide range of operations, and a dynamic model of robot vibrations is developed. Section 3 deals with selective modal testing, and some criteria for finding robot configurations with modes of vibration dominated by only one joint, or by a couple of joints, are presented and discussed. Experimental tests carried out in the selected configurations are presented in Section 4, which significantly extends a previous work of the authors [16], in which only the first three joints of the robot were considered. Section 5 deals with the identification of joint stiffness, and an optimization technique is presented to cope with the cases in which the identified mode is dominated by a couple of joints. The use of the robot dynamic model with the identified stiffness values is presented in Section 6, and the variation of the natural frequencies in the robot workspace is analyzed. Finally, in Section 7, conclusions are drawn.

Dynamic Model
In the framework of this research an Adept Viper s650 robot manufactured by Omron Adept Inc. (Pleasanton, CA, USA) is considered. It is a medium size six-DOF industrial robot, its workable space is roughly a hemisphere with a radius of about 600 mm. Figure 1 shows that this robot has a typical anthropomorphic structure. There are three main joints that define the position of the wrist and three minor joints that define the orientation of the end-effector. All joints are driven by alternate current (AC) servomotors. The vertical plane passing through the axis of joint 1, perpendicular to the axes of joint 2 and 3, is named the meridian plane of the robot arm. The vertical plane passing through the axis of joint 1, perpendicular to the axes of joint 2 and 3, is named the meridian plane of the robot arm. This paper aims at developing a dynamic model for the study of vibrations around a working configuration of the robot defined by vector q 0 of joint variables, therefore joint variables, velocities, and accelerations are defined by the following equations: where vector ∆q contains the small variations in joint variables around the selected configuration.
Joint compliance around the rotation axis is considered, whereas link and bearing compliance is neglected.
The axis of joint 6 coincides with the approach axis of the end-effector. Since, in most cases, the external force exerted on the end-effector is aligned with or intersects the approach axis, compliance about joint 6 is neglected.
If the non-linear Coriolis and centrifugal terms [28] are neglected, according to the assumption of small oscillations, the equations of free vibrations take the form [6]: where M q 0 is the mass matrix that depends on q 0 . The elements on the diagonal of M q 0 are the direct inertia terms, the off-diagonal terms are the inertial cross-coupling terms. C q is a diagonal matrix that accounts for joint damping, K is a stiffness matrix that accounts for joint compliance and configuration-depending gravity torques. Gravity can generate large restoring or unbalancing torques about joint axes only when the robot arm is close to a vertical configuration (upwards or downwards).
If configuration-depending gravity torques are negligible the equations of free vibrations become: where K q is a diagonal matrix that accounts for joint compliance. Specific calculations have been be carried out to estimate the effect of gravity torques in the selected configurations. The identification of matrices K q and C q is the main task of this research. The mass matrix M q 0 can be calculated using the inertial data provided by the computer aided design (CAD) models of the robot. These CAD files have been retrieved from the website of the manufacturer and the links are assumed to have uniform density and no free spaces within the external shell. This is due to the uncertainty of the mass distribution within the links of the robot that could not be eliminated without a complete disassembly of the robot.
Starting from the CAD files, the following inertial parameters have been retrieved: • Mass of link i (m i ) has been calculated by distributing the total mass m tot based on the volume of each link; • Center of mass of link i (G i ) has been calculated with respect to the corresponding Denavit-Hartenberg reference frame (see Figure 2); • Inertia tensor of link i (I i ) has been calculated at the center of mass aligned with the correspondent Denavit-Hartenberg reference frame [28].

Selective Modal Testing
Vibration signals contain information about mechanical systems [29], and this paper aims to use impulsive modal analysis to identify the stiffness properties of an industrial robot, and to develop a compliant robot model.
Modal analysis is a mature technology for the study of vibrations [30], which has been successfully used in many fields of aerospace [31], vehicle, and automation engineering [16,32,33].
In the field of industrial robotics, it is well recognized that joint stiffness has a larger effect on robot performance than link stiffness, therefore the main task of impulsive modal analysis is the identification of the stiffness of the joints.
A robot with assigned joint values ( ) is considered. When an impulse force is exerted on a point of the robot by means of a hammer for modal testing and where the mass, damping and stiffness matrices are diagonal, and elements * ), * ) and * ) represent modal mass, damping and stiffness associated to the th mode.
When the robot is tested in a generic configuration the modes of vibration typically involve rotations about many joints of the robot. Therefore, modal coordinate is a combination of joint On the basis of the above-mentioned parameters mass matrix M q 0 can be calculated as follows: where R i is the rotation matrix from link i to the base frame and J p,i and J O,i are the Jacobians:

Selective Modal Testing
Vibration signals contain information about mechanical systems [29], and this paper aims to use impulsive modal analysis to identify the stiffness properties of an industrial robot, and to develop a compliant robot model.
Modal analysis is a mature technology for the study of vibrations [30], which has been successfully used in many fields of aerospace [31], vehicle, and automation engineering [16,32,33].
In the field of industrial robotics, it is well recognized that joint stiffness has a larger effect on robot performance than link stiffness, therefore the main task of impulsive modal analysis is the identification of the stiffness of the joints.
A robot with assigned joint values (q) is considered. When an impulse force is exerted on a point of the robot by means of a hammer for modal testing and N acceleration components are measured by means of accelerometers, a set of N modes of vibrations can be identified. These modes of vibration transform the equations of motions (3) into a set of N independent equations in the modal coordinates η k k = 1, . . . , N: where the mass, damping and stiffness matrices are diagonal, and elements M * kk (q), C * kk (q) and K * kk (q) represent modal mass, damping and stiffness associated to the kth mode. When the robot is tested in a generic configuration the modes of vibration typically involve rotations about many joints of the robot. Therefore, modal coordinate η k is a combination of joint rotations, which is defined by the corresponding mode of vibration. Modal mass, damping and stiffness do not have a direct physical meaning, and are not very useful to identify the stiffness and damping characteristics of the robot joints.
Conversely, when the same robot is modally tested in a specific configuration in which the mass matrix shows that a joint variable i has negligible inertial cross coupling terms, one of the modes of vibration is dominated by the rotation about joint i, and the corresponding modal equation becomes: The modal coordinate η i coincides with joint variable q i , modal mass coincides with inertia about joint i (I i zz (q)), modal stiffness and damping coincide with joint stiffness (k i ) and damping (c i ), respectively. Mode i is very useful to identifying the stiffness of joint i.
Inertia about joint i (I i zz (q)) depends on the inertias of the following links and can be calculated by using the Huygens-Steiner theorem: where R ji is the rotation matrix that transforms the reference frame of link j to the reference frame of link i, and x l ji , y l ji and z l ji are the coordinates of the center of mass of link j defined in the reference frame of link i. It is worth noticing that usually mode i (dominated by joint i) includes small contributions of the other joints, and sometimes of structural compliance of the links, therefore assuming M * ii (q) = I i zz (q), K * ii (q) = k i and C * ii (q) = c i is only an approximation.
For the above-mentioned reasons, a specific modal testing method is proposed, which is based on robot configurations and excitation directions that make possible the excitation of modes of vibrations dominated by only one joint.
In order to have an insight into the properties of the mass matrix of the tested robot, the elements of the mass matrix were evaluated in random configurations. In particular 1,000,000 random combinations of q 1 . . . q 5 with values within the actual joint ranges were considered.
Since the inertial cross-coupling terms may be positive or negative, depending on the configuration, the mean values and standard deviations of the moduli of the elements of the mass matrix were calculated to avoid cancellations.
The results of this analysis are summarized in Figure 3, and show the following properties.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 19 Figure 3. Mean values and standard deviations of mass matrix terms (10 random configurations), direct inertia terms (red bars) and inertial cross-coupling terms (blue bars).
In a specific configuration, the inertial cross-coupling terms of joint ( , ≠ ) are small or negligible if the velocity of joint generates a distribution of robot arm velocities that cannot be generated by the velocities of the other joints. Even if the previous condition is satisfied, cross-coupling terms may be present, if the axes of the joints are not principal axes.
According to the results of Figure 3, the configurations for selective modal testing summarized in Figure 4 and Table 1 were chosen.  (a) The direct inertia terms of the first three joints are by far larger than the ones of the wrist joints-notice that a logarithmic scale is adopted in Figure 3. Joint 2 (J2) has the largest mean value followed by Joint 1 (J1). (b) The direct inertia term of joint i (M ii ) depends on the joint variables of the following links, hence M 11 depends on q 2 . . . q N , whereas M NN is constant. Therefore, the standard deviation decreases from J1 to J5. (c) The inertial cross-coupling terms between the first three joints are much smaller than the direct inertia terms with the exception of M 23 , which is comparable with M 33 . (d) M 23 has a large standard deviation, and this means that there are robot configurations with a minimum coupling between J2 and J3. (e) The mean values of the inertial cross-coupling terms of the first three joints with the wrist joints are very small in comparison with the direct inertia terms of the first three joints. They are more important if compared with the direct inertia terms of the wrist joints (f) There is a negligible coupling between wrist joints (M 45 ).
In a specific configuration, the inertial cross-coupling terms of joint i (M ij , j i) are small or negligible if the velocity of joint i generates a distribution of robot arm velocities that cannot be generated by the velocities of the other joints. Even if the previous condition is satisfied, cross-coupling terms may be present, if the axes of the joints are not principal axes.
According to the results of Figure 3, the configurations for selective modal testing summarized in Figure 4 and Table 1 were chosen.
In a specific configuration, the inertial cross-coupling terms of joint ( , ≠ ) are small or negligible if the velocity of joint generates a distribution of robot arm velocities that cannot be generated by the velocities of the other joints. Even if the previous condition is satisfied, cross-coupling terms may be present, if the axes of the joints are not principal axes.
According to the results of Figure 3, the configurations for selective modal testing summarized in Figure 4 and Table 1 were chosen.   Joint 1 (J1) is always perpendicular to joints 2 and 3 (J2 and J3). If joint variables q 2 , q 3 are set to 0 • and 90 • , respectively, the robot arm is horizontal, and if q 4 q 5 are set to 0 • , the approach axis is horizontal as well. In this configuration (Test 1) only J1 is able to generate relevant velocities in the horizontal plane. This configuration minimizes the inertial coupling between J1 and the other joints (which is usually small) and enhances the excitation of J1, when a lateral force (in y direction) is exerted on the end-effector.
Joints J2 and J3 are always parallel, and generate velocity components in the same plane, Figure 3 shows that the inertial cross coupling may be large. The mass matrix of a robot with two parallel joints [28] shows an inertial cross-coupling term depending on joint variable q 3 (the angle between the links). A numerical calculation of the ratios between the inertial cross-coupling term M 23 and the direct inertias of J2 and J3 (M 22 and M 33 respectively) is depicted in Figure 5. They show that the inertial cross-coupling term tends to decrease when the two links tend to overlap, notice that the two links are aligned when q 3 = 90 • . Minimum coupling configurations can be found for q 3 = −38 • and q 3 = 293 • . These theoretical values have to be compared with the actual range of the joint (−29 • ≤ q 3 ≤ 256 • ), and, for this reason, a configuration with q 3 = 243 • (Test 3) is suited to identifying the stiffness both of J2 and J3. It is worth noticing that, in this configuration, a vertical force generates a large moment about J3, thus, it is able to excite the mode dominated by this joint. Conversely, in this configuration a vertical force produces a small moment about J2, thus, the excitation of the mode dominated by this joint is poor. Figures 3 and 5 show that the inertial cross-coupling term is less important for J2 than for J3, since the direct inertia M 22 term is much larger than M 23 . For this reason, configuration Test 2 of Table 1, in which links 2 and 3 are almost perpendicular, gives a good compromise between the reduction of inertial cross coupling and the generation of a large moment about J2. large moment about J3, thus, it is able to excite the mode dominated by this joint. Conversely, in this configuration a vertical force produces a small moment about J2, thus, the excitation of the mode dominated by this joint is poor. Figures 3 and 5 show that the inertial cross-coupling term is less important for J2 than for J3, since the direct inertia 22 term is much larger than 23 . For this reason, configuration Test 2 of Table 1, in which links 2 and 3 are almost perpendicular, gives a good compromise between the reduction of inertial cross coupling and the generation of a large moment about J2.  Joint J4 is able to generate velocity components perpendicular to the vertical plane that contains the robot arm (meridian plane). J2, J3 generate velocity components in meridian plane, and the inertial cross coupling with J4 is small. Conversely, J1 is able to generate velocity components in the same direction as J4, and there is a significant cross coupling with J4 (see Figure 3). The inertial cross coupling between J4 and J5 is always negligible (see Figure 3). In order to enhance the moment that a lateral hit can generate about J4, q 5 was set to ±90 • and configuration Test 4 was defined.
Finally, when q 4 is set to zero, Joint 5 (J5) is able to generate velocity components of the last links in the meridian plane of the robot arm, but these components can be generated by J2 and J3 as well. When q 4 is set to 90 • , J5 is able to generate velocity components out of the meridian plane of the robot arm, but these components can be generated by J1 as well. Therefore, Figure 3 shows that there is always a significant coupling between J5 and the arm joints. The test configuration can be chosen only with the aim of increasing the moment about J5 generated by the hammer hit. Test configuration Test 5 satisfies this condition for a hammer hit in the longitudinal direction.
The relevance of gravity restoring (or unbalancing) torques was evaluated in the selected configurations. The calculated values resulted much smaller than reference values of joint torques caused by joint stiffness [16].

Experimental Tests
The equipment needed to carry out the modal analysis of a robot with the impulsive method is rather simple and cheap. It includes a hammer for modal testing, an accelerometer, and a data acquisition board. In the framework of this research, a PCB 086C01 hammer (with load cell sensitivity 0.2549 mV/N), a PCB 356A17 tri-axial accelerometer (sensitivity 50.5 mV/(m/s 2 )), and a NI9234 data acquisition board were used. The hammer and the accelerometer were made by PCB Piezotronics Inc., Depew, NY, USA, whereas the board was made by National Instruments Inc., Austin, TX, USA. Measured signals were analyzed in the frequency domain by means of ModalVIEW R2 (2017), which is a specific software for modal analysis developed by ABSignal-National Instruments Inc., Austin, TX, USA.
In order to characterize the vibrations of the robot 10 testing points were defined. Figure 6 shows that testing point 1 is located on the flange of the end-effector, testing point 10 is located on the base, before joint 1, and the other testing points are evenly distributed on the links. The number of testing points is a trade-off between the need for a large number of measurement locations, to accurately describe robot vibrations, and the need for a quick testing procedure. Modal analysis was carried out with the rowing response approach [34] Excitation was always exerted on the end-effector, whereas the tri-axial accelerometer was sequentially moved to the 10 testing points. The three axes of the accelerometer were always parallel or perpendicular to the robot surface.
In order to characterize the vibrations of the robot 10 testing points were defined. Figure 6 shows that testing point 1 is located on the flange of the end-effector, testing point 10 is located on the base, before joint 1, and the other testing points are evenly distributed on the links. The number of testing points is a trade-off between the need for a large number of measurement locations, to accurately describe robot vibrations, and the need for a quick testing procedure. Modal analysis was carried out with the rowing response approach [34] Excitation was always exerted on the end-effector, whereas the tri-axial accelerometer was sequentially moved to the 10 testing points. The three axes of the accelerometer were always parallel or perpendicular to the robot surface.  Since motions of the base of the robot may have a negative effect on the quality of measurements, the robot was rigidly fastened to a large steel base.
After some preliminary tests, data acquisition was carried out with a sampling frequency of 1024 Hz and 2048 samples. Since motions of the base of the robot may have a negative effect on the quality of measurements, the robot was rigidly fastened to a large steel base.
After some preliminary tests, data acquisition was carried out with a sampling frequency of 1024 Hz and 2048 samples.   Figure 8 shows the experimental results obtained in configuration Test 2 with hammer excitation in the vertical direction (z). In this case, the direct point FRF and the CMIF show a high and isolated peak at about 17 Hz. ModalVIEW made it possible to identify at 17.5 Hz a mode of vibration dominated by the rotation of the whole robot arm about J2, which is represented in Figure  8. The corresponding damping ratio is 1.9%, this value is congruent with the damping ratio of J1.  Figure 8 shows the experimental results obtained in configuration Test 2 with hammer excitation in the vertical direction (z). In this case, the direct point FRF and the CMIF show a high and isolated peak at about 17 Hz. ModalVIEW made it possible to identify at 17.5 Hz a mode of vibration dominated by the rotation of the whole robot arm about J2, which is represented in Figure 8. The corresponding damping ratio is 1.9%, this value is congruent with the damping ratio of J1. The experimental results obtained in configuration Test 3 are summarized in Figure 9. The direct point FRF shows a high peak at about 30 Hz with a large phase change, but in this case the main peak is not alone, since there is another peak at about 35 Hz. The CMIF plot confirms this scenario. ModalVIEW made it possible to identify two modes. The first mode is at 30.8 Hz with a damping ratio of 2.1%. This mode, which is represented in Figure 9, is dominated by the compliance of J3, and shows some contribution of J2, owing to the inertial-cross coupling between the two joints. The second mode at 34.7 Hz is influenced by the structural compliance. In particular, structural compliance in the direction perpendicular to the axes J2 and J3 leads to an apparent torsion of the robot arm. The large damping ratio (4.4%) is in agreement with the contribution of structural deformation. The presence of structural modes in the range of frequencies of the modes dominated by joint compliance has been detected by other researchers [6]. Figure 10 deals with experimental results obtained in configuration Test 4 with lateral excitation (y direction). Both the direct point FRF and the CMIF highlight the presence of an important resonance peak beyond the 10÷40 Hz frequency band, which contains the previously identified modes. Modal analysis made it possible to identify a mode of vibration at 56.8 Hz with a 3.0% damping ratio, which is similar to the ones of the other joints. This mode, which is represented in Figure 10, shows a large contribution of the rotation about J4, which leads to lateral displacements of the points of the following links.
Finally, Figure 11 shows experimental results obtained in configuration Test 5 with longitudinal excitation (x direction). At high frequency, there is a very important peak at about 170 Hz. Modal analysis showed that at 171.4 Hz there is a mode of vibration with large displacements localized in the last link and caused by the compliance of J5 (see Figure 11). The damping ratio of this mode is 2.8%. The experimental results obtained in configuration Test 3 are summarized in Figure 9. The direct point FRF shows a high peak at about 30 Hz with a large phase change, but in this case the main peak is not alone, since there is another peak at about 35 Hz. The CMIF plot confirms this scenario. ModalVIEW made it possible to identify two modes. The first mode is at 30.8 Hz with a damping ratio of 2.1%. This mode, which is represented in Figure 9, is dominated by the compliance of J3, and shows some contribution of J2, owing to the inertial-cross coupling between the two joints. The second mode at 34.7 Hz is influenced by the structural compliance. In particular, structural compliance in the direction perpendicular to the axes J2 and J3 leads to an apparent torsion of the robot arm. The large damping ratio (4.4%) is in agreement with the contribution of structural deformation. The presence of structural modes in the range of frequencies of the modes dominated by joint compliance has been detected by other researchers [6]. Figure 10 deals with experimental results obtained in configuration Test 4 with lateral excitation (y direction). Both the direct point FRF and the CMIF highlight the presence of an important resonance peak beyond the 10 ÷ 40 Hz frequency band, which contains the previously identified modes. Modal analysis made it possible to identify a mode of vibration at 56.8 Hz with a 3.0% damping ratio, which is similar to the ones of the other joints. This mode, which is represented in Figure 10, shows a large contribution of the rotation about J4, which leads to lateral displacements of the points of the following links.     Finally, Figure 11 shows experimental results obtained in configuration Test 5 with longitudinal excitation (x direction). At high frequency, there is a very important peak at about 170 Hz. Modal analysis showed that at 171.4 Hz there is a mode of vibration with large displacements localized in the last link and caused by the compliance of J5 (see Figure 11). The damping ratio of this mode is 2.8%.

Identification
Testing configurations Test 1 and Test 2 made it possible to identify modes of vibration clearly dominated by the stiffness properties of J1 and J2 respectively, and these stiffness properties are the most important from the point of view of robot operations (e.g., machining or assembly). Testing configurations Test 4 and Test 5 made it possible to identify modes of vibration with a large contribution of the stiffness of J4 and J5. Owing to the small moments of end-effector forces about J4 and J5, the stiffness of wrist joints is less critical than the stiffness of arm joints. Testing configuration Test 3 minimized the coupling between J2 and J3, nevertheless the identified modal shape a at 30.8 Hz showed some influence of the rotation about J2 and of structural deformation, since in this configuration the robot has a structural mode ad 34.7 Hz. Therefore, the measured frequency is chiefly influenced by the stiffness of J3, but it includes the contribution of other stiffness properties.
For the aforementioned reasons, the stiffness about J1, J2, J4 and J5 was directly calculated from equation 7, assuming modal coordinate coincident with and using a 1 DOF model: where is the natural frequency measured in testing configuration . The damping coefficient of joint ( ) was calculated with the following equation: where is the identified damping ratio of mode . The results of these calculations are summarized in Table 2.

Identification
Testing configurations Test 1 and Test 2 made it possible to identify modes of vibration clearly dominated by the stiffness properties of J1 and J2 respectively, and these stiffness properties are the most important from the point of view of robot operations (e.g., machining or assembly). Testing configurations Test 4 and Test 5 made it possible to identify modes of vibration with a large contribution of the stiffness of J4 and J5. Owing to the small moments of end-effector forces about J4 and J5, the stiffness of wrist joints is less critical than the stiffness of arm joints. Testing configuration Test 3 minimized the coupling between J2 and J3, nevertheless the identified modal shape a at 30.8 Hz showed some influence of the rotation about J2 and of structural deformation, since in this configuration the robot has a structural mode ad 34.7 Hz. Therefore, the measured frequency is chiefly influenced by the stiffness of J3, but it includes the contribution of other stiffness properties.
For the aforementioned reasons, the stiffness about J1, J2, J4 and J5 was directly calculated from Equation (7), assuming modal coordinate η i coincident with q i and using a 1 DOF model: where f ni is the natural frequency measured in testing configuration q i . The damping coefficient of joint i (c i ) was calculated with the following equation: where ζ i is the identified damping ratio of mode i. The results of these calculations are summarized in Table 2. For joint J3, the value given by Equation (11) was considered only a first estimate of the actual joint stiffness (k 3 ). The identification of this stiffness was then improved carrying out an optimization. The following penalty function was defined: This is the sum of the squared differences between the measured natural frequencies and the natural frequencies f pjk predicted by the model with assigned k 3 . Five modes of vibration (j = 1, . . . , 5) in the five testing configurations are considered (k = 1, . . . , 5). Minimization was carried out with the function fmincon of MATLAB. The optimized stiffness k 3 = 6485 Nm/rad was rather different from the first attempt value. Finally, damping coefficient was calculated according to Equation (12) and resulted c 3 = 1.41 Nms/rad.
The identified stiffness and damping properties were the implemented in the mathematical model of the robot (Equation (3)), and the natural frequencies and damping ratios were calculated in some validation configurations that were experimentally tested as well. Joint variables of validation configurations are summarized in Table 3.  Table 4 shows the comparison between numerical values and experimental values in the validation configurations. In validation configuration 1, the predicted natural frequencies of the first four modes are a bit lower than the experimental natural frequencies. The largest error is less than −10%, and it takes place in the second mode, which is a mode on the meridian plane dominated by the rotation about J2. The predicted natural frequency of Mode 5 is a bit larger than the experimental value.
In validation configuration 2, the predicted natural frequencies of Mode 1 and Mode 4 are a bit lower than the experimental values (about −3%), whereas the other natural frequencies are higher than the experimental values. The largest error (+15%) takes place in the frequency of Mode 3. This happens because modal analysis (Figure 9) showed that in this range of frequencies structural deformability affects the modes of vibration, and this effect is not taken into account by the numerical model, which is based on joint stiffness.
Looking at modal damping, there is a general agreement between numerical and experimental values in both configurations. The modes that show a large experimental damping typically show a large numerical damping. Only in the first mode the calculated value is significantly lower than the experimental value. The presence of configuration dependent friction forces, due to large gravity loads, may be the cause of this phenomenon.

Numerical Simulations
Recent studies [2,6,35,36] have highlighted that, in some robot operations (e.g., milling tasks), the dynamic properties of the robot arm are important, and that a simple optimization of robot configurations based on static stiffness does not lead to the best performance [2]. The modes of vibration of a robot are configuration-dependent, and the variation of natural frequency and modal damping in the workspace is an important feature of robot dynamics. It gives some hints useful to improve robot performance selecting robot configurations with highly damped modes having natural frequencies far from the excitation frequencies.
Experimental modal analysis is time consuming, thus, the mathematical model based on modal analysis in a small number of selected configurations is a very useful tool to predict the variation in the dynamic properties of the robot over the whole workspace.
The first numerical analysis aimed to understand the effect of robot configuration on the natural frequency and damping of the modes of vibration. A total of 10,000 random configurations of the robot corresponding to 10,000 random combinations of joint variables J2, J3, J4, and J5 (within joint ranges) were defined, and the modal properties were calculated. Since the robot model is axial-symmetric about J1, the variation of J1 was not taken into account, and this joint variable was set to zero. The results of these calculations are collected in Figures 12 and 13, which depict the mean values and standard deviations of natural frequencies and modal dampings, respectively.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 19 Looking at modal damping, there is a general agreement between numerical and experimental values in both configurations. The modes that show a large experimental damping typically show a large numerical damping. Only in the first mode the calculated value is significantly lower than the experimental value. The presence of configuration dependent friction forces, due to large gravity loads, may be the cause of this phenomenon.

Numerical simulations
Recent studies [2,6,35,36] have highlighted that, in some robot operations (e.g., milling tasks), the dynamic properties of the robot arm are important, and that a simple optimization of robot configurations based on static stiffness does not lead to the best performance [2]. The modes of vibration of a robot are configuration-dependent, and the variation of natural frequency and modal damping in the workspace is an important feature of robot dynamics. It gives some hints useful to improve robot performance selecting robot configurations with highly damped modes having natural frequencies far from the excitation frequencies.
Experimental modal analysis is time consuming, thus, the mathematical model based on modal analysis in a small number of selected configurations is a very useful tool to predict the variation in the dynamic properties of the robot over the whole workspace.
The first numerical analysis aimed to understand the effect of robot configuration on the natural frequency and damping of the modes of vibration. A total of 10,000 random configurations of the robot corresponding to 10,000 random combinations of joint variables J2, J3, J4, and J5 (within joint ranges) were defined, and the modal properties were calculated. Since the robot model is axial-symmetric about J1, the variation of J1 was not taken into account, and this joint variable was set to zero. The results of these calculations are collected in Figure 12 and 13, which depict the mean values and standard deviations of natural frequencies and modal dampings, respectively.
The natural frequencies of the second and third mode, which are usually characterized by large displacements in the meridian plane, show the largest standard deviations. Therefore, these modes are influenced by the robot configuration. The first natural frequency, which usually corresponds to a mode of vibration with large displacements in the horizontal plane, shows a mild standard deviation, since the robot configuration affects the moment of inertia about J1. The natural frequencies of the last two modes, which are chiefly dominated by J4 and J5, show very small standard deviations. This happens because the robot configuration has a small effect on the direct inertia term of J4 and no effect on the direct inertia term of J5, but it only changes the inertial cross-coupling terms (see Figure 3).   The statistical analysis of modal damping variation with robot configurations is depicted in Figure 13. This figure shows that robot configuration affects the modal damping of the first three modes, since standard deviations are comparable with the mean values. Conversely, robot configuration has a very small effect on the damping of the last modes, which chiefly involve wrist joints, because standard deviations are very small. Figure 12b and Figure 13b show the result of a similar calculation, in which only random variations in J2 and J3 were performed. The standard deviations of the natural frequencies and modal dampings are very similar to the ones of Figure 12a and Figure 13a, therefore, this result points out that the variations in J2 and J3 have the largest effect on the modal properties of the robot. Figures 12 and 13 highlight another important feature of the robot. Variations in joint variables are able to modify the natural frequencies, but they do not change the typical frequency bands of the modes. In other words, the frequencies of the modes associated to the wrist joints are always higher than the frequencies associated to the arm joints. The first two modes belong to the same frequency band (15 ÷ 30 Hz), and this band has a negligible overlap with the bands containing the other modes.
The statistical analysis showed mild variations in the properties of the first modes of vibrations of the robot, due to variations in J2 and J3. These joint variables determine the configuration of the robot in the meridian plane and in particular the location of the wrist (the origin of coordinate system 4 in Figure 2). Therefore, a further numerical analysis was carried out to study the dependence of the first three natural frequencies on the position of the wrist in the meridian plane. The results are presented in terms of contour plots, in which the darker colors represent the lower frequencies and the lighter colors represent the higher frequencies. Since some locations can be reached with two robot configurations (elbow-up and elbow-down), both configurations are considered in separate plots. Figure 14 shows that the first natural frequency decreases in a regular way when the distance between the wrist and the robot base increases, and this effect takes place both in the elbow-up and in the elbow-down configuration.
The contour plot of the second natural frequency is more complex (Figure 15). The second natural frequency reaches the lowest values when the arm is extended in order to reach locations in front of the base or on the back of the base. In these configurations, links 2 and 3 tend to align. The maximum value of the second natural frequency takes place when the arm has to reach locations above the base. This effect is present in both the considered configurations. Figure 16 shows a rather regular trend of the third natural frequency. In both configurations, this natural frequency tends to increase when the wrist location is far from the base and links 2 and 3 tend to align (joint variable J3 is close to 90°). Therefore, the angle ( ) between link 2 and 3 has opposite effects on the second and third natural frequencies that usually correspond to modes that take place in the meridian plane. The second natural frequency decreases when the two links tend to align, because the direct inertia term The natural frequencies of the second and third mode, which are usually characterized by large displacements in the meridian plane, show the largest standard deviations. Therefore, these modes are influenced by the robot configuration. The first natural frequency, which usually corresponds to a mode of vibration with large displacements in the horizontal plane, shows a mild standard deviation, since the robot configuration affects the moment of inertia about J1. The natural frequencies of the last two modes, which are chiefly dominated by J4 and J5, show very small standard deviations. This happens because the robot configuration has a small effect on the direct inertia term of J4 and no effect on the direct inertia term of J5, but it only changes the inertial cross-coupling terms (see Figure 3).
The statistical analysis of modal damping variation with robot configurations is depicted in Figure 13. This figure shows that robot configuration affects the modal damping of the first three modes, since standard deviations are comparable with the mean values. Conversely, robot configuration has a very small effect on the damping of the last modes, which chiefly involve wrist joints, because standard deviations are very small. Figures 12b and 13b show the result of a similar calculation, in which only random variations in J2 and J3 were performed. The standard deviations of the natural frequencies and modal dampings are very similar to the ones of Figures 12a and 13a, therefore, this result points out that the variations in J2 and J3 have the largest effect on the modal properties of the robot. Figures 12 and 13 highlight another important feature of the robot. Variations in joint variables are able to modify the natural frequencies, but they do not change the typical frequency bands of the modes. In other words, the frequencies of the modes associated to the wrist joints are always higher than the frequencies associated to the arm joints. The first two modes belong to the same frequency band (15 ÷ 30 Hz), and this band has a negligible overlap with the bands containing the other modes.
The statistical analysis showed mild variations in the properties of the first modes of vibrations of the robot, due to variations in J2 and J3. These joint variables determine the configuration of the robot in the meridian plane and in particular the location of the wrist (the origin of coordinate system 4 in Figure 2). Therefore, a further numerical analysis was carried out to study the dependence of the first three natural frequencies on the position of the wrist in the meridian plane. The results are presented in terms of contour plots, in which the darker colors represent the lower frequencies and the lighter colors represent the higher frequencies. Since some locations can be reached with two robot configurations (elbow-up and elbow-down), both configurations are considered in separate plots. Figure 14 shows that the first natural frequency decreases in a regular way when the distance between the wrist and the robot base increases, and this effect takes place both in the elbow-up and in the elbow-down configuration.
increases. The third natural frequency increases when the two links tend to align, because the inertial cross coupling term ( ) between J2 and J3 increases (see Figure 5). This effect takes place even if a simple 2 DOF model of the robot in the meridian plane is adopted [28].
Finally, the analysis of Figure 15 and 16 shows that it is rather difficult to maximize both the second and the third natural frequencies. Only in a limited region of the workspace above the robot it is possible to achieve large values of both frequencies.    The contour plot of the second natural frequency is more complex (Figure 15). The second natural frequency reaches the lowest values when the arm is extended in order to reach locations in front of the base or on the back of the base. In these configurations, links 2 and 3 tend to align. The maximum value of the second natural frequency takes place when the arm has to reach locations above the base. This effect is present in both the considered configurations.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 17 of 19 increases. The third natural frequency increases when the two links tend to align, because the inertial cross coupling term ( ) between J2 and J3 increases (see Figure 5). This effect takes place even if a simple 2 DOF model of the robot in the meridian plane is adopted [28].
Finally, the analysis of Figure 15 and 16 shows that it is rather difficult to maximize both the second and the third natural frequencies. Only in a limited region of the workspace above the robot it is possible to achieve large values of both frequencies.     Figure 16 shows a rather regular trend of the third natural frequency. In both configurations, this natural frequency tends to increase when the wrist location is far from the base and links 2 and 3 tend to align (joint variable J3 is close to 90 • ).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 17 of 19 increases. The third natural frequency increases when the two links tend to align, because the inertial cross coupling term ( ) between J2 and J3 increases (see Figure 5). This effect takes place even if a simple 2 DOF model of the robot in the meridian plane is adopted [28].
Finally, the analysis of Figure 15 and 16 shows that it is rather difficult to maximize both the second and the third natural frequencies. Only in a limited region of the workspace above the robot it is possible to achieve large values of both frequencies.    Therefore, the angle (q 3 ) between link 2 and 3 has opposite effects on the second and third natural frequencies that usually correspond to modes that take place in the meridian plane. The second natural frequency decreases when the two links tend to align, because the direct inertia term M 22 increases. The third natural frequency increases when the two links tend to align, because the inertial cross coupling term (M 23 ) between J2 and J3 increases (see Figure 5). This effect takes place even if a simple 2 DOF model of the robot in the meridian plane is adopted [28].
Finally, the analysis of Figures 15 and 16 shows that it is rather difficult to maximize both the second and the third natural frequencies. Only in a limited region of the workspace above the robot it is possible to achieve large values of both frequencies.

Conclusions
The proposed testing method based on selective modal testing requires simple equipment, does not take a long testing time, can be implemented in every robot configuration, and makes it possible to detect the presence of modes of vibration affected by structural and bearing compliance.
The combination of the results of selective modal testing with a mathematical model of mass distribution of the robot makes it possible a straightforward identification of joint stiffness and damping. A simple optimization procedure on a small set of parameters is carried out, to identify joint stiffness in the presence of a strong inertial cross-coupling between the joints.
The mathematical model implemented with the identified stiffness and damping properties gives information about the variation of the dynamic characteristics of the robot in the workspace. This information is useful to optimize the performance of the robot in machining tasks.
In future, the method will be improved, in order to cope with the presence of the structural mode in the band of frequency of interest.