Robust Parameter Estimation of Robot Manipulators Using Torque Separation Technique

In this paper, we propose a robust method for estimating the parameters of robot manipulators using the torque separation technique, which was developed previously by the authors, to extract the inertial, gravitational, and frictional components from the input torques of multiple sinusoidal joint motions. The separated components of the input torque produce a set of reduced linear regression equations where the dynamic parameters of the robot manipulators appear decoupled or minimally coupled. A mathematical analysis is presented to show that the set of reduced regression equations tends to yield more robust parameter estimation than the conventional way where a large dimensional full regression equation is solved simultaneously. An iterative scheme to alleviate the possibility of parameter estimation error caused by the deviation from the assumed ideal sinusoidal motions is also proposed. Experimentation with two robot manipulators is used to verify the proposed approach and related claims. The proposed method is simple and pragmatic without requiring any specialized signal filtering and technical know-hows which numerous previous methods often demanded for explicitly or implicitly.


I. INTRODUCTION
Estimating the dynamic parameters of robot manipulators requires a rational dynamic model, a persistently exciting reference trajectory, reliable signal processing for measured torque/motion data, and a judicious choice of a mathematical algorithm to solve the necessary regression problem [1]- [7]. If either of these elements is not handled properly, the resulting estimated parameters will fail to represent the actual dynamics of the robot manipulators, while these issues often require technical know-how that the practitioners or field engineers do not have easy and explicit access to. This makes the dynamic parameter estimation still remain a difficult and tricky task, even though the topic has been studied for such a long time for precise motion and force controls for various robotic applications.
In this paper, we propose a novel parameter estimation method that is not only easy to implement but also reliable in performance. The proposed method estimates the dynamic The associate editor coordinating the review of this manuscript and approving it for publication was Aysegul Ucar . parameters in a dissected manner by using the torque separation technique [8], developed previously by the authors, which enables the whole joint torque to be divided into inertial, gravitational, and frictional components through the superposition of the sinusoidal joint responses with different cyclic frequencies. By doing so, the proposed method becomes to have better robustness to the perturbations and uncertainties occurring from various sources. Besides, our method is simple and straightforward to implement in the field, as it does not depend on any particular method of signal filtering nor require any art of choosing good excitation trajectories, while many previous methods did.
The property that the robot dynamic model appears to be linear in the parameters [9], [10] is commonly used to transform a time history of dynamic equations into a set of standard linear regression equations, AX = B, where A is a long rectangular regression matrix, X denotes a vector of unknown parameters, and B is a vector of measured joint torques. Depending on whether the dynamic parameters are chosen minimally or non-minimally, the form of the regression equation could vary [5], [11]. The non-minimal VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ parameters consist of standard link parameters such as mass, the first mass moment, and the inertia tensor of all links, as well as the frictional coefficients of all joints [12]- [15]. Because of structural restrictions, some of the non-minimal parameters in the regression equation become dependent on others, resulting in the condition that the regression matrix becomes rank deficient. On the contrary, the minimal parameters, also known as the base parameters [16], imply the least number of parameters that can describe the dynamics of the robot manipulators [17]. The minimal parameters can be obtained by linear or nonlinear combinations of non-minimal parameters. Khalil and Gauthier [16] proposed a method for symbolically obtaining the minimal parameters, which was realized by a publicly open software, OpenSYMORO [18]. Minimal parameters can also be determined numerically by eliminating the dependent set of parameters from non-minimal parameters using techniques such as the QR factorization or the singular value decomposition of the regression matrix [1], [9]. With minimal parameters, the regression equation becomes a normal least squares problem that yields a unique solution [19], [20]. The regression equation, whether expressed by a set of minimal or non-minimal parameters, has hinted and inspired a numerous amount of studies on topics such as, estimating the best set of parameters from the mathematical least squares viewpoint [6], [7], [10], [21]- [23], designing optimal test trajectories that fully excite the robotic system by examining the maximum likelihood and condition number of the regression matrix [24]- [27], combining state estimators to use less noisy velocities and acceleration information [28]- [30], or characterizing the statistical effect of the noise and uncertainties to achieve favorable parameter estimation [31]- [33].
Recent progress in robot parameter estimation research has been, in particular, focused on how to estimate the nonminimal parameters while maintaining physical feasibility, such as the positivity of mass and inertia tensor. Sousa and Cortesão [11] interpreted the dynamic parameter estimation problem as a semi-definite programming problem with linear inequalities that impose the physical conditions of mass and inertia tensor. In estimating the dynamic parameters, Joukov et al. [34] used the extended Kalman filter approach, in which the physical feasibility conditions were imposed. Lee et al. [35] used a differential geometric method to effectively develop the robot parameter estimation as a convex optimization problem with physical constraints and validated their method in the parameter estimation of highly articulated and complex systems. Wensing et al. [36] developed physical consistency constraints as a linear matrix inequality in terms of the covariance of the mass distribution.
In principle, most aforementioned methods attempted to estimate the dynamic parameters of a robot manipulator all at once by solving a single linear regression equation. However, these one-shot methods will naturally lack robustness if no special care is taken. The reason is because the pure inertial torque component in most conceivable robotic motions is much smaller than the gravitational and frictional components, and bad conditioning of the regression matrix with a large condition number can happen [24], [37], [38]. Therefore, the least squares solution would be sensitive to the perturbation or noise in the regression matrix elements and can lose reliability.
Unlike the conventional one-shot approaches, this study proposes a separated approach for estimating the inertial, gravitational, and frictional parameters individually from a set of reduced regression equations, which enhances the robustness without losing accuracy of estimation. This set of reduced regression equations is established with the assistance of the torque separation technique developed by the authors [8]. A mathematical analysis is supplemented to show how the proposed method can improve the reliability of the estimated parameters in representing the dynamics. Additionally, a complementary error correction algorithm is also proposed, which further improves the estimation accuracy when the sinusoidal responses are compromised by the tracking error.
Before closing, a short review of studies on robust parameter estimation in diverse fields is given as follows. Studies in [39]- [41] addressed the robustness of the parameter estimation enhanced by modeling and compensating the structured disturbance with a predictable form. For the unstructured disturbance, the stochastic approach was used to model and compensate the disturbance; see, for example, [42]- [44]. In [45], [46], Stewart derived the sensitivity of estimated parameters in linear estimation problem against the structured and/or unstructured disturbances, which are distinguished by how the disturbances occur in relation to the range space of the regression matrix. Golub and Pereyra proposed the separable nonlinear least squares method by separating the regression matrix into the range and orthogonal space terms [47]. In addition, a robust regression equation can be used to maintain the performance against modeled or unstructured uncertainties. Many other robust linear regression methods, such as M-estimates, S-estimates, Least Median of Squares(LMS) just to name a few, have been studied by considering the stochastic properties of the disturbance and the induced output. See, for example, [48], [49]. The filtering techniques were also treated importantly in eliminating the disturbance and uncertainty in the data signal in the regression problems [50], [51].

II. PRELIMINARIES: DYNAMICS MODEL AND PARAMETER REGRESSION A. DYNAMIC MODEL
Consider the dynamic equation for an n-DoF (degree of freedom) robotic manipulator written by the tensor notion as n j=1 m ij (q(t))q j (t) + n j=1 n k=1 h ijk (q(t))q j (t)q k (t) where q(t) = [q 1 (t) q 2 (t) · · · q n (t)] T ∈ R n is the vector of joint variables with q i (t) being the i-th joint angle; τ i (t) is the i-th joint torque; m ij (q(t)) is the inertial coefficient of the j-th link applied to the i-th joint, h ijk (q(t)) := ∂m ij (q(t))/∂q k (t)− 1/2∂m jk (q(t))/∂q i (t) is the coefficient of Coriolis and centripetal acceleration; g i (q(t)) denotes the i-th component of the gravitational torque; α i and β i are the coefficients of the Coulomb's friction and the viscous friction, respectively; and sgn(q i (t)) is the signum function. Please refer to, for example, [52]- [54] for detailed information on the terms of (1). (In the sequel, the explicit time dependency in the state variables shall be omitted if no confusion arises in the context.) The fundamental dynamic parameters relevant to the rigid link i of the robot manipulator, shown in Fig.1, are defined as where elements from I i,xx to I i,zz represent the six independent components of the inertia tensor defined by the body-fixed frame i , m i is the mass, r i ∈ R 3 is the position vector of the center of mass from the origin of i . Because these parameters appear linearly in the dynamic equation of the robot manipulator [9], [10], [14], the dynamic equation in (1) can be rewritten using a linear regression form as where Y( ) ∈ R n×12n is the regression matrix with [q TqTqT ] T being the concatenated vector of joint variables and their first two derivatives; and = [ T 1 T 2 · · · T n ] T ∈ R 12n denotes the vector of dynamic parameters for all links. Extra dynamic parameters such as the rotor inertia of each motor and higher-order frictional coefficient may appear in the parameter vector , but we ignore them as their effect on the overall dynamics is insignificant.
Certainly the reliability of the dynamic model (1) or (3) is dependent on accurate knowledge of , which would then help us design various kinds of model-based controllers including inverse dynamics based control [53], [54], model predictive control [55], Lyapunov-based control [56], and optimal control [57].

B. CONVENTIONAL PARAMETER ESTIMATION
When a sufficient amount of motion data ( and τ ) is available at instances, t = k t, k = 1, 2, · · · N , where t denotes the sampling time, a large dimensional regression equation is constructed as where Y(k t) implies Y( ) at t = k t, and Y ∈ R nN ×12n is the assembled full regression matrix, which is vertically very long. This equation is posed as one of the least squares problems, defined as follows: where the weight W −1 physically implies the inverse of the noise covariance W ∈ R N ×N . Given that Y is probably rank deficient due to inherent structural mobility restriction in some proximal links, including the first link to be attached on earth, infinitely many feasible solutions inevitably exist. One may choose the minimum norm solution among the feasibles.
To briefly overview the associated mathematical nature, the common QR factorization is applied to Y as follows where P ∈ R 12n×η is the permutation matrix; η is the rank of Y; Y 1 ∈ R nN ×η and Y 2 ∈ R nN ×(12n−η) are, respectively, the decomposed base and redundant regression matrices. Matrices Q 1 ∈ R nN ×η and Q 2 ∈ R nN ×(12n−η) are the unitary sub-matrices of the column space of Y; R 1 ∈ R η×η and R 2 ∈ R η×(12n−η) are the corresponding coordinate matrices. The above permutation matrix P allows the following transformation: where s ∈ R η is a primary parameter set and r ∈ R 12n−η is a redundant parameter set. Then the set of base parameters b ∈ R η , representing the minimal set of parameters required to describe the dynamics, satisfies And b is determined by solving the following least squares problem redefined from (5): The transformed least squares problem must have a unique solution, which we treat as a de facto set of estimated parameters. However, the obtained solution b invites many possible 's that satisfy (8), which depends on the choice of r [19]. The minimum norm solution among many 's can be the VOLUME 9, 2021 most preferred candidate if no other criteria exist. Please refer to the details of the algorithm in the literature, e.g., [16].
There is another class of approaches tackling the least squares problems with rank-deficient regression matrices, which has been known as the Tikhonov's regularization method [35], [58], [59]. By this method, the best is determined by minimizing a single cost function that combines the squared residual error and squared deviation of from a certain value 0 , as follows: where the first part of the cost penalizes the residual error, whereas the second one, being the Tikhonov's regularization term, keeps suppressing the solution to remain near 0 which can either be a nominal value of from CAD data or be just null if no a priori knowledge of is available. And µ denotes a constant weight. The closed-form solution to (10) can be written as where I n denotes the identity matrix of order n. The regularization term has transformed the original rank-deficient least squares problem into the one that is computationally wellposed, saving an effort required to handle the redundancy of . The trade-off is a possible reduction in regression accuracy.
Recently an interesting strand of researches was done, focusing on combining additional equality and/or inequality constraints into the least squares problem in order to find a physically more consistent and reliable solution that accounts for the requirements such as positivity of mass and inertia, triangular inequality of the mass-inertia parameters and so on. See references, for example, [11], [14], [35]. These arrangements could be effective to some extent but the fundamental and intrinsic difficulty of selecting a particular set of dynamic parameters from an infinite number of possibilities may not be avoided, and the estimated parameters may not be effective enough to reproduce the dynamics for different reference trajectories.

III. PARAMETER ESTIMATION BY TORQUE SEPARATION TECHNIQUE
In this section, a new pragmatic approach to achieving more reliable parameter estimation is presented by using the torque separation technique, which allows us to separately identify the parameters into relevant groups, thereby robustifying the regression accuracy.
A. TORQUE SEPARATION TECHNIQUE: REVIEW [8] A brief overview of the torque separation technique begins with the following sinusoidal joint trajectory: where q 0 is the initial joint configuration of the manipulator, a = [a 1 a 2 · · · a n ] T ∈ R n is the vector of amplitudes of the sinusoidal joint trajectory, and ω is the frequency of the whole joint motion. If the robot manipulator is forced to move along the above joint trajectory, input torque i satisfies where τ ω i (t) and q ω (t) denote τ i (t) and q(t), respectively, to highlight the associated sinusoidal trajectory (12) with frequency ω. If we double and triple the frequency of the joint trajectory in (12), the joint torque becomes τ 2ω i (t) and τ 3ω i (t), respectively. Because q ω (t) = q 2ω (t/2) as depicted in Fig.2, Similarly, because q 2ω (t/2) = q 3ω (t/3), h ijk (q ω (t))a j a k sin 2 ωt + β i a i ω sin ωt.
By subtracting (14) from (15) and then dividing by 2, the inertial component of torque for trajectory (12) can be determined as follows.
The frictional torque can be separated from the full dynamics by using the symmetric properties of the sinusoidal functions, cos ωt = cos ω(T − t) and sin ωt = − sin ω( and dividing the result by 2 give us the frictional torque as And, after some algebraic manipulation, the gravitational torque is separated as follows: Thus, the inertial, frictional, and gravitational torques have been completely isolated simply by superposing three dynamic responses with a set of sinusoidal joint trajectories. And obviously these separated torque components reconstruct the whole joint torque τ ω i (t) if added.

B. SEPARATED TORQUE BASED PARAMETER ESTIMATION
As illustrated in (16) - (18), each separated component of joint torque is equivalent to a fraction of dynamics associated with a reduced number of parameters. To examine further, the whole parameters are divided into groups as where i,inr , i,grc , i,frc are the dynamic parameters of the i-th link related to the inertial, gravitational and frictional components of torque, respectively, defined as (16), the linear regression form only for the inertial torque is obtained as where τ ω inr (t) = [τ ω 1,inr (t) τ ω 2,inr (t) · · · τ ω n,inr (t)] T is the vector of inertial torques of all joints, and Y inr (·) and Y inr (·), respectively, denote the sub regression matrices related to inr and grv . From (18), we also obtain τ ω grv (t) = Y grv (q ω , t) grv (20) where τ ω grv (t) = [τ ω 1,grv (t) τ ω 2,grv (t) · · · τ ω n,grv (t)] T is the vector of gravitational torques of all joints, and Y grv (·) is the sub regression matrix related to grv . Similarly, from (17) τ ω where τ ω frc (t) = [τ ω 1,frc (t) τ ω 2,frc (t) · · · τ ω n,frc (t)] T denotes the vector of frictional torques of all joints, and Y frc (·) is the sub regression matrix related to frc .
Each of the relations (19) - (21) can be stacked up for instances at t = t, 2 t, · · · , N t, to yield a set of separated regression equations with reduced size as where the symbols( ·) denote the concatenation of the corresponding quantities for all instances. From (23) and (24), grv and frc can be independently estimated by solving the following small-sized optimizations: and arg min The solutions to the above optimizations can be solved, for example, using the QR factorization method. Once these have been solved, inr from (22) can be estimated by solving the following: where grv is supposed to be known from the previous step. Each of the weighting matrices, W inr , W grv and W frc , correspond to the noise covariance of each separated regression equation [60], [61]. These optimizations can be developed differently by taking other choices of cost functions and additional constraints, but the separated nature of these optimizations is still maintained. The main advantage of the separated sequential optimizations (25)-(27) over a one-shot optimization (5) is improved reliability against noises and/or perturbations in the measured data. The improved reliability is achieved because the coupling effect can be avoided between groups of parameters, particularly between inr and grv , which have different sensitivity under the uncertainties inherently, thereby reducing the possibility of overfitting effect that conventional one-shot parameter regression methods are likely to suffer from.
To validate our claim, we first define the order of magnitude for each sub-regression matrix in the following proposition, using its explicit symbolic form summarized in Appendix A. (For more complete derivation of the symbolic regression matrix for robot manipulators, please refer to [14].) Proposition 1: In the robot dynamics with a sinusoidal joint trajectory (12), the order of magnitude of sub-regression matrix Y i , i = {inr, inr , grv, fric} is determined as where O(·) denotes the order operator and g r is the gravitational constant.
The order of magnitude is defined based on the dominant features that make each separated regression matrix, Y i , differ from others. Looking closely into the symbolic forms, Y inr and Y inr depend linearly on the angular acceleration and the square of angular velocity, in which the square of the sinusoidal frequency ω of the joint trajectory defines their magnitude; thus it is reasonable to say Y inr and Y inr is of O(ω 2 ). (Even in other kinds of joint trajectories for parameter estimation, they can be regarded as O(ω 2 ) if the trajectories are approximated as sinusoidal ones.) In Y grv , the peculiar difference is that it has the factor of gravitational constant g r ; therefore Y grv is of O(g r ). Similarly, Y frc must be of O(ω) because it depends on the joint velocity.
Within mind the orders of magnitude of these regression matrices, consider the robustness of the conventional regression problem in (5) when perturbations δ Y and δ τ occur in the regression matrix and joint torque, respectively, due to disturbances and uncertainties. The change of solution δ , affected by δ Y and δ τ , satisfies the following: where δ Y R and δ Y N , respectively, denote the projections of δ Y onto the range space of Y and onto the orthogonal complement of the range space; similarly, τ R and τ N denote, respectively, the projections of τ onto the range space of Y and its orthogonal complement; δ τ R is the projection of δ τ onto the range space of Y; and κ(·) denotes the condition number of a matrix defined as σ M (·)/σ m (·), where σ M (·) and σ m (·) represent the maximum and minimum non-zero singular values, respectively. See [45], [46] for detailed information on how the relation in (29) has been derived. According to the above equation, δ is affected by two components related to the projection of the range and orthogonal space of δ Y. The term related to δ Y N is magnified by a factor κ 2 ( Y)|| τ N ||/|| τ R ||, which is the ratio measuring how nearly τ is on the range space of δ Y. The disturbance on the range space of δ Y is called structured (parametric) disturbance which is linked to the error of parameters. Besides, the disturbance on the orthogonal space of δ Y is called the unstructured disturbance which could not be compensated by model in regression. If the structured disturbance is dominant in the whole disturbance, || τ N ||/|| τ R || goes to zero, and the sensitivity follows κ( Y). Otherwise, if the unstructured disturbance is dominant, the sensitivity follows κ 2 ( Y). The higher condition number makes the effect of the disturbance significant and causes a greater error in estimating the parameters whether the disturbance is a structured or unstructured.
The same logic can be applied to each of the three separated sub-regression problems in (25) - (27). Therefore, if all the sub-regression matrices have much smaller condition numbers than the full regression matrix, we can conclude that separated parameter estimation is more robust to the perturbations than the conventional one-shot approach. The following assumption plays as the prerequisite for the sub-regression matrices to have much smaller condition numbers.
Assumption 1: For the sinusoidal joint motion with ω 2 g r , the relative magnitudes of the maximum and minimum non-zero singular values between Y i , i = {inr, grv, frc}, satisfy the following inequalities: Loosely speaking, Assumption 1 implies that the magnitude of Y grv is greater than those of Y inr and Y frc , which may be conjectured from Proposition 1. Because the elements of Y inr and Y frc , containing joint velocity and acceleration, are varying sinusoidally, their magnitudes in average sense become much smaller than those of Y grv which depend linearly on the gravitational constant g r .
Proposition 2: The condition number of the full regression matrix Y is always much larger than the condition number of Y i , i = {inr, grv, frc}, that is, if the speed of the sinusoidal trajectories satisfies ω 2 g r . Proof: Consider the following rearranged full regression equation as which is equivalent to (4) but with a different sequence of parameters. Then, the self product of the above Y becomes where Y grv Y grv + Y inr . On condition that ω 2 g r , certainly Y grv ≈ Y grv since Y inr ∼ O(ω 2 ) and Y grv ∼ O(g r ). Now, by Cauchy's eigenvalue interlacing theorem [62], the eigenvalues of the principal submatrix Y T i Y i , i = {inr, grv , frc} must hold the following inequalities: , and a similar analogy exists between Y T i Y i and Y i . Combining the conditions in Assumption 1 with the above inequalities yields which lead to the conclusion: So far, the better robustness of the separated parameter estimation approach has been looked at from the mathematical viewpoint. However, another practical reason for the robustness of the proposed approach can be found as well. In the proposed approach, grv is estimated on the basis of the column space of Y grv in (23), which is dependent only on measurements of joint angles that are rarely affected by noises. Therefore. using cleaner measurements, grv would be estimated more reliably, which is in return expected to enhance the reliability of inr in (22) as inr is affected by the estimated grv . In conventional methods where all the parameters are estimated simultaneously, the noises and perturbations in the regression matrix can produce a compromised solution for both inr and grv -let alone the frictional parameters frc -that simply minimizes the residual error of the coupled regression equation (4) in the brute-force manner. Given that the portion of the inertial torque in nominal operations of the manipulator is much smaller than that of the gravitational torque, any small perturbation in the motion data can significantly reduce the reliability of the estimated inr , possibly resulting in poor reproduction of inverse dynamics when highly inertial movements are required.
Practical proposed approach is that it does not any special sophisticated signal processing techniques, which frequently need complicated tuning of filter parameters [60], [63], [64]. Additionally, it is conceptually simple and general enough to be used in the parameter estimation of almost all robotic manipulators available.

C. COMPENSATION OF JOINT TRACKING ERROR
Recall that the torque separation and the ensuing separated parameter estimation took advantage of the harmonic properties of sinusoidal trajectories with an assumption that the robot manipulator is expected to track the sinusoidal trajectories perfectly. In practice, such an ideal or near-ideal tracking performance may be difficult to achieve for custommade robot manipulators before fine controller tuning, which, however, relies on a precise knowledge of the dynamic parameters. Thus, the proposed parameter estimation method should cope with the possibility of poor joint tracking. In the sequel, we present a recursive algorithm that uses the idea of the Newton-Raphson method [65] to update the estimated parameters in such a way to compensate for the effect of the tracking error.
Before going into details, we first define τ i ( , ) as the i-th inverse dynamic torque that explicitly shows the dependency on and . Certainly τ i ( , ) should be equal to the measured joint torque τ i (t) or τ i for short. Then the Taylor series expansion of the joint torque at a given time t can be written as follows.
where d denotes the desired joint trajectory for the tracking task as defined in (12);ˆ denotes the known previous estimate of parameters, which was estimated in previous iteration; and [∂τ i /∂ ] * and [∂τ i /∂ ] * represent the partial derivatives of τ i for and , evaluated at ( = d , = ). The first calculation ofˆ is done under the condition with = d , via the least squares algorithm in (25)-(27). We know from (3) However, the second term on the right-hand side of (31) does not reduce to a simpler form. Thus, we use the first-order Taylor series approximation of τ i ( ,ˆ ) for about d to obtain Substituting (32) and (33) into (31) and then rearranging the terms yield where the term on the left-hand side is the inverse dynamic torque of the true system for the ideal sinusoidal trajectory, the first term on the right-hand side must be the measured joint torque, and the rest of the two terms on the right-hand side are the reproduced inverse dynamic torques of the system with untrue parameterˆ for actual and ideal joint trajectories, respectively. Because all terms on the right-hand side are known or computable, the joint torque of the true system for the ideal sinusoidal trajectory is achievable. The last two terms in the right-hand side of (34) act as the corrective torque to produce joint torque for the ideal sinusoidal trajectory when the joint tracking error is not negligible. When (34) is applied to the joint torques under sinusoidal joint motions, corrected sets of inertial, gravitational, and frictional components of joint torques can be separated even when the joint responses deviate from the ideal sinusoidal ones. The extracted torque components can then be used to construct regression equations in (19) -(21) just as before, whence the updated estimates for inr , grv , and frc are determined sequentially. To improve the quality of torque Algorithm 1 Procedure of the proposed parameter estimation algorithm with joint error compensation Input: N sample sets of motion data including desired and actual joint trajectories ( d , ), and measured joint torques τ (t) with (t ∈ { T , 2 T , · · · , N T }) for each of three sinusoidal trajectories of cyclic frequencies ω, 2ω, and 3ω; Output: Estimates of dynamic parameters (ˆ inr ,ˆ grv ,ˆ frc ); (τ ω inr , τ ω grv , τ ω frc ) ← Enhanced torque separation using the updated joint torque; 7: Break if ||τ ( d , ) − τ ( rω , )|| < ; 8: end for End procedure separation, repeat the above steps based on the updated estimates of the partitioned parameters. The overall flow of the proposed scheme is summarized in Algorithm 1.
A potential limitation on the proposed the recursive error compensation algorithm is that the estimated solution might not converge uniformly when the tracking error is too large. Because the algorithm is based on the idea of the Newton-Raphson's method, the joint torque expression in a first-order linear polynomial form (31) is not valid under the condition of a large tracking error, and therefore, the proposed algorithm might not work as intended. A possible remedy is to slow down the speed of numerical update in (34) as where K Y < 1 is a small gain that suppresses the effect of the abrupt change.

IV. EXPERIMENTAL RESULTS
In this section, the torque separation-based parameter estimation method will be demonstrated through two case studies: one using the Frank Emika Panda (https://www.franka.de/) the other one using a custom-made robot manipulator [66]. The step-by-step procedure for the proposed method is detailed and then the performance of the method is compared with the conventional method.

A. CASE STUDY 1: FRANKA EMIKA PANDA
The Franka Emika Panda is an affordable research-oriented collaborative robot manipulator with 7 joints, and is equipped with a joint torque sensor in each joint. Our experimental setup with the Franka Emika Panda is shown in Fig.3. This robot is controlled by the Franka Control Interface (FCI) under the Robot Operating System (ROS) at 1[kHz] of sampling speed. The FCI controller receives either a joint position command or a joint torque command, depending on its operational mode, and has feedback information on the position, velocity, and torque of each joint.
Six cases of joint trajectory q(t) = q 0 + a(1 − cos(rωt)), with different q 0 's and a's, were chosen to collect motion data  for parameter estimation, as summarized in Table 1. And each was sequentially executed with three different frequencies, rω, r = 1, 2, 3 where ω = π/3. The collection of these trajectories covered a large portion of robot's workspace. Especially, the motion of the wrist joints was made as wide as possible to draw a recognizable inertial effect, because the torques of the wrist joints were usually much smaller than other joints'. The robot manipulator was made to track each  trajectory 10 times repeatedly by using the original position servo mode of FCI, and the ensemble average of the joint responses was obtained for reliability. Fig.4 shows the joint response for joint trajectory set 1 with the base frequency ω = π/3. The observed joint tracking error was in the order of 10 −3 [rad], but the magnitude of the joint tracking error tended to increase for higher cyclic frequencies. The measured joint torque was quite smooth even without any signal smoothing (see Fig.4(c) and (d)). Since joints 5-7 were placed at the wrist, less torque in those joints was produced during the motion.
Next, the joint torque data obtained from the tests with frequency rπ/3, (r = 1, 2, 3) were combined to extract τ ω inr (t), τ ω grv (t), and τ ω frc (t) as described in (16) - (18). Fig.5 shows these decomposed components for one period of cycle (T = 6[s]). A moving average filter was used to smooth the decomposed components. As shown in Fig.5, τ ω grv (t) took the largest proportion from the whole joint torque for all joints except for joints 1 and 7. (Note that joint 1 was parallel to gravitational direction and joint 7 had a link with negligible mass.) The proportions of τ ω inr (t) and τ ω frc (t) were found to be similar, which seemed unusual given that the frictional torque in most cases would be much larger than the inertial torque. This means that the mechanical condition of the Frank Emika Panda was well managed.
In order to estimate the dynamic parameters, we constructed the full regression matrix ( Y) in (4) and the reduced regression matrices ( Y inr , Y inr , Y grv , and Y frc ) in (22) -(24) as well. The necessary joint velocities and accelerations were derived using the numerical differentiation with a simple moving average that suppresses the numerical noise. These regression matrices were a conglomeration of the experimental results from all six sets of trajectories. The QR decomposition was then used to solve the related rank-deficient least squares problems. Table 2 shows the rank, the maximum/minimum non-zero singular values (σ M and σ m ) and the condition number (κ). Under the given experimental conditions, the sinusoidal frequency (ω = π/3) was slightly higher than 1, implying that ω 2 g r was satisfied. In Table 2, σ M of Y inr , Y inr , and Y frc look similar, while σ M of Y grv is approximately g r /ω 2 times higher than σ M of Y inr , Y inr , and Y frc , which means that Assumption 1 was valid in this experimentation. Note that σ M of Y was almost close to σ M of Y grv , but σ m of Y was slightly lower than σ m of Y inr . Therefore, κ of Y is higher than κ's of the separated regression matrices as stated in Proposition 2. Table 3 shows the estimated dynamic parameters obtained from the three approaches: the conventional one-shot approach (CM), the proposed separated torque approach (TS), and the proposed separated torque approach with joint error compensation (TS+EC). The listed parameters are written in the order of the magnitude of the associated singular values of the regression matrices. The base dynamic  parameters from CM were obtained by (9), while the base dynamic parameters from TS and TS+EC were obtained by (25)- (27) after reducing the dependent variable via QR factorization by (6). All noise covariance matrices in the least squares problem were set to be the identity matrices, so that the comparison between the methods could not be affected by the hyper parameters. The algorithms were computed under MATLAB R2019b on a dedicated PC with core i5-7600 3.5GHz, 32GB RAM. The computational times of CM, TS, and TS+EC were, respectively, 14.13, 7.20, and 14.37 minutes. The majority of the computational time was due to the QR factorization of the related regression matrix. However, this amount of time consumption should not be a significant issue as the algorithm is to be run off line, just to extract the parameter values.
To evaluate these sets of estimated parameters, the measured joint torques for a new experimental trajectory tracking task were compared with the corresponding estimated joint torques, as shown in Fig.6. The new test trajectory was made by combining a set of sinusoidal functions of the form: q i (t) = a i sin(ω i t)+q 0,i , i = 1, 2 · · · , 7, which is depicted in Fig.6(a). Table 4 contains the numerical values of the trajectory parameters. The result shows that the estimated torque of the CM differed significantly from the measured torque at the moments of high acceleration, which seemed to be caused by incorrect inertial parameters estimation ( Fig.6(b) -(h)). Especially, the torque estimation errors of joints 1, 3, 5 were remarkably high since these are more influenced by inertial motion. On the contrary, the proposed two approaches were highly accurate in estimating the joint torque for most of the time. A statistical summary of the torque estimation is presented in Fig.6(i). To evaluate the aptness of the estimated dynamic parameters, we used the root mean square (RMS) of torque estimation error e RMS as the performance indicator: where T f denotes the time span of evaluation. The evaluation by the RMS torque error captures the average performance on how correctly the estimated parameters represent the actual dynamics. The RMS torque error in the ascending order was 0.3405 (TS+SC), 0.4226 (TS), and finally 0.6788 (CM). Not only RMS error, but the spread of the torque estimation error was small in the same order among the approaches. These results supported our claim that the proposed separation approaches would be more robust and reliable in the parameter estimation than the conventional approach.

B. CASE STUDY 2: GRAVITY BALANCING ROBOT MANIPULATOR
In this case study, the gravity balancing robot manipulator, the so-called eZRO, developed by the Korea Institute of Machinery & Materials (KIMM) was employed to validate the effectiveness of the proposed parameter estimation method. The eZRO shown in Fig.7 is a 6-DOF robot manipulator which is capable of counterbalancing the weights of  main links through internal sliding masses with spring, even though complete gravity compensation may not be possible due to added weights in the wrist and inevitable uncertainties in the counterbalancing mechanism. The expected mechanical benefits would be that the size of the motors can significantly be reduced, thus the system becomes energy efficient and probably safer for unexpected collision. Refer to [66] for the detailed information about this robot manipulator.
To estimate the parameters of the eZRO using the proposed approach, four sets of sinusoidal joint trajectories, summarized in Table 5, were used to gather the necessary motion data. The cyclic frequencies used in this case study were rω (r = 1, 2, 3) with ω = π/3[rad/s], which were the same as the previous case study. Fig.8 shows a cycle of ensemble-averaged joint responses for the first set of joint trajectories, obtained after 10 repeated tests. In this set of joint trajectories, an identical joint travel path was imposed on each joint. The magnitude of the tracking error ranged in the order of 10 −3 [rad], and the maximum joint torque was as high as 70 [Nm] in the second shoulder joint.
The inertial, gravitational and frictional components of torques were extracted as before by combining all sets of input joint torques in the experiments, as shown in Fig.9. The gravitational torques in the second and third joints were mostly caused by the weights at the wrist joints that were not completely counter balanced. In terms of the proportion, the gravitational torque took less than a half of the overall torque, as opposed to the case with the Frank Emika Panda, where  the gravitational torque component in the load-bearing joints took as high as 90 percent of the total torque. Next, we constructed the related regression matrices and examined the rank, the maximum/minimum non-zero singular values σ M and σ m , and the condition number κ as summarized in Table 6. An extremely large condition number was obtained for the full regression matrix,Ỹ; thus, it is doubtful that the parameters of eZRO can robustly be estimated than those of the Franka Emika Panda by solving the conventional one-shot regression problem. Opposedly, the separated regression matrices have much lower condition numbers than the full regression matrix. Moreover, Table 6 shows that σ M ( Y inr ), σ M ( Y inr ), and σ M ( Y frc ) are of a similar order of magnitude, whereas Y grv is much higher by the order of gravitational constant than the others. Similar patterns were found for the minimum singular values. Each of these results verified the validity of Assumption 1 and Proposition 1 and 2. Table 7 shows the base dynamic parameters estimated by using the three approaches: CM, TS, and TS+EC. The estimated parameters in Table 7 are listed in the order of the magnitude of the associated singular values of the regression matrices. As was in the case of Franka Emika Panda, all noise covariance matrices in the least squares problem were set to be the identity matrices for a fair comparison. The computational times of CM, TS, and TS+EC were, respectively, 2.71, 1.35, and 2.81 minutes.
The measured joint torques were compared with the reconstructed joint torques from inverse dynamics for a new joint trajectory made by a mixture of the fifth polynomial functions with high and low accelerations to validate the correctness of the estimated parameters. As shown in Fig.10, the estimated torque of CM could not correctly produce the joint torques at the moments of high acceleration, which seemed to be caused by the incorrect estimation of the inertial parameters. Like the results of Franka Emika Panda, the torque errors of the joints that are more influenced by inertia (joints 1 and 4) were considerable by CM. To the contrary, the proposed two approaches (i.e., TS and TS+ES) were highly accurate in estimating the joint torques for most of the time, although the performance was not as good as the case of Franka Emika Panda, which seemed due partly to a large amount of noises in the measured torque signal of eZRO. Especially, the peak torque error of TS+EC was just a half of that from CM. A statistical summary of the RMS torque estimation error is shown in Fig.10(i). The RMS torque estimation error in the ascending order was 2.477 (TS+EC), 2.747 (TS), and finally 4.137 (CM). The spread of the estimation error was also lined up in the same order among the approaches.
These results again showed that the proposed approaches based on the separated torques can robustly estimate dynamic parameters, without the use of specialized signal filtering  techniques or other sophisticated mathematical algorithms. Furthermore, numerous additional experimental tests carried out by using different conditions and robot manipulators confirmed that the performance of the proposed methods were consistent and dependable, always excelling the performance of the conventional method.

V. CONCLUDING REMARKS
In this study, we presented a robust parameter estimation method for robot manipulators, based on the separated components of joint torques (that is, the inertial, the gravitational and the frictional components) obtained by superposing a set of sinusoidal responses with different cyclic frequencies. These separated components of joint torques enabled us to formulate a set of reduced regression problems. We showed mathematically that these reduced regression problems by nature can be more robust to the perturbations due to the associated regression matrices having much smaller condition numbers than the conventional regression matrix. Furthermore, we developed a complementary iterative error correction scheme in order to overcome the effect of trajectory deviation from the ideal sinusoidal one. The performance of the proposed parameter estimation method was verified by employing a commercial 7-jointed robot and a custom-made 6-jointed robot. In both experimental verifications, the proposed method outperformed the conventional method in the parameter estimation. We confirm that the proposed method is practical and reliable, and works in most cases of robot parameter estimation without the use of specialized signal processing techniques.

APPENDIX DERIVATION OF REGRESSION EQUATION
To better understand Proposition 1, we review the procedure to transform the Newton-Euler dynamic equation into a regression equation in terms of the generic dynamic parameters. Most of the subsequent derivation is taken from the authors' previous work in [14].
The Newton-Euler dynamic equation for the entire robot system can be written aṡ where ν ∈ R 3n is generalized momentum for all links G ∈ R 3n denotes the gravitational force vector of all links; F c ∈ R 3n and F ∈ R 3n , respectively, are the equivalent constraint forces/moments and non-conservative forces/moments acting on the center of mass. The momentum of the robot manipulator ν may be written as where M(q) = blockdiag[m 1 I 3 · · · m n I 3 | I 1 · · · I n ], where V i ∈ R 3 and ω i ∈ R 3 are, respectively, the linear and angular velocities of link i; J vi and J wi are the linear and angular Jacobian matrices defined about the center of mass of the i-th link; I 3 is the 3 × 3 identity matrix; and I i is the inertia tensor matrix of link i,. We can modify the Jacobian matrix as where J vi is the Jacobian matrix about the link's distal end. Plugging (37) into (36) and then pre-multiplying J T (q) yield By the virtual work principle, we have J T F = τ , and the constraint condition yields J T F c = 0. Thus, (39) becomes M(q)q + C(q,q)q + g(q) = τ where M = J T MJ, C = J TṀ J + J T MJ, and g = J T G.
If written in a summation form, (39) becomes: whereV i is the linear acceleration on the end of link i;ω i i = R T iω i is the angular acceleration with respect to frame of the link i; S(ω i ) := (ω i ×) is a matrix operator of vector cross product; B(d) ∈ R 3×6 is a matrix operator for any vector d ∈ R 3 satisfying I i d = B(d) i,inr .
If the joint friction is added to (41), then the resulting equation becomes arranged as a regression equation: where τ a := τ + τ frc (corresponding to the motor torque)