Modal Kinematic Analysis of a Parallel Kinematic Robot with Low-Stiffness Transmissions

: Several industrial robotic applications that require high speed or high stiffness-to-inertia ratios use parallel kinematic robots. In the cases where the critical point of the application is the speed, the compliance of the main mechanical transmissions placed between the actuators and the parallel kinematic structure can be signiﬁcantly higher than that of the parallel kinematic structure itself. This paper deals with this kind of system, where the overall performance depends on the maximum speed and on the dynamic behavior. Our research proposes a new approach for the investigation of the modes of vibration of the end-effector placed on the robot structure for a system where the transmission’s compliance is not negligible in relation to the ﬂexibility of the parallel kinematic structure. The approach considers the kinematic and dynamic coupling due to the parallel kinematic structure, the system’s mass distribution and the transmission’s stiffness. In the literature, several papers deal with the dynamic vibration analysis of parallel robots. Some of these also consider the transmissions between the motors and the actuated joints. However, these works mainly deal with the modal analysis of the robot’s mechanical structure or the displacement analysis of the transmission’s effects on the positioning error of the end-effector. The discussion of the proposed approach takes into consideration a linear delta robot. The results show that the system’s natural frequencies and the directions of the end-effector’s modal displacements strongly depend on its position in the working space.


Introduction
Parallel robots are Parallel Kinematics Machines (PKM) which have been investigated for a long time, as demonstrated by some key papers [1][2][3]. The parallel kinematic configuration gives to these robotic systems a high stiffness/inertia ratio. This characteristic is due to closed kinematic chains, where some members are mainly subject to axial loads; hence, these links can be designed to be slender and light without the risk of decreasing the system's overall stiffness. Moreover, the parallel kinematic configuration allows the positioning of the driving motors on the fixed frame of the robot, reducing, therefore, the moving masses [4,5].
The linear delta kinematic configuration [3] is characterized by linear axes instead of the rotary ones used within the well-known classical Delta configuration [1,2]. As a result, for this configuration, the working volume is spatially extended along a direction whose length depends only on the length of the linear axes. Applications in the industrial field requiring high speed and working volume extended along a direction can use linear delta robots designed with a suitable transmission that on one hand permits high linear speed and on the other has low inertia. This configuration is suitable for the primary packaging station of a production line or for transferring production goods between two lines. Linear transmissions based on a timing belt are an interesting solution due to the low mass, high speed, long travel length and low cost, especially if compared with other kinds of linear transmissions, such as ball-screws.

The Parallel Robot under Analysis
The parallel robot on which we have applied our "Modal Kinematic Analysis" approach is a 3 d.o.f. parallel kinematic machine produced by the Italian company Mechatronics and Dynamic Devices s.r.l. [31] that is characterized by the typical linear delta kinematic configuration. Figures 1 and 2 present the linear delta, whose main characteristics are the following: The distance l r between axes is 200 mm. • The length l d of the links connecting the carriages to the end-effector is 400 mm.
The constraints between the links and, respectively, the carriages and the mobile platform are realized by means of universal joints, highlighted in blue for the carriages and in red for the mobile platform in Figure 3. Hence, each link is composed of two rods realizing a four bar linkage with a parallelogram configuration.

Kinematics and Dynamics of the Parallel Kinematic Part
This section presents the model for the kinematics and the dynamics of the parallel kinematic part of the manipulator. This part has very high stiffness and low weight due to the carbon fiber composite material of the rods. For these reasons, it is possible to say that the flexibility of the parallel kinematic part can be considered negligible with respect to the compliance of the linear belt transmissions. This subsystem can be thus modeled as being composed of rigid bodies, whereas the only flexible elements are the transmissions. Moreover, the belt's mass is negligible with respect to the carriages mass; accordingly, we neglect the inertial contributions of the belt, and simply model its stiffness using configuration-dependent lumped parameters.
The rigid structure of the manipulator is responsible not only for the dynamic coupling between the linear axes, but also for the overall mass distribution of the system; the rigidbody kinematics and dynamics of this subsystem-composed of the mobile platform, of the distal links and of the trucks-must, therefore, be thoroughly accounted for.
The kinematics of the manipulator are largely analogous to those of the Delta robot designed by Reymond Clavel and analyzed, for example, in [1,2]; as such, this mechanical structure results in a mobile platform with three translational degrees of freedom [3]. Figure 2 shows the reference frame with respect to which the end effector's position is The angular positions ϕ j , j = 1, 2, 3 of the universal joints that connect the distal rods to the platform are also shown in Figure 3. The free coordinates x 1 , x 2 and x 3 that govern the kinematics of the manipulator are those relative to the longitudinal displacement of the three trucks. Given the geometric parameter l r that defines the distance between two adjacent linear guides, the truck positions e j are trivially expressed as functions of x j : The direct and inverse functions between these quantities and the center point p of the moving platform must furthermore be determined. The chief relationships that allow to do so are the vector loop equations written for each kinematic chain: In Equation (2), l d,j is the vector pointing from the j th truck e j to the position c j of the j th platform universal joint constraint, while the difference (p − c j ) is a constant vector entirely determined by the geometry of the moving platform. The vector loop equations, whose geometric meaning is depicted in Figure 4, must be satisfied for every assembled configuration of the robot. To exactly enforce them when the truck position vectors e j are known, it is sufficient to intersect three spheres of radius l d and centered in the points, indeed, the intersection point coincides with the admissible p, from which the constraint positions c j can be then found. On the other hand, if p is given the points c j can be immediately determined; the intersection of the sphere of radius l d and center c j with the segment generated by the j th linear guide allows then the determination of e j (and x j ) that satisfy the vector loop equations.
Other positional quantities of interest are the centers of mass of the distal linkages and the Euler angles describing the rotation of the linkage rods. Since the two rods of each linkage are constrained to perform the same motion, the kinematics of the entire distal link can be abstracted to those of its center-line.
Given that each distal rod is a slender, axially symmetric body, the position d j of the center of mass of the j th linkage can be expressed as a weighted average of the constraint and carriage position.
where W c is the weight of the universal joint and W t is the weight of the truck. A physically meaningful description of the distal link rotation must take into account the actual configuration of the universal joint connecting each rod pair to the end-effector. With reference to Figure 5, where the red dashed line represents the direction of the generic link j, the most straightforward solution involves a joint allowing a rotation of the rod pair first around the local x-axis by an angle α j , and then around the rotated local y -axis by an angle β j . The third Euler angle γ j around the local z -axis, directed as the distal rod itself, is kinematically constrained to be equal to zero. Given the vector (5) expressed in the j th frame of reference, these angles are easily computed as: Additionally, the rotation matrix R [j] d,j of the distal link can be constructed from the Euler angles through standard formulas for the composition of elementary rotations.
As the position kinematics of the platform are described as a sequence of geometric rather than analytical operations, the Jacobian analysis must also be developed in the same vein. The first goal is to determine the Jacobian matrix relating the time derivative of the truck coordinates to the velocity of the mobile platform.
To do so, the formula for the velocity kinematics of rigid bodies can be applied to the three kinematic chains connecting the platform to the base. Recalling that the platform can only translate and thereforeċ j =ṗ the following equations can be written: Robotics 2021, 10, 132 The yet unknown distal angular velocities ω j can be eliminated from Equations (9) and (11) by multiplying each side by the quantity (p − e j ) . In matrix notation: The truck velocities can then be expressed in terms of the time derivatives of the free coordinates by highlighting the Jacobian matrices D e,j of e j with respect to whereẋ = ẋ 1ẋ2ẋ3 . By substitution of Equation (13) in Equation (12), the following terms can be highlighted: More explicitly, The platform Jacobian matrix can finally be expressed as D p = D l −1 D r . Each Jacobian matrix related to the motion of the center of mass of the j th distal linkage is then easily derived: Given that both positions and velocities of the notable points of the structure are made explicit, it is also possible to determine analytically the time derivative of the platform Jacobian matrix as a function of x andẋ: To computeḊ r andḊ l , it is enough to differentiate, with respect to time, each element of D r and D l , which can be easily done thanks to the clear geometric meaning of these matrices.
The Jacobian analysis of the distal link rotations might then be performed. Equations (6)-(8) can be expressed in vector form as follows: where ψ j = α j β j γ j . The function f can be easily differentiated with respect to the components of u j to yield the matrix D f . On the other hand, the Jacobian matrix of u j with respect to x is where R f ,j is the constant rotation matrix due to ϕ j . It follows thaṫ A linear relationship ω [j] j = Ψ jψj subsists between the angular velocity of the distal rod and the time derivative of its Euler angles, with Ψ j being defined as As a result the angular velocity of the distal rods can be expressed as Finally, in the absolute frame of reference, Since the kinematics of the rigid subsystem is expressed as a function of x andẋ, its kinetic energy can be written as The mass matrix M ss of the manipulator, which appears inside the parentheses of Equation (24), is thus composed of terms related to the translations of the mobile platform; the translations of the trucks; the translations of the center of mass of each distal linkage; and the rotations of the distal links. Accordingly, the masses m t , m d and m p of the trucks, the distal linkages and the platform appear alongside the inertia matrix I d,j of the distal links. These quantities are set for the following numerical analyses as Assuming W t = W c , the distal link inertia is calculated in the principal and barycentric frame as Conversely, the gravitational potential of the subsystem is related only to the mobile platform and to the center of mass of the distal links, as only these are allowed to move along the vertical direction. Indicating as g the gravitational constant, the gravitational potential U g can be written as: The gradient of U g can be computed straightforwardly as follows: Robotics 2021, 10, 132 9 of 22

Belt Transmission Dynamics
As already mentioned, the mass of the belts can be neglected; therefore, the kinetic energy associated to the belt transmission is determined by the rotational inertia not only of the pulleys, but more significantly of the three motors and gearboxes, which are assumed to rotate rigidly with the actuated pulley.
To properly describe the system, six additional coordinates ϑ j,1 ϑ j,2 (with j = 1, 2, 3) associated to the rotation of the pulleys are required. Figure 6 shows a diagram of the timing belt in which the stiffnesses and the coordinates are highlighted for the j th axis. Which pulley is driven by the motor is shown in Figure 2. Denoting as I p,j,k the rotational inertia of the pulley described by the coordinate ϑ j,k , the kinetic energy of the three transmissions can be expressed in the form By introducing the array of pulley coordinates and by suitably collecting the mass terms, the kinetic energy of the transmissions can be expressed in matrix notation as In Equation (30) the transmission mass matrix M t is a 6-by-6 diagonal and constant matrix. For the subsequent analyses, pulleys having the following characteristics are considered: Inertia around the rotation axis I p,j = 130 × 10 −6 kg · m 2 . Furthermore, the motor and gearbox inertia-suitably projected on the machine side through the gearbox reduction rate-is added to the inertia of the actuated pulleys.

Overall Dynamics
The total kinetic energy of the system can then be expressed as: Introducing the array of the free coordinates of the overall system as the complete mass matrix can be easily assembled: The overall mass matrix is a function of the trucks' positions. In a similar way, the mechanical stiffness of each transmission changes with the position of the truck, as this defines the free lengths of the three segments of which the timing belt is composed. As already highlighted in Section 2.1, the specific stiffness k sp quantifies the stiffness of the timing belt for each unit length and unit width of the belt itself, and as such, it is one of the parameters commonly specified within the product datasheet. For a given configuration, the stiffness coefficients of the belt sections are therefore computed according to the following equation: in which w b is the belt width, and l f ,j,k is the free length of the considered belt segment. More explicitly, the stiffness coefficients for each belt transmission are computed as where x * j represents the undeformed j th truck coordinate, l t is the length of the truck, and min(x j ) and max(x j ), shown in Figure 6, are fixed structural parameters that correspond to the maximum and minimum displacements achievable by an ideally dimensionless cart.
If large rigid motions of the system are considered, the elastic actions do not admit a potential function, due to the non-constant stiffness coefficients. However, if small displacements around a given undeformed configuration are to be investigated, the stiffness variations can be neglected, and an elastic potential for each belt can be determined as Here, the variables ∆ϑ j,1 , ∆ϑ j,2 and ∆x j are to be interpreted as displacements around the undeformed configuration ϑ * j,1 , ϑ * j,2 , x * j : The potential function of the free system can then be obtained by summation of the gravitational and elastic potentials: The dynamic model describes a semi-definite system, which therefore can undergo rigid motions. Indeed, the motors apply torques to the system and do not set the position of the actuated pulleys. The equilibrium condition of the system, for a given position of the end-effector, can, however, be evaluated by introducing static torques τ st,1 , τ st,2 , τ st,3 applied to the actuated pulleys in order to keep the end-effector in the desired configuration.
The external generalized forces acting on the generalized coordinates of the system can therefore be expressed in vector notation as It should be noted that the vector Q defined in (43) is consistent with the location of the motors reported in Figure 2, and thus properly distinguishes the actuated degrees of freedom (namely ϑ 1,1 , ϑ 2,2 and ϑ 3,2 ) from the passive ones. In more detail, the position of the static torques within vector Q corresponds to the position of the driven pulleys' rotational angles in vector q. The other elements of vector Q are equal to zero because they correspond to the position of the non-actuated coordinates.
A static equilibrium configuration should satisfy the condition The solution of the set of Equation (44) yields the static equilibrium configuration around which we develop the modal analysis. In particular, the static torques and the pulleys' rotations are the unknowns of the problem, while the displacements of the trucks, which alone define the position of the mobile platform, are fixed. This choice is motivated by the need to express the several results of the modal analysis as functions of the mobile platform coordinates, as this can shed additional insight also on the functional characteristics of the robotic device.

Configuration-Dependent Modal Analysis
The modal analysis around a given equilibrium configuration q eq can be performed by substituting into Lagrange's equations the approximated kinetic energy and potential functions.
In particular, where H U is the Hessian matrix of the total potential. Lagrange's equations assume the form Remembering the equilibrium conditions, Equation (47) can be written in the familiar matrix notation, in which the stiffness matrix is simply defined as K = −H U . It might be observed that the Hessian matrix associated to U el, can be analytically obtained by double differentiation, while the one relative to the gravitational potential U g is computed using the following property of the time derivative of the Jacobian matriceschiefly D p -that govern the gradient ∇U g , EvaluatingḊ p (x, δ jk ), with δ jk being Kronecker's Delta, the partial derivatives ∂ ∂x j D p can be recovered, from which, in turn, the Hessian matrix of the gravitational potential can be straightforwardly constructed. The mass and stiffness matrices are constant for each position q eq , taken into account inside the investigated workspace. For each position, an eigenvalue problem for the matrix M −1 K is then set up and solved in order to highlight the configuration-dependent natural frequencies and modal vectors of the system. As expected, the first three modes represent the rigid motions of the system, while the remaining six are proper vibration modes.
Although a detailed workspace analysis of the manipulator is outside the scope of this work, the investigation refers to a significant working plane selected by the Jacobian's condition number of the rigid subsystem. The Jacobian's condition number of the chosen plane is more uniformly distributed and has the lower mean value, i.e., 3.5. Figure 7 shows the selected plane at the position z = −245 mm. An equally distributed sampling of this plane defines the investigated end-effector's pose, giving the equilibrium positions q eq . 245 working plane working volume

Results and Discussion
This section presents the results of the modal analysis developed on the end-effector's positions belonging to the plane defined in the previous section. The software package selected for the implementation of the system's model, and for the development of the modal analysis for different end-effector's positions is MathWorks ® MATLAB ® , owing to its matrix processing and data visualization capabilities. The first results reported are the natural frequencies for each of the six modes of vibration (the remaining three, related to the rigid motion of the system, are not included); f 1 − f 6 are the natural frequencies (expressed in Hz) of the modes 1 to 6, sorted in ascending order. In particular, Figure 8    This rather clean distinction between lower and higher frequency modes is also reflected As expected, all the graphs are symmetric with respect to the system's central linear axis (axis 1), while the change of the natural frequencies along the y axis depends on the specific mode of vibration. Regarding the first mode, the graph shows that the frequency has its maximum value located in the lower half of the reference plane; starting from this point, moving toward the upper or the lower end of the working plane (i.e., increasing of decreasing coordinate y p ), the frequency decreases, reaching its minimum value at the upper end. For the second mode, minimum and maximum values correspond, respectively, to the upper end and to the lower end of the working area. Moving toward lower y p coordinate values, the frequency increases with a gradient more significant in the lower half of the area. The third mode has the frequency's minimum value at approximately the middle of the working area, and the frequency increases, both increasing and decreasing coordinate y p . The gradient of the frequency becomes quite significant near the ends of the area. For the fourth mode, the behavior is almost symmetric also along the y direction, with the maximum value at the upper end and the minimum value in the middle. The fifth mode is characterized by a progressive increase of the frequency from the lower end to the higher one. The sixth mode has the minimum value located approximately in the middle of the area, and the maximum value at the lower end. Again, the frequency increases, both increasing and decreasing the position y p , but with a gradient higher in the upper half rather than in the lower one.
However, the information concerning the natural frequencies is not enough to completely understand the system's behavior. As a matter of fact, even if we know the frequency of vibration corresponding to a specific mode, we do not know how the system vibrates according to the mode. So we do not know which is the influence of that specific mode on the behavior of the end-effector. In order to reach this goal, the attention must be focused on the modal vectors and on their projection on the end-effector's coordinates; this is the key point. The values of the modal vectors for the positions corresponding to the minimum and maximum values of the frequency, taken along the symmetry axis of the plane within the workspace area not affected by boundary effects, are summarized for each mode of vibration in Tables 1 and 2. The truck translations are divided by the pulley radius to compare non-homogeneous linear and angular coordinates fairly.  In both Tables 1 and 2, it can be seen that the motion of the free pulley is dominant within the modal vectors of modes 4, 5, and 6; on the other hand, the component associated to the actuated pulleys is the most significant for modes 1, 2 and 3. This rather clean distinction between lower and higher frequency modes is also reflected in the magnitude of the displacements of the trucks; indeed the maximum absolute truck displacement over modes 1, 2 and 3 is at least one order of magnitude greater than that of modes 4, 5 and 6.
The kinematic modal analysis results are also represented from the working space point of view, highlighting, for different positions on the reference plane, the directions of the end-effector's displacement as a consequence of the projection of the modal vector components related to the carriages' displacement. In other words, they represent the projection of the modal vectors after suitable normalization on the mobile platform using the Jacobian matrix D p . The direction of the end-effector's displacements are represented, for modes 1 to 6, in Figures 9-14. of modes 4, 5 and 6.       From the analysis of these diagrams, it is evident that the direction of vibration of the end-effector, due to a specific mode of vibration, strongly depends on its position within the reference plane. For all the modes, the displacement directions are symmetric with respect to the central axis, and they progressively change by changing the y p coordinate. Figures 15-20 represent, in another way, the same modal vectors, focusing on values of in-plane and out-of-plane displacements rather than on the directions of vibration, using the same modal vector normalization for all the figures.  As regards the second mode of vibration, Figure 16 shows an opposite behavior 509 with respect to Figure 15, for both the in-plane and the out-of-plane displacements.

510
The in-plane displacement decreases as y P coordinate increases, with the minimum This approach highlights the first-order effects of the truck displacements on the displacement of the end-effector. The in-plane component was defined and represented as δp 2 1 + δp 2 2 , where δp 1 and δp 2 are the first two components of the modal vector projected on the end-effector, while the absolute value of the vertical out of plane component δp 3 is represented separately. The vibrational displacements of the mobile platform are rep-resented using the same scale, allowing a quantitative comparison between the effects of different modes.
In particular, Figures 15-17 are related to the first three modes and show that, for these modes, the direction's displacements of the end-effector, both along the reference plane and along the vertical direction, are significant (please note that the color bars of the graphs have all the same scale). Moreover, the shape of the diagrams is clearly symmetric with respect to the central axis. As far as the y direction is concerned, Figure 15, which refers to the first mode of vibration, shows that the upper half of the working plane, the one nearer to motors 1 and 2, is characterized by values of in-plane displacements lower than the ones corresponding to the lower half-plane, the one nearer to motor 1. Hence, displacements increase as the y p coordinate increases. For the out-of-plane components graph, the behavior is opposite: as y p coordinate increases, the displacement decreases. The displacement's maximum value is located in a small area just below the half of the plane. Moreover, it should be noted that along the symmetry axis, there is a position near to the middle of the area where the vertical displacement has a sudden change.
As regards the second mode of vibration, Figure 16 shows an opposite behavior with respect to Figure 15, for both the in-plane and out-of-plane displacements. The in-plane displacement decreases as the y p coordinate increases, with the minimum value located just above the central position, and the maximum one located in small areas on opposite sides just below the center of the plane. On the contrary, the vertical displacement increases as the y p coordinate increases, and there is an area around the middle of the plane where it has its maximum value and where a sudden change occurs.
For the in-plane displacement related to the third mode of vibration (Figure 17), the working plane is again divided in two equal parts: the lower one is characterized by values lower than the upper part. Moreover, no sudden changes are highlighted, just small areas where the values slightly change with respect to the values of the respective half area. The out-of-plane component is quite uniformly distributed along the plane, with a central area where the values increase and some areas are located approximately near the ends of the plane's horizontal and vertical middle lines where the values decrease.
As regards the displacements related to the fourth, fifth and sixth modes, respectively depicted in Figures 18-20, the displacements, both in-plane and out-of-plane, of the endeffector are clearly negligible. This confirms on the entire working plane that-as already seen for some notable configurations analyzed in Tables 1 and 2-the truck displacements associated to modes 4, 5 and 6 are comparatively smaller, and thus do not generate appreciable motions at the end-effector.

Conclusions
The proposed approach is applied to a parallel kinematic manipulator with driven joints characterized by low mechanical stiffness. It outlines the effects of the stiffness of the transmission, of the mass distribution and of the coupling between the joints due to the parallel kinematic chains on the end-effector's vibration direction and magnitude. The mathematical approach considers the masses of all the elements of the parallel kinematic part, the system's actual configuration, i.e., that it is subject to a motion control algorithm, where the motors are torque controlled, and the transmissions' stiffness, leading to a 9 d.o.f. system. The direction of vibration of the end-effector, which depends on the dynamic configuration of the system, is represented by the modal vector calculated from the equations of the system's model, linearized around the investigated points of the workspace.
The discussed results show that the magnitude and direction of the modal displacements at the end-effector, evaluated in suitable points of the working space, are influenced by the configuration-dependent mass distribution and transmission's stiffness.
The method, applied to a linear delta manipulator, highlights that the effects of the different vibration modes can be effectively compared, considering the amplitude of displacement of the end-effector in the x − y plane or z-direction.
Moreover, the discussion of the results outlines that the displacement's direction of the end-effector changes along the working area as a function of the kinematic and transmission configuration. For the same vibration mode (i.e., the first), in some zones of the working space, the end-effector is subject to vibration in the vertical direction, while in other zones, the displacement's direction is predominantly horizontal. In the same way, the frequency associated with the mode varies in the workspace. The proposed approach constitutes a useful support for the system's design in evaluating the end-effector vibration direction for a given vibration mode, in choosing the working area within the plane, in selecting the proper motion laws, and in synthesizing the control system.
Finally, the proposed approach introduces, in an effective and computationally efficient way, the mass distribution and the coupling effects of the parallel kinematic part as a function of the investigated end-effector positions, thanks to the used mathematical approach. Lastly, our work could open the way toward its implementation in the cases of robots where the stiffness of the parallel kinematic part cannot be neglected.