Abstract-In this paper, we address the problem of extracting a feasible set of dynamic parameters characterizing the dynamics of a robot manipulator. We start by identifying through an ordinary least squares approach the dynamic coefficients that linearly parametrize the model. From these, we retrieve a set of feasible link parameters (mass, position of center of mass, inertia) that is fundamental for more realistic dynamic simulations or when implementing in real time robot control laws using recursive Newton-Euler algorithms. The resulting problem is solved by means of an optimization method that incorporates constraints on the physical consistency of the dynamic parameters, including the triangle inequality of the link inertia tensors as well as other user-defined, possibly nonlinear constraints. The approach is developed for the increasingly popular Panda robot by Franka Emika, identifying for the first time its dynamic coefficients, an accurate joint friction model, and a set of feasible dynamic parameters. Validation of the identified dynamic model and of the retrieved feasible parameters is presented for the inverse dynamics problem using, respectively, a Lagrangian approach and Newton-Euler computations.

I. INTRODUCTION

The knowledge of accurate dynamic models is of fundamental importance for many robotic applications. It is necessary, in fact, for designing control laws with superior performance, in free motion or when interacting with the environment [START_REF] Siciliano | Robotics: Modeling, Planning and Control[END_REF], e.g., in strategies for the sensorless detection, isolation and reaction to unexpected collisions [START_REF] Haddadin | Robot collisions: A survey on detection, isolation, and identification[END_REF] or when regulating force or imposing a desired impedance control at the contact [START_REF] Magrini | Hybrid force/velocity control for physical human-robot collaboration tasks[END_REF].

In order to obtain an estimation of the dynamic model, regression techniques are widely employed for industrial [START_REF] Hollerbach | Model identification[END_REF], [START_REF] Janot | A generic instrumental variable approach for industrial robot identification[END_REF] or humanoid robots [START_REF] Jovic | Humanoid and human inertia parameter identification using hierarchical optimization[END_REF], [START_REF] Bonnet | Optimal exciting dance for identifying inertial parameters of an anthropomorphic structure[END_REF]. These techniques are hinged on a fundamental property: the linear dependence of the robot dynamic equations in terms of a set of ρ dynamic coefficients π R ∈ R ρ [START_REF] Gaz | Extracting feasible robot parameters from dynamic coefficients using nonlinear optimization methods[END_REF], also known in the literature as base parameters [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF], which are linear combinations of the dynamic parameters of each link composing the robot. In particular, each link has 10 parameters, specifying the mass, the position of the center of mass (CoM), and the 6 elements of the symmetric inertia tensor. Then, a robot with l links has a total 10 l of such parameters, denoted as p ∈ R 10 l . In addition, one may also include a number of parameters for modeling joint friction.

The regrouping of dynamic parameters, i.e., the dynamic coefficients, occurs because some parameters are not excitable during motion (they have no influence on the robot dynamics), while some others have an effect on the dynamics only in combinations (they are not separately identifiable).

The identification of the dynamic coefficients π R is often sufficient for many robotic applications, such as dynamic motion robot control and motion planning, since knowledge of π R allows for a numerical evaluation of the robot dynamic model in the Euler-Lagrange (E-L) form. However, the retrieval of a set of feasible numerical values for the dynamic parameters p is also relevant. This is the case, for instance, when performing dynamic simulations via a CAD-based robotic simulator -like V-REP [START_REF]V-REP virtual robot experimentation platform[END_REF] -or when implementing torque-level control laws (such as the feedback linearization) under hard real-time constraints. In this case, a widely adopted solution is to use the recursive numerical Newton-Euler (N-E) algorithm, which is preferred to the evaluation of the symbolic computationally more expensive E-L approach (which relies on dynamic coefficients). However, usual N-E routines require the knowledge of the dynamic parameters p of each link in the kinematic chain, and not just of the dynamic coefficients π R . In [START_REF] Gaz | Extracting feasible robot parameters from dynamic coefficients using nonlinear optimization methods[END_REF], we addressed the problem of recovering a complete set of values for the original robot parameters starting from the identified dynamic coefficients. In general, this is a nonlinear problem admitting an infinite number of solutions. However, not all solutions are physically consistent (as example, negative masses may appear). In order to discard unfeasible solutions, we considered upper and lower bounds on each component of p by solving a constrained nonlinear optimization problem. Because of the ill-conditioned nature of the solution space, we used global optimization methods, such as simulated annealing.

The physical consistency of the identified dynamic parameters, as introduced in [START_REF] Mata | Dynamic parameter identification in industrial robots considering physical feasibility[END_REF], is currently attracting more and more attention. Researchers have treated the problem of physical feasibility within the framework of linear matrix inequalities (LMIs), solving the problem of the identification of a physically consistent set of parameters by means of semidefinite programming (SDP) techniques [START_REF] Sousa | Physical feasibility of robot base inertial parameter identification: A linear matrix inequality approach[END_REF]. Recently, this framework has been enriched by the addition of the triangle inequality of the inertia tensors [START_REF] Wensing | Linear matrix inequalities for physically-consistent inertial parameter identification: A statistical perspective on the mass distribution[END_REF], [START_REF] Sousa | Inertia tensor properties in robot dynamics identification: A linear matrix inequality approach[END_REF], a constraint which was originally mentioned in [START_REF] Traversaro | Identification of fully physical consistent inertial parameters using optimization on manifolds[END_REF]. The approach presented in this paper can be considered as an alternative to the LMI-SDP framework.

All these approaches act on the parameter identification phase, obtaining from this dynamic coefficients that are typically different from the classic ordinary least squares (OLS) solution [START_REF] Hollerbach | Model identification[END_REF], [START_REF] Janot | A generic instrumental variable approach for industrial robot identification[END_REF]. In any event, the approach presented in [START_REF] Sousa | Physical feasibility of robot base inertial parameter identification: A linear matrix inequality approach[END_REF], [START_REF] Sousa | Inertia tensor properties in robot dynamics identification: A linear matrix inequality approach[END_REF] requires to express constraints as linear matrix inequalities, while the proposed optimization algorithm manages both linear and nonlinear (e.g., if-else) constraints, without any mathematical manipulation. This additional flexibility allows to handle directly nonlinear constraints coming from the geometric shape of cylindrical or spherical links (like for Universal Robots manipulators or for the sixth link of the KUKA LWR IV+), or when the use of approximate box constraints may generate solutions which are even unfeasible (e.g., a center of mass outside a convex link of the robot). On the other hand, a drawback of the proposed algorithm is that convergence to a global optimum in a finite number of steps cannot be guaranteed. The proposed approach is general and can be applied to a large class of robot manipulators. In this work, as case study, we apply it to the Franka Emika Panda robot (see Fig. 1), a manipulator that is attracting a large interest in the robotics and industrial communities due to its high usability and relatively low price among the torque-controlled manipulators. For this robot we obtain a complete identification of a set of feasible dynamic parameters.

Summarizing, the main contributions of this paper are: (i) presentation of a framework for the robot dynamic parameters retrieval that deals with linear, nonlinear and conditional constraints; (ii) identification of the dynamic model of the Franka Emika Panda robot1 ; (iii) retrieval of a feasible set of dynamic parameters.

The paper is organized as follows. In Sec. II, we briefly present the main features of the Franka Emika Panda robot. In Sec. III we recall the general procedure for identifying the dynamic coefficients, and we present the problem of physical consistency of the parameters. The identification procedure used to retrieve the dynamic model for this robot is described in Sec. IV, together with an estimation procedure of the joint FCI controller libfranka user input ( , , ) Fig. 2. Signal flows from and to the controller. The user sends a command to the libfranka interface that communicates with the FCI controller. This input is then converted to a commanded torque τ c to the robot that returns the measured joint torque τ , as well as the joint positions q and velocities q. The FCI controller computes the numerical values for the inertia matrix M (q), as well for the gravity vector g(q), the Jacobian J(q), and the Coriolis term c(q, q). These data are sent back to the user through the libfranka interface. A more detailed description of the FCI can be found at: https://frankaemika.github.io/docs/index.html. friction. Retrieval of the feasible parameters from the dynamic coefficients is presented in Sec. V. Finally, validation results are reported in Sec. VI and conclusions are drawn in Sec. VII.

II. THE FRANKA EMIKA PANDA ROBOT

Figure 1 shows the Franka Emika Panda robot and its kinematic parameters according to the modified Denavit-Hartenberg convention. This robot is equipped with n = 7 revolute joints, each mounting a torque sensor, and it has a total weight of approximately 18 kg, having the possibility to handle payloads up to 3 kg. It is possible to control the robot through the Franka Control Interface (FCI), that is able to provide, via the libfranka interface (at 1 kHz), the joint positions q and velocities q, as well as link side torque vector τ . Moreover, it returns the numerical values of the inertia matrix M (q), of the gravity vector g(q), as well as the Jacobian J (q) and the Coriolis term c(q, q) at a given joint position q and velocity q. These data will be of fundamental importance for the identification of the dynamic model of the robot (see Section IV-A for details).

The robot can be controlled in different modalities, according to the user requirements: torque-mode (by providing a vector τ d to the robot motors), position-mode (by giving a desired joint position q d ), velocity-mode (by sending a desired joint velocity vector qd ) are the most common modalities. The FCI controller is designed in such a way that the command inputs given by the user are appropriately manipulated so that the motors generate the proper torque τ c for the commanded task. Figure 2 depicts all the control signals above-mentioned.

III. PRELIMINARIES

A. Building the inverse dynamic model

In order to derive the symbolic dynamic model of a robot with elastic joints, such as the Franka Emika Panda, one may follow the procedure presented in [START_REF] Spong | Modeling and control of elastic joint robots[END_REF], separating the motor torques from the link-side torques. Nevertheless, the particular features offered by the robot controller allow us to simplify the modeling: in fact, since the FCI controller is able to return the estimations of the link-side torques (exploiting the motor position measures read from the encoders), we are able to adopt the classical model structure as for a rigid joints robot, neglecting the elasticity [START_REF] Gaz | Identifying the dynamic model used by the KUKA LWR: A reverse engineering approach[END_REF], and henceforth τ ∈ R n is the vector of the link-side torques.

From the E-L equations [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF], we can obtain the dynamic model of a n-dof robot as M (q)q + S(q, q) q + g(q) = τ ,

where q, q, q ∈ R n are, respectively, the joint positions, velocities and accelerations, M (q) ∈ R n×n is the inertia matrix, g(q) ∈ R n is the gravity vector and S(q, q) q = c(q, q) ∈ R n is the vector of the Coriolis and centrifugal forces. The dynamic model in the form (1) includes typically nonlinear functions of q, q, q and the dynamic parameters described in detail further. For each link i , i = 1, . . . , n, let m i be the mass and let

i r i,ci =   c ix c iy c iz   , i J i =   J ixx J ixy J ixz J ixy J iyy J iyz J ixz J iyz J izz   , (2)
be the position of the center of mass and the symmetric inertia tensor with respect to the i-th link frame, respectively. At this stage, if we collect the dynamic parameters of all the robot links in the three vectors

p 1 = m 1 . . . m n T , p 2 = c 1x m 1 c 1y m 1 c 1z m 1 . . . c nx m n c ny m n c nz m n T , p 3 = J T 1 . . . J T n T , (3)
with p 1 ∈ R n , p 2 ∈ R 3n , p 3 ∈ R 6n and J i = J ixx J ixy J ixz J iyy J iyz J izz T , (4)
it is possible to rearrange (1) as

Y (q, q, q) π(p 1 , p 2 , p 3 ) = τ , (5)
where the vector π(p

1 , p 2 , p 3 ) = p T 1 p T 2 p T 3 T ∈ R p .
Moreover, π appears linearly in the dynamic model ( 5), multiplied by the regressor matrix Y of known time-varying functions.

The dynamic identification procedure is performed by collecting M n p joint torque samples as well as M joint position samples, while the joint velocity and the acceleration are computed by off-line differentiation. For each numerical sample (τ k , q k , qk , qk ), with k = 1, . . . , M , we have

Y k (q k , qk , qk )π = τ k . ( 6
)
By stacking these quantities in vectors and matrices, one has

Y π = τ , (7)
with τ ∈ R M n and Y ∈ R M n×p . According to [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF], we can prune the stacked regressor Y so as to obtain a matrix with full column rank Y R , and then identify the dynamic coefficients by solving an ordinary least-squares (OLS) problem via

pseudoinversion πR = Y # R τ . (8)
With the solution πR ∈ R ρ of regrouped dynamic parameters, i.e. the dynamic coefficients, we can provide a joint torque estimate as τ = Y R (q, q, q) πR (9)

for validation on any new motion q(t). Finally, following [START_REF] Gaz | Extracting feasible robot parameters from dynamic coefficients using nonlinear optimization methods[END_REF], one can extract from the identified vector πR a feasible set of dynamic parameters p = (p 1 , p2 , p3 ) -not necessarily the true ones -such that π R (p 1 , p2 , p3 ) = πR and the upper/lower bounds on the components of p i , i = 1, 2, 3, are also satisfied. However, the triangular inequality constraint of inertia tensors is not taken into account in [START_REF] Gaz | Extracting feasible robot parameters from dynamic coefficients using nonlinear optimization methods[END_REF], as instead done in Sec. V.

B. Physical consistency of the dynamic parameters

The obtained estimation of the dynamic coefficients vector πR might be possibly physically inconsistent (e.g., a negative link mass), and this can be caused, for instance, by modeling errors or by noisy measurements. Recent works [START_REF] Wensing | Linear matrix inequalities for physically-consistent inertial parameter identification: A statistical perspective on the mass distribution[END_REF]- [START_REF] Traversaro | Identification of fully physical consistent inertial parameters using optimization on manifolds[END_REF] highlighted these physical constraints and provided frameworks to consider them during the identification phase, by solving linear constrained optimization problems using the following cost function:

min π f (π) = Y π -τ . ( 10
)
Physical constraints regard the mass of each link, which has to be positive, and the barycentric inertia tensor of each link, which has to be positive definite. That is, for each link i of the manipulator, one has:

m i > 0 (11)
and

I i =   I ixx I ixy I ixz I ixy I iyy I iyz I ixz I iyz I izz   0, (12)
where I i is the inertia tensor of link i with respect to its center of mass. Moreover, it is always possible to express the barycentric inertia tensor in a diagonal form, exploiting a particular rotation matrix Ri , such as

I i = Ri Ī i RT i , (13)
where Ī i is the diagonal inertia tensor. Since the diagonal elements ( Īi,x , Īi,y , Īi,z ) of Ī i are also the eigenvalues of I i , condition [START_REF] Sousa | Physical feasibility of robot base inertial parameter identification: A linear matrix inequality approach[END_REF] can be rewritten as

Īi,x > 0 , Īi,y > 0 , Īi,z > 0 . ( 14
)
These three inequalities are however included in the following triangle inequality:

   Īi,x + Īi,y > Īi,z Īi,y + Īi,z > Īi,x Īi,z + Īi,x > Īi,y (15)
which, with simple manipulations [START_REF] Sousa | Inertia tensor properties in robot dynamics identification: A linear matrix inequality approach[END_REF], leads to the following condition on I i :

tr(I i ) 2 -λ max (I i ) > 0, (16)
since Īi,x + Īi,y + Īi,z = tr(I i ) and denoting with tr(I i ) and λ max (I i ), respectively, the trace and the maximum eigenvalue of the inertia tensor I i . Therefore, in order to have physical consistency, conditions ( 11) and ( 16) must be satisfied for each link.

IV. IDENTIFICATION PROCEDURE A. Identifying the model used by the controller

Exploiting the features offered by the controller of the Panda robot (see Sec. II), it is possible to avoid the classical procedure involving exciting trajectories [START_REF] Swevers | Dynamic model identification for industrial robots[END_REF], and obtain the estimates of the gravitational and inertial coefficients by collecting a set of static positions only by means of a reverse engineering procedure 2 . The same procedure had been used to retrieve the dynamic coefficients of the KUKA LWR robot [START_REF] Gaz | Identifying the dynamic model used by the KUKA LWR: A reverse engineering approach[END_REF], [START_REF] Jubien | Dynamic identification of the Kuka LightWeight robot: Comparison between actual and confidential Kuka's parameters[END_REF]: in that case, the gravitational and the inertial coefficients have been estimated separately, while now a slightly different approach has been used, due to the fact that many coefficients can be retrieved both from the inertia matrix and from the gravity vector. As a first operation, one has to rearrange the symbolic inertia matrix M (q) in a vector form (the inertia stack), by exploiting its symmetry. Having for the Panda robot n = 7 joints, we obtain a vector m(q) ∈ R m , with m = n(n + 1)/2 = 28 components, containing all the lower triangular elements of M (q). Now, it is possible to obtain -as described in Sec. III -the symbolic regressor Y s (q) from the column vector s(q) ∈ R m+n in such a way that

s(q) = m(q) g(q) = Y s (q)π s , (17)
where π s ∈ R 10n is the vector containing both the gravitational and inertial dynamic parameters, which are the sameexcluding joint friction and motor inertias -as in vector π of eq. ( 7).

In order to obtain a numerical estimation of the dynamic coefficients vector, a data acquisition procedure should be carried out: in the general case, performing exciting trajectories is required in order to span all the admissible joint positions, velocities and accelerations (see eq. ( 8)). Since s(q) depends only on the joint positions q, it is just sufficient to retrieve data by imposing static positions.

The libfranka software provides the numerical evaluation of the gravity vector and the inertia matrix of the Panda robot at the current link position. Therefore, it is possible to collect a fair amount of data even only in a static way, i.e., bringing the manipulator to a desired configuration and then retrieving and storing the numerical values of the gravity vector and the inertia matrix. This acquisition procedure can be performed during a motion as well. The main advantage of imposing static joint positions is that this procedure avoids any influence of friction and uncertainty (e.g., due to measure noise or to any unmodeled phenomenon).

The data is acquired and collected in a list of M different (special and/or random) configurations, under the weak condition M n p.

For a generic configuration q k , with k = 1 . . . M , we have

s k = m k g k = Y sk π s , (19)
2 See the supplementary material accompanying this paper for a comparison of the dynamic coefficients estimated through the reverse engineering process with those obtained from the classical approach that exploits exciting trajectories.

where m k and g k represent, respectively, the numerical inertia stack vector and the numerical gravity vector, as they are retrieved from the libfranka interface at a given configuration q k , and Y sk = Y s (q k ) is the evaluated regressor.

When all data are collected, they can be stacked into a vector s such as:

s =      s 1 s 2 . . . s N      =      Y s1 Y s2 . . . Y sN      π s = Y s π s . (20)
Nevertheless, the regressor Y s is typically rank-deficient: this implies that the elements of the vector π s are not fully identifiable. Therefore, it is necessary to drop linear dependent columns of the regressor in order to reach a full (column) rank condition (i.e., by means of the Gauss-Jordan elimination technique). As a consequence, some dynamic parameters will be grouped together accordingly, in the form of dynamic coefficients [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF]. Exploiting condition [START_REF] Swevers | Dynamic model identification for industrial robots[END_REF] on the minimal number M of samples to retrieve, the ill-conditioning of the matrix is avoided. Denoting as πs,R the vector containing the regrouped parameters (a.k.a. dynamic coefficients), and as Y s,R the full rank numerical regressor, eq. ( 20) is solved using a least squares technique as

πs,R = Y T s,R Y s,R -1 Y T s,R s = Y # s,R s, (21)
where # denotes pseudoinversion. Once one has the dynamic coefficients estimation πs , from eq. ( 17), it is possible to obtain the estimates ĝ(q) and M (q) as:

ŝ(q) = m(q) ĝ(q) = Y s,R (q) πs,R ,

where Y s,R (q) is the symbolic regressor pruned of the dependent columns (according to the full-rank matrix Y s,R ) and M (q) is built from the estimated inertia stack m(q). Finally, the estimation of the Coriolis and centrifugal forces vector ĉ(q, q) is derived from M (q) using the Christoffel's symbols according to [START_REF] Siciliano | Robotics: Modeling, Planning and Control[END_REF]. The form of the inverse dynamics formula, providing an estimation of the joint torques τ needed to accomplish a given trajectory (q, q,q), is therefore: τ (q, q, q) = M (q)q + ĉ(q, q) + ĝ(q). (23)

B. Friction estimation

If a validation trajectory is executed and the measured torques τ are compared (after a filtering procedure) with the estimated torques τ generated by eq. ( 23), a difference in the two signals may be observed. This discrepancy is due to estimation errors (e.g., due to the noise affecting the measurements) or to unmodeled effects, such as joint friction. Typically, the latter effect, assumed to act separately on each joint, is expressed as an additional torque τ f,j , j ∈ {1, . . . , n}:

τ f,j ( qj ) = f v,j qj + f c,j sign( qj ) + f o,j , (24)
where f v,j and f c,j represent, respectively, viscous and Coulomb friction, while f o,j is the Coulomb friction offset.

Model (24), although being simple and effective, has the main drawback of exhibiting sudden discontinuities in the neighborhood of qj = 0. In order to attenuate this chattering, a sigmoidal friction model can be used for avoiding discontinuities for low joint velocities. In case the viscous effects are negligible and a symmetric behavior for positive and negative joint velocities is observed (as for the Panda robot), τ f,j ( qj ) can be expressed as the following function:

τ f,j ( qj ) = ϕ 1,j 1 + e -ϕ2,j ( qj +ϕ3,j ) - ϕ 1,j 1 + e -ϕ2,j ϕ3,j , (25)
which is characterized by 3 parameters for each joint (ϕ 1,j , ϕ 2,j , ϕ 3,j ). In order to estimate the 3n = 21 friction parameters for the Panda robot, trajectories spanning all possible joint velocities can be executed for each joint (possibly, keeping the others at rest). Computing the inverse dynamics for the previous trajectories according to (23), the difference ∆τ j = τ j -τj (where τ j and τj are, respectively, the measured and the estimated joint torques) can be interpreted as a measure of the friction effects. Solving a least squares problem, it is possible to find the parameters that make the curve (25) to fit data at best -e.g., by means of a Nelder-Mead routine.

The presented reverse engineering approach has the main advantage of allowing an estimation of the dynamic parameters by exploiting only static measures (which are affected by low noise) without the need of numerical derivations, which can dramatically affect the results of the identification process. Moreover, having a separate step for joint friction identification allows to repeat only this step (instead of the whole identification process) when lubrication changes or when moving to a new instance of the same robot. Conversely, we must rely on the numerical values of s k provided by the manufacturer.

V. RETRIEVAL OF FEASIBLE PARAMETERS

In [START_REF] Gaz | Extracting feasible robot parameters from dynamic coefficients using nonlinear optimization methods[END_REF] a framework has been presented for retrieving a feasible set of dynamic parameters from the previously identified dynamic coefficients, by solving a nonlinear global optimization problem.

For some purposes, in fact, the estimated dynamic coefficients (linear combination of dynamic parameters, such as masses, position of the centers of mass and inertia tensors) are not sufficient, and an estimation of the dynamic parameters themselves is needed. This occurs, for instance, if the dynamic parameters values are needed for simulating the robot behavior in a CAD software (like V-REP), or if a Newton-Euler routine is required to compute the joint torques estimation during the robot motion, with strict time constraints (in this case, a Lagrangian approach may not fit the real-time requirements), or in case one is interested in computing the wrenches acting on the robot joints.

The approach presented in [START_REF] Gaz | Extracting feasible robot parameters from dynamic coefficients using nonlinear optimization methods[END_REF] is able to return a possible and feasible set of the dynamic parameters: in general, infinite solutions exist (among which the real one) and the more constraints are given, the closer to the real solution is the returned one.

Starting from the identified vector of dynamic coefficients πR , the following transformations may be applied to the corresponding symbolic vector π R (p) (parallel axis theorem, see [START_REF] Siciliano | Robotics: Modeling, Planning and Control[END_REF]):

J ixx → I ixx + m i c 2 iy + m i c 2 iz J ixy → I ixy -c ix c iy m i J iyy → I iyy + m i c 2 ix + m i c 2 iz J ixz → I ixz -c ix c iz m i J izz → I izz + m i c 2 ix + m i c 2 iy J iyz → I iyz -c iy c iz m i (
p 1 = m 1 . . . m n T ∈ R n , p 2 = c 1x c 1y c 1z . . . c nx c ny c nz T ∈ R 3n , p 3 = I T 1 . . . I T n T ∈ R 6n , (27)
where

I i = I ixx I ixy I ixz I iyy I iyz I izz T . ( 28
)
It is possible to provide lower bounds (LB) and upper bounds (UB) to p based on a priori information. In particular, condition ( 11) is managed by assigning a lower bound of zero to each link mass. Upper bounds for the masses are assigned exploiting, for instance, data retrieved from the datasheets of the robot. Moreover, for each link, one can easily infer that the center of mass is located inside the smallest parallel box which includes the link geometry, in the most general case. The lower and upper bounds are then set in order to guarantee a physical meaning to the obtained solution. In order to retrieve a possible set of dynamic parameter p, we propose to solve the nonlinear optimization problem depicted in Algorithm 1.

       min p k f (p k ) = φ(p k ) + ξ k γ(p k ) = π(p k ) -π 2 + ξ k ι g(hι(p k )) s.t. LB ≤ p k ≤ U B ξ k+1 ← 10k (29) 6 7 end
The first two lines of Algorithm 1 are the initialization step of the algorithm: the starting point is randomly selected between the lower and the upper bounds using a uniform distribution. Moreover, ξ 1 is set to zero. Lines 3 -6 of Algorithm 1 consist in solving the constrained nonlinear optimization problem κ times: at a given step k = 1 . . . κ, the initial state is the optimal solution found at step k -1. The objective function presents also an exterior penalty function, in which ξ k is the -progressively increasing -penalty coefficient: this function provides a penalty in case one of any additive constraint h ι (p k ) is violated; function g(•) is chosen to return a measure of the constraint violation. In particular, two kinds of constraints have been included:

• for each link i , the triangle inequality (16) must be satisfied; • the total sum of the link masses must be in a given range, that is:

m rob,min ≤ i m i ≤ m rob,max . (30)
Note that the presented framework is extremely flexible, and further external constraints can be easily included. The manifold generated by the cost function f (p k ) contains multiple local minima, and therefore a global optimization method, like genetic algorithms [START_REF] Mitchell | Machine Learning[END_REF] or simulated annealing (SA) [START_REF] Russell | Artificial Intelligence: A Modern Approach[END_REF] is mandatory to address the problem (29). In the present case, we have used SA, applying a more sophisticated interior-point (IP) Nelder-Mead local optimization algorithm at the end of each SA iteration k. Moreover, Q runs have been launched having a different random initial point p 0 , in order to span as much as possible the cost function manifold. The improvement of the parameters retrieval framework (with respect to the one presented in [START_REF] Gaz | Extracting feasible robot parameters from dynamic coefficients using nonlinear optimization methods[END_REF]) has been proven necessary since the introduction of some constraints -as for instance the triangle inequality of the inertia tensorseventually led the algorithm to get stuck in local minima. This was caused by recurring abrupt changes in the cost function given by the constant penalty function adopted before. Therefore, a violation-dependent penalty [START_REF] Bertsekas | Constrained Optimization and Lagrange Multiplier Methods[END_REF] has been implemented for the present algorithm, solving this problem. Moreover, the sequence of successive runs of the algorithm could in practice help in improving the solution. The term φ(p k ) of the parameters retrieval algorithm in eq. ( 29) requires, to be computed, a previous identification step returning the coefficients estimation π. Another possible φ(p k ) function, yielding to a single-step procedure, is described in the Supplementary material accompanying this letter.

VI. RESULTS

In this section, results from the dynamic coefficients and joint friction estimations for the Panda robot are reported, in addition to results from the parameters retrieval procedure.

In order to obtain a numerical estimation of the dynamic coefficients π s,R (see eq. ( 21)), the numerical values of the inertia matrix and the gravity vector are retrieved from a set of M = 1010 static positions, spanning the whole joint space, according to the robot documentation 3 . The symbolic vector s(q) (see eq. ( 17)) has been computed according to [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF], using the modified Denavit-Hartenberg convention.

Since our robot is mounted on a table parallel to the ground, the first element of the gravity vector (relative to joint 1) is not informative, and therefore it is discarded.

Stacking all the numerical quantities of the data acquisition phase and after evaluating the corresponding regressor (see eq. ( 20)), we obtained that rank(Y s ) = 43: using the Matlab function rref, which implements Gauss-Jordan elimination 3 The joint position bounds of the Franka Emika Panda robot can be found at: https://frankaemika.github.io/docs/control_parameters.html technique, we finally obtained the dynamic coefficients estimations πs,R according to eq. ( 21), from a regressor Y s,R whose condition number is 49. Moreover, the relative error percentage of predictions (defined in eq. ( 84) of [START_REF] Sousa | Physical feasibility of robot base inertial parameter identification: A linear matrix inequality approach[END_REF]) for the identification set is 0.031%. Computing their standard deviations (see [START_REF] Janot | A generic instrumental variable approach for industrial robot identification[END_REF], [START_REF] Sousa | Physical feasibility of robot base inertial parameter identification: A linear matrix inequality approach[END_REF]), we found that two coefficients exhibit a standard deviation greater than 30%, and therefore they were discarded (since their estimations are not reliable). All the estimated dynamic coefficients are reported with their standard deviations in the Supplementary material, together with a comparison with the corresponding coefficients obtained from the classical identification procedure (see eq. ( 8)).

From the identified dynamic coefficients πs,R , we were able to reconstruct the inertia matrix M (q), the gravity vector ĝ(q) and the Coriolis and centrifugal force vector ĉ(q, q) following a Lagrangian approach (see Sec. IV).

After deriving the inverse dynamic model ( 23), we validate it by comparing the measured joint torques with the estimated ones during several motions. In particular, the robot was commanded in velocity-mode by means of sinusoidal trajectories in the joint space. In other words, each joint is commanded according to the following equation:

qi,des (t) = A i sin 2 π T i t , i ∈ [1, . . . , n] (31)
where A i is the amplitude of the velocity profile and T i is the period of the sinusoidal signal for the i-th joint. The numerical values for A i and T i , i ∈ [1, . . . , n] for a typical experiment are reported in Table I, where 5458 samples were collected. The joint torque signals were recorded (and filtered through a 4-th order zero-phase digital Butterworth filter with a cutoff frequency of 1 Hz) during this motion, and compared with our Lagrangian inverse dynamics estimations τ (q, q, q) (from eq. ( 23)), feeding that model with the measured joint positions and velocities, and with the joint accelerations obtained by numerical differentiation of the filtered velocities. The joint torque comparison is reported in Fig. 3: the torque estimations are almost perfectly superimposing the measured ones for joints 1 . . . 4, while the last joints show some discrepancies. One can notice, though, that the difference between the two signals strongly depends on joint velocities (it is more evident, e.g., for joint 7): this behavior is typical of joint friction. Therefore, we performed the estimation of the joint friction according to the procedure reported in Sec. IV-B: we collected more than 10k samples of joint velocities and torques during sinusoidal motions in (31); eventually, we found that the best friction model that fits the data was given by a sigmoidal The red lines represent the measured joint torques during the validation experiment reported in Tab I, while the green lines are the torque estimations τ computed according to eq. ( 23)). The blue lines represent the joint torque estimates comprehensive of the joint friction term (25). The dashed and the dotted lines are the errors between the torque sensors readings and, respectively, the torques estimates without and with the friction component. function (25), which is characterized by 3 parameters per joint. These were estimated by solving a nonlinear least squares problem by means of a Nelder-Mead routine (Matlab function fminsearch), using as fitting data the differences ∆τ i between the measured joint torques and the estimated torques τi for each joint i separately. The results of the fitting procedures are reported in Fig. 4, while the numerical values are reported in the supplementary material. Finally, adding the newly estimated joint friction components τ f ( q) = τf,1 ( q1 ) τf,2 ( q2 ) . . . τf,7 ( q7 ) T to the previous estimations τ , we obtained a satisfactory compensation, as shown with the blue solid lines of Fig. 3.

The estimated inverse dynamic model in Lagrangian form described so far, though, is very cumbersome for a 7 DoF robot, and computationally intensive, such that it would not be reliable under real-time constraints. For this reason, a Newton-Euler (N-E) approach would be more appropriate and effective to quickly return joint torque estimates, due to its recursive form. Nevertheless, a N-E routine requires the dynamic parameters (masses, inertia tensors and centers of mass of each link) of the robot. Exploiting the parameters retrieval algorithm (29), though, we are able to extract a feasible set of dynamic parameters, which provides the same dynamics (although there is no guarantee that the estimated robot parameters set is coincident with the real one).

In order to implement Algorithm 1, the (bounded) simulated annealing Matlab function simulannealbnd has been used, together with the IP hybrid function fmincon; the problem parameters were Q = 100 (total number of independent runs, Fig. 4. Joint friction estimates. The red dots are the differences ∆τ j = τ j -τ j for each joint j, while the blue lines are the sigmoidal fitting functions. The effect of friction is clear for joints 4 . . . 7, while its contribution to the total joint torque is slight for the other joints.

providing Q solutions -possibly coincident, since there might exist multiple minima) and κ = 10 (number of successive runs). The penalty functions g(h ι (p k )) are chosen as the distance functions of the external constraints, that is:

where g 1,i regards the triangle inequality on each inertia tensor (eq. 16), while g 2 regards the constraint on the total mass of the robot (we chose m rob,min = 16 kg and m rob,max = 20 kg).

A feasible set of dynamic parameters p has been then retrieved and its numerical values are included in the supplementary material. In order to validate this set, we inserted the retrieved dynamic parameters in a N-E routine in order to compute the joint torques τ ne necessary to perform a given validation trajectory in the joint space. In particular, we commanded a sequence of sinusoidal trajectories (with different amplitudes and periods) to the joints according to eq. (31). We then compared the measured joint torques τ with the estimations τ ne + τ f : the result of this comparison is reported in Fig. 5, showing good results.

VII. CONCLUSIONS

In this letter, we addressed the problem of extracting a feasible set of parameters that characterizes the dynamics of the robot. We identified the dynamic coefficients through a standard least squares algorithm by means of a reverse engineering approach. Thanks to an improved version of the algorithm proposed in [START_REF] Gaz | Extracting feasible robot parameters from dynamic coefficients using nonlinear optimization methods[END_REF], we retrieved a set of feasible parameters by solving a nonlinear optimization problem, taking into account several constraints including the physical bounds on the dynamic parameters (such as the total mass of the robot) and the triangle inequalities of the link inertia tensors. The proposed framework was validated by deriving the Lagrangian model of the Franka Emika Panda robot and testing it through experiments. We also validated the extracted feasible set of dynamic parameters using a Newton-Euler routine. To the best of the authors' knowledge, this is the first paper that retrieves the dynamic coefficients and a feasible parameter set for the Panda, a robot that is being more and more used in the robotics community. A major feature of the proposed framework is that it can be easily modified to include further (possibly) nonlinear constraints (e.g., based on a priori information on the robot). We are releasing both the dynamic model and the parameters retrieval framework as open-source code, such that they will be available to the robotics community at the following website: http://diag.uniroma1.it/~gaz/panda2019.html.

As a future work, the authors would like to consider other possible nonlinear constraints and use the derived dynamic model for detecting, isolating and reacting to collisions between the robot and the environment.

Fig. 1 .

 1 Fig. 1. Denavit-Hartenberg frames and table of parameters for the Franka Emika Panda. The reference frames follow the modified Denavit-Hartenberg convention. In the figure, d 1 = 0.333 m, d 3 = 0.316 m, d 5 = 0.384 m, d f = 0.107 m, a 4 = 0.0825 m, a 5 = -0.0825 m, a 7 = 0.088 m.

1 p T 2 p T 3 T

 123 26) for each link i , i = 1 . . . n. We can now rearrange the parameters vector p = p T ∈ R 10n as:

Algorithm 1 : 4 / 5 /

 145 Parameters retrieval 1 p 0 ← LB + (U B -LB)u, with u ∼ U (0, 1); 2 ξ1 ← 0; 3 for k = 1, . . . , κ do / Start the optimization from the previous step solution p k,init ← p k-1 ; / Solve the following optimization problem

Fig. 3 .

 3 Fig. 3. Comparison between the torques of the derived E-L dynamic model.The red lines represent the measured joint torques during the validation experiment reported in Tab I, while the green lines are the torque estimations τ computed according to eq. (23)). The blue lines represent the joint torque estimates comprehensive of the joint friction term (25). The dashed and the dotted lines are the errors between the torque sensors readings and, respectively, the torques estimates without and with the friction component.

g 1 , 7 i=1 m i , 7 i=1m

 177 i = -min {0, tr(I i )/2 -λ max (I i )} , i ∈ [1, . . . 7] g 2 = -min 0, m rob,maxi -m rob,min

Fig. 5 .

 5 Fig. 5. Comparison between the torques τ ne (including friction τ f ) from the N-E dynamic model and the ones retrieved from the the joint torque sensors. The overlapping of two signals shows the quality of the obtained parameters estimation p.

  Dynamic Identification of the Franka Emika Panda Robot with Retrieval of Feasible Parameters Using Penalty-based Optimization Claudio Gaz 1 Marco Cognetti 2 Alexander Oliva 3 Paolo Robuffo Giordano 2 Alessandro De Luca 1

TABLE I AMPLITUDES

 I A i AND PERIODS T i , i ∈ [1, . .. , n] OF THE SINUSOIDAL TRAJECTORIES IN EQ. (31) USED FOR VALIDATING THE FRICTION ACTING

			ON THE ROBOT JOINTS.
		1	2	3	4	5	6	7
	A i	2.21 -2.21	1.2	-2.1	-2.3	2.1	-2.5
	T i	3.68	2.04	2.98 1.75 4.43 2.749 1.06

Dipartimento di Ingegneria Informatica, Automatica e Gestionale, Sapienza Università di Roma, via Ariosto

25, 00185 Roma, Italy. E-mail: {gaz, deluca}@diag.uniroma1.it.2 CNRS, Univ Rennes, Inria, IRISA, Rennes, France, E-mail: {marco.cognetti, prg}@irisa.fr.3 Inria, Univ Rennes, CNRS, IRISA, Rennes, France, E-mail: alexander.oliva@inria.fr.

The dynamic identification procedure has been performed in this paper according to a reverse engineering approach. In the Supplementary material accompanying this manuscript, however, the reader can find also the results of the classical identification approach.


Claudio Gaz

Marco Cognetti

Alexander Oliva

Paolo Robuffo Giordano

Alessandro De Luca

Dynamic Identification of the Franka Emika Panda Robot with Retrieval of Feasible Parameters Using Penalty-Based Optimization

come

Abstract-In this document, we report a number of supplementary materials to our RA-L paper [START_REF] Gaz | Dynamic identification of the Franka Emika Panda robot with retrieval of feasible parameters using penalty-based optimization[END_REF]. I provides the numerical values of the parameters obtained with the fitting procedure described in Sec. IV.B of [START_REF] Gaz | Dynamic identification of the Franka Emika Panda robot with retrieval of feasible parameters using penalty-based optimization[END_REF]. These parameters refer to the equations τ f,j = ϕ 1,j 1 + e -ϕ2,j ( qj +ϕ3,j ) -ϕ 1,j 1 + e -ϕ2,j ϕ3,j , j ∈ [1, . . . , 7], which model the friction acting at the robot joints. II. COMPARISON OF DYNAMIC COEFFICIENTS The manufacturer does not provide any information on the dynamic parameters of the Panda robot, and we have no ground truth values at our disposal for the dynamic coefficients (a.k.a., base parameters) that have to be estimated.

I. JOINT FRICTION PARAMETERS Table

Nevertheless, we have compared the dynamic coefficients obtained with the reverse engineering approach in [START_REF] Gaz | Dynamic identification of the Franka Emika Panda robot with retrieval of feasible parameters using penalty-based optimization[END_REF] with those of a classical identification method [START_REF] Gautier | Direct calculation of minimum set of inertial parameters of serial robots[END_REF] that uses measurements of the joint positions and joint torques only. In this case, static positioning is no longer sufficient and the robot must execute sufficiently exciting trajectories in order to allow a reliable estimation of the dynamic coefficients. In particular, following [START_REF] Sousa | Physical feasibility of robot base inertial parameter identification: A linear matrix inequality approach[END_REF], we designed the reference trajectories for the joints j ∈ {1, 7} as q j (t) = L l=1 a l,j lω f sin (lω f t) -b l,j lω f cos (lω f t) + q 0,j , where L = 5 is the number of harmonics, ω f = 0.15π, and the coefficients a l,j , b l,j and q 0,j are reported in Table II. These trajectories have been commanded to the robot joints during 21 s, collecting a total of 20544 samples of positions and torques per joint. Since we noticed already that joint friction is not negligible for the robot under study, we added in the robot dynamic model suitable friction functions to each joint j ∈ {1, 7}, which are linear in the three parameters f v,j , f c,j , and f o,j -see eq. ( 8) in [START_REF] Gaz | Dynamic identification of the Franka Emika Panda robot with retrieval of feasible parameters using penalty-based optimization[END_REF]:

τ f,j ( qj ) = f v,j qj + f c,j sign( qj ) + f o,j .
These functions may suffer from sudden discontinuities in the neighbourhood of qj = 0, j = 1, . . . , 7, but are adopted here because of their linearity in the defining parameters.

The results of the identification process is summarized in Table III. The first column provides the symbolic expressions of the dynamic coefficients. The second and third columns show the values of the coefficients πR-E identified using the reverse engineering approach and, respectively, their standard deviations %σ(π R-E ). The regressor condition number is 48.9, while the relative error percentage is 0.0315%. The fifth and sixth columns show the values of the coefficients πCLS identified through the classical method with exciting trajectories and, respectively, their standard deviations %σ(π CLS ). In this case, the regressor condition number is 121.88, while the relative error percentage is 1.41%. Finally, the fourth column reports the numerical value of the dynamic coefficients π(p) computed by substituting the values of the dynamic parameters retrieved with our framework [START_REF] Gaz | Dynamic identification of the Franka Emika Panda robot with retrieval of feasible parameters using penalty-based optimization[END_REF] in their symbolic expressions. It can be noted that these values are slightly different from the original ones, showing that the ordinary least squares solutions πR-E and πCLS are usually not physically consistent [START_REF] Sousa | Physical feasibility of robot base inertial parameter identification: A linear matrix inequality approach[END_REF]. On the other hand, in order to estimate a set of dynamic coefficients which are physically consistent, the following one-step approach can be used. It is possible, indeed, to change the term φ(p k ) of the cost function in the Parameters Retrieval algorithm (29) in [START_REF] Gaz | Dynamic identification of the Franka Emika Panda robot with retrieval of feasible parameters using penalty-based optimization[END_REF], as:

φ(p k ) = Y π(p k ) -τ 2 ,
where Y is the stacked regressor, τ is the stacked measurements vector and π(p k ) is the coefficients vector computed from the current parameters vector p k . The drawback of the use of this function φ (with respect to the one reported in [START_REF] Gaz | Dynamic identification of the Franka Emika Panda robot with retrieval of feasible parameters using penalty-based optimization[END_REF]) is that it is computed from the cumbersome stacked regressor, thus affecting the computational time. The identified dynamic coefficients πOSI obtained from the previously described one-step identification procedure are reported in the seventh column of Tab. III. Note that the fourth column of Table III has two missing entries for the eight and twenty-fourth dynamic coefficients: these small coefficients are in fact discarded due to their large standard deviations (see the corresponding second column in the table).

In conclusion, the similarity between the results obtained with the classical method and with our approach indirectly confirms the validity of the set of dynamic coefficients in Table III for the Panda robot.

III. COMPARISON OF DYNAMIC PARAMETERS

The dynamic parameters of the Panda robot, i.e., the mass, center of mass, and inertia tensor of each link, have been retrieved by solving the nonlinear optimization problem presented in Sec. V of [START_REF] Gaz | Dynamic identification of the Franka Emika Panda robot with retrieval of feasible parameters using penalty-based optimization[END_REF] as well as by using the LMI-SDP framework presented in [START_REF] Sousa | Physical feasibility of robot base inertial parameter identification: A linear matrix inequality approach[END_REF], based on the Python code available at https://github.com/cdsousa/wam7_dyn_ident. This code has been slightly modified in order to include the triangular inequalities on the inertia tensors, according to [START_REF]Inertia tensor properties in robot dynamics identification: A linear matrix inequality approach[END_REF]. The results are reported in the two Tables VIII andIX.

In Tab. VIII, the value of c 1z is reported as " * " (don't care) since it does not have any influence on the dynamics (i.e., it does not appear in the E-L model, and can take any value).

The used lower (LB) and upper (UB) bounds and the final value p obtained with our optimization algorithm are shown in the second to fourth columns. The final parameters δ (including 21 friction parameters) obtained with the LMI-SDP method are reported in the fifth columns. In particular, these latter values were retrieved from the dynamic coefficients β e by means of eq. ( 49) in [START_REF] Sousa | Physical feasibility of robot base inertial parameter identification: A linear matrix inequality approach[END_REF], using the same bounds (properly manipulated) that we adopted for our algorithm. The resulting δ set was then slightly manipulated, in order to obtain for each link the center of mass from the its first moment of inertia and the barycentric inertia tensor from the inertia tensors w.r.t. the link frame.

Figure 1 shows a comparison between the torques measured by the joint torque sensors of the Panda robot during a validation trajectory and the estimated torques generated from a Newton-Euler (N-E) routine using the p parameters and the δ parameters. In order to further validate the parameter sets p and δ, we also designed 10 new validation trajectories (each lasting for 10 seconds), performing the same joint torque comparisons and then computing the mean square error (MSE) for each trajectory. The concatenated joint trajectories are reported in Fig. 2, the concatenated joint torques in Fig. 3, and the total MSE in Table IV. These results confirm that both estimates are reliable and consistent. Fig. 1. Comparison between measured joint torques during a validation trajectory (blue solid lines) and estimated torques generated from a N-E routine using the p parameters from our parameter retrieval method in [START_REF] Gaz | Dynamic identification of the Franka Emika Panda robot with retrieval of feasible parameters using penalty-based optimization[END_REF] (solid red lines) and the parameters δ from the LMI-SDP framework in [START_REF] Sousa | Physical feasibility of robot base inertial parameter identification: A linear matrix inequality approach[END_REF] (solid yellow lines). The torque errors are reported as well.

TABLE IV MEAN SQUARE ERROR (MSE) [(NM) 2 ] OF TORQUE PREDICTIONS

USING THE PARAMETERS p OF OUR METHOD [START_REF] Gaz | Dynamic identification of the Franka Emika Panda robot with retrieval of feasible parameters using penalty-based optimization[END_REF] AND THE PARAMETERS δ OF LMI-SDP APPROACH [START_REF] Sousa | Physical feasibility of robot base inertial parameter identification: A linear matrix inequality approach[END_REF]

IV. VALIDATION IN V-REP

In order to validate the set of identified parameters reported in Tables VIII and IX, we built a robot model in V-REP [START_REF]V-REP virtual robot experimentation platform[END_REF], as shown in Figure 4. We adjusted the parameters for each link, in order to be compatible with the V-REP interface (e.g., for each link, the inertia tensor in V-REP has to be expressed w.r.t. the given CoM reference frame).

The validation process follows the idea that if the identification of the robot dynamic parameters is satisfactory, the measured torques on the real robot, apart from friction which is not simulated in V-REP, should match the ones provided Fig. 3. Joint torques recorded during the 10 validation trajectories reported in Fig. 2. by the physics engines of the simulator when performing the same trajectory.

In particular, V-REP provides four different free licence physics engines. A performance comparison of these engines is reported in Figs. 5 and6. Both the real and the simulated robot are commanded through sinusoidal joint velocity inputs (spanning their maximum available range). Figure 5 shows, for each joint, a filtered version of the measured torques of the real robot (solid blue lines) vs. those provided by the different physics engines (light blue: Bullet 2.78, orange: Bullet 2.83, yellow: ODE, and purple: Newton). As it can be seen, the simulated torques for all engines have a good overlap with the measured ones. Figure 6 shows the mean and the standard deviation of the torque errors for each joint. This error is defined as the difference between the measured torque on the real robot and the one generated from the simulation. We can notice that the mean torque errors are close to zero while the peaks of the standard deviations are due to measurement noise on the real robot torque readings and to numerical inaccuracies of the physics engines. In our experiments, all engines showed a good numerical performance, as confirmed by the values in Table V, and looking at the results in joint space, ODE (although it is the most noisy among the engines) and Bullet 2.83 seem to be the best performing engines in V-REP. In order to confirm the quality of the proposed solution, we analyzed the performance also in Cartesian space. To this aim, we compared the 3D pose (as obtained from V-REP) of the end-effector of the simulated robot with the corresponding pose obtained by reading the joint positions from the real robot and passing them to the direct kinematics module. We defined then a position and an orientation error (following the Roll-Pitch-Yaw convention), whose mean values are reported in Fig, 7 for all physics engines available in V-REP. This comparison is shown also in the video accompanying [START_REF] Gaz | Dynamic identification of the Franka Emika Panda robot with retrieval of feasible parameters using penalty-based optimization[END_REF].

The results support once again the quality of the estimated dynamic parameters, since the mean position error is below 1.5 cm and the mean orientation error is below 3 degrees for all the engines. Particularly good results are obtained with the ODE engine, with mean errors of less than 5 mm in position and below 1.2 deg in orientation. Moreover, the maximum errors for the worst performing physics engine (Bullet 2.78) are 9.13 cm in position and 3.14 deg in orientation.

V. USE OF OUR METHOD WITH NONLINEAR AND CONDITIONAL CONSTRAINTS

In order to retrieve a feasible set of dynamic parameters p for the Panda robot, only linear constraints were used, since the non-convex shapes of the links would allow their centers of mass to be located even outside the link themselves. We present here a simple example of retrieval of dynamic parameters where the use of nonlinear constraints would be preferable, and which fits then in the framework of our method [START_REF] Gaz | Dynamic identification of the Franka Emika Panda robot with retrieval of feasible parameters using penalty-based optimization[END_REF].

Consider the spatial 2R robot with orthogonal joint axes in Fig. 8. Its (standard) Denavit-Hartenberg kinematic parameters are reported in Table VI. The robot body consists of two links: link 1 is a cylinder of height l 1 = 1.5 m and radius r 1 = 0.15 m; link 2 has a total length of l 2 = 0.8 m and is composed by a cylinder of length l 2 -l 2 and radius r 2 = 0.08 m, connected at the end with a truncated cone of length l 2 = 0.2 m and with radius r 3 = 0.02 m for its smaller base.

i a i α i d i θ i 1 0 π/2 l 1 q 1 2 l 2 0 0 q 2
The assumed ground truth values of the dynamic parameters p real are given in the fourth column of Table VII.

The symbolic form of the dynamic coefficients π 2R ∈ R 8 of this robot is:

π 2R =           J1yy + J2yy -0.64m2 J2xx -J2yy + 0.64m2 J2xy J2xz -0.8c2zm2 J2yz J2zz -0.8m2 0.8m2 + c2xm2 c2ym2          
. In order to estimate the values of the dynamic coefficients π 2R , a numerical simulation is performed using p real and the joint torques are recorded during an exciting motion. Then the dynamic coefficients are properly estimated with an ordinary least squares method.

At this stage, the algorithm for the retrieval of the dynamic parameters p is launched, using the bounds on physical feasibility reported in Table VII, upper and lower bounds on the total mass, 1 ≤ m 1 + m 2 ≤ 15, and linear constraints on the inertia tensors (exploiting the triangular inequality). When using box constraints on the position of the center of mass of each of the two links, the solution would be searched in an extra volume (i.e., a parallelepiped) which is V link1 = 4/π 1.27 times larger than the volume of the cylindric link 1 and

V link2 = 4/π (1 + 2l 2 /l 2 )
2 times larger than the one of link 2. On the other hand, considering that the distance of the center of mass from the major link axis must be less than its radius, the previous (approximate)

c 1z -∞ ∞ -0.02 -3.8e-04 m c 2x -0.8 0 -0.1 -0.24487 m c 2y -∞ ∞ 0 6.16e-06 m c 2z -∞ ∞ 0 -3.66e-03 m I 1xx 0 1 0.03 0.61296 kg•m 2 I 1xy -1 1 0 0.071738 kg•m 2 I 1xz -1 1 0 -0.0099797 kg•m 2 I 1yy 0 1 0.03 0.022982 kg•m 2 I 1yz -1 1 0 0.031154 kg•m 2 I 1zz 0 1 0.01 0.6172 kg•m 2 I 2xx 0 1 0.01 0.019819 kg•m 2 I 2xy -1 1 0 -4.525e-06 kg•m 2 I 2xz -1 1 0 -0.0077412 kg•m 2 I 2yy 0 1 0.04 0.35386 kg•m 2 I 2yz -1 1 0 -5.9342e-05 kg•m 2 I 2zz 0 1 0.04 0.34412 kg•m 2
box constraints on the center of mass of the two links can be replaced by the following (exact) nonlinear and conditional (if-else) constraints g CoM1 and g CoM2 , to be used within our core optimization Algorithm 1 in [START_REF] Gaz | Dynamic identification of the Franka Emika Panda robot with retrieval of feasible parameters using penalty-based optimization[END_REF]:

α 1 = r 1 -c 2 1y + c 2 1z , α 2 =      r 2 -c 2 2y + c 2 2z , if c 2x < -l 2 r 2 + r3-r2 l 2 (l 2 + c 2x ) -c 2 2y + c 2 2z , otherwise, g CoM1 = -min{0, α 1 }, g CoM2 = -min{0, α 2 }.
Taking advantage of the knowledge of the link shapes, these nonlinear constraints ensure that each center of mass will lie inside the corresponding link shape but that no feasible position is being excluded, thus yielding a complete solution.

A validation test was finally performed, comparing the joint torques on sinusoidal trajectories that are computed using the real parameters with the joint torques estimated by means of a N-E routine fed with the parameters p reported in the fifth column of Table VII. The retrieved solution generates just the same dynamics as obtained with the real parameters.

Another interesting quantity that can be estimated only when using the N-E inverse dynamics algorithm is the wrench acting at the robot joint level, i.e., the exchanged forces/moments between two successive links connected by a joint. To this aim, it is very important to retrieve a set of dynamic parameters that is as close as possible to the real one. To show this, we relaxed the constraints in order to obtain a second solution set of dynamic parameters p , which, however, contains also non-feasible parameters (i.e., parameters that have no physical meaning). VII. A N-E routine fed with the feasible parameters p and with some unfeasible parameters p returns similar estimated torques (respectively, red and yellow lines) along the motion. When comparing these estimates with the actual torques recorded along a sinusoidal joint trajectory (blue lines, almost overlapped by the red line), we notice that both sets provide a good torque estimation.

In Fig. 9, the yellow lines represent the estimated joint torques during a validation experiment coming from a N-E routine that is being fed with the unfeasible parameters p . Despite of this, reliable motion torque estimations can still be appreciated. On the contrary, a strong inconsistency if feasible (i.e., physically consistent) dynamic parameters are provided to the routine, while estimates obtained from an unfeasible set p can be even 10 6 times larger (!) than real internal forces (note that the yellow line was scaled for the ease of reading).

may be observed when estimating the internal forces acting on the joints, obtained as a byproduct of the standard N-E algorithm. As shown by the values of these estimated forces in Fig. 10, it can be seen that estimations obtained from the feasible parameters p (red lines) are close to the real forces (blue lines), while estimations retrieved from the unfeasible parameters p (yellow lines) are almost 10 6 times larger than the real values "felt" at the joints. This confirms that a physically-consistent solution which is close to the real one is mandatory for estimating the joint wrenches.

Fig. 2 .

 2 Fig.2. Joint positions recorded during 10 validation trajectories (10 seconds each). The associated joint torques are reported in Fig.3.

Fig. 4 .

 4 Fig. 4. The Franka Emika Panda robot in V-REP.

Fig. 5 .

 5 Fig. 5. Comparison between the torques generated by the real robot during an experiment (blue) and those obtained from the simulated Panda robot in V-REP using different physics engines (Bullet 2.78, Bullet 2.83, ODE, Newton).

Fig. 7 .

 7 Fig. 6. Torque errors (w.r.t. the real robot) of the physics engines Bullet 2.78, Bullet 2.83, ODE, and Newton in V-REP during an experiment with the Panda robot for each of its 7 joints. The boxes indicate mean values of the errors while the bars denote their standard deviations.

Fig. 8 .

 8 Fig. 8. A spatial 2R robot with joints along axes z 0 and z 1 .

Fig. 9 .

 9 Fig. 9. Validation of estimated parameters in TableVII. A N-E routine fed with the feasible parameters p and with some unfeasible parameters p returns similar estimated torques (respectively, red and yellow lines) along the motion. When comparing these estimates with the actual torques recorded along a sinusoidal joint trajectory (blue lines, almost overlapped by the red line), we notice that both sets provide a good torque estimation.

  with feasible pars.( 10 6 ) est. force with unfeasible pars.

Fig. 10 .

 10 Fig.10. Internal forces acting on the joints of the simulated 2R robot during a sinusoidal joint trajectory. Estimates of the internal forces obtained from a N-E routine are significantly closer to the real internal forces (blue lines) if feasible (i.e., physically consistent) dynamic parameters are provided to the routine, while estimates obtained from an unfeasible set p can be even 10 6 times larger (!) than real internal forces (note that the yellow line was scaled for the ease of reading).

  Dynamic Identification of the Franka Emika Panda Robot with Retrieval of Feasible Parameters Using Penalty-Based Optimization -Supplementary material -Claudio Gaz 1 Marco Cognetti 2 Alexander Oliva 3 Paolo Robuffo Giordano 2 Alessandro De Luca 1

TABLE I ESTIMATED

 I JOINT FRICTION PARAMETERS.

		Joint 1	Joint 2	Joint 3	Joint 4	Joint 5	Joint 6	Joint 7	units
	ϕ1	5.4615e-01	8.7224e+03	6.4068e-01	1.2794e+00	8.3904e-01	3.0301e-01	5.6489e-01	N•m
	ϕ2	5.1181e+00 9.0657e+00	1.0136e+01	5.5903e+00	8.3469e+00	1.7133e+01	1.0336e+01	s/rad
	ϕ3	3.9533e-02	2.5882e-02	-4.6070e-02	3.6194e-02	2.6226e-02	-2.1047e-02	3.5526e-03	rad/s

TABLE III DYNAMIC

 III COEFFICIENTS IDENTIFIED USING OUR REVERSE ENGINEERING APPROACH ( πR-E ) AND USING A CLASSICAL APPROACH ( πCLS ), WITH THEIR STANDARD DEVIATIONS (IN PERCENTAGE). THE COEFFICIENTS π(p) ARE COMPUTED FROM THE RETRIEVED DYNAMIC PARAMETERS p. DYNAMIC COEFFICIENTS πOSI ARE OBTAINED FROM A ONE-STEP IDENTIFICATION PROCEDURE.

	Dynamic coefficients	πR-E	%σ(πR-E)	π(p)	πCLS	%σ(πCLS)	πOSI
	J2yy + J1zz	0.0292	0.3735	0.0373	0.0204	2.3609	0.0132
	J2xx -J2yy +J3yy +0.0998m3 +0.632c3zm3 +0.1067(m3 +m4 +m5 +m6 +m7)	0.9818	0.0176	0.9896	0.9188	0.1202	0.9485
	J2xy	-0.0052	1.5625	-0.0040	-0.0033	11.0392	0.0014
	J2xz	0.0284	0.2698	0.0103	0.0042	14.3570	7.4640 × 10 -4
	J2yz	-0.0035	2.3330	0.0008	-0.0155	2.3054	-0.0113
	J3yy + J2zz + 0.09985 m3 + 0.632 c3z m3 + 0.1067 (m4 + m5 + m6 + m7)	1.0428	0.0091	1.0358	0.9346	0.0973	0.9491
	J3xx -J3yy + J4yy -0.0068 m4	0.0106	1.2664	0.0116	0.0149	5.8149	-0.0164
	J3xy + 0.0825 c4z m4	-0.0002	30.3622	----	0.0127	3.5669	-0.0020
	J3xz	-0.0104	0.4874	-0.0055	-0.0061	5.2129	0.0021
	J3yz	-0.0048	1.2379	-0.0044	-0.0138	2.3907	-0.0055
	J4yy + J3zz + 0.0068 m4 + 0.0136 m5 + 0.0136 m6 + 0.0136 m7	0.1169	0.0691	0.1245	0.0825	0.5394	0.0922
	J4xx -J4yy + J5yy + 0.768 c5z m5 + 0.1406 (m5 + m6 + m7)	0.5324	0.0202	0.5408	0.5462	0.0951	0.5381
	J4xy + 0.0317 m5 + 0.0317 m6 + 0.0317 m7 + 0.0825 c5z m5	0.1501	0.0274	0.1388	0.1202	0.1665	0.1255
	Jxz4	0.0048	0.9079	0.0039	0.0040	6.5285	-0.0020
	J4yz	-0.0027	1.6177	-0.0016	0.0078	3.4723	1.8047 × 10 -4
	J5yy + J4zz + 0.1543 m5 + 0.1543 m6 + 0.1543 m7 + 0.768 c5z m5	0.6369	0.0103	0.6325	0.6196	0.0829	0.6089
	J5xx -J5yy + J6yy + 0.0077 m7	0.0322	0.2739	-0.0075	0.0114	3.9131	0.0137
	J5xy	-0.0037	0.9578	-0.0015	0.0036	5.7718	2.3992 × 10 -4
	J5xz	-0.0061	0.4932	-0.0046	-0.0102	1.6764	-0.0117
	J5yz	0.0077	0.4345	0.0022	0.0080	2.7133	0.0077
	J6yy + J5zz + 0.0077 m7	0.0181	0.3169	0.0276	0.0216	1.6395	0.0210
	J6xx -J6yy + J7yy -0.0077 m7	-0.0065	1.0591	-0.0009	-0.0157	1.8290	-0.0097
	J6xy + 0.088 c7z m7	0.0054	0.5331	0.0055	0.0014	8.9602	0.0032
	J6xz	-1.2901 × 10 -5	218.6430	----	0.0047	3.8644	0.0039
	J6yz	0.0008	3.3084	-0.0001	0.0006	27.1286	-1.0794 × 10 -4
	J7yy + J6zz + 0.0077 m7	0.0250	0.1773	0.0304	0.0145	1.8227	0.0204
	J7xx -J7yy	-0.0018	2.7378	0.0024	0.0024	7.3234	2.6087 × 10 -4
	J7xy	0.0012	2.0320	-0.0004	0.0049	1.8700	0.0011
	J7xz	0.0018	1.2299	-0.0017	0.0023	4.8889	0.0013
	J7yz	-0.0006	3.6756	-0.0005	-0.0027	4.7174	-0.0016
	J7zz	0.0013	2.6129	0.0049	-0.0043	4.4680	0.0022
	c2x m2	-0.0054	0.1551	-0.0020	0.0288	1.8926	0.0280
	c2y m2 -0.316 m4 -0.316 m5 -0.316 m6 -0.316 m7 -0.316 m3 -c3z m3	-3.1026	0.0003	-3.1043	-3.1363	0.0105	-3.1353
	0.0825 m4 + 0.0825 m5 + 0.0825 m6 + 0.0825 m7 + c3x m3	0.6874	0.0014	0.6842	0.6803	0.0481	0.6814
	c3y m3 -c4z m4	0.0234	0.0353	0.0282	-0.0121	1.2201	-0.0103
	c4x m4 -0.0825 m6 -0.0825 m7 -0.0825 m5	-0.4888	0.0016	-0.4901	-0.4726	0.0368	-0.4742
	0.384 m5 + 0.384 m6 + 0.384 m7 + c4y m4 + c5z m5	1.7185	0.0005	1.7207	1.7517	0.0193	1.7525
	c5x m5	-0.0100	0.0644	-0.0147	-0.0005	26.5699	0.0015
	c5y m5 -c6z m6	0.0771	0.0074	0.0824	0.0817	0.1362	0.0815
	0.088 m7 + c6x m6	0.1657	0.0034	0.1650	0.1637	0.0744	0.1639
	c6y m6 -c7z m7	-0.0677	0.0083	-0.0688	-0.0623	0.1527	-0.0613
	c7x m7	0.0062	0.0744	0.0077	0.0013	6.5781	6.3142 × 10 -4
	c7y m7	-0.0005	0.9490	-0.0031	0.0048	1.6783	0.0033
	fv,1				0.0665	1.9658	0.0628
	fv,2				0.1987	0.9615	0.2088
	fv,3				0.0399	3.4971	0.0361
	fv,4				0.2257	0.9124	0.2174
	fv,5				0.1023	1.5868	0.1021
	fv,6				-0.0132	12.2607	1.6128 × 10 -4
	fv,7				0.0638	2.1388	0.0632
	fc,1				0.2450	0.5889	0.2549
	fc,2				0.1523	0.9602	0.1413
	fc,3				0.1827	0.7324	0.1879
	fc,4				0.3591	0.4152	0.3625
	fc,5				0.2669	0.4668	0.2728
	fc,6				0.1658	0.8077	0.1529
	fc,7				0.2109	0.7055	0.2097
	fo,1				-0.1073	0.6932	-0.1069
	fo,2				-0.1566	2.9549	-0.1601
	fo,3				-0.0686	1.1566	-0.0718
	fo,4				-0.2522	0.9563	-0.2562
	fo,5				0.0045	19.4633	0.0079
	fo,6				0.0910	0.8390	0.0935
	fo,7				-0.0127	6.2628	-0.0070

  IN VALIDATION EXPERIMENTS OF FIG. 2.

		Joint 1 Joint 2 Joint 3 Joint 4 Joint 5 Joint 6 Joint 7
	MSE p	0.4068	1.0282	0.5951	0.4231	0.1749	0.0403	0.0267
	MSE δ	0.4955	1.0155	0.2751	0.6627	0.0828	0.0508	0.0284

TABLE VI DH

 VI TABLE OF THE SPATIAL 2R ROBOT.

TABLE VII LOWER

 VII AND UPPER BOUNDS FOR THE DYNAMIC PARAMETERS OF THE 2R ROBOT IN FIG. 8, WITH REAL p real AND ESTIMATED p VALUES.

	parameter	LB	UB	p real	p	units
	m 1	0	10	8	4.5788	kg
	m 2	0	5	3	3.7828	kg
	c 1x	-∞	∞	0.01	-0.016124	m
	c 1y	-1.42	0.08	-0.75	-0.74408	m

TABLE IX LOWER

 IX AND UPPER BOUNDS FOR THE INERTIA TENSOR ELEMENTS OF THE PANDA ROBOT AND THEIR CORRESPONDING ESTIMATED VALUES USING OUR APPROACH (p) AND THE ONE IN [3] (δ).

							parameter LB UB	p	δ (from [3])	units
							I 1xx	0	1	7.0337e-01	0.6	kg•m 2
							I 1xy	-1	1	-1.3900e-04	-3.2497e-14	kg•m 2
				TABLE VIII			I 1xz	-1	1	6.7720e-03	-6.4280e-09	kg•m 2
	LOWER AND UPPER BOUNDS FOR MASSES AND CENTERS OF MASS OF	I 1yy	0	1	7.0661e-01	0.6	kg•m 2
	THE PANDA ROBOT AND THEIR CORRESPONDING ESTIMATED VALUES	I 1yz	-1	1	1.9169e-02	-6.4280e-09	kg•m 2
	USING OUR APPROACH (p) AND THE ONE IN [3] (δ).		I 1zz	0	1	9.1170e-03	2.0353e-06	kg•m 2
	parameter	LB	UB	p	δ (from [3]) units	I 2xx	0	1	7.9620e-03	0.028	kg•m 2
	m 1	0	10	4.970684	2.0643e-05	kg	I 2xy	-1	1	-3.9250e-03	5.5650e-05	kg•m 2
	m 2	0	10	0.646926	9.2017	kg	I 2xz	-1	1	1.0254e-02	0.0056	kg•m 2
	m 3	0	10	3.228604	1.4575	kg	I 2yy	0	1	2.8110e-02	0.0291	kg•m 2
	m 4	0	10	3.587895	4.5856	kg	I 2yz	-1	1	7.0400e-04	-2.7636e-04	kg•m 2
	m 5	0	10	1.225946	0.5231	kg	I 2zz	0	1	2.5995e-02	0.0011	kg•m 2
	m 6	0	10	1.666555	2.0604	kg	I 3xx	0	1	3.7242e-02	2.0154e-06	kg•m 2
	m 7	0	10	7.35522e-01	0.1718	kg	I 3xy	-1	1	-4.7610e-03	-5.9144e-09	kg•m 2
	c 1x	-0.05 0.05	3.875e-03	5.202e-08	m	I 3xz	-1	1	-1.1396e-02	1.1523e-08	kg•m 2
	c 1y	-0.05 0.05	2.081e-03	5.202e-08	m	I 3yy	0	1	3.6155e-02	2.0243e-06	kg•m 2
	c 1z	-0.4	0.05	*	-0.1750	m	I 3yz	-1	1	-1.2805e-02	6.0956e-09	kg•m 2
	c 2x	-0.05 0.05	-3.141e-03	0.0015	m	I 3zz	0	1	1.0830e-02	2.0164e-06	kg•m 2
	c 2y	-0.15 0.05	-2.872e-02	-0.0578	m	I 4xx	0	1	2.5853e-02	2.0242e-06	kg•m 2
	c 2z	-0.05 0.05	3.495e-03	-0.0384	m	I 4xy	-1	1	7.7960e-03	7.1405e-09	kg•m 2
	c 3x	-0.05 0.15	2.7518e-02	0.1096	m	I 4xz	-1	1	-1.3320e-03	8.8917e-10	kg•m 2
	c 3y	-0.05 0.05	3.9252e-02	0.05	m	I 4yy	0	1	1.9552e-02	2.0036e-06	kg•m 2
	c 3z	-0.1	0.05	-6.6502e-02	-0.1	m	I 4yz	-1	1	8.6410e-03	-2.2039e-09	kg•m 2
	c 4x	-0.15 0.05	-5.317e-02	-0.0485	m	I 4zz	0	1	2.8323e-02	2.0260e-06	kg•m 2
	c 4y	-0.05 0.15	1.04419e-01	0.15	m	I 5xx	0	1	3.5549e-02	2.0057e-06	kg•m 2
	c 4z	-0.05 0.05	2.7454e-02	0.0147	m	I 5xy	-1	1	-2.1170e-03	1.6453e-10	kg•m 2
	c 5x	-0.05 0.05	-1.1953e-02	-0.004	m	I 5xz	-1	1	-4.0370e-03	-3.6824e-10	kg•m 2
	c 5y	-0.05 0.05	4.1065e-02	0.045	m	I 5yy	0	1	2.9474e-02	2.0037e-06	kg•m 2
	c 5z c 6x c 6y c 6z c 7x c 7y c 7z	-0.05 0.05 -0.05 0.15 -0.05 0.05 -0.05 0.05 -0.05 0.05 -0.05 0.05 0.04 0.15	-3.8437e-02 6.0149e-02 -1.4117e-02 -1.0517e-02 1.0517e-02 -4.252e-03 6.1597e-02	-0.05 0.0732 -0.0251 -0.0276 -0.0326 0.0087 0.0396	m m m m m m m	I 5yz I 5zz I 6xx I 6xy I 6xz I 6yy I 6yz	-1 0 0 -1 -1 0 -1	1 1 1 1 1 1 1	2.2900e-04 8.6270e-03 1.9640e-03 1.0900e-04 -1.1580e-03 4.3540e-03 3.4100e-04	2.3080e-09 2.0028e-06 2.0026e-06 1.5823e-09 2.0530e-09 2.0070e-06 -4.9387e-10	kg•m 2 kg•m 2 kg•m 2 kg•m 2 kg•m 2 kg•m 2 kg•m 2
							I 6zz	0	1	5.4330e-03	2.0073e-06	kg•m 2
							I 7xx	0	1	1.2516e-02	2.0026e-06	kg•m 2
							I 7xy	-1	1	-4.2800e-04	2.6159e-10	kg•m 2
							I 7xz	-1	1	-1.1960e-03	1.4055e-09	kg•m 2
							I 7yy	0	1	1.0027e-02	2.0037e-06	kg•m 2
							I 7yz	-1	1	-7.4100e-04	-4.2047e-10	kg•m 2
							I 7zz	0	1	4.8150e-03	2.0020e-06	kg•m 2

Dipartimento di Ingegneria Informatica, Automatica e Gestionale, Sapienza Universitá di Roma, via Ariosto

25, 00185 Roma, Italy. E-mail: {gaz,deluca}@diag.uniroma1.it.2 CNRS, Univ Rennes, Inria, IRISA, Rennes, France, E-mail: {marco.cognetti,prg}@irisa.fr.3 Inria, Univ Rennes, CNRS, IRISA, Rennes, France, E-mail: alexander.oliva@inria.fr.