Copyright
SESSION WAl -
ROBOTICS
RECURSIVE LAGRANGIAN DYNAMICS OF FLEXIBLE MANIPULATOR ARMS VIA TRANSFORMATION MATRICES W.
J. Book
The Robotics Institute, Carnegie-Mellon University, Pittsburgh, PA, USA and School of Mechanical Engineering, Georgia Institute of Technology, Atlanta, GA , USA
Abstract The non linear equations of motion for flexible manipulator arms consisting of rotary joints connecting two flexible links arc developcd ,
Kinematics of both the rotary joint motion and the link
dcformation arc dcscribcd by 4x4 transformation matrices. Thc link deflcction is assumed small so that the link transformation can bc composcd of summations of assumed link shapes. Thc rcsulting equations are presented as scalar and 4x4 matrix operations ready for programming. The efficiency of this fonnulation is compared to rigid link cases reported in the literature.
Keywords: Robots: Distributed parameter systems; Models; Manipulation: Vibration control; Flexible mechanisms; Mechanical arms. The author gratefully acknowledges the support of this work by The Robotics Institute. CamegicMellon University and the Gcorgia Institute of Technology while thc author was on lcavc as Visiting Scicntist. Thc cncouragemcnt of Marc Raibert of The Robotics Institute was also important.
1. Introduction
1.1. Sketch of Prior Work
Increased ilccuracy in m(ldeling system behavior is
Much work has bccn donc in thc formulation of the
needed when improved performance is sought for most
dynamic cquations of motion for mcchanical arms with rigid
enginecring systcms.
links. The reader will find dcscribed in referenccs [21]. [26].
The evolution of the mechanical arm
from telcoperator and crane to present ddY industrial and space
[27]. [2] and in thcir bibliographics thc work on the "inverse
robots and largc spacc manipulators is no cxccption. Initial
dynamic
simplc kinel1lJtic and d ynamic models arc
[30]. [11]. and thcir bibliogr3phics thc work on the dynamic
110
longer adequate
for the most critical performance demands.
Both the
design simubtion.
Proposcd
new
in control . and in [28]. [19]. [31].
formubtion for simulation of rigid link arms. Thc cfficicncy of
mcchanical system and cOlllrol system require improved models for
fonnu~Jtion" u ~cd
thcse'
control
formulations
and
alternativcs
to
thcir
real
time
calculation is discusscd in [25]. [1] and the works referred to
algoritJlm; require dynamic model s for control calculation.
thcrein.
Planning and programming activities as well as m an-in-theloop simulation also requirc accurate modcls of the arm.
Thc limitation of thesc works is that rigid links are assumcd.
Accuracy is usually acquired at some cost.
The
With this assumption me tcchniques become at
some point sclf defcating. if their purpose is to improve
application of mechanical arms to economically sensitive
performance. Maintaining rigidity
endeavors in industry and incrcased cost consciousness in the
improvcd pcrformanee cvcn as tne modcl accuracy tends to
spacc program also gives incentive to improve the cffieieney of
improve it.
ofthc links inhibits
the formulation and simulation of the dynamic models. Control algorithms and man in the loop simulation rcquire
Considcration of flexibility of the links in arm type
"rcal timc" calculation of dynamic bchavior. Formulation of
devices and thcir control is rcportcd by 1972 in Mirro [23].
the dynamics in an casy to undcrstand conccptual approach is
This carly work considcrcd both the modeling and control of a
also important if maximum utilization of the results is to be
singlc link dcvice. Book (7] considcrcd the linear dynamics of
obtained.
spatial flcxible arms reprcsentcd as lumped mass and spring
5
w.
6
J. Book
components via 4x4 transformation matrices. This was refined
and open loop chains such as manipulator arms. This work
and later reported in [9).
Book and Whitney [3), [4] later
assumes a known nominal motion over time about which the
considered linear distributed dynamics of planar arms via
flexible arm equations arc !inearized. This falls short of a true
transfer matrices and the limitations flexibility imposed on
simulation of the flexible, non linear equations, but is an
control system performance [8]. Maizza and Whitney [22], [4]
interesting compromise in the interest of computational speed.
utilized a planar nonlinear model with modal representation of
This technique is oriented toward finite clement analysis to
the flexibility and considered modal control as a technique for
obtain modal characteristics of the links which arc then
overcoming the limitations of the flexibility. Whitney, Book,
combined using a time varying compatibility matrix. It uses
and Lynch [32], [4] considered the design implications of
4x4 matrices to represent the nominal kinematics and
flexibility.
derivation of the compatibility matrix.
Distrihuted
frequency
domain
analysis of
non planar arms using transfer matrix techniques [5], [6] has been used by Book. et.al in verifying the accuracy of truncated modal models of the non linear spatial dynamics of flexible manipulators
(the
Remote
Manipulator
Shuttle).A more classical appro
to
of the
Space
manipulator dynamics
both rigid [l7] and flexible [IS] has been undertaken by HusLOn and his coworkers.
1 . 2 . Pe rspective on This Wo rk
The Y-ork in this p
to
workers
in the field of robotics. It is unique in sCI'eral respects. It uses 4x4 matrices
to
repn'sent both the joint and deflrction motion.
The deflection transfurmation is represented in terms of a The work in flexible spacecraft has spawned a line of research pertaining to the interaction of articulated structures which has great relevance to the manipulator modeling problem. Entries into Lhis literature arc provided by the works by l.ikins [20) and Hughcs 124). This activitv produced a spatial. nonlincar. flcxiblc manipulator model repoTted by Ho el.al. (13) and corresponding computcr code for
~imulation.
The simubtion required great amounts of computer time and was unsuitable for even off line simulation. Further work for thc purposes of simulating the Space Shuttle Remote Manipulator was pcrformed bl' Hughcs. His lin carizcd modcl is reported in [J 5) and a more general model is reported in [16). The Hughes model ignores the interaction between structural
summation of modal shapes. The computations resulting from the l.agrang ian formulation of the dynamics are reduced to recursive form similar to lhat wh ieh has proven so efficient in the rigid link case. The eq uatio ns arc free from assumptions of a nominal moti on. and do not ignore the interaction of angular rates and deflections. They do assume small deflections of the link s which can be desc ribed by a sum mation of the modal shapes and a linear model of elasticity. Only rotational joints arc allowed.
The results arc quite tractable for automated
computer solution for arbitrary rotary joints. programs have been written
to
Preliminary
evaluate computational
efficiency.
deformation and angubr rate as might be appropriate for the Space Shuttle atm. This work and associated work at SPAR Aerospace, Ltd. and the Charles Stark Draper Laboratory, Inc.
2. Flexible Arm Kinematics The previous works on rigid ann dynamics make heavy
probahly represents the most intensive work on the modeling,
use of the serial nature of manipulator arms which results in
simulation and control (lf fl exible arms. Unfortun,nely, little of
multiplicative
this work has been rep0rled in the open literature. Recent
representation of flexible structure dynamics. on the other
examination of experimental results from the operation of the
hand. is a parallel or additive represe ntation of the system
Shuttle arm in space has confirmed the validity of these
behavior. One of the contributions of this paper is
models.
this difference in a concise way. As with many of the previous
terms
in
the
kinematics.
The
to
modal
resolve
works on rigid dyna mics, the 4x4 matrices of Denavit and Yet another branch of research that has found its way to
Hartenberg [10) are used.
Sunada and Dubowsky [291 used
the flexible manipulator dynamics problem is the study of
this representation for their flexible arm simulations but did
Dubowsky and Gardner (12) and
not produce a comolete nonlinear dynamic simulation. Other
Winfrey [33] will provide the reader with a bibliography on this
workers such as Hughes [16) relied on the more general
work. Sunacta and Dubowsky [29] have developed modeling
formulation provided by a vector-dyadic represcntation.
techniques applicable to both spatial closed loop mechanisms
While Silver 126], Hollerbach (14). and others have pointed out
flexible mechanisms.
Recursive Lagrangian Dynamics
the relative inefficiency of the 4x4 formulation. the conceptual fr~mework
is most ad\antagcous when tackling the complexity
of the flexible dynamics.
7
To incorporate the deflection of the link. the approach of modal analysis is used which is valid for small deflection of the link.
Define the position of a point in Cartesian coordinates by
(4) an augmented vector:
[I x component y component z componentf Define the coordinate system [x y zl; on link i with origin 0; at
where x. y .. . 1. = the x. y .. and z displacement IJ IJ IJ I I I components of mode j of link i's deflection, respectively. (j = the time varying amplitude ofmodej of link i ~)I = the number of modes used to describe the deflection of link i.
the proximal end (nearest the base) oriented so that the x axis is coincident with the neutral axis of the beam in its undeformed condition. The orientation of the remaining axes will be done so as to allow efficient description of the joint motion. A point on the neutral Jxis at x = 1) wllen the beam is undeformed is locatcd at ;h(1)) under a general condition of I deformation. in terms of system i.
The link transformation matrix must also incorporate the Here the rotations as wen as the
deflection of the link.
translations of the deflection must be represented.
By a homogeneous transformation of coordinates the position of a point can be desc ribed in any other coordinate system j if the transfonnation matrix jw; is known. The form
If one
consistently requires small rotations the dircction cosine matrix simplifies as noted in [91 and furthermore the small angles can be assumed to add vectorally. This is basic to the approach used here. The link transformation matrix can then be written
of this matrix is
as m.
'w, =
[
Yj Zj
1
(1)
Xj component ofO;
component ofO; component ofO;
E;
= [ H;
+
2:
(jij
(5)
Mu ]
j=l
where where
jR; = a 3x3 matrix of direction cosines o = a ]x3 vector of zeros.
H I. =
Thus in terms of the fixed inertial coordinates of the base the position of a point on link i is given as
[
(2) where the special case ofOW.I
= W. I
M. =
It is usefL!I to separate the
11
transformations due to the joint from the transformation due to the flexihle link as follows
where
A = the joint tran sfo rmation matrix for joint j J E 1 = the link tran sformation matrix for link j-l ) between joints j-1 and j W ] = th e clImubti\e tran sfo rmatioll from base J"coordinateslO C\.] at the distal end oflinkj.
Dj.]
0 0
1
0 0
0 0
1
0 0 0
0
1
0
0 0
-0
x. I)
0
J';)
o.
zIJ
-Oy;j
O il)
1I)
(6)
1
~"j 1 o
(7)
-0
XI)
-OXij
and where (3)
A
[
0 II
,\11 variables in hrdckets are c\alllated at 1.I O'if (j yij' 0 Zlj = the x,, Y;' and z; rotation components of link i, respectively. 1; = the length of link i. To find the velocity of a point on link i, take the time
•
is fixed to the link j-l and with no deflec:ion
parallel to [x y 1.lj.] with \1 coincident with x j.l '
[~
y ~ lj-l
deri\'3ti\e of tile position: dh = 1 dt
-
is
h.1 = \V I ;h I + W 1 ;h 1..
(8)
Due to the serial nature of the kinematic chain it is computationally efficient to rclate the position of a point and its derivatives to preceeding members in the chain.
By
w.
8
J. Book
differentiating 2 one obtains: .
;..
3. System Kinetic Energy In this section the expression for Ule system kinetic
A'
(9)
Wj = Wj_l .\ + Wj_l Aj 'Vj = '\-1 '\ + 2 \"j-l''\j + " 'j-l
(10)
Aj
energy is developed for use in Lagrange's equations.
The
kinetic energy for a differential clement is written then integrated over the link. This produces tenns which are the
where (11)
equivalent of the moment of inertia matrices of rigid link anns. Summation over all the links provides the total kinetic energy.
(12)
The kinetic energy of a point on the i-th link is
. { h.' 1 h'T} dk.I = - Idmlr I 2 where
q. = the joint variable ofjointj. J
.
(18)
..
dm is ule differential mass of the point and Tr{.} is the trace operator.
"'-
Thus W and \Vj can be computed recursively from Wj_l . J its derivatives. and the partials with respect to the variables of
Expanding 18 and using the fact that Tr{A BT} = Tr{B AT}>
link j-] and joint j. No mixed p<;rtials are explicitly present
the expression for dk i becomes
This computational approach is similar to that proposed by Hollerbach [14] for rigid link alms. needs Wj_1 and its derivatives.
Here onc additionally
These can be computed where
recursively from \Vj_1 and its derivatives: (13)
Wj=WjEj J..
•
(14)
Wj = Wj Ej + \Vj Ej
.. Wj = \Vj Ej + 2 Wj Ej + Wj Ej :;.:
..
(15)
kinetic energy. In this paper it is assumed that the links are slender beams because it makes the central development
mj
F:k =
1:
(20) Sij [ 0 Xij Yij lij ]T. j=1 By integrating over the link one can obtain the total link
L: Sjk Mjk
(16)
clearer. Other mass distributions could be used with a slight
k=l
departure here in the development. For slender beams dm =
mj
J.L
j.\ = L:
(17)
Sjk Mjk
d1) and one can integrate over 1) from 0 to li' Only the tenns
in ihJ and its derivatives are functions of 1) for this link. Thus the integration can be perfomlCd without knowledge ofW i and its derivative. Summing over an n links onc finds the system
k=1 One can see from the last two equations that the deflection transformations enter even more simply into the
kinetic energy to be n
kinematics on a per variable basis than do Ule joint variables.
K=
K=
the
terms
involving
~econd
derivatives of the joint and deflection variables will be
derivatives of the state variables. The "inverse dynamics" solution that proceeds directly from the Lagrange fonnulation has little obvious utility.
(22)
where
separated from the above expressions and included in the inertia matrix to make up the coefficient matrix of the
L:
Tr{ \', B3i \"; + 2\"i B2i W; i=l +WiBliW;}
velocity and acceleration is preserved from the rigid case. For equations
(21)
dk i
n
chosen for the transformation. The recursive nature of the
simulation
/i 0
i=1
This is due to the small deflection assumption and the form
the
L:
BJj
=.l/ 2
I. J
.. "T J.L \
\
d1) .
(23)
0
By interchanging the integration in 23 and the summations involved in the definition ohi i in 20 one obtains m.J Bli =
1:
j=l
L:
k=l
Sij Sik Cikj
(24)
Recursive Lagrangian Dynamics
where
I.
Cikj =
.l/1 2
0
9
3.1 . Oerriatives of Kinetic Energy
Jl [0 \k Yik
li/ [0 \j Yij Zij] d1).
(25)
It should be noticed that Cikj has units of an inertia matrix and serves a similar function . While shown here as a 4x4 matrix it is nonzero only in the 3x3 lower right hand corner. It can also
For construction of Lagrange's equations one needs aK IOq . , aK / a8! ' ~ ( aK ) ) dt
IOq) ) , and ~dt ( aK / (8)()
First consider oK / aq . This will involve the partials of )
all tlle terms in 22, some of which are zero. In fact, only
Wi for
be shown that Cikj = CJk ' It is possible by choosing the assumed mode shapes in an appropriate manner, to reduce the
j $ i $ n provides nonzero partials with respect to qj' The time
number of nonzero terms in 24. This matter is discussed in
derivative of the partial is then taken.
light of computational speed in thc conclusions.
following equivalences should be noted:
The other terms in equation 22 can similarly be found: li B2i
=.ll
(26)
2
1
= aw / oq.
/Oq )
1
(31)
)
d (. .) . aw /oq . = aw / aq . I) I) dt
(32)
0
m.I
mi
L:
B2i =
aw
In this respect the
8 ij Cij +
j= 1
(33)
m.
L: ~ 8 ik 8ij Cikj j=l
k=l
(27)
(34)
where I.
Cij =
.l/1 2
Also helpful in simplifying the result is that. Tr{A} = Tr{A T} Jl [1 1) 0 O]T [0 Xij Yij zij] d1) .
(28)
0
Considerable cancellation and combination when the terms in
Finally, by a similar approach:
Lagrangc's equation
Ij
B3i
=.ll 2
0
~ ( dt
+
L:
8 ij [C ik
+ Cr]
L: k=l
i=j
kinetic
energy
are
Oq) - 0 K / oq
=
Tr { 0 Wi oqj
C.
)
)
[[
+
I
(29)
8 it 8 jj Cikj
j=l
where
2
aK /
L:
2
mi
.ll
the
n
j= 1
Ci =
involving
combined . The result of this combination is mi
B3i = Ci
for any square matrix A and that B3i is symmetric.
1=1
k=l
\ Jl [11) 0 O]T [11) 00] d1).
(30)
~ik (C jk +
0
mi
one link in the system is to be considered rigid, in which mi = O. Should a link consist of a flexible member with rigid appendages the above derivation is readily extended to modify the matrices Cikj ' Cik ' and Ci with no further modifications to the succeeding development. In fact these matrices could be obtained by finite element analysis should the link shape be irregular as is often the case. Furthermore, the expression for B3i contains a term of order 8 2 which is by definition small and a candidate for later elimination.
L:
8 i1 Cilk
)]
w;
1=1
This final tenn contains tlle rigid body inertia terms. It should be noted that these terms are easily simplified if
mi
+ [2
L: k=l
mi
8ik (C ik +
L:
8 il Ci1k )]
1=1
'V;]} (35)
We note in the second line terms which are second order in the small variable 8jr and which can be ignored consistent with the assumption that the deflections are small. Noting the recurrance of certain terms above, it is convenient to define the following: m.1
Dik = Cik
+
L:
8 il Cllk
(36)
1=1
Finally, much of the
m.
complexity of the integration of the modal shape products can be done off line one time for a given link structure.
I
Gj = Ci +
L:
k=1
8 jt ( Cik
+ Cr ). (37)
W. J. Book
10
When these definitions are substituted into equation 35 one obtains:
and twisting about the longitudinal Xi axis. Compression is not initially included as it is generall y much smaller. Along an incrementa11cngth d1J the elastic potentia! is
~ ( oK / oei) - aK laq.J = dt J
dV
n
. { -aW-i [ lr
2:
2
+
1
..
2:
:07i)2 + Iy( :OYi)2] V1J
V1J
(00Xi)2}
+GI
.
2:
c'lik Dik \Vi +2
. T]}
c'lik Dix Wi
(40)
01J
x
m1
T
k=l
2
1
J
m.1
= -Ld1J{ E [IJ
G. \V.
oq.
i=j
ei
" T
.
k=l
(38)
where
oxi' 0 'i' and 0zi are the rotations of the neutral axis ofthe beam at the point1J in the x,1 y.1, and Z1 directions, respectivel y. Since def1ections are small. these directions arc essenti;]l1y parallel or perpendicular to the neutral axis of the beam,
are
The partials of K with respect to c'ljf and Sjf
considerably more complex due to the fact that Bli' B2i' and Bli
are functions of the deflection variables. The techniques of
simplification
are
similar
but
with
the
addition
of
E = Young's modulus of elasticity of the material
simplification due to the fact that if A is any antisymmetric G = The shear modulus of the material
matrix, as it would be if it were the difference of a maUix and its transpose, and if W was a compatible matrix for
I = The polar area moment of inertia of the link x cross section about the neutral axis.
multiplication, then Tr{ W A WT} = O.
I I = the area moment of inertia of the link eross ~ z . section about the Yi and zi axes, respectively,
~ ( oK laS·Jf ) - oK / oc'lfJ = dt n
22:
Tr
i=j+1
~ k=l
With a truncated modal appro ximation for U1e beam
[ "T - - L G. W . + { -o\V. oc'ljf
1
deformation the angles 0 xi' 0 yi' and O'i are represented as
1
summations of modal coefficients times the def1ection
~ k=l
variables. The x rotation, for example is
mi
~
k=l
+
\V.
J
~
k=l
2:
(41) c'likOxik' k=l where 0 xik is the angle about the Xi axis corresponding to the
0Xi =
k·th mode of link i at the point 1J. When dV ci is integrated over the link the integration can be taken inside the modal (39)
summations of equation 41 and its corresponding y and
Z
components. The following definitions thell prove useful:
4. System Potential Energy
(42)
The potential energy of the system arises from two sources considered here: clastic deformation and gravity. In both cases they are included by first writing the potential energy contribution of a differential clement, integrating over
where
00 ·1_X_I 00 ·k d1/ G I (1J)_X_' x 01/ 01/
(43)
the length of the link, and then summing over all links.
4.1. Elastic Potential Energy
(44)
Consider a point on the i-th link. undergoing small deflections. First restrict the link to be of the slender beam type.
The clastic potential is accounted for to a good
approximation by bending about the transverse Yi and zi axes
(45)
II
Recursive Lagrangian Dynamics
Note that Kikl = Kilk and that for certain special cases the orthorgonality of the modal functions can eliminate many of
Note that fik is found in the top row ofC ik . It is the distance from the undefonned center of gravity to the center of gravity
the tenns in equations 43 , 44, and 45. The clastic potential for
when all 8 are zero except 8 ik ' which is onc. The total distance
the total system, V0 can then be written as
n
Ve=..!.
m,
L L L
2 i=1
to the center of gravity from 0i Uoint i) is multiplied by the
mi
k=1
(46'
8ik 8il Kikl ·
mass to give ri.
1=1
Upon
Note that the Ve is independent of qi' the joint v:lriables.
taking
the
partial derivatives required by
Lagrange's equations we find for the joint variables n
aV g
(47)
= _gT
aq.
L
~
r .
. . aq.
)
I=J
(53)
1
J
For the deflection variables, for 1 :s j For deflection variables
av
----2..
a8 jf
av_ _ c
, (48)
a8 jf
For j
:s n-I'
n T
L
= -g
(54)
i=j+1
=n av T ----I:...=-gWt.
The fonn of equation 48 is much more general than the initial
08
assumptions made regJrding the contributions to the clastic
(55)
n nf
nf
potential energy would Jllow . Compression strain energy, and link forms other than beams can be represented in this fonn. The values of the coefficients
Kjkf
can be detennined
5. Lagrange's Equations in Simulation Form At this juncture the components of the complete
analytically or numerically, eg. by finite element methods.
equations of motion in Lagrange's formulation, except for the external forcing tenns. have been evaluated in equations 38,
4.2. Gravity Potential Energy
For a differemiJI clement on the i-th link of length d7j
47, and 53 for the joint equations: and in equations 39, 48, 54 and 55 for deflection equations. '111e external forcing tenns
the gravity potential is
arc the generalized forces corresponding to the generalized (49)
coordinates: the joint and deflection variables in this case. The generalized force corresponding to joint variable qj is the joint
where the gravity vector g has the fonn gT
= [0 .gx
f\
torque
gy gzl.
For the deflection variables the corresponding
When integrated over the length of the beam and summed
generalized j(lrCe will be zero if the corresponding modal deflections or rotatiuns ha\ e no di splacement at those locations
over all beams, the gravity potential becomes
where external forces are applied. Thus it is 3>sunled for the
n
Vg = - gT
L
Wiri
(50)
present de\c!opment that the modal functions arc selected so that is the case. This is con venient for utilizing the results as
i=1
~t
well. All motion
where r,
= M., rn. +
I:
the joint
IS
desc ribed in tenns of the joint
variable. (This is not tfue in the approach taken by Sunada and
8ik tik
(51)
Du bow ski (29» The fonn of Lagrangc's equations will then be:
k=1 M, = the total mass of link i
r,;
Thc joint equation j
= [I
rxi 0 01. a vector to the center of gravity from joint i (undefonned)
d.
av
av
aq .
aq .
- (aK/aq) - aK/aq. + _ c + ----2.. dt
)
)
)
)
= F.)
(56)
I.
tix
=
!'
o
JI. [0 \1 Yi1 lik IT d7j .
(52)
The deflection cq uation j, f d
.
- (oK/a8 f dt
)
) -
aKlOo
)f
av
av
aO
(8)f
+ _ 0 + ----2.. = jf
0
(57)
W. J. Book
12
These equations are in the "inverse dynamic" fonn.
~
..
To
W,vJ = Wv.J''} AJ
convert them to the simulation fonn one must extract the
'2
~
A.
+ 2 W AJ + WJ' 1 U2Jq
(64)
coefficients of the second derivatives of the generalized (65)
coordinates to compose an inertia matrix for the system. The second and first derivatives together make up the derivative of the state vector. which can be used in one of the available
5.2. Inertia Coefficients
integration schemes, e.g. Runga-Kutta, to solve for the state as a functiolJ of time for given initial conditions and inputs Fr
To obtain the inertia coefficienl~ that multiply the second derivatives. substitute equations 63 and 60 into the relevant parts of the equations of motion, equations 38 and 39,
5.1. Kinematics Revisited
respectively.
The purpose of this section will be to extend the kinematics to separate the second derivatives of the joint
Collecting the terms and arranging them for
efficient computation requires the steps outlined in this section.
variables and deflection variables from the expressions for Wi and
\1.;..I O~her
occurrences of these derivatives are already
5.2.1.lnertia Coefficients of Joint Variables in the
explicit in the formulation as it exists.
Joint Equations
/\11 occurances of q in equation 38 are in the expression First consider the product of transformations which
W.and two alternative ways of expressing it.
make up
for
,,,.T. I
J
When these terms are isolated, a double summation
over the indices i and h exists. Interchange the order of the
I
summation as follows: n
(59)
n
L h=1 L
(58)
L
i=j
h=l
n
L
i=max(h,j)
TIle resulting coefficient for joint variable qh in the joint eq uation j is
Carrying through the derivatives one obtains
(66)
L
WI
(Wh_} Uh h Wiqh
where
h=1
n
mh
+
L
W h Mhk h\\'i ShX)
+ WVi '
(67)
(60)
k=l
i=max(h,j) Note that if one exchanges j and h and transposes inside the trace operation an identical expression is obtained.
For the corresponding expression for \Vi write
This
indicates the symmetry of the inertia matrix which is used to Wi
= A} El A2 F.2 ... Ah Eh .. , E j _} Ai
reduce the number of computations required. The expression
=
hWho} Ah W h
(61)
for ifh can be computed recursively as will be described later
=
W h Eh hWi
(62)
to further improve the efficiency of calculation.
~
L
W. I
5.2.2. Inertia Coefficients of the Deflection
h- .. Who} U h Wiqh ~
Variables in the Joint Equations
h=l i-I
+
The deflection variables appear both in the expression mh
L k=1 L
.. W h Mhk Wi c'lhk +
for
h~
W vi
(63)
W1 and explicitly in equ
into equation38. collect terms in
h=1
After substituting
Sjr and
exchange the order of
summations as follows The value of\\' , and VI
W, VI
can be calculated recursively as ..
shown in equations 15 and 10, respectively, for only eliminating terms involving qj and
Sjx'
Wiand 'Vi
The result is
by
n
i-I
L h=l L
i=j
n-l
=
L
h=]
\V IT
n
L
i=max(h+l,j)
Recursive Lagrangian Dynamics
The resulting coefficient of 8hk in joint equation j is Jihk . The terms to be included depend on the relative values of j and
Forj = h =1. .. n-l:
2Tr{ Mif lll>i M;~ + Cikf } .
'iDk =
h. The following hold for 1 ~ k ~ mho Forh
= n,j
Forh = n;j
= 1 ... n:
Jink = Hr {
(\Vi.] Ui )i\\'n Dnk W~};
(68)
13
(73)
= 1 ... n-1:
lifnk = 2 Tr{ Wi 1\1/"'n Dnk
w~ } .
(74)
Forj = 1 ... n-l ; h = j+l ... n-1: 'ilhk = 2Tr{l\1if [i4>h + (69)
for h
= 1 ... j-l, j
where forh
(70)
iw. G
ill>h--
I
hW T
.
(76)
I
5.2.4. Recursions in the Calculation of the Inertia
(71)
i= max(h + 1, j)
Coefficients
Since the inertia matrix is a square matrix it requires the calculation of
It can be shown that the inertia coefficient for the
n;
coefficient for the joint variable qi in the deflection equation h.k. This further extends the symmetry of the inertia matrix and reduces the computation necessary.
n
=n +
L
mj i=l The fact that the matrix is symmetrical reduces the number of nt
distinct terms 5.2.3. Inertia Coefficients of the Deflection
terms where nt is the total number of
variables:
denection variable c'lhk in the joint equation j is the same as the
to
nt(n t + 1)/2. which still has a second power
dependence. Thus while the inverse dynamics computation
Variables in the Deflection Equation
In a manner much the same as for the previous two types c()efficienl~.
I
i=maxU+1, h+l)
n
of
(75)
n
= 1 .. . n-l , j = I ... n:
iFh --
iWh Dhk]WD ·
Terms in the above defined for j = 1 ... n-1 ; h = 1 ... n-1 are :
= 2 ... n:
link = Hr { ( \\:i.j U ) iF h M~kW~ } ; i
l\1~k+
the inertia coefficients of the denection
complexity can be made linear in n t. simulation requires the inertia matrix with complexity dependent on n~. Since nt can be quite large for practical arms it is important to reduce the
variables in the dcnection equations are evaluated. Symmetry
coefficient of the squared term as much as possible.
of the coefficients can be shown such the coefficient of variable
further that due to their short or even zero length it is possible
h.k in equationj.fis the same as the coefficient ofvariablcj,fin
for some links to be essentially rigid. Anthropomorphic arms,
Substituting equation 63 into equation 39,
for example. have two links which are much longer than the
equation h.k.
Note
isolating the second derivatives of the deflection variables, and
others and tend to dominate the compliance.
interchanging the order of sum mat ions enables the inertia
possible that many of the terms derived above will not be
Thus it is
coefficients 10 be identified. Further simplification is based on
needed for these links. four of the six links in the
the identity that. for Jny three square matricies A. Band C
anthropomorphic example.
calculating the terms in the equations should not require these
Tr{ AB C} = Tr{C :\ B} = Tr{B CA}. Furthermore the rotation matrices in the transformation matrices are orthogonal so that R.I RTI = I, a 3x3 identity matrix. This coupled with the zero first row and column of Cikf results in an especially simple form for two of the four cases. The following hold for 1 ~ k ~ mh and 1 ~ f ~ mj' For j
= h = n:
lnfnk = 2 Tr{ C nkf }
Any recursive scheme for
(72)
calculations as a means to get tu needed tcmlS. if possible. Consider the calculation of equations 67. 71 . and 76. Several recursi\'e schemes could be arranged for the efficient calculation of these quantities. Equation 71 is only needed if the link corresponding to the variable. link h. is flexible . That is, if mh > O. Equation 76 is only needed if both the link of the variable and the link of the equation. lin k j. is also flexible. Thus we propose the following recursive scheme for calculating iFh • iF h• and i4>h' The following hold for 1 ~ k ~
14 ill ;
h
W. J. Book 1 :s f
:s illj .
Intialization:
nF'n =
G
R = dynamics from the joint equation j (equation J 56) -:xcluding second derivatives of the generalized coordinates
(77)
n
:s n:
For j > h
R == dynamics from the deflection equation jf
jFh = EjA/+lfh
f J (equation 57) excluding second derivatives of the generalil.ed coordinates
(78)
For j = h:
hF'h --
Gh +
(79)
arranged to form the proper equations in the order described above. This order has been selected because it results in the
Ifm h > 0 calculate: .
JFh
=
The elements of .1 have just been formulated and can be
.-
T
(80)
JF h + 1 Ah .
symmetrical appearance of J . The elements of R have not been explicitly given with the second derivatives removed.
Ifmh > 0 and mj > 0 calculate:
These arc gi ven helow with some recursions (81)
jct>h=Aj+]jFh ·
10
fJcilitate their
computation. (83)
5.3. Assembly of Final Simulation Equations
The complete simulation equ,lliuns have now been It remains to assemble them in final form and to
derived.
point out some remaining recursion relations that can be used to reduce the numher of calculations.
The the second
derivatives of tile joint and deflection variables arc desired on m
the "left hand side" of the equatiun as unknowns and the
+
remaining dynamic effecLS and the inputs arc desired on the "right hand side." To carry out tllis process completely one
2W n L... ~.5 nk c nkf k=l
]
WT} n
m
would take the inverse of the inertia matrix J and premultiply
~
the vector of other dynamic effects. This inverse can only be
(85)
8 nk K nkf + gT W n fnf
k=l
evaluated numerically because of its complexity. Thus for the present purposes tile equations will be considered complete in the following form:
J
l
= R.
(82)
+ 2WJ
where
t
J == Inertia matrix consisting of coefficients preliously defined in the order for multiplication appropriate for z z = the vector of generalized coordinates = [q] 8]] 012 00' 8]m] q2 8 2]'00 8 2m2 00. qh 0h1
00.
0hk
.00
8 hmh
00.
k=1
where
~ .5 nk 0nk) w~
(87)
k=1
onm)T
qh == the joint variable of the h-th joint
m
Qn = G n W~n+ 2 (
QJ == Gj W;J + 2 (
I: .5
Jk DJk )
\VT
k=l
8 hk = the deflection variable (amplitude) of the k-thmode of link h R = vector of remaining dynamics and external forcing terms = [R] RH R]2 ... R]m] R2 R 2] 00. R 2m2 00. Rj Rj] ... Rjf .. · R jmj ... Rnm 1T n
(88)
Pn =AIn r n +
t
k=l
(89)
Recursive Lagrangian Dynamics m.
Pj
= MJ rJ +
t
Djk EJk
+
(90)
Ej AJ+l P j + 1
15
the number of calculations is approximately: Number Q[multiplications:
k=l
6n~ m2 + 17.5 n r m 2 + ll8 n; m+ 74 n nrm +
6. Conclusions Two measures of success of the ahove model are its
137.5 nrm + 84 n2 + 86 n nr + 279 n+ 126 nr" 57
accuracy and its speed. The two are somewhat related in that accuracy of the flexible representation can be improved by increasing the number of modes used
to
Number ill additions:
represent the link
deflection at tl"le expense of calculation time.
The issue is
65 . nr2 m2 + 19 nrm 2 + 115.5 nr2 m+ 68 n nrm +
further complicated by the choice of mode shapes. range of motion considered. and the arm configuration. Furthermore, limited
information
is
available
in
the
literature
for
comparison. t\ simple comparison has been used in the past and
can
be
performed
fur
calculation
complexity.
2
123 nrm + 85 n + 80 n nr + 329 n+ III nr" 91 where: n = total number of joints nr = number of flexible links m = number of modes describing each flexible link
Hollerbach [14J compares several approaches to the inverse
The above approximatioll
assume, an
dynamics problem of rigid arms by di fferent authors.
complexity mcr two common typc-s of rotary joints, the same
Walker [31J gives a similar count for four approaches to the
number of modes on c;lch flexible link, a rigid last link and a
simulation problem. Sunada [29J has given computation times
flexible first link.
"average" joint
for a given manipulator, trajectory, and computer for his flexible simulation. Comparison to the calculation counts of rigid models are given for a rough comparison of speeds in this No attempt at a quantitative comparison of the
section.
If assumed mode shapes arc restricted so that the shape
-
. I"
functions in the x. y, and z directions arc onh(){wnal. onlv C will be non"zero.
This is a stronger requirement than the
orthogonality of the set of complete mode shapes. but would
accuracy is made.
often be realized with simple mode shapes. It has not been In the determination of the number of calculations from the equations computed.
it
choice must be made on how some terms are
detennincd if this would improve the combination of speed and accuracy.
J-lollerhach took the approach that the most
straightforward implementation of the equations should be
This calculation count can be roughly compared to rigid
Obvious
link results available in the literature mentioned above. For a
used.
The attitude here is quite different.
simplifications in the multiplication of matrices with known
12 degree of freedom
constant rows, the top row of a transformation matrix for
transformation matrix formulation requires 2.66 times as many
example. are assumed in these computations. '111e 4x4 matrix
multiplies as the Newton·Euler fonnulation. Walker's method
transformation was chosen for its conceptual convenience and the calculation count will not be intentionally penalized for
rigid problem the inverse 3x3
3 (his best) for simulation requires 4.491 multiplies.
For 6
joints, and two flexible links with 3 modes each the method of this paper requires approximately 12,009 multiplies. The ratio
that choice. Furthermore, certain products appear in multiple equatiolls and are ~,s umed to be say·ed when needed later.
of these simulation methods is 2.67, almost exactly the same as
Spc-cial purpose multiply routines arc Llsed whenever they can
for the inverse dynamic methods with the same number of
capitalize on the speci al structure of a given matrix. Finally, in
degrees of freedom.
the simulation form the calculations needed to invert the
would be much more accurate than adding 6 imaginary joints
inertia matrix are not included. and no consideration is given
to represent compliance, but one could expect to use 15
A modal representation of flexibility
to the calculations of the intcgration routine. The general form
imaginary joints and 6 real joints with Walker's method with
of the modal parameters are used however. This results in all
fewer multiplies than with the method of this paper.
combinations of modes hand k in the matrix C ihk to be computed and
used and hence
introduces a squared
dependence on the number of modes on each inertial coefficient of the deflection variables. With these assumptions
16
W. J. Book
Thus it seems that in order to be competitive with possible Newton-Euler, non-transfer matrix approaches, the simplification of the assumed mode shapes will have to be made. It is not clear that the conceptual convenience of the transformation matrix approach can be justified relative to vector dyadic approa<.:hes of Hughes [16J.
10. Denavil. J , and R.S . Hartenherg. "A Kinematic Notation for Lower-Pair Mechanisms Based on ,~1atriccs." ASIIf E Juumal of Applied Mechanics 22 (June 1955),215-221.
11. Derby. Stephen J. Kinematic Fiasto-Dynamic AnalysiS and ('omputer Graphio Simulation ofGelleral Purpose Robot Manipulators. PhD. Th .. Renssc1aer Polytechnic Institute, August 1981.
Unfortunately,
wmpuwtion counts are nut available for that work .
References l. Albus.l.S. "i\ New Approach to Manipulator Control: The Cerebellar Model Articulation Controller (CMAC)." ASME
Juufllal of Dynamic S)'stel1l. Measurement and COlltrol97 (1975).270- 277. 2. nejczy. A,K. Robot Arm Dynamics and Control. NASA Technical Memorandum 33-669. Jet Propulsion Laboratory, California Institute of Technology, February, 1974. 3. nook. Wayne J. Mod('ling. Design and Contrul of Flexible Man ipulalUr Arms. Ph.]). Th .. Massachusetts Institute of Technology, ])ept. of Mechanical Engineering, April 1974. 4. Book. Way ne 1.. 0, Mai 7.za-Neto. and D.E. Whitney. "Feedback Control of Two Beam. Two Joint Systems with Distributed Flexibility." ASM /:' Journal of DYflamic Systems. IIIPasurement. and Control 97G, 4 (December 1975). 5. Book. Waync 1., Mark Majette. and Kong Ma. The Distributed Systems Analysis Package (DSAP) and Its Applic
12. Dubowsky. S.. and T.N . Gardner. "Design and Analysis of Multi-Link Flexible Mechanisms with Multiple Clearance Connections." AS.H E Joumal of Engincering for Industry 99, 1 (February 1977). 13. Ho. J .Y.L.. W.w. Hooker. G . Margulies. and T.P. Winarske. Remote Manipulator System (lUvIS) Simulation, Vol. 1. Locl:.heed Palo Alto Research Laboratory, October, 1974. 14. Hollerbach, John M. "A Recursive Lagrangian Formulation of Manipulator Dynamics and a Comparative Study of Dynamics Formulation Complexity." IEEE Transactions on S)'stems. Afan. and Cybernetics SMC-lO, 11 (November 1980). 730-736. 15. Hughes, P.e. Dynamics of a flexible Manipulator Arm for the Space Shuttle. AASI AIAA Astrodynami<.:s Conference, American Insitute of Aeronautics and Astronautics, September. 1977. Jackson Lake Lodge, Grand Teton National Park, Wyoming, USA 16. Hughes. P.e. "Dynamics ofa Chain ofl-lexib1c Bodies."
Journal of Astronautical Sciences 27. 4 (October-December 1979), 359-380. 17. Huston, R.L.. e.E. Passerello. "Multibody Structural Dynamics Including Translation Between Bodies." Computers and Structures f2, 5 (November 1980),713-720. 18. Kelly , Fred A. and Ronald I .. Huston. Modelling of Flexibility Effects in Robot Anns. Proceedings of the 1981 loint Automatic Control Conference. American Automatic Control Council. lune, 1981. Paper WP-2C 19. Liegois. A .. W. Khalil. J.M . Dumas, and M. Denaud. M
21. Luh. J.Y.S .. M.W. Walker and R.P.e. Paul. "On-Line Computational Scheme for Mechanical Manipulators."
ASM E ./ournal of DYl1alllic SystflllS, lIfea.lurcl/lcnt. and Control f02 (June 1980), 69-76. 22. Maizza-Neto.Octavio. Modal AnalySiS and COlltrol of Flexible Manipulator A rillS. Ph,]), Th" Dept. of Mechanical Engineering. Massachusetts Institute of Technology, 1974.
Recursive Lagrangian Dynamics 23. Mirro, John. Automatic Fccdback Control of a Vibrating Beam. MJSler 1'h., Depl. of Mcchanical Engincering, Massachusctl\ Institute of Technology, August 1972, Also published as Charles Stark Draper Laboratory Report No. 1'-571
17
29. Sunada, W. and S. Dubawsky. "The Application of Finite Element Methods to the Dynamic Analysis of Spatial and Coplanar l.:nkage Systems." ASM 1~· .Iournal of Mechanical Designl()3 (July 1981),643-651.
30. Thomas. M. and D . Tesar. Dyanamic Modcling of Serial 24. Nguyen. P,K, and P.e. Hughcs. Finite-Elcment Analysis ofCTS-Likc Flexible Spacecraft. Tcch. Rept. 205, Institute for Aerospace Studies, University of Toronto, June, 1976.
Manipulator Arms. Proceedings of the 1981 Joint Automatic Control Conference. American Automatic Control Conference. June. 19X1. also to appear in ASl'''/~'Journal of
Dynamic .'·')'s/cII7s Measurcmen/ and ('ol11rol 25. Raibert. M.H . and B.K.P. Horn. "Manipulator Control Using the Configuration Space Method." The Industrial Robot 5,2 (June 1978),69-73. 26. Silver, William M. On the Equivalence of Lagrangian and Newton-Euler Dynamics for Manipulators. Proc.cedings of the 1981 Joint Automatic Control Conference. The American Automatic Control Council, June. 1981. Paper no. TA-2A 27. Stephanenko. Y. and M. Vukobratovic. "Dynamics of Articulated Open-Chain Active Mechanisms." Mathematical Bi05ciences 28(1976).137-170. 28. Sturges. Robert. Te1coperatar Ann Design Program (TOAD). Tech. Rept. E-2746, The Charles Stark Draper Laboratory, February, 1973.
31. Walker. M.W. and D.E. Orin. "Efficient Dynamic Computer Sirnubtion of Robmic~\'kchani s m, ." ASM E Journal of Dyl/amic S~\ls!ems. Measurc/IICI/! and CUI//ruII04, 3 (September 1982). (to appear) 32. Whitney. D.E .• W.J. flook. and P.M. Lynch. Design and Control Considerations filr Industrial and Sp;lce l\iJnipulators. Proceedings of the 1974 Joint Automatic Control Conference, The American Automatic Control Council , 1974, pp. 591-598. Austin, TX 33. Winfrey. R.e. "Dynamic Analysis of Elastic Mechanisms by Reduction of Coordinates." ASM E Juurnal of Engineering for Indus!ry 94, 1 (May 1972).