Int. J. Mechanisms and Robotic Systems, Vol. 1, No.

1, 2013 15

A contribution for developing more efficient dynamic

modelling algorithms of parallel robots

Renato Maia Matarazzo Orsino*

Department of Mechanical Engineering,
Polytechnic School,
University of São Paulo,
Avenue Professor Melo Moraes 2231,
Sao Paulo, SP 05508-900, Brazil
E-mail: [email protected]
E-mail: [email protected]
*Corresponding author

Tarcisio Antonio Hess Coelho

Department of Mechatronics and
Mechanical Systems Engineering,
Polytechnic School,
University of São Paulo,
Avenue Professor Melo Moraes 2231,
Sao Paulo, SP 05508-900, Brazil
E-mail: [email protected]

Abstract: Parallel kinematic structures are considered very adequate

architectures for positioning and orienting the tools of robotic mechanisms.
However, developing dynamic models for this kind of systems is sometimes a
difficult task. In fact, the direct application of traditional methods of robotics,
for modelling and analysing such systems, usually does not lead to efficient and
systematic algorithms. This work addresses this issue: to present a modular
approach to generate the dynamic model and through some convenient
modifications, how we can make these methods more applicable to parallel
structures as well. Kane’s formulation to obtain the dynamic equations is
shown to be one of the easiest ways to deal with redundant coordinates and
kinematic constraints, so that a suitable choice of a set of coordinates allows the
remaining of the modelling procedure to be computer aided. The advantages of
this approach are discussed in the modelling of a 3-dof parallel asymmetric
mechanism.

Keywords: parallel mechanisms; manipulator; kinematics; dynamics;

homogeneous transformations; analytical mechanics; Kane; numerical
simulation.

Reference to this paper should be made as follows: Orsino, R.M.M. and

Coelho, T.A.H. (2013) ‘A contribution for developing more efficient dynamic
modelling algorithms of parallel robots’, Int. J. Mechanisms and Robotic
Systems, Vol. 1, No. 1, pp.15–34.

Copyright © 2013 Inderscience Enterprises Ltd.

16 R.M.M. Orsino and T.A.H. Coelho

Biographical notes: Renato Maia Matarazzo Orsino is currently studying for

his Masters degree in Mechanical Engineering from Polytechnic School of
Engineering of the University of São Paulo (EPUSP), Brazil. He received his
BEng in Mechanical Engineering from the Polytechnic School of Engineering
at the University of São Paulo (EPUSP), Brazil in 2011. His main research
areas are multibody dynamics and parallel robots.

Tarcisio Antonio Hess Coelho is an Associate Professor of the University of

São Paulo since 2005. He received his BEng, MEng and PhD in Mechanical
Engineering from the Polytechnic School of Engineering of the University of
São Paulo (EPUSP), Brazil in 1987, 1990 and 1997, respectively. From
2001–2002, he was a Post-Doctoral Fellow at Stanford University, USA. He is
also a member of two IFToMM Technical Committees (Mechatronics and
Robotics, Multibody Dynamics). His main research areas are parallel robots,
education and biomechanics.

1 Introduction

Parallel kinematic structures have shown a lot of advantages in many applications

compared to the traditional serial structures. These advantages are mainly related to the
fact that multi-degree of freedom parallel structures can be built in such a way that all the
actuators remain attached to a fixed base, avoiding the motion of heavy parts. In fact,
such mechanisms achieve larger velocities and accelerations without consuming a large
amount of energy. Nevertheless, the dynamics of a parallel structure is so complex,
compared to a serial equivalent one, that it is a quite challenging task to design a suitable
control system for it without a detailed dynamic model that considers relevant effects
affecting its own behaviour.
The complexity of the dynamics of parallel structures leads to some difficulties in the
analysis of such systems requiring the use of some methodology to obtain a mathematical
model which will be able to describe this dynamics in an adequate way. The aim of this
paper is to describe how some renowned methodologies of kinematics and dynamics can
be optimised in order to make more efficient mathematical models of parallel
mechanisms. These models may be deduced with the aid of a computational tool
(avoiding errors and saving time) and may be able to produce the most commonly desired
results in an analysis (such as inverse kinematics and inverse dynamics), being suitable to
make efficient numerical simulation algorithms. The procedures for analysis and
modelling of parallel mechanisms will be described and discussed, showing the most
appropriate way of integrating these methodologies to study kinematics and dynamics of
these systems.
First of all the kinematics of the system has to be well understood in such a way that
the motions of each part of a mechanism can be well described by an adequate selection
of a set of coordinates. The most important coordinates of the analysis (the ones which
the analysis will be designed to determine) are the coordinates related to the motions of
the actuators of the system and to the motions of its end-effector. Actually, the motions of
all other parts of such kinematic structure are generally of no interest and the definition of
coordinates to describe this kind of motions will be acceptable only if it leads to some
simplification in the mathematical model designed. After the identification and
description of the motions by an adequate selection of a set of coordinates, it is necessary
 A contribution for developing more efficient dynamic modelling algorithms 17

to choose a methodology to describe the relations among the end-effector coordinates, the
actuators coordinates and the auxiliary coordinates. A general and renowned
methodology to achieve this goal is the homogeneous transformations method (Craig,
1989). By this analytical mechanics method, coordinates of points described in different
coordinate systems can be related, which enables the development of the constraint
equations (with these equations, given an independent set of coordinates, like the
coordinates of the end-effector, for example, all the other coordinates of the model can be
obtained by the solution of these equations). With the constraint equations the complete
kinematic analysis of the system can be carried out.
With a good description of the kinematics of the system, the development of the
dynamic model is more properly done by methods of analytical mechanics such as
Kane’s and Lagrange’s (Kane and Levinson, 1985; Leech, 1971; Lanczos, 1983). In a
comparison with the classic approach of Newton-Euler, the main advantages of those
methods are that inertia forces can be methodically obtained by the kinematic analysis of
the system and the included active forces are only those which actually matters and not
the constraint ones. The description of such active forces may sometimes require the
definition of more coordinates in order to simplify this process. This can be done by
adding new constraint equations to the original model.
Finally, as proper values of the parameters of the system are suitable estimated,
numerical simulations can be performed, making it possible to assess the dynamic
behaviour of the system. Actually, all the methodology described here is applied to a
3-dof parallel asymmetric mechanism, showing the advantages of the methods selected to
develop each part of the analysis of such system.
In Section 2 some discussion about the methods applied is done; in Section 3 a simple
example of application of the methods to obtain the dynamic equations of a slider-crank
mechanism (with a non-rigid connecting rod) is presented in order to make the reader
familiar with the purposed approach; in Section 4 the 2RSS + PPaP mechanism is
presented, its dynamic model is detailed and some numerical simulations results are
shown and discussed, leading finally to the conclusions in Section 5.

2 Methodological background

2.1 Stages of the modelling process

There are several ways to make kinematic and dynamic analysis of a mechanical system,
but there are some characteristics that may hold in order to make an efficient model of
this system: the kinematic and dynamic methods must be easily integrated, most of the
modelling steps must be computer aided (saving time and avoiding human errors) and the
numerical simulation algorithms required to make these analysis must be as simple and
reliable as possible.
The methods of kinematic and dynamic analysis will be detailed on Sections 2.2 and
2.3 respectively. It is possible to use the methodologies for all kinds of mechanical
system, but they are optimised to be used to model complex systems such as parallel
mechanisms (it is an unnecessary effort to use sophisticated methodologies to simple
mechanical systems, where simpler analysis are demanded). Basically, the stages of the
modelling process presented here are shown in Figure 1 and in Table 1. This flowchart
indicates not only all the steps involved in the modelling process, but also their
18 R.M.M. Orsino and T.A.H. Coelho

interaction (the numbers of the stages indicate the order of the first occurrences of each of
them in the overall process and the arrows indicate the precedence in their interaction). It
is worth noting that of all stages, only the definition of coordinates is iterative, and also it
is not computer aided. As it will be discussed in subsequent sections, it is one of the most
important stages of the modelling processes, being the way to deal mathematically with
the motions of a mechanical system, so that an adequate choice of a set of coordinates is
determinant to obtain a model with the desired characteristics cited at the beginning of
this section.

Figure 1 Flowchart showing the steps of modelling a mechanical system, integrating the
topology, kinematic and dynamic analysis

Table 1 Description of the stages of the modelling process of a mechanical system according
to the flowchart shown in Figure 1

Stage Description Group

1 Chains and links identification Topology
2 Mobility analysis
3 Coordinates and generalised speeds definition
4 Constraint equations Kinematics
5 Velocities, angular velocities and its derivatives
6 Partial velocities and partial angular velocities
7 Inertia forces and inertia torques Dynamics
8 Generalised inertia forces
9 External forces/torques and additional effects
(gravity, friction, elasticity)
10 Generalised active forces
11 Dynamic equations of motion

The main objective of the procedure in Figure 1 is to obtain the constraint equations
(or kinematic equations), which corresponds to stage 4, and the dynamic equations,
which is the stage 11. There is then a need to develop a routine to simulate the system
constituted by these equations, which may involve techniques to deal with algebraic and
differential non-linear equations. Again, the choice of an adequate set of coordinates may
be responsible for avoiding the occurrence of non-linear functions that can present
singular points where the numerical algorithms used may fail.
 A contribution for developing more efficient dynamic modelling algorithms 19

2.2 Topology analysis and kinematic modelling: constraint equations and the
homogeneous transformations method
The first stage of the modelling process of a parallel robot is the identification of chains
and links. It is necessary to consider that such systems generally have two most relevant
parts which are the base of the robotic mechanism and its end-effector, whose motion
relative to the base is the main object of study in a modelling process. Then, it is
necessary to analyse the restrictions that each chain of mechanism impose to the motion
of the end-effector relative to the base. Such approach is a mobility analysis based on the
method of the Lie group of rigid body displacements (Hervé, 1999). After this analysis it
is possible to determine the number of degrees of freedom of the system, which is the
maximum number of independent virtual displacements that can occur in the
neighbourhood of any possible position of this system. This kind of analysis also
determines the number of degrees of freedom of the end-effector and the number of
inputs (external forces and torques) which are normally the same, so that this method
leads to the following natural definitions of coordinates:
• a first set consisting of coordinates that describe all the possible motions of the
end-effector of the system
• a second set consisting of coordinates that describe the motions imposed by the
actuators of the system
• a third set consisting of auxiliary coordinates defined to describe motions of links of
the mechanism’s chains which cannot be described trivially by the coordinates of the
two first sets (these variables have to be carefully chosen, otherwise the increase of
the number of variables might not lead to a significant benefit).
It is worth noting that this method is similar to one which has been very successful in
both kinematic and dynamic analysis of complex mechanical systems such as kinematic
structures presented by Tsai (1999) and Merlet (2006). In this methodology, coordinates
are defined to describe the motion of the main link of the mechanism and to describe the
motions done by the systems’ actuators. Calling z the vector constituted only by the
coordinates of the motions of the actuators and x the vector constituted by the coordinates
of the end-effector, Tsai (1999) and Merlet (2006) purpose that the constraint equations
can be written in the following form: J z z = J x x, where Jz and Jx are called the system’s
Jacobians with respect to the actuators coordinates vector z and with respect to the
end-effector coordinates vector x. Also noting that this relation is equally valid to
infinitesimal variations of the configuration of the system (i.e., Jzδz = Jxδx), this approach
leads to two independent sets of coordinates (z and x) which can be easily related and
interchanged, enabling an easy description of the main motions of the system (the most
relevant motions of the system are those of the end-effector and this approach leads
to an easy computation of the motions of the actuators related to any motion of the
end-effector, which is called an inverse kinematic analysis of the system). When
designing a dynamic model of such system, however, not only the description of the
motions of the main link and of the actuators are required, but also the motions of all
links whose inertia can have a significant influence in the system’s dynamics. Depending
on the complexity of the chains which constitute the kinematic structure studied, using
only the coordinates of z and x may lead to some complex descriptions, involving
20 R.M.M. Orsino and T.A.H. Coelho

difficult mathematical expressions and leading to a model numerically inefficient and

unstable. That is why it is recommended the use of a third set of coordinates in order to
describe some complex motions of links of the chains which cannot be trivially described
by the coordinates in z and x. Sometimes, when the dynamic model is being developed it
may be recommended the definition of more coordinates (a fourth set) that will be used to
make simpler descriptions of some of the active forces.
All the coordinates defined in the analysis can be ordered in a single column vector
q* consisting of n* coordinates qi (the * will be used to denote that this may not be the
final version of the vector of coordinates which will be used in the analysis, emphasising
the iterative character of the choice of coordinates). Considering that the whole system
has p degrees of freedom, among all the δqi (i = 1, … n* ), at most p can be considered
independent and sufficient to describe any virtual displacements in the neighbourhood of
any possible configuration of the system (if the system is in a singular configuration, less
than p independent virtual displacements are possible). Thus, it is necessary to obtain
n* − p constraint equations so that the model described by the coordinates defined lead
to motions consistent with the physical constraints of the system. Basically, the
constraints are classified as holonomic or non-holonomic. A constraint is said to be
holonomic if it can be expressed as a function among the coordinates of the system, not
involving any of its time derivatives, i.e., this constraint can be well described by an
equation of the form:

(
fi q1 ,… , qn* = 0 ) (1)

This text will make the option of using a more simple notation, referring to such equation
as fi (q* ) = 0 and to the full set of constraint equations as f (q* ) = 0 (here, all the
equations fi are treated as components of a column vector f). Any relation among the
system’s coordinates that cannot be described by this kind of function characterises a
non-holonomic constraint. A remarkable type is those Kane and Levinson (1985) call
simple non-holonomic constraint that can be described in the following form:
n*

∑ J ( q ,…, q ) dq
j =1
i, j 1 n* j =0 (2)

This differential form is not exact, otherwise it would be possible to integrate

this expression, finding a result of the form fi (q* ) = 0, which would characterise
a holonomic constraint. However, the most important methodologies of analytical
mechanics (Kane and Levinson, 1985) involve only the differential form of the constraint
equations, making them applicable to mechanical systems whose constraints are
holonomic or simple non-holonomic. In the most common kinematic structures, all the
constraints are holonomic and all the constraint equations can be obtained by relations
involving the position of relevant points of the main links of this mechanism.
Nevertheless, these equations are commonly non-linear and hardly ever will lead to a
single solution for the redundant coordinates given a set of independent ones.
Considering that these equations are normally solved by numerical methods, there is a
need to guide the convergence of the solution to the real configuration of the system.
Although the exact solution of the equations is not previously known, an interval of
 A contribution for developing more efficient dynamic modelling algorithms 21

values where each variable must be is well known, and an initial estimation can be done
by selecting a value in this interval. This, however, is not enough to ensure convergence
of the solution algorithm. Here comes the need for an adequate definition of redundant
coordinates. Noting that a coordinate may be related to the position and the motion of a
particular link of the system studied, the definition of such coordinate (even if it is not
independent) may define uniquely a position or a motion of this link. That is why
defining additional coordinates and making good initial estimates may contribute to a
proper convergence of the numerical solution algorithm.
In order to obtain the holonomic constraint equations of a kinematic structure, one
adequate method is the use of the homogeneous transformations (Craig, 1989). To
introduce how it is used, consider that two rigid bodies A and B of the system studied. In
the most general case, B is able to describe translation (at most three degrees of freedom)
and rotation (at most three degrees of freedom) relative to A. The description of the
translation involve at most three independent functions of the coordinates which may
represent three components of the position vector pOAOB of a point OB of B relative to a
point OA of A in an adequate vector basis. The description of the rotation may involve at
most three independent parameters which are called the Euler angles (which have to be
functions of the coordinates). Considering two vector basis {a1, a2, a3} and {b1, b2, b3}
fixed in A and B respectively, the matrix of the directional cosines ARB whose element of
the ith row and jth column is ai ⋅ bj (i = 1, 2, 3 and j = 1, 2, 3) can be written as a function
of the Euler angles. Consider, for instance, a vector u whose components in the basis {b1,
b2, b3} are given by the 3 × 1 matrix Bu. Its components on the basis {a1, a2, a3}, can be
obtained by the relation Au = ARBBu. Now, consider a point P and call A P = ApOA P ,
B
P = B pOB P and A
d= A OAOB
p ; A p OA P = A OAOB
p + ApOB P and A
pOB P = A
R BB pOB P ,
so that:

⎡ A P⎤ ⎡ A RB A
d⎤ ⎡ B P⎤ A A
⎢ ⎥=⎢ ⎥⎢ ⎥ or P= T BB P (3)
⎣ 1 ⎦ ⎣⎢ 01×3 1 ⎦⎥ ⎣ 1 ⎦

where

A
⎡ A RB A
d⎤ ⎡ A P⎤
TB = ⎢ ⎥ and A
P= ⎢ ⎥ (4)
⎣⎢ 01×3 1 ⎦⎥ ⎣ 1 ⎦
The matrix ATB is called the homogeneous transformation matrix between the coordinate
systems attached to A and B. The main property of these matrices is that ATB = ATCCTB
where C is any reference frame with a well defined coordinate system. This leads to
A BB A
T T = I, so that ATB = (BTA)–1.
The use of homogeneous transformations in order to obtain holonomic constraint
equations can be done by the following procedure: consider that the rigid bodies A and B
are connected by two kinematic chains K and L (which can be real or fictitious) so that
the links K1 and L1 constitute kinematic joints with A and the links Kk and Ll (where k and
l are the total number of links in K and L) constitute kinematic joints with B. Considering
the properties of the homogeneous transformations, and using the system’s coordinates to
describe each of these matrices (which is also a way of guiding the coordinates selection)
it can be stated that:
22 R.M.M. Orsino and T.A.H. Coelho

A
TB = A
T K1 …K k T B = A
T L1 …Ll T B (5)
These identities lead to the holonomic constraint equations among the coordinates which
constitutes the vector q* . Some of these equations will lead to trivial relations among
coordinates (i.e., two of the coordinates in q* are permanently equal or proportional)
which allows a reduction in the set of coordinates. The ones to be eliminated are those
considered as auxiliary (it is preferable to keep the coordinates defined to describe the
motions of the end-effector and the motions imposed by the actuators because they are
more relevant to the analysis). The coordinates vector obtained after this simplification
process will be called q (consisting of n variables qi), the trivial holonomic constraint
equations can be eliminated and the remaining holonomic constraint equations will have
∂f
∑
n
the form f(q) = 0. Its differential form is dq j = 0 which can be joined with the
j =1 ∂q
j

eventually present simple non-holonomic constraint equations, that reflect motion

constraints such as rolling without slipping, as discussed in Kane and Levinson (1985),
leading to a differential form (where J is a n – p × n Jacobian matrix):
J (q)dq = 0 and J (q)q = 0 (6)

The homogeneous transformations method can also be used to achieve the second
objective of the kinematic analysis which is the description of positions, velocities and
accelerations of relevant points and the angular velocities and angular accelerations of the
links of the system. Suppose that A is attached to an inertial reference frame, P is a point
in B (so that BP is a constant vector). Using the homogeneous transformation
A
P = AT BB P , A v P = A P and A a P = A v P . Also, adopting an identity presented in Kane
and Levinson (1985), the angular velocity of B relative to A can be determined by the
equation:

A ⎛
ωB = ⎜ B R A
⎝ dt
(
d A BB
) ⎞ ⎛
R b 2 ⋅ B b3 ⎟ b1 + ⎜ B R A
⎠ ⎝
d
dt
( A
) ⎞
R BB b3 ⋅ B b1 ⎟ b 2
⎠
(7)
⎛
+ ⎜ B RA
⎝ dt
(
d A BB
) ⎞
R b1 ⋅ B b 2 ⎟ b3
⎠
The corresponding angular acceleration AαB is the vector whose components in the vector
basis {b1, b2, b3} are the time derivatives of components of the angular velocity vector
A B
ω in the same basis.
Finally, its worth discussing that some modelling approaches such as Newton-Euler’s,
Kane’s (Kane and Levinson, 1985), Hamilton’s (Lanczos, 1983) or Poincaré’s (Udwadia
and Phohomsiri, 2007) make use of kinematic variables which replace the time
derivatives of the system’s coordinates, simplifying even more the models developed.
These variables can emerge from the physical description of the system by a well defined
relation as in Hamilton’s approach (where the generalised momenta are defined by partial
derivatives of the system’s Lagrangean function with respect to the time derivatives of
the system’s coordinates) and Newton-Euler’s approach (where the components of the
angular velocity vector of a rigid body are obtained as affine transformations involving
the time derivatives of the Euler angles which describe the orientation of this body) or
can be freely defined during the developing of the model as in Kane’s and Poincaré’s
 A contribution for developing more efficient dynamic modelling algorithms 23

approach. In any case, these variables constitute a column vector u of n components

(the same number of components as the coordinates vector q), and there exists an
invertible affine transformation between the variables in u and q that can be defined by
explicit (and known) functions Xi,j(q, t) and yi(q, t) such that:
n
ui = ∑X
j =1
i , j ( q, t )q j + yi (q, t ) (i = 1, … , n) (8)

This transformation can be written in matrix form as:

u = X ( q, t ) q + y ( q, t ) and q = V (q, t )u + w(q, t ) (9)

where V = X–1 and w = –X–1y. The reader will find several examples of ways to define
these new kinematic variables ui (i = 1, … , n) as affine transformations of the time
derivatives of the coordinates in Kane’s approach, where they are referred as generalised
speeds, in Kane and Levinson (1985) and Mitiguy and Kane (1996).
As the variables in u replace the time derivatives of the system’s coordinates, it is
necessary to express the constraint equations (6) in terms of the vector u. This can be
easily done using the transformation (9) leading to:
J (q, t )u + w(q, t ) = 0 (10)

where J = JV and w = Jw. These equations can be made explicit in u by means of an

algebraic artifice. The matrix J has rank (n – p) except in a singular configuration of the
system (because J has rank (n – p) and V is invertible). Thus, it can be stated that
it is possible to choose (n – p) columns of the matrix J constituting an invertible
(n – p) × (n – p) matrix Jˆ. Calling u the column vector constituted by the (n – p)
components of u with the same indexes as the columns which constitute the matrix Ĵ , J
the (n – p) × p matrix constituted by the other columns of J and defining u analogously
to u, the equation (10) can be written as:

Jˆ (q, t )u + J (q, t )u + w(q, t ) = 0 and (

u = − Jˆ −1 J u + w ) (11)

Considering also that u = I p× p u (I denoting the identity matrix), and that every variable
in u is either in u or in u, it can be stated that the constraint equations can be written in
explicit matricial form as:
u = B (q, t )u + d (q, t ) (12)

By equation (12) it can be noted that all the kinematic variables in u can be expressed in
terms of only p independent ones among then (that constitute the column vector u ). This
interpretation leads to physical a way of finding an invertible matrix Ĵ constituted by
(n – p) columns of the (n – p) × n of the matrix J : it is sufficient for the variables chosen
to constitute the column vector u to be independent, i.e., they must represent
independent motions of the system. The most common choice of an independent set of
variables includes the ones used to describe the motions of the end-effector, or
24 R.M.M. Orsino and T.A.H. Coelho

alternatively the ones used to described the motions imposed by the actuators
(more variables must be included in u if the number of degrees of freedom of the whole
system p is greater than the number of degrees of freedom of the end-effector).

2.3 Dynamic modelling: Kane’s method

Among the methods of analytical mechanics which are able to employ redundant
coordinates, there is one very suitable proposed by Kane and Levinson (1985). Kane’s
method is based on the approach previously developed by Gibbs and Appell, who gave a
modified formulation of the Gaussian principle of least constraint, which is a variant form
of d’Alembert principle, as shown in Lanczos (1983), that is suitable for a formulation of
dynamical equations based on kinematic variables (also called generalised speeds).
According to the discussion presented in the previous subsection, in a mechanical system
with p degrees of freedom where n generalised coordinates are defined, all the
generalised speeds u1, …, un can be expressed as functions of u1, …, up, q1, …, qn, t only,
according to equation (12). Moreover, q = B (q, t )u + d (q, t ) (where B = VB and
d = Vd + w ), so that, all the expressions of velocities and angular velocities can be
expressed as functions of u1 ,… , u p , q1 ,… , qn , t only. Kane and Levinson (1985) define
the ith partial velocity of a point P in a reference frame N (which in this text will be an
inertial reference frame) and the ith partial angular velocity of a body B in a reference
frame N respectively as:

N
N P
v iP = ∂ v and
N B
ωi =
∂ ( N
ωB ) (i = 1, … , p) (13)
∂u i ∂u i

After defining the partial velocities and partial angular velocities, the generalised active
forces and the generalised inertia forces can be determined. The former are defined as the
sum of the scalar products of the external applied forces and a given partial velocity of
the point of application of this force along with the scalar products of the external
applied torques and the respective partial angular velocity of the body where
this torque is applied, that is the ith generalised active force of the system can be
computed as:

Fi = ∑ (F ⋅
j
j
N P
vi j + ) ∑ (T ⋅
j
j
N Bj
ωi ) (14)

Analogously, the ith generalised inertia force can be defined as:

Fi* = ∑ ( −m
j
j
N P N
a j ⋅ vi
j P
) (15)

For a rigid body B the system of inertia forces of all its particles is equivalent to a single
force R* and a single moment T*, with:
*
R* = − m N a B and T* = − N α B ⋅ I* − N
ω B × I* ⋅ N
ωB (16)
 A contribution for developing more efficient dynamic modelling algorithms 25

*
where m is the total mass of B, I* is its central inertia tensor, N a B is the acceleration of
the centre of mass B* of B and NωB and NαB are its angular velocity and angular
acceleration, so that the ith generalised inertia force due to B can be calculated as:

(F ) B
*
N B N
i
* = R* ⋅ v i + T* ⋅ ωi (17)
B

The Kane’s dynamical equations of the system are given by:

Fi + Fi* = 0 for i = 1, … , p (18)

Noting that all the velocities and angular velocities can be written as functions of q, u and
t only, with all the components linear in u, the accelerations and angular accelerations as
its time derivatives are functions of q, u, u and t (remembering that u replaces q ), and
by the chain rule (from differential calculus), it can be stated that all its components are
linear in u and the terms where some u j is present are independent of u; moreover, as
partial derivatives of the velocities and angular velocities with respect to u, the partial
velocities and partial angular velocities are functions of q and t only. Taking all this into
consideration, it is possible to state that all the generalised inertia forces can be expressed

∑
n
as Fi* = − M i , j (q, t )u j − gi* (q, u , t ) while all the generalised active forces can be
j =1

expressed as Fi = − gi (q, u,τ , t ), where τ is a vector which expresses all external input to
the system, such as torques and forces provided by the system’s actuators. The dynamical
equations of the system can be expressed as:

M (q, t )u + g * (q, u , t ) + g (q, u ,τ , t ) = 0 (19)

If in a particular system all the inputs in the vector τ are control inputs, if the number of
inputs is equal to the number p of degrees of freedom of the system and if the system is
not in a singular configuration (when it is uncontrollable), then the equation (19) can be
used to determine numerically the inputs associated with a specified motion (when q and
u are prescribed time functions that respect the kinematic constraints) in a process called
inverse dynamics. However, sometimes the numerical integration of such equations is
needed for a given set of inputs, in order to determine the system motion. In this case,
which is called direct dynamics, the integration of the dynamical equations needs to occur
simultaneously with the integration of the constraint equations. This is done using the
time derivative of the equation (10) which is: Ju = −( Ju + w). The equations to be
numerically integrated can be expressed as:
q = Vu + w
−1
⎡M ⎤ ⎡ g* + g ⎤ (20)
u = −⎢ ⎥ ⎢ ⎥
⎣J ⎦ ⎢⎣ Ju + w ⎥⎦

Finally, it is important to mention that the main advantage of Kane’s approach compared
to Lagrange’s (Leech, 1971; Lanczos, 1983), which is one of the most used methods of
analytical mechanics, is that the number of dynamical equations is equal the number of
degrees of freedom of the system and not equal the number of coordinates defined
26 R.M.M. Orsino and T.A.H. Coelho

(which requires the use of the method of undetermined multipliers to obtain the
Lagrange’s dynamical equations). However, Kane’s method requires the evaluation of
accelerations and angular accelerations, which may be cumbersome, while Lagrange’s
approach requires only the evaluation of the kinetic energy of the system such that
knowing the expressions of velocities and angular velocities is enough.

3 Example of the application of the method

In order to illustrate the complete modelling process described on the previous section, a
simple example of application will be made using a slider-crank mechanism with a
non-rigid connecting rod as the one shown on Figure 2 (this connecting rod will be
considered as a linear spring with stiffness k and natural length a3 that can suffer
longitudinal deformations in tension and compression). The crank will be considered
actuated by a torque τ (with the imposed rotational motion described by the coordinate
q2) and the piston will be considered as the end-effector. A mobility analysis allows
noting that although the end-effector has only one degree of freedom, which is described
by the translational coordinate q1, the whole system has two degrees of freedom, because
an elongation of the connecting rod (indicated by the auxiliary coordinate q4) can happen
independently of the motion of the piston. The angle q3 indicating the orientation of the
connecting rod is also chosen as an auxiliary coordinate.

Figure 2 Representation of a slider-crank mechanism with a non-rigid connecting rod and its
coordinates and geometric parameters definitions

Having p = 2 degrees of freedom and n = 4 coordinates, there is a need to define n – p = 2

constraint equations. Calling si = sin(qi) and ci = cos(qi) (i = 2, 3), the following constraint
equations immediately arise from the system’s geometry:
q1 − a2 c2 − ( a3 + q4 ) c3 = 0 and a2 s 2 − ( a3 + q4 ) s3 = 0 (21)

Considering ui = qi for i = 1, 2, 3, 4 and assuming that u1 = u1 and that u 2 = u2

according to equations (6), (10), (11) and (12) for the given system the expressions of the
matrices J and B are the following (noting that the column vector d is null in this case):

⎡ 1 0 ⎤
− c3 ⎤ ⎢ 0 ⎥
⎡1 a2 s 2 l3s3 1
J =⎢ and B = ⎢ ⎥ (22)
⎣ 0 a2 c2 −l3c3 −s3 ⎥⎦ ⎢ −s3 l3 a2 c23 l3 ⎥
⎢ ⎥
⎣ c3 a2 s 23 ⎦
 A contribution for developing more efficient dynamic modelling algorithms 27

where l3 = a3 + q4, s23 = sin(q2 + q3) and c23 = cos(q2 + q3).

Considering that the connecting rod has a negligible inertia, the piston has mass m,
and the crank has a moment of inertia I around the actuator axis (with its centre of mass
on this axis), applying the equations (16, 17) leads to the following expression for the
generalised inertia forces: Fi* = −mu1 B1,i − Iu2 B2,i for i = 1, 2 (Bi,j indicates the element
in the ith row and jth column of the matrix B). Also considering that the piston is subjected
to a viscous friction force opposite to its motion with magnitude proportional to its speed,
and calling b this constant of proportionality, according to equation (14) the generalised
active forces (which are due to the friction in the piston, the elasticity of the connecting
rod and the input torque) have the following expression: Fi = −bu1 B1,i − kq4 B4,i + τ B2,i . It
is worth noting that the adequate choice of coordinates led to very simple expressions for
the active forces and for the corresponding partial velocities. Writing the dynamic
equations in the form of equation (19) leads to:

⎡m 0 0 0⎤ ⎡0⎤ ⎡bu1 + kc3 q4 ⎤

M =⎢ ⎥ g* = ⎢ ⎥ g=⎢ ⎥ (23)
⎣ 0 I 0 0⎦ ⎣0⎦ ⎣ ka2 s 23 q4 − τ ⎦

Knowing that V = I4×4, w = 0, J = J and w = 0, equations (20) can be used for making
direct dynamic simulations of this system when the torque τ is a known time function and
the initial conditions for q1, q2, u1 and u2 are given. The corresponding initial conditions
for q3, q4, u3 and u4 can be obtained by solving numerically the equation (21) and its time
derivative Ju = 0, with J defined in equation (22), for these variables.

4 Dynamical modelling of a 3-dof mechanism

4.1 Mechanism description

The 3-dof parallel asymmetric mechanism modelled according to the methodology
presented here is shown on Figure 3. According to Almeida and Hess-Coelho (2010), the
architecture of this mechanism may be referred as 2RSS + PPaP. A represents the
fixed base (considered as attached to an inertial reference frame) and B represents the
end-effector. The RSS chains are denoted as K and K′ and the PPaP as L. Each link of
each chain is numbered from the one joined to the fixed base to the one joined to the end
effector. The joints AK1 and AK1′ are actuated rotational joints and the joint AL1 is an
actuated prismatic joint. The joints K1 K 2 K3 , K1′K 2′ K 3′ , K 3 K 4 K 5 B and K3′ K 4′ K5′ B are
equivalent to passive spherical joints, L1L2L3 is a passive parallelogram and L3B is a
passive prismatic joint (its worth noting that the links K3 and K3′ are constituted by two
rigid bodies with a relative rotation motion around the longitudinal axis allowed, but this
motion cannot be influenced by the systems actuators, i.e., these are uncontrollable and
negligible degrees of freedom in this analysis, and that is why these links may not be
considered as constituted by different parts).
28 R.M.M. Orsino and T.A.H. Coelho

Figure 3 CAD and graph representations of the 2RSS + PPaP mechanism indicating the notation
adopted and the coordinates defined

In order to model this mechanism, the next stages (in the context of Section 2.1) are the
mobility analysis and the definition of coordinates. The Lie group of rigid body
displacements (Hervé, 1999) can be applied to study the mobility of the system. The
kinematic chains K and K′ (both RSS) do not create any constraint in the motion of the
end-effector B relative to the fixed base A, i.e., the presence of these chains do not
constraint any of the 6-dof of the rigid body motion of B relative to A. However, the
chain L (which is PPaP), does not let any relative rotational motion of B relative to A,
i.e., only translations are allowed by this chain. It is possible to conclude that the
motion of the end-effector relative to the fixed base has three degrees of freedom
(three-dimensional translations allowed), so that three rectangular coordinates are
sufficient to describe this motion. Defining a coordinate system O0x0y0z0, where O0 is the
midpoint between the centres of the active rotational joints, x0 is a downward vertical
axis, y0 is the axis which has the direction from O0 to the centre of the active rotational
joint of the chain K′, z0 is a horizontal axis parallel to the direction of the active prismatic
joint, forming a positively oriented basis according to the right-hand rule. The
coordinates of the midpoint between the spherical joints K3K4K5B and K 3′ K 4′ K 5′ B in the
coordinate system O0x0y0z0 will be denoted as q1, q2 and q3 respectively. Also defining
lines connecting the centres of consecutive joints as reference to measure angles, it is
possible to define a coordinate for each 1-dof joint of the system, as shown on Figure 3.
One important thing to be noted is that as the translation motion in the active prismatic
joint is equal the translation motion of the end-effector in the direction of the joint, the
coordinate q3 is sufficient to describe the movement of this joint; also, as L1L2L3 is a
parallelogram, a single coordinate is sufficient to completely describe the relative
motion of its parts. The constraint equations can be obtained by the method shown on
equation (5) (the homogeneous transformations may involve not only the coordinates, but
also some geometrical parameters of the mechanism), in such a way that a coordinate
system can be defined for each link with the local y axis aligned with its longitudinal
axis. It is worth emphasising that although the method proposed by equation (5) can be
computer aided, it normally will lead to a greater number of constraint equations than
 A contribution for developing more efficient dynamic modelling algorithms 29

needed. The choice of an independent set of constraint equations among all these must be
done by who is modelling.

4.2 Inclusion of actuators inputs, friction and gravitational effects

With the coordinates defined in the previous section and the homogeneous
transformations method, it is possible to describe the positions of relevant points, and by
temporal derivatives, its respective velocities and accelerations. Also the homogeneous
transformations make it possible to obtain the angular velocities of each link, according
to equation (7). There is no reason to make a different choice so that the generalised
speeds will be chosen such that u = q. All the generalised speeds can be written as
functions of u1, u2 and u3 which is a trivial set of independent variables among those
constituting u. Then, the partial velocities and partial angular velocities can be
determined.
The coordinates previously defined are also sufficient to describe the active forces in
the mechanism. According to Kane and Levinson (1985) inner forces and constraint
forces applied in fixed points do not have any effect in the system’s active forces, so that
the only contributions to these are the gravitational forces, the friction forces and the
inputs of the actuators. Noting that all the gravitational forces have the same direction of
the axis x0 and for each link the effect of all gravitational forces is equivalent to a single
force applied in the centre of mass of this link whose magnitude is equal to the total mass
of this link multiplied by the acceleration of gravity. The actuators inputs will be
considered as torques τ1 and τ2 applied in the rotational active joints of the chains K and
K′ respectively and a force τ3 applied in the prismatic active joint of the chain L.
Other important effects to be included in the model are the joints friction. Observing
that there is a coordinate to describe the relative motions in every joint of the mechanism,
the most natural way to describe the friction forces, is a constitutive relation involving
these coordinates and its derivatives. Obviously this kind of relation has to respect the
dissipative nature of these forces (i.e., in any moment these forces will restore the kinetic
energy of the system). In this analysis, the most trivial relations will be chosen, i.e., the
friction will be considered as a force opposed to the velocity of relative translational
motion in the prismatic joints and as a torque opposed to the relative rotational motion in
the rotational joints (this can be expressed as –biuiei, where bi is a positive constant,
ui = qi as previously stated and ei is the instantaneous unit vector representing the
relative translational motion or the instantaneous rotational axis, depending on the nature
of the joint).

4.3 Estimation of mechanism parameters

In order to make realistic numerical simulations of the dynamics of the 3-dof mechanism
presented, it is necessary to estimate its geometrical, inertial and frictional parameters.
This can be done with the aid of a prototype of this mechanism. In the laboratory of the
Group of Solid Mechanics and Structural Impact (from the Polytechnic School of the
University of São Paulo) a first prototype consisting only of the chain L was available to
make some experiments. With the results of these experiments, a CAD model of this
mechanism can be used to extrapolate the data, estimating all the parameters needed.
30 R.M.M. Orsino and T.A.H. Coelho

Being at disposal a prototype of the chain L, the relevant geometric parameters can be
easily measured and disassembling the links, the masses of each of them can be measured
using a scale. The inertia moments could be measured by some experiment, but in this
case they were estimated considering the geometry of the links: if they were bar shaped,
then the moment of inertia around the longitudinal axis is negligible and the moment of
1
inertia around any central axis orthogonal to the longitudinal axis is equal to ml 2 ,
12
where m is the total mass of the link and l its longitudinal dimension; otherwise, if all the
dimensions of the link are negligible compared to the total dimensions of the mechanism,
then all its inertia moments are negligible. The geometric parameters used in the model
and the links’ masses are indicated in Figure 4.

Figure 4 Masses and geometric parameters of the links of the mechanism used in the model

The friction parameters of the rotational and prismatic joints of the chain L were
measured by two simple experimental procedures. The experiment to estimate the friction
coefficients of the rotational joins of the parallelogram of L has been done considering
that disassembled of the end-effector, it can oscillate like a pendulum. For small
oscillations, a pendulum can be treated as a linear second order dynamic system and the
linear rotational friction coefficient can be estimated by the elapsed time to dissipate part
its initial mechanical energy (which is the elapsed time to not be able to achieve a given
amplitude angle any more). The value of the friction coefficient estimated for these joints
was 0.001 N ⋅ m ⋅ s/rad and by hypothesis all the other rotational joints of the mechanism
were considered to have the same friction coefficients as the rotational joints of the chain
L. In order to estimate the friction coefficient of the prismatic joint AL1, it has been done
an experiment where the chain L was put in motion in the direction of this joint and the
force and displacement time histories were acquired by a device consisting of a load cell
in series with a LVDT. Analysing in these time histories and verifying the mean value of
the friction force in time intervals where the speed is approximately constant, it can be
verified a correlation between friction force and the speed of the displacement, so that
adjusting the experimental results to the friction model using least square methods, the
 A contribution for developing more efficient dynamic modelling algorithms 31

value of the friction coefficient is approximately 15 N ⋅ s/m (by hypothesis the friction
coefficient of the joint L3B is supposed to be equal).

4.4 Numerical simulations of the dynamic model

In order to use the model to analyse the dynamic behaviour of the system modelled and
the influences of all the dynamic effects included, it is necessary to make some numerical
simulations that may be able to represent some typical use of this mechanism.
The simulations presented here have been done according to the following procedure.
A desired motion is imposed to the mechanism as shown on the plots of Figure 5; this
motion represents a rectilinear motion of the end-effector with null speeds and
accelerations at both extremities of the trajectory. The time plots of the accelerations
corresponding to this motion are sinusoidal. Then, an inverse kinematic and dynamic
analysis is done with the dynamic model, determining the time histories of all the other
coordinates of the system and the actuators inputs needed to enable this motion. The time
histories of these inputs are shown on the plots of Figures 6 and 7.
It is worth noting that the torque inputs necessary to move the whole system in the
desired trajectory are not much different from the static loads required to equilibrate the
system at the start and the end of the trajectory (values of the inputs in 0.0 s and 0.5 s),
showing one of the main advantages of this architecture of mechanisms. Nevertheless,
this kind of behaviour shows how sensitive is the motion of end-effector with respect to
the inputs and emphasises the importance of the dynamic model in the synthesis of a
control system for such kind of mechanisms. If the inputs provided by the actuators are
not close enough to the inputs corresponding to a desired trajectory, the real path
described by the end-effector will be quite different from the desired, so that the dynamic
model has an important role in the design of a good control system.

Figure 5 Trajectory imposed to the 2RSS + PPaP mechanism to make an inverse kinematic and
dynamic simulation
32 R.M.M. Orsino and T.A.H. Coelho

Figure 6 Time plots of the torques inputs τ1 and τ2 (active rotational joints of the chains
K and K′ respectively) corresponding to the imposed trajectory

Figure 7 Time plot of the force input τ3 (active prismatic joints of the chain L) corresponding to
the imposed trajectory

5 Conclusions

The methodologies discussed, although being well known in the literature (and some of
them being largely used to model the traditional serial mechanisms) can be specialised to
applications in the modelling of parallel kinematic structures, taking advantages on the
definition of redundant coordinates in order to make the mathematical model as simple as
possible. A good kinematic analysis has been shown to be a fundamental task of the
 A contribution for developing more efficient dynamic modelling algorithms 33

modelling process, because the application of Kane’s method automatically determinates

the generalised inertia forces based on this analysis and a well done selection of a set of
coordinates can significantly simplify the expressions of the generalised active forces.
Moreover, Kane’s method has the same advantages of the virtual works and Lagrange’s
approaches, automatically considering the system’s constraints (without the need of
constraint forces as in Newton-Euler traditional approach) but having the advantage of
dealing with them in an easier way, defining the generalised speeds, relating and
eliminating them by means of the constraint equations (when the constraints are
holonomic or simple non-holonomic), not requiring the definition of undetermined
multipliers or some other non-methodical elimination process.
The 2RSS + PPaP mechanism modelled, although having a quite complex dynamics
(being an asymmetric and three dimensional parallel structure) could be treated by the
methodology presented, being necessary only the definition of one coordinate for each
relative motion on each joint. This led to simple expressions for the positions, velocities
and accelerations of the main points as well as for the angular velocities and angular
accelerations of the links and the gravitational and friction forces. The mathematical
model is able to determine uniquely the configuration of each link of the system (some
coordinate sets, although sufficient to describe the dynamics of the system as a whole, are
not able do determine exactly the configurations of all the links, being possible more than
one configuration of the system corresponding to the same coordinates set; however, the
choice of redundant coordinates can eliminate this problem as done in this analysis) and
all the dynamic effects could be suitable included, showing all the advantages of the
methodology developed in the analysis of a complex parallel kinematic structure.

Acknowledgements

The authors would like to acknowledge the National Council for Scientific and
Technological Development (CNPq, Brazil) for having sponsored the undergraduate
research of Renato Maia Matarazzo Orsino (from July 2010 to July 2011), which has
been the beginning of the research project presented in this paper.

References
Almeida, R.Z.H. and Hess-Coelho, T.A. (2010) ‘Dynamic model of a 3-dof asymmetric parallel
mechanism’, The Open Mechanical Engineering Journal, Vol. 4, pp.48–55, (online) available
at http://www.benthamscience.com/open/tomej/articles/V004/SI0036TOMEJ/48TOMEJ.pdf
(accessed on 9 January 2012).
Craig, J.J. (1989) Introduction to Robotics: Mechanics and Control, Addison-Wesley Longman
Publishing, USA.
Hervé, J.M. (1999) ‘The Lie group of rigid body displacements, a fundamental tool for mechanism
design’, Mechanism and Machine Theory, Vol. 34, No. 5, pp.719–730.
Kane, T.R. and Levinson, D.A. (1985) Dynamics: Theory and Applications, McGraw Hill
Publishing Company, USA.
Lanczos, C. (1983) The Variational Principles of Mechanics, Dover Publications, Inc., New York,
USA.
Leech, J.W. (1971) Mecânica Analítica, Ao Livro Técnico S.A. e Editora da Universidade de
São Paulo, Rio de Janeiro, Brazil.
34 R.M.M. Orsino and T.A.H. Coelho

Merlet, J-P. (2006) Parallel Robots, 2nd ed., Springer, Netherlands.

Mitiguy, P.C. and Kane, T.R. (1996) ‘Motion variables leading to efficient equations of motion’,
The International Journal of Robotics Research, Vol. 15, No. 5, pp.522–532.
Tsai, L-W. (1999) Robot Analysis: The Mechanics of Serial and Parallel Manipulators, John Wiley
and Sons, USA.
Udwadia, F.E. and Phohomsiri, P. (2007) ‘Explicit Poincaré equations of motion for general
constrained systems. Part I. Analytical results’, Proceedings of the Royal Society, Vol. 463,
pp.1421–1434.