Screw and Lie group theory in multibody dynamics

Screw and Lie group theory allows for user-friendly modeling of multibody systems (MBS), and at the same they give rise to computationally efficient recursive algorithms. The inherent frame invariance of such formulations allows to use arbitrary reference frames within the kinematics modeling (rather than obeying modeling conventions such as the Denavit–Hartenberg convention) and to avoid introduction of joint frames. The computational efficiency is owed to a representation of twists, accelerations, and wrenches that minimizes the computational effort. This can be directly carried over to dynamics formulations. In this paper, recursive O(n)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$O ( n ) $\end{document} Newton–Euler algorithms are derived for the four most frequently used representations of twists, and their specific features are discussed. These formulations are related to the corresponding algorithms that were presented in the literature. Two forms of MBS motion equations are derived in closed form using the Lie group formulation: the so-called Euler–Jourdain or “projection” equations, of which Kane’s equations are a special case, and the Lagrange equations. The recursive kinematics formulations are readily extended to higher orders in order to compute derivatives of the motions equations. To this end, recursive formulations for the acceleration and jerk are derived. It is briefly discussed how this can be employed for derivation of the linearized motion equations and their time derivatives. The geometric modeling allows for direct application of Lie group integration methods, which is briefly discussed.


Introduction
The core task in computational multibody system (MBS) dynamics is to either construct the equations of motion (EOM) explicitly, that can be written for an unconstrained tree-topology MBS in the form M (q) q + C ( q, q) q = Q ( q, q, t) , in a way that is easy to pursue, or to evaluate them for given (q, q, q) and t, respectively to solve them, in a computationally efficient way for q (t).In continuation of [62] the aim of this paper is to present established O (n) formulations in a common geometric setting and to show that this setting allows for a flexible and user-friendly MBS modeling.Screw and Lie group theory provides a geometric framework that allows for achieving optimal computational performance and at the same time allows for an intuitive and flexible modeling.In particular, it gives rise to a formulation of the MBS kinematics that does not involve body-fixed joint frames.The kinematics modeling is indeed reflected in the formulation used to evaluate the EOM.A central concept is the representation of velocities (twists) as screws.Four different variants were recalled in [62].In this paper their application to dynamics modeling is reviewed.A wellknown approach, which exploits the fact that rigid body twists are screws, is the so-called 'spatial

Body-Fixed Representation
Starting from (101) the body-fixed acceleration is Vb i = J b i q + Jb i q, and explicitly in terms of the body-fixed instantaneous screw coordinates Using the matrix form of (103) the partial derivatives of the instantaneous screw coordinates are This can be evaluated with help of the POE formula (93) as and in the same way follows that Inserted into (2) yields , and noting (113), the final expression is Hence the partial derivative of the instantaneous joint screw J b i,j w.r.t. to q k is simply the screw product (114) of J b i,j and J b i,k .The final expression for the acceleration attains a very compact form Indeed the same result would be obtained using (103) in terms of Y i .This expression has been derived, using different notations, for instance in [16,51,65,69].
The equations ( 6) can be summarized for all bodies i = 1, . . ., n using the system twist (111) and system Jacobian (112).To this end, the derivative ( 5) is rewritten as Noticing that adk X k J b k,k = 0 the time derivative of the body-fixed system Jacobian factors as Jb (q, q) = −A b (q) a b ( q) A b (q) X b = −A b (q) a b ( q) J b (q) (8) with A b defined in (24) of [62] and with a b ( q) := diag ( q1 ad1 X1 , . . ., qn adn Xn ).
Hence the system acceleration is given in compact matrix form as Vb Remark 1 (Overall inverse kinematics solution) The relation (10) gives rise to a solution of the inverse kinematics problem on acceleration level, i.e. the generalized accelerations for given configurations, twists, and accelerations of the bodies.The unique solution is which is indeed the time derivative of (26) in [62].In components this gives the acceleration of the individual joints as qi A further time derivative of the twist yields the jerk of a body, which requires a further partial derivative of the Jacobian.Starting from (5), and using the Jacobi identity (116) and the bilinearity ], the non-zero second partial derivative is found as This gives rise to an explicit form for the body-fixed jerk ] qj qk qr .
Thus only computationally simple nested screw products are required to compute the terms that are quadratic and cubic in qj , qk .The same applies to higher derivatives (that for instance are required for motion planning).The explicit form of the ν-th order partial derivative was presented in [57] ∂ ν J b i,j where is the ordered sequence of the indices α 1 , . . ., α ν .Clearly the closed form expressions become very involved.Their explicit determination can be avoided by recursive evaluation [57].

Spatial Representation
Proceeding in the same way as for (2) the partial derivative of the spatial Jacobian is obtained as Since the spatial representation J s j is intrinsic to the joint j, rather than related to a body as is (103), the time derivative can be expressed as This relation reconfirms the special properties of spatial twists that are advantageous for recursive implementations.It may be considered as an extension of Euler's formula for the time derivative of vectors resolved in moving frames to screws.For this reason Featherstone [28,30] termed the Lie bracket the 'spatial cross product'.The spatial acceleration is therewith In matrix form the overall spatial acceleration can be summarized as and L being the lower triangular block identity matrix.A solution for q similar to (10) exists.
The second partial derivative of the spatial Jacobian is Therewith the spatial representation of the jerk of body i is found as The instantaneous joint screws (104), and thus their derivatives (15) and (20), are independent of a particular body.The closed form of the ν-th order partial derivative has been reported in [55] is the ordered sequence of the indices α 1 , . . ., α ν .The last form in (23) allows for a recursive determination.Moreover, a recursive formulation for the time derivative of spatial twists has been reported in [58].Together with the very concise form (16) this makes the spatial representation computationally very attractive.

Hybrid Form
The results in section 2.1 can be carried over to the hybrid twist making use of the relation (106).
As in (118), denote with v J h i,k and ω J h i,k the screw coordinate vectors comprising respectively the linear and angular part of the column of the hybrid Jacobian so that and thus The similarity to (5) is apparent.The difference is that the convective term due to the angular motion is missing, which is why only v J appears.The time derivative of the hybrid Jacobian can thus be expressed as Jh where ∆V h j−1,i := V h i − Ad ri,j−1 V h j−1 is the relative hybrid twist of body j − 1 and i as observed in the BFR on body i.A simpler relation is obtained by directly differentiating (105)
This yields the following explicit expressions for the hybrid acceleration For the second derivative it is simplest to start from (27), and a straightforward calculation yields Jh i,j = ad ri,j + 2ad ṙi,j ad ω s j + Ad ri,j (ad ωs j + ad ω s j ad ω s j ) 0 X j j .
The jerk in hybrid representation can thus be written as ... q j + 2ad ṙi,j qj + ad ri,j + 2ad ṙi,j ad ω s j qj + Ad ri,j 2ad ω s j qj + ad ωs j + ad ω s j ad ω s j qj ) 0 X j j .

Mixed Representation
With (100), employing the results for the mixed representation, yields 3 Newton-Euler Equations in Various Representations

Spatial Representation
Consider a rigid body with body-fixed BFR F b = {Ω; e b,1 , e b,2 , e b,3 } located at an arbitrary point Ω.Denote the inertia matrix w.r.t.this BFR with M b , ref. (40).The configuration of the BFR F b is described by C = (R, r).The spatial inertia matrix expressed in the IFR is then The spatial canonical momentum co-screw Π s = (L s , P s ) T ∈ se * (3), conjugate to the spatial twist, is thus The momentum balance yields the Newton-Euler (NE) equations in spatial representation, which attains the simple form Πs = W s where W s = (t s , f s ) T is the applied wrench, with spatial torque t s ≡ 0 t 0 and force f s ≡ 0 f , both measured and resolved in the IFR.The momentum balance equation (36) is the simplest form possible, which is achieved by using the spatial representation of twist, wrench, and momentum.Firstly, it does not involve any vectorial operation, e.g.cross products.Secondly, it is also numerically advantageous: any numerical discretization of the ODE (36) easily preserves the spatial momentum in the absence of external wrenches.This has been discussed already by Borri in [13].
In this context the spatial formulation is called the fixed pole equation.In a recent paper [32] the advantages of this form are exploited for geometrically exact modeling of beams.
The explicit and compact form in terms of the spatial twist is found, introducing (35) and using along with ad V s V s = 0, as Remark 2 Writing (38) as This property is called the skew symmetry of the motion equations [65].

Body-fixed Representation
Let F c = {C; e c,1 , e c,2 , e c,3 } be a body-fixed frame located at the COM.Its configuration is described by The inertia matrix w.r.t.this COM frame is denoted with (38) to the BFR F b yields the body-fixed momentum balance represented in F b in the concise form with the applied wrench W b = t b , f b T in body-fixed representation.The equations (41) are formally identical to the spatial equations (38).Written separately, this yields the NE equations expressed in an arbitrary body-fixed BFR When using the COM frame as special case, the momentum represented in the body-fixed COM frame is , and the momentum balance yields Written in components, this yields the NE equations represented in the COM frame Noticeably the angular and translational momentum equations are coupled even though the COM is used as reference.This is due to using body-fixed twists.

Hybrid Form
The hybrid twist V h c = (ω s , ṙc ) T of the COM frame is related to the body-fixed twist by Ad −1 Rc V b c , see (98), where R c is the absolute rotation matrix of F c in C c .The hybrid momentum screw is thus , where the hybrid representation of the inertia matrix is The hybrid momentum balance w.r.t. the COM follows from Πh with ω V h c = (ω s , 0) T (notice −ad T ω s = ad ω s ).Writing (48) separately for the angular and linear momentum balance (50) shows that the hybrid NE equations w.r.t. the COM are indeed decoupled.Here W h c = (t h c , f h c ) T denotes the hybrid wrench measured in the COM frame and resolved in the IFR.Now consider the arbitrary body-fixed BFR F b with configuration C = (R, r).The hybrid twist V h = (ω s , ṙ) T measured at this RFR is V h = Ad d bc V h c , with the displacement vector d bc from BFR to COM resolved in the IFR.The hybrid mass matrix w.r.t. to the BFR F b is found as The momentum balance in hybrid representation w.r.t. an arbitrary BFR is found, using Ȧd −1 (52), as Separating the angular and translational part results in These are simpler than the body-fixed equations ( 42) and (43).Finally notice that f h = f s .

Mixed Form
The mixed twists V m = ω b , ṙ T consists of the body-fixed angular velocity ω b , i.e. measured and resolved in the BFR F b , and the translational velocity ṙ measured at the BFR F b and resolved in the IFR.The NE equations for the mixed representation w.r.t. a general BFR are directly found by combining ( 42) and ( 55), with r = vb + ω b v b , If a COM frame is used, combining ( 45) and ( 50) yields

Arbitrary Representation
The NE equations of body i represented in an arbitrary frame F j are obtained by a frame transformation of the spatial momentum balance (36) as The spatial twist in terms of the twist of body i represented in with the inertia matrix of body i represented in frame j The spatial and body-fixed representations are special cases with i = j.Even more generally, the NE equations can be resolved in yet another frame F k .This is achieved by transforming the momentum balance (36) as where R k,j is the rotation matrix from F i to F k .The final equations follow from ( 60) and the with the mass matrix of body i measured at frame F j and resolved in frame The spatial and body-fixed representations are special cases with i = j = k, and the hybrid representation with i = j and k = 0.An alternative form of the NE equations in arbitrary reference frames was presented in [9].

Recursive Evaluation of the Motion Equations for a Kinematic Chain
The model-based control of complex MBS as well as the computational MBS dynamics rely on efficient recursive inverse and forward dynamics algorithms.A recursive Newton-Euler method for tree-topology MBS was presented in an abstract, i.e. coordinate-free, approach in [47].However, the various recursive methods using different representations give rise to algorithmically equivalent methods but with different computational costs.In the following, the various inverse dynamics algorithms are presented and their computational effort is estimated.A detailed analysis as well as the forward dynamics algorithms are beyond the scope of this paper.The presented discussion is nevertheless indicative also for the corresponding forward dynamics algorithms.Some results on the forward kinematics complexity can be found in [67,81,87].This depends on the actual implementation, however.A comparative study is still due, and shall be part of further research.
The inverse dynamics consists in evaluating the motion equations for given joint coordinates q, joint rates q, accelerations q, and applied wrenches W app i , and hence to determine the joint forces Q = (Q 1 , . . .Q n ).The starting point of recursive algorithms for rigid body MBS are the NE equations of the individual bodies.The MBS dynamics is indeed governed by the Lagrange equations.Consequently, summarizing the recursive steps yields the Lagrangian motion equations in closed form.This will be shown in the following.It is assumed for simplicity that the inertia properties, i.e. the mass matrices M b i , are expressed in the body-fixed BFR of body i determining its configuration, rather than introducing a second frame.

Body-fixed Representation
Forward Kinematics Recursion Given the joint variables q, the configurations of the n bodies are determined recursively by (92) or (93), and the twists by (108).Then also the accelerations are found recursively.The expression

and hence Vb
where (65b) and (65c) follow by replacing either argument in the Lie bracket using (108).
Remark 3 Notice that solving (108) for qi leads to the result in remark 9 of [62].Solving (65c) for qi yields (11).Using (65b) the latter can be expressed as qi Recursive Newton-Euler Algorithm Once the configurations, twists, and accelerations of the bodies are computed with the forward kinematics recursion, the Newton-Euler equations (41) for each individual body can be evaluated by an inverse dynamics backward recursion.The momentum balance of body i then yields the resulting body-fixed wrench W b i acting on the body due to generalized joint forces and constraint reactions forces.Projecting the resultant wrench onto the screw axis i X i of joint i yields the generalized force Q i .Summarizing the forward and backward recursions yields the following recursive algorithm: Forward Kinematics -Input: q, q, q -For i = 1, . . ., n The joint reaction wrench is omitted in (67a) since this is reciprocal to the joint screw, and does not contribute to (67b).Notice that, with (94), the body-fixed i X i as well as the spatial representation Y i of joint screw coordinates can be used.This form of the recursive body-fixed NE equations, using Lie group notation, has been reported in several publications [69,73,72,51].
Computational Effort For the kinematic chain comprising n bodies connected by n 1-DOF joints, in total, the twist recursion (66b) and acceleration recursion (66c) each requires n − 1 frame transformations.The acceleration recursion (66c) further requires n − 1 Lie brackets.The second argument of the Lie bracket can be reused from (66b).Hence the twist and acceleration recursion need 2 (n − 1) frame transformations and n − 1 Lie brackets.The backward recursion (67a) needs n − 1 frame transformations and n Lie brackets.In total, the NE algorithm needs 3 (n − 1) frame transformations and 2n−1 Lie brackets.The evaluation of the Lie bracket in (66c) can be simplified using (65b) since the screw vector i X i expressed in RFR is sparse and often only contains one non-zero entry.
Remark on Forward Dynamics Using the body-fixed representation, a recursive forward dynamics algorithm, making explicit use of Lie group concepts, was presented in [69,70,73,72,78].The kinematic forward recursion together with the factorization in section 3.1.3of [62] was used derive O (n) forward dynamics algorithms in [31,45], where the Lie group concept is regarded as spatial operator algebra.Other O (n) forward dynamics algorithms were presented in [2,7,35,36].The inverse dynamics formulation was also presented in [30,46] in the context of screw theory.

Spatial Representation
Forward Kinematics Recursion Expressing the spatial twist in terms of the spatial Jacobian, the expressions (104) lead immediately to The recursive determination of spatial accelerations thus only requires the time derivative ( 16) of the spatial Jacobian, so that The second form in (69) follows by inserting (68).This is the generalization of Euler's theorem, for the derivative of vectors resolved in moving frames, to screw coordinate vectors.Therefore the ad operator is occasionally called the 'spatial cross product'.
Recursive Newton-Euler Algorithm The momentum balance expressed with the spatial NE equations (38) together with (68) leads to the following algorithm: Forward Kinematics -Input: q, q, q -For i = 1, . . ., n In contrast to (66b), once the instantaneous screws (70b) and the spatial mass matrix (34) are computed, the recursions (70c), (70d), and (71b) do not require frame transformations of twists.Instead the spatial mass matrix is transformed, according to (71a), which is the frame transformation of a second-order tensor.Overall the spatial algorithm needs n frame transformations of screw coordinates, n frame transformation of a second-order tensor, and 2n − 1 Lie brackets.Comparing body-fixed and spatial formulation, it must be noticed that the frame transformation of the second-order inertia tensor has the same complexity as two screw coordinate transformations (if just implemented in the form (34)), and hence the computational complexity of both would be equivalent.This fact is to be expected since body-fixed and spatial representations are related by frame transformations.Nevertheless the spatial version has some interesting features that shall be emphasized: 1.The NE equations (38) form a non-linear first-order ODE system on SE (3) × se (3).Since a spatial reference is used, the momentum conservation of a rigid body can simply be written as Πs i = 0, where Π s i ∈ se * (3) is the momentum co-screw.Using the spatial momentum balance (36) has potentially two advantages.Firstly, ( 36) is a linear ODE in Π on the phase space SE (3) × se * (3).This implies that a numerical integration scheme can easily preserve the momentum, as pointed out in [13].Secondly, O (n) formulations using canonical momenta have been shown to be computationally advantageous.An O (n) forward dynamics algorithm based on the canonical Hamilton equations was presented in [66] that uses on the hybrid form.It was shown to require less numerical operations than O (n) algorithms based on the NE equations.It is also known that O (n) algorithms based on the spatial representation can be computationally more efficient than those based on body-fixed or hybrid representations [30].A further reduction of computational costs shall be expected from an algorithm using spatial momenta.2. It is interesting to notice that the hybrid as well as the spatial twists appear in the recursive O (n) forward dynamics formulation in [6], where the first is called 'Cartesian velocity' and the latter 'velocity state'.In this formulation the spatial twist plays a central role, and it was already remarked that the recursive relation of spatial twists, ref.
3. If a purely kinematic analysis is envisaged the forward recursion (70b)-(70d) is more efficient than the body-fixed and the hybrid version (see next section) [67] (disregarding possibly necessary transformations of the results to local reference frames).As pointed out in section 2.2 this advantage is retained for the higher-order kinematics (jerk, jounce, etc.) [55].
Remark on Forward Dynamics The spatial formulation is rarely used for dynamics.Featherstone [27,30] derived a forward dynamics O (n) algorithm.It was concluded that this requires the lowest computational effort compared to other methods.But this does not take into account the necessary transformations of twists and wrenches to local reference frames.Moreover, it was shown in [81] that the O (n) forward dynamics algorithm in body-fixed representation, using the body-fixed joint screw coordinates i X i and RFR at the joint axis, can be implemented in such a way that it requires less computational effort than the spatial version.The key is that when the BFR of F i is located at and aligned with the axis of joint i, then i X i becomes sparse.From a users perspective this is a restraining presumption, however.

Hybrid Form
Forward Kinematics Recursion The hybrid twists is determined recursively by (110) with 0 Taking into account that ad ṙi V h i − X h i qi = ad ṙi V h i−1 (because there is no angular part in ad ṙi ), and ad ṙi + ad ω s i = ad V h i , this can be transformed to Another form follows by solving (110) for V h i−1 and inserting this into (72), while noting that ad ṙi,i−1 Ad −1 ri,i−1 = ad ṙi,i−1 , as Comparing these three different recursive relations (72), (73), and (74) for the hybrid acceleration from a computational perspective (72) is the most efficient.
Recursive Newton-Euler Algorithm With the hybrid Newton-Euler equations ( 53) the recursive NE algorithm is as follows: Forward Kinematics The hybrid representation is a compromise between using twists and wrenches measured in body-fixed frames (as for the body-fixed representation, where twists and wrenches are measured at the RFR origin) and those resolved in the IFR (as for the spatial representation, where twists and wrenches are measured at the IFR origin).It has therefore been used extensively for O (n) inverse and forward dynamics algorithms.The essential difference between the forward recursion for kinematic evaluation in body-fixed and hybrid formulation is that the body-fixed recursion (66a-66c) requires frame transformations of screws involving rotations and translations whereas the hybrid recursion (75a-75d) only requires the change of reference point using position vectors resolved in the IFR.The attitude transformation only appears in (75b) and in the computation of the hybrid inertia matrix (76a).In total the forward kinematics needs n rotational transformations, 2n − 2 translational transformations.Further, (75d) needs n − 1 cross products of the form ad ṙi,i−1 The inverse dynamics needs the n rotational transformations (76a) of the second-order inertia tensor, n − 1 translational transformations of wrenches and n Lie brackets with ω s i in (76b).In total the hybrid NE algorithm needs 3n − 3 translational and n rotational transformations of screw coordinates, n rotational transformations of the inertia tensor, and 3n − 1 Lie brackets.Although the number of operations is equivalent to the body-fixed version the particular form of transformations is computationally very simple motivating its extensive us in O (n) forward dynamics algorithms.Moreover, the hybrid NE equations are commonly expressed in a body-fixed BFR at the COM, so that the hybrid NE equations simplify to (48) and ( 49), (50), respectively.Instead of transforming the joint screws i X i or Y i in the reference configuration, the instantaneous hybrid joint screws can be determined using the defining expression (36) in [62] with the current b j,i and e j .
Remark on Forward Dynamics The above inverse dynamics formulation was presented in [38,39,75,76] together with O (n) forward dynamics algorithms.An O (n) forward dynamics method was presented in [1,6].These algorithms are deemed efficient taking into account that the computation results do not have to be transformed to the body-fixed reference points of interest, as in case of the spatial version.An O (n) forward dynamics algorithm was developed in [66] using canonical Hamilton equations in hybrid representation, i.e. the momentum balance (52) in terms of the conjugate momenta Π h i , rather than the NE equations.It was concluded that its performance is comparable to that of Featherstone's [30] method in terms of spatial twists.

Choice of Body-Fixed Reference Frames
The Lie group formulation involves geometric and inertia properties that are readily available, e.g. from CAD data.In [62] and in the preceding sections two approaches to the description of MBS geometry (with and without body-fixed joint frames) and three versions for representing velocities and accelerations (body-fixed, spatial, hybrid) were presented, of which each has its merits.The description of the geometry is independent from the representation of twists.For instance, the geometry could be described in terms of joint screws expressed in the IFR while the kinematics and dynamics is modeled using body-fixed twists.This allows to take advantage of the low-complexity hybrid or spatial recursive NE equations while still having the freedom to use or avoid body-fixed joint frames.
The standard approach to model an MBS is to introduce 1.) an IFR, 2.) body-fixed BFRs, and 3.) body-fixed JFRs.The latter is avoided using spatial joint screws Y i , as already presented.It still remains to introduce body-fixed BFRs kinematically representing the bodies.However, even the explicit definition of RFRs can be avoided by properly placing them.Their location is usually dictated by the definition of the inertia tensors, and it is customary to relate the inertia data to the COM.If instead the body-fixed BFRs are assigned such that they coincide in the reference configuration (q = 0) with the IFR, then no reference configurations of bodies need to be determined (A i = I).This normally means that the RFR is outside the physical extension of a body.That is, the inertia properties of all bodies are determined in the assembly reference configuration w.r.t. to the global IFR.In other words they are deduces from the design drawing (corresponding to q = 0) relative to a single construction frame.This can be exploited when using CAD systems.The required kinematic data then reduces to the direction and position vectors, e i and y i , in order to compute Y i in (94).As a side effect it is Y i = i X i .This is an important result, and does apply to any of the discussed twist representations, since the representation of twists has nothing to do with the geometry description.Moreover, then the POE (93) and the Jacobian (103), (104), and thus (106), in terms of spatial screw coordinates simplify.The only computational drawback is that the hybrid Newton and Euler equations are not decoupled since the spatial IFR, to which the inertia data is related, is unlikely to coincide with the COM of the bodies in the reference configuration.Details can be found in [54].

Euler-Jourdain Equations
The body-fixed NE equations for the individual bodies within the MBS are where W b,c i is the constraint reaction wrench of joint i, and W b,app i represents the total wrench applied to body i including the applied wrench in joint i.Jourdain's principle of virtual power, using the admissible variation δV b = J b δ q of the system twist, and noting that δV b are reciprocal to the constraint wrenches (see appendix B.2), yields the system of n motion equations This form allows for a concise and computationally efficient construction of the motion equations.
The point of departure are the NE equations of the individual bodies.The body-fixed system Jacobian (112) is determined by the (constant) joint screw coordinates in X b and the screw transformations encoded in A b .The same applies to the other representations.The accelerations are determined by (6), respectively (10).Explicit evaluation of (78) leads to the recursive algorithm in section 4.1.Inserting the twists and accelerations in (78) yields the equations (1) that determine the MBS dynamics on the tangent bundle T V n with state vector (q, q) ∈ T V n .Alternatively, combining ( 78) with (111) yields a system of n + 6n ODEs in the state variables q, V b ∈ V n × se (3) n that govern the dynamics on the state space V n × se (3) n .The advantage of this formulation is that it is an ODE system of first order, and that the system has block triangular structure.Yet another interesting formulation follows with the NE equations (36) in terms of the conjugate momenta in spatial representation M s i J s i q = Π s i , i, . . ., n.This is a system of n + 6n first order ODEs in the phase space (q, Π s ) ∈ V n × se * (3) n .The system (79) can be solved for the Πs i and qi noting the block triangular structure of J s [62].From a numerical point of view the momentum formulation in phase space shall allow for momentum preserving integration schemes.Various versions of ( 78) have been published.Using the hybrid representation of twists, basically the same equations were reported in [4].There the system Jacobian is called the 'natural orthogonal complement' motivated by the fact that the columns of J b are orthogonal to the vectorial representations of constraint wrenches (although the former are screws while the latter are coscrews).In classical vector notation they were reported in [40,49] and [15].In [40] the equations ( 78) are called Euler-Jourdain equations.In [49], emphasizing on the recursive evaluation of the body Jacobian the instantaneous body-fixed joint screws J b i are called 'kinematic basic functions' as they are the intrinsic objects in MBS kinematics.In [15] the equations ( 78) are called 'projection equations' since the NE equations of the individual bodies are restricted to the feasible motion (although J b is not a projector).The equations ( 78) in body-fixed representation are equivalent to Kane's equations where J b i are called 'partial velocities' [41].The instantaneous joint screw coordinates, i.e. the columns J b i of the geometric Jacobian, were also called 'kinematic influence coefficients' and their partial derivatives (5) the 'second-order kinematic influence coefficients' [8,84].It should be finally remarked that due to the block triangular form of J b , solving (78), and using the inversion of A b (ref.(25) in [62]), leads immediately to an O (n) forward dynamics algorithm.This is the common starting point for deriving forward dynamics algorithms that applies to any twist representation.

Lagrange Equations
The MBS motion equations can be derived as the Lagrange equations in terms of generalized coordinates.For simplicity, potential forces are omitted so that the Lagrangian is simply the kinetic energy.Then the equations attain the form with generalized mass matrix M, and C ( q, q) q representing Coriolis and centrifugal forces.The vector Q stands for all other generalized forces, including potential, dissipative, and applied forces.Using body-fixed twists, the kinetic energy of body i is qT M q with the generalized mass matrix and Its time derivative is with (10) given as and a b defined in (9).From (7) follows This admits to identify the generalized mass matrix (81) and the matrix The first term on the right hand side in ( 83) can be simplified so that The concise expressions ( 81) and ( 84) allow for construction of the Lagrange equations in closed form.Similar expressions can be derived using the spatial and hybrid representation of twists.
For analytic investigations of the MBS dynamics it may be useful to write the Lagrange equations in components as where the Christoffel symbols of first kind are defined as The recursive relations (5) give rise to the closed form expressions with i < j ≤ k or j ≤ i < k, r = max (i, j) , s = min (i, j) .
This expression for the Christoffel symbols in Lie group notation was reported in [17,51], and already in [49] in tensor notation.This expression simplifies when Binet's inertia tensor ϑ i = with i < j ≤ k or j ≤ i < k.
The equations (86) were presented in [51,70,72,69,73], and (87) in [51].Prior to these publications the equations ( 86) and ( 87) have been reported in [49,50] using tensor notation rather than Lie group notation.Another publication that should be mentioned is [20] where the Lagrange equations were derived using similar algebraic operations.The above closed forms of EOM are derived using body-fixed twists.The potential benefit of using spatial or hybrid twists remains to be explored.

Derivatives of Motion Equations
In various contexts the information about the sensitivity of the MBS kinematics and dynamics is required either w.r.t.joint angles, geometric parameters, or dynamic parameters.Whereas it is know that the EOM of a rigid body MBS attain a form that is linear in the dynamic parameters, they depend non-linearly on the generalized coordinates and geometry.The POE formulation provides a means to determine sensitivity w.r.t. to kinematic parameters.

Sensitivity of Motion Equations
Gradients w.r.t.generalized coordinates are required for the linearization of the EOM (as basis for stability analysis and controller design) as well as for optimal control of MBS.Since the secondorder and higher derivatives (12) of the body-fixed Jacobian, (20) of the spatial Jacobian, and (25) of the hybrid Jacobian are given as algebraic closed form expressions in terms of screw products, the linearized EOM can be evaluated recursively as well as expressed in closed-form.Using the Lie group notation this was reported in [78].
The same results were already presented in [50] using tensor notation.Comparing the two formulations reveals once more that the matrix Lie group formulation provides a level of abstraction leading to compact expressions.A closed form for partial derivatives of the inverse mass matrix has been reported in [52], which is required for investigating the controllability of MBS.Using the body-fixed representation of twists, recursive O (n) where reported in [36,3].

Geometric Sensitivity
Optimizing the design of an MBS requires information about the sensitivity w.r.t. to geometric parameters.A recursive algorithm was reported in [36] and its parallel implementation in [3] where the partial derivatives are computed on a case by case basis.The Lie group formulation gives rise to a general closed-form expression.To this end, the POE formula (92) is extended as follows.
The geometry of the two bodies i and i − 1 connected by joint i is encoded in the constant part S i,i and S i−1,i in (91), respectively in B i in the formulation in (92).These are frame transformations, and can hence be parameterized in terms of screw coordinates.If B i depends on λ ≤ 6 geometric parameters, it is expressed as The screw coordinates U i1 and corresponding parameters π i1 account for the considered variations from the nominal geometry, represented by B i0 ∈ SE (3).The relative configuration due to joint i and the geometric variations is thus The key observation is that partial derivatives of B i (π i ) are available in closed form, as for the joint screw coordinates.Hence also the sensitivity w.r.t. the MBS geometry can be expressed in closed form [53].This fact has been applied to robot calibration [22,23] where the POE accounts for geometric imperfections to be identified.

Time Derivatives of the EOM
The design of feedback-linearizing flatness-based controllers for robotic manipulators that are modeled as rigid body MBS actuated by elastic actuators (so-called series elastic actuators) requires the time derivatives of the inverse dynamics solution Q (t) [26,68].That is, the first and second time derivatives of the EOM are necessary.Extensions of the classical recursive Newton-Euler inverse dynamics algorithms in body-fixed representations were presented in [19].As it can be expected the relation are very complicated.Using the presented Lie group formulation of the inverse dynamics algorithms gives rise to rather compact and thus fail-safe algorithm.This was presented in [63] for the body-fixed and hybrid version.

Geometric Integration
This paper focussed on the MBS modeling in terms of relative (joint) coordinates.Alternatively, the MBS kinematics can be described in terms of absolute coordinates.One of the issues that is being addressed when modeling MBS in terms of absolute coordinates is the kinematic reconstruction, i.e. the determination of the motion of a rigid body, represented by C (t), from its velocity field V (t).This amounts to solving one of the equations (see appendix A2 in [62]) together with the NE (41) or (38), respectively.Classically, the orientation is parameterized with three parameters.The problem encountered is that there is no singularity-free global parameterization of rotations with three parameters.Instead of local parameters (position and rotation angles) the absolute configurations of the rigid bodies within the MBS can be represented by C (t).Then a numerical integration step from time t k−1 to t k = t k−1 + h shall determine the incremental configuration update . The equations (88) are ODEs on the Lie group SE (3).These can be replaced by ODEs on the Lie algebra se (3).The motion increment from t k−1 to t k is parameterized as ∆C (t) = exp X (t) with an algorithmic instantaneous screw coordinate vector X.Then (88) are equivalent to the ODEs on the Lie algebra where dexp X : se (3) → se (3) is the right-trivialized differential of the exp mapping on SE (3) [13,60,61,71].This is the basic idea of the class of Munthe-Kaas integration schemes [21,37,64].This scheme has been adapted to MBS in absolute coordinates [82].The advantage of these integration methods is that no global parameterization is necessary since the numerical integration is pursued in terms of the incremental parameters X.The ODEs (89) can be solved with any vector space integration scheme (originally the Munthe-Kaas scheme uses a Runge-Kutta method) with initial value X (t k−1 ) = 0.
Recently the geometric integration concepts were incorporated in the generalized α method [42,18] for MBS described in absolute coordinates.In this case the representation of proper rigid body motions is crucial as discussed in [60,59], which is frequently incorrectly represented by SO (3) × R 3 .Also momentum preserving schemes were proposed [83].It should be mentioned that the concept of geometric integration schemes on SE (3) can be transferred to the kinematics of flexible bodies undergoing large deformations described as Cosserat continua.In this context the spatial description (referred to as fixed pole formulation) has proven to be beneficial [32].Recent results on Lie group modeling of beams can be found in [79,80].

Conclusions and Outlook
The computational effort of recursive O (n) algorithms, but also of the formalisms for evaluating the EOM in closed form, depends on the representation of rigid body motions and of the motions of technical joints.Since the geometry of finite rigid body and relative motions is described by the Lie group SE (3), and that of instantaneous motions be the screw algebra se (3), Lie group theory provides the geometric framework.As already shown in [62], Lie group formulations for the MBS kinematics give rise to compact recursive formulations in terms of relative coordinates.In this paper the corresponding recursive NE algorithms were presented and related to the various O (n) algorithms scattered in the literature.This allows for a comparative investigation of their efficiency in conjunction with the modeling procedure.For instance, whereas most O (n) algorithms used the hybrid representation, the spatial representation, as used by Featherstone [30] and Bottasso [13] (where it is called fixed point formulation), is receiving increased attention since it gives easily rise to structure preserving integration schemes [11,12,13,32].A conclusive investigation will be the subject of future research.Future research will also focus on combining the O (n) forward dynamics algorithm by Featherstone [30], based on NE equations using spatial representations with Naudet's algorithm [66] based on Hamilton's canonical equations in hybrid representation.The use of the spatial momentum balance shall allow for momentum preserving integration of the EOM and at the same time to reduce the number of frame transformations.A further important research topic is the derivation of structure preserving Lie group integration schemes for which the spatial formulation of EOM will be formulation of choice.

A Summary of basic Kinematic Relations
As prerequisite the kinematic relations derived in [62] are summarized.Denote with the absolute configuration of body i w.r.t. the inertial frame (IFR) F 0 .This is alternatively denoted with C i = (R i , r i ).The relative configuration of body i relative to body i − 1 is given as where is the reference configuration of body i w.r.t.body i − 1, i.e. for q i = 0, and i−1 Z i ∈ R 6 is the screw coordinate vector of joint i represented in the joint frame (JFR) J i−1,i on body i − 1. Successive relative configurations can be combined to where i X i ∈ R 6 is the screw coordinate vector of joint i represented in the joint frame fixed at body i, Y i ∈ R 6 is the joint screw coordinate vector in spatial representation (measured and resolved in IFR) for the reference configuration q = 0, and A i = C i (0) is the reference configuration of body i.The two representations of joint screw coordinates are related by where, in vector representation of screws, the adjoined transformation Ad corresponding to C ∈ SE (3) is given by the matrix For sake of simplicity, the following notations are used so that Ad C = Ad r Ad R .The twist of body i in body-fixed representation Here v b i = R T i ṙi the body-fixed translational velocity, i.e. the velocity of the origin of the bodyfixed reference frame (RFR) F i of body i measured in the IFR F 0 and resolved in F i , whereas v s i = ṙi + r i × ω s i is the spatial translational velocity, i.e. the velocity of the point of the body that is momentarily passing through the origin of the IFR F 0 resolved in the IFR.The body-fixed and spatial angular velocity, ω b i and ω s i , is defined by ω b i = R T i Ṙi and ω s i = Ṙi R T i , respectively.The hybrid twist is defined as V h i = (ω s i , ṙi ) T , and finally the mixed twist as V m i = (ω b i , ṙi ) T .The four representations are related as follows The twist of body i within a kinematic chain is determined in terms of the generalized velocities q as with the Jacobian J b i in body-fixed, J s i in spatial, J h i in hybrid, and J m i in mixed representation.The ith column of the Jacobian is respectively given by These are the instantaneous joint screw coordinates in body-fixed, spatial, and hybrid representation.The Jacobians are related as The representations of joint screw coordinates are related by (94) and by where r i is the current position of body i in C i .The twists admit the recursive expressions Summarizing the twists of all bodies in V b , V s , V h ∈ R 6n , respectively, admits the expressions in terms of the system Jacobians that admit the factorizations [62] This provides a compact description of the overall MBS kinematics.The explicit relations for the inverse of the matrices A is the starting point for deriving recursive forward dynamics O (n) algorithms.

B Rigid Body Motions and the Lie Group SE (3)
For an introduction to screws and to the motion Lie group SE (3) the reader is referred to the text books [5,48,65,77].

B.1 Derivatives of Screws
Let C i be time dependent.According to (95), the corresponding frame transformation of screw coordinates from F i to F 0 is X ≡ 0 X = Ad Ci i X. Assume that the screw coordinates expressed in body-fixed frame are constant.The rate of change of the screw coordinates expressed in IFR is is the Lie bracket on se (3), also called the adjoint mapping.In vector notation of screws, denoting a general screw vector with with The form (114) is known as the screw product [14,77].The matrix (115) has appeared under different names, such as 'spatial cross product' in [29,30,39], or the 'north-east cross product' [13].The Lie bracket obeys the Jacobi identity Allowing for time dependent body-fixed screw coordinates i X, the above relation gives rise to an expression for the time derivative of screw coordinates in moving frames This is the spatial extension of Euler's formula for the derivative of a vector resolved in a moving frame.
For the sake of simplicity throughout the paper, the following notations are used Then the matrices Screws are the geometric objects embodying twists, wrenches, and momenta of rigid bodies.These different physical meanings imply different mathematical interpretations of the geometric object.A wrench, defined by a force and moment, is denoted with W = (t, f ) T .The force applied at a point with position vector p generates the moment t = p × f .The dual to Chasles theorem is the Poinsot theorem stating that every system of forces can be reduced to a force together with a couple with moment parallel to the force.Geometrically a screw is determined by the Plücker coordinates of the line along the screw axis and the pitch.If e is the unit vector along the screw axis, and p is a position vector of a point on that axis, the screw coordinate vector of a twist is V = (ω, v) T = ω (e, p × e + he) T , where ω = ω is its magnitude, and h = v T ω/ω 2 is its pitch.The screw coordinate vector of a wrench, i.e. the force f producing a torque t about the axis e when the point of application is displaced according to p from the axis, is W = (t, f ) T = f (p × e + he, e) T , with pitch h = t T f / f 2 .Apparently the linear and angular components of the screw coordinates are interchanged for twists and wrenches.The different definition of screw coordinate vectors allows to describe the action of a wrench on a twist as the scalar product: W T V is the power performed by the wrench acting on twist V.A twist 2 V represented in frame F 2 transforms to its representation in frame F 1 according to 1 V = Ad S1,2 2 V.The power conservation yields that a wrench represented in F 1 transforms to its representation in F 2 according to While this notation is useful for kinetostatic formulations, it is inconsistent in the sense that it treats screw coordinates differently for twists and wrenches.In screw theory, aiming on a consistent treatment of screw entities, a screw is represented by its coordinates as defined by (67) in [62] and the so-called reciprocal product of two screws is used [5,14,77].The latter is defined for X 1 = (ξ 1 , η 1 ) T and X 2 = (ξ 2 , η 2 ) T as X 1 ⊙ X 2 = ξ T 1 η 2 + η T 1 ξ 2 .Two screws are said to be reciprocal if X 1 ⊙ X 2 = 0. Obviously, if twists and wrenches are represented consistently with the same definition of screw coordinates, a reciprocal twist and wrench screws means that they perform no work.Geometrically, for zero pitch screws, this means that the screw axes intersect.In screw theory wrench screws are called co-screws to distinguish them from motion screws and to indicate that a wrench acts on a motion screw (a twist) as a linear operator that returns work or power.As twists form the Lie algebra se (3), wrenches from the dual se * (3).

C Nomenclature
-JFR for joint i at body i, joint i connects body i with its predecessor body i − 1 J i−1,i -JFR for joint i at body i − 1 i r -Coordinate representation of a vector resolved in BFR on body i.
The index is omitted if this is the IFR: r ≡ 0 r.R i -Rotation matrix from BFR F i at body i to IFR F 0 R i,j -Rotation matrix transforming coordinates resolved in BFR F j to coordinates resolved in F i r i -Position vector of origin of BFR F i at body i resolved in IFR F 0 r i,j -Position vector from origin of BFR F i to origin of BFR F j x -skew symmetric matrix associated to the vector x ∈ R 3 C i = (R i , r i ) -Absolute configuration of body i.This is denoted in matrix form with C i C i,j = C −1 i C j -Relative configuration of body j w.r.t.body i k v j -Lie algebra of SE (3) -algebra of screws q ∈ V n -Joint coordinate vector V n -Configuration space ) with the body mass m and the inertia tensor Θ c expressed in the body-fixed COM frame F c .Let S bc = R bc , b d bc ∈ SE (3) be the transformation from the COM frame F c to the BFR F b .Here b d bc is the position vector from the BFR to the COM resolved in the BFR.Then the configuration of F c is given in terms of that of the BFR by C c = CS bc .The inertia matrix w.r.t. to the general BFR F b is

1 2
tr (Θ i ) I − Θ i is used in the mass matrix M b i .Then (39) is replaced by Mb ic = diag (ϑ i , m i I), and (40) by Mb i = Ad −T S bc Mb ic Ad −1 S bc .This leads to 119) are used to denote the matrix (115) for twists ω V and v V, respectively, which are the infinitesimal versions on (96).B.2 Wrenches as Co-Screws -se * (3)

i-
Translational velocity of body i measured at origin of BFR F j , resolved in BFR F k v b i ≡ i v i i -Body-fixed representation of the translational velocity of body i v s i ≡ 0 v 0 i -Spatial representation of the translational velocity of body i k ω i -Angular velocity of body i measured and resolved in BFR F k ω bi ≡ i ω i -Body-fixed representation of the angular velocity of body i ω s i ≡ 0 ω i -Spatial representation of the angular velocity of body i k V j i -Twist of (RFR of) body i measured in F j and resolved inF k V b i ≡ i V i i -Body-fixed representation of the twist of body i V s i ≡ 0 V 0 i -Spatial representation of the twist of body i V h i ≡ 0 V i i -Hybrid form of the twist of body i V b -Vector of system twists in body-fixed representation V s -Vector of system twists in spatial representation V h -Vector of system twists in hybrid representation V m -Vector of system twists in mixed representation W b i -Applied wrench at body i in body-fixed representation W s i -Applied wrench at body i in spatial representation W h i -Applied wrench at body i in hybrid representation M b i -Inertia matrix of body i in body-fixed representation M s i -Inertia matrix of body i in spatial representation M h i -Inertia matrix of body i in hybrid representation Ad R -Screw transformation associated with C = (R, 0) Ad r -Screw transformation associated with C = (I, r) Ad Ci,j-Transformation matrix transforming screw coordinates represented in F j to screw coordinates represented in F i ad X -Screw product matrix associated with screw coordinate vector X∈ R 6 [X, Y] -Lie bracket of screw coordinate vectors X, Y ∈ R 6 .It holds [X, Y] = ad X Y. X ∈ se (3)-4 × 4 matrix associated with the screw coordinate vectors X ∈ R 6 SE (3) -Special Euclidean group in three dimensions -Lie group of rigid body motions se(3)