Dynamic model of robots:
Lagrangian approach
Tayachew Fikire
Dynamic model
provides the relation between
generalized forces 𝑢(𝑡) acting on the robot
robot motion, i.e.,
assumed configurations 𝑞(𝑡) over time
𝑢𝑗𝑜𝑖𝑛𝑡𝑠(𝑡) 𝑢𝐶𝑎𝑟𝑡𝑒𝑠𝑖𝑎𝑛(𝑡)
𝑞(𝑡)
a system of 2nd order
differential equations
𝛷 𝑞, 𝑞̇, 𝑞̈ = 𝑢
2
Direct dynamics
direct relation
𝑢< 𝑞<
𝑢 𝑡 = ⋮ 𝑞 𝑡 = ⋮
𝑢𝑁 𝑞𝑁
input for 𝑡 ∈ [0, 𝑇] 𝑞 0 , 𝑞̇(0) resulting motion
initial state at 𝑡 = 0
experimental solution
apply torques/forces with motors and measure joint variables
with encoders (with sampling time 𝑇𝑐 )
solution by simulation 𝛷 𝑞, 𝑞̇, 𝑞̈ = 𝑢
use dynamic model and integrate numerically the differential
equations (with simulation step 𝑇𝑠 ≤ 𝑇𝑐 )
3
Inverse dynamics
inverse relation
𝑞𝑑 𝑡 , 𝑞̇𝑑 𝑡 , 𝑞̈𝑑(𝑡) 𝑢𝑑(𝑡)
desired motion required input
for 𝑡 ∈ [0, 𝑇] for 𝑡 ∈ [0, 𝑇]
experimental solution
repeated motion trials of direct dynamics using 𝑢 𝑘 (𝑡), with
iterative learning of nominal torques updated on trial 𝑘 + 1
based on the error in [0, 𝑇] measured in trial 𝑘: lim 𝑢 𝑘 (𝑡) ⇒ 𝑢 𝑑 (𝑡)
𝑘 →G
analytic solution 𝛷 𝑞, 𝑞̇, 𝑞̈ = 𝑢
use dynamic model and compute algebraically the values 𝑢 𝑑 (𝑡) at
every time instant 𝑡
4
Approaches to dynamic modeling
Euler-Lagrange method Newton-Euler method
(energy-based approach) (balance of forces/torques)
dynamic equations in dynamic equations in
symbolic/closed form numeric/recursive form
best for study of dynamic best for implementation of
properties and analysis of control schemes (inverse
control schemes dynamics in real time)
many other formal methods based on basic principles in mechanics
are available for the derivation of the robot dynamic model:
principle of d’Alembert, of Hamilton, of virtual works, Kane’s
equations …
5
Euler-Lagrange method
(energy-based approach)
basic assumption: the 𝑁 links in motion are considered as rigid bodies
(+ later on, include also concentrated elasticity at the joints)
𝑞 ∈ ℝ𝑁 generalized coordinates (e.g., joint variables, but not only!)
Lagrangian 𝐿 𝑞, 𝑞̇ = 𝑇 𝑞, 𝑞̇ − 𝑈(𝑞)
kinetic energy – potential energy
principle of least action of Hamilton
principle of virtual works
Euler-Lagrange 𝑑 𝜕𝐿 𝜕𝐿
equations − = 𝑢i 𝑖 = 1, … , 𝑁
𝑑𝑡 𝜕𝑞̇i 𝜕𝑞i
non-conservative (external or dissipative)
generalized forces performing work on 𝑞i
6
Kinetic energy of a rigid body
mass density
𝒗𝒄 body 𝐵
mass 𝑚 = z 𝜌 𝑥, 𝑦, 𝑧 𝑑𝑥𝑑𝑦𝑑𝑧 = z 𝑑𝑚
𝑚 𝐵 𝐵
1
position of 𝑟𝑐 = z 𝑟 𝑑𝑚
center of mass (CoM) 𝑚 𝐵
RFc when all vectors are referred to a body frame
RFc attached to the CoM, then
𝑟𝑐 𝑟𝑐 = 0 ⇒ z 𝑟 𝑑𝑚 = 0
𝐵
RF0 1
kinetic energy 𝑇 = z 𝑣 𝑇 𝑥, 𝑦, 𝑧 𝑣(𝑥, 𝑦, 𝑧) 𝑑𝑚
2 𝐵
(fundamental)
kinematic relation 𝑣 = 𝑣𝑐 + 𝜔 × 𝑟 = 𝑣𝑐 + 𝑆(𝜔) 𝑟
for a rigid body
skew-symmetric matrix
7
Kinetic energy of a rigid body (cont)
1 𝑇
𝑇= z 𝑣𝑐 + 𝑆 𝜔 𝑟 𝑣𝑐 + 𝑆 𝜔 𝑟 𝑑𝑚
2 𝐵
<
= ∫ 𝑣𝑐𝑇 𝑣𝑐 𝑑𝑚 + ∫ 𝑣𝑐𝑇 𝑆 𝜔 𝑟 𝑑𝑚 + < ∫ 𝑟 𝑇 𝑆 𝑇 𝜔 𝑆 𝜔 𝑟 𝑑𝑚
a 𝐵 𝐵 a 𝐵
1 1
= 𝑚 𝑣𝑐𝑇 𝑣𝑐 = 𝑣𝑐𝑇 𝑆 𝜔 z 𝑟𝑑𝑚 = 0 = z 𝜔 𝑇 𝑆 𝑇 𝑟 𝑆 𝑟 𝜔 𝑑𝑚
2 𝐵
2 𝐵
1 𝑇
= 𝜔 z 𝑆 𝑇 𝑟 𝑆 𝑟 𝑑𝑚 𝜔
translational 2 𝐵
kinetic energy rotational
(point mass + kinetic energy
1 𝑇
= 𝜔 𝐼𝑐 𝜔
Homework
provide the expression
at CoM) (of the whole body) 2 of the elements
of the inertia matrix 𝐼𝑐
König theorem body inertia matrix
(around the CoM)
8
Robot kinetic energy
𝑁
𝑇 = Ö 𝑇i 𝑁 rigid bodies (+ fixed base)
iÜ<
𝑇i = 𝑇i(𝑞j, 𝑞̇j; 𝑗 ≤ 𝑖) open kinematic chain
𝒗𝒄𝒊
König theorem
𝑚𝒊 1 1 𝑇
𝑇
𝑇i = 𝑚i 𝑣𝑐 𝑣𝑐 + 𝜔i 𝐼𝑐 i 𝜔i
2 i i
2
RFci
absolute velocity absolute
of the center of mass angular velocity
i-th link (body) (CoM) of whole body
of the robot
9
Kinetic energy of a robot link
1 𝑇 1 𝑇
𝑇i = 𝑚i 𝑣𝑐 𝑣𝑐 + 𝜔i 𝐼𝑐i 𝜔i
i
2 i
2
𝜔𝑖, 𝐼𝑐𝑖 should be expressed in the same reference frame,
but the product 𝜔𝑖𝑇𝐼𝑐𝑖𝜔𝑖 is invariant w.r.t. any chosen frame
[
[
𝜔 𝑇i [ 𝐼 𝑐 i (𝑞) [ 𝜔 i = [ i
𝑅 i (𝑞) 𝜔 i
𝑇 [
𝐼 𝑐i (𝑞) [ i
𝑅 i (𝑞) 𝜔 i i
= 𝜔 𝑇i 𝑅 𝑇 (𝑞) [ 𝐼 (𝑞) [ 𝑅 (𝑞) i
𝜔i
i 𝑐i i
= i 𝜔 𝑇i i𝐼𝑐i i 𝜔 i [𝐼
𝑐i 𝑞 = [𝑅i 𝑞 i𝐼 [ 𝑇
𝑐 i 𝑅 i (𝑞)
in frame RFci attached to (the center of mass of) link 𝑖
z 𝑦 a + 𝑧 a 𝑑𝑚 − z 𝑥𝑦 𝑑𝑚 − z 𝑥𝑧 𝑑𝑚
i
𝐼𝑐i = z 𝑥 a + 𝑧 a 𝑑𝑚 − z 𝑦𝑧 𝑑𝑚
constant! 𝑠𝑦𝑚𝑚 z 𝑥 a + 𝑦 a 𝑑𝑚
10
Examples of body inertia matrices
homogeneous bodies of mass 𝑚, with axes of symmetry
parallelepiped with sides
y 𝑎 (length/height), 𝑏 and 𝑐 (base)
𝑎
𝑐 <
𝑏a + 𝑐 a
x 𝐼𝑥𝑥 <a
𝑚
𝑏 𝐼𝑐 = 𝐼𝑦𝑦 = <
𝑚 𝑎a + 𝑐 a
<a
𝐼𝑧𝑧
z <
<a
𝑚 𝑎a + 𝑏a
empty cylinder with length ℎ,
y and external/internal radius 𝑎 and 𝑏
ℎ <
𝑏 𝑚 𝑎a + 𝑏a
x a
𝐼𝑐 = <
𝑚 3 𝑎a + 𝑏a + ℎa 𝐼𝑧𝑧 = 𝐼𝑦𝑦
<a
𝑎 𝐼𝑧
z’ z a 𝑧
ℎ
u
𝐼𝑧𝑧 = 𝐼𝑧𝑧 + 𝑚
2 (parallel) axis translation theorem
Steiner theorem
… its generalization:
𝐼 = 𝐼𝑐 + 𝑚 𝑟 𝑇 𝑟 k𝐸w×w − 𝑟𝑟 𝑇 = 𝐼𝑐 + 𝑚 𝑆 𝑇 𝑟 𝑆 𝑟 changes on body inertia matrix
due to a pure translation 𝑟 of
body inertia matrix identity skew- the reference frame
relative to the CoM matrix symmetric
10
Moment of Inertia tensor
Relation ship between Angular Momentum and angular Velocity
• The tensor of inertia gives us an idea
about how the mass is distributed in a
rigid body.
• Ixx,Iyy,Izz are moment of inertia about x,y,
and z axis.
• Ixy, Ixz,Iyz, are called product of inertia
and they are measure of the imbalance in
the mass distribution.
12
Robot potential energy
assumption: GRAVITY contribution only
𝑁
𝑈 = Ö 𝑈i 𝑁 rigid bodies (+ fixed base)
iÜ<
𝑈i = 𝑈i (𝑞j ; 𝑗 ≤ 𝑖) open kinematic chain
𝑈i = −𝑚 i 𝑔 𝑇 𝑟 [,𝑐i
typically
gravity acceleration position of the expressed
vector center of mass of link 𝑖
in RF0
dependence on 𝑞
𝑟[ ,𝑐 𝑟i ,𝑐 constant
[ < i– < in RFi
= 𝐴<(𝑞<) 𝐴a(𝑞a) ⋯ 𝐴 i (𝑞 i )
i i
1 1
NOTE: need to work with homogeneous coordinates
13
Robot potential energy
assumption: GRAVITY contribution only
𝑁
𝑈 = Ö 𝑈i 𝑁 rigid bodies (+ fixed base)
iÜ<
𝑈i = 𝑈i (𝑞j ; 𝑗 ≤ 𝑖) open kinematic chain
𝑈i = −𝑚 i 𝑔 𝑇 𝑟 [,𝑐i
typically
gravity acceleration position of the expressed
vector center of mass of link 𝑖
in RF0
dependence on 𝑞
𝑟[ ,𝑐 𝑟i ,𝑐 constant
[ < i– < in RFi
= 𝐴<(𝑞<) 𝐴a(𝑞a) ⋯ 𝐴 i (𝑞 i )
i i
1 1
NOTE: need to work with homogeneous coordinates
14
Robot dynamic model
in vector formats
1. 𝑀 𝑞 𝑞̈ + 𝑐 𝑞, 𝑞̇ + 𝑔 𝑞 = 𝑢
𝑐𝑘 𝑞, 𝑞̇ = 𝑞̇𝑇𝐶𝑘(𝑞)𝑞̇ 𝑘-th component
of ve ctor 𝑐
𝑘-th column 𝑇
1 𝜕𝑀𝑘
of matrix 𝑀(𝑞) 𝐶 (𝑞) = + 𝜕𝑀𝑘 −
𝜕𝑀 symmetric
𝑘
2 𝜕𝑞 𝜕𝑞 𝜕𝑞𝑘 matrix!
2. 𝑀 𝑞 𝑞̈ + 𝑆 𝑞, 𝑞̇ 𝑞̇ + 𝑔 𝑞 = 𝑢
NOTE:
the model
is in the form NOT a factorization of 𝑐
𝛷 𝑞, 𝑞̇, 𝑞̈ = 𝑢 symmetric 𝑠 𝑘j 𝑞, 𝑞̇ = Ö 𝑐 𝑘ij 𝑞 𝑞̇i by 𝑆 is not unique!
as expected matrix i
in general
15
… and interpretation of dynamic terms
𝜕𝑈
Ö 𝑚 𝑘j (𝑞)𝑞̈ j + Ö 𝑐 𝑘ij 𝑞 𝑞̇i𝑞̇j + = 𝑢𝑘 𝑘 = 1, … , 𝑁
j i ,j
𝜕𝑞𝑘
INERTIAL CENTRIFUGAL (𝑖 = 𝑗) GRAVITY
terms and CORIOLIS (𝑖 ≠ 𝑗) terms terms 𝑔𝑘(𝑞)
𝑚𝑘𝑘(𝑞) = inertia at joint 𝑘 when joint 𝑘 accelerates (𝑚𝑘𝑘 > 0!!)
𝑚 𝑘j (𝑞) = inertia “seen” at joint 𝑘 when joint 𝑗 accelerates (= 𝑚j𝑘(𝑞))
𝑐𝑘ii(𝑞) = coefficient of the centrifugal force at joint 𝑘 when
joint 𝑖 is moving (𝑐iii = 0, ∀𝑖 )
𝑐 𝑘ij 𝑞 = coefficient of the Coriolis force at joint 𝑘 when joint 𝑖
and joint 𝑗 are both moving (= 𝑐 𝑘ji (𝑞))
16
Dynamic model of a PR robot
𝑦 𝑇 = 𝑇< + 𝑇a 𝑈 = constant ⇒ 𝑔(𝑞) ≡ 0
(on horizontal plane)
𝑑 𝑐a 𝑝𝑐2 𝑞< − 𝑑𝑐< a = 𝑝̇𝑇 𝑝̇𝑐< = 𝑞̇a
𝑝𝑐 < = 0 𝑣𝑐< 𝑐< <
𝑑𝑐1 𝑞2 0
1
𝑝𝑐1 𝑥 𝑇< = 𝑚<𝑞̇<a
2
𝑞1
1 𝑇
1 𝑇
𝑇a = 𝑚a 𝑣𝑐 a𝑣𝑐 + 𝜔a 𝐼𝑐 a 𝜔a
2 a 2
𝑞< + 𝑑 𝑐a cos 𝑞a 𝑞̇< − 𝑑 𝑐a sin 𝑞a 𝑞̇a 0
𝑝 𝑐a = 𝑑𝑐 a sin 𝑞a 𝑣 𝑐a = 𝑑 𝑐a cos 𝑞a 𝑞̇a 𝜔a = 0
0 0 𝑞̇a
𝑇a = < 𝑚a 𝑞̇<
a + 𝑑 a 𝑞̇ a − 2𝑑 sin 𝑞 𝑞̇ 𝑞̇ + < 𝐼
𝑐a a 𝑐a a < a 𝑞̇
a 𝑐a,𝑧𝑧 a
a
a
17
Dynamic model of a PR robot (cont)
𝑚< + 𝑚a −𝑚a 𝑑𝑐 a sin 𝑞a 𝑐< 𝑞, 𝑞
𝑀 𝑞 = 𝑐 𝑞, 𝑞̇ =
−𝑚a 𝑑𝑐 a sin 𝑞a 𝐼𝑐a,𝑧𝑧 + 𝑚 a 𝑑 𝑐a 𝑐a 𝑞, 𝑞
𝑀< 𝑀a
a 𝑐𝑘 𝑞, 𝑞̇ = 𝑞̇ 𝑇 𝐶𝑘 (𝑞)𝑞
𝑇
1 𝜕𝑀𝑘 𝜕𝑀𝑘 𝜕𝑀
where 𝐶𝑘(𝑞) =
2 𝜕𝑞 + 𝜕𝑞 −
𝜕𝑞𝑘
1 0 0 0 0
𝐶< 𝑞 = + − 0 0
2 0 −𝑚a 𝑑𝑐 a cos 𝑞a 0 −𝑚a 𝑑𝑐 a cos 𝑞a 0 0
𝑐< 𝑞, 𝑞̇ = −𝑚a𝑑𝑐a cos 𝑞a 𝑞̇aa
0 −𝑚a 𝑑𝑐 a cos 𝑞a 0 0
1 +
0 0 −𝑚a𝑑𝑐a cos 𝑞a 0
𝐶a 𝑞 = =0
2 0 −𝑚a 𝑑𝑐 a cos 𝑞a
−
−𝑚a 𝑑𝑐 a cos 𝑞a 0
𝑐a 𝑞, 𝑞̇ = 0
18
Dynamic model of a PR robot (cont)
𝑀 𝑞 𝑞̈ + 𝑐 𝑞, 𝑞̇ = 𝑢
𝑚< + 𝑚a −𝑚a 𝑑𝑐 a sin 𝑞a 𝑞̈< −𝑚 𝑑 cos 𝑞 𝑞̇ a 𝑢<
+ a 𝑐 a a = 𝑢
−𝑚a𝑑𝑐a sin 𝑞a 𝐼𝑐a,𝑧𝑧 + 𝑚 a 𝑑 a 𝑞̈a a 0 a
𝑐a
NOTE: the 𝑚𝑁𝑁 element (here, for 𝑁 = 2) of 𝑀(𝑞) is always constant!
Q1: why does variable 𝑞<not appear in 𝑀(𝑞)? … this is a general property!
Q2: why Coriolis terms are not present?
Q3: when applying a force 𝑢<, does the second joint accelerate? … always?
Q4: what is the expression of a factorization matrix 𝑆? … is it unique here?
Q5: which is the configuration with “maximum inertia”?
19
Dynamics of an actuated pendulum
a first example
𝑏𝑚 𝜃̇𝑚 = 𝑛𝑟𝜃 ̇ 𝜃𝑚 = 𝑛𝑟𝜃 + 𝜃 𝑚[
viscous friction 𝐼𝑚
𝑟 = 𝑛𝑟 𝑟𝑚 =0
𝑛𝑟 ≥ 1 𝑟𝑚, 𝜃𝑚
𝑏𝑙
motor 𝑞 = 𝜃 (or 𝑞 = 𝜃𝑚 )
transmission
𝑑 (with reduction gear) (… around the ∥ axis
𝑇 = 𝑇𝑚 + 𝑇𝑙 through its base)
𝑟, 𝜃 𝑚 1 1
𝑙 𝑇𝑚 = 𝐼𝑚 𝜃̇𝑚
a 𝑇𝑙 = 𝐼𝑙 + 𝑚𝑑 a 𝜃̇ a
𝐼𝑙 𝐹𝑥 2 2
𝑚𝑔 [ motor inertia link inertia
(around its (around the 𝑧-axis through
spinning axis) its center of mass…)
1
𝑇= 𝐼 + 𝑚𝑑 a + 𝐼 𝑛a 𝜃̇ a = 1 𝐼𝜃̇ a
kinetic energy
2 𝑙 𝑚 𝑟 2
2
0
Dynamics of an actuated pendulum (cont)
𝑦 𝑟𝑚 , 𝜃𝑚 𝑈 = 𝑈[ − 𝑚𝑔[ 𝑑 cos 𝜃 potential energy
𝑥
− 𝑑 cos𝜃
1 a
𝑔[ 𝐿 = 𝑇 − 𝑈 = 𝐼𝜃̇ + 𝑚𝑔 [ 𝑑 cos 𝜃 −𝑈[
2
𝑛𝑟 𝑟𝑚 , 𝜃
𝑚 𝐹𝑥
𝜕𝐿 𝑑 𝜕𝐿 𝜕𝐿
𝑝𝑥 = 𝑙 sin 𝜃 = 𝐼𝜃 = 𝐼𝜃 = −𝑚𝑔[ 𝑑 sin 𝜃
𝜕𝜃 𝑑𝑡 𝜕𝜃 𝜕𝜃
𝑝̇𝑥 = 𝑙 cos 𝜃 k𝜃̇ = 𝐽𝑥 𝜃̇
𝑢 = 𝑛𝑟 𝑟𝑚 − 𝑏𝑙 𝜃̇ − 𝑛𝑟 𝑏𝑚 𝜃̇𝑚 + 𝐽𝑇 𝐹𝑥 = 𝑛𝑟 𝑟𝑚 − 𝑏𝑙 + 𝑏𝑚 𝑛a 𝜃̇ + 𝑙 cos 𝜃 𝐹𝑥
𝑥 𝑟
“sum” of
applied or dissipated torques equivalent joint torque non-conservative
on motor side are multiplied by 𝑛𝑟 due to force 𝐹𝑥 applied to torques
when moved to the link side the tip at point 𝑝𝑥
2
1
Dynamics of an actuated pendulum (cont)
𝑟𝑚 , 𝜃𝑚
dynamic model in 𝑞 = 𝜃
𝐼𝜃̈ + 𝑚𝑔[ 𝑑 sin 𝜃 = 𝑛𝑟 𝑟𝑚 − 𝑏𝑙 + 𝑏𝑚 𝑛𝑟a 𝜃̇ + 𝑙 cos 𝜃 k𝐹𝑥
𝑛𝑟 𝑟𝑚 , 𝜃
dividing by 𝑛𝑟 and substituting 𝜃 = 𝜃𝑚/𝑛𝑟
𝐼 𝑚 𝜃𝑚 𝑏𝑙 𝑙 𝜃𝑚
a
𝜃𝑚 + 𝑔[ 𝑑 sin = 𝑟𝑚 − a + 𝑏𝑚 𝜃𝑚 + cos k𝐹𝑥
𝑛𝑟 𝑛𝑟 𝑛𝑟 𝑛𝑟 𝑛𝑟 𝑛𝑟
dynamic model in 𝑞 = 𝜃𝑚
2
2
Example 2. Dynamic model of 2 link planar Robot
23
Example 2. Dynamic model of 2 link planar Robot
24
Example 2. Dynamic model of 2 link planar Robot
25
