Closed-loop Dynamic Parameter Identification of Robot Manipulators Using Modified Fourier Series

This  paper  concerns  the  problem  of  dynamic parameter  identification  of  robot  manipulators  and proposes  a  closed‐loop  identification  procedure  using modified Fourier series (MFS) as exciting trajectories. First, a static continuous friction model is involved to model joint friction  for  realizable  friction  compensation  in  controller design.  Second, MFS  satisfying  the  boundary  conditions are  firstly  designed  as  periodic  exciting  trajectories.  To minimize  the  sensitivity  to  measurement  noise,  the coefficients  of  MFS  are  optimized  according  to  the condition  number  criterion. Moreover,  to  obtain  accurate parameter  estimates,  the maximum  likelihood  estimation (MLE) method  considering  the  influence of measurement noise  is  adopted.  The  proposed  identification  procedure has  been  implemented  on  the  first  three  axes  of  the QIANJIANG‐I  6‐DOF  robot  manipulator.  Experiment results verify  the effectiveness of  the proposed approach, and comparison between identification using MFS and that using  finite Fourier series  (FFS)  reveals  that  the proposed method achieves better identification accuracy.


Introduction
Recently, robot manufacturers have increased development efforts towards new applications, e.g., machining, laser welding, laser cutting and multi-robot collaboration [1], which puts forward new performance requirements on robots. The requirements, which mean high speed and high precision, can be satisfied by involving advanced model-based control schemes [2,3]. Design of an advanced model-based controller needs an accurate knowledge of the dynamic model. While the model structure of robot manipulators is well known, the involved parameter values are not always available, since dynamic parameters are rarely provided by the robot manufacturers and often not directly measurable. Therefore, dynamic parameter identification of robot manipulators has aroused increasing interest from researchers [4][5][6].
A robot manipulator is a multi-input multi-output (MIMO) system with coupling effects and nonlinearities, which makes the identification problem particularly complex. Since the robot manipulator is always open-loop unstable, a closed-loop identification procedure is needed, in which the robot executes specific exciting trajectories through closed-loop feedback control. An identification procedure generally consists of modelling, exciting trajectory design, data measurement, signal processing, parameter estimation, and model validation [7], where exciting trajectory design, measurement accuracy, and parameter estimation method mainly affect the identification accuracy [8]. With a review of previous studies, several identification approaches have been proposed. Atkeson, et al. [9] developed a procedure using a fifth-order polynomial trajectory in joint space as exciting trajectory, and standard least squares estimation (LSE) was adopted as the parameter estimation method. Armstrong [10] firstly studied the optimization of the exciting trajectories. He introduced the condition number of the observation matrix as a measure of the noise immunity, and the optimization was regard as a nonlinear optimization problem. Swevers, et al. [11] proposed a periodic exciting trajectory, parameterized as finite Fourier series (FFS). The periodicity of FFS enables continuous repetition of the identification experiment and allows time domain averaging of the measured data, which improves the signal to noise ratio. This procedure has an important effect on subsequent research. However, the boundary conditions are not satisfied for FFS, which makes it difficult for robots to track. Considering the influence of measurement noise, Olsen and Petersen [12] proposed the maximum likelihood estimation (MLE) method based on a statistical framework for parameter estimation. To improve the measurement accuracy, Bingul and Karahan [8] used a motion analysis system with three cameras to measure the robot velocity and acceleration, and load-cell sensors to measure the joint torques. Different from the identification procedures mentioned above, Qin, et al. [5] suggested a sequential identification procedure. The friction, gravity, and inertial parameters were identified step by step, while an accumulation of errors may occur since identification errors were propagated from one step to the next.
In this paper, a closed-loop identification procedure using modified Fourier series (MFS) as exciting trajectories is developed. First, instead of the basic friction model that is the combination of viscous and Coulomb friction, a static continuous friction model is involved to model joint friction for realizable friction compensation in controller design. The robot dynamic model is then expressed linearly with respect to the dynamic parameter vector to be identified (Section II). Second, motivated by FFS in [11], MFS are designed, and the coefficients of MFS are optimized using genetic algorithm (GA) (Section III). Moreover, to estimate the dynamic parameters accurately, MLE is adopted as parameter estimation method (Section IV). The identification experiment and model verification process have been carried on the first three axes of the QIANJIANG-I 6-DOF robot manipulator. Experiment results verify the effectiveness of the proposed identification procedure. In addition, comparison between identification using MFS and that using FFS reveals that the proposed approach can achieve better identification accuracy (Section V).

Robot Dynamics
The dynamics of a rigid serial n-link robot manipulator with revolute joints can be written as where qR n denotes the joint angle vector, M(q)R n×n represents the inertia matrix, C(q, q  )R n×n represents the centripetal-Coriolis matrix, G(q)R n is the vector of gravity effect, τ f R n represents the friction effect and τR n is the vector of input torque.
Robot manipulators are usually actuated by motors with harmonic drive reducer or RV reducer, and for most applications, it is sufficient to regard the robot joint friction as a static nonlinear function of the joint velocity with Stribeck effect. Thus, the friction of joint where f ci is the level of Coulomb friction, f si is the level of stiction force, si q  is the Stribeck velocity,  is an empirical parameter which is dependent on the specific situation. ξ can be equal to 2, but not necessarily, and f vi is the viscous friction coefficient. As shown in Fig.1, this model is discontinuous at zero velocity. Therefore, we cannot use this model for friction compensation in controller design since it would be sensitive to quantization errors and measurement noise at zero velocity and may excite high frequency dynamics. In addition, when the velocity changes direction, the sudden jump of friction compensation can never be delivered by the motor. To solve this problem, the discontinuous friction model (2) is approximated by a continuous function (as shown in Fig.1) which can capture the essential characteristics of friction and can be realized by the motor. Thus the static continuous friction model can be obtained as 20 where F fi ( i q  ) is a known nonlinear continuous function, A fi is the unknown amplitude to be identified.
Dynamic model (1) contains kinematic parameters, inertial parameters and friction coefficients. Kinematic parameters, i.e., link lengths, twist angles, and link offsets, are either provided by the robot manufacturer or can be measured directly on the robot. Inertial parameters are rarely provided, but can be identified by proper procedures. Inertial parameters of the ith link, including link masses, first-order mass moment, and link moment of inertia, can be defined by the following vector: m m r m r m r I I I I I I  p (4) where m i is the mass of the link, m i r xi , m i r yi and m i r zi are the components of the first-order mass moment, I xxi , I xyi , I xzi , I yyi , I yzi , I xzi are the components of the inertial tensor. Friction coefficients are also unknown and need to be identified. Dynamic model (1), excluding friction, can be formulated linearly with the inertial parameters [9]: where P=[p 1 T p 2 T … p n T ] T R 10n×1 is inertial parameter vector with p i (i=1,2, … ,n) defined by Eq. (4). κ(q, q  , q  )R n×10n is so-called regress matrix, and elements of which are nonlinear functions of joint position, velocity, and acceleration. In fact, not all the inertial parameters in P have an effect on the dynamic model, which means only a part of them can be identified. The inertial parameters can be categorized into three classes [13]: uniquely identifiable, identifiable in linear combination, and unidentifiable. By eliminating those that have no effect on the dynamic model and regrouping some others, we can obtain a base parameter set (BPS) which is a set of minimum number of identifiable inertial parameters [14]. It is proved that the dynamic model excluding friction is linear with respect to BPS [15] where P b R p×1 (p <10n is the number of elements in BPS) is the vector of base parameters, ( , , )  q q q    R n×p is corresponding regress matrix. Note that the static continuous friction model (3) also can be expressed linearly with the unknown amplitude A fi and viscous coefficient f vi . Therefore, Eq. (6) can be rewritten as where is a dynamic parameter vector including BPS and friction coefficients. Φ(q, q  , q  )R n×(p+2n) , whose elements can be expressed as nonlinear functions of motion data, is the corresponding regress matrix. Thus, the robot identification can be formulated as dealing with the problem of estimating the dynamic parameter vector β from the data measured when the robot is executing the exciting trajectories. In most cases, the measured data are joint angles q and joint torques τ. Assume we have collected K samples of q and τ corresponding to time instants t 1 , t 2 , … , t K . Since Eq. (7) holds for each sample instant, we obtain following overdetermined set of equations: ,

Design of Exciting Trajectories
Measurement noise, friction phenomenon not captured by the proposed friction model, and other un-modelled dynamics (e.g., elasticity in transmissions) may cause bias of the estimated parameters. To reduce the effect of the measurement noise, the exciting trajectories must be carefully designed. During design, it is necessary to consider noise immunity, sufficiency of the excitation, motion constraints of the robot, and convergence rate of the parameter estimation. Generally, the design procedure consists of two steps. First, parameterize the trajectories. Second, optimize the trajectories with motion constraints according to suitable optimization criterion.

Modified Fourier series
Several approaches have been presented in parameterization of exciting trajectory, e.g., fifth-order polynomial trajectory [9] and FFS [11]. In [11], the FFS of the ith joint is as follows: where ω f is the fundamental frequency of the exciting trajectory and is common for all joints to guarantee the periodicity of the robot excitation. N is the number of harmonics. q i0 is the constant term. a l i , b l i (i=1, 2,…, n) are the coefficients which constitute the DOF of the trajectory optimization. FFS has advantages in terms of signal processing and noise rejection, however, there is a sudden change of velocity and acceleration at the start and end of motion, i.e., (0) ≠0 (i=1, 2, … , n; t f is the period), which may cause vibration of the robot. Therefore, it is difficult for robots to track closely with FFS, and the identification accuracy may be deteriorated. To remedy this disadvantage, the exciting trajectories are required to satisfy boundary conditions which can be described as where q init is the initial configuration of the robot manipulator. Eq. (10) means that the robot manipulator moves back to initial configuration at the end of each period, so that it can repeat the exciting trajectory period by period continuously. Eqs. (11), (12) guarantee that the velocity and acceleration are continuous at the start and end of the motion. To satisfy these requirements, we construct a quintic polynomial to substitute for the constant term q i0 in Eq. (9). Then, MFS of ith joint can be obtained as where t f is the period, j (j=1, 2, 3, …) represents the jth period and      is the top integral function. Eq. (13) specifies a periodic trajectory with period 2π/ω f . Solving Eqs. (10) (11) (12) and (13), the six coefficients c i k (k=0, 1, … ,5) can be represented in terms of coefficients a i l ,b i l and q init : int ( , , ).
Substitute (14) in (13), Eq. (13) can be expressed only containing coefficients a i l ,b i l . Therefore, we can remedy the disadvantages of FFS without adding DOF of the trajectory optimization. Coefficients a i l ,b i l can be determined through a trajectory optimization procedure in the next section.
The proposed MFS has following advantages: i) Satisfaction of the boundary conditions allows the robot manipulator to execute the exciting trajectories exactly.
ii) The periodic property allows time-domain data average, which improves the signal-to-noise ratio of the measured data.
iii) Continuous repetition of MFS allows estimation of characteristics of the measurement noise that is valuable in parameter estimation.

Trajectory optimization using GA
To find optimal exciting trajectories, several criteria have been proposed in previous literature. In [16], trajectory optimization is executed according to the d-optimality criterion. In [17], the exciting trajectories were optimized by minimizing the condition number and maximizing the smallest singular value of the observation matrix. In this work, minimizing the condition number as in [2] is adopted. The condition number of the observation matrix represents the upper bound for input/output error. It directly affects the convergence rate and noise immunity of the identification experiment. Thus, the optimization problem can be formulated as minimizing the condition number of the observation matrix. Meanwhile, the motion constraints on the robot manipulator need to be considered. The optimization can be described as where { ( ( ))} t s q calculated using forward kinematics represents the set of end-effector positions when the robot is executing the exciting trajectories. S is the available workspace. min q , max q , max q  , and max q  are the vectors containing the physical limits of joint position, velocity, and acceleration.
This constrained nonlinear optimization problem can be solved using GA with binary coding. During the GA optimization, each decision variable a i l and b i l is coded as a binary integer, respectively. Then, an individual is obtained by connecting all the binary integers into a binary string. The condition number of the observation matrix is chosen as the fitness function. At each generation, the fitness value of each individual is evaluated, and individuals are selected for reproduction based on the fitness values. The selected individuals undergo recombination under the action of the crossover and mutation operators. Moreover, in decoding phase, the nonlinear constraints are taken into account. After a certain number of generations, GA converges to an individual that is best suited.

Pre-process of data
The designed exciting trajectories are repeated by the robot manipulator period by period continuously, while joint angles and joint torques are sampled at specified frequency. Decentralized PID controllers are applied to generate the control inputs. As we know, the collected position and torque data are corrupted by measurement noise, which causes uncertainty and bias errors in parameter estimation. Therefore, the collected data should be pre-processed before using them in parameter estimation.
Since the exciting trajectories are periodic, the following data averaging can be used to improve the signal-to-noise ratio: where k represents the kth measurement in each period. M is the number of measured periods. q mj (k) and τ mj (k) are the kth measurements within the jth period.
To estimate the noise level on measured data, the variance of the measurement noise is estimated as where σ qi 2 , σ τi 2 are variance of noise on measured positions and torques of joint i. K is the number of measurements in each period. q mij (k) and τ mij (k) are the kth measurement within the jth period of joint i.
In addition, actual information of joint velocity and acceleration are required in calculation of the observation matrix in parameter estimation. Usually, joint velocities and accelerations are not available for lack of sensors and need be estimated from joint positions. Since numerical differentiation of joint positions yields very noisy joint velocities and accelerations that cannot be used in parameter estimation, we adopt an analytical approach as suggested in [11]. The measured joint positions are first approximated as MFS using the linear least square method. Then the joint velocities and accelerations are calculated as the derivatives of the obtained MFS.
Remark: this method should only be used if the robot manipulator follows the exciting trajectories closely.

Maximum likelihood parameter estimation
The most commonly used method in parameter estimation is LSE that optimizes the root-mean-square residual error of the model. However, measurement noise often deteriorates the accuracy of parameter estimation. To solve this problem, we need to consider the effect of measurement noise, and MLE based on the statistical framework is adopted.
For measurement noise, Capisani, et al. [18] analysed the noise on a COMAU SMART3-S2 industrial manipulator, which demonstrates that the noise distribution can be considered practically to be Gaussian. Also, in system identification based on the statistical framework, it is common to assume the measurement noise to be independent zero-mean Gaussian noise [6,12]. Thus, the measured joint position and joint torque can be written as where n qi (k) and n τi (k) are independent zero-mean Gaussian noise on the position and torque measurements of joint i. q i (k) and τ i (k) represent the noise-free joint position and joint torque which satisfy equation (7). The maximum likelihood estimate of parameter vector β is the value of β that maximizes the likelihood of the measurements. Since the measurement noise is independent, the maximum likelihood of β can be obtained by minimizing following cost function: where σ qi 2 and σ τi 2 are known constants calculated in Eqs. (18) and (19). Considering Eqs. (7), (20) and (21), the minimization of criterion (22) is a nonlinear least squares optimization problem. The DOF of this optimization is the number of elements in β. Solving this problem requires calculation of n qi (k) and n τi (k) for every estimate of β given the measured joint positions and torques. However, this is impossible since the noise-free joint position q i (k) and joint torque τ i (k) are unknown.
To solve this problem, Swevers, et al. [16] proposed a method that extended the parameter vector β with the coefficients of exciting trajectory. However, this method requires an initial guess of the parameters, and the convergence to global optimum is not guaranteed. Actually, if the joint position data are free of noise, the estimation can simplify significantly, and Eq.(22) can be written as This assumption can be justified by the fact that the measurement noise on joint positions is much smaller than that on joint torques, and most of the noise on joint positions is filtered through data averaging in the preprocess. Therefore, taking into account of Eqs. (7), (21) and (23), the MLE simplifies to weighted linear least squares estimation, and the maximum likelihood estimate of the parameter vector β can be calculated as , Σ is the diagonal variance matrix of the measurement noise on joint torques, and Γ m =[τ m (t 1 ) T τ m (t 2 ) T … τ m (t K ) T ] T .

Model verification
To assess the quality of the identified model, verification trajectories different from the exciting trajectories are executed by the robot manipulator. Meanwhile, joint torque is sampled. Joint torque predictions τpred are calculated from Eq. (7) using the estimated parameter vector βML and the desired motion data. Then, the torque predictions and measured torques are compared. Fig. 2 shows the scheme of experimental procedure used to verify the identified model. The quality of the identified model can be evaluated by the root-mean-square (RMS) residuals which is defined as where τ' m (k) is the kth sample of the joint torque when the robot manipulator executing the verification trajectories, τ pred (k) is the predicted torque.

Experimental Results and Comparisons
To verify the effectiveness of the proposed approach, the identification procedure was implemented on the first three axes of the QIANJIANG-I 6-DOF robot manipulator (see Fig. 3), and the remaining axes were kept locked. The six rotational joints of the robot manipulator are actuated by AC servo motors incorporating increment encoders for position measurements. The joint torque data are obtained through motor current measurements. The control unit of the robot manipulator is constituted of an industrial PC and a motion control card. The motion control card, on which control algorithms are implemented, is open for programming. The reduced dynamic model (7) for the first three axes is derived, where the continuous function F fi ( i q  ) is chosen as (2/π)arctan(950 i q  ) and the parameter vector β contains 21 base parameters (15 inertial parameters and 6 friction coefficients).
The exciting trajectories are postulated as in (13) shows the corresponding trajectory of end-effector in Cartesian space (only one period of the trajectories is shown). As can be observed, the boundary conditions as well as the motion constraints are exactly satisfied. The parameter vector β are identified through procedure using FFS [11] and that using MFS, respectively. The resulted parameters are shown in   To verify the accuracy of the identified parameters, we chose a circle with diameter 100mm in Cartesian space as the verification trajectory, since circular motion is commonly used in robot manipulator operation.
Trajectory planning was carried out in Cartesian space, and the trajectories in joint space have been calculated through inverse kinematics. During the motion, the maximum velocity of the robot end-effector is 0.5m/s in Cartesian space. Fig. 8 shows the desired circle and actual circle performed by the robot manipulator through PID control. During the motion, the joint torque has been sampled. In addition, the joint torque predictions are calculated using the identified parameters and the desired motion data. Fig. 9 compares the measured and predicted torque, which reveals that the needed joint torques can be reconstructed using the dynamic parameters obtained from the identification procedure using MFS. For the torque prediction errors, they are mainly caused by the bias of estimation, friction phenomenon that are not captured by the proposed friction model, and un-modelled dynamics.    Table 3 shows the RMS residuals of the joint torque prediction. As can be observed, identification using MFS allows smaller prediction errors than that using FFS. The reason is that MFS satisfy the boundary conditions exactly, which makes it possible for robot manipulators to track closely. Therefore, accurate estimates of joint velocity and acceleration can be obtained. Moreover, since the desired MFS have been optimized for noise immunity, tracking closely is helpful to reduce the influence of measurement noise. From table 3 we can also observe that the improvement of torque prediction accuracy of joint 1 and joint 2 is more obvious than joint 3. This is because the first two joints have a large load and severe coupling effect. For such joints, it is much easier to track with MFS than FFS. While, joint 3 has a relatively small load and is not influenced by other joints.

Conclusions
In this paper, a closed-loop identification procedure using MFS has been proposed for dynamic parameter identification of robot manipulators. A static continuous friction model is involved to model joint friction for realizable friction compensation in controller design. MFS satisfying the boundary conditions are firstly designed as exciting trajectories. Considering motion constraints on robot manipulators, the coefficients of MFS are optimized according to the condition number criterion for minimizing the sensitivity to measurement noise. To improve the identification accuracy, the effect of measurement noise is considered in parameter estimation, and MLE based on the statistical framework is adopted. The proposed approach has been implemented on the first three axes of the QIANJIANG-I robot manipulator. Experiment results show that the needed joint torques can be reconstructed using the dynamic parameters resulting from the proposed identification procedure. Moreover, comparison between identification using MFS and that using FFS reveals that the former achieves better identification accuracy than the latter.

Acknowledgments
The research work is jointly supported by the Zhejiang Provincial Natural Science Foundation of China