A Comprehensive Multibody Model of a Collaborative Robot to Support Model-Based Health Management

: Digital models of industrial and collaborative manipulators are widely used for several applications, such as power-efﬁcient trajectory deﬁnition, human–robot cooperation safety improvement, and prognostics and health management (PHM) algorithm development. Currently, models with simpliﬁed joints present in the literature have been used to evaluate robot macroscopic behavior. However, they are not suitable for the in-depth analyses required by those activities, such as PHM, which demand a punctual description of each subcomponent. This paper aims to ﬁll this gap by presenting a high-ﬁdelity multibody model of a UR5 collaborative robot, containing an accurate description of its full dynamics, electric motors, and gearboxes. Harmonic reducers were described through a translational equivalent lumped parameter model, allowing each constitutive element of the reducer to have its decoupled dynamics and mating forces through non-linear penalty contact models. To conclude, both the mathematical model and the real robot on a test rig were tested with a set of different trajectories. The experimental results highlight the ability of the proposed model to accurately replicate joint angular rotation, speed and torques in a wide range of operational scenarios. This research provides the basis for the development of a model-based PHM-oriented framework to carry out detailed and advanced analyses on the effects of manipulator degradations.


Introduction
Recent years have seen a significant growth in the adoption of collaborative robots (cobots) in industrial environments, since their low entry-point cost and relative simplicity has popularized the adoption of increased automation levels in small and medium enterprises [1]. Since these manipulators are largely employed in assembly lines [2], inspection [3], or repair operations [4,5], their malfunctioning might jeopardize the operators' safety and cause unexpected downtimes. To mitigate these risks, preventive actions, such as scheduled maintenance, are usually adopted. Nevertheless, since each robot is programmed to execute a specific task, its load history and degradation rate are not equal to the ones of the same machine performing a different task. Other possible solutions are represented by stand-by working stations, backup manipulators, and fault-tolerant control algorithms [6], which, however, are not always applicable.
Within this context, condition-based maintenance represents a key technology to ensure line availability and prevent economical losses. The ability to detect a fault at its early stage allows scheduling the replacement of the faulty components during periods of lower production, or before the robot can cause any damage to the operators with which the cobot shares its workspace or the production line itself. To do so, data-driven models (DDMs) are often used to extract health indicators from the machine and predict its remaining useful life using the obtained dataset to test machine learning algorithms [7]; for example, confronting simulated acceleration and measured acceleration by means of an accelerometer mounted on the tip of the robot [8,9]. While these algorithms can also be used to detect motor failures [10], according to [11], electric faults are likely to be autonomously detected by the built-in control logic of the joint. As such, prognostics and health management (PHM) studies on robotic systems are mainly focused on the detection and classification of mechanical faults on bearings [12] and gearboxes. By adding several sensors [13][14][15][16], it is possible to obtain information about, for example, speeds and accelerations. After the definition of the dataset, the analysis of these data can be performed using, for example, generative adversarial networks [13][14][15][16] to solve the problem of data imbalance and multiscale convolutional neural network to realize the fault diagnosis of the harmonic drive. Another approach is based on the use of deep learning [13][14][15][16] to extract, optimize, and classify different features. However, the main drawback of this approach is the need for large datasets of both healthy and faulty robot behavior, which are rarely available [17,18] due to an average mean time between failures of tens of thousands of hours [19]. To overcome this limitation, a database of simulated data in the presence of faults and failures can be built through a high-fidelity model of the robot arm. Such data, together with the experimental data, can be used to train DDMs for the proper detection and classification of the failure modes to which the machine is subjected. In [20,21], the authors proposed the application of a prognostic framework based on a long short-term memory neural network and the automatic recognition of the assembly task to detect mechanical faults in a robot for roller hemming employed in the automotive sector. Another example is represented by [20,21], in which physics-based simulation models, based on the digital twin concept, are used to obtain an estimate of the remaining useful life of machinery equipment.
With this aim, several studies focused on researching kinematic analyses and calibration of anthropomorphic robot manipulators, such as the UR5 collaborative robot from Universal Robots A/S, and on the creation of its digital twin to replicate the dynamic behavior, have obtained insights on internal quantities and studied the effects of specific degradations on both its position accuracy and overall performance. In Ref. [22], the authors present a kinematic model of a 6-DoF manipulator that utilizes the Denavit-Hartenberg convention to describe the relationships between the reference frames of each joint in the robot. The same approach was followed in [23] to obtain a simplified mass-concentrated model of a 3-DoF anthropomorphic manipulator: this kind of model can be used for kinematic calibration, as in [24]. In Ref. [25], a similar approach was used to perform forward and inverse kinematics to obtain the joint angular rotations from the desired Tool Center Point (TCP) pose and vice versa. To achieve this, a specific trajectory must be followed, and the inverse dynamic modeling approach allows calculating the torques required for such a trajectory. This requires considering the inertia of the links, Coriolis forces, gravity, and every exogenous force or torque applied on the robot, as shown in [26].
However, in high-fidelity models, particularly those involved in the creation of nonnominal behavior datasets for PHM purposes, a causal forward dynamic approach is preferred since it is dynamically consistent and satisfies the fundamental laws of physics at each point of the simulation.
The ability to replicate not only the actual joint positions, but also the torques necessary to perform the desired trajectory, is crucial for several studies, such as the ones aiming to minimize the power consumption of industrial manipulators [27,28] or the ones focused on using robot digital twins to evaluate the safety and ergonomics of specific applications [29,30]. The development of accurate models allows improvements in productivity and accuracy [31] as well. For these reasons, dynamic models of increasing complexity have been developed according to the specific objective and needs [32]. One of the main purposes related to the definition of dynamic models concerns the design and testing of the control algorithms. In Ref. [33], the PID controller of the robot was tuned based on the Ziegler-Nichols method, exploiting a classic dynamic model of a robot obtained by means of Lagrange equations [23]. In Ref. [34], the authors developed a simplified dynamic model by means of the Kane equation to study a robot design environment for the use of collaborative robots for rehabilitation purposes.
Good et al. [35] investigated the robot control algorithm and underlined the importance of taking into account the interactions between the robotic arm and the electromechanical drives. The model presented in [36] enhanced the mathematical representation of the robot structure, considering a simplified electric motor model. Tarn et al. [37] developed a thirdorder dynamic model of the robot obtained by combining the manipulator and the motor dynamics to study the effect of a nonlinear feedback control. Here, the model received as input the motor's armature voltage and returned as output the kinematic quantities of each joint.
An effective description of robot dynamics is necessary to achieve simulation models with high accuracy. For this reason, in both forward and inverse dynamics approaches, it is vital to know the correct dynamic parameters of the robot. On this topic, several studies concentrated on the development of identification techniques for these parameters [31], minimizing the error between the signal measured on the real robot, subjected to a proper excitation trajectory [38], and those obtained from the model simulation. Kovincic et al. [26] identified the dynamic parameters using a Moore-Penrose pseudo-inverse of the regressor of a UR5, obtaining a functional set of parameters to reproduce the robot's trajectory. However, this approach did not include any constraints on the values of the dynamic parameters, which can lead to a physically meaningless set of parameters. To solve this problem, in [39], a different method of identification of the dynamic parameters was proposed. Such an approach offers promising results, but requires the use of additional sensors such as accelerometers to measure the motion of each joint. An alternative way to identify the dynamic parameters involves the use of statistical techniques [36]. This method can provide excellent results in terms of overlapping of the model with the experimental results, but suffers from a very loose representation of the system physics. An alternative approach was developed by Raviola et al. [40], which applied a constrained genetic algorithm minimization, taking into account also the effect of mounting configuration for different temperature conditions. Hitherto, in the literature, the majority of dynamic models of robotic manipulators were developed with the goal of describing the macroscopic behavior of the entire system without taking into account the punctual description of each subsystem of the robot and their interaction. However, if a physics-based description of faults and degradation in the system has to be inserted for PHM analyses, a high level of detail and fidelity of the model is required, and the dynamic behavior and mutual interaction of each subsystem should be taken into account.
To solve these issues, the present paper describes a high-fidelity model of a UR5, a six-degrees-of-freedom (DoF) collaborative robot from Universal Robots A/S. Hayati's additional kinematic parameters were introduced to model assembly errors and obtain more physically representative kinematic calibration parameters. From a dynamic perspective, the dynamic parameters were identified for the UR5 under study according to [40]. Additionally, a detailed joint model was developed to accurately model the elements of each joint, including the electric motor and the reducer. For the latter, a strain wave gear dynamic model was considered together with the description of the internal interactions between its main elements, enabling a comparison between real robot currents and those generated by the model. The proper modeling and implementation of the reducer in the description of the robot dynamic behavior is of paramount importance in the robot design phases as well. The harmonic reducer is not the only element that can be exploited in the joint architecture of a robot. Another possible solution is the Zero-backlash High-Precision Roller Enveloping Reducer (ZHPRER) proposed in [41]. In their work, Wang et al. tested and verified the applicability of such a novel reducer by means of a dynamic model established using multibody dynamics. A further alternative is represented by the Rotate Vector (RV) reducer, whose applicability and performance were tested by means of a high-fidelity model as discussed in [42,43]. Finally, friction effects were included in each joint by means of the LuGre model, which allows a representation of both stiction and dynamic viscous friction depending on the operative conditions, temperature, and joint configuration [40]. The same logic can be applied to the development of digital twins of other manipulators with a different number of degrees of freedom. The major novelty of the present work regards the integration of the robot joint model in an electro-mechanical actuator architecture to obtain a physical-based description of the joint torque. This approach allows a better description of the dynamic behavior of the overall manipulator with respect to the previous literature, also obtaining a good approximation on the single joint's torque evaluation.
The present paper is organized as follows: first, the structure and characteristics of the UR5 robot under analysis are briefly recalled. Then, the mathematical model of the aforementioned robot is presented, and the strain wave gear mathematical representation is introduced, as well as that of the electric motor. Next, the performed experimental procedure is described and, to conclude, the comparison between the simulated data and experimental signals is presented in terms of position, speed, and torque for each joint. In particular, two quite different trajectories were tested in order to investigate the model ability of replicating the real robot behavior in different operational scenarios.

Mathematical Model of the UR5
The case study faced in this work concerns the development of a high-fidelity model of the UR5 robot of Universal Robots A/S depicted in in Figure 1. The UR5 is a collaborative anthropomorphic robot with six degrees of freedom, which corresponds to the same number of revolute joints. A peculiarity of the UR5 structure is that it does not have a spherical wrist; therefore, there is no point in which the three axes of the wrist joints intersect. From an implementation point of view, the first three joints are equal to each other, as are the last three. represented by the Rotate Vector (RV) reducer, whose applicability and performance were tested by means of a high-fidelity model as discussed in [42,43]. Finally, friction effects were included in each joint by means of the LuGre model, which allows a representation of both stiction and dynamic viscous friction depending on the operative conditions, temperature, and joint configuration [40]. The same logic can be applied to the development of digital twins of other manipulators with a different number of degrees of freedom. The major novelty of the present work regards the integration of the robot joint model in an electro-mechanical actuator architecture to obtain a physicalbased description of the joint torque. This approach allows a better description of the dynamic behavior of the overall manipulator with respect to the previous literature, also obtaining a good approximation on the single joint's torque evaluation.
The present paper is organized as follows: first, the structure and characteristics of the UR5 robot under analysis are briefly recalled. Then, the mathematical model of the aforementioned robot is presented, and the strain wave gear mathematical representation is introduced, as well as that of the electric motor. Next, the performed experimental procedure is described and, to conclude, the comparison between the simulated data and experimental signals is presented in terms of position, speed, and torque for each joint. In particular, two quite different trajectories were tested in order to investigate the model ability of replicating the real robot behavior in different operational scenarios.

Mathematical Model of the UR5
The case study faced in this work concerns the development of a high-fidelity model of the UR5 robot of Universal Robots A/S depicted in in Figure 1. The UR5 is a collaborative anthropomorphic robot with six degrees of freedom, which corresponds to the same number of revolute joints. A peculiarity of the UR5 structure is that it does not have a spherical wrist; therefore, there is no point in which the three axes of the wrist joints intersect. From an implementation point of view, the first three joints are equal to each other, as are the last three. The mathematical model of the UR5 collaborative robot, whose structure is reported in Figure 2, was developed in the MATLAB-Simulink environment. The mathematical model of the UR5 collaborative robot, whose structure is reported in Figure 2, was developed in the MATLAB-Simulink environment.
is provided as input to the six joints of the robot digital replica. The main elements of each joint (e.g., control logic, power electronics, electric motor, gearbox, and sensors) were modeled by finding a good compromise between simulation time and level of detail. To simulate one second, a computer with an Intel ® i7-8750H processor at 3.91 GHz equipped with 16 GB of RAM DDR4 takes about five minutes.
For each joint, by means of three nested control loops (joint position, motor velocity, and motor current), the reference voltage ( ) is calculated and provided as input to the power electronics block, where, through Pulse Width Modulation (PWM), the supply voltage ( ) is obtained. Inside the electric motor block, the motor current ( ) and velocity ( ) are computed and used as inputs to the current sensor and the optical encoder, respectively, to obtain the relative feedback signals ( , , ). The gearbox model compares the motor and the joint angular positions and velocities ( and , and ) to calculate the reaction torque to the motor/input shaft ( , ) and the joint output torque ( = , ). To this last value, the friction torque ( ) is then subtracted to obtain the useful torque ( ) acting on the joint shaft. The useful torques of each of the six joints of the manipulator are used in the robot arm block to solve the forward dynamics from which the joints' angular accelerations are obtained. Once integrated, they provide the joints' angular positions ( ), which are given as inputs to the magnetic encoder to determine the feedback signal ( ). Finally, through the forward kinematics of the robot arm, the pose (position and orientation) of the TCP is computed.

Control Logic and Power Electronics
Since the actual control algorithms implemented by the robot manufacturer on the real robot are unavailable, for the simulated manipulator a control logic based on proportional and integrative (PI) gains was adopted, similar to the ones proposed by [25,33]. The control parameters reported in Table 1 were tuned by analyzing the frequency response function derived for each joint as described in Appendix A.  As for the real manipulator, the set of the joints' angular positions (q set = [q set,1 , q set,2 , q set,3 , q set,4 , q set,5 , q set,6 ]) is provided as input to the six joints of the robot digital replica. The main elements of each joint (e.g., control logic, power electronics, electric motor, gearbox, and sensors) were modeled by finding a good compromise between simulation time and level of detail. To simulate one second, a computer with an Intel ® i7-8750H processor at 3.91 GHz equipped with 16 GB of RAM DDR4 takes about five minutes.
For each joint, by means of three nested control loops (joint position, motor velocity, and motor current), the reference voltage (V re f ) is calculated and provided as input to the power electronics block, where, through Pulse Width Modulation (PWM), the supply voltage (V sup ) is obtained. Inside the electric motor block, the motor current (i) and velocity ( q) to calculate the reaction torque to the motor/input shaft (T SWG,in ) and the joint output torque (T j = T SWG,out ). To this last value, the friction torque (T f ) is then subtracted to obtain the useful torque (T u ) acting on the joint shaft. The useful torques of each of the six joints of the manipulator are used in the robot arm block to solve the forward dynamics from which the joints' angular accelerations are obtained. Once integrated, they provide the joints' angular positions (q), which are given as inputs to the magnetic encoder to determine the feedback signal (q FB ). Finally, through the forward kinematics of the robot arm, the pose (position and orientation) of the TCP is computed.

Control Logic and Power Electronics
Since the actual control algorithms implemented by the robot manufacturer on the real robot are unavailable, for the simulated manipulator a control logic based on proportional and integrative (PI) gains was adopted, similar to the ones proposed by [25,33]. The control parameters reported in Table 1 were tuned by analyzing the frequency response function derived for each joint as described in Appendix A.
The joint control logic was structured into three nested loops. The outer one compares the desired (q set ) and the feedback (q FB ) signals of the joints' angular positions to determine the position error (e pos ) of each joint. Such value is then multiplied by a proportional gain (P), to which the desired motor angular velocity ( . q M,target ) is summed to obtain the reference motor angular velocity ( . q M,re f ). This approach adds a speed feedforward branch to enhance the control capabilities, exploiting the a priori knowledge of the entire trajectory, usually known in classic robot operations. The mismatch (e vel ) between such signal ( . q M,re f ) and the feedback value ( . q M,FB ) coming from the optical encoder determines the velocity error which, once compensated, becomes the reference current signal. The motor current control loop evaluates the reference voltage (V re f ) as a function of the mismatch (e cur ) between the reference (i re f ) and the feedback (i FB ) motor currents. This is used by the power electronics block, which regulates the network voltage (V net = 48 V) providing the PWM supply voltage (V sup ) to the electric motors.

Electric Motor
The UR5 is equipped with six Permanent Magnets AC Synchronous Motors (PMSMs), which combine the high precision in speed control typical of a synchronous motor with the simple and robust design of an asynchronous one. To reduce the model complexity, the motors of the simulated manipulator were described through their equivalent DC monophase models receiving as input the supply voltage (V sup ) and providing, as output, the motor torque T M = k t i, where k t is the torque constant and i is the motor current. The motor dynamics are described through two electrical and mechanical equilibrium equations: where R and L are, respectively, the equivalent DC motor winding resistance and inductance, k e is the motor back-EMF constant, assumed, at first approximation, to be equal to k t , . q M is the motor angular velocity, J M is the rotor inertia, and T SWG,in is the reaction torque generated by the reducer on the motor shaft and represents the disturbance torque to the electric motor. The values of the aforementioned parameters are provided in Table 2.

Gearbox
The Strain Wave Gear (SWG), also known as harmonic drive, is a reducer widely adopted in several industrial fields, ranging from industrial and collaborative robotics [44,45] to aerospace applications [46], due to its high reduction ratio, virtually zero-backlash, and compact and light-weight design. This gear can be functionally described through its three main components: Wave Generator (WG), Flexspline (FS), and Circular Spline (CS). The WG consists of an elliptical cam (wave generator plug) inserted in a deformable ball bearing which is then pushed inside the FS, a thin-walled gear with external toothing engaging with the CS, a rigid gear with internal toothing. The behavior of the gear changes according to the mounting configuration of these three components. As an example, for the UR5 manipulator, the input shaft is always integral to the WG, while the output one is integral to the CS, as the FS fixed to the joint case. This leads to a reduction ratio (τ), defined as: where Z FS and Z CS are, respectively, the number of the flexspline and the circular spline teeth. Since, for the robot under analysis, Z FS = 200 and Z CS = 202, the gear ratio reported in Equation (2) is τ = 101. Due to its elliptical shape, the rotation of the wave generator elastically deforms the flexspline along its major axis, allowing it to mesh with the circular spline. Since the number of FS teeth is two fewer than the one of CS, for a full rotation of the input shaft, the output one rotates of an angle corresponding to two teeth. For the proposed case study, according to Equation (2), for each 360 • of rotation of the WG, the CS rotates, in the same direction, by about 3.56 • .
To replicate the behavior of the strain wave gear, a translation equivalent model, based on the research in [47,48], was implemented on each joint of the UR5 manipulator. This allows describing the cyclic deformation of the flexspline and its interaction with the circular spline through equivalent tangential and radial displacements. Damping-stiffness (c − k) coefficients [49] were used to describe the interactions not only between the contact interfaces within the SWG, but also between the motor shaft and the WG, and between the CS and the joint output flange. Therefore, it is possible to decouple the dynamics of the SWG from those of the input and the output shafts, leading to a four-DoF model. The schematic representations of the proposed gear model is reported in Figure 3. Here, the curved continuous arrows are used to schematize the rotation of the gear elements, while the dashed ones represent their small oscillations due to their compliant connections to the joint housing represented through a damper and a torsional stiffness.
Robotics 2023, 11, x FOR PEER REVIEW 7 of 23 example, for the UR5 manipulator, the input shaft is always integral to the WG, while the output one is integral to the CS, as the FS fixed to the joint case. This leads to a reduction ratio ( ), defined as: where and are, respectively, the number of the flexspline and the circular spline teeth. Since, for the robot under analysis, = 200 and = 202, the gear ratio reported in Equation (2) is = 101. Due to its elliptical shape, the rotation of the wave generator elastically deforms the flexspline along its major axis, allowing it to mesh with the circular spline. Since the number of FS teeth is two fewer than the one of CS, for a full rotation of the input shaft, the output one rotates of an angle corresponding to two teeth. For the proposed case study, according to Equation (2), for each 360° of rotation of the WG, the CS rotates, in the same direction, by about 3.56°.
To replicate the behavior of the strain wave gear, a translation equivalent model, based on the research in [47,48], was implemented on each joint of the UR5 manipulator. This allows describing the cyclic deformation of the flexspline and its interaction with the circular spline through equivalent tangential and radial displacements. Damping-stiffness ( − ) coefficients [49] were used to describe the interactions not only between the contact interfaces within the SWG, but also between the motor shaft and the WG, and between the CS and the joint output flange. Therefore, it is possible to decouple the dynamics of the SWG from those of the input and the output shafts, leading to a four-DoF model. The schematic representations of the proposed gear model is reported in Figure 3. Here, the curved continuous arrows are used to schematize the rotation of the gear elements, while the dashed ones represent their small oscillations due to their compliant connections to the joint housing represented through a damper and a torsional stiffness. To properly describe the behavior of the gearbox, both the kinematic and dynamic parameters (e.g., inertia and torsional stiffness) of the gears were characterized. According to [50], the UR5 collaborative robot was equipped with the same joint/reducer size for the first three joints (base, shoulder, and elbow), while a smaller one is mounted on the last three (wrist 1, wrist 2, and wrist 3). For convenience, in the present document, these two sizes will be labeled as A and B, respectively, for the bigger and the smaller reducers.
The translation model of the SWG reported in Figure 3b is based on the possibility to correlate the angular quantities of the actual gear (e.g., input and output shaft rotations) to the equivalent translational ones through a primitive equivalent radius ( ) equal to the undeformed FS radius. To respect the actual transmission ratio, the elliptical cam To properly describe the behavior of the gearbox, both the kinematic and dynamic parameters (e.g., inertia and torsional stiffness) of the gears were characterized. According to [50], the UR5 collaborative robot was equipped with the same joint/reducer size for the first three joints (base, shoulder, and elbow), while a smaller one is mounted on the last three (wrist 1, wrist 2, and wrist 3). For convenience, in the present document, these two sizes will be labeled as A and B, respectively, for the bigger and the smaller reducers.
The translation model of the SWG reported in Figure 3b is based on the possibility to correlate the angular quantities of the actual gear (e.g., input and output shaft rotations) to the equivalent translational ones through a primitive equivalent radius (r g ) equal to the undeformed FS radius. To respect the actual transmission ratio, the elliptical cam equivalent angle (α n ) was described as a function of the gear ratio (τ) and the tooth pressure angle (α t ) as:

of 22
While the inertia values of the wave generator (J WG ), the flexspline (J FS ), and the circular spline (J CS ) for the two reducers were derived from the CAD files, the torsional stiffness (K WG , K FS , K CS ) of the single elements of the strain wave gears were approximated using the formula for hollow cylinders: where E = 206 × 10 9 N/m 2 and ν = 0.26 are the steel Young modulus and Poisson coefficient, l is the shaft length, and d ext and d int are, respectively, the external and internal diameters of the single component. On the contrary, the linear stiffnesses used to model the contact between the WG outer race and the FS inner surface (K b ) and the one associated with the FS-CS teeth meshing (K m ) were derived from the literature [51], while all the damping factors were set equal to 0.1. The values used in the proposed model are summarized in Table 3. To calculate the torque transmitted by the reducer (T j = T SWG,out equal to T FS or T CS , according to the gear mounting configuration), it is necessary to solve the dynamic equilibrium of the WG, the FS, and the CS schematized in Figure 3b. The equations describing the interactions among the single elements of the gear are [51]: where, for better clarity, the general notation c x = cos(x) and s x = sin(x) was adopted. The nine equations used to describe the kinematic and dynamic behavior of the single reducer can be classified into two groups: • Group 1: comprised of the first four equations describing, respectively, the dynamic equilibrium of the wave generator, the flexspline radial and tangential contributions, and the one of the circular spline; • Group 2: comprised of the last five equations, representing the interactions between the input shaft and the elliptical cam, through the WG bearing, between the FS and the CS teeth, and between the FS and the joint case and the CS with the output shaft.
To ensure the correct functioning of the proposed gear models, the initial conditions were defined according to the starting angular positions of the input/motor shafts and the output/joint shafts, obtained from the solution of the following set of equations:

Friction
In contrast with [48], friction was not implemented at the single contact interfaces of the SWG, but defined as an overall joint friction torque. This was chosen to implement the Coulomb ( f c ) and the viscous ( f v ) friction coefficients derived through the base dynamic parameters identification process in the proposed robot model, as described in [52]. However, the friction model used to build the regression matrix in [40] does not properly replicate the actual trend of a joint friction. To describe the behavior associated with the input and output bearings inside the robot joint, the Stribeck curve proposed in [53] was modified in the one reported in Equation (7) through the exponent associated with the joint angular velocity ( . q) [51]. The friction torque T f between the input and output of each joint was modeled through the interaction of the LuGre elastic bristles as proposed in [54], leading to: where σ 0 = 5 × 10 6 Nm/rad and σ 1 = 100 Nms/rad are, respectively, the stiffness and the damping coefficient of the bristles subjected to an average deflection z at a speed . z, .
q is the joint angular velocity (equal to the relative velocity of the two contact surfaces), g . q also depends on surface materials, lubrication regime, and temperature, . q s = 0.5 rad/s is the sliding speed threshold used to describe the Stribeck effect, K static = 0.3 is the static friction increment with respect to the Coulomb friction f c which leads to the stiction friction force, while k q and f v are the Stribeck parameter and the viscous friction coefficient, respectively, selected in order to reproduce the experimental behavior [52]. Other possible formulations of the friction curve for a robot joint can be found in [55,56]. The LuGre model allows considering micro-displacements, which are represented with a spring-like behavior.

Sensors
According to the scheme reported in Figure 2, each joint of the robot was equipped with three sensors:

1.
Optical encoder: an encoder integral to the motor/input shaft used to derive the motor angular velocity ( Magnetic encoder: an encoder integral to the joint/output shaft used to measure the joint actual angular position (q) used to close the position control loop through the feedback signal (q FB ); 3.
Current sensor: used to measure the motor current (i) necessary to close the current control loop through the feedback signal (i FB ).
Such sensors were modeled using first-order transfer functions for the current sensors, while second-order ones were implemented for both input and output encoders. For a more realistic description of the real sensors, noise was added to the feedback signals, which were then digitalized by means of a sample and hold circuit and a quantizer. The noise standard deviation was the same as the one measured from the UR5.
Moreover, in the real robot, each joint is also equipped with a thermocouple used to measure the joint temperature. However, since heat transfer was neglected in the proposed joint model, this sensor was not implemented.

Six-DoF Forward Kinematics and Dynamics
The robot arm subsystem, shown in Figure 2, contains the Simscape Multibody model with which the single joints and links of the UR5 digital replica are connected, as depicted in Figure 4. joint actual angular position ( ) used to close the position control loop through the feedback signal ( ); 3. Current sensor: used to measure the motor current ( ) necessary to close the current control loop through the feedback signal ( ).
Such sensors were modeled using first-order transfer functions for the current sensors, while second-order ones were implemented for both input and output encoders. For a more realistic description of the real sensors, noise was added to the feedback signals, which were then digitalized by means of a sample and hold circuit and a quantizer. The noise standard deviation was the same as the one measured from the UR5.
Moreover, in the real robot, each joint is also equipped with a thermocouple used to measure the joint temperature. However, since heat transfer was neglected in the proposed joint model, this sensor was not implemented.

Six-DoF Forward Kinematics and Dynamics
The robot arm subsystem, shown in Figure 2, contains the Simscape Multibody model with which the single joints and links of the UR5 digital replica are connected, as depicted in Figure 4.
As in Section 2.3, for better clarity, the general notation = cos( ) and = sin( ) was adopted. The values of the variables in Equation (8) were obtained through post-processing of the ones derived from the dual robot calibration performed by the manufacturer [57]. Nominally, each UR5 collaborative robot is geometrically described by the Denavit-Hartenberg (DH) parameters provided in [58]. However, due to production tolerances and mounting errors, it is necessary to calibrate each robot to achieve the desired position The relative position and orientation of the joint reference frame i with respect to the previous one i − 1 are described through the 4 × 4 roto-translation matrix i−1 A i as: As in Section 2.3, for better clarity, the general notation c x = cos(x) and s x = sin(x) was adopted. The values of the variables in Equation (8) were obtained through postprocessing of the ones derived from the dual robot calibration performed by the manufacturer [57]. Nominally, each UR5 collaborative robot is geometrically described by the Denavit-Hartenberg (DH) parameters provided in [58]. However, due to production tolerances and mounting errors, it is necessary to calibrate each robot to achieve the desired position accuracy and repeatability. According to [59], through the DH convention, for links between two nearly parallel joints, it is not possible to obtain, after the calibration process, corrections to the robot geometric parameters with a physical meaning. In the proposed use case, for example, the values of the DH parameter d for links 2 and 3 have an order of magnitude of kilometers. To avoid using unrealistic numbers, it was necessary to merge the Denavit-Hartemberg convention with the one proposed by Hayati in [59], leading to the roto-translation matrix shown in Equation (8) and to the kinematic parameters reported in Table 4.
To complete the description of the robot arm, the dynamic properties reported in Table 5 were associated with the single links of the manipulator. From a dynamic point of view, the description of the behavior of the manipulator was obtained by exploiting the Simulink/Simscape environment. Within it, a casual object-oriented programming language, individual system model, and its interconnections describe the underlying dynamics equation. This approach allows the design of the model to be based on the actual structure of the physical system. The kinematic and dynamic parameters reported in Tables 4 and 5 were used to build the multibody Simscape model schematized in Figure 5. As shown in Figure 2, each joint sub-model receives as input the torque ( ) provided by the reducer and calculated through Equation (5). As an output, each block provides the joint angular position ( ). Within the sub model of each joint, mass, Coriolis, and gravity matrices are automatically constructed, while the friction matrix is provided by exploiting the LuGre model as illustrated in Section 2.4, to compute the direct dynamics of the manipulator.

Experimental Testing
To examine the goodness of the proposed model, several trajectories were applied to both the UR5 and its digital replica. As an example, in the present paper, two different test trajectories were analyzed. On one hand, specific movements designed to properly excite the robot dynamic parameters, such as the ones reported in [26,61,62], were provided as inputs to evaluate the ability of the proposed model to replicate a high-dynamic As shown in Figure 2, each joint j sub-model receives as input the torque (T j ) provided by the reducer and calculated through Equation (5). As an output, each block provides the joint angular position (q j ). Within the sub model of each joint, mass, Coriolis, and gravity matrices are automatically constructed, while the friction matrix is provided by exploiting the LuGre model as illustrated in Section 2.4, to compute the direct dynamics of the manipulator.

Experimental Testing
To examine the goodness of the proposed model, several trajectories were applied to both the UR5 and its digital replica. As an example, in the present paper, two different test trajectories were analyzed. On one hand, specific movements designed to properly excite the robot dynamic parameters, such as the ones reported in [26,61,62], were provided as inputs to evaluate the ability of the proposed model to replicate a high-dynamic behavior of the real robot. On the other hand, since, according to [63], pick-and-place is the most common application for both industrial and collaborative robots, the model performance was also evaluated with a trajectory representative of such movements. In both cases, the desired joints' angular positions of the test trajectories, reported in Figure 6, were applied to the UR5 through port 29999 [64], while the joints' angular positions and velocities and the motors' current signals were acquired at a frequency of 125 Hz using port 30003.
vided by the reducer and calculated through Equation (5). As an output, each block provides the joint angular position ( ). Within the sub model of each joint, mass, Coriolis, and gravity matrices are automatically constructed, while the friction matrix is provided by exploiting the LuGre model as illustrated in Section 2.4, to compute the direct dynamics of the manipulator.

Experimental Testing
To examine the goodness of the proposed model, several trajectories were applied to both the UR5 and its digital replica. As an example, in the present paper, two different test trajectories were analyzed. On one hand, specific movements designed to properly excite the robot dynamic parameters, such as the ones reported in [26,61,62], were provided as inputs to evaluate the ability of the proposed model to replicate a high-dynamic behavior of the real robot. On the other hand, since, according to [63], pick-and-place is the most common application for both industrial and collaborative robots, the model performance was also evaluated with a trajectory representative of such movements. In both cases, the desired joints' angular positions of the test trajectories, reported in Figure 6, were applied to the UR5 through port 29999 [64], while the joints' angular positions and velocities and the motors' current signals were acquired at a frequency of 125 Hz using port 30003.  For both trajectories, to clearly plot the joints' torques and to improve the paper readability, the current signals were filtered with a low-pass filter at a cutoff frequency of 10 Hz. On the contrary, no filtering was required for the joints' angular positions.

Dynamic Parameters Excitation Trajectory
Due to the complexity of such excitation trajectory, to move the real manipulator, the robot built-in function servoj (target joint configuration, acceleration, velocity, time, lookahead time, gain) was adopted [65]. The numerical values of acceleration, velocity, time, lookahead time, and gain were set equal to 0 rad/s 2 , 0 rad/s, 0.08 ms, 0.03, and 500, respectively.
The ability of the proposed model to follow the desired trajectory can be seen in Figure 7, where both set and feedback joints' angular position signals are reported.
Since, for all the joints, the trends of the three curves well overlap, Figure 8 reports the mismatches between such set and feedback signals to better highlight the tracking error of both the real and the simulated robots. For better comparison, the axes scales are the same for all the joint plots.
By analyzing the position errors for each joint, it can be seen that they not only have the same order of magnitude for both manipulators, but they also share similar trends, proving that the proposed model is able to replicate the actual behavior of the real UR5. It should be also noticed that the position error of the proposed model is lower than that of the real robot, with a maximum value of about 0.08 • . However, the apparent better performance of the proposed model is likely related to the simplifications adopted. Among all, the strain wave gear torsional hysteresis [66] was not modeled in the translation equivalent model of the strain wave gear of Figure 3. Such a simplification leads to the gears always having the same torsional stiffness during the entire trajectory, while studies in the literature [67,68] showed that this is a function of the applied torque. This, together with the approximations of the gear stiffness values, might lead the digital replica to be stiffer than the real robot, allowing its control logic to better compensate for position errors.
For both trajectories, to clearly plot the joints' torques and to improve the paper readability, the current signals were filtered with a low-pass filter at a cutoff frequency of 10 Hz. On the contrary, no filtering was required for the joints' angular positions.

Dynamic Parameters Excitation Trajectory
Due to the complexity of such excitation trajectory, to move the real manipulator, the robot built-in function servoj (target joint configuration, acceleration, velocity, time, lookahead time, gain) was adopted [65]. The numerical values of acceleration, velocity, time, lookahead time, and gain were set equal to 0 rad/s 2 , 0 rad/s, 0.08 ms, 0.03, and 500, respectively.
The ability of the proposed model to follow the desired trajectory can be seen in Figure 7, where both set and feedback joints' angular position signals are reported. Since, for all the joints, the trends of the three curves well overlap, Figure 8 reports the mismatches between such set and feedback signals to better highlight the tracking error of both the real and the simulated robots. For better comparison, the axes scales are the same for all the joint plots. robot built-in function servoj (target joint configuration, acceleration, velocity, time, lookahead time, gain) was adopted [65]. The numerical values of acceleration, velocity, time, lookahead time, and gain were set equal to 0 rad/s 2 , 0 rad/s, 0.08 ms, 0.03, and 500, respectively.
The ability of the proposed model to follow the desired trajectory can be seen in Figure 7, where both set and feedback joints' angular position signals are reported. Since, for all the joints, the trends of the three curves well overlap, Figure 8 reports the mismatches between such set and feedback signals to better highlight the tracking error of both the real and the simulated robots. For better comparison, the axes scales are the same for all the joint plots. Nevertheless, the low tracking errors between the set and feedback joints' angular positions translate into an actual path of the TCP close to the ideal/desired one, which is of paramount importance for applications such as welding, dispensing, and gluing. The position and the orientation of the robot TCP with respect to the fixed base reference frame ( 0 A 6 ), in fact, are described by the combination of rigid body transformation matrices described in Equation (8) for each joint as: According to the control logic reported in Figure 2, the position error was used to compute the target/reference motor angular velocity whose difference with the feedback signal, once compensated, provides the motor current for each joint. By multiplying the feedback currents (i FB ) by the reduction ratio (τ) and the motor torque constants (k t ) reported in Table 2, the joint reducer output torques of the two manipulators are derived and compared in Figure 9.
of paramount importance for applications such as welding, dispensing, and gluing. The position and the orientation of the robot TCP with respect to the fixed base reference frame ( 0 A6), in fact, are described by the combination of rigid body transformation matrices described in Equation (8) for each joint as: According to the control logic reported in Figure 2, the position error was used to compute the target/reference motor angular velocity whose difference with the feedback signal, once compensated, provides the motor current for each joint. By multiplying the feedback currents ( ) by the reduction ratio ( ) and the motor torque constants ( ) reported in Table 2, the joint reducer output torques of the two manipulators are derived and compared in Figure 9.  Even though the simulated torques (Model) well approximate the measured ones (UR5), some intervals of the tested trajectory show higher mismatches between the two sets of data. Such a residual difference might be related to several causes, such as: • Implementation of simplified DC models of the joint motors instead of using threephase ones, whose parameters, such as resistance and inductance, are not fully known. Moreover, the motor back-EMF constant (k e ) was assumed to be constant for the entire trajectory, while it slightly changes with the applied load. • A possible incorrect description of the joint friction. As already mentioned, the identification algorithms provided by [40,52] only take into account viscous and Coulomb friction, while they are not able to describe more complex trends, such as the ones affecting the real robot. To do so, friction was modeled through Equation (7), whose coefficients were tuned according to several experimental campaigns aiming to reach a good overall fit between simulated and measured motor currents. • Lack of both the kinematic error and the gear torsional hysteresis in the translationequivalent models of the strain wave gears.

•
The joint control logic schematized in Figure 2 is not the same as the one implemented in the real robot, which is unknown. As an example, studies such as [69,70] highlight how manipulators from Universal Robots are equipped with vibration suppression algorithms, which were not implemented in the current model. Having two different control logics deeply affects the torque trends. Since, according to the scheme reported in Figure 2, the current control loop is the most internal one, the reference value of the motor current is affected by both the position and the velocity control loops. As a direct consequence, since, as depicted in Figure 8, the joint position error is different between the two manipulators, the velocity reference signal will also be different.
Hence, based on the adopted PI gains, the motor angular velocity error differs from the one of the real robot as well.
• While the control parameters reported in Table 1 are constant, it should not be excluded that Universal Robots A/S might adopt variable settings according to the specific operating conditions of the manipulator. such as very fast or very slow movements.
All the aforementioned major simplifications adopted in the present work increase the difficulty of building a digital shadow of the UR5 able to replicate the actual joint torques. More accurate estimates could be reached with a more detailed model of the single components of the robot arm, but this would negatively affect the computational cost of the robot digital replica, leading to simulation times high enough to nullify the advantage of higher precision in describing the actual dynamic behavior of the manipulator. Nevertheless, the obtained results on the torque signals are quite representative of the real behavior and constitute an innovation with respect to the usually employed literature models, in which a much lighter description of the internal quantities of the robot is present.

Pick-and-Place Trajectory
In contrast with the previous trajectory, the pick-and-place trajectory was obtained using the movel and movej commands [65] with different settings of speed and acceleration using the Polyscope programming interface on the robot teach pendant. An analysis, similar to the one reported in Section 3.1 for the dynamic parameters identification trajectory, was also conducted for the pick-and-place trajectory. The results are reported in Figures 10-12.
By comparing Figure 11 with Figure 8, it can be seen that, on average, the joints' angular positions errors are lower for the pick-and-place trajectory. This is most likely related to the lower dynamics involved in this task with respect to the rapid and continuously changing movements required to identify the UR5 base dynamic parameters.
The lower angular velocities involved during the pick-and-place give more time to the joint control logic to compensate for a mismatch between the desired and actual joint positions. This is particularly true for the last joint, which rotates a few tenths of degrees just to maintain the desired orientation of the TCP.
However, a position error of approximately 0.24 • , about 0.05 • higher than the one of the UR5, is visible at about 8.5 s in the second wrist. This is caused by a rapid deceleration of the joint from a speed of 180 • /s to 0 in a short interval. In so doing, the control logic is required to rapidly adjust its target from keeping the joint at its highest rotation speed to break it according to a trapezoidal velocity trend. This rapid change generates a current spike, visible in the torque peak in Figure 12. The lower angular velocities involved during the pick-and-place give more time to the joint control logic to compensate for a mismatch between the desired and actual joint positions. This is particularly true for the last joint, which rotates a few tenths of degrees just to maintain the desired orientation of the TCP.    However, a position error of approximately 0.24°, about 0.05° higher than the one of the UR5, is visible at about 8.5 s in the second wrist. This is caused by a rapid deceleration of the joint from a speed of 180°/s to 0 in a short interval. In so doing, the control logic is required to rapidly adjust its target from keeping the joint at its highest rotation speed to break it according to a trapezoidal velocity trend. This rapid change generates a current spike, visible in the torque peak in Figure 12. For the tested trajectory, the shoulder joint was affected by the highest errors. However, it should be noted that it was also the most loaded joint and its dynamics are directly affected by those of the other four following joints. It should be pointed out that, while, as shown in Figure 9, there is a good match between the measured and the simulated torques of the last joint (wrist 3), the same is not demonstrated for the pick-and-place trajectory in Figure 12. The difference in the two trends is probably caused by the possible different control logic adopted by Universal Robots A/S to handle very low movements, which would lead to extremely low angular velocity errors, as in the proposed case study where the maximum error between reference and feedback speed for all the joints is about 0.2 °/s over the entire trajectory. In this scenario, the joint is almost still, but since the velocity error, even if low, is still present, it generates a current signal which is proportional to such error.

Conclusions
The present paper describes a high-fidelity dynamic model of a UR5 collaborative robot where the main components of robot joints were described in a detailed and physically representative way. Testing trajectories of different types were applied to both the For the tested trajectory, the shoulder joint was affected by the highest errors. However, it should be noted that it was also the most loaded joint and its dynamics are directly affected by those of the other four following joints. It should be pointed out that, while, as shown in Figure 9, there is a good match between the measured and the simulated torques of the last joint (wrist 3), the same is not demonstrated for the pick-and-place trajectory in Figure 12. The difference in the two trends is probably caused by the possible different control logic adopted by Universal Robots A/S to handle very low movements, which would lead to extremely low angular velocity errors, as in the proposed case study where the maximum error between reference and feedback speed for all the joints is about 0.2 • /s over the entire trajectory. In this scenario, the joint is almost still, but since the velocity error, even if low, is still present, it generates a current signal which is proportional to such error.

Conclusions
The present paper describes a high-fidelity dynamic model of a UR5 collaborative robot where the main components of robot joints were described in a detailed and physically representative way. Testing trajectories of different types were applied to both the real robot and its digital replica, showing a good match between experimental and simulated signals in terms of position, speed, and joint torques. In particular, the latter information represents the innovation with respect to a dynamic literature model, which usually disregards the dynamics of the joint electric drives.
Future work will be focused on simulating the presence of faults and failures of different origins and magnitudes to evaluate their impact on the overall behavior of the manipulator. The digital replica of the real robot will be used as a virtual test bench to generate a dataset of nominal and degraded responses, paving the way for a more detailed and thorough model-based PHM investigation. In addition, the proposed model of the strain wave gear will be improved to take into account additional dynamic effects, such as the kinematic error, the load-dependent torsional stiffness, and its relative degradation in time, due to an in-depth analysis of the internal interactions among the gear's constitutive elements.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
The following symbols are used in the manuscript:

Appendix A
The control parameters reported in Table 1 were obtained by analyzing the joints Frequency Response Function (FRF) to ensure the stable dynamic behavior of the robot. To do so, each control loop (in the order: current, velocity, and position) was tested with sinusoidal signals of different frequencies to test the ability of the joint to follow the reference signals under different operating conditions. In detail, the frequency ranges were 10-10,000 Hz, 1-200 Hz, and 0.1-10 Hz, respectively, for the current, velocity, and position control loops. The equivalent joint model used for this study is similar to the one reported in Figure 2, and it is schematized in Figure A1. The proportional (P) and integrative (I) parameters of each control loop were tuned from the inner (current) to the outer (position) loop using the model in Figure A1 configured as: The required stability thresholds for the proposed model were a phase and a gain margin of at least 60° and 8 dB.
The joint load ( ), represented by the "robot equivalent inertia" block in Figure A1, was chosen as the mean value of the diagonal elements of the positive definite symmetric inertia matrix ( ), as a function of the robot joints' angular positions [71].  The proportional (P) and integrative (I) parameters of each control loop were tuned from the inner (current) to the outer (position) loop using the model in Figure A1 configured as: The required stability thresholds for the proposed model were a phase and a gain margin of at least 60 • and 8 dB.
The joint load (J eq ), represented by the "robot equivalent inertia" block in Figure A1, was chosen as the mean value of the diagonal elements M ii of the positive definite symmetric inertia matrix M(q), as a function of the robot joints' angular positions q [71]. Figure A2 shows how the robot equivalent inertia acting on the base (M 11 ), shoulder (M 22 ), elbow (M 33 ), and first wrist (M 44 ) joints varies as a function of the position of the two consecutive joints. margin of at least 60° and 8 dB.
The joint load ( ), represented by the "robot equivalent inertia" block in Figure A1, was chosen as the mean value of the diagonal elements of the positive definite symmetric inertia matrix ( ), as a function of the robot joints' angular positions [71]. The values for the last two wrists of the UR5 are not reported because, since the robot wrist group is quite compact, the angular positions of the second wrist almost do not affect the equivalent inertia to which the first one is subjected. The same reasoning can be applied to wrist 3, which moves link 6 and which can be approximated with a cylinder. This is clearly visible from Table A1, where the maximum, mean, and minimum values of the robot equivalent inertia are listed. For the wrist group, the equivalent inertia is almost constant for the fourth joint of the robot, while no sensitive variations are present for the last two joints.  The values for the last two wrists of the UR5 are not reported because, since the robot wrist group is quite compact, the angular positions of the second wrist almost do not affect the equivalent inertia to which the first one is subjected. The same reasoning can be applied to wrist 3, which moves link 6 and which can be approximated with a cylinder. This is clearly visible from Table A1, where the maximum, mean, and minimum values of the robot equivalent inertia are listed. For the wrist group, the equivalent inertia M ii is almost constant for the fourth joint of the robot, while no sensitive variations are present for the last two joints.