Efficient Analytical Derivatives of Rigid-Body Dynamics using Spatial Vector Algebra

An essential need for many model-based robot control algorithms is the ability to quickly and accurately compute partial derivatives of the equations of motion. State of the art approaches to this problem often use analytical methods based on the chain rule applied to existing dynamics algorithms. Although these methods are an improvement over finite differences in terms of accuracy, they are not always the most efficient. In this paper, we contribute new closed-form expressions for the first-order partial derivatives of inverse dynamics, leading to a recursive algorithm. The algorithm is benchmarked against chain-rule approaches in Fortran and against an existing algorithm from the Pinocchio library in C++. Tests consider computing the partial derivatives of inverse and forward dynamics for robots ranging from kinematic chains to humanoids and quadrupeds. Compared to the previous open-source Pinocchio implementation, our new analytical results uncover a key computational restructuring that enables efficiency gains. Speedups of up to 1.4x are reported for calculating the partial derivatives of inverse dynamics for the 50-dof Talos humanoid.


I. INTRODUCTION
Rigid-body dynamics models are widely used in the development of control algorithms for quadruped and humanoid robots, with recursive dynamics algorithms at the core of many real-time controllers. Example use cases include current approaches to controller design via optimization [1], [2]. Likewise, robotics libraries such as Crocoddyl [3] and Drake [4] require the partial derivatives of rigid-body dynamics with respect to the state and control variables. The calculation of these derivatives represents the majority of the CPU time required in many optimal control applications [5].
There are multiple general-purpose methods that can be applied to rigid-body dynamics for computing derivatives. Finite differences [7]- [9] are simple to implement and are trivially parallelized, but suffer in accuracy ( Table I). The complex-step method [12] is accurate to machine precision but requires all functions to be computed in the complex plane, leading to additional overhead. Cossette et al. [13] extended this method to matrix Lie groups, providing direct applicability to rigidbody dynamics. An approach that avoids the overhead of the complex step is automatic differentiation (AD), also called algorithmic differentiation [4], [10], [11]. AD calculates the partials of an algorithm by automating the accumulation of chain-rule expressions through operator overloading [10] or source code transformation [4], [11]. The first approach overloads basic data types to compute partials concurrently with all arithmetic. The second approach builds an expression graph from source code and augments the source directly to calculate the partial derivatives.
Without relying on AD, but similar to it, others have developed analytical methods that accumulate chain-rule expressions on top of existing algorithms. Examples include many existing descriptions of recursive algorithms for dynamics derivatives [5], [15], [16]. The derivation in Ref. [5] applies chain rule on the classical two-pass Recursive-Newton-Euler Algorithm (RNEA) for inverse dynamics (ID). The presented method has computational complexity O(N d) for both the forward and backward pass of the RNEA derivatives, where N is the number of bodies in the system and d is the depth of the kinematic connectivity tree. However, the algorithm in the C++ Pinocchio (2.6.0) code [6] associated with [5] is distinctly different from the algorithm presented in [5]. The Pinocchio code includes a more efficient O(N ) forward pass, coupled with an O(N d) backward pass.
Derivative algorithms that directly carry out chain-rule on an existing algorithm (using AD or by hand) may not always be the most efficient. Other special-purpose methods have been proposed to fully exploit the structure of the rigid-body equations of motion. Such approaches often analytically differentiate closed-form equations of motion and then design an algorithm to compute the result. An example is Ref. [17] where partial derivatives are considered of the mass matrix M (q), the Coriolis matrix C(q,q), and the gravity vector g(q) with respect to the state variables q andq, leading to an O(N 3 ) algorithm.
Ref. [18] gives a recursive method for partial derivatives of inverse and forward dynamics (FD) for serial kinematic chains with single-DoF joints. Their approach includes an O(N 2 ) algorithm for the partials of FD, and an O(N d) algorithm for the partials of ID.
The main contribution of this paper is to extend previous analytical results on partial derivatives of ID [18] to generic multi-DoF joint robots with a fixed or  1 Algo-2,3 of [5] provides a chain-rule method for ID partial derivatives. 2 Pinocchio code for ID partial derivatives accompanying [5] but different from chain-rule. 3 This paper provides a Special-Purpose Analytical Method that outperforms the state-of-the-art in terms of speed. floating base. We first derive closed-form expressions for the first-order partial derivatives of ID. To the best of authors' knowledge, these results represent the first of their kind for general rigid-body systems (fixed or floating base) with multi-DoF joints. The new expressions lead immediately to an algorithm of order O(N d) that naturally generalizes the one shown in [18]. The method is fundamentally different than the straight chain-rule approach presented in [5]. However, the distinct algorithm in the Pinocchio code [6] ends up being directly related to the computations required in our algorithm. Despite the similarity, our new closed-form expressions uncover a key restructuring that accelerates computations. We use the relationship between FD and ID [5], [18] in an efficient manner to ultimately calculate the partial derivatives of FD with complexity O(N 2 ).

II. DERIVATIVES OF RIGID-BODY DYNAMICS
Rigid-Body Dynamics: For a rigid-body system, the state variables are the configuration q and the generalized velocity vectorq, while the control variable is the generalized torque vector τ . The equation of motion, also called the Inverse Dynamics (ID), is given by where M ∈ R n×n is the mass matrix, C ∈ R n×n is a Coriolis matrix, and g ∈ R n is the vector of generalized gravitational forces, and n is the DoF of the system. For a fixed state, ID calculates τ for a givenq (Eq. 1), and Forward Dynamics (FD) computesq for a given τ : Toward supporting control applications (e.g., [1], [9]), the objective of this work is to compute the partial derivatives of FD with respect to the state (q,q) and control variables (τ ). The most efficient algorithms for calculating ID and FD are the O(N ) Recursive-Newton-Euler Algorithm (RNEA) and the Articulated-Body Algorithm (ABA) respectively [19].
Derivatives of Inverse and Forward Dynamics: As presented in [5], a direct approach to compute the partial derivatives of ID is to manually differentiate the RNEA algorithm via chain-rule, denoted as RNEACR (Table II). A similar chain rule approach for the derivatives of ABA is denoted as ABACR. Because the ABA is more computationally intensive than RNEA, a more efficient approach for derivatives of FD first applies RNEACR, and then uses the relationship [5], [18] where the variable u ∈ {q,q}. From Eq. 3, ∂F D/∂τ can be directly calculated as M −1 , enabling re-use in Eq. 5. This method for the partials of FD via RNEA Chain Rule is abbreviated as FDCR.
The state-of-the-art numerical method for computing partial derivatives of FD is provided in the Pinocchio package [6] (Pinocchio FD Derivs, Table II). The method first computes the derivatives of ID (Pinocchio ID Derivs, Table II). It then computes M −1 and uses Eq. 5 to obtain the derivatives of FD.
In the following sections, Spatial Vector Algebra (SVA) is reviewed for dynamics analysis, followed by a theoretical development of analytical expressions that lead to an algorithm for partial derivatives of ID. The new expressions result in algorithm with a key difference from the state-of-the-art [6], enabling a speedup over it.

III. SPATIAL VECTOR ALGEBRA (SVA)
Notation: Spatial vectors are 6D vectors that combine the linear and angular aspects of a rigid-body motion or net force [19]. Spatial vectors are denoted with lowercase bold letters (e.g., a), while matrices are denoted with capitalized bold letters (e.g., A). Motion vectors, such as velocity and acceleration, belong to a 6D vector space denoted M 6 . Spatial vectors are usually expressed in either the ground coordinate frame or a body coordinate (local coordinate) frame. For example, the spatial velocity k v k ∈ M 6 of a body k expressed in the body frame is given by is the angular velocity expressed in a coordinate frame  fixed to the body, while k v k ∈ R 3 is the linear velocity of the origin of the body frame. When expressed in the ground frame, the spatial velocity for body k is denoted as 0 v k . The vector is again composed of angular and linear velocity components. However, the linear velocity is associated with the body-fixed point on body k that is coincident with the origin of the ground frame 0. In this paper, when the frame used to express a spatial vector is omitted, the ground frame is assumed. Force-like vectors, such as force and momentum, belong to another 6D vector space F 6 , which is dual to M 6 . A spatial cross-product between two motion vectors (v,u), written as (v×)u, is given by (Eq. 6). This operation can be understood as providing the time rate of change of u, when u is moving with a spatial velocity v. For a Cartesian vector, the matrix ω× is the classical 3D cross product operator. A spatial cross-product between a motion and a force vector is written as (v× * )f , as defined in Eq. 7.
An operator × * is defined by swapping the order of the cross product, such that (f × * )v = (v× * )f [20]. Further introduction to SVA is provided in Appendix A.
Connectivity: An open-chain kinematic tree with serial or branched connectivity ( Fig. 1) is considered with N links connected by joints, each with up to 6 DoF. Body i's parent toward the root of the kinematic tree is denoted as λ(i). ν(i) denotes the set of bodies in the subtree rooted at body i, while ν(i) denotes the set of bodies in ν(i) excluding the body i . We define i j if body i is in the path from body j to the base. The spatial velocities of the neighbouring bodies in the tree are related by v i = v λ(i) + S iqi , where S i is the joint motion subspace matrix for joint i [19] anḋ q i the joint rates for joint i. For a single DoF revolute joint, S i ∈ R 6 , andq i is a scalar. For a 6-DoF freemotion joint, S i ∈ R 6×6 , whileq i ∈ R 6 . The velocity v i can also be written as the sum of joint velocities over predecessors as v i = l i S lql . The derivative of joint motion subspace matrix due to the axis changing with respect to local coordinates (commonly denotedS i [19]) is assumed to be zero. The quantityṠ i = v i × S i signifies the rate of change of S i due to the local coordinate system moving.

Dynamics:
The spatial equation of motion [19] is given for body k as: where f k is the net spatial force on body k, I k is its spatial inertia [19], and a k is its spatial acceleration. Instead of treating gravity as an external force, a common trick is to accelerate the base upwards opposite of the gravitational acceleration (a 0 = −a g ), providing the acceleration of body k as: For later use, we decompose a k into terms from joint accelerations, and terms from joint rates, according to In a similar fashion, the net spatial force on body k is then decomposed as: where η k = I k γ k (η ∈ F 6 ) is the spatial force on the body caused by joint accelerationsq, and ζ k = v k × * I k v k + I k ξ k (ζ∈ F 6 ) gives the Coriolis and centripetal forces on body k. From the formulation of the RNEA [19], Proceeding to consider the system overall, for any particular joint i, Eq. 1 can be written as follows.
Then, using τ i = S T i f C i , Eq. 9 and 11: where I C i denotes the composite rigid-body inertia of the sub-tree rooted at body i, given by IV. ANALYTICAL PARTIAL DERIVATIVES USING SVA In this section, closed-form expressions are derived for the derivatives of ID. The reader less interested in the derivation may skip to Eq. 30 and Eq. 32 as a summary.

A. Building Blocks
Several kinematic identities are given in Appendix C as a basis for the main derivation. We denote by n j the number of DoFs for joint j. The motion subspace matrix for the joint is then given as S j = s j,1 · · · s j,nj where each spatial vector s j,p gives the p-th free-mode of motion for joint j. With a slight liberty of notation, ∂/∂q j,p denotes an operator for the directional derivative along this p-th free mode of a joint. For example, considering the case with j i gives identity J1: which provides the rate of change in S i with respect to relative motion s j,p at joint j earlier in the chain. For revolute joints, these derivatives are just conventional derivatives with respect to joint angles. For multi-DoF joints such as a floating base, they are more formally Lie derivatives. Considering a configuration-dependent vector u, we denote by ∂u ∂q j = ∂u /∂q j,1 · · · ∂u /∂q j,n j the matrix of derivatives associated with joint j. To illustrate, we derive identity J7 for later use in the section.
Considering j i and using the definition of γ i : Using J1 and switching the order of the cross product: Collecting all DoFs of joint j, Eq. 80 becomes: Using the definition of γ l , and summing over l: The subsequent derivations employ these formulae and the identities from Appendix C to obtain the partials of the terms in Eq. 11. Since the derivations are quite lengthy, the presentation focuses on explaining the methodology. A full derivation is provided in Appendix D.
Partial Derivative of Eq. 12: We derive the partial derivative of [M (q)q] i with respect to q j for all i, j ∈ {1, . . . , N }. First, consider the case when j i. Using the product rule of differentiation in Eq. 12: Using identity J3 for the term , and cancelling terms: Upon summing over the index k, the final expression is: For the case i≺j, we have that: Using the identities J3, J4 and J7 and cancelling terms: Summing over the index k results in: where η C j collects joint-acceleration-dependent forces for the subtree, calculated as η C j = k j η k . For ease of implementation, we consider a single case where j i. The indices i and j are switched in Eq. 24, to get the expression for for the case j ≺ i. Therefore, Eq. 25 gives the two expressions formulated for the general case j i.
Partial Derivative of Eq. 13: Similarly, we use identities J2-J6 (Appendix C) to get the first-order partial derivatives of [C(q,q)q] i with respect to q j for the case j i as shown in Eq. 26, 27. In these equations, the composite of ζ C i for the subtree is defined as ζ C i = k i ζ k , and the matrix B k (v k , I k ) is a body-level Coriolis matrix [20]- [22] given by: Partial Derivative of the Gravity Term: Using the identities J3 and J4, for the case j i, the partial derivative of g i (Eq. 14) with respect to q j is: Summary for Partials of ID w.r.t. q: The partials of the individual components are now collected together. Terms in Eqs. 25, 26, 27, and 29 are added to get the total expressions for ∂τ ∂q for the case j i. The full derivation is provided in Appendix E.
The spatial quantitiesΨ j andΨ j are defined as:Ψ The partials of τ with respect toq depend only on the the Coriolis terms Cq. Using the identities J8 and J9 leads to expressions for the case when j i(see Appendix F for full derivation):

D. Algorithm for First-Order Partials of ID
Algorithm 1 returns the terms in Eq. 30 and Eq. 32 and has computational complexity O(N d). The method operates with all spatial vectors expressed in ground frame coordinates. The first pass is in the forward direction and goes from root to leaves calculating the spatial velocity, acceleration,Ṡ i ,Ψ i ,Ψ i , B i and the spatial force f i . When joint i has a single DoF,Ṡ i =Ψ i andS i =Ψ i and the steps are the same as in [18,Alg. 2].
The second pass progresses in the backward direction from leaves to the root. Index i takes all the values from 1 to N . The quantity ∂τ ∂q [ν(i), i] denotes the partial derivatives of τ for all bodies in the sub-tree of i with respect to q i , while ∂τ ∂q [i, ν(i)] is the partial derivative of τ for body i with respect to each of the bodies in the sub-tree of i, excluding body i. This backward pass again generalizes the one in [18,Alg. 2]. The differences are that the third term on line 15 can be dropped in the single-DoF case (since S T i (f C i × * )S i = 0 in that case), while lines 17-20 restructure a second inner loop in [18] into matrix multiplies.

A. Algorithm correctness
The partial derivatives of FD from ABACR, FDSVA, & FDCR are compared with derivatives calculated in a Fortran implementation using the complex-step method for accuracy. For N = 100, the ABACR method results in a term-by-term root-mean-square (rms) relative error of 10 −13 , while both FDSVA and FDCR result in an rms error of approximately 10 −12 . The rms relative error for all methods grows linearly with DoF on a log-log scale.
B. Runtime for Partials of ID/FD vs. RNEACR/FDCR as presented in [5] via Fortran implementation of both We consider an N link serial or branched kinematic tree with all revolute joints about their local z-axis.
RNEACR in body-coordinates [5, Algos. 2 & 3] (Table II) is used to calculate the partial derivatives of ID, and compared with the IDSVA (Table II) method. All the algorithms are written in Fortran 90 and implemented using the Intel Fortran compiler on a 3.07 GHz Intel Xeon processor. To calculate the average run time, each algorithm is run 10,000 times with randomized inputs for the state and control variables. Fig. 2 shows the comparison of the two methods. For N = 100, a speedup of 15× for IDSVA over RNEACR is found. Fig. 2 also shows the comparison of FDSVA with FDCR and ABACR (Table II)  Since SVA allows for coordinate-free expressions, recursive algorithms can be formulated in either bodycoordinates or ground-coordinates [19]. For the bodycoordinate algorithms, all the intermediate quantities are transformed between body coordinates using transformation matrices i X λ(i) , while for the ground-coordinates algorithms, the quantities can be left in the ground frame. A major advantage for the latter comes by avoiding the repeated transformation of the quantities between local body frames in the backward pass of the algorithm. This gain is achieved at the cost of expressing kinematic quantities (velocities 0 v i , accelerations 0 a i ), joint motion subspace matrices 0 S i , and inertia matrices 0 I i in the ground coordinate frame during the outward pass. Fig. 3 shows a comparison of IDSVA (see Table II) runtime in body coordinates and ground coordinates. Speedups for ground coordinate algorithms are between 1.3 to 1.9 for N = 2 to N = 500.

C. Comparing Runtime for C++ Partials of ID/FD with
Pinocchio Implementation Accompanying [5] IDSVA in ground coordinates is implemented in C++ within the Pinochhio framework. This strategy enables direct comparison with Pinocchio's original ID partial derivatives [6]. Fig. 4 shows the comparison of the two methods for serial and branched kinematic trees with some branching factor bf [19]. For a serial chain with N = 100, a speedup of 2× is found using the gcc-9.0  (Table II)

D. Comparing CPU Runtime for AZA vs. DMM in C++
with Pinocchio Implementation Accompanying [5] AZA is implemented within the Pinocchio [6] framework. Pinocchio's M −1 algorithm [25] is modified to include the AZA. In Fig. 7, the "crossover" N denotes the point below which the DMM performs better than the AZA. This point depends on the hardware and compiler optimization settings used to implement the algorithm, but for high N , the O(N 2 ) AZA is efficient because it avoids the expensive matrix-matrix product. Table III   gives the cross-over N values (for kinematic trees) above which AZA performs better than DMM for different compiler settings. The effect of this cross-over N can be seen in Fig. 6   calculate the inverse and forward dynamics partials for robots with multi-DoF joints and a floating base. Several algorithmic optimizations and Spatial Vector Algebra (SVA) identities are exploited to enable an efficient implementation. The method provides a 1.4× speedup for the TALOS humanoid model over the state-of-theart Pinocchio FD derivatives using the gcc compiler and a 1.2× speedup using the Clang-10 compiler. The reduction in runtime for partial derivatives enables faster optimization algorithms for both on-line and off-line applications. These improved timings can ultimately lead to better motion planning of legged and industrial robots.

A. Spatial Vector Algebra
A body k with spatial velocity k v k ∈ M 6 in the body frame is decomposed in its angular and linear components as: where k ω k ∈ R 3 is the angular velocity of the body in a coordinate frame fixed to the body, while k v k ∈ R 3 is the linear velocity of the body-fixed point at the origin of the body frame. Spatial vectors can also be expressed in the ground frame. For example, the spatial velocity of the body k in the ground frame is denoted as 0 v k . In this case, the linear velocity is associated with the body-fixed point on body k that is coincident with the origin of the ground frame. The net spatial force 0 f k ∈ F 6 defined in Eq. 36 on the body can be calculated from the spatial equation of motion (Eq.37): where 0 n k ∈ R 3 is the net moment on the body about the origin of the ground frame, 0 f k ∈ R 3 is the net linear force on body, 0 I k is the spatial inertia of the body k that maps motion vectors to force vectors, and 0 a k ∈ M 6 is the spatial acceleration of the body. The transformation matrix i X j is used to transform vectors in frame j to frame i is defined as: where i R j ∈ R 3×3 is the rotation matrix from frame j to frame i, p i/j ∈ R 3 is the Cartesian vector from origin of frame j to i, and 0 is the 3 × 3 zero matrix. p× is the 3D vector cross product on the elements of p, defined as: The spatial transformation matrix 0 X k can be used to obtain spatial velocity vector 0 v k from the vector k v k as: A spatial cross-product operator between two motion vectors (v,u), written as (v×)u, is given by (Eq. 41). This operation can be understood as providing the time rate of change of u, when u is moving with a spatial velocity v. A spatial cross-product between a motion and a force vector is written as (v× * )f , and defined in Eq. 42.
Hence, the three spatial vector cross-product operators defined above map between the motion and the force vector space as: [19] × :

B. Properties of Spatial Vectors
Assuming u, v, m, v 1 , v 2 ∈ M 6 , and f ∈ F 6 , many spatial vector properties [19] are utilized herein:

C. Multi-DoF joint Identities and Expressions
Identities are shown below for partial derivatives of common spatial quantities by perturbing a multi DoF joint position variable (q j ). These identities are used to then derive the partial derivatives of inverse dynamics for a rigid-body system with multi-DoF joints and a floating base.
otherwise Identities listed above are derived in detail as follows: J1. For directional derivative of S i along the p th freemode of the joint j, total derivative with respect to time is taken for the numerator and denominator as: Using the definition ofṠ i asṠ i = v i × S i , and the definition of v i : Only a single term remains for the case j i: J2. The directional derivative of the spatial velocity of a body i along the p th free-mode of the joint j, where j i is given as: Using J1, we get: By switching signs, Eq. 49 can also be written as: Eq 50 can be written for each p th mode of the joint j to get the partial derivative (of size 6 × n j , where n j is the number of DoF for joint j) of v i with respect to q j as: Using the definition of v l , Eq. 51 now is: Using Eq. 52, the partial derivative of v i × a with respect to q j for j i, where a ∈ M 6 is any fixed motion vector is given as: Similar to Eq. 53, the partial derivative of v i × * f , where f ∈ F 6 is any fixed force vector, is calculated as: J3. For any fixed force vector f , the partial derivative of S T i f with respect to q j for j i is calculated. Using J1: Eq. 56 can be written for each DoF of the joint j. Hence, the partial derivative with respect to q j is: J4. For any fixed motion vector a, the partial derivative of I i a is calculated with respect to q j for j i. The directional derivative of I i in the direction of p th free-mode of joint j [19] is: Multiplying a on both sides result in: Similar to Eq 56, Eq 59 can be written for each DoF of joint j collectively as: J5. The partial derivative of I i v i with respect to q j for j i is now calculated using the identities derived above. Using the product rule: For the first term in Eq. 61 on the RHS, v i is assumed to be a fixed motion vector. Hence, the identity J4 can be used. For the second term in Eq. 61, Eq. 52 is used as: Expanding terms: Upon cancellations and simplification, Eq. 63 becomes: whereΨ j is defined in Eq. 31. J6. Using identities derived above, and the definition of ξ i , the partial derivative of ξ i with respect to q j for j i is calculated as: where v J l is the joint velocity defined as: Using the product rule of differentiation: The partial derivative of the joint velocity v Ji with respect to q j for j i is calculated. Using the definition of v Ji (Eq. 68) we get: Using identity J1, we get ∂v Ji ∂q j,p = s j,p × S iqi (70) Collectively, for all DoF of joint j, Eq. 71 can be written as: Treating v J l as a fixed motion vector in the first term in Eq. 67 on the RHS, Eq. 52 is used. Eq. 53 is used for the second term.
Expanding terms, and using the definition ofΨ j (Eq. 31): Combining terms and simplifying, Using spatial vector cross property P1, definition of ξ i and simplifying, 76) Summing over the terms and simplifying, we get: J7. Using the definition of γ i , the directional derivative of γ i along the p th free-mode of the joint j for j i is calculated as: Using J1, we get: Eq. 79 can also be written as Collectively for all DoF of joint j , Eq. 80 is written as: Using the definition of γ l , Eq. 81 becomes: J8. The partial derivative of v i is calculated with respect toq j for j i. Using the definition of v i , we get: Summing over all indices l, all the terms vanish except the ones pertaining to the index j resulting in: J9. Using the definition of ξ i , the partial derivatives of ξ i with respect toq j for j i is: Using the product rule of differentiation, we get: Using J8, summing over the index l results in: Upon simplifying and using the definition ofΨ j (Eq. 31), we get: D. Partial Derivatives of ID w.r.t q: Derivations Partial Derivatives of [M (q)q] i : 1) Case when j i Using product rule of differentiation in Eq. 12, we get: Using identities J3, J4, and J7, we get: Upon cancellations and simplifications, the final expression is: 2) Case when j ≻ i For this case, identities J3, J4 and J7 are used as: Expanding: Upon cancellations and simplification, we get the following expression.
Partial Derivative of [C(q,q)q] i : 1) Case when j i Using product rule of differentiation in Eq. 13, we get: Using the property P3 to combine terms, Summing over the index k, we get: Partial Derivative of the Gravity Term:

1) Case when j i
For the case j i, the partial derivative of g i with respect to q j is: Using identities J3 and J4, we get: Cancelling terms, and summing over index k, we get: 2) Case when j ≻ i For the case j ≻ i, we follow the similar process and use identity J4 as: Expanding: Summing over the index k, we get: ∂g i ∂q j = S T i (I C j a 0 )× * S j + I C j (a 0 × S j ) (113) Re-arranging terms here UsingΨ k (Eq. 31), the final compact expression for ∂τ i ∂q j is: Re-arranging terms, we get: Using the definition of a i and f i (Eq. 10), we get: Using the definition ofΨ i (Eq. 31), we get the final expression for ∂τ j ∂q i as: F. Partial Derivatives of ID w.r.tq: Derivations 1) Case when j i Using product rule of differentiation in Eq. 13, we get: Using identities J8 and J9, we get: Simplifying and collecting terms, we get: Summing over the index k, we get: 2) Case when j ≻ i Similar to the previous case, product rule of differentiation in Eq. 13 is used: Using identities J8 and J9, we get: Simplifying and collecting terms, we get: Summing over the index k, we get: Flipping the indices i and j in Eq. 130 to get for the case j ≺ i: Hence, the partial derivatives of τ with respect toq are: