On the elastodynamics of a five-axis lightweight an- thropomorphic robotic arm

This paper presents elastodynamic modeling and analysis for a five-axis lightweight robotic arm. Natural frequencies are derived and visualized within the dexterous workspace to show the overall performances and compare them to the frequencies when the robotics is with payload. The comparison shows that the payload has a relatively small influence to the firstand second-order frequencies. Sensitivity analysis is conducted, and the system’s frequency is more sensitive to the second joint stiffness than the others. Moreover, observations from the displacement response analysis reveal that the robotics produces linear elastic displacements of the same level between the loaded and unloaded workingmodes but larger rotational deflections under the loaded working condition. Themain contribution of this work lies in that a systematic approach of elastodynamic analysis for serial robotic manipulators is formulated, where the arm gravity and external load are taken into account to investigate the dynamic behaviors of the robotic arms, i.e., frequencies, sensitivity analysis, and displacement responses, under the loaded mode.


INTRODUCTION
Lightweight robotic arms and anthropomorphic assistive robots with high payload capacity are desired for applications of industry and welfare, among other fields, such as assisted daily living [1][2][3] , pick-and-place operations [4] , etc. Pick-and-place robots are well suited for a static environment where the task is repeated and precise tolerances are demanded [5] . As a mechanical system, the dynamic characteristics of the robotic arm is of importance to account for the requirements of application, such as high precision, speed, and payload. Henceforth, higher natural frequencies and low elastic displacements of a robotic manipulator will allow higher operational speeds and working cycles for efficient productivity [6] . Natural frequencies indicate the condition in which a mechanism tends to vibrate [7,8] . Differing from a structure or element, the dynamic behavior of a mechanism usually heavily depends on its architecture and configurations [9] ; thus, it is not a trivial task to characterize the robot dynamics throughout the workspace, which calls for the kineto-elastodynamic analysis to provide the fundamentals of the modeling, design and control.
The elastodynamic modeling and analysis of a robotic manipulator have been reported previously [10,11] , and they are roughly grouped into two categories: lumped modeling [12][13][14][15] and distributed-flexibilities modeling [9,[16][17][18][19] . In general, with lumped modeling it is simpler to model the elastodynamic equation with acceptable computational accuracy, while the latter provides a more accurate model but with the high-dimensional generalized coordinate space and more complex procedure [20] . The commonly used method to study the elasticity of the robotic manipulators is the virtual joint method (VJM) as it can provide acceptable computation accuracy that is close to that of finite element analysis (FEA) [21] . Besides, VJM can be time efficient. VJM is based on pseudo-rigid body models with "virtual joints" [22][23][24][25] . Generally, the link flexibilities and linear/torsional springs take into account the bending contributions to the mechanism [26][27][28][29] . The stiffness formulated in the above approaches is limited to a subspace defined by the degrees of freedom (dofs) of the manipulator end-effector. Pashkevich et al. [30] overcame this issue by introducing a full-mobility lumpedparameter model by localizing 6-dof virtual springs to the links' ends and/or joints. In these models, the stiffness matrix is calculated in an unloaded equilibrium configuration of a robotic manipulator. On the other hand, the external loads directly influence the manipulator equilibrium configuration and, consequently, may modify the static properties. The lightweight design of the robotics accordingly decreases the link structural stiffness; thus, the robot geometry change due to external loads should be considered [31][32][33] . Consequently, elastodynamics of the robotic manipulators is an important concern in their design and applications. Based on the matrix structural analysis, Cammarata et al. [9,34] proposed an algorithm to assemble the stiffness matrix to investigate the manipulators with lower kinematic pairs. In this manner, the overall robotic manipulator inparallel architecture can be split into substructures for modeling the elastodynamics [35] . Wu et al. [36] analyzed and compared the stiffness and natural frequencies of a 3-dof parallel manipulator with/without a redundant leg, where the joint deformations are ignored in the stiffness modeling. The small-amplitude deformations of the active joints can be considered as parameter uncertainties in terms of small variations to be integrated into the dynamic model [37] . Briot and Khalil [14] used the Newton-Euler recursive approach to develop a general symbolic elastodynamic calculation model for flexible parallel robots. Taghvaeipour et al. [15] derived the posture-dependent stiffness matrix in the elastodynamic modeling by resorting to the generalized spring concept. The previous models were established in the nominal configurations; hence, the geometry changes of the manipulator in this work are considered in the elastodynamic modeling and analysis.
In this paper, the elastodynamic characteristics of a lightweight robotic arm are investigated. The arm gravity and external load are taken into account to derive the stiffness matrix. Isocontours of natural frequencies over the dexterous workspace are formulated and sensitivity analysis is conducted. The frequencies and displacement responses of the robotics with payload are analyzed and compared with the dynamic behaviors of the unloaded mode. The main contribution of this work lies in that a systematic approach of elastodynamic analysis for serial robotic manipulators is formulated, where the arm gravity and external load are taken into account to investigate the dynamic behaviors of the robotic arms, i.e., frequencies, sensitivity analysis, and displacement responses, under the loaded mode.  [38] .

KINEMATICS OF THE LIGHTWEIGHT ROBOTIC ARM
The lightweight robotic arm under study has five degrees of freedom (dof) [38] , which adopts a modular design approach, as shown in Figure 1. The revolute joints are composed of CPU series gearboxes of Harmonic Drive and Maxon motor with gearhead to enhance the torque capabilities, except Joint 4 with geared motor. The actuators of joints are controlled by Maxon EPOS controllers. The Controller Area Network (CANopen) bus is adopted to build the communications between motors and controllers, and A CAN-USB interface is used to establish the communications between CANopen bus and the PC [38] . In accordance with the Denavit-Hartenberg (D-H) convention [39] , the Cartesian coordinate systems are established accordingly.

Kinematics of robotic arm
Throughout this work, i, j, and k stand for the unit vectors of the -axis, -axis, and -axis, respectively. The transformation matrix in forward kinematics of the end-effector in reference frame is expressed as where D-H parameters are given in Table 1, and the inverse geometry problem for this robotics is well documented in the literature [8] .

Kinematic jacobian matrix
The velocities between the joints and end-effector are mapped with the Kinematic Jacobian matrix  is the velocity of the end-effector. Moreover, J is the kinematic Jacobian matrix of the robotic arm [40] , namely, where R −1 and q −1 denote the rotation matrix and position vector of the transformation matrix from the reference coordinate system to the ( − 1)th coordinate system, respectively, which can be extracted from

Dexterous workspace
The reachable workspace of the robotic arm can be visualized by considering the limitation of the joint displacements and link dimensions. To effectively perform the kinematic performance, a dexterous workspace is defined, throughout which the inverse of the condition number of the Jacobian matrix is greater than 0.2, namely −1 (J) ≥ 0.2. Since the Jacobian matrix of Equation (4) is not homogeneous, a characteristic length [41] is introduced to normalize the Jacobian matrix as follows: By constraining the condition number of the kinematic Jacobian matrix, a regular dexterous workspace is quarterly visualized in Figure 2.

ELASTODYNAMIC MODEL OF ROBOT
The elastodynamic modeling procedure pertains to the calculation of the stiffness and mass matrices of the manipulator, which is described in the following sections. Prior to the derivation of the elastodynamic model, the following assumptions are made: • The actuator stiffness is considered as an 1-dof torsional spring, while the link is considered as cantilever with a 6-dof spatial spring located at the end but treated as rigid. • The centers of mass of the regular components are coincident with their geometric centers.
• The sum of moments of inertia of the actuators and Harmonic drivers are considered as lumped.

Stiffness matrix
To derive the elastodynamic equation of the robotic arm, the stiffness matrix is calculated with the virtual spring approach [42] , based on the screw coordinates [43] . Hence, the component masses and external loads are taken into account to compute the Cartesian stiffness matrix. Figure 3 shows the VJM model of the robotic arm, where G , = 1, 2, ..., 7, stands for the gravity and F for the external loads.
Let and be the original and the deformed displacements of the virtual springs, respectively, following the principle of virtual work, i.e., the work of the auxiliary forces is equal to the work of internal forces , namely, where the virtual displacements t and t can be computed from the linearized geometrical model derived from t = J ( − ) and t = J ( − ), respectively, J and J being the Jacobians, namely, where J (:, 1 : ) stands for the first columns in J and stands for the total degrees of freedom of the virtual springs from the base to G . Moreover, (7) is rewritten as consequently, the force equilibrium equation is derived as Assuming that K is the stiffness matrix in the joint space, with the linearized force-deflection relation, the equilibrium condition can be written as with where act, is the actuation stiffness and K and K are the upper and lower link stiffness matrices, respectively.
To calculate the stiffness matrix of the loaded mode, a neighborhood in the loaded configuration in which the external loads and the joint location are supposed to be incremented by small values F and , which can still satisfy the equilibrium conditions, is considered, leading to (16) and the linearized kinematic constraint Based on Equation (14), expanding Equation (16) yields where the symbol ⊗ represents the Kronecker product between matrices and H = J / , H = J / . Combining Equations (17) and (18), the stiffness model of the robotic manipulator is reduced to with From F = K t, the Cartesian stiffness matrix K of the robotic arm is calculated as

Mass matrix
The mass matrix can be derived from the expression of the system' s kinetic energy, consisting of energies of the revolute joints, links, and end-effector. The energy of the five active joints are and where , is the moment of inertia of the th joint, , is the mass, and v , is the velocity in the Cartesian space. Let I = diag[ ,1 , ,2 , ..., ,5 ]; then, Equation (22) can be written in a compact form, namely, with The kinetic energy of the upper/lower links and the wrist link can be expressed as where the subscripted I, , and v stand for the moment of inertia, mass, and velocities in the Cartesian space, respectively, and where q and q are the position vector of the centers of the mass of the upper and lower links, respectively. Equation (27) can be cast in a matrix form as follows: with Similarly, the kinetic energy of the end-effector can be obtained as where I is the moment of inertia of the end-effector and is the mass.
From the total kinetic energy of the robotic arm = + + , the mass matrix M for the robotic arm can be expressed as

Dynamic equation and analysis
The dynamic equation of the robotic arm can be formulated as where C is the damping matrix, F is the resultant force, and u and u are the elastic displacement and acceleration, respectively. Since damping can only slightly influence the natural frequency and mode of free vibrations, the damp can be ignored to determine the natural frequencies. Simplification of Equation (35) results in the linearized elastodynamic equation below The rigidity of the system may be represented by the natural frequency. The higher is the frequency, the higher is the stiffness. From Equation (36), we get where = /2 denotes the natural frequency.
The displacement response analysis can be carried out from Equation (35) based on the initial conditions Here, the damping ratios are set to = 6% according to the manipulator structure. From Equation (37), the displacement vector u can be represented in terms of the modal contributions, namely, where Q and are the modal matrix and the vector of the displacements in each mode, respectively. Consequently, Equation (35) can be rewritten as (41b) Since the mass and stiffness matrices in Equation (35) are time-varying, the common way to solve such a problem is to divide the motion period into extremely short intervals, where the stiffness and mass matrices are considered as constant in each interval. Let denote the complete motion period that is divided into intervals, namely, Δ = / . In the th time interval ∈ [ −1 , ], the equation of motion in the th mode is expressed as Thus, the th mode contributes to the displacement response [44] is with Hence, ( ) and ( ) can be solved as long as ( −1 ) and ( −1 ) are given, and where e is the th column of the modal matrix. The total displacement response is calculated by the following addition Consequently, the natural frequency and displacement response can be obtained with numerical calculations.

NUMERICAL SIMULATION
Elastodynamic characteristics of the robotic arm are investigated in this section. The properties of the robotics components are listed in Tables 2 and 3, respectively. Moreover, according to the output shaft of the gearbox, the actuation stiffnesses are calculated and set to , = 2 · 10 4 Nm/rad, = 1, ..., 5, and the link stiffness matrices given in Appendix A are derived by means of FEA with ANSYS [45] . The numerical simulation was carried out with Matlab.

Natural frequency
To effectively measure the overall performance of the robotic arm, the distributions of natural frequencies over the dexterous workspace in Figure 2   vertical and horizontal, respectively. It can be observed that the nonsymmetric distributions of the natural frequencies in Figure 5 are different from the symmetric ones in Figure 4. This is because the robot configurations are not axisymmetric about the vertical direction with the vertical end-effector, leading to different inverse kinematic solutions of such a 5-dof robotic arm, which are different from the axisymmetric robot configurations with horizontal end-effector. As the mass and stiffness matrices of the robot are configuration dependent, non-symmetric distributions of natural frequencies in Figure 5 occur. These two figures show that the first two orders of natural frequencies increase with the increasing coordinates but with decreased and coordinates, namely both the first and second frequencies increases from the workspace boundaries to the origin of the global coordinate systems. As displayed in Figure 4, when the end-effector remains vertical, the natural frequencies have the same varying trend in any vertical cross-section of the workspace. By contrast, the first-and second-order frequencies become smaller counterclockwise within the workspace when the endeffector is in the horizontal configuration, as shown in Figure 5. Moreover, it is found that the differences among the frequencies of the manipulator in different configurations are not so large, which means that the robotic arm has close frequencies inside the overall workspace.

Sensitivity analysis
Sensitivity analysis can be used to evaluate the influence of the geometric parameters and design variables to the manipulator performances. Based on the elastodynamic equation, there exists  Figure 6 illustrates the sensitivity of the first-order natural frequency to the first two active joints with constant orientation [0, /2, 0]. It is found that the first-order natural frequency is much more sensitive to the second joint, particularly in the upper and lower workspace regions, which implies that the robot' s dynamic performance can be improved by replacing the second joint with a stiffer actuator. It is noted that the distributions of sensitivity coefficients are not symmetric, which is because the robot configurations are not axisymmetric about the vertical direction when the robot end-effector moves with some constant orientations, since the robot under study is a 5-dof robotic arm. Moreover, if a payload with more mass were exerted to the robot, it could be predicted that the sensitivity coefficients will be increased with very tiny varying trends, compared to the present results.

Dynamic analysis of loaded system
With the payload 5 kg applied to the end-effector of the robotic arm, they constitute a new dynamic system and the solved frequencies with constant-orientation [0, /2, 0] are illustrated in Figure 7, from which it is observed that the frequencies of the loaded robotic system decrease about 20% compared to Figure 5. Table 4 lists the average frequencies [46] within the constant-orientation workspace defined bȳ where Ω stands for the workspace volume. Different from the traditional industrial robots with low frequencies, the high order frequencies have large values to make the manipulator achieve high-speed motion. Compared to the average natural frequencies, the frequencies of the robotics with payload reduce 10%-40% for the six orders of frequencies. From the view of kineto-elastodynamic characteristics, the difference between the frequency of the loaded system and its natural frequency could be a consideration in the design of the mechanical system, where the smaller difference implies higher rigidity and higher payload capability.  Assuming that the motion of the robotic arm follows the trajectory (unit: mm) defined by where the end-effector keeps constant-orientation [0, , 0] and the motion period = 0.5 s is divided into 1024 intervals, Figure 8 shows the displacement responses of the end-effector, from which it is seen that the linear elastic displacement responses are close, whenever the robotic arm is under loaded and unloaded working modes. The angular displacements of the end-effector generate relatively large differences. The largest deformations appear around 0.3 s where the end-effector is located in the middle layer of the workspace, approximately = 250 mm. Figure 9 shows the comparison of the joint angular displacements between the numerical simulation and ex-  perimental measurements along previous trajectory, where the experimental data are read from the motor encoders. Due to the frictions and time-varying disturbance in the joints, the experimental curve profiles have more fluctuations and larger vibration amplitudes than the simulation ones. On the other hand, the comparison shows that the differences between these two curves are small, thus, the built analytical model can be acceptable for dynamic analysis of the robots.

CONCLUSION
This paper presents the elastodynamic characteristics of a 5-dof lightweight robotic arm. The main contribution is that a systematic approach of elastodynamic analysis for serial robotic manipulators is formulated, where the arm gravity and external load are taken into account to investigate the dynamic behaviors of the robotic arms, i.e., frequencies, sensitivity analysis, and displacement responses, with auxiliary payloads exerted to the robot. The modeling in this work eases the evaluation of elastodynamics of the manipulator at a large number of postures as the elastodynamic aspect is usually time-consuming. As the mass and stiffness matrices are posture dependent, the proposed method can effectively provide a symbolic calculation and achieve the modal analysis along an operating trajectory. Moreover, such a model can compute the additional mass or evaluate the influence of an isolator to the system more precisely to eliminate/reduce vibration in the vibration control. The developed model can be used in either performance evaluation or design optimization.
The frequencies of the loaded robotics are visualized within the representative workspace regions to show the overall dynamic performance and compare them with the natural frequencies. The comparison reveals that the studied robot keeps relatively high rigidity with high payload ratio. It is found from sensitivity analysis that the natural frequency can effectively increase by improving the second joint stiffness. Based on the displacement responses analysis, the payload has a slight influence on the translational elastic displacements of this robotic system, although it leads to reduced frequencies, while the effect on the rotation deflections cannot be ignored.
In the future, the developed model will be integrated into its control system and an optimum redesign of the robotics will be conducted.

APPENDIX A: STIFFNESS MATRICES OF ARM LINKS
The stiffness matrices of the upper and lower links K and K for the 5-dof robotic arm, computed by means of finite element analysis (FEA) with ANSYS, are given as where the blocks corresponding to rotation, translation, and coupling terms are given in Nm/rad, N/rad, and N/m, respectively.