A two-stage calibration method for industrial robots with joint and drive flexibilities

Dealing with robot calibration the neglection of joint and drive flexibilities limit the achievable positioning accuracy significantly. This problem is addressed in this paper. A two stage procedure is presented where elastic deflections are considered for the calculation of the geometric parameters. In the first stage, the unknown stiffness and damping parameters are identified. To this end the model based transfer functions of the linearized system are fitted to captured frequency responses of the real robot. The real frequency responses are determined by exciting the system with periodic multisine signals in the motor torques. In the second stage, the identified elasticity parameters in combination with the measurements of the motor positions are used to compute the real robot pose. On the basis of the estimated pose the geometric calibration is performed and the error between the estimated end-effector position and the real position measured with an external sensor (laser-tracker) is minimized. In the geometric model, joint offsets, axes misalignment, length errors and gear backlash are considered and identified. Experimental results are presented, where a maximum end-effector error (accuracy) of 0.32 mm and for 90 % of the poses a maximum error of 0.23 mm was determined (Stäubli TX90L).


Introduction
One of the main characteristics of industrial robots is their positioning accuracy, strongly depending on the sensor resolution and the geometric parameters.The calculation of the real geometric parameters is called geometric robot calibration and is crucial for accurate robot movements.In literature different calibration methods exist (see Khalil and Dombre, 2004 or Siciliano and Khatib, 2008).However, they are dealing with kinematic models only, neglecting the effects of flexibilities in the joints and drives and thus cause a systematic error in the calculation of the real geometric parameters.This systematic error limits the achievable positioning accuracy of the robot essentially.Only a few publications where flexibilities of the robot in the calibration are considered exist, see Whitney et al. (1986), Gong et al. (2000) or Khalil and Besnard (2002).In the paper Whitney et al. (1986) similar to our contribution a two stage method is proposed.They first identified a compliance model by performing multiple experiments with different external forces and then secondly the geometric calibration is performed by including displace-ments according to the identified model from stage one.The papers of Gong et al. (2000) and Khalil and Besnard (2002) include the elastic displacements in the calculation of the geometric parameters but assumed that the stiffness is known.The goal of this contribution is to present a procedure, which in the first stage identifies the main flexibilities of our robot without external hardware.To this end, the methods presented in Hardeman (2008) and Wernholt (2007) are implemented which both performed frequency domain identifications of industrial robots.In the second stage, the elastic deflections are considered in the geometric calibration leading to reliable geometric parameters.
In Sect. 2 the dynamic model of our 6-axis articulated robot with joint and drive elasticities is derived.Section 3 deals with the modeling of the geometric error parameters, like joint offsets, axes misalignment, length errors and gear backlash.Subsequently, in Sect. 4 the elasticity parameters are identified, using a frequency domain approach as in the work of Hardeman (2008) and Wernholt (2007).The identification of the robots frequency response matrix is presented, Ñ ÒØ× q 1 q 2 q 3 q 4 q 5 q 6 I r E ϕ E I O where the computation of suitable excitation signals are of prime importance.Then the matching between the transfer matrix of the dynamic model and the identified robot frequency response matrix is discussed.Finally in Sect. 5 the computation of the geometric parameters is addressed and the results for an industrial robot Stäubli TX90L are presented.

Dynamic modeling
This section deals with derivation of a dynamic model for the 6-axis industrial robot depicted in Fig. 1.Looking at the mechanical setup, the flexibilities of the first three joints and drives influence the positioning accuracy of the robots endeffector essentially.Hence, a finite bearing and drive stiffness is modeled leading to 4 degrees of freedom (DOF) for each axis i = {1, 2, 3}.The 4 DOF are composed of the motor position q iM , the small bearing distortions q ix and q iy and the arm rotation q i (see Fig. 2).The bearing distortions and the arm rotation are combined in the vector q T i = [q ix , q iy , q i ] representing the arm orientation.The modeling of flexible joints is motivated by the work of Hardeman (2008) which showed that the bearing and drive stiffness is at the same order of magnitude.Elastic effects of the last three axes (the spherical wrist) are of minor importance w.r.t. the endeffector positioning accuracy and are thus modeled without flexibilities.For the derivation of the dynamic model, the transformation from the body fixed frame of the previous arm (p) into the body fixed frame of the current arm (c) is necessary.Hence, the rotation matrix is introduced, where successive rotations about the main axes are performed.The rotation matrix for a rotation about the x, y or z axis is given by respectively.The modeling is based on the paper of Öhr et al. (2006) where consecutive rotations are suggested to model joint and drive compliance.In total, the robot is represented by the vector q T = [q T M , q T A , q T SW ] consisting of the motor coordinates q T M = [q 1M , q 2M , q 3M ] the arm coordinates q T A = [q T 1 , q T 2 , q T 3 ] and the coordinates of the spherical wrist q T SW = [q 4 , q 5 , q 6 ].For the dynamic model, the bearing stiffness is modeled with linear springs and dampers c ix , c iy and d ix , d iy respectively and the drive stiffness and damping with the parameters c i and d i for axis i = {1, 2, 3}.For further details about the dynamical modeling see Öhr et al. (2006).Finally, the dynamic model is derived with the Projection Equation, see Bremer (2008) leading to the equation of motion M(q) q + g q, q, with the mass matrix M(q) ∈ R 15×15 , the motor torques of the first three axis τ M ∈ R 3×1 and the motor torques of the spherical wrist τ SW ∈ R 3×1 .Throughout the whole paper, the identity matrix is represented by I and the zero matrix/vector by 0. The vector g(q, q, p elast ) contains the Coriolis, centrifugal, gravitation and friction forces as well as the stiffness and damping parameters which are going to be identified and are combined in the vector It is apparent from Eq. ( 4), that we are dealing with an underactuated system.

Linearized dynamic model
Basis for the calculation of the transfer matrix is the linearization of the dynamic model.A static equilibrium q 0 ( q0 = 0, q0 = 0) is used as linearization point.linearization, a Taylor series expansion of Eq. ( 4) at the equilibrium configuration q 0 up to order 1 is carried out, leading to (5) with the mass matrix M lin = M q 0 , the damping matrix D lin = ∂g q, q, p elast ∂ q q=q 0 , q=0 , the stiffness matrix K lin = ∂M(q) q + g q, q, p elast ∂q q=q 0 , q=0, q=0 and τ M = τ M0 + τ M and τ SW = τ SW0 + τ SW .The vectors τ M0 and τ SW0 represent the constant motor torques in the static configuration q 0 .Thus, if q 0 is a valid static configuration the term g(q 0 , 0, p elast ) − B M τ M0 − B SW τ SW0 vanishes, and we get the linearized equations of motion

Motor transfer matrix
For the identification of the modeled elasticity parameters, see Sect. 4, the transfer matrix from the motor torques τ M to the motor accelerations qM is necessary (Only the first three axes are included).By using the motor accelerations the double integrator behavior in the transfer matrix is avoided.
Starting with the linearized robot dynamics of Eq. ( 6), the system is reduced to the coordinates q M by using the selection matrix F T M = [I, 0] and the relation q M = F T M q.Using the principle of virtual work, we get the linearized motor equations The motor torques of the spherical wrist τ SW vanish because F T M B SW = 0. Applying the Laplace transformation to Eq. ( 7) we get with a M = qM .The vectors a M and τ M are the Laplace transformed values of a M and τ M , respectively.Thus the transfer matrix leads to where the superscript q 0 denotes the linearization point.

Determination of a static equilibrium
The calculation of a static equilibrium q 0 for our underactuated system, see Eq. ( 4) is neccessary for two different reasons.First, for the correct evaluation of the transfer matrix by using the correct linearization point, and second for the evaluation of the real robot pose and thus for the real endeffector position and orientation in a static configuration.
The motor position q M and the coordinates of the spherical wrist q SW are measured.Hence, only the static values of the arm coordinates q A have to be computed.Therefore the linear system Eq.( 5) is reduced to the coordinates of interest i.e. the arm coordinates q A , by again applying the principal of virtual work and the relation q A = F T A q with F T A = [0 9×3 , I, 0 9×3 ].Assuming a static solution -time derivatives of q are zero -we get Since the arm coordinates are not actuated 10) leads to If q 0 is a static pose of our robot, then g(q 0 , 0, p elast ) = 0 otherwise the displacement q A is calculated by according to Eq. ( 11).Thus the static solution q 0 can be calculated iteratively by with the iteration index i.Because the stiffness of our robot is quite high, the initial solution q (0) 0 = q 1M , q 2M , q 3M , 0, 0, q 1M , 0, 0, q 2M , 0, 0, q 3M , q 4 , q 5 , q 6 T ( 14) is obvious and often one iteration is sufficient to obtain a valid static solution.Having a valid static solution calculated, a geometric model is needed to evaluate the corresponding end-effector position and orientation.

Geometric model
The forward kinematics describes the end-effectors position I r E and orientation ϕ E (see Fig. 1).Furthermore, it depends on the arm and wrist coordinates, q A and q SW respectively, and also on the geometric parameters.In this paper, the geometric parameters are separated into the known nominal values p nom which describe the nominal kinematics of the robot and the unknown geometric error parameters p ge comprising joint offsets, axes misalignment, length errors, and gear backlash.To model a joint offset p 0 for a rotation around the x axis with the DOF q the rotation matrix R leads to Axes misalignment, also for the example of a rotation around the x axis are introduced by adapting the rotation matrix R to with p β and p γ as misalignment angles.To include length errors the connection vector r which describes the links of our robot, is extended by the length error parameters p lx , p ly and p lz leading to Note, the parameters l x , l y and l z are the nominal values of r.In contrast to the previous errors, gear backlash is not a geometric error and depend on the pose of the robot.In order to identify the direction in which the clearance is present, the center of mass of each body must be considered.To avoid the evaluation of the body dynamics the decision is based on the sign of the motor torque τ M in each axis leading to the rotation matrix In the final geometric model of our robot, joint offsets, axes misalignments and length errors combined in the vector p i = [p 0i , p βi , p γ i , p lxi , p lyi , p lzi ] T for all axes i = {1, . . ., 6} as well as for the end-effector tool i = E are included.Gear backlash p BLi is only modeled for the second and third axes i = {2, 3}.The offset of the inertial frame of the robot and the external sensor (laser-tracker) is modeled by the length errors p lxI , p lyI and p lzI .Thus finally, 47 error parameters combined in the vector p ge = [p lxI , p lyI , p lzI , p T 1 , p T 2 , p BL2 , p T 3 , p BL3 , p T 4 , p T 5 , p T 6 , p T E ] T describe the end-effector position I r E (q A , q SW , p nom , p ge ) and orientation ϕ E (q A , q SW , p nom , p ge ) of the robot.Since the endeffector position and orientation depend on the arm coordinates, and the arm coordinates highly depend on the elasticity parameters their identification is treated in the next section.

Identification of elasticity parameters
For the identification of elasticity parameters, first the nonparametric frequency response matrix (FRM) see Pintelon and Schoukens (2012) of the real robot is determined, and second the transfer matrix calculated in Eq. ( 9) is adapted to the real one by adapting the elasticity parameters.

Identification of the non-parametric frequency response matrix
In standard operation mode, the robot is controlled by a feedback PD controller.To improve the identification accuracy, the influence of the PD controller is reduced by using small feedback gains.To ensure that the robot maintains the position, feed forward torque calculated from the inverse dynamics is added.The measurement setup is depicted in Fig. 3.
The system is excited with periodic motor torques while simultaneously the actual motor torque and motor position is measured.The obtained measurement data is then transformed into the frequency domain and used for the FRM calculation.To reduce the effects of stick-slip transitions a sine wave with an amplitude of 3 • serves as reference signal q Md for the position controller.For the excitation of the system, normalized random multisine signals are used which leads to a separation of the FRM in three parts with G BLA (j ω k ) the best linear approximation (BLA), G S (j ω k ) the stochastic nonlinear contributions and N G (j ω k ) the errors due to the output noise (see Pintelon and Schoukens, 2012).We are interested in G BLA (j ω k ) which is used for the identification of the elasticity parameters in Sect.4.2.

Calculation of the FRM
Assuming an appropriate excitation of the system, the measured motor torques and motor positions can be transformed to the frequency domain using the fast Fourier transformation (FFT).Note a signal x is represented in the frequency  domain as x or F{x}.In order to reduce the nonlinear contributions G S , N r different excitation signals and to reduce the noise contribution N G , N p periods are recorded.The FRM for an excitation signal r in pose q is calculated by Ĝ with QM = [ qM,1 qM,2 qM,3 ] and M = [τ M,1 τ M,2 τ M,3 ] resulting from measurement data for 3 orthogonal excitations, see Sect.4.1.2.For different excitation signals, also slightly different FRM results are computed which have to be combined using averaging techniques.The choice of the averaging method depend on the signal to noise ratio (SNR) of the measurement data and the measurement setup.A comparison of different strategies is discussed in the paper of Wernholt and Moberg (2008).For our robot, the simple arithmetic averaging is useful, which is given by Ĝ

Synthesis of the excitation signal
The quality of the FRM highly depends on the excitation signal.For frequency domain identification, periodic signals are very useful because they do not suffer from leakage effects.Furthermore the frequency resolution and power spectrum can be customized for the robot resulting in a minimum measurement time by a maximum quality of the signals (see Pintelon and Schoukens, 2012).Odd random phase multisine signals feature all these properties and are given by consisting of N f different frequencies of odd multiplicity of the basis frequency ω 0 .We selected uniform distributed ran-dom phases φ k and a constant amplitude spectrum A k .Because the FRM of the first three axes should be identified, a excitation signal for each motor has to be generated.Furthermore, for the FRM calculation with Eq. ( 20) the motor torque matrix M must have full rank.To get this property, three sets of orthogonal excitation signals are calculated.Therefore, three different odd random phase multisine signals according to Eq. ( 22) are generated and combined in a diagonal matrix leading to R is generated in the frequency domain and k represents the discrete frequency index.The columns of D [r] R (k) are orthogonalized using the matrix with n u the number of input signals (n u = 3).This yields the excitation matrix Each column of [r]  exc (k) represents one set of excitation signals and is orthogonal to the other ones.Using the matrix T in combination with the multisine signals D

Measurement procedure for the FRM identification
To calculate the FRM of our robot, N r different excitation signals are calculated, each for N p periods.Note, at least the first period of the measurement data can not be used for the FRM calculation because the system has to be in steady state.To reduce nonlinear contributions each excitation signal is applied once with positive and negative sign but with the same reference trajectory q Md .Then the measurement data of these two excitations are subtracted from each other and are used for FRM calculation, as suggested in Hardeman (2008).
Figure 4 shows a graphical representation of the procedure.It starts with the calculation of three independent excitation signals τ i,exc using Eq. ( 22), and it ends with an accurate estimate of the FRM Ĝ[r,q] for one excitation r in pose q.The procedure is performed N r times, leading to N r slightly different FRM results which are in a next step averaged using Eq. ( 21).Our experiments have shown that N r = 4 different realizations each with 4 periods, two periods for the calculation of the transfer matrix and two periods to get in steady state are convenient.The basis frequency of the excitation signal (see Eq. 22) is defined to ω 0 = 2π 0.25 rad s −1 .Hence the signal has a period of 4 s.The excitation of the system for obtaining the transfer matrix for a single realization requires 96 s.This time is composed of the duration of the 3 orthogonal excitations (see Fig. 4) where each consists of a measurement with positive and negative sign with 4 periods.Thus all measurements (4 different realizations) are completed after 384 s.This yields Ĝ[q] which is shown in Fig. 5 for the pose qP 1 of Fig. 6.To get a global representation of the robot, the FRM is identified in 3 different joint configurations.The selected poses are depicted in Fig. 6.These poses represent extrema regarding the moment of inertia in the considered axes.By using poses with extrema in the inertia, the obtained transfer matrices are assumed to be sufficiently different to identify global valid parameters.For all poses the FRM is  identified and used for the elasticity parameter identification of our robot.

Identification of the elasticity parameters
To identify the elasticity parameters, the parametric transfer matrix of our robot -derived in Sect.2.1 and given by Eq. ( 9) -is fitted to the measured one, by changing the parameters p elast .The fitting procedure performs a minimization of the cost functional q∈{qP 1,qP 2,qP 3} k=1...N ω [q] ω k , p elast * W [q] (ω k ) [q] ω k , p elast (26) with the complex vector The minimization is carried out by applying a genetic algorithm followed by a gradient based minimization algorithm.
The vector [q] (ω k , p elast ) represents the complex error between the measured FRM and the transfer matrix of the linearized system see Eq. ( 9).The superscript * denotes the conjugate transpose matrix.To get a good fit in the region of interest, a weighting matrix W [q] (ω k ) is introduced.Especially the error in the vicinity of poles and zeros in the diagonal elements is amplified, leading to reasonable system parameters.
Figure 7 shows both, the non-parametric FRM Ĝ[qP 1] and the optimized parametric transfer function G [qP 1] in the first pose of the robot.It can be seen, that the coupling between the first and the other two axes is quite low and thus the measured and identified transfer functions agree not very well.Furthermore, the weighting matrix was tuned to get a good fit in the diagonal elements of the transfer matrices of poses qP 1, qP 2 and qP 3. ments of the non-parametric and optimized parametric transfer matrices of the poses qP 2 and qP 3. Unfortunately, not all poles and zeros agree exactly, which is mainly attributed to nonlinearities in the stiffness parameters.Overall, a good fit between the non-parametric and parametric transfer functions is achieved and the identified parameters are listed in Table 1.Note, if only one pose is matched the accordance of the identified and measured FRM is very high but the parameters are not reliable anymore.

Geometric calibration
With the identified dynamic model of the last section the real robot pose due to the motor positions can be calculated.This estimated poses in combination with the real end-effector po-sition, measured by an external sensor (laser-tracker), is used to calibrate the geometric model of Sect. 3. Our external sensor can not measure the end-effector orientation and thus the geometric calibration is discussed by using end-effector positions only.

Calculation of the estimated geometric parameters
The estimated end-effector position I r [q] E (q A , q SW , p nom , p ge ) can be evaluated by inserting the calculated arm coordinates q A in Eq. ( 13) and the measured spherical wrist coordinates q SW .Assuming small geometric errors we start with the initial vector of p (0) ge = 0. Thus the end-effector error for a set of N poses follows with I r [q] E,m the end-effector position measured by the lasertracker in pose q of the robot.Performing a Taylor series expansion on e geo at p ge = p .The superscript (i) indicates the current set of parameters.Unfortunately, not all geometric error parameters are independent, i.e. rank of is less than the number of parameters p ge .These linear dependencies must be eliminated.Therefore a numerical regularization algorithm which is based on a QR-decomposition is applied (see Gautier, 1991).30 independent parameters out of the 47 modeled error parameters are found.The parameters p ge are calculated by minimizing the error e geo , using e.g. the Levenberg-Marquardt algorithm.The quality of the calibration result highly depends on the poses.

Calculation of optimal poses for the calibration
In order to get a good excitation of the geometric parameters the calibration poses are very important.In the paper of Sun and Hollerbach (2008a) the choice of the observability index for the pose calculation and in Sun and Hollerbach (2008b) and Zhuang (1994) selection algorithm for the computation of optimal poses are presented.We used the algorithm presented in Sun and Hollerbach (2008b) where the minimum singular value of the covariance matrix = T is maximized.Finally, 40 optimal poses are calculated with a minimum singular value of σ min = 0.46.Exemplarily, 3 out of the 40 poses are depicted in Fig. 10.For all optimal poses the end-effector is measured with a laser tracker and the actual motor positions are measured by the motor encoders.Then, the estimated arm angles q A are evaluated using Eq. ( 13) and the error of Eq. ( 28) is minimized using the Levenberg-Marquardt algorithm.The absolute positioning accuracy for the uncalibrated and calibrated case is depicted in Fig. 11.The positioning accuracy for the calibrated case is shown in Fig. 12.For the calibration poses, a positioning accuracy of about 0.1 mm is obtained which is in the range of the repeatability of our robot, and thus a very good result.

Evaluation of the absolute positioning accuracy
With the geometric parameters p ge identified in the last section, the positioning accuracy in the whole workspace should be evaluated.Therefore 150 random poses in the whole workspace of the robot are generated and measurements with the laser-tracker.The calculated end-effector position due to the estimated robot configuration is compared to the measured position and the obtained absolute positioning accuracy is shown in Fig. 13.For 150 arbitrary poses an absolute error less than 0.32 mm and in 90 % of the poses an error less than 0.23 mm is obtained.If the elastic deflections are not considered in the calibration, the maximum error for the same poses is 0.7 mm which thus justifies the effort of the presented procedure.

Conclusions
This paper presents the identification of geometric parameters by including stiffness and damping parameters in the calibration process.Not only the inclusion of the elasticity parameters also their identification was discussed.The identification was carried out in the frequency domain, where the transfer matrix of the linearized system was adapted to the real robots transfer matrix by changing the elasticity parameters.Regarding the calculation of the real robots transfer matrix, the excitation signal is crucial and thus it is discussed in detail.The identification of globally valid parameters requires different poses.In this paper, the poses were selected manually.The calculation of optimal poses for the elasticity parameter identification is part of future work.With the identified elasticity parameters, the real robots configuration is estimated using the system dynamics and the geometric calibration is performed.Experimental results for a Stäubli TX90L are reported.Thereupon geometric model parameters were identified which lead to a positioning accuracy of 0.32 mm in the whole workspace.

Figure 2 .
Figure 2. Model of the third axis.

Figure 3 .
Figure 3. Block diagram of the control concept.
arbitrary number of different but orthogonal excitation signals can be calculated.

Figure 4 .
Figure 4. Measurement schedule for the FRM calculation of one realization.

Figure 6 .
Figure 6.Poses for the FRM calculation represented by qP 1, qP 2 and qP 3 from left to right respectively.
e geo,lin = p ge + e geo | p (i