A Novel Optimal Design of Measurement Configurations in Robot Calibration

Robot calibration highly depends on redundantmeasurement configurations to collect enough sample data for higher accuracy, but excessive measurements seem to be uneconomical and time-consuming.Thus lots of observability indexes to evaluate the goodness of the measurement configurations have been proposed. However, in some circumstances, it is of critical importance to obtain accurate kinematic parameters and estimate the end-effector pose precisely at the same time. Obviously one single observability index can hardly meet this need yet. Accordingly, after analyzing the essential constrains of robot calibration and the influence of measurement configurations on the observability indexes, an optimization model with two observability indexes to be taken into consideration is proposed in this paper, and then the existing DETMAX algorithm is modified to seek optimal design of measurement configurations, by adopting a set-constructing method and a set-shrinking method. Much better results have been obtained by simulation study, which implies that the proposed model and the modified DETMAX algorithm perform well in both kinematic parameter identification and pose estimation of the end-effector of the robot.


Introduction
Robot calibration, the purpose of which is to obtain the accurate parameters representing the robot structure and the exact estimation of the end-effector pose, can maintain greater economy and convenience than the way guaranteeing robot accuracy by improving manufacturing and assembly precision.Generally, its procedure should encompass four main steps: kinematic modeling, measurement of the actual pose of the end-effector, identification of the kinematic parameters based on the calibration model, and compensation of the motion error by the employment of the identified kinematic parameters [1,2].Judicious selection of robot configurations used for redundant measurements can guarantee robustness of robot calibration with respect to sensor noise [3].The investigation conducted on the evaluation indexes for robot configurations and their applications are elaborated briefly as follows.
In an attempt to evaluate efficacy of a certain configuration set composed of different configurations for robot calibration, several observability indexes have been proposed.[4] proposed an index known as  1 , which is based on the product of the nonzero singular values of the identification matrix.Maximization of  1 indicates more effect of kinematic parameter errors exerted on the estimate pose error.This is also true for the index  5 , proposed by Sun and Hollerbach [5] and based on the sum of the reciprocals of the nonzero singular values mentioned above.The second index proposed by Driels and Pathre [3], here termed  2 , aims to maximizing the inverse of the condition number of the identification matrix, which is equivalent to the smallest nonzero singular value divided by the largest one. 2 obtains its maximum value at 1 when all singular values are equal.The index  3 proposed by Nahvi and Hollerbach [6] is just the minimum singular value, referring to the worst observability of the parameter error as the criterion, whereas the index  4 [6] is the product of  2 and  3 , and actually there is no significant difference from  2 to  4 .Though both  1 and  5 can represent the overall performance of the identification matrix, by emphasizing the volume of a hyper-ellipsoid whose directions are determined by the singular values, these two indexes may lay more emphasis on 2 Mathematical Problems in Engineering one direction than another.Besides the other three indexes focusing on the partial features of the singular values, the global features may be ignored.Consequently it is necessary to take the partial features and the global ones of the singular values into account at the same time.

Borm and Menq
With the observability index determined, the additional challenge is to find an effective algorithm to construct a configuration set whose observability index is maximal.The algorithm most widely used is DETMAX proposed by Mitchell [7], and there have been lots of applications.Daney [8] investigated the influences of different observability indexes on a Gough platform by DETMAX.Li et al. [9] and Zhou et al. [10] adopted this algorithm for calibration of serial robots, and the calibration effect was improved.However, it is sensitive to local convergence.To circumvent this problem, Zhuang et al. [11] proposed a simulated annealing approach to obtain optimal measurement configurations.And then the genetic algorithm was also adopted to plan the robot calibration experiments [12].Wang et al. [13] proposed the modified simulated annealing algorithm for the polishing robot.These heuristic algorithms improved convergence performance but often led to a low convergence rate.
Accordingly, two main improvements are proposed in this paper for robot calibration.The first improvement is to use the combination of  1 and  3 to achieve a multicriteria optimization.The second one is to modify the DETMAX algorithm for a better design of the measurement configuration set.The remainder of this paper is organized as follows.Essential constraints of robot calibration are presented in Section 2 to conduct later analysis of observability problem.In Section 3, it is demonstrated that maximizing  1 always means an unbiased estimate of the kinematic parameters with smallest variance, while  3 is the right criterion for minimizing the possible error of the end-effector for any configuration, from the experimental design literature.DET-MAX is modified by adopting a set-constructing method and a set-shrinking method in Section 4. Simulation validation is presented in Section 5. Finally, conclusions and recommendations are given in Section 6.

Essential Constraints of Calibration
Without loss of generality, a serial robot with  actuators controlling  DOF (degrees of freedom) is taken as the research object.Essential constraints of robot calibration are investigated in this section to conduct subsequent study on observability problems.Conventions used below are introduced for convenience as follows: (i)  denotes the angle of one actuator of the robot, and thus   denotes the angle of the th actuator, where  = 1, 2, . . ., .
(ii) Θ is a  × 1 vector composed of  joint angles, also known as a robot configuration.
(iii) Φ is a set composed of different robot configurations, called a configuration set.
(iv)  denotes the permissible operating range of all the actuators, and also means a set of all possible robot configurations, called the candidate configuration set pool.
Generally the kinematic model of the robot relates the joint angles Θ * and the end-effector pose Ψ * through a function  * based on the irreducible exact kinematic parameters as illustrated in the following.
The joint angles are not known with exactitude due to the existence of sensor noise.Similarly, defects in manufacturing and assembly result in parameter errors.The difference between the nominal values Ω and Θ and the exact values Ω * and Θ * is assumed to be Ω and Θ, such as Ω = Ω * + Ω, Θ = Θ * + Θ.Moreover, for lack of completeness, the modeling function also leads to errors, set as   .In conclusion, it is the modeling error   , the angle measurement error Θ, and the kinematic parameter error Ω that yield the estimate error Ψ of the end-effector pose.To relate these factors to the estimate error of the end-effector pose, differentiate (1) with respect to all its parameters.
where J Ω is the generalized Jacobian matrix, also called the identification matrix, and J Θ is the Jacobian matrix in terms of a robot configuration Θ.The kinematic parameter error makes the greatest contribution to the estimate error of the end-effector pose, considering that either the modeling error or the measurement error is usually small [14].That is to say, the systematic error of the robot system is much greater than the random error.Hence (2) can be simplified as follows: The number of the independent parameters to describe the pose error Ψ is 6, which is much less than that of the kinematic parameter error Ω, indicating that (3) is underdetermined, especially under only one robot configuration.Thus sufficient configurations have to be measured to establish an overdetermined system.Then the definition of parameters in (3) is further extended as Ψ = [Ψ 1 , Ψ 2 , . . ., Ψ  ]  , J Ω = [J Ω1 , J Ω2 , . . ., J Ω ]  .The configuration Θ is replaced with a configuration set Φ.Here (•)  indicates a transpose of a matrix, and  denotes the amount of the measurement configurations.It is noteworthy that the identification matrix is obtained with the inaccurate parameters Ω and Θ, so that the least squares method should be applied iteratively to guarantee accuracy as shown in the following.
where  denotes the total number of iterations, and (•) + means the Moore-Penrose inverse of a matrix.The estimate pose error is obtained by Ψ = Ψ  − Ψ * , in which the exact values Ψ * are replaced with the measured values Ψ  .The process is iterated until Ω  converges to a small value.The identified kinematic parameters can finally be obtained by (5).Note that the measured values Ψ  is different from the exact values Ψ * for defects of measuring equipment and methods, and the error is set as   = Ψ * − Ψ  .
The repetitive positioning accuracy of the end-effector has been determined once the robot is designed and manufactured, which is the robot's property that cannot be modified, whereas the robot calibration process involving kinematic modeling, external measurement of the endeffector poses, and selection of measurement configurations can be optimized.Four general constraints are summed up for robot calibration, the first two of which are prerequisites for calibration and need to be validated in advance, while the last two are interrelated with each other and can be optimized for a better design of measurement configurations.
Constraint A. Sufficiently high end-effector repeatability is a necessary condition for robot calibration, and in general, the precision of the absolute positioning is lower; for example, for a robot with a repeatable positioning accuracy of ±0.05 mm, the absolute positioning accuracy can only reach 2∼3 mm [15].The factors influencing the end-effector repeatability include sensor noise of joint angles measurement, backlash caused by gear gaps, change of temperature, and other random factors.The sensor noise is the main influencing factor, which has to satisfy the inequality as follows: where ΔΘ is the maximum value of the joint angle measurement error and   denotes the permissible range of the repeat positioning error.For all the reachable robot configurations, the repeatability should meet the accuracy requirements.However, (6) only needs to be verified under enough configurations in practice.
Constraint B. Lots of kinematic modeling methods have been proposed for efficient and succinct description of the relation between the adjoining link coordinate frames, such as the widely used DH (Denavit-Hartenberg) method [16].Yet all the modeling methods do not meet three basic requirements: model completeness, parameter minimality, and model continuity [17]; for example, DH method has only four link parameters, which is not enough to achieve a complete kinematics model, and the effects of machining and assembly errors may be overlooked.This is also true for other modeling methods.
It should be noted that nonsingularity of the identification matrix is essential because a pseudo-inverse of the identify matrix is adopted in (4), such as where det(•) is the determinant of a matrix.In order to accurately identify the kinematic parameters, the redundant parameters should be removed in the process of modeling, ensuring that the identification matrix is full rank.
Constraint C. Since the actual end-effector pose of the robot is expressed by the measurement values Ψ  instead of the exact value Ψ * , the parameter identification error is mainly derived from the external measurement noise under the assumption that the linearization error can be remedied through iteration.Hence external measurement instruments including the laser tracker, the coordinate measuring machine (CMM), and the optical measuring instrument generally have a high measuring accuracy; for example, the laser tracker AT501B produced by Leica has a static positioning accuracy of up to 10 −6 m.
Constraint D. On the one hand, sufficient measurement configurations are selected to establish an overdetermined system of equation.Assume that the kinematic modeling function  * describes the adjoining link coordinate frames with  parameters, none of which is linearly dependent, so that there are  = × kinematic parameters to be identified for a  DOF robot.The end-effector pose is expressed by  parameters, and thus the minimum amount of measurement configurations can be calculated as follows: In practical application, for more robustness of robot calibration against the measurement noise, several times, or even ten times, the minimum number required of measurement configurations is selected.On the other hand, the composition of the measurement configuration set is also worth studying.Mapping of the kinematic parameter error and the estimate pose error of the end-effector, here called the identification matrix J Ω (Ω, Φ), obviously depends on the configuration set Φ. A better configuration set can reflect more effects of the kinematic parameter error on the pose error for higher identification accuracy, which also lays the foundation for further optimization of the measurement configurations.

Observability Problems
As stated above, the process of robot calibration involves a number of error factors, such as modeling errors, highorder residuals from linearization, measurement noises, and other random errors.Under the assumption that linearization errors can be remedied through iteration and modeling errors have a slight impact, (3) is mainly influenced by measurement noises   , and in the lecture of regression it can be written as where y = Ψ is a  × 1 observation vector,  = Ω an  × 1 estimate variable vector, X = J Ω (Ω, Θ) a  ×  regression Mathematical Problems in Engineering matrix, and   a  × 1 error vector, given  =  ⋅ .According to the definition of the identification matrix, where Φ = [Θ 1 , Θ 2 , . . ., Θ  ] is a set composed of  robot configurations, and Ω = [Ω 1 , Ω 2 , . . ., Ω  ] is a set composed of  kinematic parameters to be identified.For the generalized regression model ( 9), the information matrix is expressed as In statistics, optimal experiment design theory has given rise to several evaluation indices of the goodness of the information matrix in order to direct experiment design, of which the most significant are defined as follows [5]: (i) D-optimality: selection of experiment design maximizes the determinant of the information matrix (ii) A-optimality: selection of experiment design minimizes the trace of the information matrix (iii) E-optimality: selection of experiment design maximizes the minimum eigenvalue of the information matrix (iv) G-optimality: selection of experiment design minimizes the maximum variance of the response variables The relation between optimal experiment design indices and observability indexes are elaborated in this section for building an optimization model of configuration set, according to the two typical purposes of robot calibration, one of which is to minimize the variance of the estimate parameters in the robot model and the other one is to minimize the maximal possible uncertainty of the estimate end-effector pose.These two purposes are not equivalent in optimal experiment design theory.The conclusion comes out that maximizing  1 always means an unbiased estimate of the kinematic parameters with smallest variance, while  3 is the right criterion for minimizing the possible error of the endeffector for any configurations.

Optimization of Kinematic Parameter Estimation.
With least squares, the kinematic parameter error  can be estimated as where (X  X) is assumed to be invertible.Minimization the variance var( β) of the estimated variable  indicates an unbiased estimate of the kinematic parameters.Assume that each factor of the measurement error vector   is independent and identically distributed, such as where   denotes the th factor in the error vector   and  denotes the number of parameters expressing the endeffector pose.Hence the covariance matrix of the measurement error vector can be expressed as where var(•) is the variance of a vector and  denotes the amount of measurement configurations.
According to (14), the covariance matrix of the estimated variable  can be obtained as Given the definition of the information matrix as ( 11), ( 15) can be rewritten as Equation (16) indicates that the smaller the information matrix M is, the more accurate the estimate of kinematic parameters will be obtained.However, the minimization of a matrix has lots of expressions corresponding to different physical meanings.D-optimality seeks an optimal estimate of the estimated variable , and for a linearization model such as (9), the volume of the dense ellipsoid defined by the estimated variable  is adopted to evaluate the efficacy of least squares estimation.
The variance and mean of the estimated variable  are assumed to be known, and an ellipsoid called the dense ellipsoid exists in m-dimensional space, on the region surrounded by which, m-dimensional random variables have the same variance and mean as an m-dimensional uniform distribution.The volume of the dense ellipsoid indicates the concentrative degree of the estimated variable.Generally different experiment designs produce different dense ellipsoids, and the measurement configuration set plays a decisive role in robot calibration, such as  =  V (Φ).For a configuration set Φ, or called an experiment design, the volume of the dense ellipsoid can be calculated with its information matrix as where Γ(•) is a Γ function.For the same model ( 9), assume that there are two different designs such as Φ 1 and Φ 2 with (Φ 1 ) < (Φ 2 ).Under D-optimality, Φ 1 is referred as a better design than Φ 2 , indicating that det (Φ 1 ) > det (Φ 2 ).The determinant of the information matrix is finally selected to represent its magnitude.
Using the singular value decomposition method, (9) becomes where U and V are orthonormal matrices, and for  measurements and  parameters, where   ,  = 1, 2, . . .,  are the nonzero singular values of the matrix  with  1 ≥  2 ≥ ⋅ ⋅ ⋅ ≥   ≥ 0. The observability index  1 is defined as Considering that reciprocal of product of eigenvalues is actually the determinant of a matrix, the determinant can be expressed as It can be concluded from ( 21) that the determinant of the information matrix M is typically equivalent to observability index  1 , which can be used as the evaluation index of configuration set optimization, guaranteeing an unbiased estimate of the kinematic parameters.

Optimization of End-Effector Pose Estimation.
With the parameters β estimated, the response variable ỹ can be estimated as For a robot, the end-effector pose can be obtained through the kinematic parameters calibrated with the forward kinematic function.Since the estimate kinematic parameters are still not exact, the estimate end-effector pose has errors, too.In some cases, the goal of robot calibration is to build a prediction model to place the robot end-effector to a certain pose accurately.In statistics, the experiment design under G-optimality indicates that the maximum of variance of each point in the design domain with respect to the GLSE (generalized least squares estimation) estimate β is minimized; that is to say the maximum of estimate variance of the responsible variable ỹ is minimized.
For (22), minimizing the maximum possible variance of the responsible variable ỹ means the minimization of max{var(ỹ)}.For any row x of the identification matrix X, ‖x‖ ≤ 1; thus In practice, all experimental designs are exact; it is supposed to be an -point design with one observation at x.For an exact design, the information matrix is expressed as For an -point exact design, the G-optimality is According to the paper [5], if the input variable is bounded as ‖x‖ ≤ 1, the exact design G-optimality is equivalent to E-optimality, which can be written as where  min denotes the minimal nonzero singular value of the identification matrix.Since the observability index  3 is defined by the minimal nonzero singular value  min , it should be the right criterion for minimizing the possible error of the end-effector for any configurations, for the same number of measurement configurations.

Modeling of Configuration Set Optimization.
A combination of the indexes  1 and  3 is used to achieve a multicriteria optimization of measurement configurations in an attempt to obtain an unbiased estimate of the kinematic parameters and an optimal estimate of the end-effector pose at the same time.
In the experiment design lecture, selection of measurement configurations in robot calibration can be thought as an point design.The design problem can be expressed as seeking a configuration set whose size is , configurations of which are selected from the candidate configuration set pool .The indexes  1 and  3 corresponding to the configuration set designed reach their maximums. max where A is a function to achieve integration of two indicators, and there are different integration strategies in practical applications.In general, the subobjective functions of multiobjective optimization are transferred into a single objective function by mathematical transformation and then singleobjective optimization methods are used to solve the optimization problem.Among these methods, the hierarchical sequence method is an analytical decision-making method optimizing one by one in the order of the importance of the objectives while the objective planning method is to make each objective converge to a given corresponding target value under constraints as close as possible, rather than directly seeking an extremum of each objective, which has a characteristic of high solving efficiency but requires a priori knowledge.
The influence of measurement configurations on the two observability indexes was investigated to determine a Mathematical Problems in Engineering better optimization strategy.Since the complete expression of the relation between measurement configurations and observability indexes is too complex to show its influence intuitively using an analytical method, in Section 5, kinds of robots were adopted to analyze the influence of measurement configurations.As the basis for the design of the calibration configuration optimization model, conclusions were summed up as follows: (1) The indexes  1 and  3 almost have no relevance with each other for the same configuration set.
(2) With the incensement of the size of the measurement configuration set, its corresponding  1 converges to a certain system finally, while  3 always has an increasing trend.
(3) The index  1 has an upper limit according to the second conclusion.
(4) When the measurement configuration set contains more than a certain number of configurations, the influence of composition of the measurement configuration set on  3 is greatly higher than that on  1 .
Taking the above conclusions into account, a combination of the hierarchical sequence method and the objective planning method was used for an optimal experiment design, where the minimal size of the measurement configuration set as its  1 converges to a certain system, recorded as , and the converge value, recorded as  1  , should be firstly determined.And then an -point design is produced to maximize  3 with Θ  ∈ ,  = 1, 2, . . ., . (28)

Modified Experiment Design Method
It is well-known that the DETMAX algorithm and its improved versions have been widely used in experiment design of robot calibration [18][19][20][21].The traditional DETMAX algorithm relies highly on the principle of exchange, including the PoseAdded algorithm and the PoseRemoved algorithm.In calibration applications, the PoseAdded algorithm aims to pick up an optimal configuration one time from the candidate configuration set pool to be added into the original -point design in order to yield the greatest increase of the corresponding observability index; on the contrast, the PoseRemoved algorithm aims to remove one configuration from the set constructed in the previous step which affects the index value minimally.These two algorithms are used interchangeably, and the observability index of the -point experiment design increases, thus achieving the optimal design.
To settle with the continuousness of the permissible operating range of all the robot's actuators, a rasterization method is adopted by the DETMAX algorithm.Nevertheless, the amount of configurations in the candidate configuration set pool will be greatly huge if a small interval is used during rasterization, especially when the robot has plenty of DOFs.
Another flaw in this algorithm is that the original selection of the -point design has a big effect upon the final optimal design, so that the DETMAX algorithm is executed several times to choose the best solution, implying that the algorithm cannot guarantee a global optimal solution.Accordingly, the DETMAX algorithm was improved in this paper from two aspects: initial set construction and iterative optimization, for a better design performance.
4.1.Ways to Improve DETMAX Algorithm.In the original algorithm, configurations in the initial set are randomly selected from the candidate configuration set pool.And the complexity of the configuration optimization problem makes the optimization result heavily dependent on the random selection, indicating that a good initial set will lead to better optimization results.Thus the following methods are proposed to build the initial set.
The Incremental Method.The minimum amount of measurement configurations needed for robot calibration can be calculated as (8), which is much less than the actual amount of measurement configurations used to guarantee robustness against measurement noises, giving us a new idea to build the initial set.The minimum amount of configurations is randomly selected from the candidate configuration set pool, and then a configuration, which contributes most to the increment of the observability index, is picked up from the remaining configuration set one by one to be added into the initial set, until the size of the measurement configuration set reaches .
The Decreasing Method.The incremental method is to find the optimal configuration in the remaining configuration set by traversal to increase the amount of the initial configuration set.Firstly, a large amount of configurations is randomly selected from the candidate configuration set pool, and then a configuration, which contributes least to the increment of the observability index, is removed from the initial configuration set one by one, until the size of the measurement configuration set reduces to .
The set-constructing method is proposed with the incremental or the decreasing method, which makes the resulting measurement configuration set a better observability index on one hand and reduces the time required for iterative optimization on the other hand, since a better initial set is usually closer to the local optimal solution.However the result of a single construction may still be unsatisfactory and thus several tries are quite necessary.
The DETMAX algorithm replaces the inappropriate configuration in the initial set with a better configuration from the remaining set one by one until the construction of an optimal experiment design.And the rasterization makes it possible to pick up the best configuration by traversing the candidate configuration set pool, but difficult to pursue better optimization results on the continuous permissible operating range.Consequently a set-shrinking method aiming to refine the candidate configuration set pool is proposed as follows.
The Shrinking Method.By constructing the initial set and the subsequent exchange optimization, an -point design Φ  is obtained, which is composed of  measurement configurations, and each configuration is expressed by  joint angles, for a  DOF robot.The -point design Φ  can be written as where  min  and  max  denote the lower and upper limit of the th joint angle,  denotes the number of the configurations Θ in the configuration set Φ  ,  is the number of times the corresponding joint angle is iteratively optimized, and Δ is the stated width of the new permissible range.Each joint angle is supposed to be iteratively optimized no more than  times.
The new permissible ranges are rasterized equally in ten pieces to yield a discrete configuration set pool needed.For all the joint angles in the set Φ  , the best configuration will be picked up from the new configuration set pool for replacement, constructing a better -design Φ  .During the next iteration process, a new configuration set pool corresponding to the new set Φ  can be obtained as (30).The optimization procedure expressed above will be repeated until the times that all the joint angles are optimized reach the upper limit .

Application of Modified DETMAX Algorithm.
The modified DETMAX algorithm is applied to solve the optimization problem illustrated in (28) with the proposed setconstructing method and set-shrinking method applied.Compared with the foregoing measurement configurations optimization methods, the proposed one has the following characteristics: (i) The optimization of observability indexes has been taken into account while building the initial configuration set just from the beginning.
(ii) Two observability indexes are selected for evaluation of the measurement configuration set, with  1 to be the objective and  3 to be the constraint.
(iii) The set-shrinking method makes it possible to further refine the measurement configurations.
Based on the conventions given in Section 2, more is added to outline the realization of the new algorithm in detail as follows: (i) Φ  is a set of  measurement configurations, also called an -point experiment design.
(ii) Φ +1 is a set of ( + 1) measurement configurations, which is a new measurement design after the PoseAdded algorithm is used.
(iii) Φ − is the set of all the possible configurations except these in Φ  , where the added configuration is picked up.(iv) Θ  denotes a robot configuration chosen by Pose-Added algorithm.(v) Θ  denotes a configuration that should be removed from the configuration set Φ +1 by the PoseRemoved algorithm.(vi)   denotes the minimum amount of configurations needed for calibration, calculated by ( 8).(vii)   denotes the maximum amount of measurement configurations given for analyzing the influence of configurations on observability indexes.
(viii) Φ ℎ is a configuration set whose size is ℎ.
(ix)  1  denotes the converging value of the observability index  1 as the size of the measurement configuration set increases.(x)   denotes the th joint angle of the th configuration in the configuration set Φ  .
(xi)    is the new feasible domain calculated as (30).(xii)   is the number of times the joint angle   has been replaced, whose upper limit is set as .
As shown in Figure 1, the improved DETMAX algorithm can be roughly divided into four steps: determination of the optimal size of the measurement configuration set; construction of the initial set by the incremental or decreasing method; iterative optimization by the PoseAdded and PoseRemoved algorithm; and further optimization by the shrinking method.Specific steps are as follows.
(1) Calculate the observability index  1 corresponding to different amounts of measurement configurations and analyze the trend of  1 changing with respect to the amount.The optimal size of the measurement configuration set  and the lower limit  1  can be determined as follows.
Step 1. Assume the initial size of the configuration set to be ℎ =   .
Step 2. ℎ configurations are selected randomly from the candidate configuration set pool to yield Φ ℎ .
Step 3. Calculate the observability index  1 corresponding to the configuration set Φ ℎ .

Initialize the candidate configuration set pool by rasterization
Determine the optimal size of the measurement configuration set and the lower limit of O1 Step 5. Draw the curve indicating the trend of the observability index  1 with the incremental size of the measurement configuration set and choose the minimal number as  1 converges to a threshold to be the optimal size, recorded as .
Step 6.The lower limit  1  can be obtained as the mean of the observability index  1 corresponding to all the configuration sets Φ ℎ with  ≤ ℎ ≤   .
(2) Two initial set constructing methods including the incremental method and the decreasing one are proposed, but due to limited space, only the application of the incremental method is presented afterwards.It is worthy of note that since either of the methods get rid of the effects of random selection, it is still necessary to construct the initial set several times for choosing the best one.
Step 1. Assume the initial size of the configuration set to be ℎ =   .
Step 2. ℎ configurations are selected randomly from the candidate configuration set pool to yield Φ ℎ .
Step 3. one configuration is selected from the remaining configuration set Φ −ℎ , which yields the highest observability index  1 corresponding to the new measurement configuration set.
Step 4. Let ℎ = ℎ + 1.If ℎ ≤   then go to Step 2; otherwise construction of the initial configuration set Φ  is done.
(3) Since the two observability indexes  1 and  3 are taken into account, the actual application of the PoseAdded algorithm and the PoseRemoved algorithm are slightly different from the previous way, as indicated below.
Step 1.One configuration called Θ  is obtained by the PoseAdded algorithm from Φ − to combine with Φ  and create a new set Φ +1 , which maximizes the observability index  3 of the configuration set Φ +1 with  1 ≥   1 .
Step 2. The PoseRemoved method is then adopted on the configuration set Φ +1 to remove one configuration from it which affects the index value  3 of Φ +1 minimally, still with  1 ≥   1 .
Step 3. If the removed configuration Θ  is different from the added one Θ  during this iteration, go to Step 1; otherwise this step stops.
(4) Constructing an initial configuration set as well as the following iterative optimization has the need to traverse the entire configuration set pool, and thus a large interval is used in rasterization for the sake of time-saving and economy, which makes it hard to obtain a better experiment design of robot calibration.Consequently, the shrinking method helps to remedy this problem.
Step 1.Let   = 0 for all the joint angles   .
Step 2. The new feasible domain    is obtained by (30) for all the joint angles   .
Step 3. Find a joint angle   with   <  in the configuration set Φ  , which can be replaced by one selected from its own feasible domain    and yield a new configuration set with the highest observability index  3 with  1 ≥   1 , and then let   =   + 1.
Step 4. If   for each joint angle in the configuration set Φ  satisfies   ≥ , this step stops and the optimal experiment design for robot calibration is obtained.Otherwise, go on to Step 2.

Simulation Study
There are two purposes of simulation study here: the first is to investigate the influence of measurement configurations on the indexes  1 and  3 ; the second is to verify the validity of the proposed experiment design method.In simulation applications, several kinds of robots with different structures are taken into consideration for the sake of universality, while control groups are given for validation such as the design method with only  1 to be optimized and the design method by random selection.According to the analysis in Section 2, the identification precision of the kinematic parameters and the estimate accuracy of the end-effector pose are chosen to serve as performance evaluation indicators of robot calibration.

Kinematics Modeling and Simulation
Objects.The mapping relating the end-effector pose Ψ * of the robot to its kinematic parameters Ω * and configurations Θ * is illustrated as (1), where the modeling function  * decides not only the contents of the kinematic parameters but also the concrete expression of the function  * .
Various methods of kinematics modeling for robots have been proposed, including DH method [16], S-model method [22], CPC (Complete and Parametrically Continuous) model [23], and zero reference model method [24].And a modified DH method, called MDH method or Hayati method [25], is used for kinematics modeling, which describes the deviation between two adjacent parallel joints axes with an additional rotation transformation to remedy the incompleteness.The coordinate frame is established by MDH method as shown in Figure 2. The MDH method uses five parameters (let  = 5) including   ,   ,   ,   ,   to describe the two adjoining link coordinate frames, and the homogeneous transformation between them is shown as follows: where rot and trans denote the transpose and rotation matrix, respectively, and −1 A  denotes the homogeneous transformation of joint  with respect to joint  − 1.For a  DOF robot without a tool coordinate frame, the end-effector pose can be expressed by the homogeneous transformation as follows: where 0 A  is the homogeneous transformation of the endeffector with respect to the base coordinate frame, which can be transferred into -- Euler angles.So, the pose of the end-effector pose can be described by three parameters of position and three of angles, all of which are independent, such as  = 6 in (8).And then the minimal amount of measurement configurations needed for robot calibration will be obtained with the DOF determined.Finally, the identification matrix can be calculated as (10) so as to obtain the corresponding observability indexes.
As is clear from Figure 3, three kinds of typical robots with different DOF are adopted to analyze the effects of measurement configurations on observability indexes, whose Table 1: Kinematic parameters of the 3-DOF robot.
9 0 0 0 0 kinematic parameters are shown in Tables 1, 2, and 3, respectively.It should be noted that some of the kinematic parameters are redundant, marked as "-" in these tables, which will be ignored during the calibration.The determining methods of redundant parameters can be found in the paper [26].

Influence of Measurement Configurations.
For these robots in simulation, their permissible operating range of each actuator is limited with (−180 ∘ , 180 ∘ ] and then latticed in an interval of 30 ∘ on each joint angle to generate configurations for the candidate configuration set pool .
For the purpose of analyzing the influence of the size of the measurement configuration set, related parameters are determined firstly as illustrated in Table 4, where  is the number of DOF,  denotes the number of the independent parameters to be identified,   denotes the minimal size of the measurement configuration set, and   is set as the maximum size for convenience.
Measurement configuration sets of different sizes from   to   are constructed by random selection from the candidate configuration set pool , and their corresponding observability indexes  1 and  3 are then calculated.
Results are shown in Figures 4, 5, and 6, where the optimal size of the configuration set and the corresponding limit are marked by red dotted lines.It can be illustrated that, for all three robots, with the increase of the size of the configuration set the index  1 increases firstly and then converges to a threshold, while the index  3 always has an increasing trend.
Taking the 4-DOF robot as an example, the difference between  1 and  3 is studied under the same configuration set, whose size is assumed to be 5, 10, and 20, respectively.Then according to the given size, 50 configuration sets are constructed with random selection and their corresponding indexes  1 and  3 are calculated.
In order to observe the relationship between the two indexes, 50 configuration sets are sorted in order of the value of their index  1 .Results are shown in Figures 7,  8, and 9. Aiming to quantify the varying magnitude of the observability indexes in the three figures, the following expression is presented: where  max denotes the maximum value of the observability index appearing in the 50 samples,  min denotes the minimum one, and then   is an evaluation index of the magnitude of variation in the observability index.
As shown in Figures 7, 8, and 9, the varying magnitude of the indexes  1 and  3 is 27.9% and 224.3%, 9.8% and 137.2%, and 7.7% and 73.6%, respectively.It can be concluded that the varying magnitude of the index  1 is always much smaller than that of  3 , though both of them decrease as the size of the measurement configuration set increases.Moreover, the observability indexes  1 and  3 almost have no relevance with each other for the same measurement configuration set.Composition of the measurement configuration set has a greater influence on the index  3 than on the index  1 .
The conclusions drawn from the simulation study above are verified by more experiments and there seems to be similar laws on different robots to provide the basis for optimization of measurement configuration.The conclusions of the first three graphs make it possible to find out the optimal size of measurement configuration set with  1 maximal economically while the last three indicate the maximization of  3 with  1 to be acceptable so as to produce a win-win situation.

Simulation on Robot Calibration.
It is necessary to carry out simulation study on robot calibration for validation.In practical applications, the actual end-effector pose of the robot can be obtained by external measurement instruments such as the optical measuring instrument and the CMM, while the actual joint angles can be obtained by the encoder.Then errors of the kinematic parameters are calculated as illustrated in (4), whereas in situational applications, the pose measuring noise and encoder noise are added according to their respective distributions in order to imitate the influence of sensor errors.The process of calibration simulation is shown in Figure 10, where firstly the measurement configuration set is constructed by a certain design method within the permissible operating range of joints and the given amount of configurations, and then the actual end-effector poses are calculated with the actual kinematic parameters.Finally, taking the encoder noise and the measurement noise into account, least squares method is adopted to identify the kinematic parameters.It is worthy of note that the actual parameters are calculated by the nominal parameters artificially added with parameter errors, indicating that all these parameters are known for analysis and comparison.
Taking the 4-DOF robot as an example, the angle errors are set as 0.1 ∘ while the length errors are set as 0.01 m, and the actual parameters can then be obtained.As shown in Figure 5 the optimal size of the measurement configuration set is  = 10, the lower limit is  1  = 0.934, and the minimum and maximum size are 3 and 50, respectively, as presented in Table 4.Moreover, the measurement noises are defined as normal distribution random variables with (0, 10 −3 ), whose unit is mm for the length values and rad for the angle values, while the encoder noises are defined with (0, 10 −8 ), whose unit is rad for the angle values.Finally 100 configuration sets are constructed by selecting configurations randomly from the candidate configuration set pool to serve as the validation group.Simulations are carried out under three circumstances.
(1) The calibration was carried out under the measurement configuration set whose size varies from 3 to 50, with the same design method of random selection.
(2) Let the size of the measurement configuration set be 10, and one measurement configuration set was constructed, respectively, under three methods, including the method proposed in this paper, the method maximizing  1 by DETMAX, and random selection.
(3) The calibration was carried out, respectively, under the three measurement configuration set obtained above for 500 times.

Results and Analysis.
As mentioned in this paper, the purpose of robot calibration is to obtain the accurate parameters representing the robot structure and the exact estimation of the end-effector pose.Based on this, the following criteria are proposed to evaluate the calibration performance.
(1) Capability Analysis of Identified Parameters Error.Since the magnitude of the kinematic parameters and their units are inconsistent, it is hard to compare different calibration results.So the level of similarity between the identified parameter errors and the actual ones is used as the criteria for identification accuracy of kinematic parameters, expressed as follows: where Ω denotes the difference between the actual parameters Ω and the nominal parameters Ω * , and Ω  denotes the calibrated parameters Ω  , so  Ω [] is the identification accuracy of the th kinematic parameter.For one calibration, the mean [P Ω ]and var[P Ω ]of the vector  Ω are calculated to present the identification performance of all the kinematic parameters.For one measurement configuration set designed, the distribution of the mean and variance obtained above is investigated to evaluate the performance of calibration.
(2) Capability Analysis of Euclidian Error Distance.The estimate accuracy of the end-effector pose under a certain configuration is evaluated with the end-effector position error as the standard: where P  denotes the estimate end-effector position of the robot, P  denotes the actual position, and  ,, are the deviations in , , and  direction with respect to the work reference frame.The end-effector position errors of 100 validation configurations and their mean and variance are calculated, respectively, in order to evaluate the end-effector position estimate capability.The influence of the size of the measurement configuration set on the calibration performance is shown in Figures 11 and 12.It can be concluded from the figures that, with the increase of the size, either the kinematic parameter identification index or the end-effector position estimate index has an improving trend and then converges to a certain system, indicating the optimization of observability indexes help to improve the calibration performance, but excessive configurations are not necessary.
As is clear from Figures 11 and 12, sufficient measurement configurations can yield an acceptable calibration result but consume unnecessary resources and time.Therefore, an optimal amount of configurations should be worked out that can maintain both economy and observability.Based on the optimal number shown in Table 4, three measurement configuration sets are constructed with different methods as shown in Tables 5, 6, and 7.The observability indexes  1 and  3 corresponding to the three measurement configuration sets are 0.9389 and 1.0827, 0.8442 and 0.5925, and 0.7270 and 0.4051, respectively.Obviously the measurement configuration set constructed with the proposed method has the highest observability indexes among the three.Calibration under these three configuration sets was repeated 500 times, and typical results are shown in Figures 13 and 14, where the calibration performances are presented by the position estimate errors and the identified parameters errors, respectively.Just in terms to the position estimate errors, Figure 13(a) indicates that the mean errors of the three configuration sets are 0.70 mm, 1.0 mm, and 2.1 mm, respectively, while the mean errors are 1.3 mm, 1.4 mm, and 1.1 mm, respectively.It can be concluded that the configuration set with the method proposed in this paper does better in calibration than the other two sets, while the three configuration sets have the similar calibration performance in Figure 14.Actually, the result of a single calibration simulation highly relies on the noises produced randomly and may not embody the real calibration performance of the configuration set, indicating the necessity of repeated calibration simulations.For a single calibration simulation, the mean and variance of the estimate position errors under the 100 validating configuration sets were calculated and so were these of the errors of the identified kinematic parameters.Results from the 500 simulations are shown in Figures 15 and 16, and it can be concluded that the method maximizing  1 by DETMAX does better than random selection in constructing an optimal measurement configuration set, while the method proposed in this paper does the best, either in identification accuracy of the kinematic parameters or in the estimate errors of the end-effector position.

Conclusion
A novel optimal design method of measurement configurations is proposed in this paper for robot calibration, where the optimal amount of measurement configurations is firstly determined by analysis of the influence of measurement configurations on the observability indexes, and then according to the two purposes of robot calibration which are to minimize the variance of the estimate parameters in the robot model and to minimize the maximum possible uncertainty of the estimate end-effector pose, the optimizing model with  3 to be the objective and  1 to be the constrain is built.Finally, the existing DETMAX algorithm is modified to seek optimal design of measurement configurations, by adopting a set-constructing method and a set-shrinking method.In order to verify the effectiveness of the proposed method, simulation study on robots of different structures was carried out.It is demonstrated by the simulation results that the proposed method embodies an optimized algorithm to reach higher accuracies with fewer measurement configurations compared with previous achievements.

Figure 2 :
Figure 2: The coordinate frame established by MDH method.

Figure 5 :
Figure 5: The influence of the size of the configuration set on the indexes  1 and  3 for the 4-DOF robot.

Figure 6 :Figure 7 :
Figure 6: The influence of the size of the configuration set on the indexes  1 and  3 for the 5-DOF robot.

Figure 8 :Figure 9 :
Figure 8: The difference between  1 and  3 on the same configuration set with  = 10.

Figure 11 :Figure 12 :
Figure 11: The estimate errors of the end-effector position with respect to different sizes of configuration sets.

Figure 13 :Figure 14 :
Figure 13: The calibration performances in the position estimate errors (a) and the identified parameters errors (b).

Figure 16 :
Figure 16: Distribution of the means (a) and variances (b) of the identified parameters errors in 500 simulations.

Table 5 :
The configuration set constructed with the method proposed in this paper.