Inverse, forward and other dynamic computations computationally optimized with sparse matrix factorizations

TL;DR

This paper introduces an optimized algorithm for computing the dynamics of articulated rigid bodies using sparse matrix factorizations, reducing computational complexity compared to traditional methods.

Contribution

The authors develop a novel off-line optimized sparse matrix factorization approach to efficiently compute inverse and forward dynamics of articulated bodies.

Findings

Reduces floating point operations compared to recursive Newton-Euler and articulated body algorithms.

Performs well on classical dynamic problems with lower computational complexity.

Validated through tests showing numerical efficiency even without existing gold standards.

Abstract

We propose an algorithm to compute the dynamics of articulated rigid-bodies with different sensor distributions. Prior to the on-line computations, the proposed algorithm performs an off-line optimisation step to simplify the computational complexity of the underlying solution. This optimisation step consists in formulating the dynamic computations as a system of linear equations. The computational complexity of computing the associated solution is reduced by performing a permuted LU-factorisation with off-line optimised permutations. We apply our algorithm to solve classical dynamic problems: inverse and forward dynamics. The computational complexity of the proposed solution is compared to `gold standard' algorithms: recursive Newton-Euler and articulated body algorithm. It is shown that our algorithm reduces the number of floating point operations with respect to previous approaches. We also evaluate the numerical complexity of our algorithm by performing tests on dynamic computations for which no gold standard is available.

Figures

Equations

$${\mathbf{d}}_p=\left[\begin{array}{c}{d}_{a_1}\\ ⋮\\ d_{a_m}\end{array}\right],{\mathbf{D}}_{p,q}=\left[\begin{array}{ccc}{D}_{a_1,b_1}& \dots & D_{a_1,b_n}\\ ⋮& \ddots & ⋮\\ D_{a_m,b_1}& \dots & D_{a_m,b_n}\end{array}\right].$$

$$v\times =\left[\begin{array}{c}\omega \\ \dot{p}\end{array}\right]\times =\left[\begin{array}{cc}\omega \times & 0\\ \dot{p}\times & \omega \times \end{array}\right],v{\times }^{\ast }=\left[\begin{array}{c}\omega \\ \dot{p}\end{array}\right]{\times }^{\ast }=\left[\begin{array}{cc}\omega \times & \dot{p}\times \\ 0& \omega \times \end{array}\right]$$

$$I_i=\left[\begin{array}{cc}{I}_{C,i}+m_i{c}_i\times c_i^{\top }& m_i{c}_i\times \\ m_i{c}_i{\times }^{\top }& m_i{I}_{3\times 3}\end{array}\right],$$

$$v_i=\text{prescript}i{X}_{{\lambda }_i}{v}_{{\lambda }_i}+S_i\dot{q_i},$$

$$a_i$$

$${\tau }_i$$

$$f_i$$

$$b_{c_{\text{eqref}eq:ai}}$$

$$b_{c_{\text{eqref}eq:fi}}$$

$$D_{c_{\text{eqref}eq:ai},a_i}$$

$$D_{c_{\text{eqref}eq:ai},\ddot{q_i}}$$

$$D_{c_{\text{eqref}eq:fi},f_i}$$

$$D_{c_{\text{eqref}eq:fi},f_i^x}$$

$$D_{c_{\text{eqref}eq:taui},{\tau }_i}$$

$${\mathbf{D}}_{p,q}(q,\dot{q}){\mathbf{d}}_q+{\mathbf{b}}_p(q,\dot{q})=0,$$

$$\displaystyle\tau=\mbox{InvD}({\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}model,q,\dot{q},}\ddot{q},f_{1}^{x},\dots,f_{N_{B}}^{x}).$$

$${\ddot{q}}_i$$

$$f_1^x$$

$$D_{c_{\text{eqref}eq:yd2qi},{\ddot{q}}_i}$$

$$D_{c_{\text{eqref}eq:ydfxi},f_i^x}$$

$${\mathbf{D}}_{p_{id},q}(q){\mathbf{d}}_q+{\mathbf{b}}_{p_{id}}(q,\dot{q})=0,$$

$$\begin{gathered} q=[f_1^x,{\ddot{q}}_1,\dots ,f_{N_B}^x,{\ddot{q}}_{N_B},a_1,\dots ,a_{N_B}, \\ {f}_{N_B},{\tau }_{N_B},\dots ,f_1,{\tau }_1]. \end{gathered}$$

$$\begin{gathered} p_{id}=[c_{(\text{ref}eq:meas{f}_1^x)},c_{(\text{ref}eq:meas{\ddot{q}}_1)}\dots c_{(\text{ref}eq:meas{f}_{N_B}^x)},c_{(\text{ref}eq:meas{\ddot{q}}_{N_B})},c_{(\text{ref}eq:NEdyn{a}_1)}\dots c_{(\text{ref}eq:NEdyn{a}_{N_B})}, \\ {c}_{(\text{ref}eq:NEdyn{f}_{N_B})},c_{(\text{ref}eq:NEdyn{\tau }_{N_B})},\dots ,c_{(\text{ref}eq:NEdyn{f}_1)},c_{(\text{ref}eq:NEdyn{\tau }_1)}]. \end{gathered}$$

$$D_{c_{\text{eqref}eq:ydfxi},f_i^x}$$

$$D_{c_{\text{eqref}eq:fi},f_i}$$

$$D_{c_{\text{eqref}eq:ai},a_{{\lambda }_i}}$$

$$D_{c_{\text{eqref}eq:ai},\ddot{q_i}}$$

$$D_{c_{\text{eqref}eq:fi},a_i}$$

$$D_{c_{\text{eqref}eq:fi},f_i^x}$$

$$D_{c_{\text{eqref}eq:fi},f_j}$$

$$D_{c_{\text{eqref}eq:taui},f_i}$$

$$\displaystyle\ddot{q}=\mbox{FwdD}({\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}model,q,\dot{q}},\tau,f_{1}^{x},\dots,f_{N_{B}}^{x}).$$

$$\displaystyle\begin{split}p_{i}^{A}={}&\nu_{i}-f_{i}^{x}+\sum_{j\in\mu_{i}}\prescript{i}{}{X_{j}^{*}}\Big{\{}p_{j}^{A}+I_{j}^{a}c_{j}+\\ &+I_{j}^{A}S_{j}\left(S_{j}^{\top}I_{j}^{A}S_{j}\right)^{-1}\left(\tau_{j}-S_{j}^{\top}p_{j}^{A}\right)\Big{\}},\end{split}$$

$$I_j^a=$$

$$I_i^A=$$

$${\ddot{q}}_i$$

$$a_i$$

$${\tau }_i$$

$$D_{c_{\text{eqref}eq:ytaui},{\tau }_i}$$

Full text

I Previous works

Real-time control of complex robots (such as humanoids) asks for efficient ways of computing position, velocity, acceleration and applied wrenches on all the bodies composing the robot articulated chain. Within this paper, we will refer to these quantities with the term dynamic variables even though they include positions and velocities which are strictly speaking kinematic. Efficiently computing the dynamic variables of a (possibly free-floating) robot is nowadays a solved theoretical problem. Several state of the art algorithms [1] have a computational complexity which scales linearly in the number of rigid links composing the robot. All these efficient algorithms can be obtained by exploiting the problem sparsity which derives from the iterative propagation of forces, torques and accelerations across articulated structures. Classical problems have been extensively optimized from the computational complexity point of view. As discussed in [1, Table 10.1] the recursive Newton-Euler algorithm as implemented in [2, Algorithm 5.7] is the most computationally efficient solution of the inverse dynamics problem. Similarly, the articulated body algorithm as implemented in [3, method 3] is the most computationally efficient solution of the forward dynamics problem. Yet another relevant problem is the hybrid dynamics solution, which can be efficiently solved with the articulated-body hybrid dynamics [1, Section 9.2].

In this paper we consider the computationally efficient solution of dynamic problems which do not fall within the scope of the inverse, forward and hybrid dynamics even though they are of practical interest. The paper is organized as follows. Section II present the notation. Section III presents a specific formulation of the Newton-Euler equations and Section III-A discuss the matrix form of these constraints. Section IV presents the inverse dynamics problem and its solution with the recursive Newton-Euler algorithm (Section IV-A). Section V presents the forward dynamics problem and its solution with the articulated body algorithm (Section V-A). An alternative but computationally comparable solution of the inverse and forward dynamics is proposed in Section IV-D and Section V-D respectively. Finally, Section VI extends this solution to dynamic problems which do not fall within the scope of the inverse, forward and hybrid dynamics.

II Notation

Let $\mathcal{A}$ and $\mathcal{B}$ be two arbitrary sets with cardinality $|\mathcal{A}|=m$ and $|\mathcal{B}|=n$. For any element in $a\in\mathcal{A}$ let’s associate a vector $d_{a}$. Similarly, for any pair of elements $a,b$ with $a\in\mathcal{A}$ and $b\in\mathcal{B}$, let’s associate a unique matrix $D_{a,b}$. Let $p=(a_{1},\dots,a_{m})\in S_{\mathcal{A}}$ and $q=(b_{1},\dots,b_{n})\in S_{\mathcal{B}}$ be two permutations of the elements in $\mathcal{A}$ and $\mathcal{B}$ respectively111Equivalently, let $p\in S_{\mathcal{A}}$ and $q\in S_{\mathcal{B}}$, i.e. let $p$ and $q$ be two elements of the symmetric group on $\mathcal{A}$ and $\mathcal{B}$ respectively.. We define the permuted vector $\bm{d}_{p}$ and the permuted matrix $\bm{D}_{p,q}$ induced by $p$ and $q$ as follows:

[TABLE]

As for dynamic quantities, in this paper we follow the same notation used in [1] but we avoid using the bold symbols for matrices, since the bold symbols are reserved for distinguishing $D_{a,b}$, $a\in\mathcal{A}$, $b\in\mathcal{B}$ from $\bm{D}_{p,q}$, $p\in S_{\mathcal{A}}$, $q\in S_{\mathcal{B}}$; ${v}$ denotes the spatial velocity (a six dimensional vector including angular velocities in the first three components and the rest as linear velocities), ${a}$ denotes the spatial acceleration (again angular and then linear), ${f}$ denotes the spatial force (couples in the first three components and forces in the remaining), ${}^{B}{X}_{A}$ is the motion vector transformation from $A$ to $B$ coordinates and ${}^{B}{X}^{*}_{A}$ is the analogous transformation for a force vector. Given a vector $v$ and its coordinates ${}^{A}v$ in $A$ we denote with ${}^{A}\dot{v}$ its temporal derivative in the $A$ coordinates. The Euclidean cross operator $\times$ (on $\mathbb{R}^{3}$) is defined as the usual vector product, while for a spatial velocity $v$ the cross spatial cross operator and its dual $\times^{*}$ are defined as follows:

[TABLE]

In the following sections, we restrict the analysis to robots described by kinematic trees (i.e. robots with no kinematic loops) composed of $N_{B}$ rigid links numbered from [math] to $N_{B}$, [math] being the selected fixed base (world reference frame) and $1$ being the selected floating base (reference rigid body in the articulated chain). Links numbering is obtained by defining a suitable spanning tree where each rigid link is associated to a unique node in the tree. Numbers can be always selected in a topological order so that each node $i$ has a higher number than its unique parent $\lambda_{i}$ and a smaller number than all the nodes in the set of its children $\mu_{i}$. Links $i$ and $\lambda_{i}$ are coupled with the joint $i$ whose joint motion constraints are modelled with $S_{i}\in\mathbb{R}^{6\times 1}$ (therefore without loss of generality we are assuming that all joints have a single degree of freedom). For each rigid link $i$, the system is modelled also by supplying the spatial inertia tensor222The spatial inertia tensor of the rigid link $i$ has the following form in the link reference frame:

$I_{i}=\begin{bmatrix}I_{C,i}+m_{i}c_{i}\times c_{i}^{\top}&m_{i}c_{i}\times\\ m_{i}c_{i}\times^{\top}&m_{i}I_{3\times 3}\end{bmatrix},$

where $I_{C,i}$ is the spatial inertia tensor with respect to the body’s centre of mass, $m_{i}$ is the total mass, $c_{i}$ is the relative displacement between the centre of mass and the origin of the link reference frame. $I_{i}$ and the motion-vector transform from the reference frame of the rigid link $i$ to the reference frame of the rigid link $j$, denoted ${}^{j}X_{i}$. For each link $i$ the considered kinematic variables are:

• $v_{i}$: the link spatial velocity,
• $q_{i}$: the joint $i$ position,
• $\dot{q}_{i}$: the joint $i$ velocity,

Similarly, the dynamic variables associated to the link $i$ are:

• $a_{i}$: the link spatial accelerations,
• $\ddot{q}_{i}$: the joint $i$ acceleration,
• ${\tau}_{i}$: the joint $i$ torque,
• $f_{i}$: the spatial force transmitted to body $i$ from $\lambda_{i}$,
• $f_{i}^{x}$: external forces acting on body $i$.

Remark 1

All these variables are expressed in body $i$ coordinates including $f_{i}^{x}$ which is more often expressed in absolute (i.e. body [math]) coordinates.

III Dynamic constraints

Let’s assume that all kinematic quantities, i.e. those depending on $q$ and $\dot{q}$ have been precomputed. In practice, these quantities are the transformation matrices ${}^{j}X_{i}$, ${}^{j}X^{*}_{i}$ and the linear/angular velocities $v_{Ji}$, $v_{i}$. The latter can be efficiently computed with the following recursion, propagated from $i=1$ to $N_{B}$:

[TABLE]

The dynamics quantities instead ($\tau$, $\ddot{q}$, $f_{1}^{x},\dots,f_{N_{B}}^{x}$, $a_{1},\dots,a_{N_{B}}$, $f_{1},\dots,f_{N_{B}}$) have to satisfy the Newton-Euler equations:

[TABLE]

where we defined $\nu_{i}=v_{i}\times^{*}I_{i}v_{i}$ and $c_{i}=v_{i}\times S_{i}\dot{q}_{i}$.

III-A Dynamic constraints in matrix form

Let’s define the set of dynamic variables $\mathcal{D}=\{a_{i}$, $f_{i}$, $\tau_{i}$, $f_{i}^{x}$, $\ddot{q_{i}}$ $\}_{i=1}^{N_{B}}$. We have $|\mathcal{D}|=5N_{B}$. The elements in $\mathcal{D}$ can be thought as generic quantities, not necessarily expressed in a specific reference frame. As described in Section II, for any $a\in\mathcal{D}$, let’s associate a vector $d_{a}$ (e.g. $d_{a_{i}}$ might correspond to the coordinates of $a_{i}$ expressed in a specific reference frame). Given a permutation $q\in\mathcal{S}_{\mathcal{D}}$, $d_{q}$ contains an ordered sequence of the elements in $\mathcal{D}$ expressed in specific reference frame(s).

Equations (3) can be seen as a set of equations which the vector of dynamic variables $d_{p}$ have to satisfy. Let’s put these constraints in matrix form by defining the set of constraints $\mathcal{C}=\{c_{\eqref{eq:ai}}$, $c_{\eqref{eq:fi}}$, $c_{\eqref{eq:taui}}\}_{i=1}^{N_{B}}$. We have $|\mathcal{C}|=3N_{B}$. Given $c\in\mathcal{C}$ let’s define vectors $b_{c}$ as follows:

[TABLE]

and $b_{c}=0$ otherwise. Given $c\in\mathcal{C}$ and $d\in\mathcal{D}$, let’s define matrices $D_{c,d}$ as follows (where $j\in\mu_{i}$):

[TABLE]

and $D_{c,d}=0$ otherwise.

Given two permutations $p$ and $q$ of the elements in $\mathcal{C}$ and $\mathcal{D}$, the permuted matrix $\bm{D}_{p,q}$ and permuted vectors $\bm{d}_{q}$, $\bm{b}_{p}$ satisfy the following linear equation:

[TABLE]

where we explicitly indicated the dependency on $q$ and $\dot{q}$.

IV Inverse dynamics

The inverse dynamics problem consists in finding ${\tau_{1}}$, $\dots$, ${\tau_{N_{B}}}$ which satisfies (3) given $\ddot{q}_{1}$, $\dots$, $\ddot{q}_{N_{B}}$, $f_{1}^{x}$, $\dots$, $f_{N_{B}}^{x}$. In [1] the problem is formulated as the computation of the following function:

[TABLE]

In the above equations we grayed out some variables that will not play a role in the following sections, and can be assumed either to be contestant ($model$) or measured ($q$, $\dot{q}$).

IV-A Inverse dynamics solved with RNEA

An efficient solution of inverse dynamics is given by the recursive Newton-Euler algorithm (RNEA), described hereafter. Equations (2) and (3$a_{i}$) are propagated from $1$ to $N_{B}$ with initial conditions $v_{0}=0$ and $a_{0}=-a_{g}$ which corresponds to the gravitational spatial acceleration vector expressed in the body frame [math] (null in its first three components and equal to the gravitational acceleration in the last three). Equations (3$f_{i}$) and (3$\tau_{i}$) are propagated from $N_{B}$ to $1$.

IV-B Inverse dynamics solved with matrix inversion

In this Section, we propose a way to solve the inverse dynamics problem with a matrix inversion. Remarkably, (6) represents the set of linear constraints in $\bm{d}_{q}$. Additionally, certain components of $d_{p}$ are known since $\ddot{q}_{1}$, $\dots$, $\ddot{q}_{N_{B}}$, $f_{1}^{x}$, $\dots$, $f_{N_{B}}^{x}$ are given:

[TABLE]

These constraints, will extend the set of constraints $\mathcal{C}$. In particular, we should define $\mathcal{C}_{id}=\mathcal{C}\cup\{c_{\eqref{eq:yd2qi}},c_{\eqref{eq:ydfxi}}\}_{i=1}^{N_{B}}$ so that $|\mathcal{C}_{id}|=5N_{B}$. To extend the definition of $D_{c,d}$ and $b_{c}$ we define:

[TABLE]

and $D_{c,d}=0$, $b_{c}=0$ otherwise. Again, given two permutations $p_{id}$ and $q$ of the elements in $\mathcal{C}_{id}$ and $\mathcal{D}$, the permuted matrix $\bm{D}_{p_{id},q}$ and permuted vectors $\bm{d}_{q}$, $\bm{b}_{p_{id}}$ satisfy the following linear equation:

[TABLE]

which is representation of the inverse dynamics problem. In particular, inverting the matrix $\bm{D}_{p_{id},q}$ we can compute the solution $\bm{d}_{q}$ of the inverse dynamics. A more computationally efficient solution can be obtained as described in the following section.

IV-C Inverse dynamics solved with forward substitution

In this section we prove that with suitable permutations we can compute $d_{q}$ which solves (10) with a computationally efficient algorithm, which is the forward substitution as defined in [4]. The idea is to build two find permutations $p_{id}$ and $q$ of the elements in $\mathcal{C}_{id}$ and $\mathcal{D}$ which lead to a lower triangular matrix $D_{p_{id},q}$ and then to solve (10) with a forward substitution. The permutations are inspired by the RNEA described in Section IV-A. Let’s first consider the $q$ permutation of the elements in $\mathcal{D}$. We choose:

[TABLE]

Let’s also choose a permutation $p_{id}$ of the elements in $\mathcal{C}_{id}$ :

[TABLE]

Property 1

The permuted matrix $\bm{D}_{p_{id},q}$ induced by the permutations defined in (11) and (12) is lower triangular.

Proof:

Since $\bm{D}_{p_{id},q}$ is defined by blocks, we structure this proof by considering the blocks that constitute the matrix itself. First, we prove that the blocks in the main diagonal are diagonal matrices. These blocks are:

[TABLE]

which are indeed diagonal. We are left with proving that the blocks in the upper triangular part of $\bm{D}_{p_{id},q}$ are identically null. We consider the non-null blocks in (5) and prove that each block $D_{c,d}$ is positioned in the lower triangular part of $\bm{D}_{p_{id},q}$. This is equivalent to prove that if $D_{c,d}\neq 0$ then $i_{c}\geq i_{d}$, being $i_{c}$ and $i_{d}$ the position of $c$ and $d$ in the permutations $p_{id}$ and $q$, respectively. In the definitions given in (5), we can neglect the diagonal blocks, which have been previously considered. We are left with:

[TABLE]

Easy computations can show that in these case $i_{c}\geq i_{d}$ considering a numbering scheme with $\lambda_{i}<i$ and $j>i$, $\forall j\in\mu_{i}$ as described in [1]. ∎

IV-D Inverse dynamics solved with LU factorization

In Section IV-C, suitable permutations led to a lower triangular structure for the matrix $\bm{D}$ in (9). After obtaining the lower triangular structure, the underlying linear system can be solved with a forward substitution which is a computationally efficient algorithm to solve a linear system. In this section we are interested in computing these permutations by solving the following problem.

Problem 1

Given arbitrary permutations $p$ and $q$ of the elements in $\mathcal{C}_{id}$ and $\mathcal{D}$, compute permutation matrices $P$ and $Q$ to obtain a triangular matrix $P\bm{D}_{p,q}Q$.

In consideration of Property 1, Problem 1 has necessarily a solution. In other terms, $\bm{D}_{p,q}$ is triangularizable with permutations.

Definition 1** (triangularizable with permutations)**

A square matrix $A$ is triangularizable with permutations if there exists a permutation of the rows and a permutation of the columns which result in a triangular matrix.

To solve Problem 1, we resort to a classical problem in matrix analysis: the sparse $LU$ factorization with minimum filling-in [4, Section 11.1.9].

Problem 2

Given a square matrix $A$ find the permutation matrices $P$ and $Q$, a lower triangular matrix $L$ and an upper triangular matrix $U$ such that $PAQ=LU$ and the number of filling-in (new nonzeros in $L$ and $U$ that are not present in $A$) is minimum.

Remark 2

If $A$ is triangularizable with permutations, i.e. $A=\bar{P}\bar{L}\bar{Q}$ with $L$ triangular, then the solution to problem 2 is given by $P=\bar{P}^{-1}$, $Q=\bar{Q}^{-1}$, $L=\bar{L}$ and $U$ equal to the identity matrix. In this case in fact, we can obtain zero filling-in which by definition is the minimum number achievable. Under this considerations, we can apply Problem 2 to solve Problem 1

Remark 3

Deciding if a matrix is triangularizable according to definition 1 is $\mathcal{NP}$-complete [5]. Similarly, solving problem 2 is also $\mathcal{NP}$-complete [6]333Even though [6] is cited several times [7, 8] as a proof of the $\mathcal{NP}$-completeness of Problem 2, it has to be observed that it is not clear to the authors of the present paper how to extend the results in [6] to the case of non-symmetric positive definite matrices. Available numerical tools for finding their solution are not guaranteed to reach the minimum. In solving Problem 2 we will use the unsymmetric-pattern multiFrontal method, as implemented in UMFPACK [9].

Going back to Problem 1, we can take advantage of Property 1 to guarantee that $\bm{D}_{p,q}$ is always triangularizable. Therefore, remark 2 applied to $\bm{D}_{p,q}$ guarantees that solving Problem 2 will give us also a solution to Problem 1. Given the $\mathcal{NP}$-completeness of the underlying problem, we are not guaranteed to find a solution but numerical experiments444Experiments are avaialble here https://github.com/iron76/bnt_time_varying/tree/master/experiments/computationalComplexity/RNEA. conducted so far with UMFPACK shows that a solution is always found if $N_{B}\leq 100$. As to this concern, one-hundred can be considered a a practical upper-bound for robotic applications.

Remark 4

Even though Problem 2 is $\mathcal{NP}$-complete, in practice the problem is solved once in a preliminary optimization phase and its benefits can be exploited in the runtime computations. The preliminary optimization consists in solving Problem 2 for a worst-case sparsity pattern. The runtime computations instead consist in solving (10) for different positions $q$ and velocities $\dot{q}$ exploiting the optimized permutations computed previously.

A solution of Problem 2 can be used to compute permutations that allow to solve (10) efficiently. The solution uses the sparsity pattern of $\bm{D}_{p,q}(q)$, i.e. the pattern of non-zero elements in the matrix. However, in practical applications we are interested in a sparsity pattern which somehow represents the sparsity of $\bm{D}_{p,q}(q)$ for all possible values of $q$. As to this concern, we define the worst case sparsity pattern.

The basic observation is that the robot structure (e.g. number of degrees of freedom, joint types, joint positions, tree structure of the robot) does not change. As a consequence, the underlying sparsity structure of (10), i.e. the non-zero elements in $\bm{D}$, changes only for the state-dependent elements, i.e. those that depend on $q$ and $\dot{q}$. Looking at (5), the only state-dependent blocks are the transformations $\prescript{i}{}{X_{j}^{*}}$, $\prescript{i}{}{X_{\lambda_{i}}}$ which depend on $q$. These sub-matrices do have a state-dependent sparsity structure but for the purpose of this paper, we can consider the associated worst-case sparsity structure, i.e. if an element is non-zero for at least one value of $q$, then it is considered as a non-zero element in the associated sparsity pattern. For classical joint types (revolute, prismatic, helical, cylindrical, planar, spherical and free-motion) the worst-case sparsity pattern can be easily computed by observing that the only $q$ dependent elements are either sines of cosines (see [1, Table 4.1, page 79]). These functions are zero only on a countable number of configurations (and therefore on a subset whose Lebesgue measure is zero) which are easy to enumerate.

V Forward dynamics

The inverse dynamics problem consists in finding $\ddot{q}_{1}$, $\dots$, $\ddot{q}_{N_{B}}$ which satisfies (3) given ${\tau_{1}}$, $\dots$, ${\tau_{N_{B}}}$, $f_{1}^{x}$, $\dots$, $f_{N_{B}}^{x}$. In [1] the problem is formulated as the computation of the following function:

[TABLE]

Again, in the above equations we grayed out some variables that will not play a role in the following sections, and can be assumed either to be contestant ($model$) or measured ($q$, $\dot{q}$).

V-A Forward dynamics solved with the ABA

The articulated-body algorithm [1, ABA] solves the forward dynamics problem in $O(N_{B})$ computational complexity. The algorithm consists in the following steps. First, the articulated body bias forces $p_{i}^{A}$ and the articulated body inertias $I_{i}^{A}$ are recursively computed iterating with $i$ = $N_{B}$, $\dots$, $1$ the following equations:

[TABLE]

Then the following two equations are iterated with $i$ = $1$, $\dots$, $N_{B}$ and initial condition $a_{0}=-a_{g}$.

[TABLE]

V-B Forward dynamics solved with matrix inversion

In solving the inverse dynamics, the fact that $\bm{D}_{p_{id},q}$ is lower triangular follows from the specific structure of the measurement equations (8). Changing the measurement equations would compromise the lower triangularity of $\bm{D}_{p_{id},q}$. As a consequence, the associated linear system would not have the suitable structure to apply the forward substitution. Within this context, forward dynamics give a useful example. The measured variables for the forward dynamic case are $f^{x}_{i}$ and ${\tau}_{i}$. The latter in matrix notation can be expressed as follows:

[TABLE]

These constraints, will extend the set of constraints $\mathcal{C}$. In particular, we should define $\mathcal{C}_{fd}=\mathcal{C}\cup\{c_{\eqref{eq:ytaui}},c_{\eqref{eq:ydfxi}}\}_{i=1}^{N_{B}}$ and extend the definition of $D_{c,d}$ and $b_{c}$ as follows:

[TABLE]

and $D_{c,d}=0$, $b_{c}=0$ otherwise. Given two permutations $p_{fd}$ and $q$ of the elements in $\mathcal{C}_{fd}$ and $\mathcal{D}$, a solution of the forward dynamics can be computed as the unique solution $\bm{d}_{q}$ of the following liner system:

[TABLE]

V-C Forward dynamics solved with forward substitution

Similarly to what observed in the inverse dynamic case, we might try to find permutations for $p_{fd}$ of the elements in $\mathcal{C}_{fd}$ and $q$ of the elements in $\mathcal{D}$ to obtain a matrix $\bm{D}_{p_{fd},q}$ somehow simple to invert. Specifically, we might think that the articulated-body algorithm [1, the ABA] presented in Section (V-A), could be translated into suitable permutations that transform $\bm{D}_{p_{fd},q}$ into a lower triangular matrix. Unfortunately, the ABA algorithm is instead something more than a permutation as discussed in the following.

Property 2

There exist:

a permutation $p_{fd}$ of the elements in $\mathcal{C}_{fd}$;

• -

a permutation $q$ of the elements in $\mathcal{D}$;

• -

matrices $W^{R}_{q_{1},q_{2}}$ defined for $q_{1},q_{2}\in\mathcal{D}$;

• -

matrices $W^{L}_{p_{1},p_{2}}$ defined for $p_{1},p_{2}\in\mathcal{C}_{fd}$;

such that:

[TABLE]

is lower triangular. These quantities lead to the following matrix reformulation of the ABA algorithm:

[TABLE]

The proof of the preposition above is constructive and can be found in Appendix -B.

V-D Forward dynamics solved with LU factorization

Again, part of the computational optimizations of the ABA proposed in Section V-A and revisited in Section V-C consists in row and column permutations. In this section we will suggest an algorithm for simplifying the forward dynamics computations by leveraging the intrinsic sparsity of the underlying matrices (16) and pre-computing suitable row and column permutations that will reduce the underlying computational cost. The idea consists in computing the worse case sparsity pattern for $\bm{D}_{p_{fd},q}(q)$ in (16). Solving Problem 2 with this sparsity pattern will give the row-column permutations $P$ and $Q$ suitable for solving (16). At run-time regardless of the specific $q$ the idea is to perform the sparse $LU$ factorization on $P\bm{D}_{p_{fd},q}(q)Q$ followed by a backward substitution on $U$ and a forward substitution on $L$. Figure 1 shows the computational cost of this solution against the algorithms proposed in Section V-A (ABA) and Section V-C (matrix reformulation of the ABA). The computational cost of ABA is the one reported in [1, Page 202] and corresponds to the algorithmic solution proposed in [10]. The other computational costs are numerically computed with software which is available open source555The software for these computations are available here: https://github.com/iron76/bnt_time_varying/tree/master/experiments/computationalComplexity/ABA/serial/compare..

VI Dynamic equations and LU-factorization

In the previous sections, we have seen how the RNEA and the ABA are computationally efficient solutions of the inverse and forward dynamics problems. The reduction of the computational costs is obtained by suitable permutations of the matrices $\{D_{c,d}\}_{c\in\mathcal{C}_{id},d\in\mathcal{D}}$ and $\{D_{c,d}\}_{c\in\mathcal{C}_{fd},d\in\mathcal{D}}$. In this section we consider the problem of finding similar permutations for a wider class of problems. Inspired by the matrix representations (10) and (16) of the inverse and forward dynamics respectively, we consider a generic estimation problem consisting in computing the solution $\bm{d}$ of the following linear system:

[TABLE]

obtained by combining (6) with a generic measurement equation $\bm{D}_{Y}(q)\bm{d}+\bm{b}_{Y}(q,\dot{q})=0$. In Section IV-C and Section V-C respectively, suitable permutations led to a lower triangular structure for the matrix $\bm{D}$. The underlying linear system could then be solved with a forward substitution. However, the proposed solutions rely on the specific structure of the problem and therefore extending them to solve (17) is non trivial. The algorithm proposed in Section IV-D and Section V-D seems more suitable to simplify the computational complexity of solving (17) in its generic form (i.e. beyond the inverse and forward dynamics cases).

As an example, we consider a relevant case study: the inverse dynamic computation of a humanoid robot standing on two feet. The specificity of this problem is that the robot has two external wrenches applied at the feet. Differently from the inverse dynamic computations presented in Section IV-A, we assume that these external wrenches are unknown. Using the notation used in previous sections, these external wrenches can be denoted $f_{1}^{x}$ and $f_{N_{B}}^{x}$ by labelling with $1$ and $N_{B}$ the right and the left foot respectively. We therefore consider the following problem:

[TABLE]

where we explicitly indicated that $f_{1}^{x}$ and $f_{N_{B}}^{x}$ are unknown even if they act on the system. The considered mechanical system is a free-floating articulated rigid body subject to constraints. Unfortunately solving (18) is ill-posed, i.e. given $q,\dot{q},\ddot{q},f_{2}^{x},\dots,f_{N_{B}-1}^{x}$ there exist multiple $\tau$, $f_{1}^{x}$, $f_{N_{B}}^{x}$ satisfying (3). To estimate $\tau$ we could try using additional measurements. We hereafter consider quite a common sensor distribution, nominally four load cells located on each foot666We first prove that these sensors are not yet sufficient to solve the inverse dynamics. Then we propose a different set of sensors which results in a well posed inverse dynamics problem.. These sensors are available in the NAO robot [11], in the QRIO robot [12] and the Atlas [13]. These sensors correspond to a three-axes force-torque sensor, i.e. they measure the net force orthogonal to the plane of the each foot and the projection of the torque on the same plane. As it is often the case, if we choose the left and right feet reference frames with the $x$ and $y$ axis aligned to the foot plane we have:

[TABLE]

with $y_{lf},y_{rf}\in\mathbb{R}^{3}$ being the measured three-axes force-torque and $h_{lf},h_{rf}\in\mathbb{R}^{3}$ being the non-directly measured components of the contact wrench $f_{1}^{x}$. We are now considering the following problem:

[TABLE]

To understand if this problem is well-posed, it is convenient to resort to a specific form of equation (3). This formulation is detailed in [14, eq. (41)] and it combines the floating-base dynamics with the joint dynamics. We consider here a simplified version obtained by choosing the left foot as the base frame and choosing $q$ a local parametrization of the robot pose (composed by the free-floating six-dimensional configuration and the joint $n$-dimensional configuration). We have:

[TABLE]

Using (19) in (21) and grouping the unobservable/observable quantities as in problem (20), we obtain:

[TABLE]

Understanding if the inverse dynamic problem with these measurements is well-posed boils down to understanding if the following matrix is invertible:

[TABLE]

Given the upper triangular form of this matrix, we are left with proving the invertibility of $\left[\begin{smallmatrix}H,&\prescript{0}{}{X_{N_{B}}^{*}}H\end{smallmatrix}\right]$. We are interested in guaranteeing the invertibility of this matrix for any $\prescript{0}{}{X_{N_{B}}^{*}}$, i.e. regardless of the feet relative pose.

Property 3

Let $H=\left[\begin{smallmatrix}H_{f}^{\top}&H^{\top}_{\mu}\end{smallmatrix}\right]^{\top}$ with $H_{f}\in\mathbb{R}^{3\times 3}$, $H_{\mu}\in\mathbb{R}^{3\times 3}$. If $H_{f}$ is singular, then there always exists $\prescript{0}{}{X_{N_{B}}^{*}}$ which makes $\left[\begin{smallmatrix}H,&\prescript{0}{}{X_{N_{B}}^{*}}H\end{smallmatrix}\right]$ singular.

Proof:

Using the structure of $\prescript{0}{}{X_{N_{B}}^{*}}$ and $H$ we obtain:

[TABLE]

and the results follows by observing that the multiplicand matrix on the left is singular (the first three lines are linearly dependent) and by using the fact that the rank of a product is always less or equal the rank of the multiplied matrices. ∎

Remark 5

Using Property 3 with definitions (19b) we can conclude that (20) is an ill-posed problem. In other terms, the inverse dynamics problem cannot be solved in the case of a humanoid robot standing on the two feet with four load cells on each foot.

Additional assumptions are needed to compute the inverse dynamics. A first realistic assumption in certain applications is to assume that the robot is standing in a very slippery surface and therefore the tangential forces along the $x$ and $y$ axis are negligible. In this case:

[TABLE]

with $y_{lf},y_{rf}\in\mathbb{R}^{3}$ being the measured contact forces and the measured torques on the $x$-$y$ plane and $h_{lf},h_{rf}\in\mathbb{R}$ being the torques on the $z$ axis. In this case it was numerically observed777https://github.com/iron76/bnt_time_varying/tree/master/experiments/computationalComplexity/SIE that the associated inverse dynamics problem is solvable. Fig. 2 shows a comparison of the number of floating point operations necessary to solve this specific inverse dynamics problem with or without a sparse $LU$ factorization. Remarkably this case cannot be solved with classical algorithms (e.g. hybrid dynamics [1]).

-A The articulated body equation of motion

The idea consists in recursively computing ($i$ = $N_{B}$, $\dots$, $1$) the quantities $p_{i}^{A}$ (articulated body bias forces) and $I_{i}^{A}$ (articulated body inertias) which satisfy the articulated body equation of motion:

[TABLE]

For defining these quantities, let’s start by considering $i$ such that $\mu_{i}=\emptyset$. In this case (3$f_{i}$) gives $f_{i}=I_{i}a_{i}+\nu_{i}-f_{i}^{x}$ and therefore $I_{i}^{A}=I_{i}$ and $p_{i}^{A}=\nu_{i}-f_{i}^{x}$. Recursively, let’s assume that $I_{j}^{A}$ and $p_{j}^{A}$ have been defined for every $j\in\mu_{i}$ (this time non-empty) and let’s find suitable expressions for $I_{i}^{A}$ and $p_{i}^{A}$. This is achieved by a five step procedure.

Step 1 - $S_{1}$. Replace (24$f_{i}^{ab}$) in (3):

[TABLE]

The latter is almost the equation we need for the recursive definition of $I_{i}^{A}$ and $p_{i}^{A}$ if only we could write $a_{j}$ as a function of $a_{i}$. This is achieved with the following steps.

Step 2 - $S_{2}$. Substitute (25$a_{i}$) in (25$\tau_{i}$):

[TABLE]

Step 3 - $S_{3}$. Multiply the previous equation by the inverse of $S^{\top}_{j}I_{j}^{A}S_{j}$ to obtain (14$\ddot{q}_{i}$) for $\ddot{q}_{j}$:

[TABLE]

Step 4 - $S_{4}$. Substitute the last equation in (25$a_{i}$):

[TABLE]

which evaluated with the substitution $i\rightarrow j\in\mu_{i}$ leads to:

[TABLE]

Step 5 - $S_{5}$. Substitute $a_{j}$ in (25$f_{i}$) with its expression in (27):

[TABLE]

and enforcing $f_{i}=I_{i}^{A}a_{i}+p_{i}^{A}$ leads to the definitions in (13$p_{i}$), (13$I_{i}^{a}$) and (13$I_{i}^{A}$).

-B Proof of Property 2

We hereafter assume that the articulated body inertias $I_{i}^{A}$ ($i$ = $1$, $\dots$, $N_{B}$) have been computed. The idea is to follow the steps presented in Section -A to compute the sub-blocks of the matrices $\bm{W}^{L}$ and $\bm{W}^{R}$. The latter is defined as a matrix equivalent of Step 1. The former as $\bm{W}^{L}=W^{L,4}W^{L,3}W^{L,2}W^{L,1}$ with $W^{L,1}$ representing Step 2, $W^{L,2}$ representing Step 3, $W^{L,3}$ representing Step 4 and $W^{L,4}$ representing Step 5.

Step 1 - $S_{1}$. We compute a matrix $W^{R}$ which multiplied by $d_{q}$ replaces $f_{i}$ with $p_{i}^{A}=f_{i}-I_{i}^{A}a_{i}$. As usual we define $W^{R}$ by its blocks $W^{R}_{q_{1},q_{2}}$ with $q_{1},q_{2}\in\mathcal{D}$. We have:

[TABLE]

and $W^{R}_{q_{1},q_{2}}=0$ otherwise.

Remark 6

Given two permutations $p$ and $q$ of the elements in $\mathcal{C}$ and $\mathcal{D}$, if $D_{p,q}$ represents (3), $D_{p,q}W^{R}_{q,q}$ represents (25). Therefore, in the following $\{c_{\eqref{eq:ai}}$, $c_{\eqref{eq:taui}}$, $c_{\eqref{eq:fi}},\}_{i=1}^{N_{B}}$ can be read as $\{c_{\eqref{eq:ABA_ai}}$, $c_{\eqref{eq:AB_taui}}$, $c_{\eqref{eq:AB_with_aj}}\}_{i=1}^{N_{B}}$.

Step 2 - $S_{2}$. We define $W^{L,1}$ which left multiplies $D_{p,q}W^{R}_{q,q}$ to substitute (25$a_{i}$) in (25$\tau_{i}$). The exceptions to $W^{L,1}_{p_{1},p_{2}}=0$ are:

[TABLE]

Step 3 - $S_{3}$. We define $W^{L,2}$ which left multiplies $W^{L,1}_{p,p}D_{p,q}W^{R}_{q,q}$ to multiply (26) by the inverse of $S^{\top}_{j}I_{j}^{A}S_{j}$. The exceptions to $W^{L,2}_{p_{1},p_{2}}=0$ are:

[TABLE]

Step 4 - $S_{4}$. We define $W^{L,3}$ which left multiplies $W^{L,2}_{p,p}W^{L,1}_{p,p}D_{p,q}W^{R}_{q,q}$ to replace the occurrences of $\ddot{q}_{j}$ in (25$a_{i}$) with their expression obtained in the previous step. The exceptions to $W^{L,3}_{p_{1},p_{2}}=0$ are:

[TABLE]

Step 5 - $S_{5}$. We define $W^{L,4}$ which left multiplies $W^{L,3}_{p,p}W^{L,2}_{p,p}W^{L,1}_{p,p}D_{p,q}W^{R}_{q,q}$ to replace the occurrences of $a_{j}$ in (25$f_{i}$) with their expression obtained in the previous step. The exceptions to $W^{L,4}_{p_{1},p_{2}}=0$ are:

[TABLE]

We are left with the definition of the permutations. Let’s first consider the $q$ permutation of the elements in $\mathcal{D}$. We choose:

[TABLE]

Let’s also choose a permutation $p_{fd}$ of the elements in $\mathcal{C}_{fd}$ :

[TABLE]