# A Novel Optimal Design of Measurement Configurations in Robot Calibration

Authors: Qingxuan Jia; Shiwei Wang; Gang Chen; Lei Wang; Hanxu Sun
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Robot calibration highly depends on redundant measurement 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.


---

## Body

## 1. 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. Borm and Menq [4] proposed an index known as O1, which is based on the product of the nonzero singular values of the identification matrix. Maximization of O1 indicates more effect of kinematic parameter errors exerted on the estimate pose error. This is also true for the index O5, 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 O2, 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. O2 obtains its maximum value at 1 when all singular values are equal. The index O3 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 O4 [6] is the product of O2 and O3, and actually there is no significant difference from O2 to O4. Though both O1 and O5 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 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.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 ofO1 and O3 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 O1 always means an unbiased estimate of the kinematic parameters with smallest variance, while O3 is the right criterion for minimizing the possible error of the end-effector for any configuration, from the experimental design literature. DETMAX 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.


## 2. Essential Constraints of Calibration

Without loss of generality, a serial robot withd actuators controlling d 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 θi denotes the angle of the ith actuator, where i=1,2,…,d.(ii)
Θ is a d×1 vector composed of d 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 f∗ based on the irreducible exact kinematic parameters as illustrated in the following.(1)f∗Ψ∗,Ω∗,Θ∗=0The 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 εf. In conclusion, it is the modeling error εf, 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.(2)δΨ=JΩΩ,ΘδΩ+JΘΩ,ΘδΘ+εf,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:(3)δΨ=JΩΩ,ΘδΩ.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,…,δΨnT,JΩ=[δJΩ1,δJΩ2,…,δJΩn]T. The configuration Θ is replaced with a configuration set Φ. Here (•)T indicates a transpose of a matrix, and n 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.(4)δΩi+1=JΩΩ+∑j=1iδΩj,Φ+δΨ,where i denotes the total number of iterations, and (•)+ means the Moore-Penrose inverse of a matrix. The estimate pose error is obtained by δΨ=Ψm-Ψ∗, in which the exact values Ψ∗ are replaced with the measured values Ψm. The process is iterated until δΩi converges to a small value. The identified kinematic parameters can finally be obtained by (5). Note that the measured values Ψm is different from the exact values Ψ∗ for defects of measuring equipment and methods, and the error is set as εm=Ψ∗-Ψm.(5)Ωc=Ω+∑j=1iδΩj.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 end-effector 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:(6)JΘΩ,ΘΔΘ≤εr,∀Θ∈ξ,where ΔΘ is the maximum value of the joint angle measurement error and εr 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(7)detJΩT•JΩ>0,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 Ψm 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-6m.Constraint D.On the one hand, sufficient measurement configurations are selected to establish an overdetermined system of equation. Assume that the kinematic modeling function f∗ describes the adjoining link coordinate frames with p parameters, none of which is linearly dependent, so that there are m=p×d kinematic parameters to be identified for a d DOF robot. The end-effector pose is expressed by q parameters, and thus the minimum amount of measurement configurations can be calculated as follows:(8)n≥mq.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 matrixJΩΩ,Φ, 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.


## 3. Observability Problems

As stated above, the process of robot calibration involves a number of error factors, such as modeling errors, high-order 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 εm, and in the lecture of regression it can be written as(9)y=Xβ+εm,where y=δΨ is a k×1 observation vector, β=δΩ an m×1 estimate variable vector, X=JΩΩ,Θ a k×m regression matrix, and εm a k×1 error vector, given k=n·q. According to the definition of the identification matrix,(10)XΩ,Φ=∂f∗Ω,Θ1∂Ω1⋯∂f∗Ω,Θ1∂Ωm⋮⋱⋮∂f∗Ω,Θp∂Ω1⋯∂f∗Ω,Θp∂Ωm,where Φ=Θ1,Θ2,…,Θn is a set composed of n robot configurations, and Ω=Ω1,Ω2,…,Ωm is a set composed of m kinematic parameters to be identified. For the generalized regression model (9), the information matrix is expressed as(11)M=XTX.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 variablesThe 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 maximizingO1 always means an unbiased estimate of the kinematic parameters with smallest variance, while O3 is the right criterion for minimizing the possible error of the end-effector for any configurations.


### 3.1. Optimization of Kinematic Parameter Estimation

With least squares, the kinematic parameter errorβ can be estimated as(12)β~=X+y=XTX-1XTy,where (XTX) 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 εm is independent and identically distributed, such as(13)εi~N0,σ2∀εi∈εm,i=1,2,…,q,where εi denotes the ith factor in the error vector εm and q denotes the number of parameters expressing the end-effector pose. Hence the covariance matrix of the measurement error vector can be expressed as(14)varεm=σ2Iq×n,where var(•) is the variance of a vector and n denotes the amount of measurement configurations.According to (14), the covariance matrix of the estimated variable β can be obtained as(15)varβ~=σXTX-1.Given the definition of the information matrix as (11), (15) can be rewritten as(16)varβ~=σM-1.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=fvΦ. For a configuration set Φ, or called an experiment design, the volume of the dense ellipsoid can be calculated with its information matrix as(17)VΦ=m+2m/2πm/2Γm/2+1det⁡MΦ,where Γ• is a Γ function. For the same model (9), assume that there are two different designs such as Φ1 and Φ2 with VΦ1<VΦ2. Under D-optimality, Φ1 is referred as a better design than Φ2, indicating that det⁡MΦ1>det⁡MΦ2. The determinant of the information matrix is finally selected to represent its magnitude.Using the singular value decomposition method, (9) becomes(18)y=UΣVTβ+εm,where U and V are orthonormal matrices, and for n measurements and m parameters,(19)Σ=σ10⋯0⋮⋮⋮⋮00⋯σm00⋯0⋮⋮⋮⋮00⋯0,where σi,i=1,2,…,m are the nonzero singular values of the matrix X with σ1≥σ2≥⋯≥σm≥0. The observability index O1 is defined as(20)O1=σ1σ2⋯σm1/mn.Considering that reciprocal of product of eigenvalues is actually the determinant of a matrix, the determinant can be expressed as(21)det⁡M=∏i=1m1σi2=1O12mnm.It can be concluded from (21) that the determinant of the information matrix M is typically equivalent to observability index O1, which can be used as the evaluation index of configuration set optimization, guaranteeing an unbiased estimate of the kinematic parameters.


### 3.2. Optimization of End-Effector Pose Estimation

With the parametersβ~ estimated, the response variable y~ can be estimated as(22)y~=Xβ~.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 y~ is minimized.For (22), minimizing the maximum possible variance of the responsible variable y~ means the minimization of max⁡vary~. For any row x of the identification matrix X, x≤1; thus(23)minmaxvary~∣x∣x∣X=minmaxxTvary~∣xx∣x∣X=minmaxσ2xTXTX-1x∣x∣X.In practice, all experimental designs are exact; it is supposed to be ann-point design with one observation at x. For an exact design, the information matrix is expressed as(24)M=1nXTX.For ann-point exact design, the G-optimality is(25)min⁡max⁡nxTXTX-1x∣x∣X.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(26)min⁡max⁡nxTXTX-1x∣x∣X∝max⁡σmin2,where σmin denotes the minimal nonzero singular value of the identification matrix. Since the observability index O3 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.


### 3.3. Modeling of Configuration Set Optimization

A combination of the indexesO1 and O3 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 n-point design. The design problem can be expressed as seeking a configuration set whose size is n, configurations of which are selected from the candidate configuration set pool ξ. The indexes O1 and O3 corresponding to the configuration set designed reach their maximums.(27)maxOIΦ=AO1Φ,O3ΦsubjecttoΦ=Θ1,Θ2,…,Θn,Θi∈ξ,i=1,2,…,n,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 single-objective 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 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 Section5, 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 O1 and O3 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 O1 converges to a certain system finally, while O3 always has an increasing trend.( 3 ) The index O1 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 O3 is greatly higher than that on O1.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 itsO1 converges to a certain system, recorded as n, and the converge value, recorded as O1s, should be firstly determined. And then an n-point design is produced to maximize O3 with O1≥O1s.(28)maxO3ΦsubjecttoO1Φ≥O1sΦ=Θ1,Θ2,…,Θn,Θi∈ξ,i=1,2,…,n.


## 3.1. Optimization of Kinematic Parameter Estimation

With least squares, the kinematic parameter errorβ can be estimated as(12)β~=X+y=XTX-1XTy,where (XTX) 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 εm is independent and identically distributed, such as(13)εi~N0,σ2∀εi∈εm,i=1,2,…,q,where εi denotes the ith factor in the error vector εm and q denotes the number of parameters expressing the end-effector pose. Hence the covariance matrix of the measurement error vector can be expressed as(14)varεm=σ2Iq×n,where var(•) is the variance of a vector and n denotes the amount of measurement configurations.According to (14), the covariance matrix of the estimated variable β can be obtained as(15)varβ~=σXTX-1.Given the definition of the information matrix as (11), (15) can be rewritten as(16)varβ~=σM-1.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=fvΦ. For a configuration set Φ, or called an experiment design, the volume of the dense ellipsoid can be calculated with its information matrix as(17)VΦ=m+2m/2πm/2Γm/2+1det⁡MΦ,where Γ• is a Γ function. For the same model (9), assume that there are two different designs such as Φ1 and Φ2 with VΦ1<VΦ2. Under D-optimality, Φ1 is referred as a better design than Φ2, indicating that det⁡MΦ1>det⁡MΦ2. The determinant of the information matrix is finally selected to represent its magnitude.Using the singular value decomposition method, (9) becomes(18)y=UΣVTβ+εm,where U and V are orthonormal matrices, and for n measurements and m parameters,(19)Σ=σ10⋯0⋮⋮⋮⋮00⋯σm00⋯0⋮⋮⋮⋮00⋯0,where σi,i=1,2,…,m are the nonzero singular values of the matrix X with σ1≥σ2≥⋯≥σm≥0. The observability index O1 is defined as(20)O1=σ1σ2⋯σm1/mn.Considering that reciprocal of product of eigenvalues is actually the determinant of a matrix, the determinant can be expressed as(21)det⁡M=∏i=1m1σi2=1O12mnm.It can be concluded from (21) that the determinant of the information matrix M is typically equivalent to observability index O1, which can be used as the evaluation index of configuration set optimization, guaranteeing an unbiased estimate of the kinematic parameters.


## 3.2. Optimization of End-Effector Pose Estimation

With the parametersβ~ estimated, the response variable y~ can be estimated as(22)y~=Xβ~.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 y~ is minimized.For (22), minimizing the maximum possible variance of the responsible variable y~ means the minimization of max⁡vary~. For any row x of the identification matrix X, x≤1; thus(23)minmaxvary~∣x∣x∣X=minmaxxTvary~∣xx∣x∣X=minmaxσ2xTXTX-1x∣x∣X.In practice, all experimental designs are exact; it is supposed to be ann-point design with one observation at x. For an exact design, the information matrix is expressed as(24)M=1nXTX.For ann-point exact design, the G-optimality is(25)min⁡max⁡nxTXTX-1x∣x∣X.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(26)min⁡max⁡nxTXTX-1x∣x∣X∝max⁡σmin2,where σmin denotes the minimal nonzero singular value of the identification matrix. Since the observability index O3 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.


## 3.3. Modeling of Configuration Set Optimization

A combination of the indexesO1 and O3 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 n-point design. The design problem can be expressed as seeking a configuration set whose size is n, configurations of which are selected from the candidate configuration set pool ξ. The indexes O1 and O3 corresponding to the configuration set designed reach their maximums.(27)maxOIΦ=AO1Φ,O3ΦsubjecttoΦ=Θ1,Θ2,…,Θn,Θi∈ξ,i=1,2,…,n,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 single-objective 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 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 Section5, 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 O1 and O3 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 O1 converges to a certain system finally, while O3 always has an increasing trend.( 3 ) The index O1 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 O3 is greatly higher than that on O1.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 itsO1 converges to a certain system, recorded as n, and the converge value, recorded as O1s, should be firstly determined. And then an n-point design is produced to maximize O3 with O1≥O1s.(28)maxO3ΦsubjecttoO1Φ≥O1sΦ=Θ1,Θ2,…,Θn,Θi∈ξ,i=1,2,…,n.


## 4. 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–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 n-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 n-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 then-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 n.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 n.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 n-point design Φn is obtained, which is composed of n measurement configurations, and each configuration is expressed by d joint angles, for a d DOF robot. The n-point design Φn can be written as(29)Φn=Θ1,Θ2,…,Θn,Θi=θi1,θi2,…,θid,i=1,2,…,n.It is clear that ann-point design Φn is composed of n×d joint angles. The new permissible range of each joint angle can be given as(30)θijmin≤θij≤θijmaxθijmin=max⁡θjmin,θij-Δθ2kθijmax=min⁡θjmax,θij+Δθ2ki=1,2,…,n,j=1,2,…,d,k=1,2,…,K,where θjmin and θjmax denote the lower and upper limit of the jth joint angle, i denotes the number of the configurations Θ in the configuration set Φn, k 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 K 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Φn, the best configuration will be picked up from the new configuration set pool for replacement, constructing a better n-design Φn. During the next iteration process, a new configuration set pool corresponding to the new set Φn 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 K.


### 4.2. Application of Modified DETMAX Algorithm

The modified DETMAX algorithm is applied to solve the optimization problem illustrated in (28) with the proposed set-constructing 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, withO1 to be the objective and O3 to be the constraint.(iii)
The set-shrinking method makes it possible to further refine the measurement configurations.Based on the conventions given in Section2, more is added to outline the realization of the new algorithm in detail as follows:(i)
Φn is a set of n measurement configurations, also called an n-point experiment design.(ii)
Φn+1 is a set of n+1 measurement configurations, which is a new measurement design after the PoseAdded algorithm is used.(iii)
Φξ-n is the set of all the possible configurations except these in Φn, where the added configuration is picked up.(iv)
ΘA denotes a robot configuration chosen by PoseAdded algorithm.(v)
ΘR denotes a configuration that should be removed from the configuration set Φn+1 by the PoseRemoved algorithm.(vi)
nL denotes the minimum amount of configurations needed for calibration, calculated by (8).(vii)
nM denotes the maximum amount of measurement configurations given for analyzing the influence of configurations on observability indexes.(viii)
Φh is a configuration set whose size is h.(ix)
O1s denotes the converging value of the observability index O1 as the size of the measurement configuration set increases.(x)
θij denotes the jth joint angle of the ith configuration in the configuration set Φn.(xi)
θijξ is the new feasible domain calculated as (30).(xii)
kij is the number of times the joint angle θij has been replaced, whose upper limit is set as K.As shown in Figure1, 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.Figure 1
Flowchart of the modified DETMAX algorithm.( 1 ) Calculate the observability index O1 corresponding to different amounts of measurement configurations and analyze the trend of O1 changing with respect to the amount. The optimal size of the measurement configuration set n and the lower limit O1s can be determined as follows.Step 1.
Assume the initial size of the configuration set to beh=nL.Step 2.
h configurations are selected randomly from the candidate configuration set pool to yield Φh.Step 3.
Calculate the observability indexO1 corresponding to the configuration set Φh.Step 4.
Leth=h+1. If h≤nM then go to Step 2; otherwise go on to Step 5.Step 5.
Draw the curve indicating the trend of the observability indexO1 with the incremental size of the measurement configuration set and choose the minimal number as O1 converges to a threshold to be the optimal size, recorded as n.Step 6.
The lower limitO1s can be obtained as the mean of the observability index O1 corresponding to all the configuration sets Φh with n≤h≤nM.( 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 beh=nL.Step 2.
h configurations are selected randomly from the candidate configuration set pool to yield Φh.Step 3.
one configuration is selected from the remaining configuration setΦξ-h, which yields the highest observability index O1 corresponding to the new measurement configuration set.Step 4.
Leth=h+1. If h≤nM then go to Step 2; otherwise construction of the initial configuration set Φn is done.( 3 ) Since the two observability indexes O1 and O3 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ΘA is obtained by the PoseAdded algorithm from Φξ-n to combine with Φn and create a new set Φn+1, which maximizes the observability index O3 of the configuration set Φn+1 with O1≥O1s.Step 2.
The PoseRemoved method is then adopted on the configuration setΦn+1 to remove one configuration from it which affects the index value O3 of Φn+1 minimally, still with O1≥O1s.Step 3.
If the removed configurationΘR is different from the added one ΘA 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.
Letkij=0 for all the joint angles θij.Step 2.
The new feasible domainθijξ is obtained by (30) for all the joint angles θij.Step 3.
Find a joint angleθij with kij<K in the configuration set Φn, which can be replaced by one selected from its own feasible domain θijξ and yield a new configuration set with the highest observability index O3 with O1≥O1s, and then let kij=kij+1.Step 4.
Ifkij for each joint angle in the configuration set Φn satisfies kij≥K, this step stops and the optimal experiment design for robot calibration is obtained. Otherwise, go on to Step 2.


## 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 n.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 n.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 n-point design Φn is obtained, which is composed of n measurement configurations, and each configuration is expressed by d joint angles, for a d DOF robot. The n-point design Φn can be written as(29)Φn=Θ1,Θ2,…,Θn,Θi=θi1,θi2,…,θid,i=1,2,…,n.It is clear that ann-point design Φn is composed of n×d joint angles. The new permissible range of each joint angle can be given as(30)θijmin≤θij≤θijmaxθijmin=max⁡θjmin,θij-Δθ2kθijmax=min⁡θjmax,θij+Δθ2ki=1,2,…,n,j=1,2,…,d,k=1,2,…,K,where θjmin and θjmax denote the lower and upper limit of the jth joint angle, i denotes the number of the configurations Θ in the configuration set Φn, k 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 K 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Φn, the best configuration will be picked up from the new configuration set pool for replacement, constructing a better n-design Φn. During the next iteration process, a new configuration set pool corresponding to the new set Φn 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 K.


## 4.2. Application of Modified DETMAX Algorithm

The modified DETMAX algorithm is applied to solve the optimization problem illustrated in (28) with the proposed set-constructing 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, withO1 to be the objective and O3 to be the constraint.(iii)
The set-shrinking method makes it possible to further refine the measurement configurations.Based on the conventions given in Section2, more is added to outline the realization of the new algorithm in detail as follows:(i)
Φn is a set of n measurement configurations, also called an n-point experiment design.(ii)
Φn+1 is a set of n+1 measurement configurations, which is a new measurement design after the PoseAdded algorithm is used.(iii)
Φξ-n is the set of all the possible configurations except these in Φn, where the added configuration is picked up.(iv)
ΘA denotes a robot configuration chosen by PoseAdded algorithm.(v)
ΘR denotes a configuration that should be removed from the configuration set Φn+1 by the PoseRemoved algorithm.(vi)
nL denotes the minimum amount of configurations needed for calibration, calculated by (8).(vii)
nM denotes the maximum amount of measurement configurations given for analyzing the influence of configurations on observability indexes.(viii)
Φh is a configuration set whose size is h.(ix)
O1s denotes the converging value of the observability index O1 as the size of the measurement configuration set increases.(x)
θij denotes the jth joint angle of the ith configuration in the configuration set Φn.(xi)
θijξ is the new feasible domain calculated as (30).(xii)
kij is the number of times the joint angle θij has been replaced, whose upper limit is set as K.As shown in Figure1, 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.Figure 1
Flowchart of the modified DETMAX algorithm.( 1 ) Calculate the observability index O1 corresponding to different amounts of measurement configurations and analyze the trend of O1 changing with respect to the amount. The optimal size of the measurement configuration set n and the lower limit O1s can be determined as follows.Step 1.
Assume the initial size of the configuration set to beh=nL.Step 2.
h configurations are selected randomly from the candidate configuration set pool to yield Φh.Step 3.
Calculate the observability indexO1 corresponding to the configuration set Φh.Step 4.
Leth=h+1. If h≤nM then go to Step 2; otherwise go on to Step 5.Step 5.
Draw the curve indicating the trend of the observability indexO1 with the incremental size of the measurement configuration set and choose the minimal number as O1 converges to a threshold to be the optimal size, recorded as n.Step 6.
The lower limitO1s can be obtained as the mean of the observability index O1 corresponding to all the configuration sets Φh with n≤h≤nM.( 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 beh=nL.Step 2.
h configurations are selected randomly from the candidate configuration set pool to yield Φh.Step 3.
one configuration is selected from the remaining configuration setΦξ-h, which yields the highest observability index O1 corresponding to the new measurement configuration set.Step 4.
Leth=h+1. If h≤nM then go to Step 2; otherwise construction of the initial configuration set Φn is done.( 3 ) Since the two observability indexes O1 and O3 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ΘA is obtained by the PoseAdded algorithm from Φξ-n to combine with Φn and create a new set Φn+1, which maximizes the observability index O3 of the configuration set Φn+1 with O1≥O1s.Step 2.
The PoseRemoved method is then adopted on the configuration setΦn+1 to remove one configuration from it which affects the index value O3 of Φn+1 minimally, still with O1≥O1s.Step 3.
If the removed configurationΘR is different from the added one ΘA 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.
Letkij=0 for all the joint angles θij.Step 2.
The new feasible domainθijξ is obtained by (30) for all the joint angles θij.Step 3.
Find a joint angleθij with kij<K in the configuration set Φn, which can be replaced by one selected from its own feasible domain θijξ and yield a new configuration set with the highest observability index O3 with O1≥O1s, and then let kij=kij+1.Step 4.
Ifkij for each joint angle in the configuration set Φn satisfies kij≥K, this step stops and the optimal experiment design for robot calibration is obtained. Otherwise, go on to Step 2.


## 5. Simulation Study

There are two purposes of simulation study here: the first is to investigate the influence of measurement configurations on the indexesO1 and O3; 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 O1 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.


### 5.1. 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 f∗ decides not only the contents of the kinematic parameters but also the concrete expression of the function f∗.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 p=5) including αi,ai,θi,di,βi to describe the two adjoining link coordinate frames, and the homogeneous transformation between them is shown as follows:(31)Aii-1=rotx,αi·transx,ai·rotz,θi·transz,di·roty,βi,where rot and trans denote the transpose and rotation matrix, respectively, and Aii-1 denotes the homogeneous transformation of joint i with respect to joint i-1. For a d DOF robot without a tool coordinate frame, the end-effector pose can be expressed by the homogeneous transformation as follows:(32)Ad0=∑i=1dAii-1,where Ad0 is the homogeneous transformation of the end-effector with respect to the base coordinate frame, which can be transferred into Z-Y-X 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 q=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.Figure 2
The coordinate frame established by MDH method.As is clear from Figure3, three kinds of typical robots with different DOF are adopted to analyze the effects of measurement configurations on observability indexes, whose 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].Table 1
Kinematic parameters of the 3-DOF robot.
Joint α (∘) a (m) θ (∘) d (m) β (∘) 1 0 0 0 0.1 - 2 - 90 0 0 0 - 3 0 0.5 0 - 0Table 2
Kinematic parameters of the 4-DOF robot.
Joint α (∘) a (m) θ (∘) d (m) β (∘) 1 0 0 0 0.1 - 2 - 90 0 0 - 0 3 0 0.5 0 - 0 4 0 0.5 0 - 0Table 3
Kinematic parameters of the 5-DOF robot.
Joint α (∘) a (m) θ (∘) d (m) β (∘) 1 0 0 0 0 - 2 - 90 0 0 - 0 3 0 0.45 0 - 0 4 - 90 0.02 0 0.5 - 5 90 0 0 0 -Figure 3
Robots with different DOF for simulation: (a) 3-DOF robot, (b) 4-DOF and (c) 5-DOF.
(a) (b) (c)


### 5.2. 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 Table4, where d is the number of DOF, m denotes the number of the independent parameters to be identified, nL denotes the minimal size of the measurement configuration set, and nM is set as the maximum size for convenience.Table 4
Parameters given for robot calibration.
d m n L n M 3 12 2 50 4 15 3 50 5 20 4 50Measurement configuration sets of different sizes fromnL to nM are constructed by random selection from the candidate configuration set pool ξ, and their corresponding observability indexes O1 and O3 are then calculated.Results are shown in Figures4, 5, and 6, where the optimal size of the configuration set and the corresponding limit O1s 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 O1 increases firstly and then converges to a threshold, while the index O3 always has an increasing trend.Figure 4
The influence of the size of the configuration set on the indexesO1 and O3 for the 3-DOF robot.Figure 5
The influence of the size of the configuration set on the indexesO1 and O3 for the 4-DOF robot.Figure 6
The influence of the size of the configuration set on the indexesO1 and O3 for the 5-DOF robot.Taking the 4-DOF robot as an example, the difference betweenO1 and O3 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 O1 and O3 are calculated.In order to observe the relationship between the two indexes, 50 configuration sets are sorted in order of the value of their indexO1. 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:(33)VP=OImax-OIminOImin×100%,where OImax denotes the maximum value of the observability index appearing in the 50 samples, OImin denotes the minimum one, and then VP is an evaluation index of the magnitude of variation in the observability index.Figure 7
The difference betweenO1 and O3 on the same configuration set with n=5.Figure 8
The difference betweenO1 and O3 on the same configuration set with n=10.Figure 9
The difference betweenO1 and O3 on the same configuration set with n=20.As shown in Figures7, 8, and 9, the varying magnitude of the indexes O1 and O3 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 O1 is always much smaller than that of O3, though both of them decrease as the size of the measurement configuration set increases. Moreover, the observability indexes O1 and O3 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 O3 than on the index O1.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 withO1 maximal economically while the last three indicate the maximization of O3 with O1 to be acceptable so as to produce a win-win situation.


### 5.3. 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 Figure10, 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.Figure 10
Flowchart of calibration simulation.Taking the 4-DOF robot as an example, the angle errors are set as0.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 n=10, the lower limit is O1s=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 N0,10-3, whose unit is mm for the length values and rad for the angle values, while the encoder noises are defined with N0,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 O1 by DETMAX, and random selection.( 3 ) The calibration was carried out, respectively, under the three measurement configuration set obtained above for 500 times.


### 5.4. 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:(34)PΩi=Ωci-Ω∗iδΩi×100%,i=1,2,…,15,where δΩ denotes the difference between the actual parameters Ω and the nominal parameters Ω∗, and Ωc denotes the calibrated parameters Ωc, so PΩi is the identification accuracy of the ith kinematic parameter. For one calibration, the mean PΩand varPΩof the vector PΩ 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:(35)δD=Poc-Poa=δx2+δy2+δz2,where Poc denotes the estimate end-effector position of the robot, Poa denotes the actual position, and δx,y,z are the deviations in x, y, and z 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 Figures11 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.Figure 11
The estimate errors of the end-effector position with respect to different sizes of configuration sets.Figure 12
The identification accuracy of kinematic parameters with respect to different sizes of configuration sets.As is clear from Figures11 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 O1 and O3 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.Table 5
The configuration set constructed with the method proposed in this paper.
No. θ 1 (∘) θ 2 (∘) θ 3 (∘) θ 4 (∘) 1 143.4 68 - 122 90 2 96 - 60 169.5 - 36 3 - 105 150 30 - 30 4 40 90 178 70 5 134.4 180 - 60 180 6 90 -60 70 - 144 7 148.5 152 170 180 8 - 92 -40 165 0 9 - 90 180 180 180 10 148.5 0 - 158 30Table 6
The configuration set constructed with the method maximizingO1 by DETMAX.
No. θ 1 (∘) θ 2 (∘) θ 3 (∘) θ 4 (∘) 1 143.4 68 - 122 90 2 96 - 60 169.5 - 36 3 - 105 150 30 - 30 4 40 90 178 70 5 134.4 180 - 60 180 6 90 -60 70 - 144 7 148.5 152 170 180 8 - 92 - 40 165 0 9 - 90 180 180 180 10 148.5 0 - 158 30Table 7
The configuration set constructed with random selection.
No. θ 1 (∘) θ 2 (∘) θ 3 (∘) θ 4 (∘) 1 143.4 68 - 122 90 2 96 - 60 169.5 - 36 3 - 105 150 30 - 30 4 40 90 178 70 5 134.4 180 - 60 180 6 90 - 60 70 - 144 7 148.5 152 170 180 8 - 92 - 40 165 0 9 - 90 180 180 180 10 148.5 0 - 158 30Calibration under these three configuration sets was repeated 500 times, and typical results are shown in Figures13 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.Figure 13
The calibration performances in the position estimate errors (a) and the identified parameters errors (b).
(a) (b)Figure 14
The calibration performances in the position estimate errors (a) and the identified parameters errors (b).
(a) (b)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 Figures15 and 16, and it can be concluded that the method maximizing O1 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.Figure 15
Distribution of the means (a) and variances (b) of the position estimate errors in 500 simulations.
(a) (b)Figure 16
Distribution of the means (a) and variances (b) of the identified parameters errors in 500 simulations.
(a) (b)


## 5.1. 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 f∗ decides not only the contents of the kinematic parameters but also the concrete expression of the function f∗.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 p=5) including αi,ai,θi,di,βi to describe the two adjoining link coordinate frames, and the homogeneous transformation between them is shown as follows:(31)Aii-1=rotx,αi·transx,ai·rotz,θi·transz,di·roty,βi,where rot and trans denote the transpose and rotation matrix, respectively, and Aii-1 denotes the homogeneous transformation of joint i with respect to joint i-1. For a d DOF robot without a tool coordinate frame, the end-effector pose can be expressed by the homogeneous transformation as follows:(32)Ad0=∑i=1dAii-1,where Ad0 is the homogeneous transformation of the end-effector with respect to the base coordinate frame, which can be transferred into Z-Y-X 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 q=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.Figure 2
The coordinate frame established by MDH method.As is clear from Figure3, three kinds of typical robots with different DOF are adopted to analyze the effects of measurement configurations on observability indexes, whose 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].Table 1
Kinematic parameters of the 3-DOF robot.
Joint α (∘) a (m) θ (∘) d (m) β (∘) 1 0 0 0 0.1 - 2 - 90 0 0 0 - 3 0 0.5 0 - 0Table 2
Kinematic parameters of the 4-DOF robot.
Joint α (∘) a (m) θ (∘) d (m) β (∘) 1 0 0 0 0.1 - 2 - 90 0 0 - 0 3 0 0.5 0 - 0 4 0 0.5 0 - 0Table 3
Kinematic parameters of the 5-DOF robot.
Joint α (∘) a (m) θ (∘) d (m) β (∘) 1 0 0 0 0 - 2 - 90 0 0 - 0 3 0 0.45 0 - 0 4 - 90 0.02 0 0.5 - 5 90 0 0 0 -Figure 3
Robots with different DOF for simulation: (a) 3-DOF robot, (b) 4-DOF and (c) 5-DOF.
(a) (b) (c)


## 5.2. 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 Table4, where d is the number of DOF, m denotes the number of the independent parameters to be identified, nL denotes the minimal size of the measurement configuration set, and nM is set as the maximum size for convenience.Table 4
Parameters given for robot calibration.
d m n L n M 3 12 2 50 4 15 3 50 5 20 4 50Measurement configuration sets of different sizes fromnL to nM are constructed by random selection from the candidate configuration set pool ξ, and their corresponding observability indexes O1 and O3 are then calculated.Results are shown in Figures4, 5, and 6, where the optimal size of the configuration set and the corresponding limit O1s 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 O1 increases firstly and then converges to a threshold, while the index O3 always has an increasing trend.Figure 4
The influence of the size of the configuration set on the indexesO1 and O3 for the 3-DOF robot.Figure 5
The influence of the size of the configuration set on the indexesO1 and O3 for the 4-DOF robot.Figure 6
The influence of the size of the configuration set on the indexesO1 and O3 for the 5-DOF robot.Taking the 4-DOF robot as an example, the difference betweenO1 and O3 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 O1 and O3 are calculated.In order to observe the relationship between the two indexes, 50 configuration sets are sorted in order of the value of their indexO1. 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:(33)VP=OImax-OIminOImin×100%,where OImax denotes the maximum value of the observability index appearing in the 50 samples, OImin denotes the minimum one, and then VP is an evaluation index of the magnitude of variation in the observability index.Figure 7
The difference betweenO1 and O3 on the same configuration set with n=5.Figure 8
The difference betweenO1 and O3 on the same configuration set with n=10.Figure 9
The difference betweenO1 and O3 on the same configuration set with n=20.As shown in Figures7, 8, and 9, the varying magnitude of the indexes O1 and O3 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 O1 is always much smaller than that of O3, though both of them decrease as the size of the measurement configuration set increases. Moreover, the observability indexes O1 and O3 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 O3 than on the index O1.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 withO1 maximal economically while the last three indicate the maximization of O3 with O1 to be acceptable so as to produce a win-win situation.


## 5.3. 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 Figure10, 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.Figure 10
Flowchart of calibration simulation.Taking the 4-DOF robot as an example, the angle errors are set as0.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 n=10, the lower limit is O1s=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 N0,10-3, whose unit is mm for the length values and rad for the angle values, while the encoder noises are defined with N0,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 O1 by DETMAX, and random selection.( 3 ) The calibration was carried out, respectively, under the three measurement configuration set obtained above for 500 times.


## 5.4. 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:(34)PΩi=Ωci-Ω∗iδΩi×100%,i=1,2,…,15,where δΩ denotes the difference between the actual parameters Ω and the nominal parameters Ω∗, and Ωc denotes the calibrated parameters Ωc, so PΩi is the identification accuracy of the ith kinematic parameter. For one calibration, the mean PΩand varPΩof the vector PΩ 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:(35)δD=Poc-Poa=δx2+δy2+δz2,where Poc denotes the estimate end-effector position of the robot, Poa denotes the actual position, and δx,y,z are the deviations in x, y, and z 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 Figures11 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.Figure 11
The estimate errors of the end-effector position with respect to different sizes of configuration sets.Figure 12
The identification accuracy of kinematic parameters with respect to different sizes of configuration sets.As is clear from Figures11 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 O1 and O3 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.Table 5
The configuration set constructed with the method proposed in this paper.
No. θ 1 (∘) θ 2 (∘) θ 3 (∘) θ 4 (∘) 1 143.4 68 - 122 90 2 96 - 60 169.5 - 36 3 - 105 150 30 - 30 4 40 90 178 70 5 134.4 180 - 60 180 6 90 -60 70 - 144 7 148.5 152 170 180 8 - 92 -40 165 0 9 - 90 180 180 180 10 148.5 0 - 158 30Table 6
The configuration set constructed with the method maximizingO1 by DETMAX.
No. θ 1 (∘) θ 2 (∘) θ 3 (∘) θ 4 (∘) 1 143.4 68 - 122 90 2 96 - 60 169.5 - 36 3 - 105 150 30 - 30 4 40 90 178 70 5 134.4 180 - 60 180 6 90 -60 70 - 144 7 148.5 152 170 180 8 - 92 - 40 165 0 9 - 90 180 180 180 10 148.5 0 - 158 30Table 7
The configuration set constructed with random selection.
No. θ 1 (∘) θ 2 (∘) θ 3 (∘) θ 4 (∘) 1 143.4 68 - 122 90 2 96 - 60 169.5 - 36 3 - 105 150 30 - 30 4 40 90 178 70 5 134.4 180 - 60 180 6 90 - 60 70 - 144 7 148.5 152 170 180 8 - 92 - 40 165 0 9 - 90 180 180 180 10 148.5 0 - 158 30Calibration under these three configuration sets was repeated 500 times, and typical results are shown in Figures13 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.Figure 13
The calibration performances in the position estimate errors (a) and the identified parameters errors (b).
(a) (b)Figure 14
The calibration performances in the position estimate errors (a) and the identified parameters errors (b).
(a) (b)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 Figures15 and 16, and it can be concluded that the method maximizing O1 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.Figure 15
Distribution of the means (a) and variances (b) of the position estimate errors in 500 simulations.
(a) (b)Figure 16
Distribution of the means (a) and variances (b) of the identified parameters errors in 500 simulations.
(a) (b)


## 6. 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 withO3 to be the objective and O1 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.



---

*Source: 4689710-2018-05-20.xml*