Closed-form time derivatives of the equations of motion of rigid body systems

Derivatives of equations of motion (EOM) describing the dynamics of rigid body systems are becoming increasingly relevant for the robotics community and find many applications in design and control of robotic systems. Controlling robots, and multibody systems comprising elastic components in particular, not only requires smooth trajectories but also the time derivatives of the control forces/torques, hence of the EOM. This paper presents the time derivatives of the EOM in closed form up to second-order as an alternative formulation to the existing recursive algorithms for this purpose, which provides a direct insight into the structure of the derivatives. The Lie group formulation for rigid body systems is used giving rise to very compact and easily parameterized equations.


Introduction
Rigid body dynamics algorithms for evaluating the equations of motion (EOM) and their derivatives find numerous applications in the design optimization and control of modern robotic systems. The equations of motion can be differentiated with respect to state variables, control output (generalized forces), time and physical parameters of the robot (see [8] for an overview). These derivatives can be computed with several methods: 1) approximation by finite differences, 2) automatic differentiation [31], i.e. by applying the chain rule formula in an automatic way knowing the derivatives of basic functions (cos, sin or exp), 3) closed-form derivatives of the EOM and 4) recursive and analytical formulations exploiting the structure of the closed-form equations of motion. While the first two methods are generic and numerical in nature, the latter two are analytical in nature and exploit the structure of the EOM. Analytical and recursive partial derivatives of the EOM of rigid body systems with A. Müller a.mueller@jku.at S. Kumar shivesh.kumar@dfki.de respect to state variables and generalized forces have been reported in the literature [3,17]. These are useful in optimal control of legged robots (e.g. differential dynamic programming [19]) and their computational design and optimization [12]. Time derivatives of EOM are required for the model-based control and motion planning of robots with higher-order continuity [27], since for highly dynamic applications not only the actuation forces but also their derivatives must be bounded in order to ensure feasibility. Flatness-based control of robots with flexible joints, and of robots equipped with series elastic actuators (SEA) or variable stiffness actuators (VSA), necessitate the first and second time derivatives of the EOM of the robot [6,9,26]. Therefore, recursive O (n)-algorithms for the evaluation of the time derivatives were developed [1,2,10,11,21,24] extending existing O (n)-formulations for the evaluation of EOM. While O (n)-formulations are deemed computationally advantageous when dealing with large systems, formulating and evaluating the EOM in closed form remains an efficient alternative for many robotic systems and provides insights into the structure of the problem. Yet, such closed-form formulations were not reported in the literature, with the exception of [13] where first-order time derivatives of the EOM were presented within the so-called spatial operator framework. A relatively recent research topic, where higher-order derivatives of the EOM are required, is the dynamic balancing of articulated mechanisms [4,30,32]. In [5], the time derivatives of the spatial momentum were used to derive global balancing conditions. Recently, we also proposed nth order time derivatives of EOM in both recursive and closed forms [15].
In this paper, the first and second time derivatives of the EOM are presented. For the subsequent treatment, the EOM are written in the form suitable for solving the inverse dynamics problem where the vector of generalized coordinates q = (q 1 , . . . , q n ) T comprises the n joint variables, M and C is the generalized mass and Coriolis matrix, respectively, and Q grav represents generalized gravity forces. The generalized forces due to external loads (e.g. interaction/contact forces and torques) are summarized in Q ext . Finally, Q are the generalized forces (drive forces/torques) required for a prescribed motion q (t).
In the following, the derivatives of (1) are derived using the Lie group formalism reported in [23], which is equivalent to those presented in [18] and [25]. A salient feature of the Lie group formulation is that it admits model description in terms of readily available data without compromising computational efficiency (of closed-form expressions as well as O (n) algorithms). As a side contribution, we also prove the structural properties of EOM from a closed-form perspective. For the sake of simplicity, and without loss of generality a single serial kinematic chain, comprising 1-DOF joints, mounted at the ground is considered. The generalizations to systems with arbitrary tree-topologies is straightforward, but will not be considered here in order to simplify notation and make the paper easily accessible.
Organization: Section 2 presents the equations of motion of serial kinematic chain in closed form using the body-fixed representation of the twists. Section 3 and Sect. 4 presents the first-and second-order time derivatives of the equations of motion in closed form, respectively. Section 5 proves the structural properties of the EOM using the closed-form formulations. Section 6 presents the application of the proposed derivatives in evaluating secondorder inverse dynamics of two exemplary robot manipulators and a discussion on its computational performance. Section 7 concludes the paper.

Kinematics in terms of joint screws
In the following, the notation and formulation of the EOM are adopted from [23]. The configuration (pose, posture) of body i = 1, . . . , n is denoted C i ∈ SE (3), which describes the frame transformation from a body-fixed frame (arbitrarily located at the body) to the inertial frame. Bodies and joints are numbered increasing order starting from the ground, so that joint i links body i to its predecessor i − 1, while the ground is indexed with 0, and by convention C 0 = I. The relative configuration of body j with respect to body i is then where R ij ∈ SO (3) is the rotation matrix transforming coordinates expressed in the reference frame on body j to their expression in the reference frame on body i, and i r i,j ∈ R 3 is the position vector from the origin of frame i to the origin of frame j expressed in frame i. Without loss of generality, all joints are assumed to have 1-DOF. Denote with q i the joint variable (rotation angle, translation) of joint i. The configuration of body i is determined by the product of exponential (POE) as where is the configuration of body i relative to its predecessor i − 1, in the reference configuration q = 0, and i X i is the screw coordinate vector associated to joint i represented in the body-frame of body i [22]. The vectors i X i are constant due to the body-fixed representation. For 1-DOF lower-pair joints they are given as where i e i ∈ R 3 is a unit vector along the joint axis, i x i ∈ R 3 is the vector to a point on this axis, and h i ∈ R is the pitch of the joint. In particular, for a revolute and prismatic joint, the screw coordinate vector is, respectively, In terms of these screw coordinates X, the exponential map attains the explicit form [18,20] where the rotation matrix is determined by the Euler-Rodrigues formula exp(ϕ e) = I + sin ϕ e + (1 − cos ϕ) e 2 .
Denote T the twist of body i represented in the body-fixed frame. The superscript b is used to indicate the body-fixed representation [23,25]. An alternative representation is so-called spatial representation of twist, which would be indicated by a superscript s. A recursive O (n)-algorithm for evaluating the EOM and their higher-order derivatives using the spatial representation was reported in [24]. In [14] the EOM were presented in closed form in terms of the spatial representation of twists. The potential advantage of using the spatial representation to express the EOM in closed form remains to be explored, however. In this paper the (classical) body-fixed representation is used.
The individual twists of all bodies are summarized in the vector V ∈ R 6n , which is referred to as the system twist in body-fixed representation. It is determined as with the system Jacobian J (q). The latter admits the factorization in terms of the block-triangular and block-diagonal matrices The matrix Ad C i,j transforms screw coordinates represented in the reference frame at body j to those represented in the frame on body i [18,25,28]. With the relative configuration (2) this matrix is wherex ∈ so (3) is the skew-symmetric matrix associated to vector x ∈ R 3 so thatxy = x × y. A central relation for deriving the EOM in closed form is the following expression for the time derivative of the matrix A and thus of the system Jacobian [23]: where a (q) := diag (q 1 ad1 X 1 , . . . ,q n adn Xn ).
Therein, the matrix adi X i (also called as spatial cross product by Featherstone [7]) is given in terms of the joint screw coordinate vector (4) as This gives rise to the closed-form expressions for the system acceleratioṅ For calculating the derivatives, the time derivative of matrix A will be needed. To this end, the expression is used [23], with the matrix With (16), the derivative of A is theṅ Using the relationȦd [23], the derivative of D attains the closed forṁ Clearly, the derivative (12) of the system Jacobian is recovered asJ =ȦX noting that aX ≡ 0.

Equations of motion
The generalized mass and Coriolis matrix in the EOM (1) of a simple kinematic chain mounted at the ground are found via Jourdain's principle of virtual power as (or likewise as the Lagrange equations) [23] where Therein, the (constant) 6 × 6 inertia matrix of body i expressed in the body-frame is defined as where m i is the body mass, i c i is the position vector to the COM of body i measured in the reference frame at body i, and Θ b i is the inertia tensor with respect to the body-fixed frame. The latter is related to the inertia tensor with respect to the COM, denoted Θ c i , via Steiner's ,c denotes the rotation matrix of the COM frame to the body frame.
A closed form of the EOM is obtained after replacing the system twist by (8). Alternatively, first the kinematic relation (8) and then the coefficient matrices in (22) are evaluated for a given state q,q. The generalized gravity forces are given as with Here, 0 g is the vector of gravitational acceleration expressed in the inertial frame, which is transformed to the individual bodies by U. The effect of contact or external wrenches acting on the bodies is given by the generalized forces The vector W ext,i accounts for applied load at body i. For instance, W ext,n can be used to describe the wrench at the end-effector of a robotic arm. More generally, W ext and thus Q ext , may account for arbitrary (time and velocity dependent) loads at the system.

First time derivative of the equations of motion
The first time derivative of the generalized forces iṡ The time derivative of the generalized mass matrix M in (20) follows withJ in (12) aṡ with where A = A (q) and a = a (q). By the same token, the time derivative of the generalized Coriolis matrix C isĊ where now The expression for C in (22) together with (12) yieldṡ withȧ =ȧ (q) ,ḃ =ḃ(V), and thus The form (34) allows for reusing expression (22) for the matrix C, where (22) is evaluated withq,V instead with the velocities. Also in the second form (35), terms like Aa, MAa, and b T M can be reused. A direct calculation yieldṡ where the relationU = −AaU (obtained by time differentiation of (26) and using the identity in (19)) and the fact that G is constant has been used. The time derivative of (27) along with (12) yields the expression for calculating first-order derivative of generalized forces due to external wrenchesQ

Second time derivative of the equations of motion
The second time derivative of the generalized motor forces (1) is determined bÿ The second time derivative of the generalized mass matrix is (1) . Inserting relation (19) in the differentiation of (30) yieldsṀ

Equation (43) allows for reusing M (1) (q,q). It should be observed that the term M (1) (q,q)
is the relation (30) evaluated withq instead ofq. Repeated time derivative of the generalized Coriolis matrix yields C q,q,q, ...
Taking the derivative of (36) yields the second time derivative of the generalized gravity forcesQ grav (q,q,q) = J T M (2) UG. (52) The second time derivative of the generalized forces (27) due to end-effector loads is readily found to beQ

Structural properties of the EOM
There are two important and well-known structural properties of EOM: 1. The generalized mass matrix M is symmetric and positive definite. 2.q T Ṁ (q,q) − 2C (q,q) q ≡ 0 for anyq Positive definiteness and the symmetric properties of the generalized mass matrix follows directly from its definition in (20). The symmetry property indeed applies to all higher derivatives, which is also evident from (29) and (40). An important property of the EOM (1), which is crucial for proving stability of passivity-based control schemes, is that the Coriolis matrix can be formulated asC so thatṀ − 2C is skew symmetric. To show this property, notice that bJq = bV = diag ad V 1 V 1 , . . . , ad Vn V n ≡ 0, so that Cq can be rewritten as In view of the time derivative (29) of M, the so defined Coriolis matrixC satisfies the relationC With (54), it hence followṡ It should be remarked that the Coriolis matrix is not unique, and the property 2) holds true for the particular arrangement as in (22). It may not hold if the equations are arranged differently. The skew symmetry ofṀ − 2C is indeed carried over to its derivatives.

Examples
The higher-order closed-form inverse dynamics formulation presented in Sect. 3 and Sect. 4 were implemented in MATLAB. 1 This MATLAB implementation can also be used for efficient symbolic code generation in MATLAB and C languages. This section presents the Fig. 1 Planar 2R robot schematic application of this work to the computation of higher-order inverse dynamics of two robot manipulators and presents a discussion on its computational efficiency.

Planar 2R robot
A 2R serial chain consists of a base, two links, and two revolute joints. A simple 2R chain is shown in Fig. 1. First link of length L 1 is connected to the base or ground through a revolute joint. The second link of length L 2 is connected to first link through a revolute joint. The center of mass is shown with a red circle on the links and lies at the end of each link. The mass of first and second links are m 1 and m 2 respectively. The gravity acting on the system is shown as g in the figure.
The relative reference configurations for the two bodies are The mass matrices of the two bodies in the body-fixed configuration according to (24) are

Second-order inverse dynamics
With the above information about the robot, using the closed-form expressions of M(q), C(q,q) and Q grav (q) provided in Sect. 2.2 in (1), one can arrive at the analytical expression for generalized forces (τ 1 , τ 2 ) which solves inverse dynamics problem as follows: These expressions are indeed identical to those reported in [18] for instance. Using the closed-form expressions ofṀ(q,q),Ċ(q,q,q) andQ grav (q,q) provided in Sect. 3 in (28), one obtains the analytical expression for first-order time derivatives (τ 1 ,τ 2 ) of the generalized forces, which solves the first-order inverse dynamics problem, as follows:τ The correctness can be easily checked taking the analytical time derivative of (τ 1 , τ 2 ) above.
Using the closed-form expressions ofM(q,q,q),C(q,q,q, ... q ) andQ grav (q,q,q) provided in Sect. 4 in (38), one obtains the analytical expression for second-order time derivatives of generalized forces (τ 1 ,τ 2 ), solving the second-order inverse dynamics problem: They can again be verified by taking analytic derivatives of the above (τ 1 ,τ 2 ).

KUKA LBR iiwa manipulator
The closed-form formulation presented in this paper are used to compute the second-order inverse dynamics of the 7 degrees of freedom KUKA LBR iiwa robot.
The relative reference configurations B i in (3)    Dynamic model parameters The mass and inertia data used in this example is taken from the identified model of KUKA LBR iiwa robot reported in [29] as listed in Table 1 Second-order inverse dynamics The joint trajectory is described by means of cos-function as shown in Fig. 3a, and used for the inverse dynamics computation of the robotic manipulator (see Fig. 3b). Figure 4a and Fig. 4b show first-and second-order time derivatives of the generalized forces. The results computed from the closed-form expressions are compared against the numerical differentiation. As apparent from Fig. 4, the closed-form derivatives match the numerical time differentiation of generalized forces in both cases which attests the correctness of the formulation presented in this paper.

Computational performance
The computational performance of the closed-form algorithm was evaluated by measuring the total CPU time spent in 10000 evaluations of secondorder inverse dynamics on a standard laptop with Intel Core i9-7960X CPU clocked at 2.80 GHz and 128 GB RAM. It was found that the closed-form implementation takes 16.34 seconds leading to the average computational time of 1.6 milliseconds per evaluation, which is sufficient for many real-time control applications.

Conclusion and outlook
This paper presents closed-form expressions for the first and second time derivative of the EOM of a kinematic chain. Building upon the Lie group formulation of the EOM, the formulations are advantageous as they are expressed in terms of joint screw coordinates, and thus facilitate parameterization in terms of vector quantities that can be easily obtained. The computationally efficiency of these closed-form relations compared to recursive algorithms and their efficient implementation will be a topic of further research. It is already obvious from the presented expressions that they involve many repeated terms that can be precomputed and reused. Future research will also address the time derivatives of general mechanisms with kinematic loops and an efficient C++ based implementation in Hybrid Robot Dynamics (HyRoDyn) software framework [16]. Fig. 4 Results for the higher-order inverse dynamics regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.