Geometrical Derivation of Differential Kinematics to Calibrate Model Parameters of Flexible Manipulator

Abstract

In this paper differential kinematics was geometrically derived to be utilized in a calibration algorithm that improves the accuracy of the manipulation of a robot. Even though the mechanical components are manufactured and assembled precisely, small differences between the designed and the actual system always exist, due to both geometric and unmodelled errors. In order to resolve these problems, differential relationships between the model parameters and the end-effector's posture were formulated. Subsequently, a derivative based estimation algorithm, such as an EKF (Extended Kalman Filter) manner, could be adopted to update the model parameters. The proposed algorithm includes joint flexibility, so is an advanced version of previous work, where a rigid joint model was adopted [1]. The effectiveness of the proposed algorithm was verified by a computer simulation with a 6 DOF manipulator robot.

1. Introduction

Obtaining information concerning the links and the joints of a manipulator is essential for the control of robotic systems. Especially in research involving an articulated link manipulator, which consists of a series of links, such as a robotic arm, finger, or leg and where the accuracy of the model directly influences the performance of the task. Since the model parameters determine the mathematical model, it is important to gain parameters that are as accurate as possible. Although 3D CAD tools provide accurate information about a model, the geometry of the manufactured robot is not generally identical to the designed one. Consequently, the differences reduce the performance of the manipulation significantly. These inaccuracies in the model parameters of the robotic systems arise from two kinds of errors: geometrical and unmodelled errors. The geometrical errors are caused by inaccurate manufacturing, assembling, etc. The unmodelled errors are commonly caused by the approximation of the mathematical formulation, such as the joint flexibilities, thermal dynamics, etc. Therefore, in order to improve the performance and/or the stability of the manipulation, it is necessary to calibrate the model parameters as well as utilizing the advanced model.

There has been a lot of research related to the kinematical identification and calibration of manipulator robots. A methodology that uses lasers was introduced to capture robot position data in order to model the stiffness of the robot manipulator [2] and to predict kinematic parameters [3, 4, 5, 6, 7]. O'Brian et al. also employed magnetic motion capture data to determine the kinematic parameters [8]. Renaud et al. [9] and Rauf et al. [10] used a vision-based measuring device and a partial pose measurement device, respectively, for the kinematic calibration of parallel mechanisms. Gatti and Danieli demonstrated a pose-matching procedure for calibration, which provided a low cost and easy-to-use external metrology system [11]. Santolaria et al. recently demonstrated a continuous data capture technique by using a ball bar gauge and a coupling probe to estimate the kinematic parameters of articulated arm coordinated measuring machines [12]. It was concluded that the parameter error was minimized in the measured positions, whereas the error increased in very different positions because the optimization technique was only based on the position information of the end-effector. To overcome this problem, Ye et al. [13] and Iurascu and Park [14] demonstrated a kinematic calibration method using differential kinematics and iterative algorithms to determine the parameters.

In the previous research, a novel algorithm was proposed to estimate the entire kinematic parameters of the robot manipulator by using a structured laser module (SLM), a stationary camera and a screen [1]. In the algorithm, by using the Jacobian matrix, which represents the relationship between the kinematical parameters and the laser points on the screen, an iterative estimation scheme, such as the extended Kalman filter (EKF), was adopted. For the kinematic model, Denavit-Hartenberg (D-H) convection is used and all the joints and the links were modelled as rigid. In this paper, the previous method was extended by employing the flexible joint model. To improve the accuracy of the manipulation, three types of model parameters were estimated: the D-H parameters, the joint compliances and the CM positions of the links. Similarly to previous work, an SLM system with a stationary camera and a screen was adopted to estimate the model parameters. The performance of the proposed algorithm was verified with a 6 DOF humanoid leg system by computer simulation.

This paper is organized as follows. In Section 2, problems concerning a flexible joint model are defined and described. Subsequently, in Section 3 and 4, mathematical models and Jacobians are discussed. In Section 5, the effectiveness of the proposed algorithm is verified by simulation results. Lastly, concluding remarks follow in Section 6

2. Problem Definition

In order to control a robot manipulator precisely, it is essential to know the model parameters of the manipulator. However, there are always some errors in link lengths and link twists during manufacturing and/or assembling the robot system. These undesirable errors reduce accuracy during the manipulation and sometimes make the system unstable. They may cause more problems in a complex articulated linkage system such as a humanoid robot. Small errors in each joint accumulate and eventually make a significant difference between the model and the actual system.

In previous work, a novel systematic algorithm was introduced to measure the errors in the kinematical parameters by using a structured laser module (SLM) and a stationary camera system. Using the algorithm, the kinematical errors were successfully measured; thereby it enhanced the accuracy of the forward and inverse kinematics. Since the algorithm uses a lot of standstill manipulator postures in the rigid kinematical model, when the mass of the link of the manipulator is large or the stiffness of the joint is small, differences between the computed and the measured points increase significantly because of the influence of gravitational force. These are represented in Fig. 1, where θ and φ represent the controlled joint angle in the rigid model and the additional rotation in the flexible model, respectively and τ_g represents gravitational torque. Note that in the rigid model φ is assumed to be zero, while it has a significant value in the flexible model.

Figure 1.

Influence of the joint flexibility.

In order to estimate the model parameters, a similar measurement system is introduced as shown in Fig. 2. The proposed system consists of a SLM attached to a manipulator, a stationary camera and a screen. Note that the base frame is attached to the screen and the stationary camera is used to measure the position of the laser point on the screen. When all the kinematical parameters, including joint compliances and mass positions, are accurately known the computed laser points on the screen for the given robot posture must be identical to the actual points measured by the stationary camera. Since there are some model errors, this generally does not happen as mentioned above. The proposed method uses the differences between the computed laser points with the current estimated parameters and the actual measured points to update the parameters. While doing this, the differential relationships, known as the Jacobian matrix, are utilized. A detailed mathematical model and derivation of the Jacobians are described in Section 3 and 4. Note that this paper focuses on how to derive Jacobians geometrically. Its applications, including experimental analysis, will be treated in further work.

Figure 2.

System configuration to calibrate the model parameters.

3. Mathematical Model

D-H convention is adopted to describe the kinematics because it requires a minimum number of parameters to represent the coordinate system attached to each link, as shown in Fig. 3. Here α and β are the joint offset angle and the twist angle between the joint axes and d and l are joint distance and link length, respectively. Lastly, θ and φ are the controlled joint angle and the deformed twisting angle due to the flexibility.

Figure 3.

D-H convention for coordinate frame system.

From the D-H parameters and the twisting angles caused by the joint flexibility, the end-effector's posture is determined as follows:

$\begin{array}{ll}{\pi }_e\hfill & =f(\lambda \mathrm{,}\phi )\hfill \\ \hfill & =\left[\begin{array}{c}{\mathrm{r}}_e\\ {\delta }_e\end{array}\right]\mathrm{.}\hfill \end{array}$

(1)

where λ and φ are 4(n+1) × 1 vector for the D-H parameters and n × 1 vector for the twisting angles, respectively and ∂_e is a 6 × 1 vector that consists of a position (r_e) and an orientation (δ _e ) of the end-effector in 3D space. By differentiating Eq. (1), the following differential kinematic motion is obtained:

$d{\pi }_e=\frac{\partial f}{\partial \lambda }d\lambda +\frac{\partial f}{\partial \phi }d\phi \mathrm{.}$

(2)

In a stationary state, the twisting angles can be represented by the following implicit function:

$\phi =h(\chi \mathrm{,}\lambda \mathrm{,}{\mathrm{p}}_{\mathrm{m}}\mathrm{,}\phi )$

(3)

where χ and p_m represent the joint compliance vector (n × 1) and the CM position vector (3n × 1), respectively. If the twisting angles caused by the joint flexibility are small enough, joint deformation can be approximated as a linear model such as Hooke's law. From this, the twisting angles are simply represented as follows:

$\phi =D_{\chi }\gamma (\lambda \mathrm{,}{\mathrm{p}}_{\mathrm{m}}\mathrm{,}\phi )$

(4a)

$=D_{\gamma }\chi \mathrm{,}$

(4b)

where γ and D_x are the joint torque vector (n × 1) and a matrix diagonalization for vector x, respectively. By differentiating Eq. (4), the differential of the twisting angles is derived as follows:

$\begin{array}{ll}d\phi \hfill & =\frac{\partial h}{\partial \chi }d\chi +\frac{\partial h}{\partial \lambda }d\lambda +\frac{\partial h}{\partial {\mathrm{p}}_{\mathrm{m}}}d{\mathrm{p}}_{\mathrm{m}}+\frac{\partial h}{\partial \phi }d\phi \hfill \\ \hfill & =D_{\gamma }d\chi +D_{\chi }\frac{\partial \gamma }{\partial \lambda }d\lambda +D_{\chi }\frac{\partial \gamma }{\partial {\mathrm{p}}_{\mathrm{m}}}d{\mathrm{p}}_{\mathrm{m}}+D_{\chi }\frac{\partial \gamma }{\partial \phi }d\phi \mathrm{.}\hfill \end{array}$

(5)

By solving Eq. (5) algebraically for dφ (differential relationship among the compliances), the D-H parameters and the mass positions are obtained as follows:

$d\phi ={\left(I-D_{\chi }\frac{\partial \gamma }{\partial \phi }\right)}^{-1}\left(D_{\gamma }d\chi +D_{\chi }\frac{\partial \gamma }{\partial \lambda }d\lambda +D_{\chi }\frac{\partial \gamma }{\partial {\mathrm{p}}_{\mathrm{m}}}d{\mathrm{p}}_{\mathrm{m}}\right)\mathrm{.}$

(6)

By substituting Eq. (6) into Eq. (2), differential kinematic motion is finally expressed with differentials of the D-H parameters. The compliances and the mass positions are as follows:

$d{\pi }_e=\frac{\partial {\pi }_e}{\partial \lambda }d\lambda +\frac{\partial {\pi }_e}{\partial \chi }d\chi +\frac{\partial {\pi }_e}{\partial {\mathrm{p}}_{\mathrm{m}}}d{\mathrm{p}}_{\mathrm{m}}$

(7)

where

$\begin{array}{ll}\frac{\partial {\pi }_e}{\partial \lambda }\hfill & =\frac{\partial f}{\partial \lambda }+\frac{\partial f}{\partial \phi }{\left(D_{\mathrm{k}}-\frac{\partial \gamma }{\partial \phi }\right)}^{-1}\frac{\partial \gamma }{\partial \lambda }\mathrm{,}\hfill \\ \frac{\partial {\pi }_e}{\partial \chi }\hfill & =\frac{\partial f}{\partial \phi }{\left(D_{\mathrm{k}}-\frac{\partial \gamma }{\partial \phi }\right)}^{-1}{D}_{\mathrm{k}}{D}_{\gamma }\mathrm{,}\hfill \\ \frac{\partial {\pi }_e}{\partial {\mathrm{p}}_{\mathrm{m}}}\hfill & =\frac{\partial f}{\partial \phi }{\left(D_{\mathrm{k}}-\frac{\partial \gamma }{\partial \phi }\right)}^{-1}\frac{\partial \gamma }{\partial {\mathrm{p}}_{\mathrm{m}}}\mathrm{.}\hfill \end{array}$

Eq. (7) represents the relationship of the infinitesimal displacements among the kinematical parameters, the compliances and the mass positions for the end-effect posture. By applying a derivative based estimation algorithm such as Least-squares manner or Extended Kalman Filter, the model parameters can be effectively calibrated. The overall algorithm is summarized in the block diagram in Fig. 4. In forward kinematics, the coordinate frames ∊ and the laser points on the screen p_s are calculated with arbitrary selected test postures, i.e., joint variables θ and the initial guesses for the model parameters. While doing this, the twisting angles Φ are iteratively computed until the joint torques due to gravity become identical to the restoring torques due to joint stiffness.

Figure 4.

Block diagram of overall algorithm flow.

Algorithm 1 describes how to compute twisting angles iteratively. After initialization, Homogeneous Transformation Matrices are computed with the given joint coordinates, the model parameters and the k^th twisting angles. Subsequently, from the kinematical configuration, joint torques can be obtained by static equilibrium condition. Since the joint model in this paper assumes that the compliance coefficient is constant, the twisting angles are easily calculated and updated. If the norm of the difference between the kt^h and the (k + 1)^th twisting angles is small enough, the (k + 1)^th twisting angles are saved and the iteration is terminated. After the forward kinematics stage, the Jacobian matrices are computed. Note that the derivations of the Jacobians in Eq. (7) will be described in detail in Section 4. Subsequently, the Least-squares method or EKF is utilized to update the model parameters, the D-H parameters λ, the compliances χ and the CM positions P_m, until the error norm becomes lower than the bounded value or the number of iterations reaches the maximum value.

Algorithm 1.

Compute twisting angles under joint flexibility

/* initialize the twisting angles */ φ⁰ ← 0; k ← 0; while k < n do /* compute Homogeneous Transformation Matrices */ /* compute joint torques */ τ^k←τ({Tⁱ}); /* update the twisting angles */ φ^k+1←D_χτ^k; /* check convergence of the twisting angles */ If ‖ φ^k+1 – φ^k ‖≤ ∊ then break; end if φ^k ← ^k+1; end while /* save the twisting angles */ φ ← φ^k+1;

4. Geometrical Derivation of Jacobians

4.1 Jacobian for the End-Effector Posture from the D-H Parameters and the Twisting Angles.

From Eq. (7), the Jacobian matrix for the end-effector posture from the model parameters consists of partial derivatives as follows:

${\mathrm{J}}_e=\left[\begin{array}{ccc}\frac{\partial {\pi }_e}{\partial \lambda }& \frac{\partial {\pi }_e}{\partial \chi }& \frac{\partial {\pi }_e}{\partial {\mathrm{p}}_{\mathrm{m}}}\end{array}\right]$

(8)

Since the Jacobian matrix accompanies partial derivatives, it is generally hard to compute analytically. Fortunately, however, for the articulated linkage system, it can be computed geometrically. The Jacobian matrix is computed using exactly the same manner as the general kinematics, except that the parameters are the D-H parameters λ = [··· α ^k d^k l^k β^k ···]^T and the twisting angles φ = […φ^k…]^T, not the joint angles θ = [···θ^k···]^T, as shown in Fig 5.

Figure 5.

Coordinate frames for Jacobian Computing

The Jacobian is comprised by differential translations and rotations with respect to the corresponding parameters. The differential translations and rotations by translational parameters such as d^k and l^k are simply unit vectors in the direction of the translations and zero vectors. Also, the differential translations and rotations by the rotational parameters, such as α^k, β^k and φ^k, are a cross product of rotational axis vectors and distance vectors, respectively. Eq. (9) represents these formulas.

$\left[\begin{array}{ccccc}\frac{\partial {\pi }_e}{\partial {\alpha }^{\mathrm{k}}}& \frac{\partial {\pi }_e}{\partial d^{\mathrm{k}}}& \frac{\partial {\pi }_e}{\partial {\mathrm{l}}^{\mathrm{k}}}& \frac{\partial {\pi }_e}{\partial {\beta }^{\mathrm{k}}}& \frac{\partial {\pi }_e}{\partial {\phi }^{\mathrm{k}}}\end{array}\right]=\left[\begin{array}{ccccc}{\widehat{\mathrm{z}}}^{\mathrm{k}-1}\times ({\mathrm{r}}_e-{\mathrm{r}}_{\sigma }^{\mathrm{k}-1})& {\widehat{\mathrm{z}}}^{\mathrm{k}-1}& {\widehat{\mathrm{x}}}^{\mathrm{k}}& {\widehat{\mathrm{x}}}^{\mathrm{k}}\times ({\mathrm{r}}_e-{\mathrm{r}}_{\sigma }^{\mathrm{k}})& {\widehat{\mathrm{z}}}^{\mathrm{k}-1}\times ({\mathrm{r}}_e-{\mathrm{r}}_{\sigma }^{\mathrm{k}-1})\\ {\widehat{\mathrm{z}}}^{\mathrm{k}-1}& 0& 0& {\widehat{\mathrm{x}}}^{\mathrm{k}}& {\widehat{\mathrm{z}}}^{\mathrm{k}-1}\end{array}\right].$

(9)

$\frac{\partial f}{\partial \lambda }=\left[\underset{0}{\underbrace{\begin{array}{cccc}{\widehat{\mathrm{z}}}^G\times ({\mathrm{r}}_e-{\mathrm{r}}_{\sigma }^G)& {\widehat{\mathrm{z}}}^G& {\widehat{\mathrm{x}}}^0& {\widehat{\mathrm{x}}}^0\times ({\mathrm{r}}_e-{\mathrm{r}}_{\sigma }^0)\\ {\widehat{\mathrm{z}}}^G& 0& 0& {\widehat{\mathrm{x}}}^0\end{array}}}\cdots \underset{\mathrm{k}}{\underbrace{\begin{array}{cccc}{\widehat{\mathrm{z}}}^{\mathrm{k}-1}\times ({\mathrm{r}}_e-{\mathrm{r}}_{\sigma }^{\mathrm{k}-1})& {\widehat{\mathrm{z}}}^{\mathrm{k}-1}& {\widehat{\mathrm{x}}}^{\mathrm{k}}& {\widehat{\mathrm{x}}}^{\mathrm{k}}\times ({\mathrm{r}}_e-{\mathrm{r}}_{\sigma }^{\mathrm{k}})\\ {\widehat{\mathrm{z}}}^{\mathrm{k}-1}& 0& 0& {\widehat{\mathrm{x}}}^{\mathrm{k}}\end{array}}}\cdots \right]\text{for}1\le \mathrm{k}\le n\mathrm{.}$

(10)

4.1.1 Partial Derivatives for the D-H parameters

By applying Eq. (9), the partial derivatives of the end-effector posture for the D-H parameters are computed as Eq. (10).

4.1.2 Partial Derivatives for the twisting angles

In a similar manner, the partial derivatives of the end-effector posture for the twisting angles are computed as follows:

$\begin{array}{l}\frac{\partial f}{\partial \phi }=[\cdots \underset{\mathrm{k}}{\underbrace{\begin{array}{c}{\widehat{\mathrm{z}}}^{\mathrm{k}-1}\times ({\mathrm{r}}_e-{\mathrm{r}}_{\sigma }^{\mathrm{k}-1})\\ {\widehat{\mathrm{z}}}^{\mathrm{k}-1}\end{array}}}\cdots ]\text{for}1\le \mathrm{k}\le n\mathrm{.}\hfill \end{array}$

(11)

4.2 Derivatives of the Joint Torques for the D-H Parameters, the Twisting Angles and the CM Positions

The laser pointers on the screen are measured during standstill state. Consequently, the torques applied to each joint are determined from the gravitational force only. Fig. 6 illustrates the coordinate frames and the CM positions for the revolute joints. Note that the k^th joint torque, τ^k, is affected from the k^th link to the last link.

Figure 6.

Coordinate frames and CM positions for revolute joints

The effective CM position, which is the average value from the k^th link to the last link, is calculated as follows:

$\begin{array}{ll}{\bar{\mathrm{r}}}_{\mathrm{m}}^{\mathrm{k}}\hfill & =\frac{{\mathrm{m}}^{\mathrm{k}}{\mathrm{r}}_{\mathrm{m}}^{\mathrm{k}}+{\mathrm{m}}^{\mathrm{k}+1}{\mathrm{r}}_{\mathrm{m}}^{\mathrm{k}+1}+\cdots +{\mathrm{m}}^n{\mathrm{r}}_{\mathrm{m}}^n}{{\mathrm{m}}^{\mathrm{k}}+{\mathrm{m}}^{\mathrm{k}+1}+\cdots +{\mathrm{m}}^n}\hfill \\ \hfill & =\frac{1}{{\mathrm{M}}^{\mathrm{k}}}{\displaystyle \sum _{j=\mathrm{k}}^n}{\mathrm{m}}^j{\mathrm{r}}_{\mathrm{m}}^j\hfill \end{array}$

(12)

where r_m^k and m^k represent the CM position vector and the mass for the k^th link, respectively. Subsequently, by using triple product formula, the joint torque is computed as follows:

${\tau }^{\mathrm{k}}=\left(f_{\mathrm{a}}^{\mathrm{k}}\times f_{\mathrm{b}}^{\mathrm{k}}\right)\cdot f_{\mathrm{c}}^{\mathrm{k}}$

(13)

with

$f_{\mathrm{a}}^{\mathrm{k}}={\bar{\mathrm{r}}}_{\mathrm{m}}^{\mathrm{k}}-{\mathrm{r}}_{\sigma }^{\mathrm{k}-1}\mathrm{,}$

(14a)

$f_{\mathrm{b}}^{\mathrm{k}}={\mathrm{M}}^{\mathrm{k}}g\mathrm{,}$

(14b)

$f_{\mathrm{c}}^{\mathrm{k}}={\widehat{\mathrm{z}}}^{\mathrm{k}-1}$

(14c)

where f_a^k is the distance vector from the k^th joint to the effective CM position and from the k^th link to the last link and f_b^k is the gravitational force with the total mass after the k^th link. Lastly, f_c^k is the joint unit vector.

By differentiating Eq. (13), the differential of the joint torque for the D-H parameters, the CM positions and the twisting angles are derived as follows:

$\mathrm{d}{\tau }^{\mathrm{k}}=\frac{\partial {\tau }^{\mathrm{k}}}{\partial \lambda }\mathrm{d}\lambda +\frac{\partial {\tau }^{\mathrm{k}}}{\partial {\mathrm{p}}_{\mathrm{m}}}\mathrm{d}{\mathrm{p}}_{\mathrm{m}}+\frac{\partial {\tau }^{\mathrm{k}}}{\partial \phi }\mathrm{d}\phi$

(15)

Where detailed partial derivatives are given as follows:

$\frac{\partial {\tau }^{\mathrm{k}}}{\partial \lambda }={\left(f_{\mathrm{a}}^{\mathrm{k}}\times f_{\mathrm{b}}^{\mathrm{k}}\right)}^T\frac{\partial f_{\mathrm{c}}^{\mathrm{k}}}{\partial \lambda }+{\left(f_{\mathrm{b}}^{\mathrm{k}}\times f_{\mathrm{c}}^{\mathrm{k}}\right)}^T\frac{\partial f_{\mathrm{a}}^{\mathrm{k}}}{\partial \lambda }\mathrm{,}$

(16a)

$\frac{\partial {\tau }^{\mathrm{k}}}{\partial {\mathrm{p}}_{\mathrm{m}}}={\left(f_{\mathrm{b}}^{\mathrm{k}}\times f_{\mathrm{c}}^{\mathrm{k}}\right)}^T\frac{\partial f_{\mathrm{a}}^{\mathrm{k}}}{\partial {\mathrm{p}}_{\mathrm{m}}}\mathrm{,}$

(16b)

$\frac{\partial {\tau }^{\mathrm{k}}}{\partial \phi }={\left(f_{\mathrm{a}}^{\mathrm{k}}\times f_{\mathrm{b}}^{\mathrm{k}}\right)}^T\frac{\partial f_{\mathrm{c}}^{\mathrm{k}}}{\partial \phi }+{\left(f_{\mathrm{b}}^{\mathrm{k}}\times f_{\mathrm{c}}^{\mathrm{k}}\right)}^T\frac{\partial f_{\mathrm{a}}^{\mathrm{k}}}{\partial \phi }$

(16c)

From Eq. (11), (12), and (16), the Jacobian matrix, Eq. (8), can be computed sequentially.

4.3 Comparison with numerical computation method

So far, the Jacobian matrix between the model parameters and the end-effector posture is geometrically derived. An alternative approach for the Jacobian derivation is to use numerical differentiation. Using the definition of the Jacobian matrix, it can be divided into two parts: position derivative and orientation derivative, as follows:

$\mathrm{J}={\left[\begin{array}{cc}\frac{\partial {\mathrm{r}}_e}{\partial \mathrm{x}}& \frac{\partial {\delta }_e}{\partial \mathrm{x}}\end{array}\right]}^T$

(17)

where x = [λ^T χ^T p_m^T]^T is a vector for the model parameters. The position derivative is obtained by computing the displacement with small variation as follows:

$\frac{\partial {\mathrm{r}}_e}{\partial {\mathrm{x}}^i}=\frac{{\mathrm{r}}_e({\mathrm{x}}^i+\mathrm{\Delta })-{\mathrm{r}}_e({\mathrm{x}}^i-\mathrm{\Delta })}{2\mathrm{\Delta }}$

(18)

where xⁱ is the i^th value of x. The orientation derivative is slightly complicated compared with the position case. From the properties of the rotation matrix, rotational displacement can also be expressed as a rotation matrix.

${\stackrel{\sim }{R}}^i=R{({\mathrm{x}}^i-\mathrm{\Delta })}^{T}R({\mathrm{x}}^i+\mathrm{\Delta })$

(19)

where R˜ⁱ represents the rotational displacement between successive rotations with small deviations in the xⁱ. By converting the angle-axis representation, displacement angle, ξ, is derived as follows:

$\xi ={{\displaystyle \mathrm{c}o\mathrm{s}}}^{-1}\left(\frac{t\mathrm{r}({\stackrel{\sim }{R}}^i)-1}{2}\right);$

(20)

where tr(·) represents a trace operation. The axis vector is derived as follows:

$\mathrm{a}=\frac{1}{2\text{sin}(\xi )}\left[\begin{array}{c}{\stackrel{\sim }{\mathrm{R}}}_{\mathrm{z}\mathrm{y}}-{\stackrel{\sim }{\mathrm{R}}}_{\mathrm{y}\mathrm{z}}\\ {\stackrel{\sim }{\mathrm{R}}}_{\mathrm{x}\mathrm{z}}-{\stackrel{\sim }{\mathrm{R}}}_{\mathrm{z}\mathrm{x}}\\ {\stackrel{\sim }{\mathrm{R}}}_{\mathrm{y}\mathrm{x}}-{\stackrel{\sim }{\mathrm{R}}}_{\mathrm{x}\mathrm{y}}\end{array}\right]$

(21)

From Eq. (20) and (21), the orientation derivative is numerically derived as follows:

$\frac{\partial {\delta }_e}{\partial {\mathrm{x}}^i}=\frac{\xi R({\mathrm{x}}^i-\mathrm{\Delta })\mathrm{a}}{2\mathrm{\Delta }}$

(22)

In Table 1, the accuracy of the numerically derived Jacobian is compared with the geometrically derived Jacobian. For quantitative analysis, 1000 test postures were selected and Jacobian matrices were computed according to the above two methodologies. Subsequently, the results were analysed in a statistical manner. From the table, it can be shown that both approaches give similar results. This implicitly guarantees that the geometrical formula is correctly derived and thus it gives true values. Practically, numerical derivation sometimes suffers from the division by zero problem.

Table 1.

Errors of the numerical Jacobians

	Mean	STD	Max	Min
Error (×10−4)	9.3489	7.7841	87.122	1.4477

Another important benefit of using geometrical derivation is computing time. Since the numerical derivation utilizes two different postures deviated by Δ to obtain one column of Jacobian matrix, ∂∂_e / ∂xⁱ, it requires a much greater computational process. From the simulation results, it can be seen that it takes 10.0622 seconds for 1000 calculations of geometrical derivation and 249.005 seconds for 1000 calculations of numerical derivation using Matlab. Using the geometrical derivation of a Jacobian matrix reduces the run time for the calibration algorithm dramatically.

4.4 Jacobian for the laser points on the screen from the end-effector

After obtaining the relationships from the kinematical parameters, including the compliances and the CM positions of the end-effector posture, the relationship between the end-effector posture and the laser points projected onto the screen must be found. From [1], the Jacobian matrix of the laser point on the screen from the end-effector is given by:

$d{\mathrm{p}}_{\mathrm{s}}={\mathrm{J}}_{\mathrm{s}}d{\pi }_{e\mathrm{l}}$

(23)

where

${\mathrm{J}}_{\mathrm{s}}=\left[\begin{array}{cccccc}-\frac{{\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{y}}}{{\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{x}}}& 1& 0& \frac{{\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{z}}}{{\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{x}}}{\mathrm{p}}_{\mathrm{l}\mathrm{x}}& \frac{{\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{y}}{\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{z}}}{{({\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{x}})}^2}{\mathrm{p}}_{\mathrm{l}\mathrm{x}}& \frac{{({\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{z}})}^2-1}{{({\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{x}})}^2}{\mathrm{p}}_{\mathrm{l}\mathrm{x}}\\ -\frac{{\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{z}}}{{\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{x}}}& 0& 1& -\frac{{\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{y}}}{{\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{x}}}{\mathrm{p}}_{\mathrm{l}\mathrm{x}}& \frac{1-{({\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{y}})}^2}{{({\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{x}})}^2}{\mathrm{p}}_{\mathrm{l}\mathrm{x}}& -\frac{{\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{y}}{\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{z}}}{{({\mathrm{r}}_{\mathrm{l}\mathrm{z}\mathrm{x}})}^2}{\mathrm{p}}_{\mathrm{l}\mathrm{x}}\end{array}\right]\mathrm{.}(36)$

From Eq. (8) and (23), the final Jacobian matrix is derived as follows:

$\mathrm{J}={\mathrm{J}}_{\mathrm{s}}{\mathrm{J}}_e$

(24)

Note that the Jacobian matrix is the relationship of the differential motions between the model parameters and the laser points on the screen.

5. Computer Simulation

The effectiveness of the derived Jacobians was verified by adopting a computer simulation with a 6 DOF humanoid-legged robot, as shown in Fig. 7. To calibrate the model parameters, 100 test postures were selected by arbitrarily choosing joint angles between −20 and +20 degrees. Subsequently, the laser pointers on the screen corresponding to each posture were obtained. An initial guess for the D-H parameters is randomly selected between 80% and 120% of the true value and random noise between −0.02 to +0.02 is added. The rest of the model parameters for the initial guess are randomly selected between 90% and 110% of the true value and random noise is added between −0.01 and +0.01. Also, during the measurement of the laser points on the screen, random noise between −0.001 to +0.001 was added. Note that the coordinate frame for the 1^st laser in SLM, σ_l1, coincides with the coordinate frame for the 6^th link, σ₆, while the coordinate frames for the 2^nd and the 3^rd lasers, σ₁₂ and σ₁₃, are rotated by +30 and −30 degrees about the x-axis, respectively and are followed by a rotation of +10 degrees about the y-axis. Detailed model parameters used in the simulation are summarized in Table 2.

Table 2.

Model parameters (true values).

link	D-H parameters				compliance (°/kg)	CM position (m)			mass (kg)
	α(°)	d (m)	l (m)	ß (°)		X	y	z
0	180.00	1.00	1.00	0.00	−	−	−	−	−
1	180.00	0.00	0.00	90.00	5.00	0.00	0.00	0.00	1.00
2	−90.00	0.00	0.00	−90.00	5.00	0.00	0.00	0.00	1.00
3	0.00	0.00	0.30	0.00	5.00	0.00	0.00	0.00	1.00
4	0.00	0.00	0.30	0.00	5.00	−0.15	0.00	0.00	1.00
5	0.00	0.00	0.00	90.00	5.00	0.00	0.00	0.00	1.00
6	0.00	0.10	0.00	0.00	5.00	0.00	0.00	0.00	1.00

Figure 7.

Coordinate frames of a 6 DOF leg robot.

In order to verify the effectiveness of the flexible joint model, the rigid joint model (previous work) was also examined. For EKF, both the covariance matrix of system noise and measurement noise were set as 1.0×10⁴ for D-H parameters and 1.0×10⁶ for the rest of the model parameters. The simulation was carried out using the Matlab program on Intel core i5–2500 CPU at 3.30 GHz with 4.0 GB RAM. The simulation takes 10.6 seconds for the rigid model and 146.5 seconds for the flexible model, respectively. Note that the total number of model parameters to be calibrated is 28 (excluding compliances and CM positions) for the rigid model and 52 for the flexible model.

5.1 Simulation with rigid joint model

In the rigid joint model, the Jacobian matrix of the end-effector's posture is a function of D-H parameters only. Consequently, the differential kinematic motion in Eq. (7) is simplified.

$d{\pi }_e=\frac{\partial {\pi }_e}{\partial \lambda }d\lambda$

(25)

From Eq. (25), EKF was iteratively executed with the same manner as in the previous work [1], except that compliance was taken into account in the measurement of the laser points on the screen. Convergence of the EFK and the model parameters is shown in Fig. 8 and Fig. 9, respectively.

Figure 8.

Convergence of EKF (rigid model).

Figure 9.

Divergence of model parameters (rigid model).

Although the EKF converged, the estimated model parameters could not converge as expected because the algorithm could not handle the joint flexibility. In other words, the EKF reduces the number of measurement errors by manipulating the D-H parameters to compensate the compliance effect. Consequently, the D-H parameters grow away from the true values.

5.2 Simulation with flexible joint model

The proposed algorithm was simulated with the same conditions. As shown in Fig. 10 and Fig. 11, the EKF and the model parameters converged.

Figure 10.

Convergence of EKF (flexible model).

Figure 11.

Convergence of model parameters (flexible model).

In order to investigate whether or not the estimated model parameters are effective on arbitrary postures as well as the test postures, 20 new test postures were created and examined, as shown in Fig. 9. Test postures were generated as follows:

${\text{θ}}_{\mathrm{k}}=\frac{10}{180}\pi {\left[\begin{array}{cccccc}-\sin q_{\mathrm{k}}& -\cos q_{\mathrm{k}}& -\sin q_{\mathrm{k}}& \cos q_{\mathrm{k}}& \sin q_{\mathrm{k}}& \cos q_{\mathrm{k}}\end{array}\right]}^T$

where θ _k (k = 1,···,20) represents the joint angles for the test postures and q_k is the k^th angle from 0 and 2∂.

Figure 12.

Laser points on the screen with test postures.

The ‘○’ markers and the ‘×’ markers represent the laser points on the screen with the true model parameters and the initial guess, respectively. The‘+’ markers indicates the laser points after calibration using the flexible joint model (proposed algorithm), while the ‘Δ’ markers correspond with the rigid joint model (previous algorithm). As shown in the figure, every calibrated point using the flexible joint model nearly coincides with the true points, whereas the calibrated points from the rigid joint model are located further apart from the true points.

In Table 3, the distance errors between the true points and the two calibrated points were compared with the errors between the true points and the non-calibrated points with the initially guessed parameters. As shown in the table, the accuracy was improved in every category.

Table 3.

Distance errors on the screen for the test postures.

		Mean	STD	Max	Min
Initial guess		0.22844	0.04023	0.25288	0.20436
After calibration	Rigid	0.09394	0.01510	0.14518	0.03806
	Flexible	0.00025	0.00009	0.00041	0.00010

6. Conclusion

In this paper, differential kinematics of a flexible manipulator were geometrically derived to calibrate model parameters such as the D-H parameters, the joint compliances and the mass positions. From a practical point of view, the proposed geometrical manner was compared with a numerical derivation technique. The proposed manner is not only free from the division by zero problem, but also much faster. In addition, by adopting the differential kinematics into a calibration system composed of a laser and a camera, the model parameters were iteratively updated. Consequently, the difference between the measured and the computed laser points on the screen was effectively decreased. A computer simulation with a 6 DOF humanoid-legged robot was carried out and confirmed that the proposed algorithm could calibrate the model parameters automatically. Applying the algorithm to actual robotic systems is an area for further work.