1. Introduction
Understanding a manipulator dynamic model is crucial in robotics applications, including developing model-based control systems [Reference Jiang, Huang, Bi, Ma and Yu1], building Digital Twin applications [Reference Lu, Liu, Wang, Huang and Xu2], or operating the robot safely within its kinematic and dynamic limits [Reference Faroni, Beschi, Guarino Lo Bianco and Visioli3]. To reduce the energy consumption of robot-assisted manufacturing systems [Reference Duflou, Sutherland, Dornfeld, Herrmann, Jeswiet, Kara, Hauschild and Kellens4], a dynamic model is essential [Reference Carabin, Wehrle and Vidoni5], as it enables the application of many techniques that can significantly improve the energy efficiency of an industrial robot [Reference Ruiz, Santos, Croes, Desmet and da Silva6, Reference Liu, Huang, Liu, Guo, Wang and Tan7]. In addition to the dynamic model, which relates generalized forces to motion quantities, the electrical model of the robot must be investigated if a precise energy consumption prediction is required.
Motion planning significantly impacts energy consumption during task execution [Reference Garcia, Bittencourt and Villani8–Reference Vidussi, Boscariol, Scalera and Gasparetto10]. It is possible to optimize trajectories by knowing how losses depend on motion parameters [Reference Liu, Wang, Wang and Wang11, Reference Carabin and Scalera12], or from a detailed knowledge of the dynamic model, on which to search for a minimum-energy solution for a given trajectory [Reference Paes, Dewulf, Elst, Kellens and Slaets13, Reference Boscariol and Richiedei14].
Another option is to use a data-driven approach, which involves processing trajectory and consumption data using machine learning, as in ref. [Reference Zhang and Yan15]. In that work, a set of basic motion tasks are used to build a model that can predict the overall energy consumption on the basis of two parameters that describe the trapezoidal speed profile. While that approach has been proved to be accurate, such accuracy is guaranteed only within the spectrum of the operative conditions used to train the machine learning software. In other words, the prediction model needs to be rebuilt from scratch if other motion primitives or different operative conditions need to be taken into account. Developing a true dynamic model is certainly more challenging as proper measurements and a solid understanding of robot dynamics are needed; nonetheless, the outcome is a more general tool with the broadest range of possibilities.
The study of a robot dynamic model is a topic that has already been discussed extensively in the literature. Earliest examples of robot identifications techniques can be found in the works by Gautier et al., who provide a method for the determination of the minimum set of inertial parameters of serial robots [Reference Gautier and Khalil16], as well as a way to define exciting trajectories for robot identification [Reference Presse and Gautier17]. Atkeson et al. proposed the estimation of inertial parameters of manipulators loads and links [Reference Atkeson, An and Hollerbach18]. Furthermore, Grotjahn et al. used a two-step approach for the friction and rigid body identification of robot dynamics [Reference Grotjahn, Daemi and Heimann19].
The approaches for the identification are diverse [Reference Wu, Wang and You20]: their need is motivated by the fact that rarely a robot manufacturer discloses sufficient information for building an accurate dynamic model from nameplate data. Focusing on offline identification methods, the one that are more frequently employed for industrial robots are either based on the identification of the parameters by a linear regression procedure [Reference Swevers, Ganseman, Tukel, de Schutter and Van Brussel21, Reference Neubauer, Gattringer and Bremer22] or by an optimization routine [Reference Jin and Gans23–Reference Ulici, Codrean and Natsakis26]. In both cases, the target is to ‘fit’ a model of known structure to a set of measured data to reduce the error between the latter and the output of the model being tuned. Both methods have been tested on a similar robot to the one used in this work, that is, the UR5 robot by Universal Robots: the results found by means of a linear regression method in ref. [Reference Tika, Ulmen and Bajcinca27] are indeed as accurate as the ones found in ref. [Reference Kovincic, Müller, Gattringer, Weyrer, Schlotzhauer and Brandstötter28]. Furthermore, in ref. [Reference Gaz, Cognetti, Oliva, Giordano and De Luca29] the dynamic identification of a robot similar in size and payload to the UR5e, that is, the Franka Emika Panda arm, is presented and a set of feasible link parameters is retrieved for more realistic dynamic simulations or when implementing in real-time robot control laws using recursive Newton-Euler algorithms.
An alternative approach involves three iterative loops for estimating mass, inertia, and friction, respectively, according to a linear model, then to improve the friction model by including nonlinearities, and finally the third for studying the Stribeck effect [Reference Deng, Shang, Zhang, Zhen and Cong30]. The identified robot is similar, in its architecture, to an UR manipulator. A method based on deep-learning is also presented for the UR5 model in ref. [Reference Wang, Shao, Yang and Liu31]. Another option, which has been suggested in the work [Reference Gautier and Poignet32], is to use the framework of the extended Kalman filter together with a weighted least squares estimation: as recognized by authors of that work, however, the method of [Reference Gautier and Poignet32] is more suited to an online identification, rather than an offline one. Widening the perspective, other less specific to robotics solutions could be potentially used as well, such as the neurofuzzy approach to empirical data modeling for nonlinear dynamical processes described in detail in the book [Reference Harris, Hong and Gan33].
Special attention is often paid to friction, as it is often the most challenging phenomenon to be identified. Friction appears to be related to several operating conditions [Reference Zhang, Wang, Chen, Chen, Lin and Xu34], and in lightweight collaborative robots is mainly due to the losses in the strain wave gearboxes. The friction model can be arbitrarily simplified, as much as describing it only by means of the viscous and Coulomb terms without any major modeling error [Reference Kovincic, Müller, Gattringer, Weyrer, Schlotzhauer and Brandstötter28]. Nonetheless, a higher accuracy seems to be achievable by considering also the dependence to load [Reference Tadese, Yumbla, Yi, Lee, Park and Moon35, Reference Wilfrido, Gabriel, Jean-Francois, Thibaut and Philippe36] and temperature [Reference Raviola, Guida, De Martin, Pastorelli, Mauro and Sorli37]. Some extensive studies by Madsen et al. offer a rather complete characterization of friction for the UR5e robot by means of a nonlinear model [Reference Madsen, Rosenlund, Brandt and Zhang38–Reference Madsen, Rosenlund, Brandt and Zhang40], which also takes into consideration joint stiffness. Other analyses related to harmonic gearboxes are carried out independently of the robotics application, as for instance in refs. [Reference Ma, Yan, Yin and Yang41, Reference Zou, Tao, Jiang, Mei and Wu42].
The electrical model of the UR5 robot is not a popular topic of study. Joint motors are scarcely analyzed in the literature and the available models are approximate [Reference Boscariol, Caracciolo, Richiedei and Trevisani43]. One investigation that approaches an electrical analysis of the robot concerns the modeling of a possible replacement to the UR5 joint actuators [Reference Lee, Lee, Woo, Lee, Lee and Ra44], but this analysis is not focused on the actual off-the-shelf manipulator. So far, one of the most elaborate and comprehensive work on the modeling of the UR5e robot can be found in ref. [Reference Madsen, Rosenlund, Brandt and Zhang40], which however does not consider the issue of estimating the energy consumption of that robot.
The aim of this manuscript is to fill this gap by proposing an approach to define a full electro-dynamic model of a UR5e robot by Universal Robots (Fig. 1). More in detail, the main novelties of this work include:
-
the modeling of the UR5e robot considering not only the dynamics of the robot but also the energy model of the motors of the manipulator;
-
a sequential procedure for the estimation of the inertia tensors of the links, followed by the identification of the moments of inertia of the motors and reducers, and finally the simultaneous identification of the friction coefficients and motor torque constants;
-
the formulation of a novel friction model accounting for the direction of the mechanical power flow within the harmonic drive reducers.
The same sequential procedure can be applied without significant differences to any of the manipulators by Universal Robots. The suitability of the procedure to any robot of the same family is ensured as no additional sensors or measures other than the ones provided by the robot controller are needed. Furthermore, to the best of authors’ knowledge, the approach to friction modeling proposed and tested in this work is different from alternative approaches found in the literature and has been chosen to retain a good model fitting with a minimal set of parameters. The resulting model can then be used to apply one of the many model-based energy-saving approaches to reduce the energy consumption and the carbon footprint of a work cell that includes one or more of such robots.
The paper is organized as follows: Section 2 provides a description of the setup used to gather the measurements on which the modeling procedure is based. Section 3 covers the identification of the inertial properties of the manipulator, whereas the friction parameters and the motor torque constants are identified by means of the procedure outlined in Section 4. The analysis of the energy model and the motor parameters are described in Section 5. Finally, the conclusions of this work and possible future developments are given in Section 6.
2. The UR5 e-series manipulator
2.1. Kinematics
The kinematics of the UR robot family has been extensively studied and discussed [Reference Schreiber and Gosselin45]; therefore, that topic is not covered in this work. The kinematic parameters, reported in Table I as a set of
$24$
Denavit-Hartenberg parameters, are assumed to be the nominal ones as disclosed by the manufacturer in ref. [46]. In other words, no kinematic calibration procedure has been performed and the robot used for the experiment is a completely ‘stock’ one.
Denavit-Hartenberg parameters of the UR5e robot.
The architecture of the UR5e robot is very similar to the one of the other manipulators produced by the same manufacturer: all of them are serial robots with six degrees of freedom (DOFs) and a non-spherical wrist. The joint design is modular, meaning that the whole UR family is based on
$5$
joint sizes, ranging from the smallest ‘0’ size to the largest ‘4’ size [47]. All joints comprise a brushless AC motor and a strain wave gear reducer, with different sizes and ratings. For this reason, the procedure outlined in this work can be applied with minimal changes to any robot of the family or even to robots from different manufacturers with similar features and design.
Experimental setup: the UR5e robot and its controller.
2.2. Robot controller and communication
The experimental setup comprises a UR5e arm by Universal Robot (Fig. 1) and a computer equipped with 32 GB RAM and an Intel Core i9 processor running Ubuntu 18.04 used for programming and data logging. The acquisition of measured data is carried out by means of the Real-Time Data Exchange (RTDE) protocol [48]. The robot is programmed in ROS (Robot Operating System) Melodic with Python and controlled through a TCP/IP connection. The RTDE protocol allows for a real-time communication between the robot and an external client, in this case the programming/data logging machine running ROS, ensuring the capability of collecting a large amount of experimental data without using any additional sensor. Data can be collected with a sampling time equal to
$2\,ms$
. The data collected and used for the identification of the electro-mechanical model of the UR5e robot include the following:
-
actual and reference joint positions;
-
actual and reference joint speeds;
-
reference joint accelerations;
-
reference joint torques;
-
actual motor currents.
The information consists of a variety of measurements collected by the robot onboard sensors (joint positions and speeds, motor currents) and data computed by the robot controller itself, such as the reference values for joint positions, speeds and accelerations. The latter are the set-point values that are fed to the robot control loop, which takes care of minimizing the tracking error of such references.
Test trajectory A: joint positions, joint and TCP speeds, joint accelerations.
Test trajectory B: joint positions, joint and TCP speeds, joint accelerations.
The motor torques are not recorded, since no torque sensor is included in the design of the robot. However, a rather good estimation of such data is provided through current measures. The reference joint torque will be discussed later, as being related, but not equal, to the actual torque available at each robot joint. In the following, it will be demonstrated that these measurements, when combined with some data disclosed by manufacturers, gathered from the literature and organized through identification processes, are sufficient to develop an accurate energy model of the UR5e robot. A similar procedure can be applied, with minimal differences, to all the robots produced by the same manufacturer, being all based on the same kinematic structure and on a standardized joint design [Reference Schreiber and Gosselin45].
Two test trajectories have been taken into account for the experimental tests: they are referred to as test trajectory A (Fig. 2), and B (Fig. 3). The related paths, referred to the position of the end-effector, are shown in Fig. 4. The trajectories are planned in the joint space using trapezoidal speed profiles, with maximum speed equal to
$3\,rad/s$
and maximum acceleration equal to
$6\,rad/s^2$
, in compliance to the joint speed limits recommended by Universal Robots, which is set to
$180^{\circ }/s$
for all joints. No specific data are available for the maximum acceleration. The motion planning method is also listed as “Linear interpolation with continuous acceleration blend” in the book [Reference Khalil and Dombre49]. For each trajectory, six random points are defined in order to ensure the robot traverses a significant portion of its workspace with multiple directional changes, exploiting its speed and acceleration capabilities to their maximum extent. The resulting trajectory is indeed very challenging to the robot for its maximum speed: the sixth plot in Figs. 2 and 3 shows that the speed of the flange of the robot is limited to roughly
$1.3\,m/s$
, unlike the planned trajectory. This feature is the result of the re-planning of the motion profile automatically operated by the robot controller. The test trajectory A is adopted for the identification of the unknown model parameters, whereas the other one is used to test the accuracy of the data gathered from the first test without any additional modification to the estimated model. In Figs. 2 and 3, the values of the joint speeds (as recorded by the robot controller) and joint accelerations (as estimated by the robot controller) are shown for the two tested motions, whereas Figs. 5 and 6 show exemplary frames of the experiments. Furthermore, to better appreciate the experimental tests, a video is available as supplementary material and online.Footnote 1
End-effector paths: (a) Test trajectory A; (b) Test trajectory B.
Exemplary frames of the experiment for test A.
Exemplary frames of the experiment for test B.
3. Identification of inertial properties
3.1. Identification of the link’s inertial properties
The robot dynamics can be written using the recursive Newton-Euler algorithm. This formulation provides a straightforward relationship between the robot’s inertial properties and the formulation of its inverse dynamics, and it is the most useful one for identification purposes [Reference Wu, Wang and You20, Reference Gaz, Cognetti, Oliva, Giordano and De Luca29, Reference Raviola, Guida, De Martin, Pastorelli, Mauro and Sorli37]. Neglecting, at this initial stage, friction forces, the vector
$\boldsymbol{\tau } = \left [ \tau _1, \ldots, \tau _6 \right ]^T$
of the six joint torques to be exerted is a function of the robot motion (through joint positions
$\boldsymbol{\vartheta }$
, velocities
$\dot{\boldsymbol{\vartheta }}$
, and accelerations
$\ddot{\boldsymbol{\vartheta }}$
) and the vector of external forces
$\textbf{f}_{ext}$
acting on it, as:
\begin{equation} \boldsymbol{\tau } = \textbf{f} \, (\boldsymbol{\vartheta },\dot{\boldsymbol{\vartheta }},\ddot{\boldsymbol{\vartheta }},\textbf{f}_{ext}) \end{equation}
The function
$\textbf{f}(\boldsymbol{\vartheta },\dot{\boldsymbol{\vartheta }},\ddot{\boldsymbol{\vartheta }},\textbf{f}_{ext})$
represents the nonlinear relationship between joint variables
$\boldsymbol{\vartheta }$
,
$\dot{\boldsymbol{\vartheta }}$
,
$\ddot{\boldsymbol{\vartheta }}$
and the joint torques
$\boldsymbol{\tau }$
, considering the effects of inertial forces, gravitational force and external forces.
The simplified robot joint model used for the identification of the inertial parameters.
The dynamic model of Eq. (1) takes into account the link dynamics only. To account for the effects of the motors and of the reducers and their inertial properties, as depicted in the model in Fig. 7, Eq. (1) should be rewritten as:
\begin{equation} \boldsymbol{\tau } = \textbf{f}_{link} \, (\boldsymbol{\vartheta },\dot{\boldsymbol{\vartheta }},\ddot{\boldsymbol{\vartheta }},\textbf{f}_{ext}) + \textbf{i} \left ( \textbf{J}_m + \textbf{J}_r \right ) \ddot{\boldsymbol{\vartheta }}_m \end{equation}
The formulation of Eq. (2) is introduced to explicitly separate the inertial effects due inertial terms due to the mass distribution of the robot links from the inertial terms due only to the moments of inertia of the motors and the reducers. The matrix
$\textbf i$
is a
$6 \times 6$
diagonal matrix of the gear ratios, while
$\textbf{J}_m$
and
$\textbf{J}_r$
are the diagonal matrices of the motors and reducers moment of inertia, assumed to be reflected on the motor side. The two terms are separated as they will be estimated independently from each other, as it will be explained in the following. The term
$\textbf{f}_{link} \, (\boldsymbol{\vartheta },\dot{\boldsymbol{\vartheta }},\ddot{\boldsymbol{\vartheta }},\textbf{f}_{ext})$
can be evaluated using the recursive Newton-Euler algorithm, in which the force
$\textbf{F}_j$
and the moment
$\textbf{M}_j$
acting on the
$j$
-th link are computed as:
\begin{equation} \textbf{F}_j = m_j{\dot{\boldsymbol{v}}}_j + \dot{\boldsymbol{\omega }}_j \times \textbf{P}_j + \left ({\boldsymbol{\omega }}_j \times \left ({\boldsymbol{\omega }}_j \times \textbf{P}_j \right ) \right ) \end{equation}
\begin{equation} \textbf{M}_j = \textbf{J}_j \dot{\boldsymbol{\omega }}_j +{\boldsymbol{\omega }}_j \times \left ( \textbf{J}_j{\boldsymbol{\omega }_j} \right ) + \textbf{P}_j \times \dot{\boldsymbol{v}}_j \end{equation}
The formulation of Eq. (3) and (4) involves some kinematic quantities, such as the linear acceleration
$\dot{\textbf{v}}_j$
, the angular speed
$\boldsymbol{\omega }_j$
of the
$j$
-th joint and its time derivative, as well as some links inertial properties, such as the overall mass
$m_j$
of the
$j$
-th link, and its first moment
$\textbf{P}_j$
, which is computed upon the position of the center of mass
$[c_{x_j}, c_{y_j},c_{z_j}]$
of the link measured in the local link reference frame as:
\begin{equation} \textbf{P}_j = m_j \left [ \begin{array}{c@{\quad}c@{\quad}c} c_{x_j} & c_{y_j} & c_{z_j} \end{array} \right ]^T \end{equation}
Equation (4) requires the evaluation of the inertia tensor
$\textbf{J}_j$
for the
$j$
-th link, which is formulated as:
\begin{equation} \textbf{J}_j = \left [ \begin{array}{c@{\quad}c@{\quad}c} J_{xx_j} & J_{xy_j} & J_{xz_j}\\ J_{xy_j} & J_{yy_j} & J_{yz_j}\\ J_{xz_j} & J_{yz_j} & J_{zz_j} \end{array}\right ] \end{equation}
$\textbf{J}_j$
is a symmetric
$3 \times 3$
matrix, so it has six distinct values. The implementation of the Newton-Euler algorithm also includes the effect of gravity in the framework of Eq. (3) and (4), by means of a fictitious ‘upwards’ acceleration of the robot base [Reference Sciavicco and Siciliano50]. Considering the UR5e robot, as many as
$72$
inertial parameters, usually called ‘base parameters’, have to be determined, of which only some of them are disclosed by the manufacturer. The most common method of estimating the unknown values is by means of a model identification technique, generally using a linear regression method [Reference Khalil and Dombre49], as done for example in ref. [Reference Raviola, Guida, De Martin, Pastorelli, Mauro and Sorli37]. The method used here is however based on a geometrical approach: as suggested in the work [Reference Khalil and Dombre49], the values of the inertia tensors of the links can be estimated by means of some simplified geometric descriptions of the links, for examples representing the robot links as composed by cylinders or some other solids. The lower accuracy due to this approach is almost irrelevant to the overall accuracy of the dynamic model, due to the high value of gear reduction found in the robot under investigation. Each joint of the robot includes an harmonic drive with a gear ratio equal to
$i=100$
, so that the reflected moment of inertia of the link as measured on the motor side has to be divided by
$10^{4}$
. Therefore, the actual values included in the inertia tensor
$\textbf{J}_j$
are scarcely relevant to the motor joint dynamics, so that they can be roughly estimated or even neglected without altering significantly the accuracy of the model tuning procedure. The most relevant contribution to the inertial forces acting on the motor are indeed the ones due to the motors and to the gearboxes, that is, the ones appearing in second term on the right-hand side of Eq. (2).
The mass of the individual links and the positions of their center of mass in the local reference frame have been disclosed by the robot manufacturer in ref. [46], and are reported in Table II.
Link masses
$m_j$
and center of mass positions for the UR5e robot.
However, the data supplied by the manufacturer do not include any information about the link’s inertia tensors. The approach used for defining the inertia tensors of the links is geometrical: knowing the size, the mass and the center of mass of a link, the inertia tensor is computed considering each link as a union of cylinders of different masses, according to the geometrical approach suggested in ref. [Reference Khalil and Dombre49] and applied to the UR5 robot in ref. [Reference Kufieta51]. The values included in the six inertia tensors are reported in Table III. These values are used for computing the joint torques collected in
$\textbf{f}_{link} \, (\boldsymbol{\vartheta },\dot{\boldsymbol{\vartheta }},\ddot{\boldsymbol{\vartheta }},\textbf{f}_{ext})$
used in Eq. (2) by means of the recursive Newton-Euler algorithm.
Link inertia tensors of the UR5e robot.
3.2. Identification of the motors and reducers inertial properties
The determination of the additional inertia values needed to complete the conservative dynamic model of Eq. (2), that is, the moments of inertia of each motor plus the reducer
$J_{m_j} + J_{r_j}$
, requires a different approach from the one used for the robot links as no data are available, nor among the data supplied from the manufacturer, nor in the literature, at least to the best of the authors’ knowledge. The proposed approach is based on the identification of the unknown parameters using experimental data gathered from the robot controller, which are labeled as target_moment, whose description by the manufacturer is ‘target joint moments (torques)’. In the absence of further specifications, and by the analysis with the results of extensive experimental tests, it is assumed that the six values contained in target_moment are the output of the inverse dynamic model, which defines the feedforward control action within the robot controller.
Our assumption is that the control scheme used by the manufacturer is mainly composed of the mentioned feedforward action, based on the ideal inverse dynamic model of Eq. (2), and by a feedback control action that compensates the inevitable discrepancies associated with friction forces, and external disturbances. This procedure is akin to the reverse engineering approach that is sometimes used for robot model identification [Reference Gaz, Flacco and De Luca52]. As a further support to this intuition, the six target_moment are directly related to the measured data labeled as joint_control_output (described by the manufacturer as ‘joint control currents’), but apparently not directly related to the measured data described as ‘actual joint currents’ [48], which are, as in the most common practice [Reference Khalil and Dombre49], assumed to be proportional to the actual torques produced by the actuators.
Now, if the six sampled torque signals mentioned above are collected in the vector
$\boldsymbol{\tau }^{ref}$
, it can be assumed that, as explained above,
$\boldsymbol{\tau }^{ref} \approx \boldsymbol{\tau }$
, hence Eq. (2) can be rewritten as:
\begin{equation} \boldsymbol{\tau }^{ref} \approx \textbf{f}_{link} \, (\boldsymbol{\vartheta },\dot{\boldsymbol{\vartheta }},\ddot{\boldsymbol{\vartheta }},\textbf{f}_{ext}) + \textbf{i} \left ( \textbf{J}_m + \textbf{J}_r \right ) \ddot{\boldsymbol{\vartheta }}_m \end{equation}
Equation (7) can be split into
$6$
individual equations, each one related to the j-th joint, given that the only unknown values are the ones appearing in the two diagonal matrices
$\textbf{J}_m$
and
$\textbf{J}_r$
. Now, recalling that for the j-th joint the vectors
$\ddot{\boldsymbol{\vartheta }}_{m_j}$
and
$\boldsymbol{\tau }_j^{ref}$
contain data collected from the robot and that the vector of joint torques
${\textbf f}_{link_j}$
has been computed using the recursive Newton-Euler formulation, the unknown values can be isolated in the left-hand side of Eq. (7) as:
\begin{equation}{\ddot{\boldsymbol{\vartheta }}_{m_j}} \left ( J_{m_j} + J_{r_j} \right ) = \dfrac{1}{i} \left [{{\boldsymbol{\tau }}^{ref}_j } -{{\textbf f}_{link_j} } \right ] \quad \text{for} \; j=1,\ldots,6 \end{equation}
Finally, the least square estimation of the motor plus reducer inertia is performed by the pseudo inversion (
$\dagger$
) of the measurement vector
$\ddot{\boldsymbol{\vartheta }}_{m_j}$
, as follows:
\begin{equation} \left ( J_{m_j} + J_{r_j} \right ) = \dfrac{1}{i}{{\ddot{\boldsymbol{\vartheta }}_{m_j}}}^{\dagger } \left [{{\boldsymbol{\tau }}^{ref}_j } -{{\textbf f}_{link_j} } \right ] \quad \text{for} \; j=1,\ldots,6 \end{equation}
This formulation merges the inertial properties of the motor and the reducer. However, these two quantities can be estimated separately only if they are physically separated, which is not possible without disassembling the robot. Another option is to spot exactly the model of the motor and the reducer used in each joint. The exact motor catalog number is difficult to determine, whereas the same does not apply to the reducer.
As disclosed by the manufacturer, the UR5e robot is composed of three ‘size 3’ joints (the first three ones) and three ‘size 1’ joints, and the maximum torque exerted by each joint is equal to
$150\,Nm$
and to
$28\,Nm$
[47]. Moreover, we know that all joints use a
$100:1$
reduction gear in the form of a HFUS-2SH strain wave gearing [53] produced by Harmonic Drive AG. Merging the information on the harmonic drive catalog with some available photos of a disassembled UR5 robot, the reducers used in ‘size 3’ and ‘size 1’ joint reducers are identified by the codes HFUS-25-2SH and HFUS-14-2SH, respectively [53]. According to the manufacturer catalog, the two reducers are rated, respectively, for a maximum torque equal to
$158\,Nm$
and
$28\,Nm$
, respectively, and their moment of inertia is equal to
$J_{r_{1,2,3}} = 1.07 \cdot 10^{-4}\,Nm$
and
$J_{r_{4,5,6}} = 0.091 \cdot 10^{-4}\,Nm$
. These values can be subtracted from the ones estimated by means of the linear regression of Eq. (8), resulting in the data reported in Table IV.
Estimated additional inertial parameters: motors and reducers moment of inertia.
The validation of the model developed up to this point can be performed by the direct comparison of the recorded control torque
${\boldsymbol{\tau }}^{ref}$
with the one evaluated by means of the inverse dynamic model of Eq. (2), which is referred to as
${\boldsymbol{\tau }}^{sim}$
. This comparison is shown in Figs. 8 and 9, which report the six values of the recorded joint torques, shown in green, with the torques evaluated by means of the dynamic model of Eq. (2), shown in orange, for the reference trajectory A. From the figures, it can be seen that the two set of values match each other almost perfectly, as the green line is almost completely covered by the orange one. This data experimentally prove that the approximation included in Eq. (7) is indeed a very minor one. The accuracy of the model is then also measured by the statistical data reported in Table V, which lists, for each joint torque, the mean torque error, the maximum absolute value of the torque error, its RMS value and the coefficient of determination
$R^2$
. The peak values of the torque error
$\tau _{err}^{max}$
are larger than
$1\,Nm$
for just the first two joints, and the RMS torque errors
$\tau _{err}^{RSM}$
take rather small values. The coefficients of determination
$R^2$
assume values that very close to
$1$
, showcasing an overall very good match between the model and the experimental data. Furthermore, Fig. 10a and 10b report the box plot representation of the torque errors for the tests A and B, respectively. The results show not only that the model fitted by means of the data of test A can faithfully reproduce the results of the same test, but that the identified model is capable of reproducing also the data collected during test B, which is not used for the identification procedure.
Frictionless model accuracy: torque modeling errors (mean, maximum and RMS errors) and coefficients of determination
$R^2$
.
Validation of estimated values of motor and gearbox moments of inertia from control torque: test A.
Validation of estimated values of motor and gearbox moments of inertia from control torque: test B.
Box plot representation of the torque errors: (a) Test A; (b) Test B. The central mark indicates the median, the bottom and top of each box represent the first and third quartiles, whereas the whiskers extend to the most extreme data not considered as outliers.
4. Friction modeling and torque constant identification
The dynamic model of the robot defined in Sect. 3 does not include any dissipative effect, which constitutes a significant approximation, particularly for a robot like the UR5e where friction effects are significant. This is mainly due to the use of a gear reduction equal to
$100:1$
and to the use of harmonic drive reducers, which are often chosen in robotics for their favorable power density and for their high reduction ratio, despite a less-than-ideal mechanical efficiency.
Given that the robot manufacturer does not provide any data related to friction, a model identification procedure from experimental data has to be used. The inverse dynamic model outlined in Sect. 3 must be extended by means of a joint friction model, so that the output of the resulting non-conservative model matches the experimental data. However, as stated above, the actual joint torque is not directly measured by the controller and therefore cannot be accessed, so an indirect measure has to be adopted. The closest match to the joint torque is the motor current, which is assumed to be proportional to the torque generated by each motor [Reference Tika, Ulmen and Bajcinca27, Reference Deng, Shang, Zhang, Zhen and Cong30, Reference Raviola, Guida, De Martin, Pastorelli, Mauro and Sorli37, Reference Madsen, Rosenlund, Brandt and Zhang38]. The torque
$\tau _{m_j}$
generated by the motor of the
$j$
-th joint can be written as:
\begin{equation} \tau _{m_j} = k_{t_j} I_j \end{equation}
being
$k_{t_j}$
the torque constant of the motor that drives the
$j$
-th joint, and
$I_j$
the current drawn on that motor. However, instead of identifying first the torque constant and then the friction model, it is advisable to perform a single identification, since both quantities rely on the same experimental data and a sequence of identification procedures usually brings an accumulation of errors [Reference Khalil and Dombre49].
As far as the choice of the friction model, several approaches are possible and different friction models can be combined together to capture the non-trivial relationship that relates friction forces to the operative conditions. In the UR5e robot, it is expected that friction forces are mainly due to the harmonic drive, given that the friction in a well-designed and maintained brushless motors is comparably low. However, the literature on friction modeling for robots and harmonic drives does not highlight a general consensus on which model provides the most accurate description of the phenomenon. In most works, such as in ref. [Reference Huang, Ke, Zhang and Ota54], friction is simply modeled as the sum of a viscous and a non-smooth Coulomb torque term. Similar options are chosen in refs. [Reference Tika, Ulmen and Bajcinca27, Reference Raviola, Guida, De Martin, Pastorelli, Mauro and Sorli37]. In ref. [Reference Kovincic, Müller, Gattringer, Weyrer, Schlotzhauer and Brandstötter28], the Coulomb friction term is smoothed by an hyperbolic function, so as to avoid a discontinuity at null speed. The same model, with the addition of a Stribeck friction torque is the choice of works [Reference Zou, Tao, Jiang, Mei and Wu42, Reference Huang, Ke, Zhang and Ota54], with the letter requiring the determination of
$24$
parameters for modeling dissipative forces. Furthermore,
$42$
parameters are determined in ref. [Reference Abedinifar, Ertugrul and Arguz55]. In ref. [Reference Gaz, Magrini and De Luca56], the joint friction is modeled by a continuous function of speed, shaped as a sigmoid function that depends on
$5$
parameters, whose identification requires a nonlinear estimation procedure. A different but accurate model for the dissipation in strain wave gears, such as the harmonic drive, suggests that friction forces depend on speed, direction of rotation, as well as on joint load: the latter is also assumed to vary among the four quadrants of operation [Reference Wahrburg, Klose, Clever, Groth, Moberg, Styrud and Ding57].
The availability of several approaches in the literature suggests that there is no obvious choice when striving for the most accurate friction model. A sensible choice should also focus on getting the best tradeoff between model complexity and accuracy, provided that a model with a larger number of parameters generally can be better fitted than a model with less tuning parameters. In this work, the model of choice differs from all the ones found in the works referenced above: we assume that friction is modeled by a ‘standard’ friction torque plus an asymmetric smooth Coulomb torque that, rather depending on the speed sign, depends on the direction of the mechanical power flow within the harmonic drive. As such, the friction torque for the
$j$
-th joint, assumed to be acting on the motor side, is evaluated by means of the following formulation:
\begin{align} \tau _j^{fric} = & k_v \dot{\vartheta }_m + \tanh \left ( \dfrac{\dot{\vartheta }_{m_j}}{\varepsilon } \right ) T_c^* \end{align}
\begin{align} \text{with} \;\;\; &T_c^* = T_c^{dir}\;\;\; \text{if} \;\; \dot{\vartheta }_{j}\tau _{j} \gt 0 \end{align}
\begin{align} &T_c^* = T_c^{rev} \;\;\; \text{if} \;\; \dot{\vartheta }_{j}\tau _{j} \lt 0 \end{align}
The constant
$\varepsilon$
is a smoothing parameter empirically set to
$10^{-3}$
, and the Coulomb term
$T_c^*$
assumes the value
$T_c^{dir}$
when the power is flowing from the motor to the load, and the value
$T_c^{rev}$
when the direction of the power is reversed. The previous expression can be rewritten as:
\begin{equation} \tau ^{frict}_j = k_{v_j}\dot{\vartheta }_m + \tanh{\left (\dfrac{\dot{\vartheta }^{\,dir}_{m_j}}{\varepsilon }\right )} T_c^{dir} + \tanh{\left (\dfrac{\dot{\vartheta }_{m_j}^{\,rev}}{\varepsilon }\right )} T_c^{rev}\;\;\;\;\;\;\;\; \text{for} \; j=1,\ldots,6 \end{equation}
with
\begin{align} \dot{\vartheta }_{m_j}^{\,dir} &= \dot{\vartheta }_{m_j} &\text{if} \;\; \dot{\vartheta }_{j}\tau _{j} \gt 0 \nonumber \\ \dot{\vartheta }_{m_j}^{\,dir} &= 0 &\text{if} \;\; \dot{\vartheta }_{j}\tau _{j} \lt 0 \end{align}
and
\begin{align} \dot{\vartheta }_{m_j}^{\,rev} &= \dot{\vartheta }_{m_j} &\text{if} \;\; \dot{\vartheta }_{j}\tau _{j} \lt 0 \nonumber \\ \dot{\vartheta }_{m_j}^{\,rev} &= 0 &\text{if} \;\; \dot{\vartheta }_{j}\tau _{j} \gt 0 \end{align}
The dissipative forces can then be included in the model of Eq. (2), and fitted with the current measurements, which represent the available data that is most closely related to the actual motor torques. Therefore, the value of the
$j$
-th motor current,
$I_j$
, the motor speed
$\dot{\vartheta }_{m_j}$
and the simulated joint torque
${ \tau }^{sim}_j$
(according to the model developed in Sect. 3) can be related by the linear combination:
\begin{equation}{I}_{j} = \frac{1}{k_{t_j}} \left ( \frac{{\tau }^{sim}_j}{i} + \dot{ \vartheta }_{m_j} k_{v_j} + \tanh{\left (\dfrac{\dot{\vartheta }^{\,dir}_{m_j}}{\varepsilon }\right )} T_c^{dir} + \tanh{\left (\dfrac{\dot{\vartheta }_{m_j}^{\,rev}}{\varepsilon }\right )} T_c^{rev} \right ) \quad \text{for} \; j=1,\ldots,6 \end{equation}
For each joint, the sequence of sampled data can be rearranged in matrix form leading to the following regressive equation:
\begin{equation} \left [ \begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \frac{{{\boldsymbol{\tau }^{sim}_j}}}{i} &{\dot{\boldsymbol{\vartheta }}_{m_j}} &{\tanh{\left (\dfrac{\dot{\boldsymbol{\vartheta }}^{\,dir}_{m_j}}{\varepsilon }\right )}} &{\tanh{\left (\dfrac{\dot{\boldsymbol{\vartheta }}_{m_j}^{\,rev}}{\varepsilon }\right )}} \end{array} \right ] \left [ \def \arraystretch{2} \begin{array}{c} \dfrac{1}{k_{t_j}} \\ \dfrac{k_{v_j}}{k_{t_j}} \\ \dfrac{T_c^{dir}}{k_{t_j}} \\ \dfrac{T_c^{rev}}{k_{t_j}} \\ \end{array}\right ] ={I_j} \;\;\;\;\;\;\;\; \text{for} \; j=1,\ldots,6 \end{equation}
The regressive form allows us to determine, through the pseudo inversion of the measurement vector, the
$5 \times 6$
coefficients reported in Table VI. This formulation provides the least square solution to the determination of the vector of parameters in Eq. (18), that is, the values of the four parameters are computed to minimize the quadratic norm of the estimation error on the motor currents
$\textbf{I}_j$
. The accuracy of the model fitting is associated with the data in Table VI can be first evaluated by analyzing Fig. 11 and 12. The plots show the measured motor currents versus the estimated ones for each joint in orange and blue, respectively.
Estimated friction coefficients and motor torque constants.
The model is then tested against the measurements associated with test B, and the comparison is shown in Fig. 12. Also in this case, the actual and estimated currents are in a very good agreement, as the modeled and the measured overlap each other well. The agreement of the model is quantified by the data reported in Table VII, which shows the mean, maximum and RMS current error, as well as the coefficient of determination
$R^2$
for each motor current and for each of the two tests. The peak current error is generally rather large, sometimes even larger than
$1\,A$
, whereas the RMS current error is far less relevant (the maximum RMS current error is indeed smaller than
$20\,mA$
for any of the first three joints and never exceeding
$66\,mA$
for the last three joints). The coefficients of determination are set mostly around
$0.99$
, with the lowest value being equal to
$0.985$
. Furthermore, Fig. 13a and 13b report the box plot representation of the current errors for the tests A and B, respectively. It can be stated that the estimated model provides a faithful representation of the current absorption associated with the execution of a task.
Estimated current and measured current comparison: test A.
Current model accuracy: current modeling errors (mean, maximum and RMS errors) and coefficients of determination
$R^2$
.
Estimated current and measured current comparison: test B.
Box plot representation of the current errors: (a) test A; (b) test B. The central mark indicates the median, the bottom and top of each box represent the first and third quartiles, whereas the whiskers extend to the most extreme data not considered as outliers.
5. Electric power consumption modeling
In this section, the electric model of the robot is discussed. The main challenges in the identification of a precise electric model are the unavailability of experimental data other than the six motor currents, as well as the limited technical data provided by the manufacturers. The robots by Universal Robots are equipped with custom brushless AC motors manufactured by Kollmorgen GMBH, belonging to the AKM frameless motor series [58]. Browsing the AKM motor series datasheet [59] highlights that the motors used in the UR5e robot must be a custom one, as the AKM motors are operated with supply voltages rated at least
$240\,V\,AC$
, whereas the motors mounted on the robot are powered from a
$48\,V\,DC$
bus, as disclosed by the robot manufacturer [48]. However, the same motor manufacturer offers some low-voltage motors specifically targeted to collaborative robots, such as the TBM2G frameless motor series [60], which are designed to be operated from a
$48\,V$
DC bus exactly like the one sported by the UR5e robot. There are no other clues indicating that the actual motors used in the UR5e manipulator are identical to the ones available in the catalog [60], but it is likely that the motors are based on similar geometries and winding arrangements, hence they are expected to sport similar electrical properties. Their determination is useful in the definition of the electric power absorption model that is discussed here.
The Park-Clark transformation is used to represent AC brushless motors as equivalent DC motors [Reference Park61]. This transformation treats a set of sinusoidal AC currents and voltages, namely the one absorbed by the phases of the motor, as a single equivalent DC voltage and current. That current is therefore interpreted as the current drawn in the equivalent DC motor, for which a simple electro-mechanical model is available. According to it, the equivalent voltage drop across the j-th motor leads and the current absorbed by the same motor can be expressed as:
\begin{equation} V_j = R_{m_j} I_j + k_{b_j} \dot{\vartheta }_{m_j} \end{equation}
\begin{equation} I_j = \dfrac{\tau _{m_j}}{k_{t_j}} \end{equation}
The model of Eq. (19) and (20) highlights that the current
$I$
is proportional to the generated torque, and the voltage drop
$V$
is proportional to the voltage drop across the equivalent motor winding resistance
$R_m$
plus the one due to the speed-controlled voltage generator by means of the coefficient of proportionality this is usually referred to as the back-emf constant,
$k_b$
. The power absorbed by the six motors can be defined by determining the parameters of Eq. (19) and (20). However, the only parameter that can be actually established with precision from the available data is the torque constant
$k_t$
, as discussed in Sect. 4. The instantaneous motor power consumption
$W_{e_j}$
by the j-th joint motor can be determined by the voltage-current product as:
\begin{equation} W_{e_j} = V_jI_j = R_{m_j} I_j^2 + k_{b_j} \dot{\vartheta }_{m_j} I_j \end{equation}
However, the substitution of Eq. (20) in Eq. (21) provides a more practically useful relation for computing the electric energy absorption
$W_{elec_j}$
for the
$j$
-th motor:
\begin{equation} W_{{elec}_j} = \dfrac{R_{m_j}}{k_{t_j}^2} \tau _{m_j}^2 + \dfrac{k_{b_j}}{k_{t_j}} \dot{\vartheta }_{m_j} \tau _{m_j} \end{equation}
Equation (22) can be further simplified, by noticing that the ratio
${k_t}/{\sqrt{R_m}}$
is the so-called motor constant, and that the ratio for the back-emf constant and the torque constant is theoretically unitary for the convention used here. A brushless motor like the ones used in the UR5e robot is supplied by a three-phase AC current, which is modulated across three windings on the motor stator, whereas the model that underlies Eq. (22) is developed for a DC (line to neutral) voltage and a line current. Therefore, the parameters
$R_m$
,
$k_b$
and
$k_t$
that are needed for our purposes must be converted from the ones used in the manufacturer’s catalog, which usually refers to a different convention [62]. Following the indications found in ref. [62], we can obtain the correct values for the constants in Eq. (22), for each motor in the catalog [60]. In particular, the use of a proper conversion highlights that the
$k_t/k_b$
ratio is equal to roughly
$0.89$
for any motor belonging to the catalog used for reference. The ratio
$k_t/k_b$
is frequently regarded as a unitary value; however, this assumption holds merely if the only source of loss in the motor is the one due to the Joule effect, neglecting the significant contributions of core losses, stray load losses, and other dissipative effects [Reference Paryanto, Brossog, Bornschlegl and Franke63]. In a non-ideal case, the
$k_t$
value should therefore be smaller than the
$k_b$
value, as done here using the ratio
$0.89$
suggested by the manufacturer. Using the motor constant
$k_m$
, Eq. (22) can be rewritten as:
\begin{equation} W_{{elec}_j} = \dfrac{1}{k_{m_j}^2} \tau _{m_j}^2 + \dfrac{k_{b_j}}{k_{t_j}} \dot{\vartheta }_{m_j} \tau _{m_j} \end{equation}
Therefore, only the motor constant
$k_m$
, measured in
$Nm/\sqrt{W}$
, has to be identified to model the energy consumption of a motor, assuming that the ratio
$k_t/k_b$
has been measured by the manufacturer. Given the impossibility of directly measuring this parameter without disassembling the robot, a sensible estimation has been performed by identifying the most probable replacements using the motors available in the TBM2G catalog [60], which lists several motors with external diameters ranging from
$50\,mm$
to
$115\,mm$
. Judging from the external size of the motor housings on the robot, size
$1$
motors can be of either
$50\,mm$
or
$65\,mm$
external diameter, whereas motors of size
$3$
can be either of
$68\,mm$
or
$76\,mm$
in diameter. Motor lengths can be, for all possible diameters,
$8\,mm$
,
$13\,mm$
or
$26\,mm$
. The selection of the most similar motor is to be done also on the basis of the motor parameters estimated so far, such as the motor moment of inertia, the torque constant and the maximum torque generation capability (declared by the robot manufacturer). An identification of the specific motor code is however not needed, since the motor constant
$k_m$
is independent from the winding geometry and depends only on the overall motor dimensions. As a matter of fact, it identifies the ‘motor size’ and therefore the suggested procedure needs to guess a proper motor size among a list of possible motors.
Using all the available data, one possible candidate for the motors of size
$1$
(i.e., the ones for joints
$4$
,
$5$
and
$6$
of the UR5e) is the one of the three TBM2G-05008x motors (‘x’ might be A, B or C according to the winding geometry), whose nominal constant torque ranges between
$0.054\,Nm/A$
and
$0.0266\,Nm$
(indicated within a
$\pm 10\%$
range in the catalog), and their moment of inertia are all equal to
$J_m = 0.79\cdot{10}^{-5}\,kg\,m^2$
. These values match well the expected values of Tables IV and VI, with mean
$k_t = 0.0768\,Nm/A$
and mean
$J_m = 1.2032\cdot{10}^{-5}\,kg\,m^2$
. The difference between the value of the motor moment of inertia taken from the datasheet and the estimated value, which is slightly higher, can be attributed to the fact that the observed motor inertia also takes into account the moment of inertia of the coupling joint between the motor and the reducer, whose characteristics are unknown having not disassembled the robot. The relative motor constant is declared to be
$0.06632\pm 10\%\,Nm/\sqrt{W}$
: this value is again shared by all three motors BM2G-05008x, since it is independent from the winding geometry.
The mean expected values (Tables IV and VI) for the motors of size 3 are
$J_m = 6.4013\cdot 10^{-5}\,kg\,m^2$
and
$k_t=0.1061\,Nm$
: such data are fitted at best by the group of three motors labeled as motor TBM2G-06826x, whose nominal moment of inertia is equal to
$5.18\cdot 10^{-5}\,kg\,m^2$
and whose torque constant can vary between
$0.0835\pm 10\%\,Nm/A$
and
$0.2947\pm 10\%\,Nm/A$
. The motor constant is declared to be
$0.27426\,Nm/\sqrt{W}$
for all three motor samples.
This possible choice of motors as suitable equivalents to the ones actually mounted on the UR5e robot is further supported by the fact that their peak continuous service torque is equal to
$0.27\,Nm$
for the smaller motor and to
$1.54\,Nm$
for the larger ones: such values reflect almost perfectly the peak torques available at the joints declared by Universal Robots which are set to
$28\,Nm$
and
$150\,Nm$
, respectively, which are
$100$
times larger due to the gear ratio.
Motor power and energy consumption for each joint motor, test A.
Motor power and energy consumption for each joint motor, test B.
The electric model of the robot is now complete, and the inverse dynamic model can be complemented with Eq. (23) to compute the instantaneous values of the electric power absorbed by each motor: the overall energy consumption can then be estimated by computing the time integral of the electric power. The power and energy absorption for each individual motor is represented graphically in Fig. 14, which refers to the execution of the test trajectory A. The same procedure is then applied to the other test trajectory: the power and energy outlines for test B are shown in Fig. 15.
Estimated energy consumption and trajectory efficiency for each single joint and for the entire robot.
The values of estimated energy consumption associated with the two test trajectories are reported in Table VIII, which shows both the mechanical and the electrical energy consumption, for each single joint and for the entire robot. The energy measures are then used to evaluate the reported values of the trajectory efficiency
$\eta$
: this value is evaluated as suggested in ref. [Reference Park61], as the ratio between the mechanical energy
$E_{mech}$
and the electric energy
$E_{elec}$
expenditures:
\begin{equation} \eta\; :\!=\; \dfrac{E_{mech}}{E_{elec}} \end{equation}
with
\begin{equation} E_{mech} = \int \limits _{0}^T \sum _{j=1}^{6} \dot{\vartheta }_{m_j} \tau _{m_j} \, dt \end{equation}
and
\begin{equation} E_{elec} = \int \limits _{0}^T \sum _{j=1}^{6} W_{{elec}_j} \, dt \end{equation}
6. Conclusions
The electro-mechanical modeling and identification of the UR5 e-series manipulator has been presented in this paper. The model is defined upon both experimental data acquired from the robot, using model identification procedures, and datasheet provided by the manufacturers of the manipulator itself, of the motors and of the gearboxes. The model validation is supported by the precise reproduction of the motor currents by means of the developed inverse dynamic model. Friction effects have been modeled by taking into account the non-symmetric behavior sported by the strain wave gears when operating in direct or reverse power flow conditions. The electric model identification relies mostly on the information provided by the motor manufacturer, given the unavailability of specific electric measurements. Nonetheless, the identified model provides an accurate estimation of the effective energy absorption by the robot during its operation. In future developments of this work, we plan to adopt the proposed electro-mechanical model of the UR5e robot for the model-based optimization of energy-efficient trajectories to improve sustainability and reduce the power consumption in real manufacturing scenarios.
Author contributions
Methodology: EC, LS, and PB; Validation: EC, LS, and PB; Investigation: EC, LS, and PB; Writing – Original draft: EC, LS, and PB; Writing – Review Editing: LS, PB, and AG; Supervision: PB and AG. All authors have read and commented on the final version of the manuscript.
Financial support
This research has been partially developed within the Laboratory for Big Data, IoT, Cyber Security funded by Friuli Venezia Giulia region and the Laboratory for Artificial Intelligence for Human-Robot Collaboration (AI4HRC) funded by Fondazione Friuli and partially supported by iNEST-Interconnected NordEst Innovation Ecosystem, funded by PNRR (Mission 4.2, Investment 1.5), NextGeneration EU – Project ID: ECS 00000043.
Competing interests
The authors declare no competing interests.
Ethical approval
None.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1017/S0263574724000833.
Open access
