Dynamic Identification of the KUKA LBR iiwa Robot With Retrieval of Physical Parameters Using Global Optimization

This paper focuses on the problem of extracting the physical dynamic parameters which are fundamental for computing the positive-definite link mass matrix. To solve this problem, a minimal set of dynamic parameters were firstly identified by the standard least squares method. In order to simplify the dynamics model, a new set of essential dynamic parameters were calculated by eliminating the poorly identified parameters with an iterative approach. Based on these dynamic parameters with better identification quality, a universally global optimization framework was proposed here to retrieve the set of physical dynamic parameters of a serial robot, in which parameter bounds, linear and nonlinear constraints with physical consistency can be easily considered, such as the triangle inequality of the link inertia tensors, the total link mass limitations, other user-defined constraints and so on. Finally, validation experiments were conducted on the KUKA LBR iiwa 14 R820 robot. The results show that the proposed optimization framework is effective, and the identified dynamic parameters can predict the robot joint torques accurately for the validation trajectories.


I. INTRODUCTION
Accurate dynamic parameters [1] of a robotic manipulator is important in many robotic applications [2]. Dynamic parameters are critical in the design of control laws based on dynamics model [3], in simulating the robot motion in some software [4] or in implementing some human-robot interaction algorithms, such as collision detection and reaction control [5], impedance control [6] and so on.
There are three main methods for obtaining the dynamic parameters of a robotic manipulator: physical experiments, computer aided design (CAD) techniques and dynamic The associate editor coordinating the review of this manuscript and approving it for publication was Christopher Kitts . parameter identification [7]. In physical experiment method, each link of the robot needs to be isolated and the dynamic parameters are obtained by using special measurement devices. For example, the mass of links can be weighted directly, the position vector of the center of mass can be estimated by determining counterbalanced points of the link and the diagonal elements of the inertia tensor can be obtained by pendular motions [1], [7]. However, the accuracy of this method depends on the measurement devices and the procedure is very tedious; hence, it should be conducted by the manufacturer before assembling the robot [8]. In the CAD method, the dynamic parameters of a robot can be obtained by analysing the 3-dimensional (3D) model in the CAD software. However, 3D models are usually not provided by the manufacturer and the parameters estimated with this method are inaccurate [7]. In addition, the friction parameters cannot be estimated by the first two methods. In the dynamic parameter identification method, the dynamic parameters are estimated by minimizing the errors between a function of the robot variables and an identification model. This method has been extensively used [9]. Compared to the above two methods, dynamic parameter identification can obtain more accurate dynamic parameter estimates. In general, an identification procedure ( Figure 1) usually consists of modeling, experiment design, data acquisition and signal processing, parameters estimation and model validation [10]. This paper presents a global nonlinear optimization framework for identifying the physical parameters of a KUKA LBR iiwa 14 R820, which can be generalized to any serial manipulator. In the nonlinear problem, the parameters bounds and the physical constraints can be easily added. To reduce the error accumulation, the objective function was calculated by the errors between the measured joint torque and the estimated torque. In addition, instead of using the base parameters to calculate the objective function, the essential parameters with better identified quality were used.
Atkeson et al. [11] first derived the equations of linearizing the inverse dynamics model of a robotic arm and identified its minimal set of inertia parameters, the so-called base parameters which are a set of independent identifiable parameters, by the standard least squares method. However, the universal derivation for analytical form of the base parameters set was not given in that work. To solve this problem, Gautier and Khalil [12] and Khalil and Bennis [13] proposed a direct method for determining and identifying the base parameters of serial manipulator, namely, regrouping parameters by means of closed-form functions of the geometric parameters of the robot. In these two works above, the applied identification methods are all based on standard least squares method, which is sensitive to the singularity of observation matrix. To avoid the accumulation of estimation errors and allow the calculation of the confidence intervals, the weighted least squares method, adding a weighted matrix multiplier, was proposed in [14]. However, a disadvantage of both the standard least squares method and the weighted least squares method is sensitivity to measurement noise. To overcome this problem, Swevers et al. [15] proposed an optimal excitation trajectories method of using finite Fourier series. In addition, the maximum-likelihood estimation method [16] was also formulated because of its asymptotically unbiased and efficient property.
Based on these previously excellent works, robotics researchers nowadays continue to explore and expand the methods of dynamic identification of the robot. A reverse engineering approach for identifying the dynamics model of KUKA LBR iiwa robot was presented in [17]. It can calculate the dynamic parameters inversely by using the known mass matrix and gravity vector information in KUKA Fast Research Interface (FRI) [18]. However, this method is obviously not a universal approach. For instance, not all robot companies can provide the mass matrix and gravity vector information to users in advance. Recently, a lot of parameter identification methods for the robot dynamics have also been proposed, such as the convex programming approach [19], adaptive control algorithm [20]- [22], extended Kalman filter method [23], neural networks method [24]- [26] and so on. However, the common problem for these methods are that the identification accuracy can not be guaranteed. To further simplify the dynamics model of the robot and improve the quality of identified parameters, a new set of dynamic parameters, namely essential parameters, was identified by eliminating the poorly identified base parameters in an iterative way [27], [28]. The identification of the base parameters and essential parameters are usually sufficient for many robotic applications, however, the retrieval of a set of feasible dynamic parameters, physical parameters [29], is also needed. This is the case, for example, in conducting the feedback linearization control laws under hard real-time constraints by the Newton-Euler method [30] or calculating the mass matrix of links used for robot collision detection. In this case, the symmetrical and positive-definite property of the mass matrix must be satisfied. However, the mass matrix calculated by the identified base parameters and essential parameters can only satisfy the symmetrical property but the positive-definite property can not be guaranteed. Thanks to the physical parameters, the positive-definite property of the mass matrix can be ensured by using them.
To obtain the physical parameters, the physical consistency of the identified parameters was considered in [31]. Then, the work [32] considered a nonlinear optimization problem to guarantee the physical feasibility of the identified parameters. However, the constraints considered in that work, such as only considering the positive link mass and the positive diagonal elements of the inertia tensor, were too simple. To solve this problem, some researchers proposed the framework of linear matrix inequalities to obtain the physically consistent set of parameters by semi-definite programming algorithm [33]. Afterwards, this method was enhanced by considering the triangle inequality of the inertia tensors [34], [35], which was originally presented in [36]. However, these methods require expressing constraints as linear matrix inequalities [30]. Recently, a local optimization identification framework has been proposed to identify the physical parameters of the KUKA iiwa robot dynamics in [37] by considering nonlinear constraints. The base parameters and the physical parameters of the robot were eventually identified and given. But the drawback of this method is that the optimal solution was based on a local optimization. In addition, the objective function of the nonlinear optimization was computed based on the base parameters, not the essential parameters which are with better identification quality, in that work. Gaz et al. [30] proposed a global nonlinear optimization framework, based on the inverse engineering method presented in their previous work [38], to identify the physical parameters of the iiwa robot dynamics. However, the objective function described in that paper was the errors between the identified based parameters and the one computed by the physical parameters. Compared with the objective function calculated by the errors between the measured joint torque and the torque estimated by the base parameters, the method in [30] will introduce the accumulation errors because the variables for calculating the errors were all obtained by estimation. In addition, no previous work can provide the identified results of all the three kinds of dynamic parameters simultaneously for a robot, namely the base parameters, essential parameters, and the physical parameters.
The main contributions of this work are described as follows. (1) A universally global optimization framework for identifying the physical parameters of serial manipulators is presented, that account for the parameters bounds, linear and nonlinear constraints. (2) The essential parameters with better identification quality were first used to calculate the objective errors instead of the base parameters. (3) The base parameters, essential parameters and the physical parameters of KUKA LBR iiwa 14 R820 robot were simultaneously obtained by the method presented in this paper.
The paper is organized as follows. The identification of essential parameters was shown in Section II. Section III describes the calculation of the physical parameters. In Section IV, the experiment results were given. Finally, discussion and conclusion were given in Section V.

A. INVERSE DYNAMICS IDENTIFICATION MODEL
By considering the rigid link dynamics with Lagrange method or Newton-Euler algorithm [39], the iiwa robot motion equations can be described as (1) where q,q,q ∈ n×1 , with n is the degrees of freedom of the robot, are the vectors of joint position, velocities and acceleration, respectively. M (q) ∈ n×n is the symmetric and positive-definite inertia matrix of links, C(q,q) ∈ n×n is the centrifugal and Coriolis matrix of links, and g(q) ∈ n×1 is the gravity vector of links. F v and F c ∈ n×n are the diagonal matrices of the viscous and Coulomb friction parameters. τ J ∈ n×1 is the joint torque readings of the torque sensors built in the robot.
To derive the linear model below of dynamics equation for dynamic parameters identification, the modified Newton-Euler method [1] is used in our case.
where Y (q,q,q) ∈ n×12n is the observation matrix which depends only on the motion data, π = [π T 1 , . . . , π T n ] T ∈ 12n×1 is the dynamic parameters of links, where for each link i, i = 1, . . . , n, the parameter vector π i is defined as Because some inertial parameters are completely unidentifiable and some others can only be identified in linear combinations, the set of parameters π to be identified can be reduced to a minimal set of parameters, which are referred to as base parameters π b ∈ b×1 , b is the base parameters number and it is 57 in this paper. The regrouping relationships between π and π b can be found in [12]. Therefore, (2) can be modified to where Y b (q,q,q) ∈ n×b is the observation matrix with full rank after deleting the linearly correlated columns in Y (q,q,q). In order to identify the dynamic model, the measured data of joint positions q and joint torques τ J are needed. The joint velocitiesq and accelerationsq are calculated from the measured joint positions q. After K data samplings of the robot during the experiments, one can apply the standard least-squares method to solve the base parameters by . . .
where Y b ∈ Kn×b and τ J ∈ Kn×1 are the stacked matrix of Y b and the stacked vector of measured joint torque vectors τ J , respectively. In order to reduce high frequency noise, the measured data of joint position and joint torque is filtered. Details about the data acquisition and signal processing can be found in Section II-C and Section IV-B. In order to understand the identified quality of the individual parameter in (6), an indicator, which is called the relative standard deviations (RSD) [37], is adopted. The calculation process of the RSD for the base parameters follows the equations below.
where σ π j,r% is the RSD of the j-th identified parameter in π b , σ π j is the standard deviation of the estimation error of the parameter j, j = 1, . . . , b. C π is the variance-covariance matrix of the estimation error, and σ 2 ρ is the variance of the parameter estimation.
If the RSD of an identified parameter is very large and its value is close to zero, then this parameter can be considered as a poorly identifiable parameter because its contribution to the dynamics model is negligible. Therefore, this parameter can be removed in order to simplify the dynamics model and a new set of parameters can be defined, which is called as essential parameters denoted by π e . The essential parameters are calculated using an iterative procedure starting from the base parameters estimation. At each iteration the base parameter which has the maximum RSD is cancelled. Then, a new set of parameters with new RSD is repetitively identified. The procedure ends when all parameters with absolute values less than 0.01 and corresponding RSD values larger than 40% are removed. It should be noted that the values 0.01 and 40% are chosen that work best for this problem. The readers should consider the trade-off between torque prediction accuracy by essential parameters and better identification quality of essential parameters for their own case to determine the values. The identified essential parameters will be used to calculate the objective function of the proposed nonlinear optimization problem of extracting the physical parameters in this paper. More details and results about the identified essential parameters are presented in Section IV.

B. OPTIMAL EXCITATION TRAJECTORIES
In order to improve the accuracy of the least-squares solution of (6), an optimal excitation trajectory is needed to obtain a well-conditioned observation matrix. In this paper, the classical parameterization Fourier series [15] are selected as optimal excitation trajectories and the optimal criterion to be minimized is the condition number of the observation matrix. The optimal excitation trajectory equation and the criterion rule with its constraints are shown as, min where L is the number of sine and cosine terms, a i,l and b i,l are optimized trajectory coefficients, the fundamental frequency is ω f 2π . The function cond(.) is used to solve the condition number of observation matrix Y b . The first constraint in (10) is to ensure the zero initial joint positions, velocities and accelerations of optimal excitation trajectories. The second one and the third one are to make the optimal trajectories satisfy the joint positions limits and joint velocities limits of KUKA LBR iiwa robot. More details and results about the optimal trajectories are shown in Section IV.

C. DATA ACQUISITION AND SIGNAL PROCESSING
The joint position and joint torque data are collected from the sensors in iiwa robot when the optimal excitation trajectories are run. Unfortunately, the joint torque data contains measurement noise. In addition, it is necessary to select an appropriate method to process the joint velocity and acceleration estimates used for dynamic parameters identification based on the measured joint position. Therefore, it is necessary to conduct signal processing to clean up the measured data, improve the signal-to-noise ratio of joint torque and joint position measurements and have an exact calculation of joint velocity and acceleration estimates.
Since the measured data are periodic, the signal-to-noise ratio can be improved by data averaging without using a lowpass noise filter. Averaging filter improves the quality of the data with the square root of the number of measured periods. The averaging of a signal x consisting of M periods of K samples is given by where x m (k) indicates the kth sample within the mth period and x(k) denotes the average of x. For robot dynamic identification, the signal x corresponds to the measured joint torque signal and joint position signal. Calculation of the observation matrix requires the estimates of the joint velocityq and joint accelerationq. Although the common numerical differentiation method and a low-pass filter can help to obtain the joint velocity and joint acceleration with less noise in off-line case, however, to obtain good filtering results, the filter coefficients of the low-pass filter must be tuned very carefully. The tuning process is time-consuming and the accuracy of estimated joint velocity and acceleration would be reduced if the filter coefficients were not selected well. To avoid this, VOLUME 8, 2020 a frequency-domain differentiation method originally proposed in [10] is introduced here to estimate the joint velocity and acceleration. Firstly, the averaged joint position data q(t) are transformed to the frequency domain data Q(k) by using the discrete Fourier transform (DFT). Next, a rectangular window multiplier is applied, where the spectrum is set to zero at all but the selected frequencies, changes the data to Q f (k). Afterwards, the resulting spectrum is then multiplied by the continuous-time frequency-domain representation of differentiators to velocity and acceleration by jω(k) and −ω(k) 2 , respectively, at the selected frequency, where ω(k) = 2πkf s /K , K is the number of samples of the signal in one sampled period, k is the number of the kth sample and f s is the sample frequency. Finally, a transformation back into the time domain using the inverse discrete Fourier transform (IDFT) yields the estimates of the first and second time derivative of the joint position, that are the velocity and acceleration. The velocity signal and acceleration signal are almost free of noise, that is, noise is removed from all frequencies except the selected ones. The measured averaged joint torques τ J are filtered by the low-pass Butterworth filter.

III. RETRIEVAL OF PHYSICAL PARAMETERS
The physical parameters are suitable combinations of geometric and inertial data of the robot bodies, which specify the mass, the position of the center of mass, and the symmetric inertia matrix for each robot link. For link i, the physical parameters are defined as (12) where m i is the mass of link i, and r C i x , r C i y , r C i z T is the vector from the center of the i-th link frame to the center of mass of link i. f vi and f ci are the viscous and Coulomb friction coefficients of link i, respectively. I ixx , I ixy , I ixz , I iyy , I iyz , I izz T form the link i inertia tensor matrix I C i with respect to the center of mass of link i The relationship between the dynamic parameters π i in (3) and the physical parameters µ i in (12) are shown below. The relationship of mass of link i between (3) and (12) is M i = m i , the relationship of product of the mass and centroid vector is and the relationship between inertia tensor I i related to the origin of the i-th link frame and I C i related to the center of mass of link i is In fact, only the estimated dynamic parametersπ are not sufficient and an estimation of the physical parametersμ is needed at some time. For example, the inertia matrix M (q) is needed for calculating the residual torque errors during robot collision detection. The physical parameters are needed because the inertia matrix computed by the dynamic parameters by Newton-Euler method can not satisfy the property of positive-definite matrix. To obtain consistent physical parameters, a global optimization method and framework with constraints guaranteeing physical feasibility, i.e., positive masses and positive definite inertia tensors, is proposed in this paper. Instead of using the identified base parameters, the objective function f (µ) for the constrained nonlinear optimization problem, namely the modulus of predicted torque errors vector between the averaged measurement torque and the estimated torque by the identified essential parameters, is defined as where µ = µ T 1 , . . . , µ T n T is physical parameters set of robot. It should be noted that all the measured data in the right hand side of (15) have been averaged by (11). Y e is the corresponding stacked observation matrix for the essential parameters π e , which can be obtained by eliminating the column of Y b with respect to the ignored poorly identification parameters in π b . π e (µ) is the nonlinear relationship between π e and µ. Firstly, the relationship π(µ) can be established by (14), then the regrouping relationship π b (π ) can be established by the method described in [12]. Finally, the relationship π e (π b ) can be established by eliminating the poorly identified parameters in π b .
To obtain a feasible solution, the lower bounds LB and upper bounds UB should be considered, which are referred to bounds described in [29].
In order to guarantee a physically meaningful result, constraints on the physical parameters µ are introduced. Physical constraints regard the mass of each link to be positive, namely m i > 0, which can be implemented by the lower bounds LB. In addition, the total sum of the link masses must be in a given range, which is described as Besides, the inertia tensor of each link should be positive definite, which impose the triangular inequality of the diagonal elements of the inertia tensor matrix [37] and prevent big differences between them. The constraints of the inertia tensor are described as In (18), the biggest diagonal term of the inertia can not exceed to the smallest one too much. Moreover, the smallest diagonal term of the inertia of links one to five is forced to correspond to the axis parallel to the link length. It should be noted that the fifth constraint in (18) is only used to make the identified physical parameters compatible with the iiwa robot model in Gazebo in our case. In fact, it is not necessary in the proposed framework. If the reader's robot is not an iiwa robot, this constraint can be ignored. Finally, the nondiagonal elements of the inertia tensor matrix should be small compared to the diagonal ones. In order to describe the retrieval algorithm of physical parameters more clearly, the pseudo-code of the global optimization framework proposed in this paper for the optimization problem is shown in Table 1. The first step of the algorithm is to load the averaged and stacked vector τ J and Y e , which are collected and calculated by the sensors reading of iiwa robot. Then, the constrained nonlinear optimization problem is solved by κ times, at a given step k = 1, . . . , κ, in a loop. At every step in the loop, the initial value for the optimization iteration is randomly selected and updated between the lower and the upper bounds using a uniform distribution. Afterwards, the nonlinear optimization problem is solved by a global optimization algorithm with considering the lower and the upper bounds and the constraints mentioned above. In this paper, the applied global optimization method is the multi-start algorithm because it is a global optimization method and it contributes to obtaining better optimal solution in our case compared with other global optimization algorithms in MATLAB, such as GlobalSearch. Multi-start approach can find a global solution or multiple minima solutions, which starts a local solver (such as fmincon) from multiple start points. The iteration process will end when the termination conditions are satisfied, such as the maximum iteration numbers. The proposed parameter retrieval framework makes contribution to finding a global optimization solution compared with those local optimization methods. In addition, it is easy to introduce some constraints, such as the triangle inequality of the inertia tensors, eventually ensuring the algorithm to obtain a meaningful solution. More experiment results can be found in next section.

IV. EXPERIMENT RESULTS AND VALIDATION A. KUKA LBR IIWA 14 R820 ROBOT
The KUKA LBR iiwa 14 R820 robot is equipped with n = 7 revolute joints and a torque sensor is mounted in each joint after each of the gearboxes. Moreover, the robot is equipped with joint position encoders. Figure 2 and Table 2 show the iiwa robot and its kinematic parameters according to the modified Denavit-Hartenberg (MDH) convention, respectively. It is possible to control the iiwa robot by KUKA Sunrise Toolbox [40], iiwa stack [41] or other software packages. In this paper, the authors control the robot through the iiwa stack, which can support Robot Operating System (ROS) programming and help the users to control the robot easily. This software package is able to provide functions to read the sensors data of joint position q and joint torque τ J . These data will be of fundamental importance for the identification of the dynamic model of the robot. Moreover, it supports lots of control modes, such as joint position control mode, joint velocity control mode, and impedance control mode in both joint space and Cartesian space, which can be selected by the users based on their requirements.

B. ROBOT EXCITATION, DATA ACQUISITION AND SIGNAL PROCESSING
In the parametrized Fourier series (8), the number of sine and cosine terms is set to L = 5 and the fundamental frequency to f f = 0.05 Hz here. Therefore, the duration of the trajectory is 20 s . The joint position and joint velocity limits of (10) used in this paper are listed in Table 3. By considering the optimization problem described by (8)∼(10), the optimal excitation trajectories are designed and shown in Figure 3, which gives a condition number of cond(Y b ) = 156.3. In this paper, our iiwa robot is fixed on a wooden desk, therefore, large motion velocity of the robot will make the desk vibrate and deteriorate the collected data. So, for a trade-off between good data measurement and low condition number, the velocity limits are reduced compared with the ones in iiwa manual, which makes the condition number more than 100.   described in (11) is used for the measured position and torque. From the averaged position data, the joint velocity and joint acceleration are calculated by the frequency filtering method described in the end of Section II. In addition, the computed joint acceleration of the robot for the excitation trajectories is shown in Figure 4 by using the frequency filtering method and the central differentation method, respectively. The results show that the obtained acceleration signal by the frequency filtering method is more clear than the one by the classical central differentation method. The averaged joint torque is filtered by a low-pass Butterworth filter with an order of 4 and a cutoff frequency of 6 Hz.

C. ESTIMATION OF ESSENTIAL PARAMETERS
By using the averaged and filtered joint position and torque and the estimated joint velocity and acceleration calculated above, the observation matrix Y b can be computed and    Table 4 together with their RSD, which are indicators of the quality of the individual parameter estimates. However, from the Table 4, it can be found that the RSD of the four parameters, XY 5 , XY 6 , MX 6 , FV 5 , are large and the values of the four parameters are close to zero. Therefore, their contribution to the dynamics model is negligible and the four parameters can be removed when calculating the essential parameters. The calculation procedure has been described at the end of Section II-A. Finally, the identified 53 essential parametersπ e are given in Table 5 together with their RSD, which are all no more than 40%. Figure 5 compares the averaged measured torque for the excitation trajectories with the predicted torque based on the identified base parameters and essential parameters. The curves indicate that the joint torque can be predicted accurately for all joints of iiwa robot FIGURE 5. Model validation based on comparison between mesaured torque filtered by a low-pass filter (black lines) and the predicted torque by the identified essential parametersπ e (red lines), predicted torque by the identified base parametersπ b (blue dashed lines) and predicted torque by the identified physical parametersμ (green dashed lines) for the excitation trajectories. except for the seventh joint. One reason for this is that the mass and absolute range of torque of joint seven are much smaller compared to the other joints, which makes it hard to predict the torque of seventh joint accurately.

D. ESTIMATION OF PHYSICAL PARAMETERS
The physical parameters µ in (12) of the iiwa robot have been retrieved by solving the nonlinear optimization problem presented in Table 1. In order to solve the nonlinear optimization problem, the algorithm multi-start and its related functions in MATLAB are used. The optimization algorithm in Table 1 is launched for κ = 5 times. In every iteration loop, the maximum iteration times of the objective function are set to 100000. The total optimization time is 9854.5 s in VOLUME 8, 2020 TABLE 6. Lower bounds (LB) and upper bounds (UB) of the link mass, position of center of mass and the friction coefficients of iiwa robot and the optimal solutionμ obtained in this paper and the oneμ [37] reported in [37].
MATLAB on a laptop computer with i5-8250U-1.60Ghz and 16 GB memory. For the bounds of link mass, a non-negative mass for each link is considered and a common upper limit, such as 6 kg, is assumed for simplicity. By considering the light seventh joint with a payload in our case, the upper mass limit for it is increased to 5 kg. The total sum of the link masses m rob,min and m rob,max are set to 16 kg and 26 kg, respectively. In addition, for each link, the center of mass is simply assumed to be located inside the smallest parallel box which includes the link geometry [29]. For the bounds of the elements in the inertia matrices, the mass distribution and the position of the center of mass for each link is evaluated with the aid of CAD tools. By assuming the friction coefficients as non-negative, the upper limits are set to 1 and the lower limits of viscous friction coefficients are set to 0.1. Trading off between the desire of strict bounds, so as to limit the feasible set when searching for a candidate solution, and the need of larger bounds, not to cut off any potential good candidate, the used LB and UB are eventually selected in Table 6 and  Table 7. In order to cover the largest surface of f (µ) within the feasible set, the starting point µ 0 at each search iteration is randomly selected by µ 0 = LB + (UB − LB)u, u ∈ (0, 1).
The final optimization solutionμ obtained with the proposed optimization framework is reported in Table 6 and  Table 7. The objective function value f (μ) at the optimal solution is 75.06 N .m. From the table, it can be seen the solution is in the range of the LB and UB. In addition, the mass and the elements of inertia matrix in the optimal solution satisfy the constraints (17) and (18) completely. The predicted torque estimated by the identified physical parametersμ for TABLE 7. Lower bounds (LB) and upper bounds (UB) of the elements of the link inertia matrix I of iiwa robot and the optimal solutionμ obtained in this paper and the oneμ [37] reported in [37].
the excitation trajectories are plotted in Figure 5. Compared with the measured torque, the physical parameters can also predict the joint torque accurately like the identified base and essential parameters. In general, the reported solution µ cannot be qualified as the 'true' physical parameters. Nonetheless, it satisfies the dynamics model equations, and can be safely adopted, e.g., to compute the inverse dynamics by means of a Newton-Euler algorithm.

E. IDENTIFIED RESULTS VALIDATION
In order to validate the joint torque prediction accuracy by the identified base parameters, essential parameters and physical parameters, three reference trajectories q ref ,j (t), j = 1, 2, 3, are designed and the comparison is done between the measured joint torque and the predicted torque calculated by the estimated parameters based on robot tracking the reference trajectories. The excitation trajectories in Figure 3 is selected as the first reference trajectory, q ref ,1 (t), for a direct validation and the comparison between the measured torque and predicted torque is shown in Figure 5. The results show the estimated joint torque by the identified parameters, which includes base parameters, essential parameters and physical parameters, can be predicted accurately compared with the actual measurement joint torque. The second reference trajectory is designed as a cosine trajectory q ref ,2 (t) = A (cos ωt − B), in which A = [ π 12 , π 6 , π 12 , π 6 , π 12 , π 6 , π 12 ], B = 0.1 and ω = 0.03. The third reference trajectory is designed with very slow joint velocity. In order to quantify the errors between the measured torque and predicted torque, the root mean square (RMS) errors are computed and the resulting RMS values for the three reference trajectories are given in Table 8. From the table, it can be found that the total RMS values of predicted torque errors are approximate calculated by the identified parameters for each reference trajectory, which indicates high precision of the joint torque prediction by the identified parameters. In addition, the RMS values in the third reference trajectory are largest and the ones in the second reference trajectory are smallest. One reason is that the viscous and Coulomb friction model used in this paper is not accurate enough. The final reference trajectory with lower motion velocity and acceleration leads to more complex joint friction torque, such as inclusion of the Stribeck effect after motion just starting. However, the simple viscous and Coulomb friction model used in this paper cannot model it completely. Therefore, the joint torque prediction accuracy for the final reference trajectory is worse than that for the first two reference trajectories. The reason for the second reference trajectories with smallest RMS is that the measured torques of all joints except for the seventh joint in the second reference trajectories are the smallest ones compared with other two trajectories.

F. COMPARISON WITH OTHER METHODS AND AN APPLICATION FOR THE IDENTIFIED PHYSICAL PARAMETERS
In this subsection, the results of comparison of estimated torque is given in Figure 6 and Table 9 between the proposed method and the methods in [30] and [37]. The predicted torques in Figure 6 are calculated by the identified physical parameters obtained by the three different methods based on the excitation trajectories data and the same friction model in our case. From the figure and table, it can be seen that the torque estimates from the proposed method is better than that of previous methods in [30] and [37]. For example, the total RMS value and the objective function value of the nonlinear optimization in the proposed method are all the smallest one according to Table 9.
In addition, an application for calculating the momentum observer used for robot collision detection by the identified physical parameters is given in Figure 7 and Table 10 based on the three different methods. For more details about the momentum observer and collision detection refer to [5]. The data for calculating the residual signal is based on the second reference trajectories. It should be noted that there is no actual FIGURE 6. Estimated torque comparison between measured torque filtered by a low-pass filter (black lines) and the predicted torque by the identified physical parameters from the proposed method (red lines), the method in [30] (blue dashed lines) and the method in [37] (green dashed lines) for the excitation trajectories, respectively. collision occurring during the robot motion in the second reference trajectories. In ideal condition, the residual signal in Figure 7 should be zero when there is no collision for robot. VOLUME 8, 2020 FIGURE 7. Comparison of momentum observer output calculated by the identified physical parameters from the proposed method (red lines), the method in [30] (blue dashed lines) and the method in [37] (green dashed lines) for the second reference trajectories, respectively. However, in fact, the residual signal is not zero even without collision because of the inaccurate dynamics. To detect a possible collision, a constant detection threshold for each joint should be set which should be larger than the corresponding maximum residual value in Table 10. Otherwise, false positives would happen. For example, the thresholds of joint 4 can be set to 0.5 (>0.4140), 0.7 (>0.6224) and 0.9 (>0.7976) according to the values in Table 10 calculated by the proposed method and the methods in [30] and [37], respectively. It should be known that larger thresholds will reduce detection sensitivity. Therefore, the residual signal estimated by the parameters obtained from our method can contribute to higher detection sensitivity than the ones by other two methods. In general, collisions usually occur at the last four joints, therefore, more attention should be focused on them.
Compared with the method in [30], the main difference is that the objective function, predicted torque errors, is calculated based on the identified essential parameters and the corresponding observation matrix in our case. However, the objective function was calculated based on the identified base parameters in the supplementary material of [30]. Compared with using the base parameters to calculate objective function, one direct advantage of using essential parameters with high identification quality is that the predicted accuracy of joint torque is significantly improved, such as the RMS value and objective function value in Table 9. The reason that the predicted results shown in [30] are accurate may be that it benefits from using the reverse engineering method which can contribute to obtaining base parameters with lower RSD values, and using the more complete friction model. Thanking to the method in [30], it gave us inspiration to propose our framework in this paper.
Compared with the method in [37], one difference of the proposed method is using a global optimization instead of a local optimization. In addition, the objective function of the nonlinear optimization is computed based on the essential parameters and the corresponding observation matrix, not the base parameters in [37]. One advantage of the proposed method is that the predicted accuracy of joint torque is improved compared to that of the method in [37]. In addition, another advantage is that the identified physical parameters by the proposed method can be used to obtain higher collision detection sensitivity compared with the ones calculated by other two methods in [30] and [37]. Both the identified physical parameters in this paper and the ones reported in [37] are shown in Table 6 and 7. Although the values of two sets of parameters are different, both of them can be used to predict the joint torques. It should be noted that, however, compared with the iiwa robot used in [37], there is a difference that our iiwa robot is fixed with a payload. Therefore, it can be seen that m 7 in our case is larger than the one in [37] in Table 6. Of course, payload contribution can be also removed by the method in [42] to obtain comparable parameters for the mass of link 7. In addition, although the used LB and UB are not shown in [37], one point can be certain that our LB and UB are different from theirs. For example, m 7 is equal to 0.3543 in [37], which is obviously not in the range (1 ∼ 5) of LB and UB in our case. By exploring the relationship between the identified elements of the link inertia matrix in Table 7 and the constraints in (18), it can be found that both the optimal ones in this paper and the ones of [37] are satisfied with the triangular inequality constraints.

V. DISCUSSION AND CONCLUSION
In this paper, the problem of extracting the physical parameters, which described the dynamics of the robot, was addressed by the proposed global optimization framework. The base dynamic parameters with RSD values were firstly identified by the standard least squares method. Based on the results of identified base parameters, the essential parameters were calculated and obtained by ignoring the parameters close to zero with high RSD values. Afterwards, based on the identified essential parameters with better identification quality, a set of physical parameters were retrieved by solving a nonlinear optimization problem, taking into account several constraints including the physical bounds on the total mass of the robot and the triangle inequalities of the link inertia tensors. Finally, the proposed framework for physical parameters extraction and all the identified parameters were validated by comparison experiments. In addition, a comparison between the proposed method and the methods in [30] and [37] is given. Then, an application of the identified physical parameters for calculating momentum observer is shown. The experiment results show that the identified physical parameters by the proposed method can obtain better prediction accuracy of joint torque and contribute to obtaining higher collision detection sensitivity compared with the methods in [30] and [37].
In our proposed framework, the physical constraints can be easily added and avoids using complex penalty functions like [30]. A major feature of the proposed framework is that it is a universal algorithm of the physical parameters retrieval for serial manipulators. In addition, it can be easily modified to include further nonlinear and physical constraints. In the nonlinear problem, the objective function was computed by using the identified essential parameters with better identification quality, instead of the base parameters. The total RMS value and the objective function value at the optimal solution by the proposed method were all the smallest when compared with the methods in [30] and [37]. This means that the proposed framework obtains more accurate torque prediction and better optimization solution by using the essential parameters. The identified physical parameters in Table 6 and Table 7 can satisfy both the set bounds and physical constraints. Therefore, the optimal solution is a set of physically meaningful solution. The solution can be regarded as the 'optimal' one, but it is not the 'true' one. It is possible to get the true solution only if the strictly accurate bounds and constraints can be added. Nonetheless, the optimal solution in this paper satisfies all dynamics model equations and it can be safely adopted, such as to compute the inverse dynamics, collision detection and so on by means of Newton-Euler method or Lagrange method. The base parameters, essential parameters and physical parameters of iiwa robot are all identified and given in this paper. To the best of the authors' knowledge, this is the first paper that show all the three kinds of dynamic parameters of an iiwa robot, which can benefit to the researchers in robotics to select the preferred one for their research. The total RMS values in Table 8 are approximate for each reference trajectory, which indicates the high precision of the joint torque prediction of iiwa robot by all the three sets of identified parameters.
The bounds used in this paper are not strictly accurate, which limits the possibility of accessing to the true physical parameters solution. In addition, more meaningfully physical constraints need to be further explored. Therefore, in the future work, the authors would like to explore more accurate bounds and consider other possibly physical constraints, and use the computed dynamics by the identified parameters for robot collision detection and reaction control. In 1989, he joined the Department of Mechanical Engineering, National University of Singapore, where he is currently an Associate Professor and the Acting Director of the Advanced Robotics Centre. He teaches both the graduate and undergraduate levels in the following areas: robotics; creativity and innovation, applied electronics, and instrumentation; advanced computing; and product design and realization. He is also active in consulting work in these areas. His research interests include robotics, mechatronics, and applications of intelligent systems methodologies. In addition to academic and research activities, he is actively involved in the Singapore Robotic Games as its Founding Chairman and the World Robot Olympiad as a member of the Advisory Council.