External Kinematic Calibration of Hybrid Kinematics Machine Utilizing Lower-DOF Planar Parallel Kinematics Mechanisms

This paper presents the implementation of nonlinear least squares and iterative linear least squares algorithms for external kinematic calibration of a hybrid kinematics machine composed of two 3PRR planar parallel kinematics mechanisms by utilizing a laser tracker. First the hand-eye and robot-world transformations were obtained by a separable closed-form solution and refined by the nonlinear least squares. Subsequently, the geometric parameters of the machine’s mechanisms were estimated using the two algorithms. Due to the rank deficiency, we implemented the nonlinear least squares algorithm through a subset selection approach in which we performed the estimation in two steps. We iterated the closed-form solution of the linear least squares until the solution converges to the actual values. We have shown that the nonlinear least squares algorithm successfully refined the hand-eye and robot-world transformations and outperformed the iterative linear squares algorithm in the estimation of the geometric parameters of the mechanisms.


Introduction
The parallel kinematics mechanisms (PKMs) which consist of a base, some legs, and a moving platform in a closedchain configuration have been proposed and developed to achieve more accuracy in the precision manipulation.Furthermore, hybrid kinematics machines (HKMs) have been introduced and developed to integrate the advantages of both the serial kinematics mechanism (SKM) and PKM.Any of the three schemes can compose an HKM: (1) serially connecting two or more PKMs, (2) serially connecting a PKM with an SKM, or (3) serially combining two or more PKMs, or a PKM with an SKM, through a rigid connection.In the third scheme, the rigid connection is considered a type of serial connection.In some applications such as machining, a PKM is typically used together with an SKM or another PKM and therefore becomes an HKM.Moreover, more attention has been given to the use of lower degree-of-freedom (DOF) PKMs due to less complexity in their modeling.
For a precision manipulation, an accurate pose of the moving platform is the aim to be achieved during its controlled motion.Since the control is commonly performed in the joint space, forward kinematics of the mechanism is required to transform the active joint positions to the pose of the moving platform.As the trajectory is commonly planned in the task space, inverse kinematics is needed to transform the planned trajectory to the joint space in which the control is applied.Consequently, the accuracy of the kinematics affects positioning accuracy.A kinematic calibration is commonly conducted to estimate the kinematic parameters accurately and subsequently, a compensation is performed by correcting the software, which is, in this case, the kinematics of the machine.
In general, the kinematic calibration can be either an external calibration, a constrained calibration, or a self-calibration [1].The external and constrained calibration can be conducted offline, and the number of calibration poses is usually limited.Here the number and the choice of the calibration poses become an issue.The number of the calibration poses should be adequate, whereas the selection of the calibration poses is shown to be vital as it affects the quality of the calibration.According to Bai and Theo [2], the optimal number of calibration poses is 10 or more.Moreover, a full pose measurement typically requires less number of calibration poses to get an accurate estimation [1].
In a nutshell, the kinematic calibration is performed by firstly determining the parameters to be calibrated, followed by modeling, measurement (data collection), identification (estimation), and compensation (correction) [3,4].Various algorithms in the identification stage are used to minimize the residual error.To mention some of them in the kinematic calibration of PKMs and HKMs, Vischer and Clavel [5] used the Levenberg-Marquardt algorithm for kinematic calibration of a Delta robot.Yang et al. [6] used a linear least squares algorithm for kinematic calibration of a spatial 3-DOF PKM.Patel and Ehmann [7] as well as Iurascu and Park [8] applied the total least squares algorithm for kinematic calibration of PKMs.Yu [9] used nonlinear least squares for kinematic calibration of a hexapod.Wang et al. [10] used minimal linear combinations of error parameters for kinematic calibration of a spatial 3-DOF PKM in a fiveaxis HKM milling machine.Wang et al. [11] applied the Markov Chain Monte Carlo method for kinematic calibration of an HKM.Liu et al. [12] proposed the use of the genetic algorithm for kinematic calibration of a Stewart platform while Fan et al. [13] proposed it for an HKM polishing machine.
While researchers [14,15] investigate the kinematic calibration of five-axis SKMs, this paper discusses the kinematic calibration of a five-axis HKM.This paper uses an external kinematic calibration of an HKM machine tool utilizing two different identification algorithms, namely a nonlinear least squares algorithm and an iterative linear least squares algorithm.In particular, a subset selection approach is proposed to overcome the rank deficiency in the nonlinear least squares algorithm.Furthermore, an iteration is introduced to be used in the linear least squares to achieve a solution convergence.This scheme is the first novelty of this paper.The implementation of the two algorithms as well as a comparison of the estimation accuracy provided by both algorithms is presented.This comparative study is another contribution of this paper.The HKM machine tool was built by conjugating two lower-DOF, planar PKMs having a similar topology but different geometric parameters.The planar state of the PKMs as will be shown results in a configuration degeneracy which indicates unidentified parameters in the off-plane direction of the mechanism.A physical adjustment should be conducted to impose a certain value, typically zero for convenience, of the unidentified parameters.The implementation of the mentioned algorithms to this new topology of HKM composed of planar PKMs is another novelty of this paper.
Laser tracker, which is capable of providing a full pose measurement, was used to conduct the kinematic calibration.The laser tracker offers a high accuracy pose measurement with easy use, although it is relatively expensive [16].It was mentioned earlier that a minimum of ten calibration poses should be used while a full pose measurement can reduce the required number of calibration poses.Accordingly, a full pose measurement of twelve unique poses is taken in this work.The unrepeated coordinate values for each DOF of the poses provide the uniqueness of the poses.The geometric parameters of the PKMs, which include the length of all the legs and the moving platform as well as the position of the joints connecting the legs with the moving platform, and the actual relative position between the two conjugated PKMs are the parameters to be calibrated.The model function to be used in the calibration is the forward kinematics of the machine's mechanisms.Finally, the accuracy improvement was evaluated after performing the compensation.
This paper is organized into several sections.Section 1 gives an introduction to the paper and describes its novelty.Section 2 provides the kinematics of the individual mechanism and subsequently presents the solution of its forward kinematics by using both Euler angles and quaternions representations.The quaternions are used to overcome the formulation singularity.Section 3 presents the hand-eye and robot-world calibration which should be performed before the calibration of the mechanism geometric parameters.Section 4 presents the kinematics of the complete mechanism with the external measurement device in the setup.Section 5 presents the use of nonlinear least squares to refine the hand-eye and robot-world calibration and to estimate the mechanism geometric parameters.Section 6 presents the use of iterative linear least squares to estimate the mechanism geometric parameters.In Sect.7, the estimates of the mechanism geometric parameters obtained by both the nonlinear least squares and the iterative linear least squares are used for compensation.The accuracy of the position tracking corresponding to the estimates obtained by the two algorithms is compared.Finally, Sect.8 concludes the paper.

Kinematics Equations
The HKM calibrated in this work is used for five-axis machining and composed by the third scheme of hybridization.A PKM which serves as the workpiece platform (lower platform) is conjugated with another PKM as the spindle platform (upper platform).Both the lower and upper platforms use planar PKMs because they are easier to design than spatial PKMs and more geometrical constraints on the DOFs will minimize the error due to joints' errors.The lower and upper platforms use an identical PKM topology called the 3PRR, planar PKM as depicted in Fig. 1a and optimized to a simpler topology as shown in Fig. 1b.This PKM has three DOFs which adequately enable required general motion consisting of translation in X-axis, translation in Y-axis, and rotation about the Z-axis, where the three axes are orthogonal and follow the right-hand rule.It consists of a fixed base, three legs, a moving platform, and several joints.The joints used in an order starting from the fixed base are prismatic (P) joints, revolute (R) joints which connect the sliders with the legs, and other revolute (R) joints which connect the legs with the moving platform.Therefore, the mechanism is called the 3PRR PKM.The prismatic joints are actuated, have translation along X-axis, and are implemented by using sliders moving along the base through guideway.The use of the sliders creates more workspace and enables the moving platform to move upward and downward as well as to tilt.The actuation is provided by linear motors which give more accuracy compared with pneumatic/hydraulic pistons and lead screws.
In this proposed mechanism, three legs are required to constrain the mechanism fully.Furthermore, two neighboring legs are maintained in a crossing configuration to avoid singularity as well as maximize the useful workspace and the stiffness of the mechanism.Point (x, y) in the moving platform is used to evaluate the general motion of the moving platform.Any other point in the moving platform can also be used as the moving platform is rigid.The selection of point (x, y) is mainly for more convenience.Besides, the orientation of the platform is defined as the angle made by the platform with respect to the fixed base.Upon the dimensional optimization of the mechanism by considering its workspace, stiffness magnitude, and stiffness uniformity, the mechanism is simplified to a topology with two coincident upper joints as shown in Fig. 1b.As a result, x p2 = x p3 = x p .The optimization procedure is discussed in [17].Figure 2 shows the hybrid scheme of the machine.The conjugation of the two planar PKMs in perpendicular directions easily creates the required spatial workspace and provides six DOFs.
The tool spindle is oriented horizontally to avoid redundancy between the rotation of the spindle and the rotation of the upper mechanism.This results in the proposed five-axis machine configuration.The five-axis machine has translational DOFs in x, y, and z directions as well as rotational DOFs about the Z and Y axes.Moreover, there is a redundancy between the two mechanisms in the X translational motion.This redundancy is advantageous as it provides an additional workspace in the x-direction.However, this also results in more complexity in determining which mechanism to move when a translation in the x-direction is required.In this case, one can pre-determine which mechanism is fixed while another mechanism is moved.A more systematic approach is a determination using an optimization scheme based on a defined objective to be optimized, such as energy consumption (control effort), accuracy, or stiffness.
The kinematics equations of the proposed parallel mechanism can be obtained easily by the algebraic approach.According to the geometry of the mechanism shown in Fig. 1, the following geometric relations are found for the legs: (1) Fig. 2 The prototype of the hybrid kinematics machine In this paper, only the forward kinematics of the mechanism is discussed as it is required in the kinematic calibration.The kinematics of the mechanism has the following constraints: where y max and x max are the vertical and horizontal limits of the machine volume having the numerical values of 1200 mm and 1700 mm, respectively.Those constraints indicate that the mechanism should always be inside the specified machine volume, the sliders should not interfere with each other, and two adjacent legs must always be crossing each other.

Forward Kinematics with Known θ
The complete forward kinematics of the mechanism at hand which includes the determination of the angle θ can be seen in [17].In this paper, only the determination of x and y given the angle θ is presented for the kinematic calibration purpose.In the case of a known θ, the value of x and y can be obtained through a more straightforward derivation.Let us rearrange the kinematics Eqs.(1), (2), and (3) in the following set of equations: Collecting x and y in ( 5)- (7) in the left-hand sides yields: (2) Furthermore, subtracting (9) from ( 8) and ( 10) yields a more compact set of equations as follows: where Equations ( 11) and (12) have to be solved simultaneously as x and y exist in both of the equations.The solution is given by:

Forward Kinematics Using Quaternion
Observing the solution of the forward kinematics given in ( 13) and ( 14), we realize that division by zero occurs when θ = 0 and therefore b = 0.In this situation, x still can be solved by using (13) by converting zero θ to a minimal value near (above or under) zero.The effect of the closeness of θ to zero can be evaluated by observing the value of the solution, i.e. x and y, as well as the corresponding links' dimension if we apply the solution.The latter shows a more obvious effect.After scanning some values of θ close to zero, it is shown that θ should not be larger than 1 × 10 −5 .This will give an accurate solution to four decimal places.Using a closer value of θ to zero increases further the accuracy.Table 1 shows the values of the link lengths corresponding (9) to different values of θ near zero, given the expected actual values of L 1 = 600 mm, L 2 = 600 mm, and L 3 = 750 mm.On the other hand, division by zero in ( 14) is very sensitive.Even a minimal value of θ replacing an exact zero θ leads to a significant error in the solution, i.e. the value of y.Therefore, ( 14) cannot be used to solve for y in the case of zero tilting angle of the moving platform.In such a case, either one of the following two ways can be used.First, using an iterative technique by spanning over possible values of y and checking with the inverse kinematics.For any set of x, y, and θ, it should be observed if it retrieves a known set of x 1 , x 2 , and x 3 .However, the accuracy of this technique is affected by the interval of y used in the iteration.Besides, this technique is computationally intensive.Second, solving the forward kinematics by using quaternion.The reason for using the quaternion is because a formulation singularity causes this problem due to the use of Euler angle θ.Therefore, avoiding using the Euler angle by using the quaternion will avoid the formulation singularity.
Recall the mechanism at hand defined by some position vectors as shown in Fig. 3.The geometric relations in the mechanism can also be written using a vectorial approach as follows: (15) Using Euler angles, (15) can be expanded to the following: From ( 16), we know that we have six equations with nine variables, i.e. x, y, θ, α 1 , α 2 , α 3 , x 1 , x 2 , and x 3 .However, three variables among the nine variables are known.In other words, we have six equations with six unknowns.
A quaternion can be written as: where i, j, and k are the unit vectors in x, y, and z directions, respectively.
A rotation in a three-dimensional space can be represented by a rotation matrix R in terms of quaternions as follows: The rotation of the legs and the moving platform of the mechanism, which is about the Z-axis, can be written in the quaternion notation in terms of the Euler angles α 1 , α 2 , α 3 , and θ as follows: Substituting any of ( 19)-( 20) into (18) as they have a similar form to get the corresponding rotation matrix, the rotation about the Z-axis can be written in a rotation matrix in terms of the quaternions as follows: (16) Table 1 The link lengths corresponding to values of θ near zero when solving for x L 1 (mm) 600.0025 600.0003 600.0000L 2 (mm) 600.0000 600.0000 600.0000L 3 (mm) 750.0040 750.0004 750.0000 Fig. 3 The position vectors defined in the mechanism Now we can express the geometric relations ( 16) by utilizing the rotation matrices in terms of the quaternions ( 21): Furthermore, we can write the unit norm constraints of the quaternions as follows: Since in this case: then we can rewrite ( 28)-( 29) as follows: To this point, we can see that we have ten equations with thirteen variables, i.e. x, y, q 0 (α 1 ), q 3 (α 1 ), q 0 (α 2 ), q 3 (α 2 ), q 0 (α 3 ), q 3 (α 3 ), q 0 (θ), q 3 (θ), x 1 , x 2 , and x 3 .However, three variables among the 13 variables are known.In other words, we have ten equations with ten unknowns.
For ease, x can be computed using (13).When θ is zero, we can replace it with a small value near zero to accurately compute x.Once θ and x are obtained, we can use any pair of ( 22)-( 23), ( 24)-( 25), or ( 26) - (27).Using ( 22)-( 23) is easier since we do not involve q 3 (θ) for which we need to convert θ to get q 3 (θ).From ( 22), we have: Substituting (33) into (31), we obtain the following: Now we can get y in terms of the knowns by substituting (34) into ( 23): It is worth mentioning that (35) applies for all values of θ, not only for θ equal or close to zero.Therefore, it is more practical to use (35) instead of ( 14) to compute y for all cases.

Hand-Eye and Robot-World Calibration
The kinematic calibration presented in this paper utilized a Leica absolute laser tracker in combination with a T-Mac TMC30-F reflector to measure the platform poses.The laser tracker head is automatically controlled to track the pose of the reflector rigidly attached on the platform.Therefore, the laser tracker measures the pose of the reflector instead of the platform.The use of a similar laser tracker for the external kinematic calibration of a serial kinematics mechanism can be found in [18].The laser tracker in combination with the reflector used in this work is capable of measuring all rigid body DOFs which include both position and orientation, i.e. three translational DOFs and three rotational DOFs. Figure 4 shows the local frame of the reflector.
As shown in the figure, the origin of the reflector frame takes place at the center of the reflector prism whereas the three axes of the reflector frame are aligned with the geometry of the reflector.The accuracy of the laser tracker is ± 15 μm for the position and 0.01 degree for the rotation.
Since the laser tracker is an external device, the reflector pose is provided with respect to the measurement frame fixed to the laser tracker frame.Figure 5 depicts the transformations among the machine base frame F B , the platform frame F P , the measurement frame F M , and the reflector frame Fig. 4 The T-Mac TMC30-F reflector and its local frame 1 3 F r .T 1 is the transformation from the machine base frame to the platform.This is given by the forward kinematics of the mechanism.T 2 is the transformation from the measurement frame to the reflector frame.X 1 is the transformation from the reflector frame to the platform, whereas X 2 is the transformation from the measurement frame to the machine base frame.The platform or end-effector is commonly called a hand whereas the reflector is called an eye.A calibration to estimate the transformation X 1 is accordingly called the hand-eye calibration.The measurement frame is also commonly called the world frame and consequently a calibration to estimate the transformation X 2 is usually called the robotworld calibration.Since the laser tracker at hand is capable of providing the six-DOF pose of the reflector, T 2 is wholly known.Unfortunately, X 1 and X 2 are unknown.As a result, they should be estimated before calibrating the mechanism.
To solve for X 1 and X 2 , at least three different poses of the mechanism should be measured and T 1 should be available.In this case, T 1 is obtained based on the nominal values of the mechanism geometry.Hence, the estimation of X 1 and X 2 assumes that the nominal values of the mechanism geometry are good enough.These nominal values are to be calibrated later based on the obtained X 1 and X 2 .The estimation of X 1 and X 2 can be classified into two broad categories: separable (sequential) and simultaneous techniques.The former techniques include [19][20][21][22][23] whereas the latter techniques include [24][25][26][27][28][29][30][31].In the former techniques, X 1 is first obtained and subsequently, X 2 is computed based on the obtained X 1 .This approach is more straightforward but gives a less precise solution.On the other hand, the latter techniques estimate both X 1 and X 2 simultaneously.In this work, the former approach is applied only to obtain the initial values for the next algorithm proposed to refine the solution.
Let T 1 , T 2 , X 1 , and X 2 be homogeneous transformation matrices which contain both rotation and translation components.The estimation of X 1 can be derived based on the following relation among the transformation matrices: For a pair of mechanism poses as illustrated in Fig. 6, we have the following: It can be observed from Fig. 6 that X 1 and X 2 are constant, whereas T 1 and T 2 change with the mechanism pose.
Upon manipulating (37a) and (37b), we have: where Notice that A and B are homogeneous transformation matrices.R A and R B are rotation matrices whereas t A and t B are position vectors.Equation (38), commonly known as AX = XB problem, is the hand-eye calibration problem.
A separable closed-form hand-eye calibration following [19] is employed here.Instead of describing the hand-eye calibration procedure in a derivation format, the procedure is described briefly here in a sequential manner as the following: Step 1 Find the transformation T 1 of two different poses.This is obtained by the forward kinematics of the mechanism using the nominal values of the mechanism geometry, namely L 1 , L 2 , L 3 , x p2 , and Step 2 Find the transformation T 2 of the two poses.This is given directly by the measurement of the full pose of the reflector by using the laser tracker.In this work, the measurement of the reflector orientation is provided in quaternions.
Step 3 Compute the matrices A and B for the two poses by using (39) and (40).
Step 4 Convert the rotation matrices R A and R B into an axis-angle representation defined by a rotation angle α and a rotation axis (n 1 , n 2 , n 3 ).
Step 5 Define PA and PB as follows: Step 6 Define υ and S as follows: Step 7 We have the following linear system: With only a pair of poses, this linear system is singular.Therefore, at least two pairs of poses, i.e. three different poses, are required.Perform Step 1 to Step 7 for all n poses and accordingly stack (45) from all pairs of poses to compose the following overdetermined system: Solve for PX 1 by using linear least squares: Step 8 Compute PX 1 as follows: Step 9 Convert PX 1 to the rotation matrix of the transfor- mation X 1 , namely R X 1 : where I is an identity matrix whereas skew PX 1 is a skew-symmetric matrix defined in a similar fashion to (44).To this point, we solve for the rotation matrix of the transformation X 1 .
Step 10: Define a linear system involving the translation vector of the transformation X 1 , namely t X1 , as follows: Solve the linear system (50) by using linear least squares: Equations ( 50) and (51) hold in general for three-dimensional Euclidean space.Since the mechanism at hand is planar, the dimensions of (50) and (51) should be reduced to a two-dimensional Euclidean space.Having the mechanism on the XY plane, the elements corresponding to the Z-axis in (50) and (51) should be suppressed.If it is not suppressed and accordingly x, y, and z elements of t X1 are all to be solved, the system will be rank deficient.This indicates a configuration degeneracy.Physically, this means that there is no reference in the z-direction to which the z element of t X1 should be defined.As a result, an infinite number of possible values of z exist.
Equations ( 49) and (51) completely define the transformation matrix X 1 : where, in the case of planar mechanism working on the XY plane: Equation (53) also implies that the zero coordinate of the Z-axis is aligned with the reflector.After the hand-eye transformation X 1 has been estimated, the robot-world transformation X 2 can be easily estimated by rearranging (36): Using n different mechanism poses in the measurement, we accordingly have n different T 1 and T 2 .As a result, there will be n different X 2 .To obtain a single X 2 , which is expected to be constant, we average the elements of T 1 and T 2 in an element-wise manner to get an averaged, constant X 2 .
It is worth mentioning that the selected poses should be completely different from each other in all coordinates in order to get the best estimates.In this case, the x, y, and θ values should be completely different among the poses.In this work, twelve experimental measurement data was carefully taken in which the x, y, and θ values of the poses are entirely different.

Kinematics of the Combined Mechanism with the External Measurement Device
The pose measurement is always provided with respect to the measurement frame.Since the kinematic calibration uses an external pose-measurement device, an extended kinematic model should be established to include the handeye and robot-world transformations as described earlier.A homogeneous transformation is very convenient to be used to represent this kinematics.Figure 7 shows all the homogeneous transformation matrices in the combined mechanism.For convenience, a similar matrix notation is used for the upper mechanism.The asterisk superscripts indicate transformation matrices belonging to the upper mechanism.
The transformation among the laser tracker, the reflector, and the lower mechanism has been described earlier in (36).Similarly, the homogeneous transformation of the laser tracker, the reflector, and the upper mechanism can be written as: Alternatively, the following homogeneous transformation can also be used: where X 0 denotes the homogeneous transformation from the machine base frame to the base frame of the upper mechanism.Using (56), we have: Another way to express the kinematics involving the laser tracker and the reflector is by using vector notation.Figure 8 shows the position vectors in the combined mechanism.The superscript indicates the frame in which a position vector r is expressed, whereas the subscript indicates a point defined by the position vector.The subscripts and superscripts M, B, O, P, and r respectively represent the measurement frame, machine base frame, base frame of (55) Fig. 7 Homogeneous transformations in the combined mechanism with a laser tracker and reflectors, namely T 1 , T 2 , X 1 , X 2 (for the lower mechanism), T * 1 , T * 2 , X * 1 , X * 2 (for the upper mechanism), and X 0 which connects both mechanisms The additional subscripts L and U indicate the lower and upper mechanism, respectively.The transformation between frames may involve a rotation matrix R which represents the rotational transformation between two frames.Similarly, a notation R M B , as an example, denotes the rotation matrix from the machine base frame B to the measurement frame M.
For more clarity, mainly when dealing with rotational transformation, Fig. 9 shows only the orientation of all the frames.The drawn axes indicate the positive directions of the axes.The orientation of the axes shown in this figure follows the right-hand rule in which the cross product between the X and Y axes gives the direction of the Z-axis.For more convenience, the base frame of the upper mechanism is oriented such that its X and Y axes create the planar workspace of the upper mechanism.Accordingly, the positive Z-axis is pointing downward.Using such an orientation, the inverse and forward kinematics of both the lower and upper mechanisms can be written on the XY plane.
Since the mechanisms are planar, by referring to Fig. 6, their rotation about their base frames can be represented by an elementary rotation matrix about the Z-axis, i.e.: (59) where R z is an elementary rotation matrix about Z-axis whereas θ L and θ U are the rotation angles of the lower and upper mechanisms, respectively.The reflectors are attached to the platforms, as shown in Fig. 6, such that: where I and R x are an identity matrix and an elementary rotation matrix about the X-axis, respectively.This means that the lower reflector frame is aligned with the lower platform frame whereas the upper reflector frame is also aligned but rotated about the X-axis by -π/2 rad.This arrangement can be made with the aid of measurement devices such as digital calipers, digital depth gage, and precision square.
The position of the reflector mounted on the lower platform can be expressed in the measurement frame as follows: where R M B denotes the rotation matrix from the machine base frame to the measurement frame whereas R B r,L denotes the rotation matrix from the reflector frame to the machine base frame.
Since the reflector frame is aligned with the platform frame, i.e. the orientation transformation between them is represented by an identity rotation matrix R P,L r,L , the rotation of the reflector frame in the measurement frame is given by: The position of the platform end-effector in the base frame r B P,L is given by the forward kinematics, i.e.: where f 1 and f 2 are the forward kinematics equations and P L contains all the geometric parameters of the lower mechanism.
For convenience, the Z-axis of the machine base frame is assumed to be aligned with the Z-axis of the reflector frame.In other words, the z coordinate of the reflector frame origin with respect to the machine base frame is zero: Similarly, the position of the reflector mounted on the upper platform can be expressed in the measurement frame as follows: The rotation of the reflector frame in the measurement frame, which is measured by using the laser tracker, can be written as: Solving for R O P,U from (71) yields: Because we have (62), then: Substituting (73) into (70), we obtain: In fact, both the position of the reflector r M r,U and the orientation of the reflector R M r,U , expressed in the measurement frame, are measured by using the laser tracker.The rotation matrix R M O and the position vector r M O can be estimated by using the robot-world calibration, whereas the position vector r P,U r,U can be estimated by using the hand-eye calibration.The position of the platform end-effector in the base (68) frame of the upper mechanism r O P,U is given by the forward kinematics, i.e.: where f 1 and f 2 are the forward kinematics equations and P U contains all the geometric parameters of the upper mechanism.
The transformation from the machine base frame B to the base frame of the upper mechanism O can be obtained from the following relations: From ( 76) and (77), we obtain: Figure 10 shows the laser tracker fixed sitting beside the prototype machine.Figure 11 depicts the reflector mounted on the lower and upper platforms.Since only one reflector is available, the reflector was first mounted on the lower platform and subsequently mounted on the upper platform, or vice versa.In each setup, the estimation of the hand-eye transformation, robot-world transformation, and the geometric parameters of the lower mechanism was conducted.As mentioned earlier, the reflector should (75) Fig. 10 The laser tracker with a reflector installed on the machine be installed to satisfy (61) and (62).To verify this alignment, one should observe the estimated rotation matrix of the reflector frame with respect to the platform frame, namely R X 1 ,L = R P,L r,L and R X 1 ,U = R P,U r,U .After the estimation was conducted for both the mechanisms, the transformation from the machine base frame B to the base frame of the upper mechanism O can be estimated.

Nonlinear least squares algorithm
In this section, Gaussian least squares differential correction (GLSDC) algorithm [32][33][34] (also called the Gauss-Newton nonlinear least squares) is applied for the refinement of the hand-eye calibration as well as the calibration of the mechanism geometric parameters.Using this algorithm, the nonlinearity of the system is taken into consideration.
Let P and Y be the parameters to be estimated and the model function, respectively, then the measured function values are given with a noise υ as follows: The noise υ is typically Gaussian with zero mean and a standard deviation is given by the uncertainty of the measurement device.
The model function Y for the lower mechanism is given by (67), i.e.: whereas the model function Y for the upper mechanism is given by (74), i.e.: The residual error is given by: Notice that the tilde (~) mark on a parameter indicates a measured parameter whereas a hat mark indicates an estimated parameter.Given m measurements, the measured function values Ỹ , the estimated function values Ŷ , and accordingly the residual errors e should be stacked in a single vector.
For the calibration of the lower mechanism, the residual error is defined as: The measured position vector rM r,L is provided by the laser tracker, whereas the estimated position vector rM r,L is evaluated by using (67) or (81).Similarly, the residual error for the calibration of the upper platform is defined as: While the laser tracker provides the measured position vector rM r,U , the estimated position vector rM r,U is evaluated by using (74) or (82).
Besides the position vectors rM r,L and rM r,U , the other parameters given by measurements are the orientation matrices R M r,L and R M r,U , as well as the active joint positions which contribute to the forward kinematics solution of both mechanisms.
The estimation is conducted to find the estimates of the parameters P which minimizes the following cost function: and accordingly minimizes the residual error in (83)-( 85).The cost function F in (86) indicates the norm or the square of the residual error.The matrix W denotes the weight matrix which is, in this work, chosen to be an identity matrix.Using GLSDC, the parameter estimates are computed iteratively as follows: where and the subscript k indicates the k-th iteration.
A modification to (88) called the Levenberg-Marquardt algorithm (also called the damped least squares) had been proposed to speed up the convergence in the case of initial values far from the "true" values.Since the initial values assigned in this work are considered close enough to the "true" values, the GLSDC is sufficient to give a fast convergence.
The Jacobian matrix H is obtained by differentiating the model function with respect to the estimated parameters: where Y x , Y y , and Y z denote the x, y, and z components of the model function Y.
Given m measurements, the Jacobian matrix H should be stacked as follows: The squared system Jacobian matrix in each iteration is the inverted part in (88), i.e.: The matrix J should have a full rank in order to give a trusted solution.This full rank indicates that all the parameters are independent and therefore fully identified.In other words, a rank deficiency by n indicates that n parameters are dependent and therefore unidentified.The determination of the unidentified parameters can be done mathematically or through knowledge on the parameter dependency of the physical system.In a system with a rank deficiency, the dependent parameters should be eliminated until a full-rank system is obtained.In such a case, the nominal parameter values can be used.The dependent parameters can be estimated in the next step employing the estimates of the independent parameters having been obtained previously.This is commonly called the subset (87) selection which results in a sequential estimation, i.e. the estimation is conducted in several steps.This approach assumes that the nominal parameter values should be good enough.Figure 12 shows how the subset selection is implemented in the GLSDC algorithm.As shown in the flowchart, if J does not have a full rank N then only M parameters are firstly estimated where M < N is the rank of J. This, as will be discussed soon, includes the following sequential estimation stages: The iterations in the GLSDC keeps running until the norm of ∆F is less than a specified threshold.This is the stopping criteria of the GLSDC algorithms.In this work, the norm (∆F) of 1.0 × 10 −10 is used.

Refinement of the Hand-Eye and Robot-World Calibration
After the procedure described in Sect. 1 has been performed to obtain the rough estimates of the hand-eye and robot-world transformations, a refinement is conducted by using the obtained estimates as initial values in an iterative estimation using GLSDC.For the lower mechanism, since the orientation of the reflector frame is aligned with that of the platform frame, the remaining parameters which should be estimated to establish the hand-eye and robotworld transformations are the position and orientation of the machine base frame with respect to the measurement frame, namely r M B and R M B , as well as the position of the reflector with respect to the lower platform frame, namely r P,L r,L .As the rotation matrix R M B is defined in terms of quaternions, namely q 0 , q 1 , q 2 , and q 3 , the position vector r M B is defined by the components x B , y B , and z B , and the position vector r P,L r,L is defined by the components x r and y r while z r = 0, the parameters to be estimated are: Similarly for the upper mechanism, since the orientation of the reflector frame is aligned with that of the platform frame with a rotation of -π/2 about X-axis, the remaining parameters which should be estimated to establish the handeye and robot-world transformations are the position and orientation of the machine base frame with respect to the measurement frame, namely r M O and R M O , as well as the position of the reflector with respect to the lower platform frame, namely r P,U r,U .As the rotation matrix R M O is defined in terms of quaternions, namely q 0 , q 1 , q 2 , and q 3 , the position vector r M O is defined by the components x O , y O , and z O , and the position vector r P,U r,U is defined by the components x r and z r while y r = 0, the parameters to be estimated are: The evaluation of the squared system Jacobian matrix J for both the lower and upper mechanisms shows that its rank is 9 (full rank) and therefore all the parameters in (93) and (94) can be estimated in a single step for each mechanism.
(93) P = q 0 q 1 q 2 q 3 x B y B z B x r y r T (94) P = q 0 q 1 q 2 q 3 x O y O z O x r z r T Tables 2 and 3 show both the estimates obtained using the closed-form solution as described in Sect. 1 and those obtained through the refinement using the GLSDC, for the lower and upper mechanisms, respectively.
As a partial benchmark, the available nominal value of y r is 185 mm.Since the reflector is mounted well on the lower platform, the deviation from the nominal value should not be too large.The estimated value of y r = 186.087565mm is considered acceptable.This indicates that the refinement works well.
After the hand-eye and robot-world transformations for both the lower and upper mechanisms have been performed, the transformation from the machine base frame B to the base frame of the upper mechanism O can be easily computed by using (78) and (79).The obtained rotation matrix is: whereas the position vector from B to O is: This shows that the origin of the X-axis of the base frame of the upper mechanism O is shifted by − 2.084870 mm from that of the machine base frame B. In the design, x B O is expected to be zero.However, the estimated value shows that it is not zero.This occurs because there may be an error in the manufacturing of the machine base structure.

Estimation of the Mechanism Geometric Parameters
For each mechanism, there are five geometric parameters to be estimated, namely l 1 , l 2 , l 3 , x p2 and x p3 .In the case of x p2 = x p3 = x p , which is the case of the prototype machine, there are four geometric parameters to be estimated, i.e. l 1 , l 2 , l 3 ,and x p .Therefore, the parameters to be estimated are: The observation shows that the squared system Jacobian matrix J gives only a rank of 4 in the case of five parameters to be estimated, i.e. when x p2 ≠ x p3 , and a rank of 3 in the case of four parameters to be estimated, i.e. when x p2 = x p3 = x p .This is because one of all the geometric parameters to be estimated is dependent on the other four parameters.In other words, all the geometric parameters but one are known then the one parameter remaining can be calculated by using the closed-form geometric relation of the mechanism.
Another approach is to estimate all the geometric parameters in two steps.In the first step, L 2 , L 3 , x p2 , and x p3 are estimated.In this step, the nominal value of L 1 is used in the estimation iterations.In the second step, the remaining parameter, i.e.L 1 , is estimated.Although the nominal value of L 1 is used in the first step, L 1 is subsequently estimated in the second step.Using twelve unique calibration poses for each of the lower and upper mechanisms, the estimates of the geometric parameters of the lower and upper mechanisms along with their nominal values are shown in Tables 4 and  5.The errors are defined as the difference between the estimates and the nominal values.It can be observed in Tables 4  and 5 that the estimate of L 1 is similar to its nominal value. (97 This is because the nominal value of L 1 is used in the first step and subsequently it is estimated after the other parameters have been estimated.In other words, the nominal values of L 1 are retrieved in this two-step estimation.

Estimation of the Geometric Parameters
Using Iterative Linear Least Squares In this section, a linear error model is used for the kinematic calibration.Accordingly, the closed-form linear least squares estimation can be employed.From the kinematic relations given in ( 1)-( 3), a small perturbation can be introduced to all the parameters which include the geometric parameters P, the active joint positions q, and the platform pose X.This small perturbation which represents a small error is a first-order approximation.Following ( 1)-( 3), the perturbed kinematics can be written by differentiating the kinematics with respect to all parameters: In a matrix form, (98)-(100) can be written compactly as: where (98) + (y + x p3 sin( ))( y + x p3 sin( ) + x p3 cos( ) ) (101) J x X + J q q = J p P (102) Let all the terms in the left-hand side be called A, i.e.: and B = J p , then we can rewrite (101) as: Having m measurements, we can stack (110) into the following overdetermined system: where The estimation of P is aimed at minimizing the fol- lowing cost function: As the problem is an overdetermined system, a closedform solution for δp is given by the following linear least squares: Since (115) denotes the errors in the geometric parameters, the estimated geometric parameters are given by: where P and P denote the estimated and nominal values of the geometric parameters, respectively.
Assuming that the measurement is given with respect to the machine base frame, the following is the sequential procedure to implement the linear least squares estimation based on the abovementioned linear error model: Step 1 Assign the poses to be visited, X.
Step 2 Measure the poses, X measured .
Step 4 Perform the inverse kinematics to obtain the active joint positions q = f(X,P).
Step 5 Measure the active joint positions, q measured .
Step 8 Compose the vector B(X,q,P).
Step 9 Repeat Step 1 to Step 8 for m measurements.Compose Ā and B as in (113).
Step 10 Check the rank of BT B.
Step 11 If BT B has full rank, then compute the linear least squares solution (115).Otherwise, eliminate the dependent parameters either mathematically, such as through SVD, or based on knowledge on the physical system.
In real practice, the measured poses X measured is often given by an external measurement device such as a camera or a laser tracker as in this work.As a result, X measured is given with respect to a measurement frame which is not aligned with the machine base frame.Accordingly, the assigned poses X should be defined with respect to the measurement frame.The values of X with respect to the measurement frame can be derived from (67) or (81) for the lower mechanism and from (74) or (82) for the upper mechanism.
From (67) or (81), we obtain: In a similar fashion, from (74) or (81) we have: The sequential procedure to implement the linear least squares estimation based on the linear model when an external pose measurement device is used can be summarized as follows: Step 1 Estimate the hand-eye and robot-world transformations.
Step 2 Assign X with respect to the machine base frame F B for the lower mechanism and with respect to the base frame of the upper mechanism F O for the upper mechanism.
Step 3 Measure the reflector position with respect to the measurement frame, namely r M r,L for the lower mechanism and r M r,U for the upper mechanism.
Step 4 Calculate X measured by utilizing (117) for the lower mechanism and (118) for the upper mechanism. (117) Step 5 Perform Step 3 to Step 11 as described for the pose measurement in the machine base frame.
Using the same measurement data as used in the GLSDC, the resulting overdetermined linear system (112) with all the geometric parameters to be estimated has a full-ranked BT B .Therefore, the estimation can be performed straight- forwardly in a single step.However, the estimates given by a single run of the algorithm are not satisfying.They are inconsistent and not guaranteed to be the optimum solution.Some coordinates may have an error less than that of the uncalibrated one, but some other coordinates may have an even more significant error than that of the uncalibrated one.This implies that the estimation algorithm fails to minimize the cost function, i.e. the residual errors.
To overcome this problem, the linear least squares algorithm as described above needs to be iterated until a minimum norm of the residual errors ε is achieved.Since this iterative linear least squares algorithm requires higher computational cost compared to the iterative nonlinear least squares, one can also stop the algorithm after a certain amount of time or a certain number of iterations and evaluate the obtained estimates.The flowchart of the iterative linear least squares is depicted in Fig. 13.It can be seen in the flowchart that the estimated parameter values are updated and subsequently supplied to the new linear system which will be solved iteratively until the stopping criteria as discussed is achieved.Furthermore, the norm of the residual is evaluated between any two consecutive iterations.If the norm of the residual is increasing, the new parameter values should be subtracted by the parameter errors.Otherwise, they should be summed.Using twelve unique calibration poses for each of the lower and upper mechanisms, the iterative linear least squares (124) provides the estimated geometric parameters of the lower and upper mechanisms as shown in Tables 6 and  7, respectively.In a similar manner with the two-step estimation using the nonlinear least squares, the errors are defined as the difference between the estimated and nominal values.
To this point, it still cannot be judged whether the estimates of the geometric parameters obtained by using the nonlinear least squares or those obtained by the iterative linear squares are more accurate as the true values of the geometric parameters are unknown.To evaluate and compare the accuracy of the estimates obtained by using both the algorithms, the pose errors will be evaluated in the next section upon the compensation of the kinematic parameters.

Compensation
After the estimation of the geometric parameters has been done, compensation should be conducted to improve the accuracy of the machine.This is performed by replacing the nominal geometric parameter values in the kinematics by the estimated ones.Since two estimation techniques have been implemented, the estimates from both the techniques  were used for the compensation.To evaluate the accuracy of the mechanism pose, pose errors are defined for all coordinates of the pose, i.e. x, y, and θ at the twelve different poses.Figure 14 shows a comparison between the pose errors of the uncalibrated and calibrated lower mechanism at the twelve different poses.For conciseness, the pose errors of the uncalibrated and calibrated upper mechanism are not shown in this paper as their behavior is similar to those of the lower mechanism.The pose errors are presented for each coordinate of the pose, i.e. x, y, and θ.It can be seen that the both the iterative nonlinear least squares (NLS, i.e.GLSDC) and the iterative linear least squares (LLS) algorithms successfully suppress the pose errors to values very close to zero.Some of the plots in this figure cannot show clearly the difference between the pose errors corresponding to both the algorithms as the error curves look coincident although they have a small difference.It is shown that the average position accuracy before the calibration is around 0.9 mm whereas the average orientation accuracy is less than 0.2 degree.The calibration using GLSDC provides an average accuracy of less than 0.004 mm for the position and less than 0.0002 degrees for the orientation, whereas the iterative linear least squares algorithm provides an average accuracy of less than 0.01 mm for the position and a similar average accuracy in the orientation to the GLSDC.In general, both the calibration algorithms can improve the position accuracy by around 0.2 mm.To illustrate more clearly, the plots of the pose errors of the lower mechanism corresponding to only the calibrated parameters at the twelve different poses are shown in Fig. 15.The plots show clearly that the GLSDC performs better in the position accuracy than the iterative linear least squares.However, both algorithms do not perform differently in the orientation accuracy.It can be seen in these plots that the iterative nonlinear least squares algorithm outperforms the iterative linear least squares in the pose error suppression.In addition to that, the former algorithm also outperforms the latter in the computational cost as the former converges faster than the latter.The improved accuracy mentioned above is position and orientation accuracy in the case of position tracking.To evaluate the effect of calibration in the accuracy of contour tracking, a test contour was executed using both the nominal and estimated parameters.Since the estimated parameters obtained by using the nonlinear least squares are more accurate than those obtained by the iterative linear least squares, they are used in the test contour.In this work, a full circle was selected as the test contour.Figure 16 shows the contour performed by the lower mechanism and therefore lying on the XY plane.The circle is centered at (400, 300) mm and has a radius of 50 mm.As shown in the figure, the contour corresponding to the nominal parameters was at a glance very close to that corresponding to the calibrated (estimated) parameters.In order to compare the accuracy, contour errors which are defined as the difference between the nominal coordinates and the calibrated coordinates were plotted as shown in Fig. 17.Since the planar contour has two coordinates, i.e. x and y, the contour errors should be evaluated for each coordinate.The figures depict the contour errors along the whole trajectory of the contour.It is shown in the figures that the contour corresponding to the calibrated parameters is more accurate by around 0.2 mm than that corresponding to the nominal (uncalibrated) parameters.This is consistent with the position accuracy of around 0.2 mm provided by both the calibration algorithms.The remaining position accuracy is basically due to the posture and/or dynamics of the mechanism which is a task for the control system to overcome.Along the trajectory of the circle, the mechanism posture keeps changing.As can be seen in Fig. 17, at some postures the error is zero while at other postures the error starts increasing.

Conclusion
It is shown that the nonlinear least squares algorithm successfully refined the hand-eye and robot-world transformation previously obtained by using a simple separable technique.The evaluation of the contouring error also shows that the nonlinear least squares algorithm estimates more accurately the geometric parameters of the machine's mechanisms than the iterative linear least squares algorithm.In return to this higher accuracy, the nonlinear least squares algorithm should be performed in two steps by employing the subset selection approach in which m parameters should be firstly fixed at their nominal values while the rest of the parameters are estimated (where N is the total number of parameters to be estimated, r is the rank of the system Jacobian matrix, and m = N − r).Subsequently, the previously fixed parameter(s) should be estimated in the next estimation step.This implies that the fixed parameter(s) in the first step should have a considerably accurate value.Among all the estimated parameters, the parameter(s) which is believed to have the most accurate value(s) can be chosen as the fixed parameter(s).Such dependency to the nominal values of the parameters does not only occur in this situation; it also occurs in the hand-eye and robot-world calibration which assumes some degree of accuracy of the nominal values of the mechanism geometric parameters.On the other hand, the closedform solution of the linear least squares should be iterated until it converges to the actual values.Otherwise, it may give estimates of the parameters that are worse than the nominal values due to divergence from the actual values.Finally, using the compensated kinematics, it is shown that both the position tracking error and the contour tracking error were reduced.

Fig. 1
Fig. 1 Schematic of the a general and b optimized mechanism topology (with the Z-axis perpendicular to the paper through the origin)

Fig. 8
Fig. 8 Position vectors in the combined mechanism with a laser tracker and reflectors, where the subscripts indicate the points of interest whereas the superscripts indicate the reference frames

Fig. 11
Fig. 11 The reflector mounted a on the upper platform and b on the lower platform

Stage 1 . 4
Estimation of the hand-eye and robot-world transformation parameters to refine the closed-form solution obtained previously.Stage 2 Estimation of the geometric parameters of the lower mechanism.This consists of two steps.Stage 3 Estimation of the position of the base frame of the upper mechanism,r B O Stage Estimation of the geometric parameters of the upper mechanism.Similar to the case of the lower platform, this consists of two steps.

Fig. 12
Fig. 12 Flowchart of the subset selection in the implementation of GLSDC

Fig. 13
Fig. 13 Flowchart of the iterative linear least squares

Fig. 14 Fig. 15
Fig. 14 Pose errors of a x coordinates, b y coordinates, and c angles θ of the lower mechanism at twelve data points (poses) before calibration and after calibration using NLS and iterative LLS

Fig. 16 Fig. 17
Fig. 16 Test contour performed by the lower mechanism (blue: using nominal values, red: using calibrated values).(Color figure online)

Table 2
The estimated hand-eye and robot-world transformation parameters for the lower mechanism

Table 3
The estimated hand-eye and robot-world transformation parameters for the upper mechanism

Table 4
The estimates of geometric parameters of the lower mechanism (x p2 = x p3 = x p )

Table 5
The estimates of geometric parameters of the upper mechanism (x p2 = x p3 = x p )

Table 6
Estimated geometric parameters of the lower mechanism obtained by using iterative linear least squares

Table 7
Estimated geometric parameters of the upper mechanism obtained by using iterative linear least squares