2.1. Manipulator stiffness modeling

A typical serial articulated robot is depicted in Fig. 1(a) and its generic elastostatic architecture (Fig. 1(b)) is comprised of a fixed base and a serial chain of flexible links joined by rotational flexible actuated joints and an end-effector, over which a wrench w (force/torque) is acting on.

Figure 1. (a) Manipulator KUKA^® KR-500 L480 and (b) its complete VJM representation.

The VJM establishes that, besides the 1-d.o.f. virtual springs commonly included in the actuated joints to model the joint compliance [Reference Dumas, Caro, Garnier and Furet20], the original rigid model should be extended by adding virtual joints (localized springs with up to 6 d.o.f.) to describe the elastic deformations of the links (Fig. 1(b)) [Reference Pashkevich, Klimchik and Chablat23]. As in refs. [Reference Pashkevich22, Reference Klimchik, Furet, Caro and Pashkevich36], the equivalent extended kinematic chain is described by:

1. (i) a rigid link between the robot’s base and the first actuated joint modeled as the constant homogenous transformation matrix T _Base;

2. (ii) flexible actuated joints modeled as a homogeneous matrix function ${\bf{T}}_{{\rm{Joint}}}^i\left( {{q^i} + \theta _{Ac}^i} \right)$ , which depends on the actuated joint variable q ⁱ and the virtual joint variable $\theta _{Ac}^i$ that takes into account the joint compliance;

3. (iii) a set of rigid links, described by the constant homogenous transformation matrices ${\bf{T}}_{{\rm{Link}}}^i$ ;

4. (iv) a set of 6-d.o.f. virtual joints which model the link flexibility and are described by the homogeneous matrix function ${{\bf{T}}_{{\rm{VJM}}}}\left( {{\boldsymbol{\theta }}_{{\bf{Link}}}^{\bf{i}}} \right)$ , which depends on the virtual joint variables ${\boldsymbol{\theta }}_{{\bf{Link}}}^{\bf{i}} = \left( {\theta _x^i,\theta _y^i,\theta _z^i,\theta _{\varphi x}^i,\theta _{\varphi y}^i,\theta _{\varphi z}^i} \right)$ corresponding to the translation/rotation deflections in/around the axis x, y, z;

5. (v) a rigid link from the last joint to the end-effector, modeled as a constant homogenous matrix transformation T _Tool.

The final robot pose is modeled by the homogeneous transformation T shown in Eq. (1), which defines the end-effector location subject to the variations of all joint coordinates (virtual and actuated joints) [Reference Klimchik, Furet, Caro and Pashkevich36]:

(1) \begin{align}{\bf{T}} = {{\bf{T}}_{{\rm{Base}}}} \cdot \left[ {\prod \nolimits_{i = 1}^n {\bf{T}}_{{\rm{Joint}}}^i\left( {{q^i} + \theta _{{\rm{Ac}}}^i} \right) \cdot {{\bf{T}}_{{\rm{Link}}}} \cdot {{\bf{T}}_{{\rm{VJM}}}}\left( {\theta _{{\rm{Link}}}^i} \right)} \right] \cdot {{\bf{T}}_{{\rm{Tool}}}}\end{align}

where n is the number of links/joints (6, in the present case). For the sake of simplification, the rotation and translation components can be extracted from the homogeneous matrix T [Reference Craig41], so that the equivalent geometric model can be rewritten as shown in Eq. (2).

(2) \begin{align}{\bf{X}} = {\bf{g}}\left( {{\bf{q}},{\rm{\;}}{\boldsymbol{\theta }}} \right)\end{align}

where g( $ \cdot $ ) denotes a relevant vector function and the vector X = (p, $\boldsymbol{\varphi}$ ) ^T defines both the end-effector tool frame position p = (x, y, z) ^T and orientation $\boldsymbol{\varphi}$ = ( ${\varphi}$ _x , ${\varphi}$ _y , ${\varphi}$ _z ) ^T in relation to the robot world frame. The vector q = (q ₁ , q ₂ , …,q _n ) ^T is composed by all actuated coordinates defined by the robot controller, whereas the vector $\boldsymbol{\theta}$ = ( ${{\boldsymbol{\theta }}_1}$ , ${{\boldsymbol{\theta }}_2}$ , …, ${{\boldsymbol{\theta }}_{n\theta }}$ ) ^T consists of all virtual joint coordinates from all the n _θ virtual joints whose values are dependent on the external loading/wrench applied to the end-effector.

The representation of the manipulator stiffness behavior considers the end-effector displacement ΔX and the external wrench w as the input and output of the model, respectively. Thus, the manipulator stiffness equation is written as

(3) \begin{align}{\bf{w}} = {{\bf{K}}_{\rm{C}}} \cdot \Delta {\bf{X}}\end{align}

where K _C is the manipulator Cartesian stiffness matrix for a given configuration q, subjected to the geometric constraints given in Eq. (2) [Reference Pashkevich22].

An alternative expression for the Cartesian stiffness matrix can be developed by assuming the so-called generalized Hooke’s law for the manipulator [Reference Salisbury42]. According to that, variations in the virtual joint variables $\boldsymbol{\theta}$ are linearly related to the virtual force/torque reactions applied to the corresponding links/joints. This means τ _θ = K _θ $\cdot$ Δ $\boldsymbol{\theta}$ , where τ _θ = ( $\boldsymbol{\tau}$ _θ1 , $\boldsymbol{\tau}$ _θ2 , …, $\boldsymbol{\tau}$ _θnθ ) ^T is the vector of the virtual joint reactions and K _θ = diag (K _θ1 , K _θ2 , …, K _θnθ ) is the virtual spring stiffness matrix; K _θi represents the elastostatic property of the corresponding link/joint. Moreover, assuming small arbitrary displacements of the virtual joints Δ $\boldsymbol{\theta}$ around the equilibrium neighborhood, the corresponding deviations of the end-effector are given by ΔX = J _θ $\cdot$ Δ $\boldsymbol{\theta}$ , where J _θ = ∂g(q, $\boldsymbol{\theta}$ )/∂ $\boldsymbol{\theta}$ is the Jacobian with respect to the virtual variables $\boldsymbol{\theta}$ , which may be computed from Eq. (2) analytically or semi-analytically [Reference Pashkevich22, Reference Pashkevich, Klimchik and Chablat23].

Thus, applying the principle of virtual work and taking into account the previous considerations, and solving the resultant statics equations for w, it yields the so-called congruence transformation between the manipulator Cartesian stiffness matrix and the joint stiffness matrix [Reference Pashkevich, Klimchik and Chablat23],

(4) \begin{align}{{\bf{K}}_{\rm{C}}} = \left( {{{\bf{J}}_{\rm{\theta }}}^{ - {\rm{T}}} \cdot {{\bf{K}}_{\rm{\theta }}} \cdot {{\bf{J}}_{\rm{\theta }}}^{ - 1}} \right)\end{align}

The reciprocal relation describing a better manipulator compliance behavior results from Eqs. (3) and (4), assuming the external wrench w and the end-effector displacement ΔX as the input and output of the model, respectively, yielding

(5) \begin{align}\Delta {\bf{X}} = {\bf{C}} \cdot {\bf{w}}\end{align}

(6) \begin{align}{\bf{C}} = \left( {{{\bf{J}}_{\rm{\theta }}} \cdot {{\bf{C}}_{\rm{\theta }}} \cdot {{\bf{J}}_{\rm{\theta }}}^T} \right)\end{align}

where C = K _C ^–1 is the aggregated Cartesian compliance matrix for a given configuration q and C _θ = K _θ ^–1 is the aggregated joint compliance matrix of the manipulator. Since the matrices {K _θi , i = 1, 2, …, n _θ } (or their reciprocal joint compliance matrices {C _θi =(K _θi )^–1, i = 1, 2, …, n _θ }) are usually unknown, in order to calculate Eq. (4), they must be identified [Reference Klimchik, Furet, Caro and Pashkevich36]. This is addressed in Sections 3 and 4. In the sequence, some relevant physical and mathematical aspects related to the elastostatic models are reviewed.

2.2. Properties of the elastostatic model

The stiffness model given by Eq. (4) is based upon the following conservative assumption from the elasticity theory: (i) the force due to the stiffness matrix is integrable and conservative and (ii) the work done by this force is conservative, that is, it depends only on the initial/end potential energy states [Reference Timoshenko and Goodier43].

While the joint stiffness matrix K _θ can be manually prescribed in relation to the conservative nature of the stiffness modeling, Chen and Kao [Reference Chen and Kao44] demonstrated that Eq. (4) for computing K _C is only valid at the unloaded equilibrium configuration; otherwise, the compliance error compensation violates the conservative properties. They showed that it is a special case of the conservative congruence transformation, written as

(7) \begin{align}{{\bf{K}}_{\rm{C}}} = {\bf J_\theta }^{ - T} \cdot \left( {{{\bf{K}}_\theta } - {{\bf{K}}_w}} \right) \cdot {{\bf{J}}_\theta }^{ - 1}\end{align}

where K _w = ∂² [g ^T (q, $\boldsymbol{\theta}$ )w]/∂ $\boldsymbol{\theta}^{2}$ is a Hessian-derivative term that accounts for the incremental changes in the robot arm geometry due to the application of the external wrench. Despite of not being practical for identification purposes, Eq. (7) evinces that, while applying Eq. (3), the robot stiffness behavior is non-linear in relation to w.

Numerically and more practically, an n $ \times $ n stiffness (or compliance) matrix K is said to be conservative if it is (i) symmetric (k _ij = k _ji ) and (ii) positive-definite, where k _ij is the (i, j)-th element of K [Reference Chen and Kao44–Reference Kao and Ngo46]. Moreover, a definition and an equivalent assertion (whose proof is omitted herein) may apply for the positive-definiteness of a matrix [Reference Horn and Johnson38]:

These properties are extremely important to be considered for both the Cartesian (K _C) and the joint (K _θ ) stiffness matrices. In special, even though K _θ is supposed to be a fixed elastostatic property of the manipulator’s mechanical components, it is usual for an inattentive identification process to produce a viable numerical solution that violates positive-definiteness or symmetry [Reference Klimchik, Caro and Pashkevich35]. Considering this, a useful mathematical property of the stiffness matrices is presented in the next subsection.

2.3. Gramian nature of the stiffness matrices

Definition of the Gramian matrix, named after his author Jörgen Pedersen Gram (1850–1916) [Reference Horn and Johnson38]:

Several engineering fields use this definition, among them, the systems theory [Reference Sun, Yu and Anderson39, Reference Bianchin and Pasqualetti40]. In the framework of the lumped elastostatic stiffness modeling, the following theorem applies directly to the method proposed herein:

Since the symmetry of G is straightforwardly derived from Definition 3 and the commutativity of the inner product, the proof for the equivalence to the positive semi-definiteness is given as ref. [Reference Horn and Johnson38],

1. ✓ First part: given ${\bf{G}} = \left[ \langle{{{\bf{v}}_i},{{\bf{v}}_j}}\rangle \right] \in {\mathbb{R}^{n \times n}}$ and a not null vector $\textbf{u} \in {\mathbb{R}^{n \times 1}}$ , it follows

   (8) \begin{align}{{\bf{u}}^T}{\bf{Gu}} = \mathop \sum \nolimits_{i,j} \left\langle{{\bf{v}}_i},{{\bf{v}}_j}\right\rangle{{\rm{u}}_i}{{\rm{u}}_j} = \mathop \sum \nolimits_{i,j} \left\langle{{\bf{v}}_i}{{\rm{u}}_i},{{\bf{v}}_j}{{\rm{u}}_j}\right\rangle = \left\langle\mathop \sum \nolimits_i {{\bf{v}}_i}{{\rm{u}}_i},\mathop \sum \nolimits_j {{\bf{v}}_j}{{\rm{u}}_j}\right\rangle = \left\|\mathop \sum \nolimits_i {{\bf{v}}_i}{{\rm{u}}_i}^2 \ge 0\right\|\end{align}
   where it has been used the bi-linearity (2nd and 3rd equalities) and positive-definiteness (4th equality) of the inner product to evince the positive semi-definiteness of G.

2. ✓ Second part: given that G is symmetric and positive semi-definite

(9) \begin{align}{\bf{G}} = {\bf{RD}}{{\bf{R}}^{\textrm{T}}} = {\bf{R}}\left( {\sqrt {\bf{D}} {{\sqrt {\bf{D}} }^{\textrm{T}}}} \right){{\bf{R}}^{\textrm{T}}} = \left( {{\bf{R}}\sqrt {\bf{D}} } \right){\left( {{\bf{R}}\sqrt {\bf{D}} } \right)^{\textrm{T}}}\therefore {\bf{G}}\;{\rm{is\;the\;gramian\;of\;the\;}}{\bf{R}}\sqrt {\bf{D}} {\rm{\;rows}}\end{align}

because of the matrix spectral decomposition. A corollary (whose proof is found in ref. [Reference Horn and Johnson38]) from Theorem 1 regarding the positive-definite aspect is:

Theorem 1 implies that every conservative n $ \times $ n stiffness matrix is the Gramian matrix of some set of vectors or, equivalently and more usefully, that there should exist some n-dimensional set of vectors $\{{{\bf{v}}_j} \in {\mathbb{R}^m}\}$ able to generate it. This is used in this work to propose a novel identification method for C _θ (and reciprocally K _θ ).

4.1. Basic assumptions and definitions

Let us assume that, according to Theorem 1, there should exist an n-dimensional set of vectors { ${\bf{v}}_j^{\left( i \right)} \in {R^m},\;j = 1,\; \ldots ,\;n$ } able to generate each of the i-th n×n conservative compliance matrices {C _θi = (K _θi )^–1 , i = 1, 2, …n _θ }. Then, according to the Definition 3

(16) \begin{align}{{\bf{C}}_{\theta i}} = gramian\left( {{\bf{v}}_1^{\left( i \right)},{\bf{v}}_2^{\left( i \right)}, \ldots ,{\bf{v}}_n^{\left( i \right)}} \right)\end{align}

where the operator ‘gramian’ denotes the Gramian matrix generated over the i-th set { ${\bf{v}}_j^{\left( i \right)} \in {R^m},\;j = 1,\; \ldots ,\;n$ }. By rearranging Eq. (5)–(6) together and assuming that all the vectors from all the sets (the desired parameters in this case) can be collected together into the single ( $n\cdot m\cdot n$ _θ ) $ \times $ 1 optimization vector b = [ ${\bf{v}}_1^{\left( 1 \right)};{\bf{v}}_2^{\left( 1 \right)}; \ldots ;{\bf{v}}_n^{\left( {{n_\theta }} \right)}$ ], and further adapting the related equations for the N dataset points {q _u, w _u |ΔX _u, $u = 1, \ldots N$ },

(17) \begin{align}\Delta {{\bf{X}}_u}\left( {{{\bf{q}}_u},{{\bf{w}}_u}{\rm{|}}{\bf{b}}} \right) = \left[ {{{\bf{J}}_\theta } \cdot {{\bf{C}}_\theta }\left( {\bf{b}} \right) \cdot {{\bf{J}}_\theta }^T} \right] \cdot {{\bf{w}}_u} + {{\boldsymbol{\varepsilon }}_u}\;\;,\;\;u = 1, \ldots N\end{align}

where $\boldsymbol{\varepsilon}$ _u still denotes the vector of measurement errors. It has been made clear that C _θ (b) = diag (C _θ1 , C _θ2 , …, C _θnθ ) represents a non-linear mapping between the aggregated joint compliance matrix and the parameter optimization vector b according to Eq. (16). A convenient error vector can be derived for each u-th equation as

(18) \begin{align}{{\bf{e}}_u}\left( {{{\bf{q}}_{\rm{u}}},{{\bf{w}}_{\rm{u}}}{\rm{|}}{\bf{b}}} \right) = \Delta {{\bf{X}}_{\rm{u}}} - \left[ {{{\bf{J}}_{\rm{\theta }}}\left( {{{\bf{q}}_{\rm{u}}}} \right) \cdot {{\bf{C}}_{\rm{\theta }}}\left( {\bf{b}} \right) \cdot {{\bf{J}}_{\rm{\theta }}}{{\left( {{{\bf{q}}_{\rm{u}}}} \right)}^{\rm{T}}}} \right] \cdot {{\bf{w}}_{\rm{u}}}\;\;,\;\;{\rm{u}} = 1, \ldots {\rm{N}}\end{align}

where the error vector e _u = [δΔx _u, δΔy _u, δΔz _u, δΔ ${\varphi}$ _xu, δΔ ${\varphi}$ _yu, δΔ ${\varphi}$ _zu]^T stands for the accomplished error (δ) in the estimate of the end-effector deviation ΔX _u. Departing from this, the identification can be translated into a non-linear minimization problem defined over whole measurement dataset {q _u, w _u|ΔX _u, ${\rm{u}} = 1, \ldots {\rm{N}}$ } as

(19) \begin{align}\min_{\bf b}\left\{ {{{\rm{F}}_M} = \frac{1}{{\sqrt 2 {\sigma _p}}}\sqrt {\frac{{\mathop \sum \nolimits_{u = 1}^{\rm{N}} \left( {{\rm{\delta \Delta }}{x_u}^2 + {\rm{\delta \Delta }}{y_u}^2 + {\rm{\delta \Delta }}{z_u}^2} \right)}}{{\rm{N}}}} + \frac{1}{{\sqrt 2 {\sigma _\varphi }}}\sqrt {\frac{{\mathop \sum \nolimits_{u = 1}^{\rm{N}} \left( {{\rm{\delta \Delta }}{{\rm{\varphi }}_{xu}}^2 + {\rm{\delta \Delta }}{{\rm{\varphi }}_{yu}}^2 + {\rm{\delta \Delta }}{{\rm{\varphi }}_{zu}}^2} \right)}}{{\rm{N}}}} } \right\}\end{align}

where the cost function F _M is modified to weight differently between the translation and rotation RMS distance-based error portions by considering the measurement uncertainty: the former portion is weighted by the inverse of the position uncertainty σ_p, and the latter is weighted by the inverse of the orientation uncertainty σ $_{\varphi}$ . These quantities could be a specification from the measuring device manufacturer or empirically estimated [Reference Wu47]. Since most calibration processes require two measurements (with and without loading), these uncertainties are increased by the factor $\;\sqrt 2 $ [Reference Klimchik, Furet, Caro and Pashkevich36]. Besides of being a homogenized error-cost function, Eq. (19) leads to a theoretical minimum threshold limit for the optimization, namely 2. This would mean that the more the identification approaches this value, the more it is averagely close to the limit accuracy of the experimental data ( $\sqrt 2 $ σ_p and $\sqrt 2 $ σ $_{\varphi}$ ). Besides of being difficult to achieve, reaching a value too lower than 2 would imply on the possibility of the model getting biased by the noise. Therefore, it is enough to solve Eq. (19) for an objective limit slightly higher or close to 2.

4.2. Gramian-constrained optimization process

Under the light of the aforementioned definitions, Fig. 2 depicts the Gramian-constrained optimization (GCO) process for the stiffness model identification of industrial manipulators that is proposed herein.

Figure 2. Gramian-constrained optimization process for the manipulator stiffness model identification.

It should be stressed from Fig. 2 that the very beginning step of this process is the proper selection and parametrization of the VJM model, which means that at least some physical model reduction considerations should be applied prior to the GCO-based identification. Even though the technique proposed herein will inherently lead to a symmetric conservative compliance matrix, the application of some reduction steps to simplify the model could help to avoid convergence problems. For instance, the actuated joint compliances before which there is an elastic link should be eliminated from the set of model parameters and included in the compliance matrix of the previous link, a physical reduction step called aggregation [Reference Klimchik, Caro and Pashkevich35]. Moreover, some heavy-payload robots present some quite rigid links (comparatively to the others) and its compliances could be firstly disregarded. These simple considerations may help to reach a solution that is numerically and physically meaningful due to the Gramian constraint.

The remaining of the proposed approach is comprised of, using some relevant non-linear optimization technique (NLOT), iteratively changing the optimization vector b = [ ${\bf{v}}_1^{\left( 1 \right)};{\bf{v}}_2^{\left( 1 \right)}; \ldots ;{\bf{v}}_n^{\left( {{n_\theta }} \right)}$ ] until the desired threshold for F _M is reached. Some hybrid NLOT that applies heuristic/stochastic and deterministic methods may lead to more effective results. In this work, a genetic algorithm [Reference Goldberg49] is used to get the best initial guess b ₀, from which the pseudo-deterministic Global Search algorithm [Reference Ugray, Lasdon, Plummer, Glover, Kelly and Marti50] finds the optimal result b ^*, both from the Optimization Toolbox of the Matlab^® 2019b. The final compliance matrices are obtained from Eq. (16).

4.3. Order of the vector space

For all the n×n conservative compliance matrices C _θi , (n.m.n _θ ) parameters should be determined within the vector b. A guidance to select the order m of the vector space ${\bf{V}} \subset {\mathbb{R}^m}$ from which the vectors ${\bf{v}}_j^{\left( i \right)}$ are taken should be useful. About 6N scalar equations can be derived from the N measurement dataset points so that an overdetermined system is on ( $n\cdot m\cdot n$ _θ ) ≤ 6N. On the other hand, considering the symmetry property, each of the n $ \times $ n conservative compliance matrices C _θi presents in fact n $\cdot$ (n+1)/2 independent parameters. Thus, reasonably ( $n\cdot m$ ) ≥ (1+ n) $\cdot$ n/2, yielding to

(20) \begin{align}\left( {n + 1} \right)/2 \le m \le 6N/\left( {n \cdot {n_\theta }} \right)\end{align}

5.1. Motivation example

Let us use a numerical example extracted from ref. [Reference Klimchik, Caro and Pashkevich35] to illustrate the potential of the presented method. It is the elastostatic identification of a single link of a given robot, whose compliance matrix has been determined as

(21) \begin{align}C = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}{4.50 \cdot {{10}^{ - 8}}}&0&0&0&0&0\\ \\[-7pt] 0&{8.01 \cdot {{10}^{ - 5}}}&0&0&0&{3.98 \cdot {{10}^{ - 4}}}\\ \\[-7pt] 0&0&{{{3.64.10}^{ - 5}}}&0&{ - 1.71 \cdot {{10}^{ - 4}}}&0\\ \\[-7pt] 0&0&0&{3.76 \cdot {{10}^{ - 3}}}&0&0\\ \\[-7pt] 0&0&{ - 1.71 \cdot {{10}^{ - 4}}}&0&{1.09 \cdot {{10}^{ - 3}}}&0\\ \\[-7pt] 0&{3.98 \cdot {{10}^{ - 4}}}&0&0&0&{2.65 \cdot {{10}^{ - 3}}}\end{array}} \right]\end{align}

In this example, this compliance matrix is estimated according to a calibration process that follows both the classical least-squares and the GCO-based approaches. As in ref. [Reference Klimchik, Caro and Pashkevich35], data are generated by means of virtual experiments, where the link is assumed to be fixed on one side while an external wrench w _u is applied to the other side; six loading conditions (N = 6) were run, limiting the force and torque magnitude by 10 N and 10 N.m, respectively; the corresponding deflection is computed as ΔX _u = C $\cdot$ w _u + $\boldsymbol{\varepsilon}$ _u, where $\boldsymbol{\varepsilon}$ _u is the measurement noise, whose magnitude has been defined as σ _p = 25 μm for the position and as σ ${\varphi}_{\varphi}=0.25$ mrad for the orientation.

The obtained compliance matrices are presented in Eq. (22) and in Eq. (23) for the weighted least-squares and the GCO-based identification approaches, respectively.

(22) \begin{align}{C_{LS}} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} \mathbf{-7.32\cdot 10^{-7}}&-2.15\cdot 10^{-6}&2.05\cdot 10^{-6}&1.29\cdot 10^{-6}&-2.25\cdot 10^{-6}&\mathbf{-6.89\cdot 10^{-7}}\\ \\[-7pt] -1.58\cdot 10^{-6}&7.58\cdot 10^{-5}&-2.10\cdot 10^{-6}&6.52\cdot 10^{-7}&-9.35\cdot 10^{-7}&3.91\cdot 10^{-4}\\ \\[-7pt]-3.00\cdot 10^{-6}&1.79\cdot 10^{-6}&3.43\cdot 10^{-5}&-1.95\cdot 10^{-6}&-1.73\cdot 10^{-4}&1.03\cdot 10^{-6}\\ \\[-7pt] 8.39\cdot 10^{-6}&-1.74\cdot 10^{-5}&-3.59\cdot 10^{-5}&3.80\cdot 10^{-3}&-1.80\cdot 10^{-5}&2.50\cdot 10^{-5}\\ \\[-7pt] -3.16\cdot 10^{-5}&-1.19\cdot 10^{-5}&-1.79\cdot 10^{-4}&-1.72\cdot 10^{-5}&1.10\cdot 10^{-3}&1.85\cdot 10^{-5}\\ \\[-7pt] \mathbf{2.26\cdot 10^{-5}}&4.33\cdot 10^{-4}&2.66\cdot 10^{-5}&6.69\cdot 10^{-5}&3.21\cdot 10^{-5}&2.60\cdot 10^{-3} \end{array}} \right]\end{align}

(23) \begin{align}{C_{GCO}} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}\mathbf{1.12\cdot 10^{-7}}&-1.31\cdot 10^{-6}&3.06\cdot 10^{-7}&1.87\cdot 10^{-6}&-9.50\cdot 10^{-7}&\mathbf{5.49\cdot 10^{-7}}\\[3pt]-1.31\cdot 10^{-6}&7.62\cdot 10^{-5}&-8.24\cdot 10^{-7}&-8.45\cdot 10^{-7}&-4.65\cdot 10^{-7}&3.92\cdot 10^{-4}\\[3pt]3.06\cdot 10^{-7}&-8.24\cdot 10^{-7}&2.95\cdot 10^{-5}&-2.68\cdot 10^{-6}&-1.75\cdot 10^{-4}&6.63\cdot 10^{-7}\\[3pt] 1.87\cdot 10^{-6}&-8.45\cdot 10^{-7}&-2.68\cdot 10^{-6}&3.80\cdot 10^{-3}&-1.79\cdot 10^{-5}&6.47\cdot 10^{-5}\\[3pt] -9.50\cdot 10^{-7}&-4.65\cdot 10^{-7}&-1.75\cdot 10^{-4}&-1.79\cdot 10^{-5}&1.05\cdot 10^{-3}&3.45\cdot 10^{-5}\\[3pt] \mathbf{5.49\cdot 10^{-7}}&3.92\cdot 10^{-4}&6.63\cdot 10^{-7}&6.47\cdot 10^{-5}&3.45\cdot 10^{-5}&2.60\cdot 10^{-3}\end{array}} \right]\nonumber\\ \\[-3pt]\nonumber\end{align}

From Eq. (22), it is observed that the matrix obtained by the classical least-squares approach violates the conservative assumption: it is not symmetric and presents (upon Assertion 1 ) a negative eigenvalue; therefore, it is not positive-definite. It is important to emphasize that this result was obtained regardless of the full-rank observation matrix of this problem. On the other hand, Eq. (23) evinces that, for the same problem and dataset, the compliance matrix obtained by the GCO-based identification is conservative (it is symmetric and presents only positive eigenvalues).

5.2. Comparative illustrative example

Let us consider another example extracted from ref. [Reference Klimchik, Furet, Caro and Pashkevich36], which is a typical serial manipulator with six rotational actuated joints that presents elastic deformations on both joints and links. The adopted kinematic representation is depicted in Fig. 3. Each manipulator link has a regular hollow circular cross-section and all actuated joints have the same compliance coefficients that are equal to 10^–6 rad/N $\cdot$ m. The remaining parameters of this robot can be found in ref. [Reference Klimchik, Furet, Caro and Pashkevich36].

Figure 3. Kinematic representation of a typical 6 d.o.f. robot used in the illustrative example.

As in ref. [Reference Klimchik, Furet, Caro and Pashkevich36], position and orientation of the end-effector are assumed to be gathered with an accuracy of 25 $\mu$ m and 250 $\mu$ rad, respectively. In order to identify the elastostatic parameters of this robot, 18 measurement poses were randomly generated, each of which used for six calibration experiments with different external exciting forces/torques directed along the Cartesian axes with magnitudes set to 1000 N and 1000 N·m, driving to 108 virtual experimental test-runs.

In the frame of the GCO-based approach, it is enough to reduce the VJM model by some physical considerations only (‘Setup the VJM model’ in Fig. 2). It is worth pointing out that, from the intrinsic symmetry property associated with the Gramian constraint, the model’s parameters are inherently reduced to 153 (21 parameters for each of the 7 links and 6 joint compliances). Thus, considering just the physical aggregation of the joint compliances, the first model to be identified (MGCO147) presents 147 parameters. A second model presented in this analysis (MGCO69) departs from the assumption that, besides the joint compliances, just links 0 (base), 2 and 3 are subject to relevant compliance deviations (they are the longest links), reducing the number of parameters down to 69. This kind of assumption could be considered as reasonable as sparsing, which departs from some template assumption for the links’ compliance matrices, a similarly difficult condition to be verified in real. Nevertheless, for the sake of the completeness of this analysis, a third model (MGCO56) with 56 parameters is identified considering the full classical physical steps of aggregation and sparsing. For the weighted least-squares-based identification approach, two reduced models are considered: model MLS56, which was previously reduced according to the physical approach (symmetrization, aggregation, and sparsing) to 56 parameters, and model MLS32, which is further reduced according to the algebraic approach to the non-redundant model of the manipulator with 32 parameters [Reference Klimchik, Caro and Pashkevich35].

Table I summarizes all the identified models with their related characteristics. It is convenient to mention that, among all the models presented in Table I, only the GCO-based models led to physically meaningful models with conservative compliance matrices. The least-squares-based models did not achieve that, since this technique is expected to be sensible to the noise that overlays data in this virtual experiment.

Table I. Identification of elastostatic models for a 6 d.o.f. manipulator.

The evaluation of the validity of the GCO approach is carried out by the comparison of the identified models’ accuracy on the prediction of the compliance deviations experimented by the manipulator for 120 runs with different randomly generated poses under different loading conditions that were not used for identification. In this dataset, relevant end-effector deflections range from 9 to 41 mm. The position distance-based error (the Euclidean distance in mm) achieved by each model in each test-run is used as a metric to evaluate the models’ accuracy. To facilitate comparisons, these errors are also evaluated on a percentage basis of the end-effector deflection computed in each of these virtual test-runs.

Table II presents the results over the whole dataset, with the 95%-percentile (i.e. in 95% of all the 120 test-runs, the distance-based error of that model is lower than or equal to this value) and the RMS, on an absolute basis (mm), and with 95%-percentile and the maximum, on a percentage relative basis.

Table II. Numerical comparison of the GCO-based and least-squares-based models.

Table II shows that, besides of leading to physically meaningful models, the GCO-based approach achieved a quite acceptable accuracy. The best model was the MGCO56, with a maximum relative error of 0.62%, on a percentage basis, and a 95%-percentile of 0.072 mm, on an absolute basis. It is nearly the same as MLS32 (the best least-squares-based model), which is slightly worse with a 95%-percentile of 0.078 mm, on an absolute basis, and a maximum relative error of 0.80%. The GCO-based model MGCO69 presented a 95%-percentile of 0.084 mm, on an absolute basis, while the MGCO147 presented a maximum relative error of 0.87%, on a percentage basis. The model MLS56 performed worse than MLS32 (and the other GCO models), with a maximum relative error of 0.94%, on a percentage basis, and a 95%-percentile of 0.085 mm, on an absolute basis. It is interesting to note that, while the least-squares elastostatic identification is known to be prone to the measurement noise, the results of this section suggest that the noise did not significantly impair the accuracy of the GCO approach.