Geometric methods and formulations in computational multibody system dynamics

Multibody systems are dynamical systems characterized by intrinsic symmetries and invariants. Geometric mechanics deals with the mathematical modeling of such systems and has proven to be a valuable tool providing insights into the dynamics of mechanical systems, from a theoretical as well as from a computational point of view. Modeling multibody systems, comprising rigid and flexible members, as dynamical systems on manifolds, and Lie groups in particular, leads to frame-invariant and computationally advantageous formulations. In the last decade, such formulations and corresponding algorithms are becoming increasingly used in various areas of computational dynamics providing the conceptual and computational framework for multibody, coupled, and multiphysics systems, and their nonlinear control. The geometric setting, furthermore, gives rise to geometric numerical integration schemes that are designed to preserve the intrinsic structure and invariants of dynamical systems. These naturally avoid the long-standing problem of parameterization singularities and also deliver the necessary accuracy as well as a long-term stability of numerical solutions. The current intensive research in these areas documents the relevance and potential for geometric methods in general and in particular for multibody system dynamics. This paper provides an exhaustive summary of the development in the last decade, and a panoramic overview of the current state of knowledge in the field.


Introduction
Space kinematics is solely based on screw theory, and consequently so is the kinematics of multibody systems (MBS).Even though such concepts are not so widely known in the MBS community, which is an obstacle hindering a fruitful exploitation of such concepts for computational MBS dynamics.On the other hand, a recent trend in MBS dynamics is to employ the terminology and certain concepts of Lie groups noting that rigid body motions, i.e., finite frame transformations, form such a group possessing certain desirable properties.The link between the theory of screws and motion groups is the fact that the screw algebra is nothing but the Lie algebra of the Lie group in question.The aim of this review paper is to discuss the significance of the Lie group concept, to summarize the recent work and to discuss open research problem in this area.It is shown that such an approach gives rise to compact formulations of the basic kinematic relations and to advanced geometric integration schemes.
Spatial kinematics rests on the theory of screws as originally discussed by Ball [10].From today's perspective, this seems to be most natural, but it would not have happened without the major overhaul of geometry in all its variants triggered by Klein's Erlanger Program [83].The basis of screw theory is the line geometry introduced by Plücker [133,134] together with the projective geometry of Grassmann [68] and Cayley [36].The projective theory of screws was advanced in the context of kinematics by Study [55], von Mises [156], Hunt [77], Blaschke [15], Gibson and Hunt [65,66], Dimentberg [53], and many others, and now provides a consistent and self-contained geometric framework.Klein's work was accompanied by Lie's development of continuous transformation groups.Even more, the discovery that transformations of (Euclidean) spaces form what is called today Lie groups is an essential aspect in kinematics.Furthermore, the inextricable connection of Lie groups and spatial motions is a cornerstone of modern kinematics, and consequently of robotics, rigid body dynamics, and so forth.It is even more remarkable that screw and Lie group theory is underrepresented in computational MBS dynamics.It is even more remarkable that Liu [91] already presented an approach to MBS modeling based entirely on screws.The significance of screw theory for mechanism kinematics and dynamics cannot be too highly praised.As a curiosity a poem entitled "Song of the Screw" was published in Nature [149] during the time when Ball worked on the theory of screws.
A MBS is a set of solid (rigid or flexible) bodies kinematically or physically interconnected.If a MBS consists of rigid bodies, it is called a discrete mechanical system [1,5].Approaches to MBS modeling have traditionally been split into two main categories, namely the so-called absolute and relative coordinate approaches.The motivation to apply Lie group methods to the absolute and relative coordinate formulations is different.The former aims to apply geometric Lie group integration methods, whereas the latter aims to employ the geometry of screw motions for the modeling purposes.This will be discussed in the paper.Lie group methods are the bridge between differential geometric approaches to MBS dynamics and geometric modeling concept applicable to computational MBS dynamics.
The area of Lie group and geometric methods is a relatively young discipline.Various approaches together with different notations have merged.A review of this upcoming area is hence timely.In Sect.2, the geometry of the rigid body motion is discussed and the configuration space is identified with a Lie group.The differences of the semidirect product group SE (3) and the direct product SO (3)×R 3 are discussed.Velocities are introduced by left or right trivialization.The rigid body dynamics is reviewed in Sect.3. The governing equations are written in the form of Euler-Poincaré equations.Therewith the equations governing the kinematics and dynamics are consistently written in Lie group setting.Section 4 recalls the modeling of lower pair joints using canonical coordinates on the motion subgroups defined by lower pair joints.This is then the basis for modeling MBS with tree structure in terms of relative coordinates in Sect. 5.In this section, recursive relations for the configuration, velocity, and acceleration are recalled.Emphasis is given on the various forms that were reported in the literature.Finally, the motion equations are expressed algebraically in closed form.The Lie group modeling in terms of absolute coordinates is presented in Sect.6.The numerical integration of these model equations is discussed in Sect.7. In each section, the current trend and open issues are discussed.

The geometry of rigid body motions
A proper representation of rigid body motions is essential for any Lie group MBS modeling.Rigid body motions form a Lie group-the group if isometric and orientation-preserving transformations of Euclidean space.While the geometry of rigid body motions is well understood, several representations of rigid body motions have been introduced.

Rigid body motions as frame transformations
A rigid body is regarded as an Euclidean space equipped with a mass density function.Kinematically, a rigid body is represented by a body-fixed reference frame with respect to which the mass density is defined.The configuration of a rigid body is then represented by a rotation matrix R ∈ SO (3) and a position vector r ∈ R 3 .The configuration is denoted with C = (R, r) ∈ SE (3).The Lie group SE (3) = SO (3) R 3 of rigid body motions is the semidirect product of the rotation group SO (3) and the translation group with the group multiplication These matrices are called homogenous transformation matrices [119] since they are used to transform homogenous point coordinates.The group multiplication, i.e., the concatenation of frame transformations, is then the matrix multiplication C 1 C 2 .The inverse of the transformation (2) is , respectively,

Instantaneous screws and canonical coordinates
In matrix representation, the Lie algebra se (3) is the vector space of matrices of the form where ξ ∈ so (3) is the skew symmetric matrix associated to the vector ξ ∈ R 3 .Notice that occasionally the 'hat' notation is used instead of the tilde, as for instance in [119].Here the hat denotes isomorphism The Lie bracket on se (3) is given by the matrix commutator With the isomorphism (5), this Lie bracket reads with The Lie bracket ( 6) is also called the screw product [24,144].The matrix (7) has appeared under different names, such as 'spatial cross product' in [61,62,82], or the 'north-east cross product' [22].The closed form of the exp mapping on SE (3) is found evaluating the matrix exponential with the Euler-Rodrigues formula and dexp is the right-trivialized differential of the exp mapping on SO (3) The expressions (10) and ( 13) are preferable from a computational point of view [114].
The important point for application to MBS modeling is that se (3) is isomorphic to the algebra of screws [119,144].It actually took some time until their relation was fully understood since Ball [8][9][10] introduced the concept of screws.An excellent review on this topic can be found in [47,48].With a chosen reference frame, a general screw is represented by a screw coordinate vector X = (ξ , η), as in (5), where ξ ∈ R 3 is the angular and η ∈ R 3 is the translational component.While Ball introduced the notion of screws basically as a tool for static analysis, it was first Mozzi [109] who proved that any rigid body motion possesses an instantaneous screw axis about which the body is performing an instantaneous screw motion.Chasles [38] later proved that any finite rigid body displacement can be achieved by a screw motion about a constant screw axis.
Chasles theorem gives rise to the explicit construction of the screw coordinates for a finite screw motion.Denote with e ∈ R 3 the unit vector along the axis, and with p ∈ R 3 the position of any point on that axis.Then, the unit screw coordinate vector is [102,119,144] Here, h ∈ R is the pitch of the screw, i.e., h = 0 for pure rotations, and h = ∞ for pure translations.For the case h = ∞, the screw coordinates ( 14) are the Plücker coordinates of a line at infinite.With this particular form, the exp mapping (8) becomes The exponential mapping in its closed form is the foundation for the Lie group modeling of MBS kinematics.The six elements of the screw coordinate vector X ∈ R 6 in (8) represent canonical coordinates of the first kind on SE (3) (also called exponential coordinates).Therewith the configuration of a rigid body is determined as C (X) = exp( X).The angular part ξ ∈ R 3 is indeed the instantaneous rotation axis.A rigid body motion can also be described as combination of successive one-dimensional motions.This corresponds to expressing SE (3) as product of one-parameter subgroups.Then, the configuration is parameterized as C (ϕ) = exp(ϕ 1 e 1 ) exp(ϕ 2 e 2 ) • • • exp(ϕ 6 e 6 ) with constant basis vectors e i ∈ R 6 .The ϕ i , i = 1, . . ., 6 are canonical coordinates of the second kind.

Frame transformations of screws
Frame transformations are described by the conjugation map.Let C ∈ SE (3) be the configuration of a body w.r.t. a chosen frame.Now let B ∈ SE (3) be the transformation to another reference frame.Then, where L B and R B are left and right translation map on SE (3), respectively, and Here Ad B : se (3) → se (3) is the adjoint transformation (111).Using screw coordinates, the adjoint mapping is X = Ad B X with for B = (R, r).

Rigid body velocities: twists
Using operations on the Lie group, there are two ways to introduce the twist of a rigid body: the spatial twist and thus the twist screw vectors are Here ω s i = Ṙi R T i and ω b i = R T i Ṙi are the body-fixed and spatial angular velocity, respectively.The vector v b i = R T i ṙi is the velocity of the body-fixed frame measured in the inertial frame and resolved in the body-fixed frame.The spatial translational velocity, v s i = ṙi + r i × ω s i , is the velocity of the point of the body that is momentarily passing through the origin of the inertial frame and resolved in this frame.

Kinematic reconstruction
The twist of a rigid body does not reveal its actual finite motion.The determination of the actual motion of the rigid body amounts to solving the Eq. ( 19) that can be rewritten as These are the right and , respectively ,left Poisson equations on the Lie group SE (3).They are referred to as the kinematic reconstruction equations [76,98].These are nonlinear ODEs in the Lie algebra SE (3).That is, their solution must preserve the nonlinear invariants defining the group.More precisely, these are the invariants of SO (3), namely the orthogonality conditions.

Direct product representation of rigid body motions
The configuration of a rigid body relative to a chosen frame can also be regarded as an element of the direct product of the rotation and translation group: C = (R, r) ∈ SO (3) × R 3 .However, then the multiplication in the direct product group is Apparently, this multiplication law does not describe frame transformations.The direct product group does hence not represent rigid body motions.Nevertheless, it is frequently used as geometric model for MBS kinematics and for Lie group integration methods (see Sect. 7).The difference becomes significant if the direct product multiplication is employed in the position update of time integration schemes as will be shown later.
The inverse element is (R, r) −1 = (R T , −r).If desired, SO (3) × R 3 can be represented by the group of 7 × 7 matrices of the form [116] The multiplication is then represented by C 1 C 2 .

Canonical coordinates
The Lie algebra of the direct product SO (3)×R 3 is so (3)×R 3 .The latter is isomorphic to R 6 , and its elements can be represented as vectors X = (ξ , r) ∈ R 6 .The exponential mapping on the direct product group is with the exponential mapping (9) on SO (3).Apparently, X is not a screw coordinate vector, but rather consists of a scaled rotation axis and a translation vector.The parameters are inherited from the rotation group SO (3) and translation group R 3 .Moreover, ξ ∈ R 3 is the scaled rotation vector in (9), and r ∈ R 3 is the translation vector.These serve as canonical coordinates on the direct product group.

Rigid body velocities
Also for the direct product representation rigid body velocities can be defined using left or right trivialization. For , the left-and right-trivialized derivatives yield, in vector notation, These are not proper twists, i.e., screws.V m is commonly referred to as the mixed velocity (it contains a mix of body-fixed angular velocity ω b and translational velocity ṙ resolved in the inertial frame), whereas V h is referred to as the hybrid velocity [119].
The dexp mapping corresponding to ( 27) is with the dexp mapping on SO (3) in (11).Thus, the mixed and hybrid velocities are determined as for X = (ξ , r).The dexp mapping reveals once more the problematic point of the direct product representation namely that the angular and linear velocities are decoupled.
The Lie bracket on this algebra is This can be expressed as ad X 2 X 1 with matrix

Kinematic reconstruction
The reconstruction equations split according to the direct product representations into the integration of ṙ and These reconstruction equations are seemingly easier to solve than those for the semidirect product representation in (25).It is, however, obvious that independent integration of angular and translational part does generally not produce a screw motion, i.e., a rigid body motion.

Rigid body dynamics on a lie group
In a geometric setting, the motion of a rigid body is a curve in a Lie group.Accordingly, the rigid body dynamics is governed by the differential equations on this Lie group.These are the left-and, respectively, right-trivialized Euler-Poincaré equations.Therewith the motion equations are defined in a coordinate-free way.A particular choice of representation then yields a particular form of the motion equations.

Euler-Poincaré equations on SE (3)
The dynamics of a rigid body is governed by the well-known Newton-Euler equations.In a geometric setting, the Newton-Euler equations of a free rigid body in terms of body-fixed twists are the left-trivialized Euler-Poincar é equations on SE (3) Here b ∈ se * (3) (se * (3) is the dual space to se (3) as a vector space) is the body-fixed momentum screw defined as b := M b V b where is the (constant) body-fixed inertia matrix given in terms of the body-fixed inertia tensor b and the displacement vector b d bc from the origin of the body-fixed reference frame to the center of gravity resolved in the body-fixed frame.
In spatial representation, the Newton-Euler equations attain the form Here, M s is the spatial inertia matrix, and s ∈ se * (3) the spatial momentum screw.The Eqs. ( 38) and ( 41) are formally identical.The Eqs. ( 37) and ( 40) describe the conservation of momentum.They govern the dynamics on se * (3) after left or right trivialization, respectively.Application of Hamilton equations and use of canonical momenta are not common MBS dynamics.Rather the Eqs.( 38) and ( 41) are used that govern the dynamics on se (3).
Remark 1 The Eq. ( 37) can be derived invoking the concept of Lagrange reduction as it has been advanced by Marsden and Scheurle [100].Lagrange and Poisson reductions for general mechanical systems are powerful tools in geometric mechanics.In the general situation on a Lie group G with algebra g with given left-and right-trivialized Lagrangians L l (g, ξ) and L r (g, ξ), respectively, the reduced equations are with ad * ξ : g * → g * being the dual to the adjoint operator on g.The last equation in ( 42) and ( 43), respectively, is indeed the respective kinematic reconstruction equation.For Lagrangians that are G independent, i.e., L does not depend on g ∈ G (as for the free rigid body in body-fixed description), these equations reduce to the Euler-Poincaré equations.More on reductions can be found in [17,98,100], and also on the extension to continuum mechanics [99].

Rigid body equations of motion on SO (3) × R 3
If the direct product SO (3) × R 3 is used as configuration space, the motion equations cannot be derived as the Euler-Poincaré equations.This is a hint that the configuration space does not capture the geometry of rigid body motions.In other words, a different geometric model is used for the kinematics and for the dynamics.
In terms of mixed velocities, the Newton-Euler equations are When a reference frame at the center of gravity is used they simplify to It is important to notice that only in the latter case the angular and translational momentum balance are decoupled, and that only then the geometric model of a direct product group is applicable.The consequences will be discussed in Sect.7. The (46) can also be written as left-trivialized Euler-Poincare equations where M m is the mass matrix (39) with b d bc = 0, and the matrix (34).

Note on the kinematic reconstruction
Solving the Euler-Poincaré equations yields the time evolution of the spatial and body-fixed twist, respectively.In summary, the Euler-Poincar é equations ( 37) and ( 41) together with the Poisson equations ( 25) describe the dynamics and kinematics of a free rigid body, either in spatial or in body-fixed representation.The dynamic equations ( 37) and ( 41) are ODEs on se (3).
The solution of the kinematic reconstruction equations may seem to be a rather obvious task but these equations have always been a source for problems, mainly for two reasons: (1) there is no singularity-free parameterization of rotations, and (2) the solution must respect nonlinear invariants intrinsic to SE (3).Classically the remedy consists in using Euler parameters (unit quaternions) [143].This itself leads to numerical issues since the Euler parameters must satisfy a unit length constraints.The latter is a quadratic form and deemed to be easy to satisfy within numerical integration schemes.Yet, the necessary constraint stabilization, i.e., renormalization, is a potential cause for errors.This is classically attributed to Poisson, and ( 25) and (35,36) are called the kinematic Poisson equation [142].More precisely, (35) is referred to as the left Poisson equation and (36) as the right Poisson equation.However, the problem of determining the rotation when the angular velocity is given has also been addressed by Darboux [49].For this reason, Condurache [46] called them the kinematic Poisson-Darboux equation.Attempts to solve them explicitly can still be found in the literature.As such, exact solution for special situations was proposed [44].As special case, the Kepler problem of relative orbital motion was solved in [46].Also the problem of solving the kinematic Poisson problem in a stochastic setting was reported [69].

Relative motions as screw motions
Most technical systems are built with the so-called lower kinematic pairs [2,3,78] that are characterized by surface contact.The contact surface is thus an invariant subspace of the relative motion two bodies connected by a lower pair joint.This was studies already by Reuleaux [137,138], which is why lower pairs are also called Reuleaux pairs.Lower pair joints correspond to subgroups of SE (3).These are the isotropy groups leaving the contact surface invariant.The correspondence of motion subgroups and lower pairs has been discussed in [144].As the space of relative motions (disregarding restrictions due to contact) is a subgroup, the vector space of relative twists is a subspace of the corresponding subalgebra, i.e., it is involutive.Motivated by this observation, Stramigioli [150] uses this as the defining property of lower pairs.Recently, Selig [145] and Wu et al. [158] seized on a concept by Brockett [28] called Lie triple systems.The latter refers to systems of vector fields that are closed under the triple Lie bracket (whereas a Lie algebra is closed under the Lie bracket).Wu [158] pointed out that this concept allows for studying the kinematic of non-lower pair joints, such as constant velocity joints.
The fact that the motion of any lower pair joint can be parameterized by canonical coordinates of second kind, i.e., by 1-DOF screw motions, proves that any lower pair joint can be modeled as a combination of revolute and prismatic joints, which was the approach when computational MBS dynamics started [143,157].This notion also goes back to Denavit and Hartenberg [50] (Table 1).
Not all subgroups correspond to lower pairs.But they can be generated by combination of lower pairs

Partitioning and canonical coordinates
The standard approach in MBS modeling is to introduce body-fixed joint frames whose relative motion is describe by a screw motion.A detailed exposition of this approach can be found in [154,155].
Variant 1 Consider two bodies, indexed with i − 1 and i.Let S i−1,i be the fixed transformation from the joint frame to the reference frame on body i − 1, and with S i,i the constant transformation from the joint frame to the reference frame on body i.In the classical MBS approach, they are parameterized in terms of Denavit-Hartenberg parameters.The joint motion is then determined by the relative motion of the two joint frames.For a lower pair joint, this is expressed as exp(Z i q i ) where is the associated screw coordinate vector, and q i the joint variable (angle, displacement).The notation i x indicates that the vector is resolved in a frame on body i.The relative configuration Variant 2 In the above formulation, the joint screw coordinate vector is represented in the joint frame on body i − 1.Since the Lie group formulation is frame invariant, any other frame can be used.In particular, the reference frame on body i can be used.This leads to where the constant part is the reference configuration of body i w.r.t.body i − 1 for q i = 0.This is referred to as the zero reference configuration of the joint.The screw coordinate vector of joint i represented in the reference frame of body i is Therein, i x i,i is the position vector of a point on the axis of joint i, and i e i is the unit vector along the joint axis, both resolved in the frame in body i.
Variant 3 Instead of using the reference frame on body i, the reference frame on body i − 1 can be used.This yields the alterative form of ( 50) (53) expressed in the frame on body i − 1.Now i−1 xi−1,i is the position vector of a point on the axis of joint i measured in the frame on body i − 1.
The last two forms ( 50) and ( 52) bear important features that make them superior for efficient MBS modeling.This clearly demonstrates the advantage of using the Lie group, respectively, screw, approach.The remarkable point here is that no joint frames need to be introduced, and that no convention (like Denavit-Hartenberg parameters) needs to be respected.

Open issues
There are still several open issues regarding the modeling of general joints.The historically prevailing problem is the singularity-free description of spatial rotations, i.e., of spherical joints.This is not problematic within the Lie group framework since this is a coordinate-free approach.The problem arises as soon as an explicit canonical parameterization is introduced.Classically this is remedied by use of Euler parameters.However, in the Lie group setting, parameters are only required during the time integration.To this end, local parameters are introduced to describe the configuration increment.Therewith a singularity-free description is achieved (see Sect. 7).
Beyond lower pairs, complex mechanical systems comprise higher kinematic pairs (e.g., curve contacts) as well as the so-called non-holonomic joints.A non-holonomic joint is a higher kinematic pair imposing non-holonomic constraints to relative motion of the respective two links.That is, a non-holonomic joint is characterized by a non-integrable constraint distribution defining feasible relative velocities.Such mechanisms were proposed and analyzed by Grosch et al. [70,71], di Gregorio [51,52], and Duindam and Stramigioli [54].In [54], a modeling approach was proposed that makes use of local exponential coordinates on the covering Lie subgroup defining the relative joint motions.A more holistic approach was taken by Chhabra [42,43] who introduced a generalized exponential formula.

Relative coordinate modeling of multibody systems in Lie group setting
The term 'relative coordinate' modeling refers to the use of local parameters to describe joint motions and thus the relative motion of adjacent bodies as in Sect. 4. Historically, the relative coordinate modeling of mechanisms comprising ideal joints has been the subject of analytical mechanics, ref. to Papastavridis [124] and Lurie [92] 5.1 Configuration of a kinematic chain using the product of exponentials For illustration purpose, consider a kinematic chain with n rigid bodies interconnected by n lower pair 1-DOF joints.The central relation for the Lie group description of a serial kinematic chains is the Product-of-Exponentials (POE) formula.Although Brockett presented it in a conference paper [27], the POE formula is a key relation in the modeling of rigid body mechanisms.Later, Park [125] discussed its computational aspects and presented two variants as follows.
Denote with q ∈ V n the vector of generalized coordinates (canonical joint coordinates).The relative motions due to lower pair joints, as described in the preceding section, can be concatenated as It will be called the body-fixed Product-of-Exponentials (POE) formula in body-fixed description since the joint kinematic is expressed by exponentials of joint screws.The first version of (54) was reported in [132], and both forms in [125,126].Common to both formulations (54) is the exclusive use of body-fixed reference frames, while no joint frame is necessary.The difference though is at which body the reference frame is used.This resembles the different forms of the Denavit-Hartenberg conventions as discussed in [3].
The body-fixed POE formulation possess two advantages compared to the classical MBS modeling.Firstly it does not require introduction of joint frames.Secondly, the screw coordinates X i of joint i represented in the reference frame on body i, respectively, Xi represented in the frame on body i − 1, are given in terms of readily available geometric data as in (51), respectively (53).
Yet the body-fixed POE formula can be reformulated in terms of joint screw coordinates represented in the spatial inertial frame as follows where is the absolute 'zero reference' configuration of body i.Now Y j = Ad A j j X j = e j y j × e j + h j e j (57) is the screw coordinate vector of joint j represented in the inertial frame in the reference configuration with q = 0, where e j is the direction unit vector, and y j the position vector of any point on the joint axis, both resolved in the inertial frame.This is also the relation of the spatial representation and body-fixed representation of joint screw in the reference configuration.The spatial POE formula (55) is the original form introduced by Brockett [27] for the kinematic description of robotic manipulators.The initial transformation A i is referred to as the 'zero reference.'Without resorting the exponential mapping, a 'zero reference formulation' was already reported by Gupta [72] in terms of frame transformation matrices.The POE formulae (55) was first used for MBS modeling in [30].
The main advantages of the spatial POE formulation is that all required data are represented in the spatial inertial frame.It thus allows for formulating the MBS kinematics without body-fixed joint frames and without referring to body-fixed reference frames.This aids the MBS modeling since no joint kinematics, i.e., S i,i , S i−1,i or B i , need to be described.Rather the absolute reference configurations A i w.r.t. to the inertial frame and the screw coordinates (57), i.e., e i and p i are required.

Velocity of tree topology MBS
The configuration of body i is given recursively as C i (q) = C i−1 (q) B i exp(X i q i ).Inserting this in the definition of body-fixed and spatial velocity (20) yields the recursive relation for the body-fixed and spatial twist Resolving the recursion, the POE formula leads to where are the instantaneous screw coordinates that constitute the columns of the body-fixed and spatial Jacobian J b i and J s i .The body-fixed version of the recursive relations ( 58) and ( 62) was reported in [125,126,132].The spatial version was presented in [115].
The body-fixed Jacobians have been used already by Maißer in [94].The remarkable fact that the Jacobians are given algebraically motivated to call them 'kinematic basic functions' since they are the intrinsic objects required for MBS kinematics.Angeles [4] termed the Jacobian the 'natural orthogonal complement.' Denote with T the overall body-fixed twist vector of the kinematic chain.This is determined as where the 6n × n matrix J b = A b X b is the body-fixed system Jacobian, and Analogously denote with V s = (V s 1 , . . ., V s n ) T the overall spatial twist vector.This is given as where the spatial system Jacobian possesses the factorizations with

Acceleration of tree topology MBS
The central relation that allows for the recursive and derivative-free generation of motion equations is the following closed-form expression for the partial derivatives of the body-fixed Jacobian This only involves the computationally simple Lie bracket (6).It yields the closed-form expression for the acceleration in a very concise and compact form The expression (71) has been derived, using different notations, for instance in [27,115,119,125].
Therewith the time derivative of the Jacobian is This gives rise to the time derivative of the body-fixed system Jacobian where This leads to a compact matrix form of the system acceleration An analogous derivation yields ∂J s and thus for the time derivative of the spatial Jacobian and the spatial acceleration The overall spatial acceleration is thus expressed as The expressions for the body-fixed acceleration were presented in [125,126,132].The role of the Jacobian was, however, not discussed.It is clear that the Jacobian is the essential object in the MBS kinematics description.Its derivatives are given simply in terms of the Lie brackets, i.e., simple algebraic operations.It is to be mentioned that the following closed-form expression for the higher-order partial derivatives exists [112,113] where is the ordered sequence of the indices α 1 , . . ., α ν .

Motion equations in closed form
Finally the above algebraic relations allow for constructing the motion equations of a MBS in closed form.They are commonly written in the form M (q) q + C ( q, q) q = Q ( q, q, t) (81) where M denotes the generalized mass matrix, C ( q, q) q represents the Coriolis and centrifugal forces, and Q stands for all other generalized forces, including potential, dissipative, and applied forces.
Using body-fixed twists, the kinetic energy of the MBS comprising n rigid bodies is qT M q where the generalized mass matrix is with 74) the time derivative of the conjugate momentum vector The matrix representing Coriolis and centrifugal terms can thus be identified as The expressions ( 82) and ( 84) are given algebraically, and the Lagrangian motion equations are thus algebraically determined in closed form.The only information needed are the instantaneous joint screw coordinates in body-fixed representation as it was discussed in [125,126,132].
It is common to rephrase the Lagrange equations as the equations of geodesic in the configuration space V n equipped with the mass matrix as Riemannian metric n j=1 M i j (q) q j + n j,k=1 Here = Γ ik j are the Christoffel symbols.These possess a closed-form expressions, obtained using the relations ( 62) and ( 70), The expression (86) for the Christoffel symbols was reported in [30,115] using Lie group notation.Prior to these publications, an equivalent expression was presented in [94] using tensor notation, where the bodyfixed Jacobians were called 'kinematic basic functions.'Also Burdick [35] presented a similar, although not equivalent, algebraic form of the Lagrange equations.
Analyzing the sensitivity of the motion of a MBS requires linearization of the motion equations, which amounts to taking partial derivatives.It is clear from the preceding derivation that this shall be possible in a purely algebraic way.Accordingly, Sohl and Bobrow [148] presented a closed-form algebraic formulation in Lie group setting.Also this is equivalent to the one presented in [95] using tensor notation.In [110], a closed form for the partial derivatives of the inverse mass matrix was presented, and in [64] an explicit form of the time derivatives of the motion equations.

Open issues
The computationally efficient evaluation and solution of the motion equations resorts to recursive O (n) formulations.In the Lie group context, such have been proposed by Ploen et al. [130,132].The point of departure is the inverse of the matrix (65) that can be expressed as with the nilpotent matrix This is the common starting point for various O (n) methods [60][61][62]80,81] .Thus, the basic operations are the same.However, the actual number of operations depends on the particular transformations involved.That is, the representation of twists is decisive about the number of operations.When using spatial twists, for instance, there are no frame transformations in the recursion (58).This is the motivation behind the algorithm introduced by Featherstone [62].On the other hand, the instantaneous screw coordinates as well as the spatial mass matrix must be transformed.Several established O (n) algorithms use the hybrid twist, notably the spatial operator algebra introduced by Rodriguez [139][140][141].Orin et al. [121,122] analyzed the computational effort when using the different representations.This was made under the premise of modeling convention.The actual effort remains to be investigated when the general screw description is used.The body-fixed representation of twists is still dominant in MBS dynamics, while the use of spatial twists is dominant in mechanism theory [103].
In certain applications, it is necessary to determine the motion equations in closed form.The above closed form is derived using body-fixed twists.As for the recursive O (n) formulation, the potential benefit of using spatial, hybrid, or mixed twists remains to be explored.
Most Lie group formulations so far have been focused on lower pair joints that allow for application of the exponential mapping to express relative motions.Recently, it was attempted in [42,43] and [54] to generalize this approach to non-holonomic systems.Along this line, an approach to model kinematic couplings that do not generate a motion subgroup was presented in [158].
Most MBS possess kinematic loops.Within the relative coordinate modeling, this means that the loop constraints must be formulated in terms of POE formula.Moreover, their derivatives must be determined algebraically.There is yet no systematic approach to this important issue.
Rigid body motions can be represented by dual quaternions.This was advanced by Chevallier [39][40][41] for application to MBS.Also recently Shoham proposed dynamics formulations using dual quaternions [31,45].The motivation behind this is that this is deemed to be more compact than the screw formulation.
Non-canonical parameterizations such as Cayley-Rodrigues parameters allow for algebraic description of motions on the expense that they are not globally valid.Like Euler parameters, such algebraic parameterizations, also called vector parameters, are deemed advantageous.Lie group formulation in terms of vector parameters has been presented in [106][107][108].Their benefits are to be explored.

Ambient space of an MBS
A rigid body is kinematically represented by an element of a 6-dimensional Lie group G (either SE (3) or SO (3) × R 3 ).The starting point for a Lie group formulation is to describe the motion of a MBS as a curve in a product Lie group.For a MBS comprising n rigid bodies the, ambient configuration space is the 6n-dimensional direct product Lie group where G i is the Lie group used to describe the motion of body i w.r.t. to the inertial frame.The absolute configuration of the MBS is then g and the identity element is I = I G 1 , . . ., I G n with the identity element The Lie algebra of the ambient space G is g := g 1 × • • • × g n .A typical element has the form X = (X 1 , . . ., X n ) ∈ g.A Lie bracket on this direct product Lie algebra g is inherited componentwise from there exists a X i ∈ g i , serving as canonical coordinates, such that C i = exp X i , with the exp mapping on G i .The induced exp mapping on the direct product group is exp : Correspondingly, its right-trivialized differential is dexp In this setting, X ∈ g represents absolute canonical coordinates of the MBS.
The tangent space of G is isomorphic to the Lie algebra g via left-trivialization Accordingly, the velocity of the MBS can be introduced via left-trivialization as V := T e L −1 g • ġ = g −1 ġ ∈ g or right-trivialization as V := T e R −1 g • ġ = ġg −1 ∈ g.Most frequently, the left-trivialized form is used where G i is either SE (3) or SO (3) × R 3 .

Configuration space of an MBS as subvariety
The bodies of an MBS are geometrically constrained.It is assumed for sake of simplicity that the MBS is subjected to a system of m scleronomic geometric (and no additional non-holonomic) constraints of the form 0 = h (g) (91) with the constraint mapping h : G → R m .This accounts for lower or higher pair joints as well as bilateral contacts.The geometric constraints define the MBS configuration space This is an subvariety in the ambient manifold G. Time differentiation of the geometric constraints (91) yields the corresponding velocity constraints.This can be expressed in terms of the left-and right-trivialized derivatives, respectively, as where B l (g) , B r (g) : R 6 n → R m are the Jacobians in vector representation of g.When the semidirect product group is used, then V l is the body-fixed system twist, and V r is the spatial system twist.When the direct product group is used, then V l is the mixed, and V r is the hybrid twist.

Equations of motion in lie group setting
The state space of an MBS is T G ∼ = G × g.The motion equations in absolute coordinate representation can be written where Q ∈ g * represents all Coriolis, centrifugal, potential, and applied forces.This is an index 3 differentialalgebraic system.Using the time derivatives of the velocity constraints (93), it can be resolved in the form This is a standard approach in MBS dynamics [11,143].

Geometric integration of MBS models in absolute coordinates
In the realm of geometric MBS modeling, recent research addressed the solution of the motion equations index 3 system (94) as well of (95) directly on the Lie group.Common to all Lie group integration methods, using canonical coordinates, is that a coordinate vector X i ∈ g is determined so that the configuration at time step i + 1 is given as g i+1 = g i exp X i [74].
In a series of publications, the solution of the differential-algebraic equation (DAE) system (94) has been addressed with amended generalized α schemes that are known from computational structural mechanics.Brüls et al. [32][33][34] introduced generalized α schemes adopting the concept of Lie group integration using canonical coordinates.In [34], the method was extended to MBS containing flexible bodies.It was observed that the direct and semidirect product representations lead to different numerical behavior.A first analysis of this phenomenon was reported in [33].Recently, Arnold et al. [6] presented a consistent theory for the error analysis for generalized α schemes for constrained MBS on Lie groups.
Classical Lie group integration methods were developed for ordinary differential equations.The best known are the Munthe-Kaas methods [74,79,117,118].This integration scheme has been applied to the rigid body problem from the outset [37,57,59,96,101,123].Recently Engø [56] introduced partitioned Runge-Kutta schemes as a special type of Munthe-Kaas integration schemes.In a recent paper [151], Munthe-Kaas method was extended to be used for DAE problems and applied for solving MBS forward dynamics modeled as DAE index 1 in MBS state space, introduced as Lie group.
Krysl [86] proposed a version of a generalized α scheme for the solution of ODE on a Lie group.This was one of the first formlations that initiated the developments for DAE.These ODE integration schemes can be applied to the formulation (95).Maekinen [96] analyzed the applicability of generalized α schemes.
For both approaches, solving the motion equations in the form of the DAE ( 94) or in the ODE form (95), the configuration space can be the direct or the semidirect product Lie group.For the latter, this was investigated in [116].It was concluded that the SE (3) formulation leads to better convergence in particular cases, but both formulations perform equal for general MBS.It will be the subject of future study to investigate the error and convergence of the proposed schemes when using different Lie groups to model the MBS.
Eventually the Lie group integration solves the problem of singularity-free parameterization MBS model.It is common to use Euler parameters as a remedy.Then, however, the unit length constraint must be imposed.In a recent paper [153], a Lie group integration method was proposed that allows for singularity-free integration of the quaternion description without explicitly imposing the unit length constraint.To this end, quaternions are treated as elements of the Lie group Sp (1).In this way, the standard quaternion integration based on a system of four DAEs (which also requires stabilization of the unit length constraint) is transformed to the problem of integrating a system of three ODEs on the Lie algebra sp (1) ∼ = so (3).More precisely, this ODE describes the time evolution of the incremental rotational vector.This integration can be performed by using any standard vector space ODE integration scheme.The Munthe-Kaas method used a Runge-Kutta method, for example.The quaternion is reconstructed with the exponential mapping on the quaternion group Sp (1), which assures satisfaction of quaternion unit length.This approach can be adopted to integration of dual quaternions that describe rigid body motions [114].
In the field of MBS geometric integration, special attention is devoted to structure preserving methods that exploit rich geometric structure of rigid body rotational dynamics (see [7,16,25,58,67,73,84,85,87,89,90,97,104,135,146,147,152] and references cited therein).To this end, rigid body rotational dynamics is studied most conveniently as Lie-Poisson system that is defined on so * (3) (the dual space of so( 3)).Beside free rigid body rotation, Lie-Poisson systems arise in ideal fluid mechanics (Euler equations for an inviscid fluid), heavy top equations [59], Vlasov-Poisson equations, and other settings.
As expressed via non-canonical coordinates, Lie-Poisson formalism is generalization of Hamiltonian formulation that is expressed in canonical form [58,88]. Actually, Lie-Poisson system can be regarded as a reduced system resulting from Lie-Poisson reduction of a canonical Hamiltonian system [97,98].It is linked to Euler-Poincare system (reduced Lagrangian being defined on pertinent Lie algebra) via reduced Legendre transform.In the case of rigid body motions modeled on SE(3), the Euler-Poincare system is introduced by ( 38) and (41).
With the aim of deriving geometric methods for free rigid body rotation, several authors used algorithm that preserves coadjoint orbits in dual space SO * (3) as a departure point toward design of the structure preserving schemes that respect integrals of motion.Therefore, by following [58,84,89,152] a geometric scheme can be introduced that extends the coadjoint orbit-preserving integration method for SO(3).
We start from the Euler equation of free rigid body rotation given as Lie-Poisson system [98] in the form π = − ω b π (96) where ω b is the body angular velocity tensor and π ∈ R 3 , which can be identified with so * (3) [75], represents angular momentum in the body attached frame.By following, for example [89], (96) can be expressed as the coadjoint operator on the dual so * (3) as where ad * is the dual of the ad operator ad a b = ab, which is commutator in the Lie algebra so(3) [79], identified here with R 3 .Within each integration step, the solution of ( 96) can be expressed as an action of SO(3) on R 3 in the form which leads to solving of ODE on the Lie group that reads The update step (98) can be written as coadjoint action [75] , denoted Ad * , of SO(3) on R 3 in the form [58] where Ad * Q(t) π = Q T (t) π is valid and the coadjoint orbit is given as Since for the step initial condition π n the coadjoint orbit O π n is a sphere of radius π n , it follows from (100) that the magnitude of free body angular momentum in the body attached frame is exactly preserved during the step, independently of the accuracy of the integration method for determining Q in (99).As it is explained in [25,58,84,85,89,147,152], by starting from the geometric scheme that preserves coadjoint orbits, further algorithmic upgrades can be made to allow for numerical preservation of additional integrals of motion, such as kinetic energy of the free rotating body or heavy top constants of motion.Another route to derive algorithm that preserves rigid body integrals of motion is treating body as constrained mechanical system, which leads to Hamiltonian formalism in canonical form.To illustrate this, we introduce diagonal matrix D = diag (a 1 , a 2 , a 3 ) with the coefficients defined as θ 1 = a 2 + a 3 , θ 2 = a 1 + a 3 , θ 3 = a 1 + a 2 , a k = B x 2 k dm (x), where θ 1 , θ 2 , θ 3 are eigenvalues of the rigid body inertia tensor.The kinetic energy of the body can be written as [98] T B = we obtain system Hamiltonian in the form where we suppose a potential U (R).Then, the equation of motion for rigid body, modeled as a constrained Hamiltonian system [88,105], can be written in canonical form as where we use the notations ∇ R U = ∂U/∂R i j , ∇ R H = ∂ H/∂R i j , and similarly for ∇ P H . Here, the independent coefficients of the symmetric matrix correspond to the six Lagrange multipliers associated with the constraint equation 0 = R T R−I.The corresponding velocity constraint equation reads as 0 = R T Ṙ+ ṘT R, which yields R T PD −1 + D −1 P T R = 0.These two constraints imply that the equations (105) constitute a Hamiltonian system constrained on the manifold which, however, is not the cotangent bundle T * SO (3) of the manifold SO(3).Discretization of ( 105) that is based on standard Störmer-Verlet integration algorithm [73,74,136] leads to symplectic RATTLE scheme for rotational rigid body dynamics that exhibits excellent structure preserving properties [88].structure preserving scheme for rigid body dynamics that also treats rigid body as constrained mechanical system is proposed in [13].In the recent paper [152], Störmer-Verlet integration algorithm is a departure point for derivation of another structure preserving integration algorithm for rigid body rotational dynamics.This algorithm is formulated in Lie group setting and exactly preserves free rotating body angular momentum and kinetic energy, along an excellent approximation of heavy top integrals of motion.
Future research will extend Lie group integration schemes to MBS with flexible bodies in particular for flexible bodies undergoing large deformations.A Cosserat rod, for instance, is a curve in a Lie group.

Conclusions
The motion of a MBS evolves on a Lie group.It is hence most obvious to describe MBS models in a Lie group.This poses new challenges but at the same time opens up new opportunities.It was the aim of this review article to summarize these concepts and to point out some of the major issues.It has been shown that from modeling perspective, Lie group methods offer more flexibility compared to classical modeling approaches.It has also been discussed that it will be imperative to pursue research on the numerical integration using Lie group schemes.This must, in particular, include flexible bodies undergoing large deformations.In this regard, Lie group formulations were already (at least implicitly) used in [18][19][20][21][22][23]63] .
Acknowledgements Open access funding provided by Johannes Kepler University Linz.Andreas Müller acknowledges that this work has been partially supported by the Austrian COMET-K2 program of the Linz Center of Mechatronics (LCM).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/),which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix: Notations
Let G be a Lie group.The Lie algebra, denoted with g, is identified with the tangent space of G at the group identity, equipped with the Lie bracket of left-invariant vector fields.In case of matrix groups, the Lie bracket is the matrix commutator.The left and right translation mappings is denoted with L g and R g , respectively.For a matrix group G, they are L g (h) = gh and R g (h) = hg, for gh ∈ G. Let exp : g → G be the exponential mapping on G.The right-trivialized differential of the exp mapping dexp X : g → g is defined via with the tangent mapping T X exp (Y ) = d dt t=0 exp (X + tY ) : g → T exp x G. On a matrix Lie group G ⊂ G L (n), for g (t) = exp X (t), the following relation holds dexp X Ẋ = ġg −1 (109) dexp −X Ẋ = g −1 ġ. (110)

3 )
is a matrix Lie group.Throughout the paper, the notations C ∈ SE (3) and C ∈ SE (3) are used interchangeably.

Table 1
Subgroups of SE(3) and the corresponding lower pair joint 2 Planar joint + screw joint with axis normal to plane 3 S O(2) R 2 = SE (2) Planar joint 4 S O(2) R 3 = SE (2) R Planar joint + prismatic joint with axis normal to plane 6 S E(3) 'Free joint'