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

Alexandr Klimchik, Benoit Furet, Anatol Pashkevich

► To cite this version:

Alexandr Klimchik, Benoit Furet, Anatol Pashkevich. Algebraic technique for the stiffness model reduction in elastostatic calibration of robotic manipulators. 2014 IEEE/ASME International Conference on Advanced Intelligent Mechatronics, Jul 2014, Besançon, France. pp.1710-1715, 10.1109/AIM.2014.6878330. hal-01206894

HAL Id: hal-01206894

https://imt-atlantique.hal.science/hal-01206894

Submitted on 6 Jul 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Algebraic technique for the stiffness model reduction in elastostatic calibration of robotic manipulators

Alexandr Klimchik, Benoit Furet, Anatol Pashkevich

**Abstract**— 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.

**Keywords**— Robotic manipulator, elastostatic calibration, parameter identifiability, model reduction, stiffness modeling,

## I. INTRODUCTION

Application area of robotic manipulators in industry is progressively increasing, they become attractive in many technological processes, including precise machining. In such a process, the robot is usually subject to essential external loadings caused by the machining forces. These forces may lead to significant deviation of the end-effector position. This feature becomes extremely important when the accuracy requirements are very high. In this case, stiffness modeling of robotic manipulator and related compliance error compensation become the key issues allowing user to improve the positioning accuracy [1-4].

At present, there are three major methods for the manipulator stiffness modeling: the Finite Element Analysis (FEA), the Matrix Structural Analysis (MSA), and the Virtual Joint Method (VJM). Among them, the VJM-based technique provides reasonable trade-off between the model accuracy and computational complexity, it is able to describe the elastic properties of the manipulator elements by 6x6 compliance matrices. Because of its advantages it will be further used in this paper. However, evaluation or identification of the stiffness model parameters is not a trivial problem that has not been paid enough attention yet.

The main difficulty of elastostatic calibration is that direct application of the VJM-based method gives excessive number of parameters whose impact in the robot positioning accuracy essentially differ (mainly because of their

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 non-trivial. 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

$$ \Delta t = \sum_{i=1}^{n} (\mathbf{J}_{0}^{(i)} \mathbf{k}_{0}^{(i)} \mathbf{J}_{0}^{(i)T}) \cdot \mathbf{w} \quad (1) $$

where $\Delta t$ is the vector of the end-effector displacements under the loading $\mathbf{w}$, the matrices $\mathbf{k}_{0}^{(i)} = (\mathbf{K}_{0}^{(i)})^{-1}$ denote the link/joint compliances that should be identified by means of calibration, and the matrices $\mathbf{J}_{0}^{(i)}$ are corresponding sub-Jacobians obtained by the fractioning of the aggregated Jacobian $\mathbf{J}_{0}^{T} = [\mathbf{J}_{0}^{(1)T}, \mathbf{J}_{0}^{(2)T}, \dots]$. For the identification purposes, this expression should be transformed into the form, where all desired parameters (elements of the matrices $\mathbf{k}_{0}^{(i)}, i = 1, 2, \dots$) are collected in a single vector $\pi = (k_{011}^{(1)}, k_{012}^{(1)}, \dots, k_{066}^{(n)})$. It yields the following linear equation

$$ \Delta t = A(\mathbf{q}, \mathbf{w} | \pi) \cdot \pi \quad (2) $$

$$ A(\mathbf{q}, \mathbf{w} | \pi) = [\mathbf{J}_1 \mathbf{J}_1^T \mathbf{w}, \mathbf{J}_2 \mathbf{J}_2^T \mathbf{w}, \dots, \mathbf{J}_m \mathbf{J}_m^T \mathbf{w}] \quad (3) $$

is so-called observation matrix that defines the mapping between the unknown compliances $\pi$ and the end-effector

<footnote>*Research supported by ANR (Project ANR-2010-SEGI-003-02-COROUSSO), France and FEDER ROBOTEX, France.</footnote>

<footnote>Klimchik and A. Pashkevich are with Ecole des Mines de Nantes, France and with Institut de Recherches en Communications et Cybernétique de Nantes (IRCCyN), France (Tel.+33-251858300; fax.+33-251858349; e-mail: alexandr.klimchik@mines-nantes.fr, anatol.pashkevich@mines-nantes.fr). B. Furet is with University of Nantes, France and with IRCCyN (e-mail: benoit.furet@univ-nantes.fr)</footnote>
displacements $\Delta\mathbf{t}$ under the loading $\mathbf{w}$ for the manipulator configuration $\mathbf{q}, \mathbf{J}_i$ is the column of Jacobian matrix $\mathbf{J}_0$.

Taking into account that the calibration experiments are carried out for several manipulator configurations defined by the actuated joint coordinates $\mathbf{q}_j, j=1,m$, the system of basic equations for the identification can be presented in the following form

$$ \Delta\mathbf{t}_j = \mathbf{A}(\mathbf{q}_j, \mathbf{w}_j | \boldsymbol{\pi}) \cdot \boldsymbol{\pi} + \boldsymbol{\varepsilon}_j; \quad j = \overline{1,m} \quad (4) $$

where $\boldsymbol{\varepsilon}_j$ denotes the vector of measurement errors. For further convenience, let us also present equation (4) in a the matrix form

$$ \Delta\mathbf{t}_a = \mathbf{A}_a(\mathbf{q}_a, \mathbf{w}_a | \boldsymbol{\pi}) \cdot \boldsymbol{\pi} + \boldsymbol{\varepsilon}_a \quad (5) $$

where the subscript "a" indicates that matrices/vectors aggregate corresponding components for $m$ configurations. Below, for the presentation convenient the matrix $\mathbf{A}_a(\mathbf{q}_a, \mathbf{w}_a | \boldsymbol{\pi})$ will be also referred to as $\mathbf{A}_{\pi}$, 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

$$ F = \sum_{j=1}^{m} (\mathbf{A}_{\pi j} \boldsymbol{\pi} - \Delta\mathbf{t}_j)^T \boldsymbol{\eta}^T \boldsymbol{\eta} (\mathbf{A}_{\pi j} \boldsymbol{\pi} - \Delta\mathbf{t}_j) \rightarrow \min_{\boldsymbol{\pi}} \quad (6) $$

where $\boldsymbol{\eta}$ is the matrix of weighting coefficients that normalizes the measurement data, $\mathbf{A}_{\pi j} = \mathbf{A}(\mathbf{q}_j, \mathbf{w}_j | \boldsymbol{\pi})$. This minimization problem yields the following solution

$$ \hat{\boldsymbol{\pi}} = \left( \sum_{j=1}^{m} \mathbf{A}_{\pi j}^{T} \boldsymbol{\eta}^{T} \boldsymbol{\eta} \mathbf{A}_{\pi j} \right)^{-1} \cdot \left( \sum_{j=1}^{m} \mathbf{A}_{\pi j}^{T} \boldsymbol{\eta}^{T} \boldsymbol{\eta} \Delta\mathbf{t}_{\pi j} \right). \quad (7) $$

If the measurement noise is Gaussian (as it is assumed in conventional calibration techniques), expression (7) provides us with a unbiased estimates for which $E(\hat{\boldsymbol{\pi}}) = \boldsymbol{\pi}$.

It can be proved [8] that the best results in terms of the identification accuracy are achieved if $\boldsymbol{\eta} = \boldsymbol{\Sigma}^{-1}$, where the matrix $\boldsymbol{\Sigma}^2 = E(\boldsymbol{\varepsilon} \cdot \boldsymbol{\varepsilon}^T)$ describes the statistical properties of the measurement errors. It leads to the following covariance matrix of the manipulator compliance parameters

$$ \mathrm{cov}(\hat{\boldsymbol{\pi}}) = \left( \sum_{j=1}^{m} \mathbf{A}_{\pi j}^{T} \boldsymbol{\Sigma}^{-T} \boldsymbol{\Sigma}^{-1} \mathbf{A}_{\pi j} \right)^{-1} \quad (8) $$

Such assignment of the weighing coefficients $\boldsymbol{\eta}$ 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 $\{\mathbf{q}_j, j=1,m\}$ that is also known as the calibration experiment planning [9], which directly influences on the observation matrices $\mathbf{A}(\mathbf{q}_j, \mathbf{w}_j | \boldsymbol{\pi})$ and on the covariance matrix (8).

It is clear that expression (7) gives reliable estimates of the parameters $\boldsymbol{\pi}$ if and only if the matrix $\sum_{j=1}^{m} \mathbf{A}_{\pi j}^{T} \boldsymbol{\Sigma}^{-T} \boldsymbol{\Sigma}^{-1} \mathbf{A}_{\pi j}$ 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 sub-sections, 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 $\boldsymbol{\pi}$ 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 $\mathbf{A}(\mathbf{q}_j, \mathbf{w}_j | \boldsymbol{\pi})$ are computed for certain set of measurement configurations $\{\mathbf{q}_j\}$ and loadings $\{\mathbf{w}_j\}$. Depending on the matrix set $\{\mathbf{A}_{\pi j}\}$, corresponding system of linear equations can be solved for $\boldsymbol{\pi}$ either uniquely or may have infinite number of solutions. In general, if the information matrix is rank-deficient, a general solution of the system (4) can be presented in the following form

$$ \hat{\boldsymbol{\pi}} = \mathbf{A}_{\Sigma}^{+} \cdot \mathbf{B}_{\Sigma} + (\mathbf{I} - \mathbf{A}_{\Sigma}^{+} \mathbf{A}_{\Sigma}) \cdot \boldsymbol{\lambda}. \quad (9) $$

where $\mathbf{A}_{\Sigma} = \sum_{j=1}^{m} \mathbf{A}_{\pi j}^{T} \boldsymbol{\eta}^{T} \boldsymbol{\eta} \mathbf{A}_{\pi j}$, $\mathbf{B}_{\Sigma} = \sum_{j=1}^{m} \mathbf{A}_{\pi j}^{T} \boldsymbol{\eta}^{T} \boldsymbol{\eta} \Delta\mathbf{t}_{\pi j}$, the superscript $+^m$ denotes the Moore-Penrose pseudoinverse and $\boldsymbol{\lambda}$ is an arbitrary vector of the same size as $\boldsymbol{\pi}$. Using the later expression, all desired parameters contained in the vector $\boldsymbol{\pi}$ can be divided into the following groups [6]:

**G1:** Identifiable parameters that are independent from the arbitrary vector $\boldsymbol{\lambda}$ 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

$$ v_i = |\hat{\pi}_i| / \sigma_i, \quad i = 1, 2, \dots \quad (10) $$

where $\sigma_i$ is the standard deviation of the parameter estimate $\hat{\pi}_i$ extracted from the diagonal of the covariance matrix (8).
It is clear that $v_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: $v_i > v_0^+$;

**G1-: Practically non-identifiable parameters**, for which the accuracy indicator is low: $v_i < v_0^-$;

**G1~: Practically semi-identifiable parameters**, for which the accuracy indicator is intermediate: $v_0^- \le v_i \le v_0^+$.

An open question however is related to justified assigning of the upper and lower bounds $v_0^+$ and $v_0^-$. From practical point of view, it is reasonable to use $v_0^+ = 5$ and $v_0^- = 2$, 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: PHYSICAL APPROACH ($\pi \to \pi'$)

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 $\{k_{ij}, k_{ji}\}$ by a single one $k_{ij}, i < j$.

For each link, this procedure is equivalent to re-definition of the model parameters vector in the following way

$$ \pi = M \cdot \pi' \qquad (11) $$

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 $k_{ij}$ corresponding to zeros in the stiffness matrix template $k^0$ 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

$$ k^0 = \begin{bmatrix} * & 0 & 0 & 0 & 0 & 0 \\ 0 & * & 0 & 0 & 0 & * \\ 0 & 0 & * & 0 & * & 0 \\ 0 & 0 & 0 & * & 0 & 0 \\ 0 & 0 & * & 0 & * & 0 \\ 0 & * & 0 & 0 & 0 & * \end{bmatrix} \qquad (12) $$

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 ($\pi' \to \pi''$)

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 non-redundant and the number of parameters is equal to the
observation matrix rank). In practice, the following inequality is often satisfied: $rank(A_π) < dim(π')$. 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:

$$ \mathrm{rank}(\mathbf{A}_{\pi}) = \mathrm{rank}(\mathbf{A}_{\pi'}) = \mathrm{dim}(\pi'') \quad (13) $$

They are based on the partitioning of the parameters set $π'$ into three non-overlapping groups (identifiable, non-identifiable 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 $A_a(q_a, w_a | π')$ as the product of three matrices $UΣV^T$. Here $U=[U_1,...U_m]$ and $V=[V_1,...V_v]$ are orthogonal matrices of the size $m×m$ and $n×n$ respectively whose columns are denoted as $U_i$ and $V_j$; the second factor $Σ$ is an rectangular diagonal matrix of the size $m×n$ containing $r$ positive real numbers $σ_1, σ_2,..., σ_r$ in descending order; $m=dim(Δt_a)$ is the number of rows in the observation matrix (i.e. number of equations used for the identification), $n=dim(π')$ is current number of the model parameters, and $r$ is the rank of the aggregated observation matrix $A_a(q_a, w_a | π')$. It is clear that $r$ defines the maximum number of parameters that can be identified using given set of configurations {$q_i$} and wrenches {$w_i$}.

Further, after substitution of this decomposition into (5) and left-multiplication by $U^T$, the system of $m$ identification equations (5) can be rewritten as $ΣV^T·π' = U^T·Δt_a$ or, in more detailed form,

$$ \left[ \begin{array}{cc|c} \sigma_1 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \sigma_r \\ \hline 0 & \cdots & \sigma_{r-1} \times r \\ \hline 0 & \cdots & \sigma_{r-1} \times (n-r) \end{array} \right] \cdot \left[ \begin{array}{c} \mathbf{V}_1^T \\ \vdots \\ \mathbf{V}_n^T \end{array} \right] \cdot \pi' = \left[ \begin{array}{c} \mathbf{U}_1^T \\ \vdots \\ \mathbf{U}_m^T \end{array} \right] \cdot \Delta t_a \quad (14) $$

where the number of equation is equal to $n$ and corresponds to the dimension of vector $π'$ (it is obvious that $n≤m$). 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

$$ \begin{aligned} \sigma_i \mathbf{V}_i^T \cdot \pi' &= \mathbf{U}_i^T \cdot \Delta t_a; & i &= 1, 2, \dots, r \\ \mathbf{0} \cdot \pi' &= \mathbf{U}_i^T \cdot \Delta t_a; & i &= r+1, \dots, n \end{aligned} \quad (15) $$

where the second group of $m-r$ 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 $U_i^T·Δt_a = \mathbf{0}$ for $i>r$ if the measurement vector $Δt_a$ 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

$$ \pi'_{0} = \sum_{i=1}^{r} \frac{1}{\sigma_{i}} \mathbf{V}_{i} \cdot \mathbf{U}_{i}^{T} \cdot \Delta \mathbf{t}_{a} \quad (16) $$

This allows us to present general solution (9) as the sum of this partial solution and arbitrary vector from the subspace with the basis $V_1, V_2, ..., V_n$

$$ \hat{\pi}' = \pi'_{o} + \sum_{i=r+1}^{n} \lambda_{i} \mathbf{V}_{i} \quad (17) $$

where $λ_i$ are arbitrary real values.

It is clear that all equalities in (15) should be satisfied for any displacements $Δt_a$ (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

$$ \mathbf{V}_i^T \cdot \pi' = 0; \quad i = \overline{1, r} \quad (18) $$

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 $V_i$) all model parameters $π'$ can be partitioned into three groups G1, G2 and G3. Parameters of group G1 (identifiable) uniquely defines by the equation (17) and do not depend on the arbitrary values $λ_i$, for these parameters corresponding row of the sub-matrix $[V_{r+1}, ..., V_m]$ is equal to zero. Parameters of group G2 (non-identifiable) do not effect on the $Δt_j$, for them corresponding column of the sub-matrix $[V_1, ..., V_r]$ is equal to zero. Parameters of group G3 (semi-identifiable) effect on the $Δt_j$ but cannot be identified uniquely, couplings between the parameters of this group is defined by the vectors $V_i$, 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$, $i=1, dim(π')$:

Rule 1: Include parameter $π'_i$ into the group G1 if the $i$th row of the sub-matrix $[V_{r+1}, ..., V_m]$ is equal to zero;

Rule 2: Include parameter $π'_i$ into the group G2 if the $i$th column of the sub-matrix $[V_1, ..., V_r]$ is equal to zero;

Rule 3: If the parameter $π'_i$ is not included in G1 or G2, include it in the group G3.

**M4b:Elimination.** Eliminate from the set of unknowns (model parameters) non-identifiable parameters that correspond to group G2.

As a result of these techniques the set of model parameters $π'$ will be reduced to the set $π' \setminus π_{G2}$ that does not reduce the rank of observation matrix $rank(A_{π' \setminus π_{G2}}) = rank(A_{π'})$, but which still is redundant, i.e. $rank(A_{π' \setminus π_{G2}}) < dim(π' \setminus π_{G2})$. It should be noted that $rank(A_{π' \setminus π_{G2}}) = dim(π_{G1})$ while $rank(A_{π_{G3}}) < dim(π_{G3})$. So, another and the most difficult problem that arises after M4, is to define among the set of parameters $π_{G3}$ 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 $\pi_{G3}$ 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 $\pi_{G3}$ into the non-overlapping groups of coupled parameters $\pi_{G3}$ for which the following conditions are satisfied:

$$ (a) \quad \pi_{G3} = \pi_{G3}^1 \cup \pi_{G3}^2 \cup \dots \cup \pi_{G3}^m, \quad \pi_{G3}^i \cap \pi_{G3}^j = \emptyset \quad \forall i \neq j; $$

$$ (b) \quad \mathrm{rank}(\mathbf{A}_{\pi_{G3}^i}) = \mathrm{rank}(\mathbf{A}_{\pi_{G3}^i, \pi_{G3}^j}), \quad \forall j=1:\dim(\pi_{G3}^i) $$

$$ (c) \quad \mathrm{rank}(\mathbf{A}_{\pi_{G3}^i}) < \mathrm{rank}(\mathbf{A}_{\pi_{G3}^i, \pi_{G3}^j}), \quad \forall i \neq j, \\ k = 1:\dim(\pi_{G3}^i) $$

In practice, when grouping is not evident it is possible to use numerical technique, which is based on the SVD-decomposition of the observation matrix $\mathbf{A}_a(\mathbf{q}_a, \mathbf{w}_a | \pi_{G3})$

$$ \mathbf{U} \cdot \mathbf{\Sigma} \cdot \mathbf{V}^T = \mathbf{A}_a (\mathbf{q}_a, \mathbf{w}_a | \pi_{G3}) \quad (19) $$

where the matrix $\mathbf{V}$ can be presented as $\mathbf{V}=[\mathbf{V}_{r-}, \mathbf{V}_{r+}]$ in accordance with the rank of $\mathbf{A}_a(\mathbf{q}_a, \mathbf{w}_a | \pi_{G3})$. Here, matrix $\mathbf{V}_{r-}$ 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

$$ \mathbf{L} = (\mathbf{V}_{r-}^*)^{-1} \mathbf{V}_{r-}^T \quad (20) $$

where matrix $\mathbf{V}_{r-}^*$ is full-rank square matrix composed of the column of the matrix $\mathbf{V}_{r-}^T$. It should be noted that matrix $\mathbf{V}_{r-}^*$ is not unique, nevertheless, any of them allows getting matrix of couplings between the model parameters $\mathbf{L}$ in a sparse form. The groups of elements can be easily detected after transformation of the matrix $\mathbf{L}$ into the block-diagonal form. Here, the sub-set of parameters $\pi_{G3}^j$ corresponds to jth block of the matrix $\mathbf{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 $\pi_{G3}^j$ and froze remainder:

**M5b:Selection.** In each group of parameters $\pi_{G3}^j$ specify $n_j = \mathrm{rank}(\mathbf{A}_a(\mathbf{q}_a, \mathbf{w}_a | \pi_{G3}^i))$ parameters that will be treated as identifiable

**M5c:Assigning.** In each group of parameters $\pi_{G3}^j$ fix $m_j = \dim(\pi_{G3}^j) - \mathrm{rank}(\mathbf{A}_a(\mathbf{q}_a, \mathbf{w}_a | \pi_{G3}^j))$ parameters to some constants; these parameters will be treated as non-identifiable

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 $\pi_{G3}^j$. 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 $\pi_{G3}$ will be split into two subsets: the subset of the parameters that will be treated as identifiable $\pi_{G3}^i$ and subset that will be

treated as non identifiable ones $\pi_{G3}^i$ and will be assigned to some constant values ($\pi_{G3}^i = \mathrm{const}$); i.e. $\pi_{G3}^i \cup \pi_{G3}^i = \pi_{G3}$, $\pi_{G3}^i \cap \pi_{G3}^i = \emptyset$.

As the output of the algebraic approach, the complete set of parameters $\pi'$ will be reduced to $\pi''$ that includes all parameters from group G1 and identifiable one from the group G3, i.e. $\pi'' = \pi_{G1} \cup \pi_{G3}^i$. Nevertheless, this model reduction will not violate completeness of the model i.e. $\mathrm{rank}(\mathbf{A}_{\pi'}) = \mathrm{rank}(\mathbf{A}_{\pi'})$.

## VI. MODEL REDUCTION: STATISTICAL APPROACH($\pi'' \rightarrow \pi'''$)

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 $\hat{\pi}'$ and their covariance matrix $\mathrm{cov}(\hat{\pi}')$ by means of (7) and (8);

**Step 2:** Using the parameters estimates $\hat{\pi}'$ and the diagonal elements of the covariance matrix $\mathrm{cov}(\hat{\pi}')$, compute the parameter-to-noise ratios $v_i$ in accordance with expression (10);

**Step 3.** For all compliance matrices $\mathbf{k}$ to be identified, eliminate from the set of unknowns the parameters $k_{ij}$ for which parameter-to-noise ratios $v_{ij}$ is lower the user defined threshold: $v_{ij} < v_0^*$.

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].

<figure><img src="image_2.png" /><figcaption>Figure 1. Experimental setup for the elastostatic calibration</figcaption></figure>

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).

<table><caption>TABLE I. SEQUENTIAL REDUCTION OF ELASTOSTATIC MODEL FOR 6 DOF INDUSTRIAL ROBOT KUKA KR-270</caption><thead><tr><th>Step</th><th>Model description</th><th>Number of parameters</th></tr></thead><tbody><tr><td>Original model</td><td>6 joints +7 links (36 par. for each)</td><td>258</td></tr><tr><td>M1: Symmetrisation</td><td>6 joints +7 links (21 par. for each)</td><td>153</td></tr><tr><td>M2: Sparcing</td><td>6 joints +7 links (8 par. for each)</td><td>62</td></tr><tr><td>M3: Aggregation</td><td>7 links (8 par. for each)</td><td>56</td></tr><tr><td>M4a:Partitioning</td><td>G1: Identifiable parameters – 4</td><td rowspan="3">52</td></tr><tr><td>M4b:Elimination</td><td>G2: Non-identifiable par. - 4</td></tr><tr><td></td><td>G3: Semi-identifiable par. – 48</td></tr><tr><td>M5a:Splitting</td><td rowspan="3">28 parameters in one group of 48 semi-identifiable parameters</td><td rowspan="3">32</td></tr><tr><td>M5b:Selection</td></tr><tr><td>M5c:Assigning</td></tr><tr><td>M6: Neglecting</td><td>Results depend on step M5</td><td>25-30</td></tr></tbody></table>

<table><caption>TABLE II. PRINCIPAL ELASTIC PARAMETERS OF ROBOT KUKA KR-270</caption><thead><tr><th>Parameter</th><th>Value</th><th>CI</th><th>Parameter-to-noise ratio</th></tr></thead><tbody><tr><td>k<sub>2</sub>, [rad×µm/N]</td><td>0.302</td><td>±0.004 (1.3%)</td><td>231</td></tr><tr><td>k<sub>3</sub>, [rad×µm/N]</td><td>0.406</td><td>±0.008 (2.0%)</td><td>150</td></tr><tr><td>k<sub>4</sub>, [rad×µm/N]</td><td>3.002</td><td>±0.115 (3.8%)</td><td>79</td></tr><tr><td>k<sub>5</sub>, [rad×µm/N]</td><td>3.303</td><td>±0.162 (4.9%)</td><td>61</td></tr><tr><td>k<sub>6</sub>, [rad×µm/N]</td><td>2.365</td><td>±0.095 (4.0%)</td><td>75</td></tr></tbody></table>

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.

## ACKNOWLEDGMENT

The work presented in this paper was partially funded by ANR (Project ANR-2010-SEGI-003-02-COROUSSO) and FEDER ROBOTEX, France

## REFERENCES

[1] J. Kövecses and J. Angeles, "The stiffness matrix in elastically articulated rigid-body systems," Multibody System Dynamics, vol. 18, pp. 169-184, 2007.

[2] G. Alici and B. Shirinzadeh, "Enhanced stiffness modeling, identification and characterization for robot manipulators," Robotics, IEEE Transactions on, vol. 21, pp. 554-564, 2005.

[3] A. Pashkevich, A. Klimchik, and D. Chablat, "Enhanced stiffness modeling of manipulators with passive joints," Mechanism and machine theory, vol. 46, pp. 662-679, 2011.

[4] A. Klimchik, A. Pashkevich, D. Chablat, and G. Hovland, "Compliance error compensation technique for parallel robots composed of non-perfect serial chains," Robotics and Computer-Integrated Manufacturing, vol. 29, pp. 385-393, 4// 2013.

[5] W. Khalil and E. Dombre, Modeling, identification and control of robots: Butterworth-Heinemann, 2004.

[6] A. Pashkevich, "Computer-aided generation of complete irreducible models for robotic manipulators," in The 3rd Int. Conference of Modelling and Simulation. pp. 293-298.

[7] A. Klimchik, S. Caro, and A. Pashkevich, "Practical identifiability of the manipulator link stiffness parameters," ASME 2013 International Mechanical Engineering Congress & Exposition, pp. 1-10, 2013.

[8] J.A. Rice. Mathematical Statistics and Data Analysis. Second Edition: Duxbury Press, 1995.

[9] J.-H. Borm and C.-H. Menq, "Determination of Optimal Measurement Configurations for Robot Calibration Based on Observability Measure," The Int. J. of Robotics Research, vol. 10, pp. 51-63, 1991.

[10] S. Timoshenko, J. N. Goodier, Theory of elasticity, 3d ed.: McGraw-Hill, New York, 1970.

[11] W. Khalil, M. Gautier, and J. Kleinfinger, "Automatic generation of identification models of robots," International Journal of Robotics and Automation, vol. 1, pp. 2-6, 1986.

[12] A. Klimchik, Y. Wu, A. Pashkevich, S. Caro, and B. Furet, "Optimal Selection of Measurement Configurations for Stiffness Model Calibration of Anthropomorphic Manipulators," Applied Mechanics and Materials, vol. 162, pp. 161-170, 2012.