Identification of geometrical and elastostatic parameters  
of heavy industrial robots 
A. Klimchik, Y. Wu, C. Dumas, S. Caro, B. Furet and A. Pashkevich  
 
Abstract— The paper focuses on the stiffness modeling of 
heavy industrial robots with gravity compensators. The main 
attention is paid to the identification of geometrical and elas-
tostatic parameters and calibration accuracy. To reduce impact 
of the measurement errors, the set of manipulator configura-
tions for calibration experiments is optimized with respect to 
the proposed performance measure related to the end-effector 
position accuracy. Experimental results are presented that il-
lustrate the advantages of the developed technique.  
 
Keywords— industrial robot, stiffness modeling, elastostatic 
calibration, gravity compensator, design of experiments. 
I. INTRODUCTION 
Further developments in aeronautic industry require high-
precision and high-speed machining of huge aircraft compo-
nents where industrial robots successfully replace conven-
tional CNC-machines. For these applications, robots are quite 
attractive due to their large workspace (that can be easily ex-
tended) and high-speed motion capacity, as well as capability 
to process the parts with complex shape and geometry. How-
ever, because of high forces required for processing of mod-
ern aeronautic materials, the influence of robot stiffness be-
comes one of the key factors that essentially affect the posi-
tion accuracy. To enhance robot stiffness, manufacturers tend 
to increase the manipulator link cross-sections that obviously 
leads to the augmentation of the robot mass. So, the gravity 
forces applied to the manipulator components become non-
negligible and also contribute to the position errors. To over-
come this difficulty, the robots are usually equipped with dif-
ferent types of gravity compensators, which, however, con-
siderably complicate the stiffness modeling of those heavy 
manipulators.  
The problem of stiffness modeling for the heavy manipu-
lators with gravity compensators has been in the focus of 
rather limited number of works. In contrast, for conventional 
serial manipulators without gravity compensators, the prob-
lem has been studied by a number of authors that considered 
both industrial and medical robots with essential compliance 
in the links and joints. [1-3]. Relevant works are mainly 
based on the virtual joint method (VJM), which was firstly 
proposed by Salisbury [4] and lumped elastostatic properties 
of robot components in virtual springs. To our knowledge, 
the stiffness modeling for the manipulators with gravity com-
pensators has not been studied in detail yet. Currently, the 
main activity in this area focuses on the gravity compensator 
design, which differs in kinematics [5] and/or may also em-
ploy some software tools embedded in the robot controller 
[6]. On the other hand, since the considered robots include 
closed loops induced by the compensators, some technique 
developed for the parallel manipulators can be adopted[7-9].  
 
*The work presented in this paper was partially funded by the Agence 
Nationale de la Recherche (ANR), France (Project ANR-2010-SEGI-003-
02-COROUSSO).  
A. Klimchik, A. Pashkevich and Y. Wu are with Ecole des Mines de 
Nantes, 4 rue Alfred-Kastler, Nantes 44307, France and with Institut de Re-
cherches en Communications et Cybernétique de Nantes (IRCCyN), 44321 
Nantes, France (phone: Tel.+33-251-85-83-00; fax.+33-251-85-83-49; e-
mail: 
alexandr.klimchik@mines-nantes.fr, 
anatol.pashkevich@mines-
nantes.fr, yier.wu@mines-nantes.fr). 
S. Caro is with Centre National de la Recherche Scientifique (CNRS), 
France and with IRCCyN (e-mail: stephane.caro@irccyn.ec-nantes.fr). 
C. Dumas,  and B. Furet are with University of Nantes, France and with 
IRCCyN 
(e-mail: 
claire.dumas@univ-nantes.fr, 
benoit.furet@univ-
nantes.fr). 
This paper proposes a VJM-based stiffness model for a 
serial manipulator with a compensator attached to the second 
joint. It is assumed that the gravity compensation torque is 
generated by a spring incorporated in an additional link, 
which creates the closed loop to be included in the stiffness 
model. The main attention is paid to the identification of the 
model parameters and calibration experiment planning. The 
developed approach is confirmed by the experimental results 
that deal with compliance error compensation for robotic cell 
employed in manufacturing of large dimensional aircraft 
components. 
To address these problems the remainder of the paper is 
organized as follows. Section 2 presents the stiffness model-
ing background. In Section 3, the geometrical and elastostatic 
models of the gravity compensator are presented. Sections 4 
and 5 are devoted to the geometrical and elastostatic calibra-
tion. Section 6 summarizes the main contributions.  
II. STIFFNESS MODELING BACKGROUND 
Stiffness of a serial robot highly depends on its configura-
tion and is defined by the Cartesian stiffness matrix that 
should be computed in the neighborhood of the loaded equi-
librium configuration. Therefore, let us first present technique 
for computing the static equilibrium configuration and after 
focus on the computing of the stiffness matrix.  
Using the VJM-based approach adopted in this paper, the 
manipulator can be presented as the sequence of rigid links 
separated by the actuators and virtual flexible joints incorpo-
rating all elastostatic properties of compliant elements [10]. 
Since for the considered robots weights of the link is not neg-
ligible, it should be also incorporated in the model. In the 
frame of the VJM-based technique, it is convenient to replace 
the link weights 
i  by the equivalent pair of the forces ap-
plied to the both link ends 
i  and 
i . Further, after aggre-
gation of all weight components applied to the same virtual 
joint, the external loading caused by the link weights can be 
described by the set of forces 
 concentrated in the virtual 
P

P

P
i
G
 
 
joint centers. For computational convenience, they can be 
collected in the matrix 


1...
n

G
G
G
 that will be referred to 
as "the auxiliary loading". This allows us to distinguish im-
pact of the gravity from the influence of the conventional 
loading 
 applied to the robot end-effector (that can be 
caused by the cutting forces in milling, for instance). Using 
those notations and the above defined assumptions the elas-
tostatic model of the heavy serial robot can be presented as 
shown in Fig. 1. 
F
 
Figure 1.  Virtual joint model of a serial industrial robot 
For this serial chain, the end-effector location  and loca-
tions of the node centers 
t
jt
,
g
q θ
 can be defined by the vector 
functions 
 and 
j
, where 
 denote the vec-
tors of the actuator and virtual joint coordinates respectively. 
It can be proved that applying the principle of virtual work 
and linearizing the geometrical model, the static equilibrium 
equations can be written in a matrix form as [11] 
,
g q θ

θ
,
q θ



(G)T
(F)T
θ
θ





J
G
J
F
K
θ
where 
θ
θ
θ
 and 
θ
θ  denote the Jaco-
bians computed with respect to the relevant nodes where the 
forces are applied to, i.e. 
(G)
(1)
( )
T
T
T
...
n




J
J
J
(F)
( )
n

J
J


( )j

1
n
θ
θ
,...,
di
K
K
K
θ
j
; and the matrix 
 aggregates stiffnesses of all virtual 
springs.  
,


J
g
q θ
θ

θ
ag
1
1
i

θ

To find the desired static equilibrium configuration, the 
above system should be solved subject to the geometrical 
constraint 
, where the end-effector location is as-
sumed to be given and the force 
 is treated as an unknown. 
It is obvious that it is a dual problem compared to the tradi-
tional one where F  is given and the vector 
 is unknown, 
but this approach essentially reduces the computational ef-
fort. Relative iterative algorithm for solving this system can 
be presented as  
,

t
g q θ
i

F
t











1
(F)
1
(F)T
θ
θ
(F)
(F)
1
(G)T
θ
θ
θ
1
(G)T
(F)T
θ
θ
1
θ
1
,
·
i
i
i
i
i
i


















F
J
K
J
t
g q θ
J
θ
J
K
J
G
K
J
G
J
F
This algorithm allows us to find the vectors F , 
 defining 
deflections in the virtual springs under the loading F , which 
will be required for computing the stiffness matrix.  
θ
To find the corresponding stiffness matrix 
C
K , the force-
deflection relation (1) should be linearized in the neighbor-
hood of the equilibrium configuration 
. After rele-
vant transformations [11], one can get the following expres-
sion 
( ,
F θ, )t




1
(F)
1
(F)T
C
θ
θ
θθ
θ
(
)




K
J
K
H
J
which includes both the first and second order derivatives 
(Jacobians and Hessians) of the functions 


,
j
g
q θ  describ-
ing 
the 
manipulator 
geometry: 


(
)
θ
F
,


q θ
θ
J
g
, 
(F)
(G)
(G)T
θθ
θθ
θθ
θ





H
H
H
J
G
θ , 
(G)
2
θθ
1
2
T
j
j
j
n





H
g G
θ , 
(F
θθ
2
)
2
T


H
g F
θ . 
It should be noted that the above presented expressions 
have been derived for the case of manipulators without grav-
ity compensators, where the matrix 
θ  is constant. In the 
following Section, to take into account the compensator in-
fluence, these expressions will be modified by using the con-
figuration dependent joint stiffness matrix 
 describing 
properties of the virtual springs. 
K
θ( )
K
q
III. MECHANICS OF GRAVITY COMPENSATOR 
The mechanical structure of the gravity compensator un-
der study is presented in Fig. 2. The compensator incorpo-
rates a passive spring attached to the first and second links, 
which creates a closed loop that generates the torque applied 
to the second joint of the manipulator. This design allows us 
to limit the stiffness model modification by incorporating in it 
the compensator torque 
c
M  and adjusting the virtual joint 
stiffness matrix 
 that here depends on the second joint 
variable 
 only. 
θ
K
2q
xa
y
a
2q



s
L
0P
0P
2P
1P
01
P
02
P
2P
1P
a
ck
2
qk
sF
cF
lF
O
 
Figure 2.  Gravity compensator and its model 
The compensator geometrical model includes three node 
points P0, P1, P2, where two distances 
1
2
,
P P , 
0
2
,
P P  are 
constants and the third one 
0
1
,
P P  varies and depends on 
2 . 
Let us denote them 
q
1
2 , 
,
L
P P

1
2 , 
,
a
P P

1
2 . Be-
sides, let us introduce the angles 
,
s
P P

,  and the distances 
xa  
and 
y , whose geometrical meaning is clear from Fig. 2. Us-
ing these notations, the variable 
a
s  describing the compensa-
tor spring deflection can be computed from the expression  

2
2
2
2
· · ·co
2
s(
a
L
)
s
a L
q







 
 
which defines the function 
.  


2
s q
This mechanical design allows to balance the manipulator 
weight for any given configuration by adjusting the compen-
sator spring preloading. It can be taken into account by intro-
ducing the zero-value of the compensator length 
0s  corre-
sponding to the unloaded spring. Under this assumption, the 
compensator force applied to the node P1 can be expressed as 
follows 

0
·(
)
s
c
s
F
K
s




where 
 is the compensator spring compliance.  
ck
Further, the angle  between the compensator links P0P1 
and P1P2 (see Fig. 2) can be found from the expression 

2
sin
s·
(
ni
q
a s
)






which allows us to compute the compensator torque 
c
M  ap-
plied to the second joint 



0
·(1
)· · ·
(
/
sin
c
c
s
s
L
M
a
K




2)
q
2
q
Upon differentiation of the latter expression with respect to 
2 , the equivalent stiffness of the second joint (comprising 
both the manipulator and compensator stiffnesses) can be ex-
pressed as: 
q

2
2
0
·
·
c
K
K
K
L
a







where the coefficient  

2
2
0
2
2
2 sin
cos
co
·
(
s
)
(
)
(
q
s
a L
q
q
s
s















2)
q


highly depends on the value of joint variable 
2  and the ini-
tial preloading in the compensator spring described by 
q
0s . To 
illustrate this property, Fig. 3 presents a set of curves 
2
(
)
q

 
obtained for different values of 
0s  (the remaining parameters 
, 
, 
 correspond to robot KUKA KR-270 studied in the 
experimental part of the paper, see Section IV). 
a
L
-150
-120
-90
-60
-30
0
-0.05
0
0.05
0.1
0.15
0
0 [m]
s 
0
0.2 [m]
s 
0
0.2 [m]
s 
0
0.4 [m]
s 
0
0.6 [m]
s 
0
0.8 [m]
s 
2,[deg]
q

 
Figure 3.  Variation of the gravity compensator impact on the equivalent 
stiffness of the second joint  
Hence, using expression (8), it is possible to extend the 
classical stiffness model of the serial manipulator presented 
in Section II by modifying the virtual spring parameters in 
accordance with the compensator properties. While in the pa-
per this approach has been used for the particular compensa-
tor type (spring-based, acting on the second joint), the similar 
idea can be evidently applied to other compensator types. 
Summarizing this Section, it is worth mentioning that the 
geometrical and elastostatic models of a heavy manipulator 
with a gravity compensator should include some additional 
parameters (, 
, 
 and 
c
a
L
K , 
0s  for the presented case) 
that are usually not included in datasheets. For this reason, 
the following Sections focus on the identification of the ex-
tended set of manipulator parameters.  
IV. IDENTIFICATION OF COMPENSATOR GEOMETRY  
In contrast to the serial manipulator that can be treated as 
a principal mechanism of the considered robots (whose ge-
ometry is usually defined in datasheets and can be perfectly 
tuned by means of calibration [12][13]), geometrical data 
concerning gravity compensators are usually not included in 
the technical documentation provided by the robot manufac-
turers. For this reason, this Section focuses on the identifica-
tion of the geometrical parameters for the described above 
compensator mechanism (see Fig. 2). 
A.  Methodology 
The geometrical structure of the considered gravity com-
pensator is presented in Fig. 4. Its principal geometrical pa-
rameters are denoted as 
, 
L
xa , 
y , where 
x
a
·cos
a
a


, 
·sin
y
a
a


 (see notation in Section III). As follows from 
the figure, the identification problem can be reduced to the 
determination of relative locations of points P0 and P1 with 
respect to P2.  
It is assumed that the measurement data are provided by 
the laser tracker whose "world" coordinate system is located 
at the intersection of the first and second actuated 
manipulator joints. The axes Y, Z of this system are aligned 
with the axes of joints #1 and #2 respectively, while the axis 
X is directed to ensure right-handed orthogonal basis. To 
obtain required data, there are several markers attached to the 
compensator mechanism (see Fig. 4). The first one is located 
at point P1, which is easily accessible and perfectly visible 
(the center of the compensator axis P1 is exactly ticked on the 
fixing element). In contrast, for the point P0, it is not possible 
to locate the marker precisely. For this reason, several 
markers P0i are used that are shifted with respect to P0, but 
located on the rigid component of the compensator 
mechanism (these markers are rotating around P0 while the 
joint coordinate 
2  is actuated). It should be noted that for 
the adopted compensator geometrical model (which is in fact 
a planar one), the marker location relative to the plane XY is 
not significant, since the identification algorithm presented in 
the following sub-section will ignore Z-coordinate.  
q
Using this setting, the identification problem is solved in 
two steps. The first step is devoted to the identification of the 
relative location of points P1 and P2. Here, for different 
values of the manipulator joint coordinates 
, 
2
{
i
q
1, }
i
m

, 
the laser tracker provides the set of the vectors 
1
ip
 describ-
ing the points that are located in an arc of the circle. After 
matching these points with a circle, one can obtain the de-
sired value of 
 (circle radius) and the Cartesian coordinates 
2  of the point P2 (circle center) with respect to the laser 
tracker coordinate system. 
L
p
 
 
xa
y
a
2q

L
0P
2P
1P
a
O
02
P
01
P
1
R
2
R
2

1


,
,
x
y
L a
a
03
P
04
P
3
R
4
R
 
Figure 4.  Geometrical parameters of the gravity compensator and location 
of the measurement points  labeled with markers  
The second step deals with the identification of the rela-
tive location of points P0 and P2. Relevant information is ex-
tracted from two data sets 
01  and 
02  that are provided 
by the laser tracker while targeting at the markers P01 and P02. 
Here, the points are matched to two circle arcs with the same 
center (explicitly assuming that the compensator model is 
planar), which yields the Cartesian coordinates 
0  of the 
point P0, also with respect to the laser tracker coordinate sys-
tem. Finally, the desired values 

ip

ip
p
xa , 
 are computed as a 
projection of the difference vector 
 on the corre-
sponding axis of the coordinate system. 
y
a
2


a
p
0
p
As follows from the presented methodology, a key nu-
merical problem in the presented approach is the matching of 
the experimental points with a circle arc. It looks like a clas-
sical problem, however, there is a particularity here caused by 
availability of additional data 
2i  describing relative loca-
tions of the points 
. This feature allows us to reformu-
late the identification problem and to achieve higher accuracy 
compared with the traditional approach.  

q

1
ip
B.  Identification algorithm 
The above presented methodology requires solution of 
two identification problems. The first one aims at approxi-
mating of a given set of points (with additional arc angle ar-
gument) with an arc circle, which provides the circle center 
and the circle radius. The second problem deals with an ap-
proximation of several sets of points by corresponding num-
ber of circle arcs with the same center. Let us consider them 
sequentially. 
To match the given set of points 
 with additional set 
of angles 

ip

iq
 with a circle arc, let us define the affine map-
ping 

i
i



p
Ru
t 

where 
i
i
i
 denotes the set of reference 
points located on the unit circle whose distribution on the arc 
is similar to 
i , 
[cos
, sin
, 0]T
q
q

u
p
 is the scaling factor that defines the de-
sired circle radius, R  is the orthogonal rotation matrix, t  is 
the vector of the translation that defines the circle center. It 
worth mentioning that such a formulation has an advantage 
(in the sense of accuracy) comparing to a traditional circle 
approximation and it is a generalization of Procrustes prob-
lem known from the matrix analysis. 
Using equation (10), the identification can be reduced to 
the following optimization problem 





1
,
,
min
T
i
i
i
i
i
m
F











R t
p
Ru
t
p
Ru
t
which should be solved subject to the orthogonality con-
straint 
T

R R
I . After differentiation with respect to t , the 
latter variable can be expressed as  



1
1
1
1
i
m
i
m
m








t
p
R
i
m
iu
That leads to the simplification of  (11) to 




1
,min
m
i
r
T
i
i
i
i
F








R
p
Ru
p
Ru






where  

1
1
1
1
;
i
i
i
i
m
i
i
i
i
m
m
m










p
p
p
u
u
u



Further, differentiation with respect to  yields to 

1
1
m
T
T
i
i
i
i
m
i





p Ru
u ui




So, finally, after relevant substitutions the objective func-
tion can be presented as   


1
1
1
1
2
min
m
m
m
i
i
i
T
T
T
i
i
i
i
i
i
F






















R
p p
u u
p Ru




where the unknown matrix 
 must satisfy the orthogonality 
constraint 
R
T

R R
I . Since the matrix R  is included in the 
second term only, the problem can be further simplified to  



1
1
max
m
m
T
T
i
i
i
i
i
i
F
trace







R
p Ru
R
u p






and can be solved using SVD-decomposition of the matrix  

1
m
i
i
T
i



u p
U ΣV




where the matrices 
 are orthogonal and 
 is the diago-
nal matrix of the singular values. Further, using the same ap-
proach as for the Procrustes problem, it can be proved that 
the desired rotation matrix can be computed as 
,
U V
Σ



T

R
V U
which sequentially allows to find the scaling factor  defin-
ing the arc radius and the vector  defining the arc center.  
t
The second problem aims at approximating of several 
point sets 

1 ,...,
k
p
p
i
i
 by corresponding number of concen-
tric circle arcs with the same center 
0
p . It should be noted 
that here the data set 
i  is not useful, since the required 
angles 
q

i
 are not measured directly and cannot be com-
puted without having exact compensator geometry. In this 
case, the objective function can be written in a straightfor-
ward way 






0
2
0
2
1
1
,
0
min
j
k
m
j
j
j
i
i
j
i
R
T
F
R







p
p
p
p
p


 
 
But it can be proved that this optimization problem does not 
lead to a unique solution (in fact, it gives the rotation axis 
passing through the desired center). For this problem, differ-
entiation of the objective function F  with respect to 
2
j
R  
yields  






0
2
1
1
0
T
m
j
j
j
i
i
i
R
m




p
p
p
p

min

which after substitution into (20) allows us to rewrite the 
problem in the following way  



0
1
1
2
0
2
k
m
j
j
i
i
j
i
T
F
s





p
p p



where  
1
1
;
m
1
1
m
j
j
j
j
j
j
i
i
l
i
i
i
i
i
l
T
m
s
m







p
p
p
p p
p p
j
j
l
T







Further, after differentiation with respect to 
0
p , one can get 
the underdetermined system of linear equations 



1
1
0
1
1
1
2
k
m
k
m
j
j
j
i
i
i
i
j
i
j
i
T
j
s







p p
p
p





whose solution  

0
c



p
p
n 

is expressed via the rotation axis vector 
 and a point be-
longing to this axis 
 (here 
n
c
p
 is an arbitrary scalar factor). 
To solve this ambiguity, an additional objective should be 
defined  



2
1
min
j
k
j
j
R
R



which leads to the following solution for the scalar parameter 


1
1
1
1
k
m
c

j
i
j
i
T
k m









n
p
p


and for the vector  



1
1
1
1
1
2
k
m
k
m
1
j
j
c
i
i
j
i
j
i
T
j
j
i
i
s








p
p p
p





So, the desired arc center is expressed as follows 



1
1
0
1
1
T
T
c
k
m
j
i
j
i
k m








p
I
nn
p
nn
p 

It should be mentioned that practical application of the latter 
expression is essentially simplified by the adopted assump-
tion concerning orientation of the reference coordinate sys-
tem (see previous sub-section), the direction of the identified 
rotation axis is close to Z-direction.  
Hence, the developed algorithms allows us to identify the 
compensator geometrical parameters 
, 
L
xa , 
y  that are di-
rectly related to the above mentioned rotation center points 
P0, P2 and corresponding radii. Below they will be applied to 
the processing of the experimental data.  
a
C.  Experimental results 
To demonstrate efficiency of the developed technique, the 
experimental study has been carried out. The experimental 
setup employed the robot KR-270 and the Leica laser tracker, 
which allowed us to measure the Cartesian coordinates of the 
markers attached to the compensator elements (see Figs 2, 4). 
Six different manipulator configurations where considered 
that differed in the value of the joint angle 
2  and three 
markers has been used. The experimental data are presented 
in Table I.  
q
These data has been processed using the identification al-
gorithm presenting in the previous sub-section. The obtained 
values for the parameters of interest 
, 
L
xa , 
y  are given in 
Table II. It also includes the confidence intervals computed 
as 
a
3

, where the standard deviation  has been evaluated 
using the Gibbs sampling. In the next section, the obtained 
model will be used for some transformations required during 
elastostatic calibration 
TABLE I.  
EXPERIMENTAL DATA FOR GEOMETRICAL CALIBRATION 
P1 
P01 
P02 
q2 
[deg] 
x, [mm] 
y, [mm] 
x, [mm] 
y, [mm] 
x, [mm] 
y, [mm] 
-0.01 
-31.84 
183.86 
-872.10 
-125.38 
-813.50 
-255.59 
-30 
-118.44 
143.42 
-872.30 
-126.07 
-813.33 
-256.18 
-60 
-173.30 
65.12 
-872.50 
-109.90 
-825.09 
-244.64 
-90 
-181.76 
-30.14 
-868.43 
-78.20 
-844.66 
-219.04 
-120 
-141.45 
-116.82 
-858.90 
-47.60 
-859.43 
-190.44 
-145 
-78.10 
-165.47 
-852.53 
-33.68 
-864.66 
-176.01 
TABLE II.  
GEOMETRICAL PARAMETERS OF GRAVITY COMPENSATOR 
 
L, [mm] 
ax, [mm] 
ay, [mm] 
value 
184.72 
685.93 
120.30 
CI 
±0.06 
±0.70 
±0.69 
 
V. ELASTOSTATIC PARAMETERS IDENTIFICATION   
In contrast to geometrical calibration, where the manipu-
lator and compensator can be considered independently, in 
elastostatic calibration the corresponding equations cannot be 
separated and the model parameters should be identified si-
multaneously. This section gives general ideas of the devel-
oped methodology and relevant identification algorithms, and 
also presents experimental results validating the proposed 
technique. 
A.  Methodology 
In the frame of the adopted VJM-based modeling ap-
proach the desired stiffness parameters describe elasticity of 
the virtual springs located in the actuated joints of the ma-
nipulator, and also elasticity and preloading of the compensa-
tor spring (see Section III for details). Let us denote them as 
,
1,6
k
j 
j

 for the manipulator joint compliances and 
0
,
ck
s  
for the compliance and preloading of the compensator. 
To find the desired parameters, the manipulator sequen-
tially passes through several measurement configurations 
where the external loading is applied to the special end-
effector described in Fig. 5 (it allows to generate both forces 
and torques applied to the manipulator). Using the laser 
tracker, the Cartesian coordinates of the reference points are 
 
 
measured twice, before and after loading. To increase identi-
fication accuracy, it is preferable to use several reference 
points (markers) and to apply the loading of the maximum al-
lowed magnitude. It is worth mentioning that in order to 
avoid numerical singularities, the direction of the external 
loading should not be the same for all experiments (in spite of 
the fact that the gravity-based loading is the most attractive 
from the practical point of view). Thus, the calibration ex-
periments yield the dataset that includes values of the ma-
nipulator joint coordinates 
i , applied forces/torques 

q

iF
 
and corresponding deflections of the reference points 

i
p
. 
Using these data, the elastostatic parameters of 
,
1,6
j
k
j 
 
and 
0
,
ck
s  should be identified. 
x
z
460 mm
75.5 mm
Joint #6
Tool
F
Markers
 
Figure 5.  End-effector used for elastostatic calibration experiments 
B.  Identification algorithm 
To take into account the compensator influence while re-
taining our previous approach developed for serial robots 
without compensators, it is proposed to include in the second 
joint an equivalent virtual spring with non-linear stiffness de-
pending on the joint variable 
2  (see eq. (8)). Using this idea, 
it is convenient to consider several independent parameters 
2i
q
k corresponding to each value of 
2q . This allows us to ob-
tain linear form of the identification equations that can be 
easily solved using standard least-square technique. 
Let 
us 
denote 
the 
set 
of 
desired 
parameters 
1
21
22
3
6  as the vector k  and linearize (1) with 
respect to this vector. This allows us to present the relevant 
force displacement relations in the form  
,(
,
...),
,...,
k
k
k
k
k



(
)
p
i
i


p
B
k
where matrices 
 are composed of the elements of the ma-
trix  
(
)
p
i
B

1
1
,...
1,
)
,
(
T
T
i
i
i
i
ni
ni
i
i
m




F
A
F
J J
J J



that is usually used in stiffness analysis of serial manipulators 
[14]. Here, 
ni  denotes the manipulator Jacobian column, 
iF  
is the applied external force, and superscript '(p)' stands for 
the Cartesian coordinates (position without orientation). It is 
clear that transformation from 
i  to 
J
A
(
)
p
i
B
 is rather trivial and 
is based on the extraction from 
i  the first three lines and 
inserting in it several zero columns.  
A
Using these notations, the elastostatic parameters identifi-
cation can be reduced to the following optimization problem 

0
,
,
( )
( )
1(
) (
)
j
c
m
p
T
p
i
i
i
i
k
k
i
F
which leads to the following solution  



1
(
)
(
)
( )
1
1
·
T
m
m
p
p
p
i
i
i
i
i
i
















k
B
B
B
p
T


where the parameters 
1
3
6  describe the compliance of 
the virtual joints #1,#3,...#6, while the rest of them 
21
22
 
present an auxiliary dataset allowing to separate the compli-
ance of the joint #2 and the compensator parameters 
c
,
,...,
k k
k
,
...
k
k
0
,
k
. 
Using eq. (8), the desired expressions can be written as  




2
2
1
0
0
1
1
·
q
q
i
T
c
c
i
m
m
T
T
i
i
i
i
K
K
K
K
s











C C
C

where 
q  is the number of different angles 
 in the ex-
perimental data, 
m
2q


2
2
1
·
/ ·
/
·
cos
sin
cos
i
i
i
L
L
a
s
a
L
a
s
i










C

here 
2
i
i
q




. 
Thus, the proposed modification of the previously devel-
oped calibration technique allows us to find the manipulator 
and compensator parameters simultaneously. An open ques-
tion, however, is how to find the set of measurement configu-
rations that ensure the lowest impact of the measurement 
noise.  
C.  Design of calibration experiments 
The main idea of the calibration experiment design is to 
select a set of robot configurations 
i  (and corresponding 
external loadings 

q

i ) that ensure the best identification ac-
curacy. The key issue here is the ranging of different plans in 
accordance with the prescribed performance measure. This 
problem has been already studied in the classical regression 
analysis [15], however the results are not suitable for the elas-
tostatic calibration and require additional efforts.  
F
In this work, it is proposed to use the industry oriented 
performance measure that evaluates the calibration plan qual-
ity. Its physical meaning is the robot positioning accuracy 
(under the loading), which is achieved after compliance error 
compensation based on the identified elastostatic parameters. 
It should be noted that usual approach (used in the classical 
design of experiments) evaluates the quality of experiments 
via covariance matrix of the identified parameters, which is 
does not have sense for our application.   
Assuming that each experiment includes the additive 
measurement error 
i , the covariance matrix for the desired 
parameters  
 can be expressed as 
ε
k





1
( )
( )
1
1
( )
( )
( )
( )
1
1
cov( )
E
T
T
T
m
p
p
i
i
i
m
m
p
T
p
p
p
i
i
i
i
i
i
i
i











k
B
B
B
ε ε B
B
B


min







B
k
p
B
k
p


Following also usual assumption concerning the measure-
ment errors (independent identically distributed, with zero 
expectation and standard deviation 
2
 for each coordinate), 
the above equation can be simplified to 
 
 




1
2
( )
( )
1
cov( )
T
m
p
p
i
i
i





k
B
B 
Hence, the impact of the measurement errors on the accuracy 
of the identified parameters k  is defined by the matrix 
1
( )
( )
T
m
p
p
i
i

B
B
i
 (in regression analysis it is known as the in-
formation matrix). 
It is evident that in practice the most essential is not the 
accuracy of the parameters identification, but the accuracy of 
the robot positioning achieved using these parameters. Tak-
ing into account that this accuracy highly depends on the ro-
bot configuration (and varies throughout the workspace), it is 
proposed to evaluate the calibration accuracy in a certain 
given "test-pose" provided by the user. For the considered 
application, the test pose is related to the typical machining 
configuration 
0  and corresponding external loading 
0  re-
lated to the technological process. Let us denote the mean 
square value of the mentioned positioning error as 
0
q
F
2
 and 
the matrix 
( )
p
i
A
 (see eq. (31)) corresponding to this test pose 
as 
. 
( )
0
p
A
It should be noted that that the proposed approach oper-
ates with a specific structure of the parameters included in the 
vector 
, where the second joint is presented by several 
components 
21
22
 while the other joints are described by 
a single parameter 
1
3
6 . This motivates further re-
arrangement of the vector k  and replacing it by several vec-
tors 
1
2
3
6
j
j
 of size 
. Using this notation, 
the above mentioned performance measure can be expressed 
as 
k

,
..
k
k
( ,
,
,.
k k
k
.

( )
,
...
k k
k
..
)
k
k
6 1



2
( )
0
0
0
1E
T
q
T
p
p
j
j
m
j




k
A
A
k



where 
j
k  is the elastostatic parameters estimation error 
caused by the measurement noise for 
2 j
q
T
. Further, after sub-
stituting 
 and taking into account that 

trace
T




p
p
p p
cov(
)

E(
)
T
j
j
j , the performance measure 


k
k

k
2
0
 can 
be presented as 



1
( )
( )
2
0
2
)
0
1
1
(
)
trace
T
q
T
m
j p
j p
m
p
i
i
i
j






0
(
p





A
A
A
A

Based on this performance measure, the calibration ex-
periment design can be reduced to the following optimization 
problem 




1
( )
( )
0
0
1
{
,
}
( )
( )
1
trace
min
T
q
T
i
i
m
j p
j p
i
p
p
j
i
m
i











q F
A
A
A
A
subject to 

max,
1.
i
.
F
i
m


F


whose solution gives a set of the desired manipulator con-
figurations and corresponding external loadings. It is evident 
that its analytical solution can hardly be obtained and a nu-
merical approach is the only reasonable one. More detailed 
description of the developed technique providing generation 
of the optimal calibration plans is presented in our previous 
publication [14], where the problem of gravity compensator 
modeling had not been addressed yet.  
D.  Experimental results 
Similar to the previous section, the developed technique 
has been applied to the robot KR-270. The parameters of in-
terest were the compliances 
jk  of the actuated joints and the 
elastostatic parameters 
0
,
ck
s  of the gravity compensator. To 
generate elastostatic deflections, the gravity forces 250-280 
kg have been applied to the robot end-effector (see Fig 6). 
The Cartesian coordinates of three markers located on the 
tool (see Fig, 5) have been measured before and after the 
loading.  
Gravity compensator
Tool
F
Force censor
Test loading
 
Figure 6.  Experimental setup for the identification  
of the elastostatic parameters 
To find optimal measurement configurations, the design 
of experiments has been carried out for five different angles 
2  that are distributed between the joint limits. For each 
2  
three optimal measurement configurations have been found 
taking into account physical constraints that are related to the 
joint limits and the possibility to apply the gravity force 
(work-cell obstacles and safety reasons). The results of the 
calibration experiment design are presented in Table III. 
Here, values of 
1  have been chosen to ensure good visibility 
of the markers for the laser tracker. To ensure identification 
accuracy for each configuration, the experiments were re-
peated three times. In total, 405 equations were considered  
for the identification, from which 7 physical parameters have 
been obtained (because of page limit, the experimental data 
cannot be presented in the paper ).  
q
q
q
TABLE III.  
OPTIMAL MEASUREMENT CONFIGURATIONS 
Joint angles, [deg] 
q1 
q2 
q3 
q4 
q5 
q6 
79.20 
-5.57 
51.00 
-97.52 
-91.67  
63.00 
-12.22 
-56.49 
41.42 
150.55 
63.00 
-0.01 
-47.98 
-70.04 
-61.55 
177.16  
95.00 
33.00 
129.69 
-98.10 
90.57 
95.00 
-107.01 
109.95 
-61.19 
174.21 
105.00 
-25.24 
14.30 
55.21 
41.26 
-152.97 
56.60 
44.54 
-55.11 
41.90 
152.06 
56.60 
64.73 
-129.65 
-98.260 
-90.55 
144.80 
-56.9 
104.49 
-69.41 
61.67 
-6.33  
-41.00 
-91.68 
55.12 
41.53 
-152.48 
-143.00 
-32.64 
110.31 
-61.47 
-6.29 
-143.00 
-99.85 
-72.01 
129.65 
-98.09 
90.82 
133.00 
147.68 
129.64 
-97.90 
90.99 
-60.00 
7.59 
-110.09 
-61.36 
-174.09 
-60.00 
-140 
-52.00 
-124.89 
-41.62 
27.78 
 
 
TABLE IV.  
ELASTO-STATIC PARAMETERS OF ROBOT KUKA KR-270  
Parameter  
value 
CI 
kc, [rad×μm/N] 
0.144 
±0.031 (21.5%) 
s0, [mm] 
458 
±27 (5.9%) 
k2, [rad×μm/N] 
0.302 
±0.004 (1.3%) 
k3, [rad×μm/N] 
0.406 
±0.008 (2.0%) 
k4, [rad×μm/N] 
3.002 
±0.115 (3.8%) 
k5, [rad×μm/N] 
3.303 
±0.162 (4.9%) 
k6, [rad×μm/N] 
2.365 
±0.095 (4.0%) 
 
The obtained experimental data have been processed us-
ing the identification algorithm presented in sub-section V.B. 
Corresponding values of the gravity compensator and the 
manipulator elastostatic parameters are presented in Table IV. 
It also includes the confidence intervals computed as 
3

, 
where the standard deviation  has been evaluated using 
Gibbs sampling. Using the obtained results it is possible to 
identify an equivalent non-linear spring 
2
2  (see Fig, 7) 
that is used in the stiffness modeling of the manipulator with 
the gravity compensator (see Section II). 
(
)
k
q
The identified joint compliances can be used to predict 
robot deformations under the external loading. The identifica-
tion errors of the joint compliances vary from 1.3% to 4.9%, 
where the lowest errors have been achieved for the joints #2 
and #3. The highest errors are in the joints #4-#6. So, higher 
precision has been achieved for the joints that are closer lo-
cated to the manipulator base. This fact is due to the different 
joint errors contribution to the total positioning error, which 
have 
been 
minimized 
while 
design 
of 
calibration 
experiments. Comparatively low identification accuracy for 
the compensator spring is caused by limited number of 
different angles 
2q . Fig. 7 shows that due to the gravity 
compensator the equivalent compliance of the second joint 
has been decreased comparing to the compliance of the 
second joint of the corresponding serial manipulator.  
-150
-120
-90
-60
-30
0
2.75
2.8
2.85
2.9
2.95
3
3.05 x 10
-7
2k
0
2k
2q
rad·m/N
,[deg]
[11] A. Klimchik, D. Bondarenko, A. Pashkevich, S. Briot, B. Furet, 
"Compensation of tool deflection in robotic-based Milling," the 9th In-
ternational Conference on Informatics in Control, Automation and 
Robotics (ICINCO 2012), July 28-31, 2012, Rome, Italy, pp. 113-122. 
 
Figure 7.  Compliance of equivalent non-linear spring  in the second joint 
VI. CONCLUSIONS 
The paper presents a new approach for the identification 
of the elastostatic parameters of heavy industrial robots with 
the gravity compensator. It proposes a methodology and data 
processing algorithms for the identification of both geometri-
cal and elastostatic parameters of gravity compensator and 
manipulator. To increase the identification accuracy, the de-
sign of experiments has been used aimed at proper selection 
of the measurement configurations. In contrast to other 
works, it is based on a new industry oriented performance 
measure that is related to the robot accuracy under the load-
ing. 
The advantages of the developed techniques are illus-
trated by experimental study of the industrial robot Kuka KR-
270, for which the joint compliances and parameters of the 
gravity compensator have been identified. 
ACKNOWLEDGMENT 
The work presented in this paper was partially funded by 
the 
ANR, 
France 
(Project 
ANR-2010-SEGI-003-02-
COROUSSO). The authors also thank Fabien Truchet, Guil-
laume Gallot, Joachim Marais and Sébastien Garnier for their 
great help with the experiments.  
REFERENCES 
[1] 
M.A. Meggiolaro, S. Dubowsky, C. Mavroidis, "Geometric and elastic 
error calibration of a high accuracy patient positioning system," 
Mechanism and Machine Theory, vol. 40, 415–427, 2005. 
[2] 
J. Kövecses, J. Angeles, "The stiffness matrix in elastically articulated 
rigid-body systems," Multibody System Dynamics (2007) pp. 169–
184, 2007. 
[3] 
B.-J. Yi, R.A. Freeman, "Geometric analysis antagonistic stiffness re-
dundantly actuated parallel mechanism", Journal of Robotic Systems 
vol. 10(5), pp. 581-603, 1993 
[4] 
J. Salisbury, “Active stiffness control of a manipulator in Cartesian 
coordinates,” in Proc. 19th IEEE Conf. Decision Control, 1980, pp. 
87–97. 
[5] 
N. Takesue, T. Ikematsu, H. Murayama and H. Fujimoto, "Design and 
Prototype of Variable Gravity Compensation Mechanism (VGCM)," 
Journal of Robotics and Mechatronics, Vol. 23, No. 2, pp. 249-257, 
2011. 
[6] 
A. De Luca, F. Flacco, ''A PD-type Regulator with Exact Gravity Can-
cellation for Robots with Flexible Joints," in Proc. 2011 IEEE Interna-
tional Conference on Robotics and Automation, pp. 317-323 
[7] 
O. Company, F. Pierrot, J.-C. Fauroux, "A Method for Modeling Ana-
lytical Stiffness of a Lower Mobility Parallel Manipulator," in: Pro-
ceedings of IEEE International Conference on Robotics and Automa-
tion (ICRA), 2005, pp. 3232 - 3237 
[8] 
J.-P. Merlet, C. Gosselin, Parallel mechanisms and robots, In B. Sicili-
ano, O. Khatib, (Eds.), Handbook of robotics, Springer, Berlin, 2008, 
pp. 269-285. 
[9] 
B.C. Bouzgarrou, J.C. Fauroux, G. Gogu, and Y. Heerah, "Rigidity 
analysis of T3R1 parallel robot with uncoupled kinematics," Proc. of 
the 35th International Symposium on Robotics, Paris, France, 2004. 
[10] A. Pashkevich, A. Klimchik, D. Chablat, "Enhanced stiffness model-
ing of manipulators with passive joints," Mechanism and Machine 
Theory 46(2011) 662-679. 
[12] D. Daney, N. Andreff, G. Chabert, Y. Papegay, "Interval method for 
calibration of parallel robots: Vision-based experiments," Mechanism 
and Machine Theory, vol. 41, 2006, pp. 929-944. 
[13] Hollerbach J., Khalil W., Gautier M., Springer Handbook of robotics, 
Springer, 2008, "Chapter: Model identification," pp. 321-344. 
[14] A. Klimchik, Y. Wu, A. Pashkevich, S. Caro, B. Furet, "Optimal Se-
lection of Measurement Configurations for Stiffness Model Calibra-
tion of Anthropomorphic Manipulators," Applied Mechanics and Ma-
terials, Volume 162, 2012, pp. 161-170 
[15] Atkinson A., DoneA., Optimum Experi-ment Designs. Oxford Univer-
sity Press, 1992