New Indices for Evaluating Vibration Characteristics of Flexible-Joint Robots

: Structural vibration is a signiﬁcant consideration for robotic applications such as machining where the robot is subject to large dynamic loading. Aiming at providing an e ﬃ cient means to evaluating the vibration characteristics of industrial robots for these applications, this work proposes two new indices to quantify the elastic displacement of the tool mounted on the robot caused by the vibrations induced by external process loading for ﬂexible-joint robots. For this purpose, a structural dynamic model is ﬁrst developed to derive the frequency responses of the tool displacement. Then, the displacement-force and displacement-torque frequency response ratios are deﬁned, which represent the mapping from the amplitudes of an external harmonic force and torque to the amplitude of tool displacement respectively. The upper bounds of the two ratios are used as evaluation indices for the vibration characteristics of the robot, which represent the worst situation of the tool displacement due to harmonic excitation with amplitude of unit force and unit torque respectively. With these indices, an e ﬃ cient method is provided to predict whether the tool misalignment caused by periodic loading is acceptable for process quality requirement. Numerical simulation demonstrates the e ﬀ ectiveness of the proposed method for a robotic riveting system being developed for aerospace assembly. natural frequencies of the robot. These peak values show that only considering the elastostatic performance is insu ﬃ cient for dynamic processes. An eigenvector analysis can provide the amplitudes and phases of the force and torque components under which the worst situation of the tool vibration occurs. For a given process, the maximum tool misalignment can be easily evaluated by combining the robot-dependent indices and the process-dependent Fourier series, so that the robot users can consider


Introduction
Although industrial robots have been successfully used in production automation for a wide range of tasks including painting, welding, packaging, and material handling [1], they still face challenges for more complex and dynamic applications like machining that involve large time-varying loading during the processes and require high positioning accuracy. Recent development of these robotic applications include the European projects HEPHESTOS and COMET that aimed at developing robotic systems for small-batch machining [2,3] and Boeing's initiative to use robots for automatic drilling and riveting in aircraft assembly [4]. The challenges for these robotic processes result from inherent disadvantages of industrial robots compared with machine tools, such as low positioning accuracy and being prone to vibrations [5]. Modern metrology techniques such as laser tracking [6,7] and photogrammetry based measurement [8,9] can be adopted to improve the accuracy of industrial robots through kinematic calibration [10,11] and positioning error compensation [12][13][14]. However, the influence of robot vibration on the processes is difficult to mitigate by these metrology techniques.
These robotic applications involve dynamic loading in a way that the process generates force and torque acting on the robot whose magnitudes and directions dramatically change with time. The process force and torque can excite vibrations into robot structure and lead to much more complicated dynamics than conventional robotic tasks like welding and material handling where only static or low contact forces act on the robot. Since the natural frequencies of an industrial robot strongly depend on its configuration [15] and its stiffness varies significantly in different directions [16,17], the vibration characteristics depend on the robot configuration and the direction of external force. It has been found that structural vibrations of the robot can cause tool misalignment [18], damage the surface quality [19], and reduce the fatigue lives of the tool and the robot [20,21]. Furthermore, the system can become instable due to the interaction between the robot vibrations and the process forces [22,23] and the stability condition can be affected by some process parameters such as feed rate and spindle speed [24,25]. All these have indicated the complexity of vibration characteristics of robots under dynamic loading and the possible serious effect on process quality. A first question for developing an automation system for such processes is whether to use and how to select a robot so that the forced vibrations excited by the process are acceptable for quality requirement. This brings the need for an efficient means to evaluating the vibration characteristics of robots, which, however, is usually not provided by robot manufacturers.
A variety of indices can be found in literature to evaluate robot performance [26]. For instance, the condition number of Jacobian matrix has been proposed to measure the robot dexterity [27], and the manipulating force ellipsoid can quantify the static torque-force transmission from robot joints to its end-effector [28]. The dynamic manipulability ellipsoid and the radius of acceleration have been derived for evaluating robots' acceleration performance respectively. The former represents the achievable end-effector accelerations in a configuration for given joint torque limits [29], and the latter is a global performance measure of the achievable acceleration from any joint position and velocity in the operating region [30]. The inertial matching ellipsoid is a performance index combining dynamic manipulability and manipulating force ellipsoids which represents the dynamic torque-force transmission efficiency between the joint torque and the static force applied to a part mounted on the robot end-effector [31]. These performance indices are derived based on the assumption that the robot is rigid and are usually proposed for dimensional synthesis. Stiffness is an important structural consideration for robotic applications and represents the effect of applied static forces and torques on the end-effector deformation. Similar to the dexterity index, Liu et al. [32] used the condition number of the stiffness matrix to evaluate the stiffness of 3-DOF (degrees of freedom) spherical parallel robots. Courteille et al. [33] chose the maximum and minimum eigenvalues of stiffness matrix to calculate the upper and lower bounds of the robot stiffness which occurs in the corresponding eigenvector directions. However, none of the aforementioned indices can be used to evaluate the structural vibration characteristics when the robot end-effector affords dynamic loading. Aiming at addressing this issue, the main contribution of this work is to derive two indices which can be used for efficient evaluation of the vibration characteristics of flexible-joint robots under a dynamic force and torque, respectively.
For typical industrial articulated robots, the joints have been considered as most significant flexible components, due to the deformations of transmission components such as gears, shafts, belts and harmonic drives [1,26]. Therefore, we limit our scope to the robots with rigid links and flexible joints in this work. It is also assumed that the influence of the robot deformation on the external force and torque can be neglected. Two indices are defined for evaluating the vibration characteristics of robots under a dynamic force and torque, respectively. To this end, a structural dynamic model is first established for flexible-joint robots that undergo forced vibrations, and obtain the frequency response of the tool displacement due to the process excitation. Then, the displacement-force and displacement-torque frequency response ratios are derived from the frequency-domain analysis, which represent the mapping from the amplitudes of external force and torque to the amplitude of tool displacement due to forced vibrations, respectively. The upper bounds of the two ratios are used as indices to describe the worst situation of robot vibrations under a harmonic excitation with an amplitude of unit force and unit torque, respectively. The development leads to an efficient method to evaluate the tool misalignment caused by periodic process excitation.

System Description
A prototype mobile robotic assembly riveting system at Shanghai University for the milling, drilling and percussive riveting of sheet metals for aerospace assembly. As displayed in Figure 1, an M-20iA industrial robot from Fanuc China in Shanghai with six revolute joints is mounted on a mobile platform driven by four Mecanum wheels. The platform can move in an arbitrary direction on the ground, and the base of the robot can also move vertically for tasks at different heights. Time-varying machining forces and torques act on the machining robot in different directions during the processes. In addition, the excitation frequencies can be as low as 5 Hz for percussive riveting and as high as 50 Hz for drilling and milling. Thus, the problem under study can be stated as to evaluate whether the elastic displacement of tool due to the forced vibrations during these processes is within certain tolerance to guarantee the process quality.

System Description
A prototype mobile robotic assembly riveting system at Shanghai University for the milling, drilling and percussive riveting of sheet metals for aerospace assembly. As displayed in Figure 1, an M-20iA industrial robot from Fanuc China in Shanghai with six revolute joints is mounted on a mobile platform driven by four Mecanum wheels. The platform can move in an arbitrary direction on the ground, and the base of the robot can also move vertically for tasks at different heights. Timevarying machining forces and torques act on the machining robot in different directions during the processes. In addition, the excitation frequencies can be as low as 5 Hz for percussive riveting and as high as 50 Hz for drilling and milling. Thus, the problem under study can be stated as to evaluate whether the elastic displacement of tool due to the forced vibrations during these processes is within certain tolerance to guarantee the process quality.

Equations of Motion
The kinematic scheme of a robot with n links and n revolute joints is plotted in Figure 2. The position vector of the ith joint in the inertial frame {OXYZ} is denoted by Pi. A local body frame {oixiyizi} is established on each of the links with the origin at the ith joint. A tool is mounted on the robot's end-effector to complete the process and the tool tip is denoted as E, with its position vector written in the inertial frame as E. Then the external force and torque acting on the tool from the workpiece are denoted as F and τ respectively. Additionally, to describe the process force and torque, a work frame {o'x'y'z'} is established and fixed on the part. In this work, it is assumed that the robot links are rigid and that the joints can be modeled as linearly elastic torsional springs with small deformations in series with the joint actuators [34]. To

Equations of Motion
The kinematic scheme of a robot with n links and n revolute joints is plotted in Figure 2. The position vector of the ith joint in the inertial frame {OXYZ} is denoted by P i . A local body frame {o i x i y i z i } is established on each of the links with the origin at the ith joint. A tool is mounted on the robot's end-effector to complete the process and the tool tip is denoted as E, with its position vector written in the inertial frame as E. Then the external force and torque acting on the tool from the workpiece are denoted as F and τ respectively. Additionally, to describe the process force and torque, a work frame {o'x'y'z'} is established and fixed on the part.

System Description
A prototype mobile robotic assembly riveting system at Shanghai University for the milling, drilling and percussive riveting of sheet metals for aerospace assembly. As displayed in Figure 1, an M-20iA industrial robot from Fanuc China in Shanghai with six revolute joints is mounted on a mobile platform driven by four Mecanum wheels. The platform can move in an arbitrary direction on the ground, and the base of the robot can also move vertically for tasks at different heights. Timevarying machining forces and torques act on the machining robot in different directions during the processes. In addition, the excitation frequencies can be as low as 5 Hz for percussive riveting and as high as 50 Hz for drilling and milling. Thus, the problem under study can be stated as to evaluate whether the elastic displacement of tool due to the forced vibrations during these processes is within certain tolerance to guarantee the process quality.

Equations of Motion
The kinematic scheme of a robot with n links and n revolute joints is plotted in Figure 2. The position vector of the ith joint in the inertial frame {OXYZ} is denoted by Pi. A local body frame {oixiyizi} is established on each of the links with the origin at the ith joint. A tool is mounted on the robot's end-effector to complete the process and the tool tip is denoted as E, with its position vector written in the inertial frame as E. Then the external force and torque acting on the tool from the workpiece are denoted as F and τ respectively. Additionally, to describe the process force and torque, a work frame {o'x'y'z'} is established and fixed on the part. In this work, it is assumed that the robot links are rigid and that the joints can be modeled as linearly elastic torsional springs with small deformations in series with the joint actuators [34]. To In this work, it is assumed that the robot links are rigid and that the joints can be modeled as linearly elastic torsional springs with small deformations in series with the joint actuators [34]. To evaluate the elastic displacement of the tool, we develop a structural model of the robot to analyze the forced vibrations of the joints in a given configuration denoted by the input joint angles as q = [θ 1 , θ 2 , . . . , θ n ] T as follows [28]: where ∆q is the elastic joint deflection from their equilibrium values q e and represents the torsional deflections of the joints, M(q) is the n × n symmetric generalized inertial matrix and is configuration dependent, K is the diagonal stiffness matrix with the diagonal entries equal to the torsional stiffness of the corresponding joints, and C is a damping matrix. For a complete formulation, matrix C should include the effect of the Coriolis and centrifugal terms related to the rigid-body motion on the elastic deformations. But this work is mainly focused on the robotic percussive riveting operations where the robot moves to a required spot and then stops there during the riveting. In other words, the robot has no rigid-body motion when its end-effector is subject to dynamic loading. Thus, it can be assumed that the coupling between the vibration and the rigid-body motion can be neglected, and then matrix C here only includes the structural damping. The 3 × n matrices J v and J w can be obtained from the Jacobian matrix J that describes the mapping from joint rates to the twist of the end-effector as below: q is the joint rates, v is the velocity of the tool tip, and ω is the angular velocity of the tool. We assume that the joint deflections are small, then the linear displacement of the tool tip can be obtained as: From the perspective of process analysis, Equation (3) gives the tool misalignment due to the structural deformation caused by the process force and torque acting on the robot.

Frequency Response Functions
We first analyze the frequency response of the displacement of the tool tip due to an external force. To do this, we assume that the end-effector is subject to a simple harmonic force with a frequency ω which can be written in vector form as: where F is a 3 × 1 complex vector whose norm and phase angle of each component represent the amplitude and phase of the force in the corresponding direction. Then the steady-state response of the joint deflection ∆q can be expressed as: where ∆q v is an n ×1 complex vector which contains the amplitude and phase of the steady-state response of the joint deflections. Substituting (4) and (5) into (1) leads to: The steady-state response of the tool displacement due to the harmonic force F is then computed from (3) as below: where H v (ω) is the 3 × 3 displacement-force frequency response function as follows: Similarly, when the end-effector is subject to a simple harmonic torque with a frequency ω as: where τ is a 3 × 1 complex vector which contains the amplitudes and phases of the three torque components, then the steady-state response of joint deflection can be written as: Now, the steady-state response of the tool displacement due to the harmonic torque τ can be computed as: where H w (ω) is the 3 × 3 displacement-torque frequency response function as follows: The frequency response functions H v (ω) and H w (ω) represent the elastodynamic compliance of the robot under simple harmonic force and torque, respectively. They not only depend on the robot configuration, but also depend on the excitation frequency. When ω is equal to zero, i.e., the robot is subject to a static force or torque, then (8) and (12) are reduced to: which represent the elastostatic compliance of the robot under static force and torque respectively.

Frequency Response Ratios
Now we investigate the forced vibration characteristics of the robot under a harmonic force. From (7), we have: where the superscript H represents the conjugate transpose of a complex matrix. Then we define the displacement-force frequency response ratio as: Hermitian matrix. The ratio η v represents the mapping from the amplitude of external force to the amplitude of tool displacement due to structural vibration, and η v depends on the robot configuration and the force frequency. Furthermore, the ratio is related to the amplitudes and phases of the three force components in the complex vector F. When the robot configuration and force frequency are given, the upper and lower bounds of η v with different F can be computed as [35]: where λ vmin and λ vmax are the minimum and maximum eigenvalues of the Hermitian matrix A v , respectively. The force vector F corresponding to the upper and lower bounds of η v are: where e vmin and e vmax are the eigenvectors of A v corresponding to λ vmin and λ vmax respectively. Equations (16) and (17) give the range of the displacement-force frequency response ratio η v and the amplitudes and phases of the force components for the boundary values of η v . Similar analysis can be carried out on the forced vibration characteristics of the robot under a simple harmonic torque. We define the displacement-torque frequency response ratio as: is also a 3 × 3 Hermitian matrix. The ratio η w represents the mapping from the amplitude of external torque to the amplitude of tool displacement due to structural vibration. When the robot configuration and torque frequency are given, the upper and lower bounds of η w with different τ can be obtained as [34]: where λ wmin and λ wmax are the minimum and maximum eigenvalues of the Hermitian matrix A w respectively, and the corresponding torque vector are given as: where e wmin and e wmax are the eigenvectors of A w corresponding to λ wmin and λ wmax respectively. Equations (19) and (20) give the range of the displacement-torque frequency response ratio η w and the amplitudes and phases of the torque components for the boundary values of η w . The maximum frequency response ratios η vmax and η wmax provide a measure of the maximum amplitudes of tool displacement due to harmonic loading with an amplitude of unit force and unit torque respectively, and the corresponding eigenvectors represent the amplitudes and phases of the force and torque components under which the worst vibration situation occurs. The ratios η vmax and η wmax are inherent characteristics of the robot and independent of the application processes, therefore they can be used as indices to evaluate the vibration characteristics of the robot under harmonic excitation.
More importantly, the eigenvectors e vmax and e wmax represents the amplitudes and phases of the force and torque components corresponding to η vmax and η wmax . The two indices can be used to predict the worst situation of forced vibrations of the robot for given processes. If the rotation matrix from the work frame (see Figure 2) to the inertial frame is written as R, then the force and torque acting on the tool can be written as: where F' and τ' are the process-dependent force and torque written in the work frame. In the case of machining, when the effect of robot deformation on F' and τ' can be neglected, these force and torque can be predicted via machining force modeling or via machining tests on Computer Numerical Control (CNC) machines. Usually, the process force and torque are periodic functions about time and can be approximated as an expansion of 2N + 1 Fourier series as follows: where ω v and ω w are the force and torque frequencies, F k and τ k are complex vectors containing the Fourier coefficients and: where F * k and τ * k are the conjugates of F k and τ k respectively. Then, the maximum amplitude of tool displacement due to the kth Fourier series in F and τ are respectively computed as: In (24), the Fourier coefficients F k and τ k are process dependent and the indices η vmax and η wmax are robot dependent. Then, once η vmax and η wmax had been provided, for example, by robot manufacturers and F k and τ k are modeled or measured by robot users, the largest tool misalignment due to each Fourier series can be easily evaluated using (24).

System Parameters
In this case study, the Fanuc M-20iA robot for machining and assembly shown in Figure 1 is simulated. The reference frames of the robot are drawn in Figure 3. Only the deflections of the first three joints are considered since the moment arms for the external force acting at the Tool Center Point (TCP) about the last three joints are much shorter than those about the first three joints. Note that the modeling approach presented in Section III is generic and can be adopted to robots with over three flexible revolute joints. The inertial frame {oxyz} is established at the robot base with its positive Z axis vertically upward and its Y axis perpendicular to the plane containing the second and third links. The z i axes of the local frames are in the rotation directions of the ith joints; when the robot is stationary, the x 1 axis is parallel to the X axis, and the x 2 and x 3 axes are along the link length directions. The directions of other reference axes can be determined by the right-hand rule.

System Parameters
In this case study, the Fanuc M-20iA robot for machining and assembly shown in Figure 1 is simulated. The reference frames of the robot are drawn in Figure 3. Only the deflections of the first three joints are considered since the moment arms for the external force acting at the Tool Center Point (TCP) about the last three joints are much shorter than those about the first three joints. Note that the modeling approach presented in Section III is generic and can be adopted to robots with over three flexible revolute joints. The inertial frame {oxyz} is established at the robot base with its positive Z axis vertically upward and its Y axis perpendicular to the plane containing the second and third links. The z i axes of the local frames are in the rotation directions of the ith joints; when the robot is stationary, the x 1 axis is parallel to the X axis, and the 2 and x 3 axes are along the link length directions. The directions of other reference axes can be determined by the right-hand rule.  Table 1. The torsional stiffness of the joints is given as: = 2.4 × 10 5 N-m/rad, k 2 = k 3 =10 5 N-m/rad. A damping ratio of 0.06 is used to compute the damping matrix in this simulation based on the modal test data for industrial robots given in [36]. The ranges of the joint angles are θ 1 ∈[-π, π], θ 2 ∈[-π/2, π/2], and θ 3 ∈[-π/3, π/3]. For this assembly application, the robot moves within a cuboid shaped operating workspace that can be represented by the width a = 1.8 m, the height h = 0.8 m, and the depth b = 0.5 m shown in Figure 3. The coordinates of the center of the operating workspace in the inertial frame are x = 0.75 m, y = 0 m, and z = 0.65 m. The Denavit-Hartenberg parameters of the robot are listed in Table 2. Lastly, the mass of the rivet  Table 1. The torsional stiffness of the joints is given as: k = 2.4 × 10 5 N-m/rad, k 2 = k 3 = 10 5 N-m/rad. A damping ratio of 0.06 is used to compute the damping matrix in this simulation based on the modal test data for industrial robots given in [36]. The ranges of the joint angles are θ 1 ∈ [−π, π], θ 2 ∈ [−π/2, π/2], and θ 3 ∈ [−π/3, π/3]. For this assembly application, the robot moves within a cuboid shaped operating workspace that can be represented by the width a = 1.8 m, the height h = 0.8 m, and the depth b = 0.5 m shown in Figure 3. The coordinates of the center of the operating workspace in the inertial frame are x = 0.75 m, y = 0 m, and z = 0.65 m. The Denavit-Hartenberg parameters of the robot are listed in Table 2. Lastly, the mass of the rivet gun is 5.8 kg. More information about this robotic application can be in [37].

Performance Indices for Given Positions
The proposed evaluation model has been implemented in the MATLAB environment. We first investigate the variation of performance indices η vmax and η wmax with different frequencies when the robot moves to given positions. As an example, the values of η vmax and η wmax with frequencies from 0 to 200 Hz at Point A (x = 0.1 m, y = 0, and z = 1 m) and Point B (x = 1 m, y = 0, and z = 1 m) are compared in Figures 4 and 5 respectively. It is seen that when the amplitude of external force/torque is equal to unit force/torque, the maximum amplitude of TCP displacement strongly depends on the excitation frequency and can be much higher than the largest static displacement when a unit static force/torque acts on the end-effector (represented by the values of η vmax or η wmax when the frequency is equal to zero). This implies that only considering the elastostatic performance is insufficient for dynamic robotic processes.

Performance Indices for Given Positions
The proposed evaluation model has been implemented in the MATLAB environment. We first investigate the variation of performance indices ηvmax and ηwmax with different frequencies when the robot moves to given positions. As an example, the values of ηvmax and ηwmax with frequencies from 0 to 200 Hz at Point A (x = 0.1 m, y = 0, and z = 1 m) and Point B (x = 1 m, y = 0, and z = 1 m) are compared in Figures 4 and 5 respectively. It is seen that when the amplitude of external force/torque is equal to unit force/torque, the maximum amplitude of TCP displacement strongly depends on the excitation frequency and can be much higher than the largest static displacement when a unit static force/torque acts on the end-effector (represented by the values of ηvmax or ηwmax when the frequency is equal to zero). This implies that only considering the elastostatic performance is insufficient for dynamic robotic processes.     At Point A, it is found that the first two modes of the robot have natural frequencies of 54 rad/s and 103 rad/s respectively, and both mode shapes are related to the deflections of Joint 2 and Joint 3. The third mode of the robot has a natural frequency of 155 rad/s, and its mode shape is related to the deflection of Joint 1. From Figures 4 and 5, we observe two peak values of η υmax and η wmax at frequencies close to the first two natural frequencies. As the moment arm about Joint 1 is small at Point A, we do not observe a third peak value related to the third mode on both curves of η υmax and η wmax , which indicates that the TCP vibration in the y direction is less obvious than those in the oxz plane. We further investigate the eigenvectors of matrices A v and A w corresponding to η υmax and η wmax at different frequencies.  N-m). Both indicate that the worst case of the TCP vibration due to a harmonic force or torque with this frequency is caused by the deflection of Joint 1.
The above eigenvalue and eigenvector analysis procedure can be used to analyze the worst situation of forced vibrations of the robot in other positions and at other frequencies.

Distributions of Performance Indices
Now we analyze the distributions of the two performance indices η υmax and η wmax over the robot workspace. As shown in Figures 6 and 7, these distributions are plotted in the oxz plane due to the symmetry about Joint 1, with force and torque frequencies of 20, 40, 60, and 80 rad/s. It is found that these dynamics performance indices dramatically change with the excitation frequencies.
Generally speaking, the structural dynamics performance at 40 rad/s is lower than those at other three frequencies. Moreover, at each frequency, the maximum frequency response ratios η υmax and η wmax in certain positions are significantly higher than those in other positions, and a further observation finds that the natural frequencies in these positions are close to the excitation frequency. A detailed investigation on the eigenvectors can give the amplitude and phase of each force/torque component under which η υmax and η wmax occurs.       The mean value and maximum value η υmax and η wmax of over the robot workspace with different excitation frequencies are also listed in Table 3. Again, it is seen that both indices strongly depend on the frequencies of the external force or moment, and that the worst structural dynamics performance occurs with an excitation frequency of about 40 rad/s.

Performance Evaluation for Percussive Riveting
For the application scenario as illustrated in Figure 1, the robot can be subject to forces and torques in different directions and due to different processes within its operating workspace. As an example, we demonstrate the evaluation of vibration characteristics for robotic percussive riveting. This operation is completed by the robot with a percussive rivet gun and the principle of percussive riveting is illustrated in Figure 8. A rivet is first inserted into pre-drilled mating holes on the panels. Then the simulated robot drives a percussive rivet gun and another supporting robot (see Figure 8) drives a bucking bar to the riveting spot. Once the rivet gun and the bucking bar are aligned to the normal direction of the panels, the robot stops and activates the gun, and then the hammer of the gun repeatedly hits the rivet head. The impact forces cause the rivet to deform. When the diameter of the rivet becomes equal to the diameter of the pre-drilled holes, the panels are joined. For motion planning of the robotic riveting operations, operational space trajectory planning is used to avoid interference during the motion. Note that since the robot is stationary when its end-effector is subject to dynamic loading during its riveting operations, the trajectory between every two riveting spots does not affect the vibration characteristics if the forced vibration can be damped out during the motion from one spot to the next.

Performance Evaluation for Percussive Riveting
For the application scenario as illustrated in Figure 1, the robot can be subject to forces and torques in different directions and due to different processes within its operating workspace. As an example, we demonstrate the evaluation of vibration characteristics for robotic percussive riveting. This operation is completed by the robot with a percussive rivet gun and the principle of percussive riveting is illustrated in Figure 8. A rivet is first inserted into pre-drilled mating holes on the panels. Then the simulated robot drives a percussive rivet gun and another supporting robot (see Figure 8) drives a bucking bar to the riveting spot. Once the rivet gun and the bucking bar are aligned to the normal direction of the panels, the robot stops and activates the gun, and then the hammer of the gun repeatedly hits the rivet head. The impact forces cause the rivet to deform. When the diameter of the rivet becomes equal to the diameter of the pre-drilled holes, the panels are joined. For motion planning of the robotic riveting operations, operational space trajectory planning is used to avoid interference during the motion. Note that since the robot is stationary when its end-effector is subject to dynamic loading during its riveting operations, the trajectory between every two riveting spots does not affect the vibration characteristics if the forced vibration can be damped out during the motion from one spot to the next.

Rivet gun
Bucking bar Panels Hammer Bucking bar Panels Repetitive impacts During percussive riveting, the end-effector is subject to repetitive impact forces normal to the panel. If the work frame is established on the panel at the riveting spot and with the ' z axis normal to the panel, then the process force and torque can be written in the work frame as below: The percussive impact force F(t) is plotted in Figure 9a, where the amplitude Fm is 550 N, the force frequency ων is about 62.8 rad/s (10 Hz), and Td/T is 0.08. The force F(t) is approximated by Fourier expansion using (22), with the Fourier coefficients drawn in Figure 9b Note that the force in (26) can be in different directions relative to the inertial frame for different parts with transformation (21). During percussive riveting, the end-effector is subject to repetitive impact forces normal to the panel. If the work frame is established on the panel at the riveting spot and with the z axis normal to the panel, then the process force and torque can be written in the work frame as below: The percussive impact force F(t) is plotted in Figure 9a, where the amplitude F m is 550 N, the force frequency ω ν is about 62.8 rad/s (10 Hz), and T d /T is 0.08. The force F(t) is approximated by Fourier expansion using (22), with the Fourier coefficients drawn in Figure 9b for k = 0, ±1, ±2, . . . , ±20, i.e., Note that the force in (26) can be in different directions relative to the inertial frame for different parts with transformation (21). The robot is used for the assembly in a cuboid-shaped operation workspace as drawn in Figure  10. The operation workspace can be determined by Points Pk with k from 1 to 8. The coordinates of P1 is [1, 0.5, 1] T (in m) and the length, height, and width of the operation workspace are P1 P2 = 1 m, P1 P4 = 0.8 m and P1 P5 = 0.5 m, respectively. The robot configurations corresponding to Points P1 to P4 are also plotted in Figure 10. To demonstrate the vibration of the robot joints due to the percussive impact force, the time and frequency responses of the tool displacement are analyzed for Point P1 when the percussive impact is along the -X direction (see Figure 10). The frequency response function for the tool displacement due to a force along the-X direction can be obtained from (7), and the magnitude of the frequency response function with different frequencies are plotted in Figure 11. If the dynamic loading is approximated as an expansion of Fourier series as (26). Then, the amplitude of the tool vibration depends on the Fourier coefficients and the frequency response function at the corresponding frequencies (  v ,  2 v ,  3 v as shown in Figure 11). It can be seen that the first two Fourier series have larger influence on the tool vibration than the other Fourier series. In addition, although the riveting force is along the -X direction, the tool displacement in the Z direction has a larger amplitude than those in the other two directions. Furthermore, the time response of the tool displacement can be computed by solving (1) and (3). The resulting tool displacement along different directions are plotted in Figure 12 for five periods of the external force (about 0.5 s). Again, the tool displacement in the Z direction has the largest amplitude among the three directions.  The robot is used for the assembly in a cuboid-shaped operation workspace as drawn in Figure 10. The operation workspace can be determined by Points P k with k from 1 to 8. The coordinates of P 1 is [1, 0.5, 1] T (in m) and the length, height, and width of the operation workspace are P 1 P 2 = 1 m, P 1 P 4 = 0.8 m and P 1 P 5 = 0.5 m, respectively. The robot configurations corresponding to Points P 1 to P 4 are also plotted in Figure 10. To demonstrate the vibration of the robot joints due to the percussive impact force, the time and frequency responses of the tool displacement are analyzed for Point P 1 when the percussive impact is along the -X direction (see Figure 10). The frequency response function for the tool displacement due to a force along the-X direction can be obtained from (7), and the magnitude of the frequency response function with different frequencies are plotted in Figure 11. If the dynamic loading is approximated as an expansion of Fourier series as (26). Then, the amplitude of the tool vibration depends on the Fourier coefficients and the frequency response function at the corresponding frequencies (ω v , 2ω v , 3ω v as shown in Figure 11). It can be seen that the first two Fourier series have larger influence on the tool vibration than the other Fourier series. In addition, although the riveting force is along the -X direction, the tool displacement in the Z direction has a larger amplitude than those in the other two directions. Furthermore, the time response of the tool displacement can be computed by solving (1) and (3). The resulting tool displacement along different directions are plotted in Figure 12 for five periods of the external force (about 0.5 s). Again, the tool displacement in the Z direction has the largest amplitude among the three directions. The robot is used for the assembly in a cuboid-shaped operation workspace as drawn in Figure  10. The operation workspace can be determined by Points Pk with k from 1 to 8. The coordinates of P1 is [1, 0.5, 1] T (in m) and the length, height, and width of the operation workspace are P1 P2 = 1 m, P1 P4 = 0.8 m and P1 P5 = 0.5 m, respectively. The robot configurations corresponding to Points P1 to P4 are also plotted in Figure 10. To demonstrate the vibration of the robot joints due to the percussive impact force, the time and frequency responses of the tool displacement are analyzed for Point P1 when the percussive impact is along the -X direction (see Figure 10). The frequency response function for the tool displacement due to a force along the-X direction can be obtained from (7), and the magnitude of the frequency response function with different frequencies are plotted in Figure 11. If the dynamic loading is approximated as an expansion of Fourier series as (26). Then, the amplitude of the tool vibration depends on the Fourier coefficients and the frequency response function at the corresponding frequencies (  v ,  2 v ,  3 v as shown in Figure 11). It can be seen that the first two Fourier series have larger influence on the tool vibration than the other Fourier series. In addition, although the riveting force is along the -X direction, the tool displacement in the Z direction has a larger amplitude than those in the other two directions. Furthermore, the time response of the tool displacement can be computed by solving (1) and (3). The resulting tool displacement along different directions are plotted in Figure 12 for five periods of the external force (about 0.5 s). Again, the tool displacement in the Z direction has the largest amplitude among the three directions.    Now, we proceed to analyze the potential risk of the tool vibration when the percussive force acts in different directions on the robot end-effector. Once the maximum displacement-force frequency response ratio η vmax is given, then the maximum amplitude of TCP displacement due to the kth Fourier series of the periodic impact force ΔE v max,k can be easily calculated using (24). For example, the values of ΔE v max,k are plotted over the potential workspace when x is 0.5, 0.6, 0.7, 0.8, 0.9 and 1 m in Figure 13a,b, for k = ± 1 and ± 2 respectively. It can be found that the distributions of ΔE v max,k are symmetric about the oxz plane. When only considering the Fourier series with k = ± 1, the amplitude of the TCP displacement can be as high as 2 × 10 −3 m in certain positions where the first natural frequency of the robot (related to deflections of Joint 2 and Joint 3) is close to ω v . When only considering the Fourier series with k = ± 2, the amplitude of the TCP displacement can be as high as 2.3 × 10 −3 m in the positions where the second natural frequency (related to the deflection of Joint 1) is close to 2ω v . It is also observed that when k > 4, kω v is beyond the robot's natural frequencies, and the values of ΔE v max,k are much smaller than those for k = ± 1, ± 2, ± 3 and ± 4.  Now, we proceed to analyze the potential risk of the tool vibration when the percussive force acts in different directions on the robot end-effector. Once the maximum displacement-force frequency response ratio η vmax is given, then the maximum amplitude of TCP displacement due to the kth Fourier series of the periodic impact force ΔE v max,k can be easily calculated using (24). For example, the values of ΔE v max,k are plotted over the potential workspace when x is 0.5, 0.6, 0.7, 0.8, 0.9 and 1 m in Figure 13a,b, for k = ± 1 and ± 2 respectively. It can be found that the distributions of ΔE v max,k are symmetric about the oxz plane. When only considering the Fourier series with k = ± 1, the amplitude of the TCP displacement can be as high as 2 × 10 −3 m in certain positions where the first natural frequency of the robot (related to deflections of Joint 2 and Joint 3) is close to ω v . When only considering the Fourier series with k = ± 2, the amplitude of the TCP displacement can be as high as 2.3 × 10 −3 m in the positions where the second natural frequency (related to the deflection of Joint 1) is close to 2ω v . It is also observed that when k > 4, kω v is beyond the robot's natural frequencies, and the values of ΔE v max,k are much smaller than those for k = ± 1, ± 2, ± 3 and ± 4. Now, we proceed to analyze the potential risk of the tool vibration when the percussive force acts in different directions on the robot end-effector. Once the maximum displacement-force frequency response ratio η vmax is given, then the maximum amplitude of TCP displacement due to the kth Fourier series of the periodic impact force ∆E vmax,k can be easily calculated using (24). For example, the values of ∆E vmax,k are plotted over the potential workspace when x is 0.5, 0.6, 0.7, 0.8, 0.9 and 1 m in Figure 13a,b, for k = ± 1 and ± 2 respectively. It can be found that the distributions of ∆E vmax,k are symmetric about the oxz plane. When only considering the Fourier series with k = ±1, the amplitude of the TCP displacement can be as high as 2 × 10 −3 m in certain positions where the first natural frequency of the robot (related to deflections of Joint 2 and Joint 3) is close to ω v . When only considering the Fourier series with k = ±2, the amplitude of the TCP displacement can be as high as 2.3 × 10 −3 m in the positions where the second natural frequency (related to the deflection of Joint 1) is close to 2ω v .
It is also observed that when k > 4, kω v is beyond the robot's natural frequencies, and the values of ∆E vmax,k are much smaller than those for k = ±1, ±2, ±3 and ±4.  A detailed investigation on the eigenvectors of matrix A v (kω v ) can give force components under which the maximum amplitude of TCP displacement due to the kth Fourier series occurs. The evaluation of ∆E vmax,k calculated from (24) is conservative in this case, since the three force components in the riveting force (25) have identical phases which can be insufficient to meet the requirement for the phases of the components in the eigenvectors of A v (kω v ). Note that the kth Fourier series cannot excite multiple vibration modes simultaneously since different mode shapes are orthogonal to each other, therefore. As a result, the maximum values of ∆E vmax,k for k = ±1, ±2, . . . , ±20 are drawn over the operating workspace in Figure 14. This gives an efficient evaluation of the maximum amplitudes of TCP displacement due to the periodic impact force in (25). The robot users should consider whether the resulting tool misalignment is acceptable for process quality requirement within the operating workspace when developing the robotic system.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 17 components in the riveting force (25) have identical phases which can be insufficient to meet the requirement for the phases of the components in the eigenvectors of A v (kω v ). Note that the kth Fourier series cannot excite multiple vibration modes simultaneously since different mode shapes are orthogonal to each other, therefore. As a result, the maximum values of ΔE v max,k for k = ±1, ± 2, …, ± 20 are drawn over the operating workspace in Figure 14. This gives an efficient evaluation of the maximum amplitudes of TCP displacement due to the periodic impact force in (25). The robot users should consider whether the resulting tool misalignment is acceptable for process quality requirement within the operating workspace when developing the robotic system. Figure 14. Maximum value of ΔE v max,k over the operating workspace.

Conclusions
In this paper, displacement-force and displacement-torque frequency response ratios are defined for flexible-joint robots, which represent the mapping from the amplitudes of external harmonic force and torque to the amplitude of tool displacement respectively. The upper bounds of these ratios are used as indices for structural dynamics performance of the robot. These indices give the worst situation of the tool displacement due to harmonic excitation with amplitude of unit force and unit torque respectively. The development leads to an efficient method to evaluate the maximum tool misalignment due to periodic external force and torque.
Numerical simulation is carried out for a robotic system for aerospace assembly with large dynamic loading. It is found that the dynamics performance indices dramatically change with the excitation frequencies for a given robot configuration. They reach peak values at frequencies close to the natural frequencies of the robot. These peak values show that only considering the elastostatic performance is insufficient for dynamic processes. An eigenvector analysis can provide the amplitudes and phases of the force and torque components under which the worst situation of the tool vibration occurs. For a given process, the maximum tool misalignment can be easily evaluated by combining the robot-dependent indices and the process-dependent Fourier series, so that the robot users can consider whether the forced vibration behavior of the robot is acceptable for process quality requirement when developing the robotic system.
Although the structural dynamics model in this work is derived for flexible-joint robots and the case study is given for manufacturing applications, the two performance indices are applicable for other robots and other applications, once the frequency response functions are available, for example,

Conclusions
In this paper, displacement-force and displacement-torque frequency response ratios are defined for flexible-joint robots, which represent the mapping from the amplitudes of external harmonic force and torque to the amplitude of tool displacement respectively. The upper bounds of these ratios are used as indices for structural dynamics performance of the robot. These indices give the worst situation of the tool displacement due to harmonic excitation with amplitude of unit force and unit torque respectively. The development leads to an efficient method to evaluate the maximum tool misalignment due to periodic external force and torque.
Numerical simulation is carried out for a robotic system for aerospace assembly with large dynamic loading. It is found that the dynamics performance indices dramatically change with the excitation frequencies for a given robot configuration. They reach peak values at frequencies close to the natural frequencies of the robot. These peak values show that only considering the elastostatic performance is insufficient for dynamic processes. An eigenvector analysis can provide the amplitudes and phases of the force and torque components under which the worst situation of the tool vibration occurs. For a given process, the maximum tool misalignment can be easily evaluated by combining the robot-dependent indices and the process-dependent Fourier series, so that the robot users can consider whether the forced vibration behavior of the robot is acceptable for process quality requirement when developing the robotic system.
Although the structural dynamics model in this work is derived for flexible-joint robots and the case study is given for manufacturing applications, the two performance indices are applicable for other robots and other applications, once the frequency response functions are available, for example, through vibration tests. Therefore, a future research topic is the validation experiments for the proposed method on a real prototype.