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

Abstract

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.

Keywords

Closed-loop identification robot dynamics friction exciting trajectories

1. 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-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).

2. Robot Dynamics

The dynamics of a rigid serial n-link robot manipulator with revolute joints can be written as

M (q) \ddot{q} + C (q, \dot{q}) \dot{q} + G (q) + τ_{f} = τ,

(1)

where q ∈Rⁿ denotes the joint angle vector, M ( q )∈R^n×n represents the inertia matrix, C ( q , $\dot{q}$ )∈R^n×n represents the centripetal-Coriolis matrix, G ( q )∈Rⁿ is the vector of gravity effect, τ _f ∈Rⁿ represents the friction effect and τ ∈Rⁿ 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 i (i = 1,2 … n) is

τ_{f i} = [f_{c i} + (f_{s i} - f_{c i}) e^{- {| {\dot{q}}_{i} / {\dot{q}}_{s i} |}^{ξ}}] \cdot sign ({\dot{q}}_{i}) + f_{v i} {\dot{q}}_{i},

(2)

where f_ci is the level of Coulomb friction, f_si is the level of stiction force, ${\dot{q}}_{s i}$ 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

Figure 1.

Joint friction model

τ_{f i} = A_{f i} F_{f i} ({\dot{q}}_{i}) + f_{v i} {\dot{q}}_{i},

(3)

where F_fi ( ${\dot{q}}_{i}$ ) 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:

p_{i} = {[m_{i} m_{i} r_{x i} m_{i} r_{y i} m_{i} r_{z i} I_{xx i} I_{xy i} I_{xz i} I_{yy i} I_{yz i} I_{zz i}]}^{T},

(4)

where m_i is the mass of the link, m_ir_xi, m_ir_yi and m_ir_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]:

τ = κ (q, \dot{q}, \ddot{q}) P + τ_{f},

(5)

where P =[ p ₁^T p ₂^T … p _n^T]^T∈R^10n×1 is inertial parameter vector with p _i (i=1,2, …,n) defined by Eq. (4). κ ( q , $\dot{q}$ , $\ddot{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], thus we have

τ = κ^{'} (q, \dot{q}, \ddot{q}) P_{b} + τ_{f},

(6)

where P _b ∈R^p×1(p <10n is the number of elements in BPS) is the vector of base parameters, $κ^{'} (q, \dot{q}, \ddot{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

τ = Φ (q, \dot{q}, \ddot{q}) β,

(7)

where β =[ P _b^TA_f1f_v1 … A_fnf_vn]^T∈R^(p+2n)×1 is a dynamic parameter vector including BPS and friction coefficients. Φ ( q , $\dot{q}$ , $\ddot{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₁, t₂, …, tK. Since Eq. (7) holds for each sample instant, we obtain following over-determined set of equations:

W β = Γ,

(8)

where $W = [\begin{matrix} Φ (q (t_{1}), \dot{q} (t_{1}), \ddot{q} (t_{1})) \\ Φ (q (t_{2}), \dot{q} (t_{2}), \ddot{q} (t_{2})) \\ ⋮ \\ Φ (q (t_{K}), \dot{q} (t_{K}), \ddot{q} (t_{K})) \end{matrix}]$ is so-called observation matrix and Γ =[ τ (t₁)^T τ (t₂)^T … τ (t_K)^T]^T.

3. 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.

3.1. 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:

q_{i} (t) = \sum_{l = 1}^{N} (\frac{a_{l}^{i}}{ω_{f} l} \sin (ω_{f} l t) - \frac{b_{l}^{i}}{ω_{f} l} \cos (ω_{f} l t)) + q_{i 0},

(9)

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ⁱ, b_lⁱ (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., ${\dot{q}}_{i} (0)$ ≠0, ${\dot{q}}_{i} (t_{f})$ ≠0, ${\ddot{q}}_{i} (0)$ ≠0, ${\ddot{q}}_{i} (t_{f})$ ≠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

q (0)= q_{i n i t}, q (t_{f})= q_{i n i t},

(10)

\dot{q} (0)= 0, \dot{q} (t_{f})= 0,

(11)

\ddot{q} (0)= 0, \ddot{q} (t_{f})= 0,

(12)

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

\begin{array}{c} q_{i} (t) = \sum_{l = 1}^{N} (\frac{a_{l}^{i}}{ω_{f} l} \sin (ω_{f} l t) - \frac{b_{l}^{i}}{ω_{f} l} \cos (ω_{f} l t)) + \sum_{k = 0}^{5} c_{k}^{i} {(t - (j - 1) t_{f})}^{k}, \\ t_{f} = \frac{2 π}{ω_{f}}, \\ j = ⌈ \frac{t}{t_{f}} ⌉, \end{array}

(13)

where t_f is the period, j (j=1, 2, 3, …) represents the jth period and $⌈ \cdot ⌉$ 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ⁱ_k (k=0, 1, …,5) can be represented in terms of coefficients aⁱ_l,bⁱ_l and q _init:

c_{k}^{i} = F (a_{l}^{i}, b_{l}^{i}, q_{int}).

(14)

Substitute (14) in (13), Eq. (13) can be expressed only containing coefficients aⁱ_l,bⁱ_l. Therefore, we can remedy the disadvantages of FFS without adding DOF of the trajectory optimization. Coefficients aⁱ_l,bⁱ_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.

3.2. 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

\begin{array}{c} \min_{a_{l}^{i}, b_{l}^{i}} cond(W), \\ s.t. {s (q (t))} \subset S, \\ q (0) = q_{i n i t}, \dot{q} (0) = 0, \ddot{q} (0) = 0, \\ q (t_{f}) = q_{i n i t}, \dot{q} (t_{f}) = 0, \ddot{q} (t_{f}) = 0, \\ q_{\min} \leq q (t) \leq q_{\max}, - {\dot{q}}_{\max} \leq \dot{q} (t) \leq {\dot{q}}_{\max}, - {\ddot{q}}_{\max} \leq \ddot{q} (t) \leq {\ddot{q}}_{\max}, \end{array}

(15)

where ${s (q (t))}$ calculated using forward kinematics represents the set of end-effector positions when the robot is executing the exciting trajectories. S is the available workspace. $q_{\min}$ , $q_{\max}$ , ${\dot{q}}_{\max}$ , and ${\ddot{q}}_{\max}$ 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ⁱl and bⁱ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.

4. Parameter Estimation

4.1. 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:

\bar{q} (k) = \frac{1}{M} \sum_{j = 1}^{M} q_{m j} (k),

(16)

\bar{τ} (k) = \frac{1}{M} \sum_{j = 1}^{M} τ_{m j} (k),

(17)

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

σ_{q i}^{2} = \frac{1}{K M - 1} \sum_{j = 1}^{M} \sum_{k = 1}^{K} (q_{m i j} (k) - {\bar{q}}_{i} (k))^{2},

(18)

σ_{τ i}^{2} = \frac{1}{K M - 1} \sum_{j = 1}^{M} \sum_{k = 1}^{K} (τ_{m i j} (k) - {\bar{τ}}_{i} (k))^{2},

(19)

where σ_qi², σ_τi² 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.

4.2. 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

q_{m i} (k) = q_{i} (k) + n_{q i} (k),

(20)

τ_{m i} (k) = τ_{i} (k) + n_{τ i} (k),

(21)

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:

C_{ML} = \frac{1}{2} \sum_{i = 1}^{n} \sum_{k = 1}^{K} (\frac{n_{q i}^{2} (k)}{σ_{q i}^{2}} + \frac{n_{τ i}^{2} (k)}{σ_{τ i}^{2}}),

(22)

where σ_qi² and σ_τi² 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

C_{ML} = \frac{1}{2} \sum_{i = 1}^{n} \sum_{k = 1}^{K} \frac{n_{τ i}^{2} (k)}{σ_{τ i}^{2}} .

(23)

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 pre-process. 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

β_{ML} = {(W_{m}^{T} \sum^{- 1} W_{m})}^{- 1} W_{m}^{T} \sum^{- 1} Γ_{m},

(24)

where $W_{m} = [\begin{matrix} Φ (q_{m} (t_{1}), {\dot{q}}_{m} (t_{1}), {\ddot{q}}_{m} (t_{1})) \\ Φ (q_{m} (t_{2}), {\dot{q}}_{m} (t_{2}), {\ddot{q}}_{m} (t_{2})) \\ ⋮ \\ Φ (q_{m} (t_{K}), {\dot{q}}_{m} (t_{K}), {\ddot{q}}_{m} (t_{K})) \end{matrix}]$ , Σ is the diagonal variance matrix of the measurement noise on joint torques, and Γ _m=[ τ _m(t₁)^T τ _m(t₂)^T … τ _m(t_K)^T]^T.

4.3. 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.

Figure 2.

Model verification scheme

The quality of the identified model can be evaluated by the root-mean-square (RMS) residuals which is defined as

ε_{RMS} = \sqrt{\frac{1}{K} \sum_{k = 1}^{K} {({τ^{'}}_{m} (k) - τ_{pred} (k))}^{2}},

(25)

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.

5. 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 ( ${\dot{q}}_{i}$ ) is chosen as (2/π)arctan(950 ${\dot{q}}_{i}$ ) and the parameter vector β contains 21 base parameters (15 inertial parameters and 6 friction coefficients).

Figure 3.

QIANJIANG-I robot manipulator

The exciting trajectories are postulated as in (13), with N=5 and ω_f=0.1Hz. This specifies a cycle period of 10s. The number of periods repeated is 10. The motion constraints of the QIANJIANG-I robot manipulator are listed in table 1. In the table, Z_end is the height of the end effector in Cartesian space. To avoid collision, Z_end is limited between 400 and 1500. The coefficients aⁱ_l,bⁱ_l of MFS were found using GA for constrained nonlinear optimization which minimizes the condition number of the observation matrix. The condition number of the resulting trajectories is 145. Figs. 4-6 show the resulting optimized exciting trajectories in joint space, and Fig. 7 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 table 2. In this table, ai (i=2,3) is the length of link i, d₂ is the offset of link 2.

Table 1.

Motion constraints of robot manipulator

	Joint	Min	Max
Position(deg)	1	−40	220
	2	−190	60
	3	−126	86
Velocity(deg/s)	1	−100	100
	2	−95	95
	3	−100	100
Acceleration(deg/s²)	1	−45	45
	2	−40	40
	3	−75	75
Z_end (mm)	–	400	1500

Figure 4.

Joint position of the optimized exciting trajectories

Figure 5.

Joint velocity of the optimized exciting trajectories

Figure 6.

Joint acceleration of the optimized exciting trajectories

Figure 7.

Optimized exciting trajectories in Cartesian space

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 2.

Identified parameter values using FFS and MFS

	Parameter	Unit	FFS	MFS
β ₁	I_xx3- I_yy3	Kg.m²	2.3016	2.4006
β ₂	I _xy3	Kg.m²	1.0597	0.3805
β ₃	I _xz3	Kg.m²	−0.1276	−0.0624
β ₄	I _yz3	Kg.m²	−0.2639	0.2106
β ₅	I _zz3	Kg.m²	2.2872	2.6921
β ₆	m ₃ r _x3	Kg.m	2.6275	2.6013
β ₇	m ₃ r _y3	Kg.m	3.9644	3.9539
β ₈	I_xx2- I_yy2-m₃a₃²	Kg.m²	−14.8817	−14.4210
β ₉	I _xy2	Kg.m²	−0.7683	0.6260
β ₁₀	I_xz2+a₃m₃ r_z3	Kg.m²	−2.2521	−1.4897
β ₁₁	I _yz2	Kg.m²	−0.5026	−0.0428
β ₁₂	I_zz2+ m₃a₃²	Kg.m²	13.8072	14.6066
β ₁₃	m₂ r_x2 + m₃a₃	Kg.m	28.4567	28.4915
β ₁₄	m ₂ r_y2	Kg.m	−0.0050	0.0492
β ₁₅	I_zz1+ I_yy2+ I_yy3+m₂(a₂²+ d₂²) +m₃(a₂²+ a₃²+ d₂²) +2d₂(m₂ r_z2 +m₃ r_z3)	Kg.m²	21.9545	21.2018
β ₁₆	A _f1	Nm.s	15.5693	15.4669
β ₁₇	f _v1	Nm.s	5.0485	6.0369
β ₁₈	A _f2	Nm.s	15.3187	12.0470
β ₁₉	f _v2	Nm	17.7341	17.7218
β ₂₀	A _f3	Nm	8.1208	8.1302
β ₂₁	f _v3	Nm	3.0682	3.2264

Figure 8.

Circular motion in Cartesian space

Figure 9.

Comparison of measured and predicted joint torque

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.

Table 3.

Comparison of RMS residuals of torque prediction

		Joint 1	Joint 2	Joint 3
RMS	FFS	11.5657	15.6386	4.0360
(Nm)	MFS	5.5254	10.8921	4.0168

6. 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.

Footnotes

7. Acknowledgments

The research work is jointly supported by the Zhejiang Provincial Natural Science Foundation of China (No : Y1100693), the Fundamental Research Funds for the Central Universities and the open foundation of the State Key Laboratory of Fluid Power Transmission and Control.

References

Brogardh

, “Robot Control Overview: An Industrial Perspective,” Modeling Identification and Control, vol. 30, pp. 167–180, 2009.

Kostic

, “Modeling and identification for high-performance robot control: An RRR-robotic arm case study,” IEEE Transactions on Control Systems Technology, vol. 12, 2004.

Bona

, “Rapid prototyping of a model-based control with friction compensation for a direct-drive robot,” IEEE-ASME Transactions on Mechatronics, vol. 11, pp. 576–584, 2006.

Benimeli

, “A comparison between direct and indirect dynamic parameter identification methods in industrial robots,” Robotica, vol. 24, pp. 579–590, 2006.

Qin

Z. K.

, “A new approach to the dynamic parameter identification of robotic manipulators,” Robotica, vol. 28, pp. 539–547, 2010.

Calanca

, “MIMO Closed Loop Identification of an Industrial Robot,” IEEE Transactions on Control Systems Technology, vol. 19, pp. 1214–1224, 2011.

, “An overview of dynamic parameter identification of robots,” Robotics and Computer-Integrated Manufacturing, vol. 26, pp. 414–419, 2010.

Bingul

Karahan

, “Dynamic identification of Staubli RX-60 robot using PSO and LS methods,” Expert Systems with Applications, vol. 38, pp. 4136–4149, 2011.

Atkeson

C. G.

, “Estimation of inertial parameters of manipulator loads and links,” International Journal of Robotics Research, vol. 5, pp. 101–119, 1986.

10.

Armstrong

, “On finding exciting trajectories for identification experiments involving systems with nonlinear dynamics,” International Journal of Robotics Research, vol. 8, pp. 28–48, 1989.

11.

Swevers

, “Experimental robot identification using optimised periodic trajectories,” Mechanical Systems and Signal Processing, vol. 10, pp. 561–577, 1996.

12.

Olsen

M. M.

Petersen

H. G.

, “A new method for estimating parameters of a dynamic robot model,” IEEE Transactions on Robotics and Automation, vol. 17, pp. 95–100, 2001.

13.

Khosla

P. K.

, “Categorization of parameters in the dynamic robot model,” IEEE Transactions on Robotics and Automation, vol. 5, pp. 261–268, 1989.

14.

Gautier

Khalil

, “Direct calculation of minimum set of inertial parameters of serial robots,” IEEE Transactions on Robotics and Automation, vol. 6, pp. 368–373, 1990.

15.

Lin

S. K.

, “Minimal linear-combinations of the inertia parameters of a manipulator,” IEEE Transactions on Robotics and Automation, vol. 11, pp. 360–373, 1995.

16.

Swevers

, “Optimal robot excitation and identification,” IEEE Transactions on Robotics and Automation, vol. 13, pp. 730–740, 1997.

17.

Antonelli

, “A systematic procedure for the identification of dynamic parameters of robot manipulators,” Robotica, vol. 17, pp. 427–435, 1999.

18.

Capisani

L.M.

, “MIMO identification with optimal experiment design for rigid robot manipulators,” IEEE/ASME International Conference on Advanced Intelligent Mechatronics, vols 1–3, pp.877–882, 2007.