Sage Journals: Discover world-class research

Abstract

In this article, a parameter identification algorithm of kinematic calibration in parallel manipulators is proposed and compared with two conventional algorithms. By analyzing the identification matrix with singular value decomposition, the mathematical properties of the identification algorithms are investigated, and a novel algorithm based on singular value decomposition is proposed. Two parallel manipulators, which have been successfully applied in practice, are introduced and their system models are established, respectively. In order to show the reliability and feature of algorithms, the calibration simulations of these two parallel manipulators are performed with different parameter identification algorithms. Since these parallel manipulators are needed to guarantee 5-degree-of-freedom motion accuracy at least, the parameter identification simulations are performed based on both 6-degree-of-freedom measurement and 5-degree-of-freedom measurement. Compared with two conventional algorithms, the novel algorithm has the best solution in all simulation situations. The investigation in this article is very helpful for the parameter identification of kinematic calibration in parallel manipulators.

Keywords

Parallel manipulators kinematic calibration parameter identification accuracy

Introduction

Compared with a serial mechanism, a parallel manipulator has the good performances of high speed, heavy load handling, and low moving inertia due to its structural feature that its end effector connects to its base by several parallel chains.^1,2 Due to the merits, parallel manipulators have developed a mass of applications in the industry successfully, such as motion simulators, solar trackers, and tool heads of hybrid machine tools.^3–5 However, the large number of passive links and joints assembled in parallel manipulators leads to a mass of geometrical errors, which decrease the accuracy of parallel manipulators. There is no doubt that the methods to improve accuracy are extremely needed for parallel manipulators. The main methods can be divided into two categories: accuracy design and calibration.

Different from the method of accuracy design, the method of calibration is easy to perform in practice, and the cost is usually low. The robot calibration is usually divided into two categories: kinematic calibration for compensating geometry errors and dynamic calibration for compensating non-geometry errors. This article focuses on geometry errors, since the geometry errors are the main error source of parallel manipulators. The target of kinematic calibration is to make the theoretical output of kinematic model approach the real output of a parallel manipulator by correcting geometrical error parameters in kinematic model.^6,7 Usually, the process of kinematic calibrations includes four steps:^8,9 error modeling, measurement, parameter identification, and compensation. Among them, parameter identification is one of the most important steps in kinematic calibration, since it affects the whole process of kinematic calibration and determines the final calibration results. It is known that the geometrical error parameters which need to identify determine which error parameters should be given in the error model. The parameter identification model determines the measuring degree-of-freedom (DOF) number and influences the selection of the measurement method. The way to compensate the end-effector accuracy is determined by the geometrical error parameters which are identified. Furthermore, if the parameter identification algorithm adopted is not suitable, the identified geometrical error parameters cannot guarantee the end-effector accuracy due to the nonlinear ill-conditioning problem of the identification matrix. If wrong algorithms are adopted, even the geometrical error parameters cannot be obtained, and kinematic calibration cannot be performed. Although the error measurement methods are more important than the identification algorithms, this article is investigated in computer simulations. It is not needed to develop the measurement method which is used in practice. Hence, this article just focuses on the issue about the parameter identification algorithm of kinematic calibration in parallel manipulators.

The main task of parameter identification is to analyze the nonlinear relation between the output error parameters and geometrical error parameters and finally find a suitable numerical algorithm to obtain geometrical error parameters. The intelligent optimization method and the Newton iterative method are used most in practice. In intelligent optimization methods, radial basis function (RBF) neural network,¹⁰ fuzzy interpolation method,¹¹ and genetic algorithm¹² had been used widely. In Newton iterative methods, a mass of scholars^13–15 used this method based on the least square (LS) algorithm in the kinematic calibrations of different parallel manipulators successfully. A part of scholars,^16–18 respectively, proposed some mortified Newton iterative methods, key algorithms of which are mainly derived from LS such as the most used ridge estimation (RE) algorithm. Compared with intelligent optimization methods, Newton iterative methods, owning the explicit physical meaning, less calculating time, and the more reliable solutions, are simple and used more widely in practice.

However, some important issues about Newton iterative method are needed to be deeply investigated for its further applications. The first one is about the number of measuring DOF. Not all parallel manipulators are needed to guarantee 6-DOF motion accuracy, sometimes just needed to guarantee 5-DOF motion accuracy. Different motion accuracy requirements correspond to the different measuring DOF selections. Compared with 6-DOF measurement, 5-DOF measurement totally changes the mathematical property of an identification matrix. At this situation, the Newton iterative methods based on LS and RE are needed to consider which one is more suitable. The second one is about the determination of the algorithm parameter. The algorithm parameter of RE cannot be directly obtained by an explicit theoretical formula, just relying on trying over and over again, which results in a very tedious and time-consuming routine. The reason is that the mathematical property of an identification matrix is rarely considered when RE is used, which leads to lack of understanding of the identification process and does not improve nonlinear ill-conditioning problem of an identification matrix. According to the above-mentioned problems, in this article, the singular value decomposition (SVD) is employed to deeply investigate the mathematical properties of the identification matrix in LS and RE, a process in which an algorithm based on SVD is proposed and compared with LS and RE.

This article is organized as follows. By investigating the mathematical properties of the identification matrices by the SVD method, LS and RE algorithms are analyzed, and a novel algorithm based on SVD is proposed in section “Analyses of parameter identification algorithms.” In order to show the reliability and feature of algorithms, section “Parallel manipulators’ description” gives two parallel manipulators and establishes their kinematic models and error models according to the closed vector loop kinematic equations and the linear perturbation method, respectively. Section “Numerical example” deals with the parameter identification simulations of using these three algorithms, which are performed on the parallel manipulators of section “Parallel manipulators description,” with 6-DOF and 5-DOF measurements, respectively. Section “Discussions” discusses the simulation results. Conclusions are organized in section “Conclusion.”

Analyses of parameter identification algorithms

In this section, in order to analyze the mathematical properties of two algorithms (LS and RE), the parameter identification model is introduced, and the method of SVD is employed. Based on the analysis result, truncated singular value decomposition (TSVD) algorithm is proposed as an available selection for parameter identification.

Parameter identification model

The parameter identification model is carried out based on the error model. The first step is to set up an equation between the output error parameters and geometrical error parameters, namely, the error model. The error model can be obtained in the following expression

$Δ C = J Δ p$ (1)

where Δ C represents the output error vector, Δ p represents the geometrical error vector, and J represents the error transfer matrix.

Although equation (1) establishes an equation for output error parameters and geometrical error parameters, more equations are needed since the number of equations is too less to solve all unknown geometrical error parameters. Measuring output poses can add equations. Hence, a large number of measurement pose data are used, and the number sets as n. With the measurement pose data, the identification matrix W stacking by all J _i (i = 1–n) in row can be obtained as follows

$W = [\begin{matrix} J_{1} \\ ⋮ \\ J_{i} \\ ⋮ \\ J_{n} \end{matrix}]$ (2)

In order to identify the geometrical error parameters Δ p *, the residuals of the output error Δ C * should be calculated as follows

$Δ C * = [\begin{matrix} f_{1} (q_{1}, p + Δ p *) - C_{M 1} \\ ⋮ \\ f_{i} (q_{i}, p + Δ p *) - C_{Mi} \\ ⋮ \\ f_{n} (q_{n}, p + Δ p *) - C_{Mn} \end{matrix}]$ (3)

where f _i represents the forward kinematic function, q _i represents the input data corresponding to each measurement pose C _Mi, p is the normal kinematic parameter, and Δ p * is a zero vector in the first iterative calculation.

At last, based on equations (2) and (3), the parameter identification model can be shown as follows

$Δ C * = W Δ p *$ (4)

It should be pointed out that Δ p in error model is different from Δ p * in parameter identification model. Δ p represents the geometrical error parameters that should be identified. Δ p * represents the geometrical error parameters that can be identified. Usually, the parameter number of Δ p is more than the one of Δ p *.

Equation (4) is a linear equation in which the geometrical error parameters are solved by different algorithms. If the solved Δ p * does not satisfy the equation |Δ C *| < ε (ε is a judging criteria), then Newton iteration will go on with the following equation

$p_{m} = p_{m - 1} - Δ p *$ (5)

where m represents the iterative number, p _m represents the geometrical error vector in the m-iterative step, and p _m−1 represents the geometrical error vector in the m−1-iterative step. At last, the parameter identification can be ended when the equation |Δ C *| < ε is achieved.

It is seen that each iterative step makes the Δ p * approach the real solution, and the key of Newton iterative method is the algorithm which solves Δ p *. Hence, the different algorithms to solve Δ p * are needed to investigate deeply.

Analyses of two algorithms (LS and RE) with SVD

LS and RE are the two frequently used algorithms to solve Δ p * in Newton iterative method. They are effective for different situations.

First, LS is introduced and analyzed. The LS expression to solve geometrical error parameters is written as follows

$Δ p * = {(W^{T} W)}^{- 1} W^{T} Δ C *$ (6)

In order to analyze mathematical property of this algorithm, SVD is employed. W is written in SVD format

$W = US V^{T}$ (7)

It is a normal SVD format of W . U and V are, respectively, composed of a set of standard orthogonal basis, and S is a diagonal matrix. The singular values of W locate at the diagonal line of S according to the order from big to small. Equation (6) is rewritten as follows

$Δ p * = {({(U S V^{T})}^{T} U S V^{T})}^{- 1} {(U S V^{T})}^{T} Δ C * = V S^{- 1} U^{T} Δ C *$ (8)

Furthermore

$V^{T} Δ p * = S^{- 1} U^{T} Δ C *$ (9)

Equation (9) reveals that the vector Δ p * after performing V ’s orthogonal transformation reveals a linear relation with the vector Δ C * after performing U ’s orthogonal transformation. The linear relation is determined by S ⁻¹.

In order to figure out the linear relation, equation (9) is rewritten as follows

$Δ p * = \sum_{i = 1}^{j} \frac{u_{i}^{T} Δ C *}{τ_{i}} v_{i}$ (10)

where u _i is a set of standard orthogonal basis of U , v _i is a set of standard orthogonal basis of V , τ_i is the singular value of W , and j represents the dimension number of W .

In equation (10), it is found that solving Δ p * is similar to processing it in a solving system, solving process of which is shown in Figure 1. Δ p * is the solving system output. Δ C * is the solving system input, and the solving system is determined by j and τ_i. In general, the more complete data Δ C * are obtained, the more accurate the Δ p * can have in the given solving system, since LS has high robustness. However, once the data Δ C * are not complete, namely, not all motions of the end effector are measured, LS is not suitable, since the solving system with the given j and τ_i does not match with the incomplete system input, which leads to the nonlinear ill-conditioning problem.

Figure 1.

The solving process.

Second, RE is introduced and analyzed. Basically, RE is derived from LS. The RE expression is written as follows

$Δ p * = {(W^{T} W + λ_{n} I)}^{- 1} W^{T} Δ C *$ (11)

where λ_n is a vector and called as the RE algorithm parameter. Each element of λ_n should be bigger than 0. I is a unit matrix with the same size of W ^T W .

Equation (11) is rewritten in SVD format

$Δ p * = \sum_{i = 1}^{j} \frac{u_{i}^{T} Δ C *}{τ_{i} + λ_{i}} v_{i}$ (12)

where λ_i is the element of λ_n.

In the past investigations, λ_i is usually replaced by a constant number λ for convenience. This article also adopts this way as follows

$Δ p * = \sum_{i = 1}^{j} \frac{u_{i}^{T} Δ C *}{τ_{i} + λ} v_{i}$ (13)

From the above analysis, λ is introduced to modify τ_i of the solving system. With adjusting λ, namely, changing the solving system, RE can match with the different measuring DOFs. λ needs to be determined before the parameter identification is performed. However, no available formulas can give λ. It is hard to obtain a suitable λ, which just relies on trying different values and always consumes a lot of time.

Analysis of TSVD

LS uses the original solving system with the given j and τ_i, and RE changes τ_i to modify the solving system. Obviously, another algorithm to modify j of the solving system can be obtained, which is the TSVD. Before presenting TSVD, a deep investigation on how RE works is given.

From the point of mathematics, τ_i influences the variance of output solutions in the solving system. If the τ_i is not changed, the solution has no variance, and the real solution can be obtained. However, when not all motions of the end effector are measured, the real solution will be nonexistent. At this situation, an approximate solution which should approach the real solution is needed. Hence, the variance of solutions should be allowed. The direct way to allow the variance is to modify the singular value τ_i, as shown in equation (13). The suitable λ is used to overcome the nonlinear ill-conditioning problem, which introduces the variance of solutions and makes the modified solving system approach the original solving system. Since parts of motions of the end effector cannot be measured, the real solution cannot be obtained in fact. The solution of RE approaches the real solution as close as possible by adjusting λ. If λ is too small, RE may be performed like LS that cannot overcome the nonlinear ill-conditioning problem. If λ is too big, RE may destroy the linear relation between the output error parameters and geometrical error parameters (namely, the variance of solutions is allowed too big), which makes the solution of RE be far away from the real solution. Hence, in RE, the key point is to determine λ.

By adjusting τ_i, an approximate solution solved by RE approaches the real solution. Obviously, another algorithm to modify j also can achieve an approximate solution. From the point of mathematics, the dimension number j of W influences the resolution of solutions. The larger the dimension number, the higher the resolution of solutions, which leads to the higher accurate solution. It should be pointed out that the dimension number j cannot be too large since it may make the identification matrix W singular. In general, the largest dimension number is determined in the error modeling process, which makes the identification matrix W nonsingular. Hence, the only way to modify the dimension number j is to reduce it, namely, the method of TSVD. The expression is written as follows

$Δ p * = {(W_{k}^{T} W_{k})}^{- 1} W_{k}^{T} Δ C *$ (14)

where W _k is a part of W , and k is called as the algorithm parameter of TSVD. Equation (14) is rewritten in SVD format

$Δ p * = \sum_{i = 1}^{k} \frac{u_{i}^{T} Δ C *}{τ_{i}} v_{i}$ (15)

Similarly, the solution of TSVD approaches the real solution as close as possible by adjusting k. If k is too large, TSVD may be performed like LS that cannot overcome the nonlinear ill-conditioning problem. If k is too small, TSVD may make the resolution of solutions too low, which makes the solution of TSVD to not approach the real solution in a small range. The advantage of TSVD is the way to determine its algorithm parameter k, process of which is simple. In the above analysis, k should approach, but should not be equal to the dimension number j of W . Furthermore, k is an integer. It can be seen that the changing range of k is limited and regular. On the contrary, λ is an arbitrary number, and the only limit for λ is that it should be bigger than 0.

From the above analysis, the selection of the algorithm parameters (λ and k) and the number of the error parameters affect each other. But in the standard process of kinematic calibration, the number of the error parameters is determined first. After that the parameters (λ and k) should be selected to match with the error parameters. Hence, it is thought that the number of the error parameters affects the selection of the parameters (λ and k).

Parallel manipulators’ description

Two parallel manipulators are described in this section. The one is a 2-DOF parallel manipulator that is applied in a solar tracker. Another one is a 3-DOF parallel manipulator that is used as a parallel module in a 5-axis hybrid machine tool.

Construction of the 2-DOF parallel manipulator and its system modeling

The 2-DOF parallel manipulator is shown in Figure 2. It is composed of a moving platform, a base platform, and three chains. The base platform is fixed, and the moving platform has two rotational motions. The first chain contains a prismatic (P) joint, a universal (U) joint, and a spherical (S) joint. The second chain contains a P joint, a revolute (R) joint, and an S joint. The third chain contains an R joint and a U joint. Among them, these two P joints are actuated. According to the naming rule for a general parallel manipulator, this 2-DOF parallel manipulator is named as RU-PUS-PRS parallel manipulator.

Figure 2.

Construction of RU-PUS-PRS parallel manipulator.

A fixed Cartesian reference coordinate system O-XYZ is fixed at the base platform OB₁B₂. OB₁ is perpendicular to OB₂. The Z axis is perpendicular to the base platform. The X axis is pointed from B₁ to O, and Y axis satisfies the right-hand rule. A moving Cartesian reference coordinate system N-UVW is fixed at the moving platform NA₁A₂. NA₁ is perpendicular to NA₂. The W axis is perpendicular to the moving platform. The U axis is pointed from A₁ to N, and V axis satisfies the right-hand rule. C₁ and C₂ locate at U joint of the first chain and R joint of the second chain, respectively. O and N locate at the centers of R joint and U joint of the third chain, respectively. The guide rails, B₁C₁ and B₂C₂, are perpendicular to the base platform, respectively. Each P joint moves along the guide rails. Four Cartesian reference coordinate systems, axial directions of which are the same with O-XYZ, locate at B₁, B₂, C₁, and C₂, respectively. They are named as B₁-X_B₁Y_B₁Z_B₁, B₂-X_B₂Y_B₂Z_B₂, C₁-X_C₁Y_C₁Z_C₁, and C₂-X_C₂Y_C₂Z_C₂.

In order to simplify the expression of the rotational matrix, T–T rotational angle in Figure 3 is used to describe the orientation of the moving platform. The rotational matrix is expressed as follows

$R_{T T} (φ, θ, ψ) = [\begin{matrix} \cos φ & - \sin φ & 0 \\ \sin φ & \cos φ & 0 \\ 0 & 0 & 1 \end{matrix}] [\begin{matrix} \cos θ & 0 & \sin θ \\ 0 & 1 & 0 \\ - \sin θ & 0 & \cos θ \end{matrix}] [\begin{matrix} \cos φ & \sin φ & 0 \\ - \sin φ & \cos φ & 0 \\ 0 & 0 & 1 \end{matrix}] [\begin{matrix} \cos ψ & - \sin ψ & 0 \\ \sin ψ & \cos ψ & 0 \\ 0 & 0 & 1 \end{matrix}]$

Figure 3.

T–T rotational angle.

The angle ψ is around the normal direction of the moving platform. When the moving platform of a solar tracker follows the sun in the sky, the changing of angle ψ will not affect the tracking accuracy. Hence, angle ψ can be ignored and only two angles are kept in the rotational matrix.

Based on the vector loops OB₁C₁A₁N and OB₂C₂A₂N, the vector loop equations can be written as follows

$H + R_{TT} a_{1} = b_{1} + R_{B 1} q_{1} + R_{B 1} R_{C 1} l_{1}$ (16)

$H + R_{TT} a_{2} = b_{2} + R_{B 2} q_{2} + R_{B 2} R_{C 2} l_{2}$ (17)

where a ₁, b ₁, l ₁, a ₂, b ₂, and l ₂ represent the position vectors of NA₁, OB₁, C₁A₁, NA₂, OB₂, and C₂A₂, respectively. H represents the position vector of the moving platform. q ₁ and q ₂ represent the input vector of each P joint. R _B1, R _B2, R _C1, and R _C2 represent the rotational matrices of the coordinate systems B₁-X_B₁Y_B₁Z_B₁, B₂-X_B₂Y_B₂Z_B₂, C₁-X_C₁Y_C₁Z_C₁, and C₂-X_C₂Y_C₂Z_C₂.

Considering the kinematic constraints of joints and links, the kinematic model can be obtained as follows

${(R_{B 1} R_{C 1} ω_{1})}^{T} [H + R_{TT} a_{1} - (b_{1} + R_{B 1} q_{1})] = ‖ R_{B 1} R_{C 1} l_{1} ‖$ (18)

${(R_{B 2} R_{C 2} ω_{2})}^{T} [H + R_{TT} a_{2} - (b_{2} + R_{B 2} q_{2})] = ‖ R_{B 2} R_{C 2} l_{2} ‖$ (19)

${(R_{B 1} R_{C 1} e_{y})}^{T} [H + R_{TT} a_{1} - (b_{1} + R_{B 1} q_{1})] = 0$ (20)

${(R_{B 2} R_{C 2} e_{x})}^{T} [H + R_{TT} a_{2} - (b_{2} + R_{B 2} q_{2})] = 0$ (21)

${(e_{y})}^{T} H = 0$ (22)

$‖ H ‖ = r$ (23)

where r represents the length of the link ON, e _x and e _y are the unit vectors of X and Y axes, and ω ₁ and ω ₂ are the unit direction vectors of l ₁ and l ₂, respectively.

Based on equations (18)–(23), the error model can be obtained with the linear perturbation method

$\begin{matrix} {(ω_{1})}^{T} [Δ H + Δ α \times R_{TT} a_{1} + R_{TT} Δ a_{1} \\ - (Δ b_{1} + Δ θ_{B 1} \times q_{1} + e_{z} Δ q_{1}) \end{matrix}] = Δ l_{1}$ (24)

$\begin{matrix} {(ω_{2})}^{T} [Δ H + Δ α \times R_{TT} a_{2} + R_{TT} Δ a_{2} \\ - (Δ b_{2} + Δ θ_{B 2} \times q_{2} + e_{z} Δ q_{2}) \end{matrix}] = Δ l_{2}$ (25)

${(e_{y})}^{T} [Δ H + Δ α \times R_{T T} a_{1} + R_{T T} Δ a_{1} - (Δ b_{1} + Δ θ_{B 1} \times q_{1} + e_{z} Δ q_{1} + (Δ θ_{B 1} + Δ θ_{C 1}) \times l_{1})] = 0$ (26)

${(e_{x})}^{T} [Δ H + Δ α \times R_{T T} a_{2} + R_{T T} Δ a_{2} - (Δ b_{2} + Δ θ_{B 2} \times q_{2} + e_{z} Δ q_{2} + (Δ θ_{B 2} + Δ θ_{C 2}) \times l_{2})] = 0$ (27)

${(Δ θ_{O} \times e_{y})}^{T} Δ H = 0$ (28)

$‖ Δ H ‖ = Δ r$ (29)

where e _z = [0 0 1]^T. Δ H represents the position error vector of the output. Δ α represents the orientation error vector of the output. Δ a ₁, Δ b ₁, Δ a ₂, and Δ b ₂ represent the geometrical error vectors of a ₁, b ₁, a ₂, and b ₂. Δq₁, Δl₁, Δq₂, and Δl₂ represent the length error parameters of q ₁, l ₁, q ₂, and l ₂. Δ θ _B1, Δ θ _C1, Δ θ _B2, and Δ θ _C2 represent the angle error vectors of R _B1, R _C1, R _B2, and R _C2. Δ θ _O represents the angle error vector of the coordinate system O-XYZ. Δr represents the length error parameter of the link ON. It should be pointed out that R _B1, R _B2, R _C1, and R _C2 are omitted in equations (24)–(29) since they are the unit matrices.

In order to make the identification matrix nonsingular, redundant error parameters should be eliminated. By analyzing equations (24)–(29) in scalar expressions, the redundant error parameters which own the same coefficients with others in the error model are eliminated. It is noted that the same analysis is performed in the following parallel manipulator and not repeated again. After the analysis, in total, 25 error parameters can be identified. They are Δa_1x, Δa_1y, Δa_1z, Δb_1x, Δb_1y, Δ θ _B1x, Δ θ _B1y, Δq₁, Δ θ _C1x, Δ θ _C1z, Δl₁, Δa_2x, Δa_2y, Δa_2z, Δb_2x, Δb_2y, Δ θ _B2x, Δ θ _B2y, Δq₂, Δ θ _C2y, Δ θ _C2z, Δl₂, Δ θ _oy, Δ θ _oz, and Δr.

Construction of the 3-DOF parallel manipulator and its system modeling

The 3-DOF parallel manipulator is shown in Figure 4. It is composed of a moving platform, a base platform, and three supporting chains with the identical kinematic structure. The base platform is fixed, and the moving platform has two rotational motions and one translational motion along Z axis. Each PRS chain contains a P joint, an R joint, and an S joint. Among them, the P joint is actuated. According to the naming rule for a general parallel manipulator, this 3-DOF parallel manipulator is named as 3-PRS parallel manipulator.

Figure 4.

Construction of 3-PRS parallel manipulator.

A fixed Cartesian reference coordinate system O-XYZ with Z axis being vertically placed is fixed at the center of the base platform B₁B₂B₃. The X axis is pointed from O to B₁, and Y axis satisfies the right-hand rule. A moving Cartesian reference coordinate system N-UVW is fixed at the center of the moving platform A₁A₂A₃. W axis is perpendicular to the moving platform. The U axis is pointed from N to A₁, and V axis satisfies the right-hand rule. B_iC_i on the base platform for i = 1, 2, and 3 represents the guide rail of each chain. The guide rail of each chain is perpendicular to the base platform. Each P joint moves along B_iC_i. P joint connects the base platform with R joint. R joint connects the P joint with the C_iA_i link. S joint connects the moving platform with C_iA_i link. All the joints attached to the base or moving platform are symmetrically distributed at vertices of the equilateral triangles. The Cartesian reference coordinate system B_i-X_BiY_BiZ_Bi locates at the P joint of each chain. Z_Bi axis is perpendicular to the base platform. X_Bi axis is pointed along the direction from O to B_i, and Y_Bi axis satisfies the right-hand rule. The Cartesian reference coordinate system C_i-X_CiY_CiZ_Ci locates at the R joint of each chain. The directions of Z_Ci axis and X_Ci axis are the same with the ones of Z_Bi axis and X_Bi axis. Y_Ci axis satisfies the right-hand rule. B_i-X_BiY_BiZ_Bi and C_i-X_CiY_CiZ_Ci are used to describe the orientations of P joints and R joints, respectively.

In order to stay the same with RU-PUS-PRS parallel manipulator, T–T rotational angle is also introduced to describe the orientation of the moving platform. Furthermore, the accuracy of the angle ψ is not necessary to assure when the five-axis machine tool machines a part. Hence, angle ψ also can be ignored.

Based on the vector loop OB_iC_iA_iN, the vector loop equations are written as

$H + R_{TT} R_{i} a_{i} = R_{i} b_{i} + R_{i} R_{Bi} q_{i} + R_{i} R_{Bi} R_{Ci} l_{i}$ (30)

where R _i represents the rotational matrix of each chain. Furthermore, a _i, b _i, and l _i represent the position vectors of NA_i, OB_i, and C_iA_i. H represents the position vector of the moving platform. q _i represents the input vector of each P joint. R _Bi and R _Ci represent the rotational matrices of the P joints and the R joints, respectively.

Considering the kinematic constraints of joints and links, the kinematic model can be obtained as

${(R_{i} ω_{i})}^{T} (H + R_{TT} R_{i} a_{i} - R_{i} b_{i} - R_{i} R_{Bi} q_{i}) = ‖ R_{i} R_{Bi} R_{Ci} l_{i} ‖$ (31)

${(R_{i} e_{y})}^{T} (H + R_{TT} R_{i} a_{i} - R_{i} b_{i} - R_{i} R_{Bi} q_{i}) = 0$ (32)

where ω _i represents the unit direction vector of l _i. e _y = [0 1 0]^T.

Based on equations (31) and (32), the error model can be obtained with the linear perturbation method

$\begin{matrix} {(R_{i} ω_{i})}^{T} (Δ H + Δ α \times R_{TT} R_{i} a_{i}) \\ = {(R_{i} ω_{i})}^{T} (R_{i} Δ b_{i} + R_{i} Δ θ_{Bi} \times e_{z} q_{i} + R_{i} e_{z} Δ q_{i} \\ + R_{i} ((Δ θ_{Bi} + Δ θ_{Ci}) \times) ω_{i} l_{i} + R_{i} ω_{i} Δ l_{i} - R_{TT} R_{i} Δ a_{i}) \end{matrix}$ (33)

$\begin{matrix} {(R_{i} e_{y})}^{T} (Δ H + Δ α \times R_{T T} R_{i} a_{i}) \\ = {(R_{i} e_{y})}^{T} (R_{i} Δ b_{i} + R_{i} Δ θ_{B i} \times e_{z} q_{i} + R_{i} e_{z} Δ q_{i} + R_{i} ((Δ θ_{B i} + Δ θ_{C i}) \times) ω_{i} l_{i} + R_{i} ω_{i} Δ l_{i} - R_{TT} R_{i} Δ a_{i}) \end{matrix}$ (34)

where Δ H represents the position error vector of the output. E is a unit matrix. Δ α represents the orientation error vector of the output. Δ a _i and Δ b _i, respectively, represent the geometrical error vectors of a _i and b _i. Δ θ _Bi and Δ θ _Ci, respectively, represent the angle error vectors of R _Bi and R _Ci. q_i and l_i, respectively, represent the original length parameters of q _i and l _i. Δq_i and Δl_i, respectively, represent the length error parameters of q _i and l _i. R _Bi and R _Ci are omitted in equations (33) and (34) since they are the unit matrices. e _z = [0 0 1]^T.

The redundant error parameters can be eliminated by analyzing equations (33) and (34) in scalar expressions until no error parameters are coupling with others. A total of 11 error parameters on each chain are kept in this error model. They are Δa_ix, Δa_iy, Δa_iz, Δb_ix, Δb_iy, Δq_i, Δl_i, Δθ_Bix, Δθ_Biy, Δθ_Cix, and Δθ_Ciz. In total, there are 33 error parameters for this parallel manipulator.

Numerical example

In order to show the validity of each algorithm, the parameter identification simulations based on LS, RE, and TSVD are performed, respectively, on the RU-PUS-PRS parallel manipulator and 3-PRS parallel manipulator with 6-DOF measurement. At the same time, the simulations with 5-DOF measurement are also performed, since 6-DOF measurement is hard to achieve in practice, and 5-DOF measurement is enough to guarantee the motion accuracy of these parallel manipulators.

Simulation process of parameter identification

The parameter identification simulations are performed in MATLAB. The detailed process is organized as follows and shown in Figure 5:

Step 1. Initialize the load initial kinematic parameters. Δ p _real is given for simulations.

Step 2. Establish the forward kinematic model f ( q , p ) = C .

Step 3. Select measurement poses. Since the selection of measurement poses can affect the identification result, the distribution of measurement poses keeps in the same in each simulation. The data of measurement poses are given in Appendix 1.

Step 4. Determine the real measurement data. Calculate the real measurement data f ( q , p +Δ p _real) = C _real in the given measurement poses. In order to simulate more reliably, C _real + σ is used as the measurement data. σ represents the measurement error. The measurement errors for two parallel manipulators are given in Appendix 1.

Step 5. Obtain the error transfer matrix J . Error model is established based on the corresponding kinematic model, and the error transfer matrix J can be obtained from the error model.

Step 6. Generate the identification matrix W . The parameter identification matrix W can be generated by stacking J in row with the given measurement poses.

Step 7. Identify geometrical error parameters Δ p *. In this step, the different algorithms are applied to solve geometrical error parameters Δ p *.

Step 8. Judge and iterate. If the solved Δ p * does not satisfy the equation |Δ C *| < ε (ε is a judging criteria), then the iteration will go on with p _m + 1 = p _m−Δ p * and back to Step 5.

Step 9. Save the result of geometrical error parameters Δ p _cal.

Figure 5.

Simulation process of parameter identification.

Obtaining the geometrical error parameters is not the end of kinematic calibration. The compensation results still should be simulated for testing these parameter identification results of each algorithm. In general, the motion accuracy of the parallel manipulator should be evaluated by the output error of the moving platform before and after calibration in a given workspace. The original output error before calibration is obtained by f ( q , p + Δ p _real) − f ( q , p ) = Δ C *_orig. The output error after calibration is obtained by f ( q , p + Δ p _real) − f ( q , p + Δ p _cal) = Δ C _cal*. Both Δ C *_orig and Δ C _cal* represent the position and orientation errors of the moving platform, respectively. Since the orientation error evaluation in practice is influenced by the length of the normal direction vector of the moving platform and in order to simplify the evaluation index of output error, just the length of the maximal position output error vector in the given output range is adopted as the output error value to evaluate these parameter identification algorithms.

Simulation results

In this section, the parameter identification simulations in different situations are performed based on the above process, and the data are presented in Tables 1 –4.

Table 1.

Normal structure parameters of RU-PRS-PUS parallel manipulator.

Parameters	Symbol	Unit	Value
Length of link, NA₁	\| a ₁\|	mm	220
Length of link, NA₂	\| a ₂\|	mm	200
Length of link, OB₁	\| b ₁\|	mm	290
Length of link, OB₂	\| b ₂\|	mm	290
Length of link, C₁A₁	\| l ₁\|	mm	560
Length of link, C₂A₂	\| l ₂\|	mm	520
Length of link, ON	\| r \|	mm	685

Table 2.

Normal structure parameters of 3-PRS parallel manipulator.

Parameters	Symbol	Unit	Value
Length of link, NA_i	\| a _i\|	mm	201
Length of link, OB_i	\| b _i\|	mm	220
Length of link, C_iA_i	\| l _i\|	mm	487

Table 3.

Output range for evaluating RU-PRS-PUS parallel manipulator.

Parameters	Symbol	Unit	Value range
Direction angle of R _TT rotational matrix	φ	rad	0–6.2831852
Tilting angle of R _TT rotational matrix	θ	rad	0–0.1745329

Table 4.

Output range for evaluating 3-PRS parallel manipulator.

Parameters	Symbol	Unit	Value range
Position along Z axis	Z	mm	700 (unchanged)
Direction angle of R _TT rotational matrix	φ	rad	0–6.2831852
Tilting angle of R _TT rotational matrix	θ	rad	0–0.1745329

It should be noted that the λ in RE and k in TSVD both greatly affect the identification results. From the above analysis, the k depends on the dimension number of the identification matrix and is determined from big to small. But no investigations give the determination of λ. Considering the mathematical properties, λ is determined according to the smallest singular value of the identification matrix. In the simulations, λ is basically adjusted from the smallest singular value to 1 with the same intervals.

From the analysis in section “Analyses of parameter identification algorithms,” it is known that both λ and k cannot be too big or small. Hence, the curve reflecting the relation curve between the output error and λ or k should be a “V-” or “U”-shaped curve. It can be thought that the best solution can be obtained when a “V-” or “U”-shaped trend emerges. In the following simulations, the best solution of each algorithm is determined by the bottom data of the “V-” or “U”-shaped curve.

The relations between the output error and algorithm parameters (λ and k) in RU-PRS-PUS parallel manipulator with 6-DOF measurement are plotted in Figure 6. The corresponding iteration processes of LS, RE, and TSVD are plotted in Figure 7. The relations between the output error and algorithm parameters (λ and k) in RU-PRS-PUS parallel manipulator with 5-DOF measurement are plotted in Figure 8. Since LS is not suitable with 5-DOF measurement, just the corresponding iteration processes of RE and TSVD are plotted in Figure 9. The corresponding simulation results of 3-PRS parallel manipulator with RE and TSVD are plotted in Figures 10 and 11, respectively. The iteration processes of each algorithm in 3-PRS parallel manipulator are similar with the ones of RU-PRS-PUS parallel manipulator and are not shown again. Since the identified error parameters on each chain are similar with others, the identified error parameters of one chain are organized in Tables 5 and 6. The output errors of each algorithm are organized in Tables 7 and 8. After calibration simulations, the output errors after calibration corresponding to each algorithm parameter are plotted in Figures 6, 8, 10 and 11. In Figures 7 and 9, the horizontal axis represents the iteration number, and the vertical axis represents the value of every error parameter.

Figure 6.

The relations between the output error and algorithm parameters (λ and k) in RU-PRS-PUS parallel manipulator with 6-DOF measurement: (a) RE and (b) TSVD.

Figure 7.

The iteration processes of LS, RE, and TSVD in RU-PRS-PUS parallel manipulator with 6-DOF measurement.

Figure 8.

The relations between the output error and algorithm parameters (λ and k) in RU-PRS-PUS parallel manipulator with 5-DOF measurement: (a) RE and (b) TSVD.

Figure 9.

The iteration processes of RE and TSVD in RU-PRS-PUS parallel manipulator with 5-DOF measurement.

Figure 10.

The relations between the output error and algorithm parameters (λ and k) in 3-PRS parallel manipulator with 6-DOF measurement: (a) RE and (b) TSVD.

Figure 11.

The relations between the output error and algorithm parameters (λ and k) in 3-PRS parallel manipulator with 5-DOF measurement: (a) RE and (b) TSVD.

Table 5.

The identified error parameters of one chain in RU-PRS-PUS parallel manipulator.

Error parameters	Real solution	6-DOFs			5-DOFs
Error parameters	Real solution	LS	RE	TSVD	RE	TSVD
Δb_1x (mm)	0.3	−1.722901	0.211841	0.113960	0.000139	−0.001614
Δb_1y (mm)	0.2	0.143312	0.474392	−0.083124	−0.005608	−0.042556
Δa_1x (mm)	0.4	0.310149	0.216371	0.291129	−0.000347	−0.024176
Δa_1y (mm)	−0.1	−0.102147	−0.096061	−0.146733	0.014761	0.124699
Δa_1z (mm)	0.2	0.283112	0.095323	0.221615	−0.005513	0.097093
Δθ_B_1x (rad)	0.00017	0.000035	0.000533	−0.000221	0.006439	0.001692
Δθ_B_1y (rad)	0.00034	0.001803	0.00261	0.002059	−0.000519	−0.000692
Δθ_C_1x (rad)	0.00022	0	0	0	0.002354	0.001258
Δθ_C_1z (rad)	0.00015	−0.000099	0.000077	0.000057	−0.002463	−0.000099
Δq₁ (mm)	0.3	−6.523845	−3.451253	0.163966	0.043	0.144353
Δl₁ (mm)	0.2	7.174132	3.600782	0.201847	0.043229	0.150479

DOFs: degrees of freedom; LS: least square; RE: ridge estimation; TSVD: truncated singular value decomposition.

Table 6.

The identified error parameters of one chain in 3-PRS parallel manipulator.

Error parameters	Real solution	6-DOFs			5-DOFs
Error parameters	Real solution	LS	RE	TSVD	RE	TSVD
Δb_1x (mm)	3	2.546116	2.320435	2.262294	2.317067	2.260464
Δb_1y (mm)	2	−0.230475	0.267064	2.143774	0.244967	2.143368
Δa_1x (mm)	4	4.055073	4.028317	4.015571	4.020576	4.007861
Δa_1y (mm)	−1	−1.020147	−1.017515	−1.001506	2.022574	2.034316
Δa_1z (mm)	2	1.930542	1.924592	1.924262	1.923355	1.923018
Δθ_B_1x (rad)	0.00017	−0.000246	−0.000118	−0.000093	−0.000124	−0.000098
Δθ_B_1y (rad)	0.00034	0.002392	0.002868	0.003126	0.002879	0.003132
Δθ_C_1x (rad)	0.00022	−0.003694	−0.002868	0.000894	−0.002903	0.000902
Δθ_C_1z (rad)	0.00015	−0.000153	−0.000111	0.000136	−0.000112	0.000136
Δq₁ (mm)	3	2.798887	2.847576	2.58509	2.846039	2.584480
Δl₁ (mm)	2	2.130641	2.070091	2.330329	2.070737	2.329989

DOFs: degrees of freedom; LS: least square; RE: ridge estimation; TSVD: truncated singular value decomposition.

Table 7.

The output errors of RU-PRS-PUS parallel manipulator in different situations.

Measurement DOF	6-DOFs			5-DOFs
Algorithm	LS	RE	TSVD	RE	TSVD
Output error after calibration (10⁻³ mm)	5.8	5.657	5.541	5.977	3.554

DOFs: degrees of freedom; LS: least square; RE: ridge estimation; TSVD: truncated singular value decomposition.

Table 8.

The output errors of 3-PRS parallel manipulator in different situations.

Measurement DOF	6-DOFs			5-DOFs
Algorithm	LS	RE	TSVD	RE	TSVD
Output error after calibration (10⁻³ mm)	6.9	6.232	6.018	6.227	6.005

DOFs: degrees of freedom; LS: least square; RE: ridge estimation; TSVD: truncated singular value decomposition.

It should be pointed out that the qualitative relationship between the output error and algorithm parameters (λ and k) is suitable for any other parallel manipulators. The quantitative relationship is specific for other parallel manipulators and its key is to obtain the suitable algorithm parameters. For λ, it is needed to find the biggest and smallest values until the solution does not diverge. Then, some suitable values of λ are interpolated between the biggest and smallest values. Compared with choosing λ, k is easy to determine. Since the dimension number of the identification matrix W is j, k should be searched in the order that it is decreased with integers from j−1. In simulation process, one algorithm parameter (λ and k) corresponds to one output error after calibration. After that the quantitative relationship can be plotted.

Discussions

According to the above-mentioned simulation results, several issues about the identified error parameters, LS, RE, and TSVD, are discussed.

In Tables 5 and 6, not all the acceptable identified error parameters approach the real error parameters given in the simulations, no matter which kind of measurement is used or which kind of parallel manipulator is performed. The reason is that with the measurement errors, the real error parameters have not been able to match with the output measurement data. Since the objective of parameter identification is to approach the output measurement data without restraining the changing range of error parameters, some identified error parameters should deviate from the real error parameters in order to fit for the output measurement data. And, different algorithms have different methods to restrain measurement errors, which lead to the different identified error parameters. Furthermore, the accuracy of some local error parameters is affected by the measurement errors. But the accuracy of the global error parameter vectors is effected by the measuring DOF. The 6-DOF measurement can obtain the higher accuracy than the 5-DOF measurement.

In Tables 7 and 8, LS is available with 6-DOF measurement and is totally failed with 5-DOF measurement, the solution of which is emanative. Compared with RE and TSVD, the solution of LS is not the best. But it should be pointed out that LS is still a good algorithm with 6-DOF measurement, since no algorithm parameters are needed to determine, and the solution of LS is usually stable.

In Figures 6(a), 8(a), 10(a), and 11(a), RE is available with both 6-DOF measurement and 5-DOF measurement, no matter which parallel manipulator is simulated. However, the suitable λ in different parallel manipulators is totally different, which ranges from 1 to 10⁻⁵. This leads to an obstacle to use RE, since λ is hard to achieve, as this algorithm parameter varies in an unlimited range and is not determined with explicit rules. Furthermore, in practice, this process to change λ will lead to a very tedious and time-consuming routine since every set of identified error parameters with a given λ should be performed on the real parallel manipulator to judge whether the parameter identification result is suitable. It should be noted that λ is chosen in several reasonable values in this article but may be not the best. If it goes on with changing λ in a careful searching, the better identification solution may be obtained in theory. However, it cannot be performed with such a careful searching in practice unless the absolute high-precision identification solution is needed without considering the time cost of searching λ.

In Figures 6(b), 8(b), 10(b), and 11(b), TSVD is available with both 6-DOF measurement and 5-DOF measurements, no matter which parallel manipulator is simulated. In Tables 7 and 8, it is seen that TSVD has the smallest output error values in many situations. Moreover, the algorithm parameter k of TSVD is easy to determine, since k is an integer and its value varies in a limit range with an explicit rule, which means that the error parameters can be identified with a great solution after trying several different k. Another advantage of TSVD is the algorithm convergence. In Figures 7 and 9, the convergence of TSVD is usually faster than the ones of RE and LS.

Furthermore, the method of restraining the measurement error is obtained. The key of this method is to find the best algorithm parameters (λ and k). For λ, it is needed to find the biggest and smallest values until the solution does not diverge. Then, some suitable values of λ are interpolated between the biggest and smallest values. These values can be used for λ in the simulations. Then, the best value of λ among them can be found when a “V-” or “U”-shaped trend like in Figures 6(a), 8(a), 10(a), and 11(a) emerges. Compared with choosing λ, k is easy to find. Since the dimension number of the identification matrix W is j, k should be searched in the order that it is decreased with integers from j−1. Then, the best value of k among them can be found when a “V-” or “U”-shaped trend like in Figures 6(b), 8(b), 10(b) and 11(b) emerges.

Conclusion

The parameter identification algorithms of kinematic calibration in parallel manipulators are investigated in this article. By analyzing the identification matrix with SVD, the mathematical properties of LS and RE are investigated, and TSVD is proposed. Parameter identification is to obtain Δ p * in a solving system. Δ p * is the solving system output. Δ C * is the solving system input, and the solving system is determined by j and τ_i. LS keeps the system constant with the given j and τ_i and just can be available in 6-DOF measurement. RE mortifies the system with changing τ_i, which is helpful for restraining the measurement error in 6-DOF measurement and obtaining the approximate solving system to approach the real solution in 5-DOF measurement. TSVD mortifies the system with changing j, which has a similar ability with RE. In the parameter identification simulations with LS, RE, and TSVD, TSVD has the smallest output error values, the best algorithm convergence, and the easier way to determine the algorithm parameter than RE. RE is available in many situations and its solution is better than LS. But the algorithm parameter of RE is hard to achieve. LS is just available in 6-DOF measurement.

It should be pointed out that each algorithm has the respective advantages and can be used in the different situations. LS without determining algorithm parameters can be performed in 6-DOF measurement within a short time. RE may obtain the best solution without considering the cost of searching λ. TSVD can obtain a suitable solution with changing k in a limited range.

This article gives the theoretical basis of LS, RE, and TSVD and shows the respective advantages of each algorithm, which are very helpful for the parameter identification in kinematic calibration of parallel manipulators.

Footnotes

Academic Editor: Jianqiao Ye

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research,authorship,and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research,authorship,and/or publication of this article: This work is supported by the National Natural Science Foundation of China (grant nos 51575307 and 51225503),the Science and Technology Major Project-Advanced NC Machine Tools & Basic Manufacturing Equipments (2013ZX04004021 and 2014ZX04002051),and Top-Notch Young Talents Program of China.

References

Piras

Cleghorn

Mills

JK.

Dynamic finite-element analysis of a planar high-speed, high-precision parallel manipulator with flexible links. Mech Mach Theory 2005; 40: 849–862.

Wang

. Dynamics and control of a planar 3-DOF parallel manipulator with actuation redundancy. Mech Mach Theory 2009; 44: 835–849.

Xie

Liu

You

. Type synthesis of 2T1R-type parallel kinematic mechanisms and the application in manufacturing. Robot CIM: Int Manuf 2014; 30: 1–10.

Chen

Wang

Design and dynamics of a novel solar tracker with parallel mechanism. IEEE/ASME T Mech 2016; 21: 88–97.

Zhang

Wang

Optimum design and performance comparison of a redundantly actuated solar tracker and its nonredundant counterpart. Sol Energy 2016; 127: 36–47.

Ramesh

Mannan

Poo

AN.

Error compensation in machine tools—a review: part I: geometric, cutting-force induced and fixture-dependent errors. Int J Mach Tool Manu 2000; 40: 1235–1256.

Roth

Mooring

Ravani

An overview of robot calibration. IEEE J Robot Autom 1987; 5: 377–385.

Duan

Qiu

Duan

. Calibration and motion control of a cable-driven parallel manipulator based triple-level spatial positioner. Adv Mech Eng 2014; 6: 368018.

Majarena

Santolaria

Samper

. An overview of kinematic and calibration models using internal/external sensors or constraints to improve the behavior of spatial parallel mechanisms. Sensors 2010; 10: 10256–10297.

10.

Ding

. RBF neural network based kinematic calibration of a planar parallel robot. In: Proceedings of the international conference on intelligent computing, Kunming, China, 16–19 August 2006. New York: IEEE.

11.

Bai

Zhuang

Wang

. Calibration of parallel machine tools using fuzzy interpolation method. In: Proceedings of the 2008 IEEE international conference on technologies for practical robot applications, Woburn, MA, 10–12 November 2008. New York: IEEE.

12.

Gao

Cheng

Zhao

Kinetic self-calibration of novel 5-DOF parallel machine tool. Comput Integr Manuf 2007; 13: 738–743.

13.

Fan

Zhao

. Calibration of a parallel mechanism in a serial-parallel polishing machine tool based on genetic algorithm. Int J Adv Manuf Tech 2015; 81: 27–37.

14.

Ibaraki

Yokawa

Kakino

. Kinematic calibration on a parallel kinematic machine tool of the Stewart platform by circular tests. In: Proceedings of the 2004 American control conference, Boston, MA, 30 June–2 July 2004. New York: IEEE.

15.

Weck

Staimer

Accuracy issues of parallel kinematic machine tools. Proc IMechE, Part K: J Mul Dyn 2002; 216: 51–57.

16.

Zhang

Online robot calibration based on hybrid sensors using Kalman Filters. Robot CIM: Int Manuf 2015; 31: 91–100.

17.

Besnard

Khalil

. Identifiable parameters for parallel robots kinematic calibration. In: Proceedings of the IEEE international conference on robotics and automation, Seoul, South Korea, 21–26 May, 2001. New York: IEEE.

18.

Fang

Dang

A step identification method of joint parameters of robots based on the measured pose of end-effector. Proc IMechE, Part C: J Mechanical Engineering Science. Epub ahead of print 3 February 2015. DOI: 10.1177/0954406215569589.