Algebraic technique for the stiffness model reduction in elastostatic calibration of robotic manipulators

The paper deals with elastostatic calibration of a serial industrial robot. In contrast to other works, all compliance sources associated with both links and joints elasticity are taken into account. Particular attention is paid to the model parameters identification using end-point measurements only. For such experimental setup, the model is transformed into the form suitable for calibration with the sufficient rank of the corresponding observation matrix. The main contributions are in developing algebraic, physical and heuristic techniques that allow user to obtain complete model with minimal number of parameters. The advantages of the developed approach are confirmed by an experimental study that deals with identification of the elastostatic model parameters for a 6 dof serial industrial robot.

magnitude).Furthermore, straightforward application of this technique produces redundant model that is not suitable for calibration.Particularly, the whole set of the elastostatic parameters for 6 dof manipulator contains 258 values and leads to the fail of the numerical routines used in calibration that is caused by singularity of the relevant observation matrix.Similar problem arises in geometric case where the concept of complete-irreducible-continues model has been introduced.To obtain such a model relevant algebraic tools for the model reduction have been developed [5][6].This idea can be also used in the elastic case, however, here complex couplings between the model elements make this step nontrivial.Besides, in elastostatic calibration, an additional difficulty is caused by essential difference in the parameter magnitudes.As follows from our experience, corresponding identification results could violate fundamental physical properties of the stiffness matrices.For this reason, a new notion of practical identifiability has been introduced in our previous work [7] that allows obtaining reliable results in real industrial environment.This paper further develops this research and pays particular attention to generation theoretically identifiable elastostatic model.In the case of elastostatic calibration, this step is more complicated comparing with the geometric one (because of complex coupling between the model elements) and requires additional investigations.

II. METHODOLOGY OF ELASTOSTATIC IDENTIFICATION
To estimate the matrices describing elasticity of the manipulator components (i.e., compliances of the virtual springs presented in the VJM-based modeling [3]), the elastostatic model can be written as [ , ,...] T T T  J J J .For the identification purposes, this expression should be transformed into the form, where all desired parameters (elements of the matrices

., ]
T T m m T  A q w J J w J J w J J w π   is so-called observation matrix that defines the mapping between the unknown compliances π and the end-effector

Algebraic technique for the stiffness model reduction in elastostatic calibration of robotic manipulators
Alexandr Klimchik, Benoit Furet, Anatol Pashkevich displacements t under the loading w for the manipulator configuration q , i J is the column of Jacobian matrix θ J .Taking into account that the calibration experiments are carried out for several manipulator configurations defined by the actuated joint coordinates , 1, j jm  q , the system of basic equations for the identification can be presented in the following form where j ε denotes the vector of measurement errors.For further convenience, let us also present equation ( 4) in a the matrix form where the subscript "a" indicates that matrices/vectors aggregate corresponding components for m configurations.Below, for the presentation convenient the matrix ( , | ) a a a A q w π will be also referred to as  A , where the subscript defines the parameters set for which the observation matrix is computed.Further, using these notations and assigning proper weights for each equation, the identification can be reduced to the following optimization problem where η is the matrix of weighting coefficients that normalizes the measurement data, ( , | )   A q w A π .This minimization problem yields the following solution If the measurement noise is Gaussian (as it is assumed in conventional calibration techniques), expression (7) provides us with a unbiased estimates for which   Ê  ππ .It can be proved [8] that the best results in terms of the identification accuracy are achieved if Such assignment of the weighing coefficients η also allows us to avoid the problem of different units in the objective function (6), which arises in straightforward application of the least-square technique to the robot parameters identification if the measurement system provides both position and orientation data.It should be noted that this particularity is usually omitted in conventional robot calibration.Another way to improve the identification accuracy is related to the proper selection of manipulator measurement configurations   , 1, j jm  q that is also known as the calibration experiment planning [9], which directly influences on the observation matrices ( , | ) A q w π and on the covariance matrix (8).
It is clear that expression (7) gives reliable estimates of the parameters π if and only if the matrix 1 1 m TT j j j     A Σ Σ A is invertible.It leads to the problem of the parameter identifiability that have been studied by a number of authors for the problem of geometrical calibration [5][6].Relevant techniques are based on the information matrix rank analysis (via either SVD-or QR-decomposition).However, in real industrial environment where the measurement not is non-negligible, the identifiable parameters are not equivalent in terms of accuracy (both absolute and relative) and expression (7) can give rather surprising results for some of them.This motivates revision of the above mentioned notion (parameter identifiability) and its extension taking into account the identification accuracy defined by the covariance matrix (8).In the following subsections, the notion of practical identifiability is introduced.Besides, existing model reduction techniques (algebraically approach) have been developed mainly for geometrical calibration.Although, the main idea from geometric calibration can be also used in the case of elastostatic one, numerical routines should be revised since in elastic calibration couplings between the model elements are more complex.

III. BASIC ASSUMPTIONS AND TERMINOLOGY
Let us assume that the vector of desired elastostatic parameters π should be identified from the set of the linear equations (4) whose least square solution is defined by the expression (7), where the observation matrices ( , | ) A q w π are computed for certain set of measurement configurations   j q and loadings   j w .Depending on the matrix set   j  A , corresponding system of linear equations can be solved for π either uniquely or may have infinite number of solutions.In general, if the information matrix is rankdeficient, a general solution of the system (4) can be presented in the following form A t , the superscript "+" denotes the Moore-Penrose pseudoinverse and λ is an arbitrary vector of the same size as π .Using the later expression, all desired parameters contained in the vector π can be divided into the following groups [6] : G1: Identifiable parameters that are independent from the arbitrary vector λ and can be obtained uniquely from (9); G2: Non-identifiable parameters that can take on any value without influence on the right-hand side of the equation ( 2) and cannot be computed from ( 9) in a unique way; G3: Semi-identifiable parameters that have influence on the right-hand side of the equation ( 2), but cannot be computed uniquely from (9).
In this paper, in contrast to previous works, this classification is enhanced taking into account practical issues related to the limited precision of the measurement system.The main idea is to compare the absolute value of the estimated parameter with the range of possible fluctuations of the estimate caused by the measurement noise.For computational reasons, it is convenient to introduce a numerical indicator similar to the signal-to-noise ratio in communication, which is defined as follows where i  is the standard deviation of the parameter estimate ˆi  extracted from the diagonal of the covariance matrix (8).
It is clear that i  can be treated as the inverse of the relative accuracy, which allows us to avoid the problem of division by zero.In the following sections this indicator will be referred to as parameter-to-noise ratio.
Using the above defined indicator, the set of parameters belonging to the group G1 (theoretically identifiable) can be further divided into three subgroups: G1+: Practically identifiable parameters, for which the accuracy indicator is high: G1-: Practically non-identifiable parameters, for which the accuracy indicator is low:  ; G1~: Practically semi-identifiable parameters, for which the accuracy indicator is intermediate:  .An open question however is related to justified assigning of the upper and lower bounds 0   and 0   .From practical point of view, it is reasonable to use 0 5 which is in a good agreement with the quantiles of the normal distribution.However, the user may modify these values in accordance with the specificity of the problem of interest.
The above presented definitions allow us to revise the concept of "suitable-for-calibration" model that in previous works included all parameters of the group G1.In this work, this model is limited to include only parameters of the subgroup G1+ that can be estimated with reasonable accuracy and provide good approximation of the original complete model.The following subsections address different aspects of model reduction allowing us to obtain the desired model suitable for the elastostatic calibration.

IV. MODEL REDUCTION
Straightforward approach to the manipulator stiffness modelling leads to the exhaustive but redundant number of parameters to be identified.For instance, each links is described by a 6 6  matrix that includes 36 parameters that are treated as independent ones.However, as follows from physics, number of the pure physical and independent parameters is essentially lower.Hence, there are strong relations between these 36 parameters but this fact is usually ignored in elastostatic calibration.Besides, due to fundamental properties of conservative system, the desired compliance matrices should be strictly symmetrical and positive-definite.In addition, for typical manipulator links, the compliance matrices are sparsed due to the shape symmetry with respect to some axis, but this property is also not taken into account in identification of the elastostatic parameters.
To take advantages of the compliance matrix properties and to increase the identification accuracy, three simple methods can be applied that allows us to reduce the number of parameters to be computed in the identification procedure (7).They can be treated as the physics-based model reduction techniques and formalized in the following way.
M1:Symmetrisation.For all compliance matrices k to be identified, replace the pairs of symmetrical parameters kk by a single one , ij k i j  .For each link, this procedure is equivalent to re-definition of the model parameters vector in the following way where the binary matrix M of size 36 21  describes the mapping from the original to reduced parameter space.It is clear that this idea allows us to reduce the number of links compliance parameters from 36 to 21 (and from 258 to 153 for the entire 6 d.o.f.manipulator).
M2:Sparcing.For all compliance matrices k , eliminate from the set of unknowns the parameters ij k corresponding to zeros in the stiffness matrix template 0 k derived analytically for the manipulator link with similar shape.
To derive desired template matrix it is convenient to use any realistic link-shape approximation.For example using the trivial beam [10], the desired template can be presented as  0 *00000 0*000* 00*0*0 000*00 00*0*0 0*000* where the symbol "*" denotes non-zero elements.It allows further reducing the number of the unknown parameters from 21 to 8, taking into account only essential ones from physical point of view.It can be also proved that the template ( 12) is valid for any link whose geometrical shape is symmetrical with respect to three orthogonal axes.But it is necessary to be careful if this property is not kept strictly.
All joint compliances cannot be identified separately and should be included in the compliance matrix of previous link by means of modification corresponding diagonal elements.These modifications touch the values that are identified and elimination corresponding columns in the observation matrix M3:Aggregation.Eliminate from the set of model parameters the ones that correspond to joint compliances before which there is a compliant link.In terms of parameters identifiability the compliance of such joints cannot be separated from the links.
Summarizing these methods, it should be mentioned that the above presented methods essentially reduce the number of parameters to be identified, but they do not violate the mode completeness, i.e. the ability to describe any deflection caused by the loading.However, the reduced model may still have some redundancy in the frame of entire manipulator.
From the geometrical calibration it is known that in spite of the fact that redundant model is suitable for direct and inverse computations it cannot be used in identification since the observation matrix does not have sufficient rank.Similar problem arises in elastostatic calibration where some stiffness coefficients of adjacent links cannot be identified separately.The problem of construction complete and irreducible model has been widely studied in geometrical calibration but have never been in the focus of elastostatic calibration.
V. MODEL REDUCTION: ALGEBRAIC APPROACH (    ππ ) The physical approach described in previous sub-section allows us essentially reducing the number of model parameters.However, it does not guarantee that the obtained model is suitable for calibrations (i.e. that the model is nonredundant and the number of parameters is equal to the observation matrix rank).In practice, the following inequality is often satisfied: To overcome the problem, this sub-section presents some algebraic tools aimed at further reduction of the model parameter set from  π to  π , which ensures full identifiability: They are based on the partitioning of the parameters set  π into three non-overlapping groups (identifiable, nonidentifiable and semi-identifiable), which either eliminated from the model or reduced to ensure the equality (13).
To introduce relevant algebraic technique, let us apply the group singular value decomposition (SVD) and present the aggregated observation matrix   dim a m  t is the number of rows in the observation matrix (i.e.number of equations used for the identification), π is current number of the model parameters, and r is the rank of the aggregated observation matrix

 
,| a a a  A q w π .It is clear that r defines the maximum number of parameters that can be identified using given set of configurations   i q and wrenches   i w .Further, after substitution of this decomposition into (5) and left-multiplication by T U , the system of m identification equations ( 5) can be rewritten as where the number of equation is equal to n and corresponds to the dimension of vector  π (it is obvious that nm  ).Taking into account particularities of the sparse matrix Σ (with r non-zero elements only) it is possible to rewrite the system (14) in the form where the second group of mr  equations should be excluded from further consideration because relevant residuals do not depend on the parameters of interest  π (they are multiplied by zero matrix).It can be proved that • a T i   t U 0 for ir  if the measurement vector a t does not contain noise.It is also worth mentioning for real identification problems (with the identification noise) the second group of equations produces constant residuals in the least-square objective (6) that cannot be minimized by varying the vector of unknown parameters  π .Hence, for the identification of n parameters comprising the vector  π we have obtained a system of r linear equations that obviously cannot be solved uniquely.Its partial solution can be found as This allows us to present general solution (9) as the sum of this partial solution and arbitrary vector from the subspace with the basis 12 , ,...
where i  are arbitrary real values.
It is clear that all equalities in (15) should be satisfied for any displacements a t (that obviously have effect also on the matrices U , Σ and V ).To analyze the model parameters identifiability, special interest represents assigning all displacements to zeros.This allows us to rewrite equation (15) as that highlights the relative impacts of model parameters on the end-effector displacement and can be used to obtain coupling between the unknowns.So, in accordance with the properties of the second term of equation ( 17) (values of vectors i V ) all model parameters  π can be partitioned into three groups G1, G2 and G3.Parameters of group G1 (identifiable) uniquely defines by the equation ( 17

V
V is equal to zero.Parameters of group G2 (non-identifiable) do not effect on the j t , for them corresponding column of the sub-matrix 1 [ ,..., ] r VV is equal to zero.Parameters of group G3 (semi- identifiable) effect on the j t but cannot be identified uniquely, couplings between the parameters of this group is defined by the vectors i V , which have more than one non- zero elements.Based on these decomposition algebraic-based model reduction techniques can be formalised in the following way: M4a:Partitioning.Divide the reduced set of the model parameters  π into three non-overlapping groups G1, G2 and G3 in accordance with the following rules applied to all i   , So, another and the most difficult problem that arises after M4, is to define among the set of parameters 3 G π the one that will be treated as identifiable and which will be fixed.
The problem of selection identifiable parameters among the set of semi-identifiable has infinite number of solutions and usually is solved numerically, without attempts to solve in a general way.Here, let us presents an algorithm that is able to split parameters 3 G π into the non-overlapping groups of coupled parameters in general way and then choose identifiable one from the group based on their physical scene M5a:Splitting.Split the set of semi-identifiable model parameters 3 G π into the non-overlapping groups of coupled parameters 3 j G π for which the following conditions are satisfied: (a) ) k  π In practice, when grouping is not evident it is possible to use numerical technique, which is based on the SVDdecomposition of the observation matrix where the matrix V can be presented as [ , ] A q w π .Here, matrix r  V defines coupling between the elements.To highlight the groups of elements it is required to sparse it.One of the most efficient ways to do this is to multiply it on the inverse of where matrix * r  V is full-rank square matrix composed of the column of the matrix T r  V .It should be noted that matrix * r  V is not unique, nevertheless, any of them allows getting matrix of couplings between the model parameters L in a sparse form.The groups of elements can be easily detected after transformation of the matrix L into the block-diagonal form.Here, the sub-set of parameters 3 j G π corresponds to jth block of the matrix L .In each diagonal block the number of lines defines the number of parameters from the block that should be included in the complete model suitable for the identification.
On the next step it is required to select parameters that will be identifiable in each set of parameters 3 j G π and froze remainder: M5b:Selection.In each group of parameters A q w π parameters that will be treated as identifiable M5c:Assigning.In each group of parameters A q w π parameters to some constants; these parameters will be treated as nonidentifiable It should be noted that the order of methods M5b and M5c is not strict, identifiable and non-identifiable parameters can be selected and fixed iteratively addressing to the methods M5b and M5c several times in the frame of one set of parameters 3 j G π .The sequence of steps and final set of parameters highly depend on the user, who can define priority for different parameters [11].As the result of methods M5a, M5b and M5c the set of parameters 3 G π will be split into two subsets: the subset of the parameters that will be treated as identifiable 3 id G π and subset that will be treated as non identifiable ones 3 ni G π and will be assigned to some constant values ( 3 As the output of the algebraic approach, the complete set of parameters  π will be reduced to  π that includes all parameters from group G1 and identifiable one from the group G3, i.e.As follows from relevant study, rigorous reduction methods based on the physical and mathematical properties of the compliance matrix are rather limited if the measurement noise is non-negligible.This gives us reasons to develop some heuristic rules that take into account the measurement noise impact on the identification accuracy.It is clear that extremely low accuracy is not acceptable, but often corresponding parameters are so small that their influence on the end-effector deflections is almost negligible.This supports an idea for heuristic reduction of small model parameters but leaving an open problem of their further reconstruction in the VJM-model using some empirical or semi-empirical relations induced by mathematical relations between the stiffness matrix elements. To take into account the relative accuracy of the parameter estimates, it is convenient to use a simple indicator showing parameter-to-noise ratio (10).It is evident that it should be applied only to those parameters that belong to the group G1.Using this index, a heuristic model reduction technique allowing us to distinguish the practically identifiable parameters from the hardly-identifiable ones can be formalised as follows: M6:Neglecting.
Step 1: Using complete but non-redundant model, compute estimates of the desired parameters ˆ π and their covariance matrix   ĉov  π by means of ( 7) and (8); Step 2: Using the parameters estimates ˆ π and the diagonal elements of the covariance matrix   ĉov  π , compute the parameter-to-noise ratios i  in accordance with expression (10); Step 3.For all compliance matrices k to be identified, eliminate from the set of unknowns the parameters ij k for which parameter-to-noise ratios ij  is lower the user defined threshold: This method allows us to eliminate from the model the parameters whose identification accuracy is comparable with the noise and cannot be considered as reliable estimates.
As follows from our experience, it a very powerful method with two useful features: (i) elimination of small (but theoretically non-zero) parameters, and (ii) detection of the elements corresponding to zeros in the matrix template (see method M2), if the latter has been defined rather carefully.

VII. APPLICATION EXAMPLE
The developed model reduction techniques have been applied to the generation sophisticated elastostatic model of the industrial robot KUKA KR-270 (Figure 1) that is suitable for the identification in real industrial environment.The application area of this robot requires high precision while performing the technological process, which generates essential deflections of the end-effector.If the elastic properties are known, these deflections can be compensated on the control level via adjusting a target trajectory [4].The considered manipulator contains 7 links separated by 6 actuated joints.Taking into account that in general elastostatic properties of each link are defined by 6x6 stiffness matrix, the complete but obviously redundant model contains 258 parameters.As a result of application model reduction techniques (M1-M5), the number of parameters to be identified has been reduced down to 32.More details on each step are given in Table I.It should be noted that the number of parameters in practically identifiable model depends on the step M5 and vary from 25 to 30.Using relevant model the compliance error compensation for the robotic based milling essential improvement of the precision has been achieved.The obtained elastic model allowed us to compensate more than 95% of deflections caused by external loading and to guarantee precision about 0.1 mm (that is comparable with the robot repeatability 0.06 mm).Additional analysis of the reduced model (that takes into account the actuators compliance only) has been done based on the proposed parameter-to-noise ratio.Relevant results have been presented in [12], for which parameter-to-noise ratios have been computed (Table II).As follows from them, in such a model all parameters are practically identifiable, however it is not able to describe all deflections (especially in x and y directions) because of the model limitations.Hence, the complete model presented in this paper has obvious advantages over the reduced model studied before.

VIII. CONCLUSION
The paper deals with the problem of elastostatic calibration of serial industrial robots, which is extremely important for robotic-based machining when high position accuracy is required.Particular attention is paid to the model parameters identification using end-point measurements only.The main contributions are in developing algebraic, that allow user to obtain complete model with minimal number of parameters, which is suitable for elastostatic calibration.The advantages of the developed approach have been illustrated by an application example that deals with the elastostatic parameters identification of industrial robot.
are orthogonal matrices of the size m m  and n n respectively whose columns are denoted as i U and j V ; the second factor Σ is an rectangular diagonal matrix of the size m n  containing r positive real numbers 12 , ,... r    in descending order;

1 [
) and do not depend on the arbitrary values i  , for these parameters corresponding row of the sub-matrix

Figure 1 .
Figure 1.Experimental setup for the elastostatic calibration

TABLE I .
SEQUENTIAL REDUCTION OF ELASTOSTATIC MODEL FOR 6 DOF INDUSTRIAL ROBOT KUKA KR-270