A data-efficient geometrically inspired polynomial kernel for robot inverse dynamics

TL;DR

This paper presents a novel polynomial kernel for Gaussian Process Regression that improves data efficiency and generalization in robot inverse dynamics estimation, reducing reliance on prior models.

Contribution

Introduction of the Geometrically Inspired Polynomial Kernel (GIP) for more data-efficient and less biased inverse dynamics estimation in robotics.

Findings

GIP kernel defines a finite-dimensional RKHS containing the inverse dynamics function.

The approach outperforms other data-driven estimators in data efficiency and generalization.

The method requires less prior information than traditional model-based estimators.

Abstract

In this paper, we introduce a novel data-driven inverse dynamics estimator based on Gaussian Process Regression. Driven by the fact that the inverse dynamics can be described as a polynomial function on a suitable input space, we propose the use of a novel kernel, called Geometrically Inspired Polynomial Kernel (GIP). The resulting estimator behaves similarly to model-based approaches as concerns data efficiency. Indeed, we proved that the GIP kernel defines a finite-dimensional Reproducing Kernel Hilbert Space that contains the inverse dynamics function computed through the Rigid Body Dynamics. The proposed kernel is based on the recently introduced Multiplicative Polynomial Kernel, a redefinition of the classical polynomial kernel equipped with a set of parameters that allows for a higher regularization. We tested the proposed approach in a simulated environment, and also in real experiments with a UR10 robot. The obtained results confirm that, compared to other data-driven estimators, the proposed approach is more data-efficient and exhibits better generalization properties. Instead, with respect to model-based estimators, our approach requires less prior information and is not affected by model bias.

Figures

Equations

$$\mathbf{\tau }=B\left(\mathbf{q}\right)\ddot{\mathbf{q}}+C\left(\mathbf{q},\dot{\mathbf{q}}\right)\dot{\mathbf{q}}+\mathbf{g}\left(\mathbf{q}\right)\text{,}$$

$$\mathbf{\tau }$$

$${\mathbf{y}}_i=\left[\begin{array}{c}{f}_i\left(\mathbf{x}\left(t_1\right)\right)\\ ⋮\\ f_i\left(\mathbf{x}\left(t_{N_{TR}}\right)\right)\end{array}\right]+\left[\begin{array}{c}{e}_i\left(t_1\right)\\ ⋮\\ e_i\left(t_{N_{TR}}\right)\end{array}\right]={\mathbf{f}}_i\left(X\right)+{\mathbf{e}}_i\text{,}$$

$${\hat{f}}_i({\mathbf{x}}_{\ast })=\sum _{h=1}^{N_{TR}}{\alpha }_i{k}_i({\mathbf{x}}_{\ast },\mathbf{x}(t_h))\text{,}$$

$$\mathbf{\alpha }={\left(K(X,X)+{\sigma }_{n_i}^{2}I\right)}^{-1}{\mathbf{y}}_i\text{.}$$

$$k_{RBF}(\mathbf{x}\left(t_h\right),\mathbf{x}\left(t_j\right))={\lambda }_{RBF}{e}^{||\mathbf{x}\left(t_h\right)-\mathbf{x}\left(t_j\right)|{|}_{{\Sigma }_{RBF}^{-1}}^2}\text{,}$$

$${\mathbf{m}}_i(X)={\Phi }_i(X)\bar{\mathbf{w}}$$

$$k_i(\mathbf{x}\left(t_h\right),\mathbf{x}\left(t_j\right))={\mathbf{\varphi }}_i\left(\mathbf{x}\left(t_h\right)\right){\Sigma }_w{\mathbf{\varphi }}_i^T\left(\mathbf{x}\left(t_j\right)\right)\text{,}$$

$$k_i(\mathbf{x}\left(t_h\right),\mathbf{x}\left(t_j\right))=$$

$${\mathbf{\varphi }}_i\left(\mathbf{x}\left(t_h\right)\right){\Sigma }_w{\mathbf{\varphi }}_i^T\left(\mathbf{x}\left(t_j\right)\right)+k_{RBF}(\mathbf{x}\left(t_h\right),\mathbf{x}\left(t_j\right)),$$

$${\mathbf{q}}_c=\left[\cos \left(q_{i_{r_1}}\right),\dots ,\cos \left(q_{i_{r_{N_r}}}\right)\right]\in {\mathbb{R}}^{N_r}\text{,}$$

$${\mathbf{q}}_s=\left[\sin \left(q_{i_{r_1}}\right),\dots ,\sin \left(q_{i_{r_{N_r}}}\right)\right]\in {\mathbb{R}}^{N_r}\text{,}$$

$${\mathbf{q}}_p=\left[q_{i_{p_1}},\dots ,q_{i_{p_{N_p}}}\right]\in {\mathbb{R}}^{N_p}\text{.}$$

$$\left\{{\dot{q}}_i{\dot{q}}_j,1\le i\le n,i\le j\le n\right\},$$

$$\bar{\mathbf{x}}=\left[{\mathbf{q}}_c,{\mathbf{q}}_s,{\mathbf{q}}_p,{\dot{\mathbf{q}}}_v,\ddot{\mathbf{q}}\right]\text{.}$$

$$deg\left(q_{c_b}\right)+deg\left(q_{s_b}\right)\le 2\text{ .}$$

$$k_{pk}^{(p)}\left(\mathbf{x}\left(t_h\right),\mathbf{x}\left(t_j\right)\right)={\left({\sigma }_p^2+{\mathbf{x}}^T\left(t_h\right){\Sigma }_p\mathbf{x}\left(t_j\right)\right)}^p\text{,}$$

$$k_{mpk}^{(p)}\left(\mathbf{x}\left(t_h\right),\mathbf{x}\left(t_j\right)\right)=\prod _{s=1}^p\left({\sigma }_s^2+{\mathbf{x}}^T\left(t_h\right){\Sigma }_s\mathbf{x}\left(t_j\right)\right)\text{,}$$

$$k_i\left(\bar{\mathbf{x}}\left(t_h\right),\bar{\mathbf{x}}\left(t_j\right)\right)$$

$$k_{cs}\left({\mathbf{q}}_{cs}\left(t_h\right),{\mathbf{q}}_{cs}\left(t_j\right)\right)=\prod _{b=1}^{N_r}{k}_{mpk}^{(2)}\left({\mathbf{q}}_{c{s}_b}\left(t_h\right),{\mathbf{q}}_{c{s}_b}\left(t_j\right)\right)\text{,}$$

$$k_p\left({\mathbf{q}}_p\left(t_h\right),{\mathbf{q}}_p\left(t_j\right)\right)=\prod _{b=1}^{N_p}{k}_{mpk}^{(2)}\left(q_{p_b}\left(t_h\right),q_{p_b}\left(t_j\right)\right)\text{.}$$

$$MSE({\mathbf{y}}_i,{\hat{\mathbf{y}}}_i)=\sum _{j=1}^N{\left(y_i(t_j)-{\hat{y}}_i(t_j)\right)}^2/N\text{,}$$

$$nMSE({\mathbf{\tau }}_i,{\hat{\mathbf{y}}}_i)=MSE({\mathbf{\tau }}_i,{\hat{\mathbf{y}}}_i)/Var\left({\mathbf{\tau }}_i\right)\text{.}$$

$$K_{eq}\mathbf{i}=B_{eq}\left(\mathbf{q}\right)\ddot{\mathbf{q}}+C\left(\mathbf{q},\dot{\mathbf{q}}\right)\dot{\mathbf{q}}+\mathbf{g}\left(\mathbf{q}\right)+F_v\dot{\mathbf{q}}+F_{c}sign\left(\dot{\mathbf{q}}\right)\text{,}$$

$$B\left(\mathbf{q}\right)=\sum _{i=1}^{i=n}{m}_i{J}_i^T{J}_i+J_{{\omega }_i}^T{R}_i^0{I}_i^i{R}_0^i{J}_{{\omega }_i}\text{,}$$

$$R_i^{i-1}=R_z\left({\theta }_i\right)R_x\left({\alpha }_i\right)\text{,}$$

$${\mathbf{l}}_i^{i-1}={\left[0,0,d_i\right]}^T+R_z\left({\theta }_i\right){\left[a_i,0,0\right]}^T\text{,}$$

$${\mathbf{\omega }}_i$$

$$J_{{\omega }_i}=\left[R_0^0\left[\begin{array}{c}0\\ 0\\ {\lambda }_1\end{array}\right],\dots ,R_{i-1}^0\left[\begin{array}{c}0\\ 0\\ {\lambda }_i\end{array}\right],0\left(3,n-i\right)\right]\text{.}$$

$$R_0^i{J}_{{\omega }_i}=\left[R_0^i\left[\begin{array}{c}0\\ 0\\ {\lambda }_1\end{array}\right],\dots ,R_{i-1}^i\left[\begin{array}{c}0\\ 0\\ {\lambda }_1\end{array}\right],0\left(3,n-i\right)\right]\dot{\mathbf{q}}\text{.}$$

$$deg\left(q_{c_b}\right)+deg\left(q_{s_b}\right)\le 2\text{,}$$

$$deg\left(q_{c_j}\right)+deg\left(q_{s_j}\right)\le 1\text{.}$$

$$\sum _{j=1}^{j=n}{c}_{ij}{\dot{q}}_j=\sum _{j=1}^{j=n}\sum _{k=1}^{k=n}\left(\frac{\partial b_{ij}}{\partial q_k}-\frac{1}{2}\frac{\partial b_{ik}}{\partial q_i}\right){\dot{q}}_k{\dot{q}}_j\text{.}$$

Full text

I INTRODUCTION

Learning the inverse dynamics model of a robot directly from data is still a challenging task in robotics, worth of investigation, as demonstrated by several important applications. For instance, by learning such a model, it is possible to design robot controllers based on feed-forward strategies [1] and on more complex Model Predictive Control approaches [2], or to provide robots with proprioceptive sensing capabilities [3], [4].

Learning models directly from data has several advantages. Firstly, the derivation of a model is not always an easy task, and, even when a model is available, its use introduces a bias, due to uncertainties on the values of parameters which are assumed known, or to assumptions which are just a rough approximation of the real behavior of the robot. Secondly, data-driven approaches are not platform-dependent, namely, the same learning technique can be applied to different physical platforms, leading to considerable advantages in terms of design time and costs.

Several data-driven strategies to learn inverse dynamics have been developed. In [5], the authors proposed a locally weighted projection of different linear models. A significant number of approaches rely on neural networks; for instance, the authors in [6] resort to the use of a recurrent neural network, while in [7] an LSTM network has been proposed. Another wide class of solutions is based on Gaussian Process Regression (GPR), [8], [9], [10] and [11]. Differently, from neural networks, GPR provides also a bound on the uncertainty of the estimates; this additional information can be exploited in different ways, see for instance in Reinforcement Learning the PILCO algorithm [12].

Although data-driven modeling techniques have been applied successfully in several control applications, see for example [12], [13], [14], they are still not able to guarantee the same generalization properties of model-based learning techniques. Indeed, data-driven approaches capture only similarity between data, without exploiting important features, like causality or the presence of constraints imposed by physics and geometry. This fact results in a considerable data inefficiency, which is particularly evident in systems with a high number of degrees of freedom (DOF). The typical huge amount of data required by standard data-driven approaches poses serious limitations on their applicability, mainly due to the high computational burden needed to process all the available information, in addition to the difficulty of guaranteeing good generalization properties.

In this paper, we investigate the possibility of developing data-driven estimators of robot inverse dynamics exhibiting good generalization properties and high data efficiency. The main contribution of the paper is the design of a data-driven inverse dynamics estimator based on GPR, more precisely on a novel kernel function, named Geometrically Inspired Polynomial Kernel (GIP). The main idea supporting our approach is related to the existence of a suitable transformation of the standard inputs, that are, positions, velocities, and accelerations, of the generalized coordinates, into an augmented space, where the inverse dynamics map derived with the Lagrangian equations is a polynomial function. Inspired by this property, we propose a model based on the Multiplicative Polynomial Kernel (MPK), recently introduced in [15], which is a re-parameterization of the standard polynomial kernel. As shown in [15], compared to the standard polynomial kernel, the MPK parametrization allows for greater flexibility in neglecting eventual unnecessary basis functions of the corresponding Reproducing Kernel Hilbert Space (RKHS), leading to higher generalization performance.

The idea of mapping the standard inputs in a space where the inverse dynamics is a polynomial function has been explored also in [16]. In [16], the authors have introduced the same augmented space we consider in this paper but have modeled only the elements of the inertia matrix and the potential energy as polynomials in this augmented space. However, instead of exploiting the encoding properties of polynomial kernels, they have designed a structural network inspired by the Lagrangian equations, where each element of the inertia matrix and the potential energy are a linear combination of monomials belonging to a particular class of polynomials. As consequence, compared to our approach, in the model proposed in [16] the number of parameters to be identified grows quickly w.r.t. the robot DOF, i.e., w.r.t. the number of possible monomials, increasing the risk of overfitting.

The polynomial-based strategy we introduce is tested both in a simulated environment and with data acquired from real experiments on a UR10 robot. Despite the GIP estimator requires minimal prior information compared to model-based estimators, the obtained results show that the proposed approach exhibits comparable performance in terms of accuracy and generalization. Additionally, compared to data-driven approaches, our learning algorithm is more data-efficient and exhibits better generalization properties.

The paper is organized as follows. In Section II, we provide an overview of the main strategies based on GPR adopted in inverse dynamics learning. In Section III we describe the approach we propose. Firstly, we identify an input transformation that leads to a description of the rigid body dynamics equations in terms of polynomial functions. Secondly, we briefly review MPK. Thirdly, we define the GIP kernel. Finally, in Section IV, we test the proposed estimator in a simulated environment, representing a SCARA robot, and on data coming from real experiments performed with a UR10 robot.

II ROBOT INVERSE DYNAMICS: LEARNING STRATEGIES

In this section, we briefly review the dynamics model of robot manipulators and the main approaches proposed to deal with the inverse dynamics problem.

Consider a robot manipulator with $n+1$ links and $n$ joints, and let $\boldsymbol{q}=\left[q_{1},\ldots,q_{n}\right]^{T}\in\mathbb{R}^{n}$ be the vector collecting the generalized coordinates associated to the joints; accordingly, let $\boldsymbol{\dot{q}}$ and $\boldsymbol{\ddot{q}}$ be the velocity and acceleration vectors, respectively. The inverse dynamics problem consists in estimating the function mapping the triple $\left(\boldsymbol{q},\boldsymbol{\dot{q}},\boldsymbol{\ddot{q}}\right)$ into the vector of generalized torques, denoted by $\boldsymbol{\tau}\in\mathbb{R}^{n}$. The estimation is typically performed starting from a set of input-output observations, which is composed by the set of input locations $X=\left\{\boldsymbol{x}(t_{1}),\dots,\boldsymbol{x}(t_{N_{TR}})\right\}$, where $\boldsymbol{x}(t)=\left[\boldsymbol{q}(t)\,,\boldsymbol{\dot{q}}(t)\,,\boldsymbol{\ddot{q}}(t)\right]$, and the set of outputs $Y$, corresponding to the the noisy measures of $\left\{\boldsymbol{\tau}(t_{1}),\dots,\boldsymbol{\tau}(t_{N_{TR}})\right\}$; $N_{TR}$ is the total number of observations. In the following, when there is no risk of confusion, we will omit the dependence on time $t$.

II-A Rigid body dynamics estimators

Several approaches which have been proposed to deal with the inverse dynamics problem are based on the rigid body dynamics (RBD) assumption. Under this assumption, the robot dynamics is described as

[TABLE]

where $B\left(\boldsymbol{q}\right)\in\mathbb{R}^{n\times n}$ and $C\left(\boldsymbol{q},\boldsymbol{\dot{q}}\right)\in\mathbb{R}^{n\times n}$ are, respectively, the inertia matrix and the Coriolis matrix, and $\boldsymbol{g}\left(\boldsymbol{q}\right)$ is the vector accounting for the gravitational contributions, see [17]. The previous equation depends on two sets of parameters, the kinematic and dynamics parameters. The first set is composed of the geometric quantities (i.e., lengths, angles) that, together with $\boldsymbol{q}$, define the forward kinematic. The second set, instead, contains the masses, centers of mass, and inertia components of the links. Remarkably, it is possible to show that (1) is linear w.r.t. the dynamics parameters, see [17]. Specifically, denoting by $\boldsymbol{w}$ the vector collecting all the dynamics parameters, (1) can be rewritten as

[TABLE]

for a suitable matrix $\Phi\in\mathbb{R}^{n\times N_{par}}$, which depends only on the kinematic parameters. Then, assuming the kinematic parameters to be known, the inverse dynamics problem boils down to the computation of $\hat{\boldsymbol{w}}$, an estimate of $\boldsymbol{w}$.

In several solutions, $\hat{\boldsymbol{w}}$ is computed relying on Fisherian techniques, see for example [18]. When the model is accurate, and the signal to noise ratio is sufficiently high, these estimators achieve accurate estimates, together with good generalization properties. However, besides the presence of noise, several aspects could limit the performance of such approaches. Indeed, it is worth stressing that errors about the kinematic parameters introduce model bias. Moreover, there are situations where it is hard to derive (2), or where the rigid body assumption is a too rough approximation of the real robot behaviors.

II-B Gaussian Process Regression for robot inverse dynamics

To overcome the limitations characterizing estimators based on the RBD assumption, several Bayesian approaches have been proposed in the last decade. Most techniques are based on GPR, see [19] for a detailed description. Typically, in GPR approaches, each joint is treated individually and modeled as a single Gaussian Process. More precisely, when considering the i-th joint, it is assumed that the output measurements $\boldsymbol{y}_{i}=\left\{\tau_{i}\left(t_{1}\right)\,,\dots\,,\tau_{i}\left(t_{N_{TR}}\right)\right\}$ are generated by the following probabilistic model

[TABLE]

where $\boldsymbol{e}_{i}$ is i.i.d. Gaussian noise with standard deviation $\sigma_{n_{i}}$, and $\boldsymbol{f}_{i}\left(X\right)$ is an unknown function defined as a Gaussian Process, namely, $\boldsymbol{f}_{i}\left(X\right)\sim N\left(\boldsymbol{m}_{i}\left(X\right)\,,K_{i}\left(X,X\right)\right)$, being $\boldsymbol{m}_{i}$ and $K_{i}\left(X,X\right)$, respectively, mean and covariance. In particular, the matrix $K_{i}\left(X,X\right)$, called also kernel matrix, is defined through a kernel function $k_{i}\left(\cdot,\cdot\right)$, i.e., the element in $h$-th row and $j$-th column is equal to $k_{i}\left(\boldsymbol{x}\left(t_{h}\right),\boldsymbol{x}\left(t_{j}\right)\right)$.

In this probabilistic framework, the maximum a posteriori probability (MAP) estimator of $f_{i}$ is given in closed form. Let $\boldsymbol{x}_{*}$ be a general input location. Then, as proved in [19] (see Chapter 2), the MAP estimate of $f_{i}(\boldsymbol{x}_{*})$ is given by

[TABLE]

where $\alpha_{i}$ is the i-th element of the vector $\boldsymbol{\alpha}$, defined as

[TABLE]

Before providing a brief overview about GPR based solutions, we comment about the computational cost of (4) and (5). Due to the matrix inversion in (5), the number of operations required to derive the $f_{i}$ estimator is proportional to the cube of $N_{TR}$. However, observe that $\boldsymbol{\alpha}$ depends only on the training data, then its computation can be performed offline. As consequence, once $\boldsymbol{\alpha}$ is computed, the cost of evaluating (4), grows linearly w.r.t. $N_{TR}$, thus allowing the use of GPR estimators in real time applications, see for example [14] and [11].

When no prior knowledge about the model is available, the GPR prior can be defined in a data-driven way. An option consists in assuming that the distribution of the outputs is stationary with respect to the locations of the inputs, and then, considering $\boldsymbol{m}_{i}\left(\cdot\right)=0$, and adopting a Radial Basis Function (RBF) as kernel. The RBF kernel is defined as

[TABLE]

where $\lambda_{RBF}$ and $\Sigma_{RBF}$ are the kernel hyperparameters, which are typically tuned from data by Marginal Likelihood (ML) maximization. In general, estimators based on RBF kernels well approximate the inverse dynamics only in a neighborhood of the training input locations. Several strategies have been designed to limit the computational complexity and to increase the generalization, see [11]. Despite that, when considering robots with a considerable number of DOF, it is still hard to design inverse dynamics estimators based on RBF kernel with remarkable generalization properties.

In case a RBD model is given, starting from (2), and modeling $\boldsymbol{w}$ as a Gaussian variable, with mean $\bar{\bf{w}}$ and covariance $\Sigma_{w}$, it is possible to derive a linear kernel that inherits all the positive aspects of the RBD estimators, but acts in a Bayesian framework. Consider joint i, then

[TABLE]

where $\Phi_{i}(X)$ is the matrix collecting all the rows $\boldsymbol{\phi}_{i}\left(\boldsymbol{x}\left(t_{j}\right)\right)$, $j=1,\ldots,N_{TR}$, and where, in this case, $\bar{\bf{w}}$ and $\Sigma_{w}$ are the kernel hyperparameters. The above kernel can be used alone, leading to the so-called Parametric Prior (PP) estimators, or together with a data-driven kernel, leading to the so-called Semiparametric Prior (SP) estimators. In the latter case, when adopting the RBF kernel as data-driven kernel, we have

[TABLE]

see for example [9], [10] and [20]. The rationale behind the use of kernel in (7) is the following: the first term allows exploiting the prior knowledge coming from the RBD, providing generalization, while $k_{RBF}\left(\cdot,\cdot\right)$ improves the estimate in a neighborhood of the training locations, compensating for model bias or un-modeled behaviors. We remark that estimators based on (6) and (7) are model-based estimators, since their kernel functions are derived starting from (2).

III PROPOSED APPROACH: GEOMETRICALLY INSPIRED POLYNOMIAL KERNEL

In this section, we propose a novel kernel that allows estimating the inverse dynamics without requiring prior knowledge of the model, preserving good generalization and high accuracy. This section is organized as follows. Firstly, we state Proposition 1, which characterizes the inverse dynamics from the functional analysis point of view. Given the type of each joint, i.e. prismatic or revolute, Proposition 1 defines a transformation of the input $\boldsymbol{x}$ where the inverse dynamics is a polynomial function. Then, we briefly review the Multiplicative Polinomial Kernel, and, finally, we define the proposed kernel function, named Geometrically Inspired Polynomial kernel.

III-A Polynomial characterization of the rigid-body model

In the following, we restrict our study to manipulators where each joint is either revolute or prismatic. Let $N_{r}$ and $N_{p}$ be the number of revolute and prismatic joints, respectively, where $N_{r}+N_{p}=n$, and let us denote by $I_{r}=\left\{i_{r_{1}},\dots,i_{r_{N_{r}}}\right\}$ and $I_{p}=\left\{i_{p_{1}},\dots,i_{p_{N_{p}}}\right\}$ the sets containing, respectively, the revolute and prismatic joints indexes. We start our analysis by defining

[TABLE]

In the following, we denote by $q_{c_{b}}$ the element in $\boldsymbol{q}_{c}$ associated to joint $i_{r_{b}}$, i.e. $\cos(q_{i_{r_{b}}})$ (similar definitions hold for $q_{s_{b}}$ and $q_{p_{b}}$). For later convenience we define also $\boldsymbol{q}_{cs}\in\mathbb{R}^{2N_{r}}$, the vector obtained concatenating $\boldsymbol{q}_{c}$ and $\boldsymbol{q}_{s}$. In addition, we denote by $\boldsymbol{\dot{q}}_{v}$ the vector collecting the elements of the set

[TABLE]

that is the set containing all the possible pairwise products of components of $\boldsymbol{\dot{q}}$. Notice that $\boldsymbol{\dot{q}}_{v}\in\mathbb{R}^{n\left(n+1\right)/2}$.

Finally, we introduce a compact notation to identify a particular set of inhomogeneous polynomial functions. Let $\boldsymbol{a}$ be the vector containing the $m$ variables $a_{1},\ldots,a_{m}$. We denote by $\mathbb{P}_{[p]}\left(\boldsymbol{a}_{[p_{1}]}\right)$ the set of polynomial functions of degree not greater than $p$ defined over the variables in $\boldsymbol{a}$, such that each variable $a_{i}$ appears with degree not greater than $p_{1}$. Similar definitions hold in case the inputs set accounts for more input vectors.

Now we consider the transformation $F:\mathbb{R}^{3n}\to\mathbb{R}^{\gamma}$, with $\gamma=2N_{r}+N_{p}+n(n+1)/2+n$, which maps the input location $\boldsymbol{x}$ into the element $\bar{\boldsymbol{x}}\in\mathbb{R}^{\gamma}$, defined as

[TABLE]

We have the following result.

Proposition 1

Consider a manipulator with $n+1$ links and $n$ joints, divided in $N_{r}$ revolute joints and $N_{p}$ prismatic joints, subject to $n=N_{r}+N_{p}$. Then, the inverse dynamics of each joint obtained through the rigid body model in (1) belongs to $\mathbb{P}_{(2n+1)}\left(\boldsymbol{q}_{c(2)},\boldsymbol{q}_{s(2)},\boldsymbol{q}_{p(2)},\boldsymbol{\dot{q}}_{v(1)},\boldsymbol{\ddot{q}}_{(1)}\right)$. Namely, each $\tau_{i}\left(\cdot\right)$ is a polynomial function in $\bar{\boldsymbol{x}}$, of degree not greater than $2n+1$, such that: (i) each element of $\boldsymbol{q}_{c}$, $\boldsymbol{q}_{s}$ and $\boldsymbol{q}_{p}$ appears with degree not greater than $2$, and (ii) each element of $\boldsymbol{\dot{q}}_{v}$ and $\boldsymbol{\ddot{q}}$ appears with degree not greater than $1$. Moreover, for any monomial of the aforementioned polynomial, the sum of the $q_{c_{b}}$ and $q_{s_{b}}$ degrees is equal or lower than two, namely, it holds

[TABLE]

The proof is reported in the Appendix.

Remark 1

*The result stated in Proposition 1 is related to the modeling property used in [16], though some important differences are present. Indeed, in [16], the authors have modeled with polynomial functions the potential energy and the elements of the inertia matrix, while Proposition 1 establishes that the whole inverse dynamics is a polynomial function in the augmented space. Moreover, Proposition 1 provides more strict constraints on the degrees of $\cos(q_{i})$, $\sin(q_{i})$, thus better characterizing the maximum degree with which each variable can appear in the different monomials, and, in turn, decreasing the number of possible monomials. *

III-B Multiplicative Polynomial Kernel

From a functional analysis point of view, Proposition 1 states that the inverse dynamics derived through the RBD belongs to the finite dimensional space of polynomial functions. A suitable set of basis functions for this space is given by the set of all the monomials in $\mathbb{P}_{(2n+1)}\left(\boldsymbol{q}_{c(2)},\boldsymbol{q}_{s(2)},\boldsymbol{q}_{p(2)},\boldsymbol{\dot{q}}_{v(1)},\boldsymbol{\ddot{q}}_{(1)}\right)$, heareafter denoted by vector $\boldsymbol{\phi}_{\mathbb{P}}\left(\boldsymbol{x}\right)\in\mathbb{R}^{N_{\mathbb{P}}}$, being $N_{\mathbb{P}}$ its cardinality. Unfortunately, $N_{\mathbb{P}}$ grows rapidly with the number of joints. To provide a couple of examples, when considering a SCARA robot, i.e. $n=4$, we have that $N_{\mathbb{P}}=1647$, while, for a standard six DOF robot, like the UR10, we have that $N_{\mathbb{P}}=302615$. Clearly, when considering GPR based approaches, the computational and memory requirements induced by the dimension of $\boldsymbol{\phi}_{\mathbb{P}}\left(\boldsymbol{x}\right)$ prevent the possibility of adopting the linear kernel $k\left(\boldsymbol{x}\left(t_{k}\right),\boldsymbol{x}\left(t_{j}\right)\right)=\boldsymbol{\phi}_{\mathbb{P}}\left(\boldsymbol{x}\left(t_{j}\right)\right)\Sigma_{w}\boldsymbol{\phi}^{T}_{\mathbb{P}}\left(\boldsymbol{x}\left(t_{j}\right)\right)$, as done in (6).

A compact solution that allows overcoming this problem consists in assuming that the target function $\tau_{i}\left(\cdot\right)$ belongs to the Reproducing Kernel Hilbert Space (RKHS) associated to a polynomial kernel, see [19]. More precisely, when considering the space of inhomogeneus polynomials defined on the components of $\boldsymbol{x}\in\mathbb{R}^{d}$, with maximum degree $p$, the polynomial kernel is classically defined as

[TABLE]

where $\sigma^{2}_{p}>0$ and $\Sigma_{p}>0$ are the kernel hyperparameters, see [19]. Unfortunately, as highlighted in [19] (chapter 4.2.2), the kernel function in (9) is not widely used in regression problems, since it is prone to overfitting, in particular when considering high dimensional inputs and $p>2$, that is exactly the situation identified in Proposition 1.

A valid alternative to (9) is represented by the MPK, recently introduced in [15]. When considering the space of inhomogeneous polynomial with maximum degree $p$, the MPK is defined as the product of $p$ linear kernels,

[TABLE]

where the $\Sigma_{s}>0\in\mathbb{R}^{d\times d}$ matrices are distinct diagonal matrices. The diagonal elements, together with the parameters $\sigma^{2}_{s}$, compose the hyperparameters set of the MPK.

Observe that the RKHSs identified by (10) and (9) contains the same basis functions. However, as discussed in [15], (10) is equipped with a richer set of hyperparameters, that can be tuned by ML maximization, and allows a better selection of the monomials that highly influence the system output.

III-C Geometrically inspired polynomial kernel

In this subsection, we describe the GIP kernel we propose to model the robot inverse dynamics. Our approach requires minimal information since we assume to know only the joints type. We model each joint torque with a zero-mean Gaussian process, and, driven by Proposition 1, we assume that the inverse dynamics is a polynomial in the input space $\bar{\boldsymbol{x}}$, defined in (8). To comply with the constraints on the maximum degree of each term, we adopt a kernel function given by the product of $N_{r}+N_{p}+1$ kernels of the type defined in equation (10), where

• •

$N_{r}$ kernels have $p=2$ and each of them is defined on a $2$-dimensional input space given by $\boldsymbol{q}_{cs_{b}}=\left[q_{c_{b}},q_{s_{b}}\right]$, with $b\in I_{r}$;

• •

$N_{p}$ kernels have $p=2$ and each of them is defined on a $1$-dimensional input, given by one of the $\boldsymbol{q}_{p}$ components;

• •

a single kernel with $p=1$ defined on the input vector $\boldsymbol{q}_{av}=\left[\boldsymbol{\ddot{q}}\,,\,\boldsymbol{\dot{q}}_{v}\right]$.

The resulting kernel for the $i$-th joint is

[TABLE]

with

[TABLE]

In Figure 1 we reported a schematic representation of the GIP kernel.

IV EXPERIMENTAL RESULTS

We tested the proposed approach both in a simulated environment and in a real environment. Regarding technical aspects, we implemented all the considered algorithms in Python. To speed up algebraic operations, we largely exploited the functionalities provided by Pytorch [21]. The code111https://bitbucket.org/AlbertoDallaLibera/gip_kernel and the datasets222https://mega.nz/#F!fbBBSCCK!NwRs60ace05mTe2Ot5fz-Q are publicly available.

IV-A Simulated SCARA robot

To evaluate the benefits of the GIP kernel, we first tested the proposed approach in a simulated environment. We considered a SCARA robot, more precisely an AdeptOne Robot. The SCARA is a $4$ DOF robot manipulator, with three revolute joints (joint $1$, $2$ and $4$) and a prismatic joint (joint $3$). As far as data generation is concerned, joint torques were computed through (1), assuming complete knowledge of the joint trajectories. Equation (1) was derived using the python package Sympybotics333https://github.com/cdsousa/SymPyBotics.

IV-A1 Estimation accuracy

In the first experiment, we tested the estimators’ accuracy. The proposed approach is compared with both model-based and data-driven estimators. As far as model-based estimators are concerned, we implemented three solutions. The first estimator is a classical Fisherian estimator (FE), based on (2); in particular, we considered the implementation provided by Sympybotics, see [18] for details. The other two model-based estimators are the PP and SP kernel-based estimators, where the model-based component is defined as in (6). For all the three solutions, we computed the $\Phi$ matrix in (2) assuming the nominal kinematic parameters provided by the manufacturer. To account for behaviors due to model bias, we varied the kinematic parameters of the model generating data around the nominal values, so that the $\Phi$ matrix used by the FE and the PP and SP kernels is different from the one generating data. Parameters perturbations are uniform random variables, ranging $\left[-0.05,0.05\right]m$ for lengths, and $\left[-5,5\right]deg$ for angles. Instead, as far as data-driven approaches are concerned, we tested an RBF kernel-based estimator and a neural network. The neural network is a fully connected network with two hidden layers, each of which is composed of $400$ sigmoids. Recurrent architectures have not been considered since typically they are applied to deal with the on-line adaptation problem, which is out of the scope of this work.

To obtain statistically significant results, we performed a Monte Carlo analysis, composed of $20$ experiments. A single experiment consists in sampling a model perturbation and performing two simulations, one for the training set and one for the test set. In each simulation, joints follow a distinct random trajectory, given by the sum of $200$ sinusoids with random amplitudes and angular velocities sampled in the range $[-2,2]\,rad/sec$. A zero-mean Gaussian noise with standard deviation $0.01\,Nm$ was added to the output of (1), resulting in a high signal to noise ratio. Both the training and test dataset are composed of $2000$ samples.

The hyperparameters of the GPR based estimators were trained by ML maximization (see [19], chapter 2.2, for a detailed discussion). As concerns $\boldsymbol{\bar{w}}$ and $\Sigma_{w}$, i.e., the prior distribution of the model-based kernels, we considered $\boldsymbol{\bar{w}}=0$, and $\Sigma_{w}$ equals to a diagonal matrix with distinct diagonal elements. Instead, as concerns the optimization of neural network parameters, the Mean Squared Error (MSE) was considered as loss function, defined as

[TABLE]

with $N$ equals to the number of samples. Both for GPRs and the neural network, we used Adam as optimizer [22]. Further details about the optimization, size of the batch and number of epochs, are provided in the source code.

Performance are compared by Normalized Mean Squared Error ($nMSE$) in the test set, defined as

[TABLE]

In Figure 2, we have plotted the obtained $nMSE$s through a boxplot. Results show that the proposed approach outperforms other data-driven estimators, which are not able to learn accurately the inverse dynamics of the SCARA robot using just $2000$ samples. Indeed, except for joint 4, the $nMSE$s of the RBF kernel-based estimator and the neural network estimator are in most of the trials higher than $10\%$. Instead, the GIP kernel-based estimator provides accurate estimates, as proven by $nMSE$s values, that are always below $1\%$, with the exception of joint 4, where two outliers are present, probably due to training inputs not sufficiently exciting. Moreover, the GIP kernel-based estimator performs similarly to the model-based approaches. Actually, in joint $2$ and $3$, the proposed approach outperforms FE and the PP kernel-based estimator, whose performances are affected by model bias. Results confirm also the validity of semiparametric schemes, proving that the addition of a data-driven component can compensate for model bias, given that the SP kernel-based estimator outperforms PP. Anyway, we highlight that in hybrid schemes the RBF component might not be effective in compensating for model bias, in particular when the performance of the data-driven estimator is low, as proven by the $nMSE$s in joint $3$, where the GIP kernel-based estimator is more accurate than SP. Finally, a comparison of the FE and PP performance suggests that Fisherian approaches are more sensitive to model bias; notice, in particular, the $nMSE$s obtained in joint 1 and 4.

IV-A2 Data efficiency

In the second experiment, we tested the data efficiency of different estimators. Since our focus is on comparing data-driven approaches, model bias was not considered, in favor of greater results interpretability. The GIP kernel-based estimator is compared with the other data-driven estimators, and also with the PP kernel-based approach. In this ideal scenario, where data are generated with the robot nominal parameters, the performance of PP might be considered as the baseline of an optimal solution.

The experiment is composed of training and test simulation, with joints trajectories generated as in the previous experiment; each dataset accounts for $4000$ samples. Results are reported in Figure 3, where we have plotted the evolution of the Global Mean Squared Error ($GMSE$), i.e., the sum of the $MSE(\boldsymbol{\tau}_{i},\boldsymbol{\hat{y}}_{i})$s of the four joints, as function of the number of training samples used to train and derive estimators. The evolutions of the errors show that the proposed solution outperforms the other data-driven estimators, both in terms of accuracy and data efficiency, given that its $GMSE$ is lower and decreases faster. As in the previous experiment, the GIP kernel-based estimator behaves more similarly to the model-based approach than to the other data-driven solutions, proving its data efficiency.

IV-B UR10 robot

We used a Universal Robots UR10 to test the proposed approach in a real setup. The UR10 robot is a $6$ DOF collaborative manipulator, where all the joints are revolute. This robot is not equipped with joint torque sensors, but one can directly measure the motor currents $\boldsymbol{i}$. Assuming that the behaviors due to elasticity are negligible, i.e. $\boldsymbol{\theta}=K_{r}\boldsymbol{q}$, where $\boldsymbol{\theta}$ contains the motor angles, and $K_{r}$ is the diagonal matrix of the gear reduction ratio, the inverse dynamics in (1) can be rewritten as

[TABLE]

where $F_{v}+F_{c}sign\left(\boldsymbol{\dot{q}}\right)$ accounts for the motors frictions, and $B_{eq}\left(\boldsymbol{q}\right)=B\left(\boldsymbol{q}\right)+K_{r}^{2}B_{m}$, with $B_{m}$ equals to the diagonal matrix of the rotor inertias; the $K_{eq}$ matrix is defined as $K_{i}K_{r}$, where $K_{i}$ is the diagonal matrix containing the torque-current coefficients of the motors.

The interface with the robot is based on ROS [23], through the ur_modern_driver444https://github.com/ThomasTimm/ur_modern_driver, and data are acquired with a sampling time of $8\cdot 10^{-3}\,sec$. The driver provides joints positions, velocities, and currents, while accelerations are computed through causal numerical differentiation. The dataset collected is publicly available, and it has been designed to stress generalization properties. The training set accounts for $40000$ samples, collected through a random exploration of the robot workspace, requiring the end-effector to reach a series of random points with variable velocity. Instead, the test dataset is composed of two types of trajectories, for a total number of $25312$ points. $22324$ points have been collected through a random exploration, similar to the one described for the training dataset. The remaining samples come from the trajectory obtained requiring the end-effector to track a circle of radius $30\,cm$ at a tool speed of $30\,mm/s$.

The optimization procedures and the considered estimators are the same of the previous experiment. Due to space constraints, we neglect the FE, which performs similarly to PP. Given the higher complexity of the UR10 inverse dynamics, the number of hidden units of the neural network was increased to $600$. Concerning the derivation of the GPR based estimators, i.e. the computation of (4) and (5), we used a subset of the training data, in order to reduce the computational burden induced by (5); in particular, we selected $4000$ samples, downsampling with constant step the $40000$ training samples. The kinematic parameters considered in the derivation of the model-driven components are the nominal values provided by the manufactures.

The results obtained in the real setup, and reported in Figure 4, confirm the behaviors obtained in the simulative setup. The proposed approach outperforms the other data-driven estimators in all the joints, confirming that data efficiency is crucial to derive inverse dynamics estimators with good generalization properties. GIP performance is close to the ones of the model-based estimators, and in joints $5$ and $6$ the proposed approach slightly outperforms the PP estimator, that, as explained before, might be affected by model errors. SP performance confirms that model errors can be compensated by the data-driven component, even though, as proven by the $nMSE$ in joint 6, the improvement might not be so significant when data-driven estimates are not accurate.

V CONCLUSION

In this paper, we have introduced a novel polynomial kernel to deal with the data-driven inverse dynamics identification. Compared to other data-driven approaches, the proposed kernel, called GIP kernel, is more data-efficient. As proven by experiments in a simulated environment and also in a real system, this property allows deriving accurate inverse dynamics estimators without the need of prior knowledge about the model. Numerical results show that the GIP kernel-based estimator exhibits behaviors similar to the ones of model-based approaches, in terms of accuracy, generalization and data efficiency. However, compared to model-based solutions, our solution is not affected by model bias. Additionally, our algorithm is convenient from an implementation point of view, given its generality and hence the possibility of applying the same approach to different physical systems.

Appendix A Proof of the Proposition 1

We prove Proposition 1 by inspection, analyzing individually all the terms in (1), i.e., the $B\left(\boldsymbol{q}\right)\boldsymbol{\ddot{q}}$ and $C\left(\boldsymbol{q},\boldsymbol{\dot{q}}\right)\boldsymbol{\dot{q}}$ contributions and the gravity term $\boldsymbol{g}\left(\boldsymbol{q}\right)$. Firstly, we provide a characterization of the elements of $B\left(\boldsymbol{q}\right)$ as polynomials in $\boldsymbol{q}_{c}$, $\boldsymbol{q}_{s}$ and $\boldsymbol{q}_{p}$. The inertia matrix is given by

[TABLE]

where $m_{i}$ and $I_{i}^{i}$ are the i-th link mass and inertia matrix, expressed in a reference frame (RF) solidal with the i-th link. $J_{i}$ and $J_{\omega_{i}}$ are the linear and angular Jacobians of the i-th RF, i.e., $\boldsymbol{\dot{c}}_{i}=J_{i}\boldsymbol{\dot{q}}$ and $\boldsymbol{\omega}_{i}=J_{\omega_{i}}\boldsymbol{\dot{q}}$, where $\boldsymbol{c}_{i}$ is the position of the center of mass of the i-th link, while $\boldsymbol{\omega}_{i}$ is the angular velocity of the i-th RF. To expand the $J_{i}$ and $J_{\omega_{i}}$ expressions, we introduce some notions regarding the kinematics. Adopting the Denavit-Hartenberg (DH) convention, the $R_{i}^{i-1}$ and $\boldsymbol{l}_{i}^{i-1}$ variables, which denote the i-th RF orientation and translation with the respect to the previous RF, are

[TABLE]

where $R_{x}$ and $R_{z}$ are the elementary rotation matrices around the $x$ and $z$ axis, while $a_{i}$ and $\alpha_{i}$ are two constant geometrical parameters, see [17]. The definitions of $d_{i}$ and $\theta_{i}$ depend on the joint interconnecting the i-th link with the previous link. When the joint is revolute, $d_{i}$ is constant and $\theta_{i}=\theta_{0_{i}}+q_{i}$, and the only terms that depend on $\boldsymbol{q}$ are $\cos\left(q_{i}\right)$ and $\sin\left(q_{i}\right)$ contained in $R_{i}^{i-1}$. Referring to the polynomial notation previously introduced, we can write that the elements of $R_{i}^{i-1}$ are functions in $\mathbb{P}_{(1)}\left(\cos\left(q_{i}\right)_{(1)},\sin\left(q_{i}\right)_{(1)}\right)$. In case the joint is prismatic, $\theta_{i}$ is constant, and $d_{i}=d_{0_{i}}+q_{i}$. Consequently the only $\boldsymbol{q}$ dependent terms are in $\boldsymbol{l}_{i}^{i-1}$. In particular the elements of $\boldsymbol{l}_{i}^{i-1}$ belong to $\mathbb{P}_{(1)}\left(q_{i_{(1)}}\right)$.

The $J_{\omega_{i}}$ matrix relates $\boldsymbol{\dot{q}}$ with the angular velocity of the i-th link. Adopting the DH convention, $\boldsymbol{\omega}_{i}^{i-1}=\lambda_{i}\left[0\,,\,0\,,\,\dot{q}_{i}\right]^{T}$, with $\lambda_{i}=1$ if the joint is revolute, and $\lambda_{i}=0$ if it is prismatic. Then, summing all the angular velocities projected in the base frame through the $R_{j}^{0}=\prod_{b=j}^{b=0}R_{b}^{b-1}$ matrices, and remarking that $\boldsymbol{\omega}_{i}=\sum_{j=1}^{j=i}\lambda_{j}R_{j-1}^{0}\boldsymbol{\omega}_{j}^{j-1}$, we obtain

[TABLE]

where $\boldsymbol{0}\left(3,n-i\right)$ is a $3\times\left(n-i\right)$ matrix containing only zero elements. The last equation implies

[TABLE]

Exploiting the properties of the rotation matrices, we obtain

[TABLE]

Let $\left\{I_{r}\leq i\right\}$ be the set of revolute joint indexes lower or equal than i, and let $\boldsymbol{q}_{c}\left(\left\{I_{r}\leq i\right\}\right)$ be the corresponding subset. Recalling that the elements of $R_{i}^{i-1}$ are functions in $\mathbb{P}_{(1)}\left(\cos\left(q_{i}\right)_{(1)},\sin\left(q_{i}\right)_{(1)}\right)$, with maximal degree one, and that $R_{j}^{k}=\prod_{b=j}^{b=k}R_{b}^{b-1}$, with $j>k$, it follows that the $J_{\omega_{i}}^{T}R_{i}^{0}I_{i}^{i}R_{0}^{i}J_{\omega_{i}}$ elements belong to $\mathbb{P}_{(2\left|\left\{I_{r}\leq i\right\}\right|)}\left(\boldsymbol{q}_{c}\left(\left\{I_{r}\leq i\right\}\right)_{(2)},\boldsymbol{q}_{s}\left(\left\{I_{r}\leq i\right\}\right)_{(2)}\right)$, where in each monomial the following constraint holds

[TABLE]

To derive a similar characterization of the $J_{i}$ elements, we analyze the $\boldsymbol{c}_{i}$ expression. The position of the i-th center of mass in the base frame is $\boldsymbol{c}_{i}=\sum_{j=1}^{j=i-1}R_{j-1}^{0}\boldsymbol{l}_{j}^{j-1}+R_{i}^{0}\boldsymbol{c}_{i}^{i}$. Then the $\boldsymbol{c}_{i}$ elements are functions in $\mathbb{P}_{(i)}\left(\boldsymbol{q}_{c}\left(\left\{I_{r}\leq i\right\}\right)_{(1)},\boldsymbol{q}_{s}\left(\left\{I_{r}\leq i\right\}\right)_{(1)},\boldsymbol{q}_{p}\left(\left\{I_{p}\leq i\right\}\right)_{(1)}\right)$, and in each monomial the following inequality holds

[TABLE]

Since $\boldsymbol{\dot{c}}_{i}=J_{i}\boldsymbol{\dot{q}}$, and since the derivative of $\cos\left(q_{j}\right)$, $\sin\left(q_{j}\right)$ and $q_{j}$ does not increase the degree of these terms when inequality (13) holds, it follows that the $J_{i}$ elements belong to the same functional space of $\boldsymbol{c}_{i}$. Consequently, the elements of $J_{i}^{T}J_{i}$ are functions in $\mathbb{P}_{(2i)}\left(\boldsymbol{q}_{c}\left(\left\{I_{r}\leq i\right\}\right)_{(2)},\boldsymbol{q}_{s}\left(\left\{I_{r}\leq i\right\}\right)_{(2)},\boldsymbol{q}_{p}\left(\left\{I_{p}\leq i\right\}\right)_{(2)}\right)$; as before, in each monomial the $q_{c_{j}}$ and $q_{s_{j}}$ degrees are subject to inequality (12).

Given the characterization of $J_{i}^{T}J_{i}$ and $J_{\omega_{i}}^{T}R_{i}^{0}I_{i}^{i}R_{0}^{i}J_{\omega_{i}}$, we obtain that the $B\left(\boldsymbol{q}\right)$ elements are functions in $\mathbb{P}_{(2n)}\left(\boldsymbol{q}_{c_{(2)}},\boldsymbol{q}_{s_{(2)}},\boldsymbol{q}_{p_{(2)}}\right)$, where in each monomial the $q_{c_{j}}$ and $q_{s_{j}}$ degrees are subject to inequality (12). Then the $B\left(\boldsymbol{q}\right)\boldsymbol{\ddot{q}}$ are functions in $\mathbb{P}_{(2n+1)}\left(\boldsymbol{q}_{c_{(2)}},\boldsymbol{q}_{s_{(2)}},\boldsymbol{q}_{p_{(2)}},\boldsymbol{\ddot{q}}_{(1)}\right)$.

As reported in [17], the i-th element of the $C\left(\boldsymbol{q},\boldsymbol{\dot{q}}\right)\boldsymbol{\dot{q}}$ product is equal to

[TABLE]

Since the $B\left(\boldsymbol{q}\right)$ elements belong to $\mathbb{P}_{(2n)}\left(\boldsymbol{q}_{c_{(2)}},\boldsymbol{q}_{s_{(2)}},\boldsymbol{q}_{p_{(2)}}\right)$ and (12) holds true, also the $b_{ij}$ partial derivatives belong to $\mathbb{P}_{(2n)}\left(\boldsymbol{q}_{c_{(2)}},\boldsymbol{q}_{s_{(2)}},\boldsymbol{q}_{p_{(2)}}\right)$, with (12) satisfied. Indeed, for each monomial in $b_{ij}$, the derivation respect to $q_{k}$ with $k\in I_{p}$ decreases the degree by one, while the derivation respect to $q_{k}$ with $k\in I_{r}$ does not alter the monomial degree, being that each $q_{k}$ appears as $\cos(q_{k})$ or $\sin(q_{k})$. Then, we obtain that the elements of $C\left(\boldsymbol{q},\boldsymbol{\dot{q}}\right)\boldsymbol{\dot{q}}$ are functions in $\mathbb{P}_{(2n+1)}\left(\boldsymbol{q}_{c_{(2)}},\boldsymbol{q}_{s_{(2)}},\boldsymbol{q}_{p_{(2)}},\boldsymbol{\dot{q}}_{v_{(1)}}\right)$.

Regarding $\boldsymbol{g}\left(\boldsymbol{q}\right)$, we observe that the i-th element is given by $-\partial\emph{U}/\partial q_{i}$, where by definition the potential energy $\emph{U}=\sum_{j=1}^{j=n}\boldsymbol{g}_{0}^{T}\boldsymbol{c}_{j}$, with $\boldsymbol{g}_{0}$ denoting the vector of the gravitational acceleration. Then, the elements of $\boldsymbol{g}\left(\boldsymbol{q}\right)$ are functions in the same space of the $J_{n}$ elements.

To conclude the proof, we just need to sum all the contributions and to note that, for each joint, the torque is a function in $\mathbb{P}_{(2n+1)}\left(\boldsymbol{q}_{c_{(2)}},\boldsymbol{q}_{s_{(2)}},\boldsymbol{q}_{p_{(2)}},\boldsymbol{\dot{q}}_{v_{(1)}},\boldsymbol{\ddot{q}}_{(1)}\right)$, with each monomial satisfying (12).