System Identification and Control of Valkyrie through SVA–Based
Regressor Computation

arXiv:1608.02683v3 [math.OC] 12 Sep 2016

Shishir Kolathaya1 , Benjamin J. Morris1 , Ryan W. Sinnet1 and Aaron D. Ames2
Abstract— This paper demonstrates simultaneous identification and control of the humanoid robot, Valkyrie, utilizing Spatial Vector Algebra (SVA). In particular, the inertia, Corioliscentrifugal and gravity terms for the dynamics of a robot are
computed using spatial inertia tensors. With the assumption
that the link lengths or the distance between the joint axes
are accurately known, it will be shown that inertial properties
of a robot can be directly evaluated from the inertia tensor.
An algorithm is proposed to evaluate the regressor, yielding
a run time of O(n2 ). The efficiency of this algorithm yields
a means for online system identification via the SVA–based
regressor and, as a byproduct, a method for accurate modelbased control. Experimental validation of the proposed method
is provided through its implementation in three case studies:
offline identification of a double pendulum and a 4-DOF robotic
leg, and online identification and control of a 4-DOF robotic
arm.

I. I NTRODUCTION
The field of system identification has been a subject of
much attention in the field of robotic systems, since the
1980’s [8], [2], [17], [12]. One of the most important reasons
is that achieving good tracking performance in robotic systems, specifically exponential convergence, without knowing
the complete model of the robot has not been shown.
Asymptotic convergence in tracking without knowing the
model parameters has been shown in [16], [4] by using a
special property of the Lagrangian dynamics of the robot,
that is, the linearity of the parameters in the dynamics.
Exponential convergence specific to a particular task using
machine learning was also shown in [5] by using the concept
of persistence of excitation [11], which requires running a
series of trials for the controller to learn.
As a deviation from the methods shown above, it can be
argued that identification of physical systems is required to
realize good tracking performance. Currently, many system
identification procedures have been implemented mainly by
either measuring the parameters of the robot part by part
[7], or dynamically by using the acceleration, velocity and
angles of the robot [12], [14]. [1] showed that the parameters
that are identified as linear combinations can be consistently
set to zero to determine the set of identifiable parameters.
Specifically, since the parameters are affine, the equation of
motion can be expressed as a matrix (regressor) multiplied by
*This research is supported by NASA grant NNX11AN06H, NSF grants
CNS-0953823 and CNS-1136104, and NHARP award 00512-0184-2009.
1 R. Sinnet and B. Morris are with the department of Mechanical
Engineering, Texas A & M University, 3128 TAMU, College Station, USA

shishirny,rsinnet,bmorris@tamu.edu
2 S. Kolathaya and Prof. A. D. Ames are with department of Mechanical
Engineering, Georgia Institute of Technology, 85 5th St, Atlanta, USA

shishirny,ames@gatech.edu

Fig. 1: Figures showing the 4-DOF robotic arm on the left
and the 4-DOF leg on the right of the Valkyrie robot on
which identification was conducted.

a vector of unknown parameters (base inertial parameters).
[18] used a novel method of selectively designing a robot
such that the resulting inertia distribution linearizes the
manipulator dynamics. In other words, the inertia parameters
are made affine in the equations of motion of the rigid body
robot. [15] used a similar method, where the affineness of
the parameters in the dynamics is leveraged to compute the
regressor.
Computing the regressor is primarily done by solving for
the dynamics for an n-DOF, b-body robot and collecting
the unknown parameters into a vector. [9] approached this
problem using an energy based approach by using the
Lagrangian formulation of robot dynamics as a starting point.
[9] also showed a second approach where the Newton-Euler
recursion method is reformulated using vector analysis–type
techniques. [1] also used Newton-Euler equations in which
the the acceleration data are obtained through a least squares
estimation and the applied torques and forces are substituted
to evaluate the regressor.
This paper uses the method adopted from [13] to evaluate
the regressor, i.e., use Spatial Vector Algebra (SVA) to
compute the regressors, an elegant way of representing the
Newton-Euler equations. If the distance between the axes are
known prior to the experiment, the parameters to be identified
are effectively the contents of the spatial inertia tensor.
These tensors are computed by shifting the inertial elements
to the joint axes. The contents of the tensor are also similar
to the D-H parameters of each link in [12]. Given a n-DOF

robot, this technique specifically lists the parameters to be
identified directly from the spatial inertia tensor. Contents of
this tensor are not the minimum representation, and therefore
will not be unique. But, it will be shown that these nonunique parameters obtained are sufficient for realizing the
model based controller, computed torque, on the robot. This
method is demonstrated on a double pendulum as well as
the leg and arm of the Valkyrie robot (Figure 1). Online
identification is done on the pendulum and the leg, and online
model based control combined with identification is done on
the arm.
We start with a brief introduction to spatial vectors in
Section II which is extracted from [3]. Representation of
kinetic energy, spatial momentum, forces and rigid body
transformations in terms of spatial vectors are also explained.
Section III shows how to use Spatial Vector Algebra (SVA)
to extract the base inertial elements and the regressor of the
robot conveniently. The derivation of this regressor and the
base inertial elements are explained in detail in Section IV.
The resulting algorithm to compute the regressor is explained
in the same section and applications to control are considered
in Section V. This is finally implemented and parameters are
identified for three models in Section VI.
II. S PATIAL V ECTOR A LGEBRA FOR A
R IGID B ODY
This section will introduce the concept of spatial vectors
and Spatial Vector Algebra. This is primarily derived from
[3] and many of the equations in this section are variants of
the equations found in the same.
Rigid body motions and forces are normally described
in two separate entities: 3D linear vectors and 3D angular
vectors. Computing linear and rotational dynamics separately
has been the normal practice in describing the equations of
motion of bodies. But, if the two 3D vectors are combined
together to form a 6D vector, a new vector space can be
described for such systems. This vector space, of course
has different rules and regulations when performing the
standard mathematical operations. The 6D vectors are formed
formally by what are called the Plücker coordinates.
A rigid body with origin located at point O, linear velocity
v, and angular velocity ω, about an axis passing through O
is shown in Figure 2. The spatial velocity of the rigid body
y
z

x

C
P
y

z
x

O
Fig. 2: Figure showing the rigid body with the origin O, an
arbitrary point P and the center of mass located at C.

can be represented as:

T
v̂O = ωx , ωy , ωz , vx , vy , vz

(1)

Similarly, the spatial force which consists of the linear force
acting at point O and the moment about the axis passing
through O can be represented as:

T
.
(2)
fˆO = τx , τy , τz , fx , fy , fz
The order of angular and translational vectors considered
is not important. Spatial vectors can also be considered
with translational vector considered first and followed by the
rotational vector.
It is important to note that the spatial vectors and forces
are independent of the origin considered and depend solely
on the bases chosen. The reason behind using spatial vectors
is that both rotations and translations can be represented in
one vector. In addition, the properties of these spatial vectors
are different and they have a different algebra. More details
about the algebra of spatial vectors can be found in [3].
Coordinate Transforms. It is important to have coordinate
transformations in order to realize rotations and translations
of coordinate frames of robotic systems using spatial vectors.
Translation. Spatial forces and velocities have the following
forms of translation from point O to an arbitrary point P .




1
0
1 −r× ˆ
v̂P =
v̂O , fˆP =
fO
(3)
−r× 1
0
1
These can be derived based on the fact that the linear velocity
vP = vO + ω × r and torque τP = τO + f × r. Here r is the
~ (see Figure 2).
position vector directed from O to P , i.e., OP
r× denotes the matrix equivalent of the cross product.
Rotation. Spatial forces and velocities have the following
forms of rotation about a point O:




ν 0
ν 0 ˆ
v̂P =
v̂O , fˆP =
fO
(4)
0 ν
0 ν
where ν indicates the rotation about the axes x, y or z or
two of them or even all of them at once. If both rotations
and translations are involved then the following relationship
is obtained:




ν
0
ν −νr× ˆ
ˆ
v̂P =
v̂O , fP =
fO . (5)
−νr× ν
0
ν
Note that the translation operation from O to P is done
first, and the the rotation about point P is carried out.
If it is required that the rotation be done first, then the
rotation matrix from (4) is multiplied with vO , and then
the translation is carried out. In doing this, it is important
to remember that the vector r also gets rotated by ν. The
resulting coordinate transformation will look like:


ν
0
v̂P =
v̂O ,
−(νr)× ν ν


ν −(νr)× ν ˆ
fˆP =
fO ,
(6)
0
ν
where (νr) is the vector r rotated by the matrix ν. Note
that (νr)× = νr×ν −1 . When this is substituted in (6) we

effectively get (5). This result will be useful in reconstructing
the spatial inertia tensors in order to conveniently evaluate
the regressor algorithm.
Momentum of a Rigid Body. If the body has mass m,
rotational inertia I¯C about its center of mass, the following
spatial momentum is described:

 

I¯C ω
I¯C
0
v̂C ,
(7)
ĥC =
=
mvC
0 m1
{z
}
|

Having the expression for momentum, the kinetic energy
can now be computed as:
T

1
1 T
ĥ v̂ = ĥTO v̂O
2
2

T

1 T ν 0
ν
v̂O
IO
0 ν
0
2

=
=

0
ν


v̂O .

(14)

III. L AGRANGIAN DYNAMICS FOR AN
N -DOF, B -B ODY ROBOT

IC

which is the product of the spatial inertia IC and the
spatial velocity v̂C . The momentum defined in (7) was w.r.t.
the center of mass. To compute the momentum about an
arbitrary point O, we have to do the transformation. So the
transformation from point O to the center of mass C is given
by:



ν 0
1 −r×
ĥC =
ĥO ,
(8)
0 ν
0
1

Equation of motion of an n-DOF manipulator is explained
in detail in this section.
A robot can be modeled as an n-link manipulator. Given
the configuration space Q ⊂ Rn , with the coordinates q ∈ Q,
and the velocities q̇ ∈ Tq Q, the Lagrangian of the n degree
of freedom robot can be defined as:

which is obtained since mvO = mvC + mr×ωC . vO , ωO
represent the spatial vectors at point O, and vC , ωC represent
the spatial vectors at point C. r is new position vector from
point O to the center of mass C (instead of P ). Expressing
ĥO in terms of ĥC , the transformation matrices get inverted:


−1
1 r×
ν 0
ĥO =
ĥC ,
(9)
0 1
0 ν

where D(q) ∈ Rn×n is the mass matrix of the robot, V (q) ∈
Rn is the potential energy of the robot. Specifically. The
equations of motion of the n-link robot can be derived as:

and using (7) in (8) and substituting for v̂C , we have:
 −1

ν
r× ν −1
ĥO =
IC v̂C
(10)
0
ν −1
 −1



ν
r× ν −1
ν
0
=
IC
v̂O .
0
ν −1
−νr× ν
If the center of mass of the rigid body is c, then at zero
rotation angle, let r = c. In other words, let r be the position
of the center of mass such that rotation of the coordinate
frame results in negative rotation of r. In other words r =
ν −1 c. Applying the trick used in (6) and substituting for
r× = (ν −1 c)× = ν −1 c × ν, we have the following result:

−1


ν 0
ν 0
ĥO =
IO
v̂O ,
(11)
0 ν
0 ν

L(q, q̇) =

1 T
q̇ D(q)q̇ − V (q),
2

D(q)q̈ + C(q, q̇)q̇ + G(q)

(15)

= BU,

(16)

where C(q, q̇) ∈ Rn×n is the matrix of Coriolis and centrifugal forces and G(q) ∈ Rn is the gravity matrix, U ∈ Rm is
the torque input with m being the number of actuators, and
B ∈ Rn×m is the mapping from actuator torques to joint
torques, often the identity map.
The first step is to determine the inertia D(q), Corioliscentrifugal C(q, q̇) and gravity G(q) matrices of the robot
via spatial vectors. We define the body spatial velocity for
each link or body i of the manipulator about the joint axis
Oi :
vOi = Ji (q)q̇,

(17)

Ji (q) is the body Jacobian of link i for the ith joint axis.
Declare the spatial inertia tensor about the joint axis Oi for
the ith link or body as IOi which is obtained from (13):


I¯Ci + mi ci × ci ×T mi ci ×
IOi =
,
(18)
mi ci ×T
mi 1

where IO :

IO =

1
0

c×
1




IC

1
c×T

0
1


,

(12)

forms the spatial inertia tensor, which is purely a function of
the parameters of the robot and independent of the orientation
of the coordinate frame considered. This equation will be
used for a general n-DOF b-body robot where the spatial
inertia tensor is effectively utilized to compute the unknown
parameters. Simplifying the spatial inertia tensor results in:


I¯C + m c× c×T m c×
I0 =
.
(13)
m c×T
m1

which is primarily obtained from the momentum equation
of (7). I¯Ci is the inertia matrix taken w.r.t. the center of
mass and mi is the mass for link i. ci is the center of mass
location of the same link w.r.t. the joint axis Oi . Accordingly,
the kinetic energy of the ith link is given by using (17).
Substituting for vOi in (17) results in:
Ti =

1 T T
q̇ Ji
2



νi
0

0
νi

T


IOi

νi
0

0
νi


Ji q̇,

(19)

where νi denotes the rotation of the ith link w.r.t. the joint
axis. The inertia matrix, D(q), can thus be expressed from

the total energy of the b bodies:
T

=

b
X

IV. T HE R EGRESSOR AND THE PARAMETERS

Ti

(20)

i=1

=

b
X

1 T
q̇
2

JiT

i=1



νi
0

0
νi

T


IOi

νi
0

0
νi

!


Ji

{z

|

q̇,
}

D(q)

where the inertia tensor IOi is obtained from (18). Note that
the Jacobian Ji is purely a function of the joint angles and the
link lengths. Therefore, assuming that the distances between
the joint axes are known (which are easy to measure), all the
other terms, namely, center of mass position, inertia, masses
are in the spatial inertia tensor IOi . This fact will be utilized
in later sections to compute the regressor.
C(q, q̇) can also be derived as a linear function of the same
elements of the tensor, by utilizing the Christoffel symbols:
Γijk is obtained from the inertia matrix D(q):


1 ∂Dij (q) ∂Dik (q) ∂Dkj (q)
+
−
Γijk =
2
qk
qj
qi
b
X
Ci (q, q̇)q̇ =
Γijk q̇j q̇k .
(21)
j,k=1

The potential energy function V (q) can be computed as
sum of the potential energies of the individual links:
V (q) =

b
X
i=1

Vi (q) =

b
X

mi ghi (q),

(22)

i=1

where mi is the mass of the individual links, g is the gravity
and hi is the vertical position of the center of mass for each
link. Therefore, hi is the sum of heights of the ith joint axis,
hOi , and the vertical height of the CoM w.r.t. the joint axis:
hi (q) = hOi (q) + hCi (q),

(23)

since the link lengths are assumed to be known, hOi in (23)
is known. The height of the CoM hCi is the dot product of
the center of mass location ci = [ci,1 , ci,2 , ci,3 ]T rotated by
a transformation matrix, ν, and the vertical axis. Assuming
that z axis is along the vertical axis, we have:

T
0
hCi (q) =  0  ν(q)ci .
(24)
1
(24) and (23) can be substituted in (22) to obtain:
 T
n
b
0
X
X
V (q) =
mi ghOi (q) +
g  0  ν(q) mi ci ,
i=1
i=1
1
{z
}
|

(25)

κ

where the unknown parameters are linear in the
and are already present in the inertia tensor IO .
in perspective, all the unknown parameters are
collected in IO , which motivates the path that
takes to compute the regressor for any general
body system.

expression
Therefore,
effectively
this paper
n-DOF b-

Since the parameters are not perfectly known, the equation
of motion, (16) computed with the given set of parameters will be henceforth have ˆ over the symbols. Therefore,
Da , Ca , Ga are the actual inertia, motor inertia, Coriolis and
gravity matrices of the robot, and D̂, Ĉ, Ĝ are the assumed
inertia, Coriolis and gravity matrices of the robot.
Consider the equation of motion of an n-link robot which
is obtained from the Lagrangian (15) and is restated here as:
K(q, q̇, q̈) = BU,

(26)

where K = D(q)q̈ + C(q, q̇)q̇ + G(q) obtained from (16).
It is a well known fact that the parameters of a robot, like
the inertia, masses, position of center of mass are affine in
(26) (see [17]). Therefore, it is possible to write (16) in the
form:
K = Y(q, q̇, q̈)Θ = U,

(27)

where Y(q, q̇, q̈) is called the regressor in [17], and Θ ∈ RnP
is called the set of base inertial elements (parameters). It
is important to note that Θ need not be unique, and if the
set is unique, then it is called the Base Parameter Set [10].
We can write, Θ = [θ1 , θ2 , θ3 , . . . ]T , where each θi is a
function of the unknown parameters of the robot. nP is the
size of the parameter set. Accordingly, Ka = Y(q, q̇, q̈)Θa ,
and K̂ = Y(q, q̇, q̈)Θ̂, where Θa is the actual set of base
inertial parameters, and Θ̂ is the assumed set of base inertial
parameters.
Determining the base inertial parameters (Θ). Since we
know the inertia of link i, I¯C ∈ R3×3 is symmetric, and the
square of a skew-symmetric matrix results in a symmetric
matrix, I¯Ci + ci ×ci ×T is symmetric. We can assign the
inertial parameters of the ith link Θi = [θi,1 , θi,2 , . . . ] to
the elements of the spatial inertia tensor IOi given in (18)
in the following manner:


θi,1
θi,2
θi,3
0
−θi,4 θi,5
 θi,2
θi,6
θi,7
θi,4
0
−θi,8 


 θi,3
θ
θ
−θ
θ
0 
i,7
i,9
i,5
i,8

,
IOi = 
θi,4 −θi,5 θi,10
0
0 
 0

 −θi,4
0
θi,8
0
θi,10
0 
θi,5 −θi,8
0
0
0
θi,10
where 10 parameters are obtained for each link. This is
similar to how the parameters were categorized in [6],
which uses the newton-euler method directly. The major
difference is that (28) is obtained from the tensors and are
used directly in online identification and control, which gets
tedious without SVA. Obtaining D(q) and C(q, q̇) from IOi
are straightforward from (20) and (21). Consider the gravity
vector, which is the partial derivative w.r.t q of the potential
energy, V (q):


n
X
∂κ
∂hOi
mi +
ci mi
(28)
G(q) =
g
∂q
∂q
i=1
Therefore, even the gravity vector G(q) is a linear function
of the inertial parameters present in the tensor.

Algorithm 1 Regressor Pseudocode
for i = 1 to n do
for j = 1 to 10 do
θi,j = 0
end for
Update IOi with the value θi,j
end for
for j = 1 to n do
for j = 1 to 10 do
θi,j = 1
Update IOi with the value θi,j
Yi+j−1 = D(q)q̈ + C(q, q̇)q̇ + G(q)
θi,j = 0
Update IOi with the value θi,j
end for
end for

Given s vectors of torque U = [u1 , u2 , . . . , un ]T , generalized
configuration q = [q1 , q2 , . . . , qn ]T , generalized velocity
data q̇ = [q̇1 , q̇2 , . . . , q̇n ]T , and generalized acceleration
q̈ = [q̈1 , q̈2 , . . . , q̈n ]T , choose model parameters Θ̂ ∈ R10b
such that1
kUC − YC Θk2

Θ̂ = argmin

(29)

Θ∈R10b

where UC is the collection of torque vector inputs and YC
is the collection of regressor matrices for s samples of angle,
velocity and acceleration data. UC and YC are given as:




Y(q[1], q̇[1], q̈[1])
U[1]
 Y(q[2], q̇[2], q̈[2]) 
 U[2] 




UC =  .  , YC = 

..


 .. 
.
Y(q[s], q̇[s], q̈[s])
U[s]
and

Comparing with [1] and [12], the inertial parameters chosen were different than the one chosen here. Specifically, the
inertial parameters in I¯C form the base inertial parameters;
whereas here the parameters in the matrix I¯Ci + ci ×ci ×T
make the unknown parameter set. Besides, it is also possible
to directly compute the regressor while evaluating the dynamics of the robot, by just picking the coefficient of every
parameter θi one by one. In fact, this becomes the basis for
a very simple algorithm shown in Algorithm 1. For a robot
having b rigid bodies, the number of unknown parameters
will be 10b.
It is shown in [3] that it is possible to compute inverse
dynamics of n-DOF robot in O(n). This is achieved by
deploying Newton-Euler recursive method. This is a standard
technique used in numerical computation of the dynamics,
which was used initially in the 1980’s. Therefore, assuming
that it is possible to compute the eom in just n recursive
iterations, we propose Algorithm 1 which calls in the current
state and acceleration of the robot and computes the regressor
from the data. The number of iterations for evaluating the
regressor is 10b, which is proportional to the number of rigid
bodies present in the robot. And the maximum number of
degrees of freedom for each rigid body is 6, which implies
that b ≤ n ≤ 6b. Accordingly in each iteration the equation
of motion is computed which takes n iterations, and the
resulting algorithmic complexity will be O(n2 ).




q[i] = 


q1 [i]
q2 [i]
..
.
qn [i]












q̇[i] = 


q̇1 [i]
q̇2 [i]
..
.









 q̈[i] = 



q̇n [i]

q̈1 [i]
q̈2 [i]
..
.




.


q̈n [i]

A vector Θ̂ that minimizes kUC − YC Θ̂k2 can be found
using the Moore-Penrose pseudoinverse:
Θ̂ = pinv(YC )UC .

(30)

An estimated (or modeled) set of computed torques ÛC can
be calculated using the parameter vector Θ̂,
ÛC = YC Θ̂.

(31)

Referring to the definition of ÛC above, a set of estimated
torques can be found for each actuator. These estimates will
be denoted [û1 , û2 , . . . , ûn ]. The coefficient of determination, R2 , can be used to describe the similarity between
UC and ÛC (or equivalently between [u1 , u2 , . . . , un ] and
[û1 , û2 , . . . , ûn ]),
e = UC − YC Θ̂
R2 = 1 − (eT e)/(UTC UC ).

(32)

Parameter Identification with an Initial Guess. An initial
guess or nominal parameter value can be incorporated into
the parameter identification problem by modifying the cost
function in (29).

V. E STIMATION AND C ONTROL
The regressor of the previous section has two main applications 1) to facilitate the identification of unknown model
parameters and 2) to enable straightforward calculation of
computed-torque controllers. This section states these problems in the general case, with specific examples to follow in
Section VI.
Parameter Identification. The problem of parameter identification can be stated as follows: Suppose we are given a
fully actuated robotic linkage with b rigid bodies, n degrees
of freedom and unknown inertial parameters Θa ∈ R10b .

Θ̂ = argmin

αkUC − YC Θk2 + (1 − α)kΘ − Θ0 k2

Θ∈R10b

(33)
where Θ0 ∈ R10b is a nominal parameter vector, and α ∈
(0, 1) is a factor that can be used to vary the relative effects
of the problem data (UC and YC ) and the initial guess (Θ0 ).
1 If necessary, the statement of the parameter identification problem (29)
can be modified to include a requirement that Θ > 0. The resulting problem
will, in general, no longer have a closed form solution and could instead
be solved using constrained quadratic programming.

A least squares solution to (33) can be found that is similar
in structure to the solution of (29),

q1

lc 1

fC )U
fC
Θ̂ = pinv(Y

(34)

m 2, IC 2

l1

q3

g

where

fC =
U

q2

m1, IC 1

α UC
(1 − α) Θ0




fC =
, Y

α YC
(1 − α) I10b



m3, IC 3

q1

.

(35)
z

Rank Properties of the Regressor. It is not necessary
to consider all of the samples to compute the parameters
because samples leading to low eigen values of the matrix
YC form a computation burden. Furthermore, it is possible to
obtain p∗ samples which give the parameter estimate Θ̂ = Θ∗
from the optimization problem (29) such that:
∗

YC (q, q̇, q̈)Θ = YC (q, q̇, q̈)Θa

(37)

The following lemma will introduce the relationship between
the desired and actual acceleration of the robot.
Lemma 1: For the fully actuated robot, i.e., B = In×n ,
if Θ̂ = Θ∗ is evaluated from the optimization problem (29)
with p∗ samples, and if the control law used is (37), then
q̈ = q̈cmd .
Note that the inertial parameters obtained from above will
not yield true parameters of the robot, but will give the same
value for the computed torque as described by Lemma 1.
This property will be used in implementing online model
based controllers and eliminate the identification of true
parameters of the robot.
VI. E XPERIMENTAL R ESULTS
This section presents three experimental studies which use
the regressor of Section IV to solve problems of identification
and the control of identified systems. The first application
will be in offline identification of a planar double pendulum,
the second in offline identification of a 3D robotic leg of
the Valkyrie robot (see Figure 1), and the third in online
identification and control of a 3D robotic arm of the Valkyrie
robot (see Figure 1).
Offline Identification of a Planar Double Pendulum.
Consider a double pendulum, such as the one illustrated
in Figure 3(a) with inertia tensors for both link 1 and 2
which are computed from (18). Note that since the pendulum
considered is planar, the number of parameters used from the
inertia tensor effectively being used is 3.
Figure 4 shows the outcome of an experiment where
angles, velocities, accelerations, and torques were measured
at each controller instant while the linkage executed a
sinusoidal motion at each joint. An estimate of the parameter
vector Θ was made without reference to an initial guess, and

y

lc 2

l2

x

(a) Illustration of
the double pendulum.

(36)

Once a parameter vector is identified, the regressor can
be used in a computed torque controller to bring about a
desired acceleration of the joints. Let q̈cmd ∈ Rn represent a
desired vector of joint accelerations, then a unique controller
to induce these accelerations is given by:
Ucmd = Y(q, q̇, q̈cmd )Θ̂.

q2

m 4, IC 4

q4

(b) Figure showing the 4-link leg.

q1

q2
m 2, IC 2

m1, IC 1

q3
z

m3, IC 3
y
x

m 4, IC 4
q4

(c) Figure showing the robot arm.

Fig. 3: Robot models considered in this paper.

thus the estimation scheme of (30) was used. The tracking
plots of Figure 4 and the scatter plots of Figure 5 illustrate
that the offline identification procedure (using all collected
data for a single parameter fit) produced a parameter set
leading to very high correlation. As noted in Section IV, the
regressor YC will not necessarily be full-rank, and thus any
estimate Θ̂ should not be expected to converge to the actual
parameters Θa , even for very large data sets.
Offline Identification of a 4-Link Robotic Leg. 4-DOF
leg shown in Figure 3 and Figure 1 was identified offline.
In the experiment a series of position, velocity, acceleration,
and torques were measured at each of the four joints of the
robot. An online estimate was made of the data, where the
parameter identification was updated at a lower, decimated
rate. The fit algorithm used here for online identification
is very similar to (34) used in offline identification of the
double pendulum. The two differences for online application
are that 1) the data vectors UC and YC grow throughout the
experiment, and as such the quality of fit improves the longer
the experiment is run, and 2) the initial parameter guess Θ0
is used to ensure bounded behavior at the beginning of the
experiment. Results of the fit are shown in Figure 6.
Online Identification of a 4-Link Robotic Arm. Online
identification was conducted on a robotic arm with four
degrees of freedom (see Figure 3(c) and Figure 1). To draw
an analogue to human physiology, the joints q1 , q2 , q3 , and
q4 can be thought of as the shoulder extensor, the shoulder

Fig. 8: Figure showing the tile of 4-DOF arm swing experiment used for identification.

Computed Torque û1 [Nm]

0.6

Angular position (rad)

5

q1
q2

0.4
0.2
0
−0.2

0
−5
−10
−15
−20
−25
−30
−30

−0.4

6

R2 J1) = 0.98

Computed Torque û2 [Nm]

0.8

−20

−10

0

Measured Torque u1 [Nm]
−0.6
0
15

50

100

150

200

Time (s)

Torque at Joint 1 [Nm]

10

250

300

350

Measured Torque u1
Computed Torque u
b1

2
0
−2
−4
−6
−6

−4

−2

0

2

4

6

Measured Torque u2 [Nm]

Fig. 5: Plotting measured torque vs. computed torque provides an

illustration of quality of fit. The coefficient of determination, R2 ,
is shown for each joint.

5
0
−5
−10
−15
−20
−25
−30
0
8

Torque at Joint 2 [Nm]

R2 (J2) = 0.81

4

50

100

150

200

Time (s)

6

250

300

350

Measured Torque u2
Computed Torque u
b2

4
2
0
−2
−4
−6
0

50

100

150

200

250

300

350

Time (s)

Fig. 4: Top: Joint trajectories as a function of time experienced
by the experimental setup of Figure 3(a). Middle, Bottom: A
comparison between experimentally measured torque vectors u1 ,
u2 and the computed torque values û1 , û2 corresponding to the
identified system.

parameters of the model. The initial guess was purposefully
chosen to be incorrect to show the convergence properties
(see Figure 7(b)) of the procedure as the nominal system
model – that estimated from engineering software – was
known with reasonable accuracy.
In order to avoid bad and potentially dangerous behavior,
a threshold of .95 was set on the coefficient of determination,
R2 , below which the nominal model would be used in place
of the identified model. The validity of the identified model
(and thereby its use in the controller) was also contingent
upon the number of data points recorded. Specifically, it
was required that the historical data buffers be full – in this
experiment, 50 historical data were used in the buffer. The
identification procedure was performed at 3 Hz and thus it
took about 16–17 seconds for the buffer to be filled. It is quite
apparent from Figure 7(c), as noted in the caption, when the
controller began to use the identified model instead of the
nominal model.
VII. C ONCLUSIONS

adductor/abductor, the upper arm pronator /supinator, and
the elbow extensor, respectively. The joint angles and even
velocities, which are necessary for computing the regressor,
can be measured from encoders but accelerations must also
be obtained and a common procedure is to filter accelerations
– in this experiment we used an exponential moving average
filter.
The experiment was run by controlling the arm to move
from one position to another and back and so forth, following
the trajectories in Figure 7(a). Tiles of the experiment are
shown in Figure 8. The system identification procedure was
brought online and was able to quickly identify the regressor

Identification and control of a physical system, in particular an n-DOF robotic system with the implementation
of an efficient computational mechanism was shown and
demonstrated on three rigid body manipulators. The regressor
involved in this implementation required a run time of O(n2 )
and computational errors resulting from this algorithm are
solely due to the error in measurement of the states of
the robots. This is a numerical method for computing the
regressor and does not use symbolic expressions which are
important for a robot like Valkyrie which has 44 degrees of
freedom. Since SVA is required for computed torque control,
evaluating the regressor through Algorithm 1 requires no
extra computational overhead. In other words, this procedure

Desired Joint Angle (rad)

5

0

Joint Torque (Nm)

−5

−10

Computed Torque û1
Computed Torque û2
Computed Torque û3
Computed Torque û4
Measured Torque u1
Measured Torque u2
Measured Torque u3
Measured Torque u4

−15

−20

−25

10

20

30

0.5
0
−0.5
q1

40

50

60

70

80

90

Time (s)

0.8

0.8
R 2 of fit

1

2

R

0.4

0.2

0.2
0
10

20

30

40
50
Time (s)

60

70

80

90

Fig. 6: The plots above show the quality of fit improves during an
online identification exercise of the robot leg pictured in Figure 3(b).
The top figure shows the instantaneous computed torque as a
function of time, along with the measured data that were used
for the parameter fit. Note that both the R2 plot at bottom and
the torque tracking plot at top both indicate that the output of the
identification algorithm leads to noticeably higher quality of fit as
the length of the available data vectors increases.

can be directly integrated within the Rigid Body Dynamics
Library [3].

60

80

5

10

15

20
25
Time (sec)

30

35

40

(b) The coefficient of determination, R2 , can be seen to converge
to right below unity. With perfect sensing and no dependence on
an incorrect initial guess, this quantity will converge precisely to
unity. Due to noise, the evolution is non-monotonic.
0.4
Joint Error (rad)

0

40
Time (sec)

0.6

0.4

0

20

(a) The arm was made to exhibit a tick tock behavior where the
joints followed the trajectories shown.

1

0.6

q3 /q4

q2

−1
0

−30
0

1

0.2
0
−0.2
q1

q2

q3

q4

−0.4

R EFERENCES
[1] Christopher G Atkeson, Chae H An, and John M Hollerbach. Estimation of inertial parameters of manipulator loads and links. Intl. J.
Robotics Research, 5(3):101–119, 1986.
[2] Björn Bukkems, Dragan Kostić, Bram de Jager, and Maarten Steinbuch. Online identification of a robot using batch adaptive control.
Proc. 13th IFAC Symp. System Identification SYSID, pages 953–958,
2003.
[3] Roy Featherstone. Rigid Body Dynamics Algorithms. Springer
Science+Business Media, LLC, 2008.
[4] Fathi Ghorbel, John Y Hung, and Mark W Spong. Adaptive control of
flexible-joint manipulators. Control Systems Magazine, IEEE, 9(7):9–
13, 1989.
[5] Roberto Horowitz, William Messner, and John B Moore. Exponential
convergence of a learning controller for robot manipulators. IEEE
Trans. Automatic Control, 36(7):890–894, 1991.
[6] Pradeep K Khosla. Categorization of parameters in the dynamic robot
model. Robotics and Automation, IEEE Transactions on, 5(3):261–
268, 1989.
[7] Pradeep K Khosla and Takeo Kanade. Parameter identification of robot
dynamics. In 24th IEEE Intl. Conf. Decision and Control, volume 24,
pages 1754–1760. IEEE, 1985.
[8] Krzysztof Kozlowski. Modelling and identification in robotics.
Springer-Verlag, New York, 1998.

0

10

20

30
40
Time (sec)

50

60

70

(c) The peak-to-peak amplitudes of the joint errors appear essentially unchanged when the controller switches to using the identified
model around 35 seconds, but the bias is clearly shifted closer to
zero.

Fig. 7: Experimental results for a 4-DOF robotic arm.

[9] W-S Lu and Q-H Meng. Regressor formulation of robot dynamics:
computation and applications. IEEE Trans. Robotics and Automation,
9(3):323–333, 1993.
[10] H. Mayeda, K. Yoshida, and K. Oshun. Base parameters of manipulator dynamic models. IEEE Trans. Robotics and Automation, 6(3):312–
321, 1990.
[11] John B Moore, Roberto Rorowitz, and William Messner. Functional
persistence of excitation and observability. In American Control
Conference, pages 308–314. IEEE, 1990.
[12] Marcelo H. Ang Jr. Ngoc Dung Vuong. Dynamic model identification
for industrial robots. Journal of Applied Sciences, 6(5):51–68, 2009.
[13] Günter Niemeyer and Jean-Jacques E Slotine. Performance in adaptive
manipulator control. The International Journal of Robotics Research,

10(2):149–161, 1991.
[14] Nicola Pedrocchi, Enrico Villagrossi, Federico Vicentini, and
Lorenzo Molinari Tosatti. On robot dynamic model identification
through sub-workspace evolved trajectories for optimal torque estimation. In Intelligent Robots and Systems (IROS), 2013 IEEE/RSJ
International Conference on, pages 2370–2376. IEEE, 2013.
[15] Shih-Ying Sheu and Michael W Walker. Identifying the independent
inertial parameter space of robot manipulators. Intl. J. Robotics
Research, 10(6):668–683, 1991.
[16] Jean-Jacques E Slotine and Weiping Li. On the adaptive control of
robot manipulators. The International Journal of Robotics Research,
6(3):49–59, 1987.
[17] Mark W Spong, Seth Hutchinson, and Mathukumalli Vidyasagar.
Robot modeling and control. John Wiley & Sons, Hoboken, NJ, 2006.
[18] DCH Yang and SW Tzeng. Simplification and linearization of
manipulator dynamics by the design of inertia distribution. Intl. J.
Robotics Research, 5(3):120–128, 1986.

