On the rotational equations of motion in rigid body dynamics when using Euler parameters

Many models of three-dimensional rigid body dynamics employ Euler parameters as rotational coordinates. Since the four Euler parameters are not independent, one has to consider the quaternion constraint in the equations of motion. This is usually done by the Lagrange multiplier technique. In the present paper, various forms of the rotational equations of motion will be derived, and it will be shown that they can be transformed into each other. Special attention is hereby given to the value of the Lagrange multiplier and the complexity of terms representing the inertia forces. Particular attention is also paid to the rotational generalized external force vector, which is not unique when using Euler parameters as rotational coordinates.


Introduction
Multibody systems allow the modeling of entire systems of rigid and flexible bodies connected by joints and driven by forces and actuators. The most common and widely used type of rotational coordinates in multibody systems that treat the orientation of either rigid or flexible bodies are Euler angles. An alternative nonminimal representation of the orientation of a body is given by a so-called quaternionic parametrization. In recent years, quaternionic parametrizations found new attraction because they substantially simplify the mathematical formulation and they lead, in contrast to parametrization groups using three variables, to a singularity-free description of the rotation.
The present paper deals with the different forms of the rotational equations of motion of an unconstrained rigid body in multibody systems using Euler parameters. The latter parameters are unit quaternions, that is, a collection of four real parameters, which are not independent. In terms of Euler's theorem [1], three of the four Euler parameters represent a rotation axis and the fourth is considered as single rotation about the latter axis. Since the four Euler parameters are not independent, one has to consider the quaternion constraint in the equations of motion. This is usually done by using the Lagrange multiplier technique [2] where the derivative of the unity of the Euler parameters is related to the angular velocity. A numerical time integration of the latter equations, however, would lead to errors in the unity of the quaternions. There-fore, it is essential to explicitly extend the equations by the unity of the quaternions yielding a set of DAEs instead of ODEs. In this case, the Lagrange multiplier has no physical meaning. Alternatively, one can find several approaches avoiding the Lagrange multiplier technique, e.g., in the work by Möller and Glocker [3], in which unconstrained quaternions are used. An alternative derivation of the quaternion equations of motion for a rigid body has been presented by Udwadia and Schutte [4]. In [5], the equations of motion have been derived firstly in terms of quaternions using Lagrange's formalism in order to analyze constrained mechanical systems including non-redundant holonomic constraints. Moreover, the latter paper shows as well how the physical torque vector can be expressed as the generalized quaternion torque vector. Alternatively, in [6], a different form for the generalized quaternion torque vector is presented, and by analytically solving the equations of motion, it has been demonstrated that the Lagrange multiplier according to the unity of the quaternion constraint is equal to twice the rotational kinetic energy, see also [7].
Similar to the work in [6], Vadali [8] has demonstrated that Lagrangian equations can be properly reformulated, such that the value of the Lagrange multiplier is zero. Although the Lagrange multiplier is zero within the latter formulation, it would lead to an error in the numerical time integration if the unity of the Euler parameters is no longer explicitly considered in the equations as a result of the drift phenomenon. Therefore, one has to re-normalize the Euler parameters after each numerical time integration step on position level in order to fulfill the unity of the Euler parameters, see e.g., [3,7,9].
Although one can find different forms of the generalized quaternion torque vector and as well suggestions on the amplitude of the Lagrange multiplier when using Euler parametrization, the connection between the generalized quaternion torque vector and the Lagrange multiplier has not yet been presented in the open literature. In the present paper, particular attention is paid on various different forms of the generalized quaternion torque vector and their influence on the amplitude of the Lagrange multiplier in rotational dynamics.
Furthermore, in the present paper, the generalized force vector is split into two parts, one part takes into account the rotational dynamics of the body, while the second part does not influence the rotation of the body at all. This fact has been mentioned as well in the work by Udwadia and Schutte [4], in which the connection between the physically applied torque and the generalized torque associated to the Euler parameters has been derived. To be more precise, the derivation in [4] reveals that the generalized torque vector includes a component in the direction of the Euler parameter vector which does not contribute to the rotational dynamics of a rigid body.
The present paper shows that all the forms of the equations of motion presented herein lead to the same rotational dynamics of the body. The difference between the presented forms is given by the amplitude of the Lagrange multiplier and the mathematical complexity in the equations of motion. The Euler's equations are the starting point of consideration in order to show that the Lagrange parameter concerning the Euler parameter constraint gets zero, contrary to the derivation in [8], in which the Lagrangian equations are reformulated leading to the Euler's equations. Moreover, the derivation of the equations of motion for a single rigid body from d'Alembert's variational principle in the present paper reveals various forms when using Euler parameters for the description of the rotation. Of particular importance is the presented derivation of the generalized quaternion torque vector in terms of the physically applied force vector, which can also be related to the approach presented in [6,10]. In addition to the work in the latter mentioned publications, in the present paper the different forms of the generalized quaternion torque vector are analyzed leading to new insights and their influence on the Lagrange parameter will be emphasized herein.

Kinematic description
This section defines the kinematic variables of interest, and recalls a number of fundamental kinematic relations, which are needed in the present paper. Fundamental derivations concerning Euler parameters and a collection of many useful quaternion identities are presented, see also [6].
In order to describe the position vector r p of an arbitrary point P on a rigid body in the floating frame of reference formulation [11] two different types of coordinate systems.
One is a time-and space-fixed inertial coordinate system and the second one is attached to the rigid body and, therefore, translates and rotates with the body.
Thus, the absolute position vector of a point P can be written as [11][12][13] where R is the position vector of the origin of the spacefixed coordinate system to the origin of the body coordinate system; the orthogonal rotation matrix A transforms quantities expressed in the body coordinate system into their counterparts in the space-fixed frame, and u p is the position vector of the point defined in the body coordinate system. Unless otherwise stated, throughout the present paper, quantities with an overbar are expressed in the body reference system. Note that any vector can, at any instant of time, be decomposed in either the space-fixed or the body coordinate system. While the three translational coordinates of the vector R are needed to fix the location of a body, in the present paper a four-component vector p containing Euler parameters is used for the parametrization of the orientation of the body. The four quantities e 0 , e 1 , e 2 and e 3 represent the Euler parameters and are not independent, since they fulfill the following relation i.e., In terms of the latter quaternions, the orthogonal rotation matrix A, which transforms quantities from the body coordinate system into the global coordinate system, is given by [14] A = 2e 2 where I is an order-three identity matrix, and the skewsymmetric matrix e, with the axial vector e, satisfies for an arbitrary vector v. The latter rotation matrix can also be written as [14] A = GL T with the pair of 3 × 4 matrices G and L given by [14] G = A direct calculation reveals that while where I 4 denotes the 4 × 4 identity matrix. In addition, G, L and p satisfy the following identities Due to the orthogonality of the rotation matrix A, we have The differentiation of the previous equation verifies that the matrix A TȦ is skew-symmetric, since Let us, therefore, define where the axial vector of ω may be interpreted as the angular velocity vector expressed in the body coordinate system. The skew-symmetric matrix ω can be alternatively written in terms of the matrix L and its first time derivative in the form The kinematic equation that relates the time derivatives of the Euler parameters to the components of the angular velocity vector ω can be written in the form [15] ω = 2Lṗ (17)

Equations of motion
In this section, we recapitulate the derivation of the equations of motion for a single rigid body from d'Alembert's variational principle and show that we end up with various forms of these equations if Euler parameters are used for the description of the rotation.
Assuming that the origin of the body coordinate system is attached to the center of mass of the body, d'Alembert's principle for an unconstrained rigid body states that [16] δR where m denotes the mass of the rigid body and J its tensor of inertia, referred to the body coordinate system. The sum of all forces acting on the body is denoted by f and the sum of the moments induced by f is defined by n. Any other pure moments acting on the rigid body are defined by the vector m. Without loss of generality, we assume, from now on, that no pure moments act on the body. Equation (18) must hold for any virtual displacement vector δR and for any rotation vector δθ , which may be expressed by Euler parameters in a similar way as ω in Eq. (17): Hence, by substituting Eqs. (16), (17) and (19) into d'Alembert's principle and considering only the rotational part (i.e., we set δR = 0), we obtain However, due to the fact that Euler parameters are not independent quantities, we may not equate the expression in brackets to zero. Because of Eq. (4), the variation δp is constrained by the relation Hence, the zero term δp T pλ may be added to Eq. (20), yielding The Lagrange multiplier λ may now be defined by the condition that the expression in brackets vanishes. This results in the following four equations for the four Euler parameters and the Lagrange multiplier λ: The missing fifth equation is the constraint equation, see Eq. (4). The term pλ represents the constraint force, whereby p is the transpose of the constraint Jacobian of the kinematic constraint p T p = 1. It can be easily demonstrated that the Lagrange multiplier λ in Eq. (23) is zero, see also [8]. For that purpose, we pre-multiply Eq. (23) with p T and make use of Lp = 0. This results in which can only be true if λ = 0, since p T p = 1. At this point, it should be mentioned that Eq. (23) could be obtained as well by inserting Euler parameter identities of Eqs. (16) and (17) into the well-known Newton Euler equations and adding the constraint force pλ.
Another form of the equations of motion is obtained by eliminating the Lagrange multiplier from Eq. (23) by transforming the equation back to the threedimensional space. This is achieved by pre-multiplying Eq. (23) with the 3 × 4 matrix L. Since Lp = 0 and by making use of LL T = I, we get Because we have now reduced the number of equations from four to three, the constraint equation, i.e., Eq. (4), must still be considered to determine the four Euler parameters. It should be noticed that Eq. (25) can also be derived from the original form of d'Alembert's principle. Since Eq. (18) must hold for any virtual rotation δθ , we obtain the Newton Euler equations and by substituting Eqs. (16,17) for ω and ω, we end up with Eq. (25) again.
One can obtain a third form of the equations of motion by starting again with d'Alembert's principle of Eq. (20). The latter equation is constrained by the relation δp T p = 0. Therefore, we substitute the identity of Eq. (11), i.e., L T L = I 4 − pp T , in the second term on the left side of Eq. (20), which actually represents the gyroscopic moment. This results in Since δp T p = 0, the gyroscopic moment can be simplified to 8L T JLṗ and, thus, the latter equation can be rewritten in a simpler form as As it is the case in Eq. (20), one may not equate the expression in brackets to zero because the Euler parameters are not independent quantities. First one has to add the zero term δp T pλ to Eq. (28) and then the expression in brackets can be set to zero which leads to the following form of the equations of motion: which has been presented as well in [6]. The latter equation is constrained by p T p = 1, which always has to be satisfied. However, the latter form is simpler than Eq. (23), since it contains only third-order terms in the Euler parameters, whereas we have to deal with fifthorder expressions in Eq. (23). This is mainly based on the simplification of the gyroscopic moment.
Since we still have four equations for p but no additional variable λ, it seems that the constraint equation is not needed anymore. However, this is not the case, because Eq. (32) also results from multiplying Eq. (25) with the 4 × 3 matrix L T . Thus, by linear combination, we have just produced four from three equations, which, of course, are linear dependent. Nevertheless, knowing that λ = 0 in Eq. (23) might be helpful to control the quality of a numerical solution of Eq. (23). Summing up, we found three different forms of the equations of motion for Euler parameters, given by Eqs. (23), (25) and (29). All of them must be solved in combination with the constraint equation p T p = 1. It seems that the most complicated form is Eq. (23) since it contains fifth-order terms in the Euler parameters and a Lagrange multiplier. In Eq. (25), the multiplier has been eliminated and the terms are of fourth order in the Euler parameters. Finally, Eq. (29) contains only third-order terms in the Euler parameters and a Lagrange multiplier.
There exist different strategies for solving the constrained rotational equations of motion when using Euler parameters. One approach is based on differentiating the kinematic Euler parameter constraint twice with respect to time and on using the augmented form of the equations of motion. This solution procedure is also known as index reduction and will be explained in brief in the remainder of this section. Another solution strategy is based on the elimination of the Lagrange multiplier and is explained in detail in [17]. If the kinematic Euler parameter constraint, see Eq. (4), is differentiated twice with respect to time, we obtain p Tp +ṗ Tṗ = 0 (33) With Eq. (23), we end up with the following set of equations for p and λ: The inverse of the augmented mass matrix on the left side reads [6] 4L T JL p Finally, Eqs. (34) and (35) yield If we use the relation Lp = 0, we notice from the second line of Eq. (36) that showing once again that the Lagrange multiplier λ in Eq. (23) is zero. Moreover, the first line of Eq. (36) is a set of four non-redundant equations determining the Euler parameters: where advantage has been used of LL T = I. In order to complete this section, in a last step, we solve the augmented form of the equations of motion of Eq. (29), i.e., Using the inverse of the augmented mass matrix of Eq.
(35), we obtain where the fact has been used that LL T = I and that Lp = 0. The Euler parameter accelerations of Eq. (40) are the same as in Eq. (38). One can also find the formula for the Euler parameter accelerations in [17]. However, the value of the Lagrange multiplier is not zero when using the equations of motion of Eq. (29), see Eq. (30) for the nonzero value of the Lagrange multiplier.

A detailed discussion on the generalized external force vector
In the last section, different forms of the rotational equations of motion of a rigid body have been presented by manipulating the gyroscopic moment. There, less attention has been paid to the generalized external force vector. The present section focuses on this issue in detail. It will be demonstrated that another form of the equations of motion is obtained by keeping the gyroscopic moment in its original form given by the second term on the left-hand side of Eq. (23), and by manipulating the generalized external force vector. This also will end up in a Lagrange multiplier value different from zero, see Sect. 4.3 for this purpose.
In the last section, we assumed that no pure moments m act on the body and this assumption will hold on in the present section. Thus, the only moment that acts on the rigid body results from an externally applied force vector f, i.e., where the position of application point of the force vector f is given by The generalized rotational force vector is represented by the right-hand side of Eq. (23), i.e., Substituting Eq. (41) into Eq. (43) yields The latter form of the generalized torque vector can also be found in [5]. In the following subsection, an alternative approach for computing the generalized external force vector will be presented.

Derivation of the generalized force vector based on the principle of virtual work
One can also obtain the generalized rotational force vector resulting from the physical force vector f by using the principle of virtual work. The virtual work of an externally applied force f is defined as which can be written in terms of the generalized coordinates R and p as where Q R represents the generalized force vector associated with the translational coordinates and Q p is the generalized force vector associated with the generalized rotational coordinates. For the evaluation of vectors Q R and Q p , the total differential of Eq. (42) has to be computed. Due to the fact that the body under consideration is a rigid body and, therefore,u p = 0, the first time derivative of Eq.
(42) is given bẏ Using the fact thatȦ = A ω, which follows from Eq. (15), one can write the above equation also aṡ where the fact has been used that for two arbitrary vectors y and z the following relation holds Substituting Eq. (17) into Eq. (48) leads tȯ where the 3 × 4 matrix B 1 is defined as In a similar way asṙ p , we can express the virtual displacement δr p as An alternative form of the total differential of Eq. (42) can be obtained by usinġ We may define the 3 × 4 matrix B 2 to be With the use of this equation and Eq. (53), (47) leads tȯ or in terms of virtual displacements A direct comparison of the matrices B 1 and B 2 shows that At this point, it should be noted that the matrices B 1 and B 2 are identical if three-dimensional parametrization groups, e.g., Euler angles, are used as rotational coordinates. However, when using Euler parameters, the matrices are not identical. Without any further knowledge on the latter two matrices, we obtain two different forms of the rotational generalized force vector when using the principle of virtual work. The first form is obtained by substituting Eq. (52) into Eq. (45) and by considering only the rotational part (i.e., we set δR = 0) which leads to Since B 1 = −2A u p L, one can show that which can be written under the usage of Eqs. (43) and (44) also as Alternatively, when using Eq. (56), we may write the rotational part of the virtual work of the externally applied force as such that Thus, we end up with two forms for the generalized external force vector associated with the generalized rotational coordinates. The first form is given by Eq.
(44) and is exactly the same vector as the vector Q p,1 defined in Eq. (58). The second form is given by the vector Q p,2 of Eq. (62). A direct comparison shows that The two different forms of the generalized torque vectors have been implemented for an arbitrarily chosen force vector, and they lead to the same orientation of the body in space. However, the value of the Lagrange multipliers is different for the two formulations. This issue will be discussed in detail in the next subsection.

Some investigations on the matrices B 1 and B 2
In order to be able to show that both generalized torque vectors Q p,1 and Q p,2 lead to the same orientation of the body but to different values of the Lagrange multiplier, the matrices B 1 and B 2 are carefully examined in the present section. Let us start with the matrix B 2 .
The matrix B 2 is defined as the partial derivative of the matrix product Au p with respect to the Euler parameter vector p, see Eq. (54), and can be expanded by using the definition of the rotation matrix of Eq. (5) as [14]: Since for any two vectors y and z one can verify by direct calculation that we can write Eq. (64) also as We may define the 4 × 4 skew-symmetric matrix H to be With the use of this definition and Eqs. (2,8), (66) can be expressed in compact form as Thus, the generalized external force vector associated with the rotational coordinates of Eq. (62) can be written as By taking the first time derivative of Eq. (4), one can verify that the second part of the right-hand side of the previous equation is orthogonal toṗ. It can be noted that this part does not contribute to the orientation of the body but has an influence on the Lagrange multiplier if implemented. The fact that its implementation has no influence on the orientation of the body can be demonstrated by substituting the expression of Q p,2 into Eq. (61), i.e., The last term on the right-hand side of the previous equation vanishes due to the orthogonality of p and δp. Thus, the second part of the right-hand side of Eq. (69) does not play any role in the rotational dynamics of the body. Now, we try to split the matrix B 1 of Eq. (51) into two components in order to show that the generalized force vector associated with the rotational coordinates given by Eq. (58) consists of a part that is orthogonal toṗ and therefore does not influence the rotational dynamics, and a second part that causes a rotation of the body.
or in a more compact form as Using Eqs. (7,72), the matrix B 1 of Eq. (51) can be rewritten as Substituting Eq. (11) into Eq. (73) yields where advantage has been taken of the fact that Gp = 0, see Eq. (12). Thus, we can also write the rotational generalized external force vector of Eq. (58) as The latter equation informs us that-as it is the case for Q p,2 -the generalized torque vector Q p,1 consists of a component in the direction of p and a second component which plays a role in the rotational dynamics of a rigid body. The form of Q p,2 is much simpler than the form of Q p,1 since it just depends linearly on the Euler parameters and, thus, yields to a Jacobian which does not depend on the Euler parameters. In contrast, Q p,1 depends cubically on the Euler parameters.
The reason for doing so is that Eq. (20) leads to the most common formulation of the rotational equations of motion where the Lagrange multiplier is zero and the gyroscopic moment is given in its original form. The focus in the present section lies on the incorporation of the external generalized moment vector into the equations of motion and its effect on the value of the Lagrange multiplier.
Using the definition of n of Eq. (41), we get for Eq. (20) whereby the variation δp is constrained by the relation δp T p = 0. Equation (77) can be rewritten by using Eqs. (58) and (59) as Without any further knowledge on the matrix B 1 and by adding the zero term δp T pλ to Eq. (78), we end up with the following equations of motion Clearly, the latter equation has to be constrained by the kinematic Euler parameter constraint. Since B T 1 f = 2L T n, see Eqs. (58) Since δp T p = 0, the latter equation can be rewritten in a simpler form as As it has been explained in Sect. 3, it is not possible to equate the expression in brackets to zero, because the Euler parameters are not independent quantities. Therefore, we first have to add the zero term δp T pλ to Eq. (81), which finally results in By differentiating the kinematic Euler parameter constraint twice with respect to time, the augmented form of the equations of motion can be written as Using the inverse of the augmented mass matrix of Eq.
(35) we finally get for the Lagrange parameter and the Euler parameter accelerations: Since LH T G T f = u p A T f = n, the Euler parameter accelerations in Eq. (84) are the same as in all other presented forms of the equations of motion, but the Lagrange multiplier takes in this case the value 2p T H T G T f. This is based on the manipulation of the external force vector.
Following the same procedure as demonstrated above for Q p,1 , one gets for the second form of the generalized force vector, i.e., Q p,2 , also two different forms of the equations of motion. This is based on the fact that Q p,2 can be separated into two parts, see Eq. (69), whereby the part 2p u p T f can be eliminated already at d'Alembert's principle level since δp T p is zero.
At this point it should be pointed out once again that all presented forms of the equations of motion must be solved in combination with the constraint equation p T p = 1. Moreover, all forms lead to the same Euler parameter acceleration, thus to the same movement of the body, but the value of the Lagrange multiplier is different. The Lagrange multiplier is only zero, when the most common form, i.e., Eq. (23), is implemented. For this special formulation, the Lagrange multiplier is zero because all generalized moments (inertia, gyroscopic, external) in the equations of motion are orthogonal to the Euler parameter vector p, and only the moment resulting from the kinematic Euler parameter constraint, i.e., pλ, is parallel to the Euler parameter vector. However, by using some Euler parameter identities, it is possible to eliminate some terms at the level of d'Alembert's principle (manipulating either the gyroscopic moment or the generalized external force vector), and for the resulting equations of motion, the Lagrange multiplier is no longer zero then.

Conclusion
Various forms of the rotational equations of motion of a rigid body have been presented for the case of using Euler parameters for the rotation. They must be solved in combination with the constraint equation of the Euler parameters. The crucial difference between the various types of equations of motion is given, on the one hand, by their mathematical complexity and, on the other hand, by the value of the Lagrange multiplier. In some cases, the Lagrange multiplier takes the value of zero, which might be helpful to control the quality of the numerical solution. Although it might seem paradoxical at first sight, for the case of using Euler parameters as rotational coordinates, both the gyroscopic moment and the generalized torque vector could be represented in various different forms. They all lead to the same rotation of the body but to different values of the Lagrange multiplier. Although the present paper deals with rigid body dynamics, exactly the same problem occurs if the body under consideration is flexible.