# Input Decoupling of Lagrangian Systems via Coordinate Transformation: General Characterization and Its Application to Soft Robotics

Pietro Pustina, Cosimo Della Santina, Frédéric Boyer, Alessandro de Luca, Federico Renda

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Pietro Pustina, Cosimo Della Santina, Frédéric Boyer, Alessandro de Luca, Federico Renda. Input Decoupling of Lagrangian Systems via Coordinate Transformation: General Characterization and Its Application to Soft Robotics. IEEE Transactions on Robotics, 2024, 40, pp.2098-2110. 10.1109/tro.2024.3370089. hal-04591947

HAL Id: hal-04591947

https://imt-atlantique.hal.science/hal-04591947

Submitted on 29 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Input Decoupling of Lagrangian Systems via Coordinate Transformation: General Characterization and Its Application to Soft Robotics

Pietro Pustina, Cosimo Della Santina, Frédéric Boyer, Alessandro de Luca, Federico Renda

<figure><img src="image_1.png" /></figure>

## To cite this version:

Pietro Pustina, Cosimo Della Santina, Frédéric Boyer, Alessandro de Luca, Federico Renda. Input Decoupling of Lagrangian Systems via Coordinate Transformation: General Characterization and Its Application to Soft Robotics. IEEE Transactions on Robotics, 2024, 40, pp.2098-2110. 10.1109/tro.2024.3370089. hal-04591947

HAL Id: hal-04591947

https://imt-atlantique.hal.science/hal-04591947v1

Submitted on 29 May 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Input Decoupling of Lagrangian Systems via Coordinate Transformation: General Characterization and Its Application to Soft Robotics

Pietro Pustina (1), Graduate Student Member, IEEE, Cosimo Della Santina (1), Senior Member, IEEE, Frédéric Boyer, Member, IEEE, Alessandro De Luca (1), Life Fellow, IEEE, and Federico Renda (1), Member, IEEE

**Abstract**—Suitable representations of dynamical systems can simplify their analysis and control. On this line of thought, this article aims to answer the following question: *Can a transformation of the generalized coordinates under which the actuators directly perform work on a subset of the configuration variables be found?* We not only show that the answer to this question is yes but also provide necessary and sufficient conditions. More specifically, we look for a representation of the configuration space such that the right-hand side of the dynamics in Euler–Lagrange form becomes $[I O]^T u$, being $u$ the system input. We identify a class of systems, called *collocated*, for which this problem is solvable. Under mild conditions on the input matrix, a simple test is presented to verify whether a system is collocated or not. By exploiting power invariance, we provide necessary and sufficient conditions that a change of coordinates decouples the input channels if and only if the dynamics is collocated. In addition, we use the collocated form to derive novel controllers for damped underactuated mechanical systems. To demonstrate the theoretical findings, we consider several Lagrangian systems with a focus on continuum soft robots.

Manuscript received 12 June 2023; revised 9 November 2023; accepted 23 January 2024. Date of publication 27 February 2024; date of current version 12 March 2024. This paper was recommended for publication by Associate Editor Fabio Ruggiero and Editor Paolo Robuffo Giordano upon evaluation of the reviewers' comments. The work of Pietro Pustina and Federico Renda was supported in part by the US Office of Naval Research Global under Grant N62909-21-1-2033 and in part by the Khalifa University of Science and Technology under Grant CIRA-2020-074 and Grant RC1-2018-KUCARS. The work of Cosimo Della Santina was supported by the Horizon Europe Program from Project EMERGE under Grant 101070918. The work of Alessandro De Luca was supported by the PNRR MUR Project under Grant PE0000013-FAIR. (Corresponding author: Pietro Pustina.)

Pietro Pustina is with the Department of Computer, Control and Management Engineering, Sapienza University of Rome, 00185 Rome, Italy, and also with the Department of Cognitive Robotics, Delft University of Technology, 2628CN Delft, The Netherlands (e-mail: pustina@diag.uniroma1.it).

Cosimo Della Santina is with the Department of Cognitive Robotics, Delft University of Technology, 2628CN Delft, The Netherlands, and also with the Institute of Robotics and Mechatronics, German Aerospace Center (DLR), 82234 Oberpfaffenhofen, Germany (e-mail: c.dellasantina@tudelft.nl).

Frédéric Boyer is with the Institute Mines Telecom Atlantique, 44307 Nantes, France (e-mail: frederic.boyer@imt-atlantique.fr).

Alessandro De Luca is with the Department of Computer, Control and Management Engineering, Sapienza University of Rome, 00185 Rome, Italy (e-mail: deluca@diag.uniroma1.it).

Federico Renda is with the Khalifa University Center for Autonomous Robotics System (KUCARS) and the Department of Mechanical and Nuclear Engineering, Khalifa University of Science and Technology, Abu Dhabi 127788, United Arab Emirates (e-mail: federico.renda@ku.ac.ae).

Digital Object Identifier 10.1109/TRO.2024.3370089

**Index Terms**—Control, dynamics, learning for soft robots, modeling, motion control, underactuated robots.

<table><thead><tr><th></th><th>NOMENCLATURE</th></tr></thead><tbody><tr><td>Symbol</td><td>Description.</td></tr><tr><td>R<sup>n</sup></td><td>Euclidean space of dimension n.</td></tr><tr><td>S<sup>n</sup></td><td>Unit sphere of dimension n.</td></tr><tr><td>R<sup>n×m</sup></td><td>Space of n × m matrices over R.</td></tr><tr><td>X</td><td>Manifold of interest with X = {M, N}.</td></tr><tr><td>R<sup>≥0</sup></td><td>Positive real numbers n.</td></tr><tr><td>T<sub>v</sub>X</td><td>Tangent space of manifold X at v ∈ X.</td></tr><tr><td>B(v)</td><td>Neighborhood of v ∈ X.</td></tr><tr><td>s<sub>e</sub>(3)</td><td>Special Euclidean algebra of dimension 3.</td></tr><tr><td>s<sub>o</sub>(3)</td><td>Special orthogonal algebra of dimension 3.</td></tr><tr><td>I<sub>n</sub> ∈ R<sup>n×n</sup></td><td>Identity matrix of dimension n.</td></tr><tr><td>O<sub>n×m</sub></td><td>Zero matrix of dimension n × m.</td></tr><tr><td>P</td><td>Symmetric positive-definite matrix.</td></tr><tr><td>S<sub>i</sub> ∈ R<sup>n</sup></td><td>Column i of matrix S ∈ R<sup>n×m</sup>.</td></tr><tr><td>S<sub>ij</sub> ∈ R</td><td>Element in row i and column j of S.</td></tr><tr><td>[v]<sub>i</sub> ∈ R<sup>i</sup></td><td>Vector containing the first i components of v ∈ R<sup>n</sup>, with i ≤ n.</td></tr><tr><td>||v||</td><td>Euclidean norm of v.</td></tr><tr><td>r&#773; ∈ s<sub>o</sub>(3)</td><td>Skew-symmetric matrix defined by r ∈ R<sup>3</sup>.</td></tr><tr><td>&alpha;&#770; = (β&#770; &gamma;&#770;)</td><td>Tensor representation of &alpha; = (&beta;<sup>T</sup> &gamma;<sup>T</sup>)<sup>T</sup> with &beta;, &gamma; ∈ R<sup>3</sup>.</td></tr><tr><td>J<sub>f</sub>(x) = &part;<sub>f</sub>/&part;x ∈ R<sup>l×h</sup></td><td>Jacobian of the vector function f(x) : R<sup>h</sup> → R<sup>l</sup>.</td></tr><tr><td>tanh(v)</td><td>Vector obtained by applying tanh(·) componentwise to v.</td></tr></tbody></table>

## I. INTRODUCTION

ELECTRICAL, hydraulic, and mechanical systems, or their combinations, are Lagrangian systems that usually exhibit complex behavior. However, their physical nature displays special properties, such as symmetry and passivity, which have been
# Input Decoupling of Lagrangian Systems via Coordinate Transformation: General Characterization and Its Application to Soft Robotics

Pietro Pustina (1), Graduate Student Member, IEEE, Cosimo Della Santina (1), Senior Member, IEEE, Frédéric Boyer, Member, IEEE, Alessandro De Luca (1), Life Fellow, IEEE, and Federico Renda (1), Member, IEEE

**Abstract**—Suitable representations of dynamical systems can simplify their analysis and control. On this line of thought, this article aims to answer the following question: *Can a transformation of the generalized coordinates under which the actuators directly perform work on a subset of the configuration variables be found?* We not only show that the answer to this question is yes but also provide necessary and sufficient conditions. More specifically, we look for a representation of the configuration space such that the right-hand side of the dynamics in Euler–Lagrange form becomes $[I O]^T u$, being $u$ the system input. We identify a class of systems, called *collocated*, for which this problem is solvable. Under mild conditions on the input matrix, a simple test is presented to verify whether a system is collocated or not. By exploiting power invariance, we provide necessary and sufficient conditions that a change of coordinates decouples the input channels if and only if the dynamics is collocated. In addition, we use the collocated form to derive novel controllers for damped underactuated mechanical systems. To demonstrate the theoretical findings, we consider several Lagrangian systems with a focus on continuum soft robots.

Manuscript received 12 June 2023; revised 9 November 2023; accepted 23 January 2024. Date of publication 27 February 2024; date of current version 12 March 2024. This paper was recommended for publication by Associate Editor Fabio Ruggiero and Editor Paolo Robuffo Giordano upon evaluation of the reviewers' comments. The work of Pietro Pustina and Federico Renda was supported in part by the US Office of Naval Research Global under Grant N62909-21-1-2033 and in part by the Khalifa University of Science and Technology under Grant CIRA-2020-074 and Grant RC1-2018-KUCARS. The work of Cosimo Della Santina was supported by the Horizon Europe Program from Project EMERGE under Grant 101070918. The work of Alessandro De Luca was supported by the PNRR MUR Project under Grant PE0000013-FAIR. (Corresponding author: Pietro Pustina.)

Pietro Pustina is with the Department of Computer, Control and Management Engineering, Sapienza University of Rome, 00185 Rome, Italy, and also with the Department of Cognitive Robotics, Delft University of Technology, 2628CN Delft, The Netherlands (e-mail: pustina@diag.uniroma1.it).

Cosimo Della Santina is with the Department of Cognitive Robotics, Delft University of Technology, 2628CN Delft, The Netherlands, and also with the Institute of Robotics and Mechatronics, German Aerospace Center (DLR), 82234 Oberpfaffenhofen, Germany (e-mail: c.dellasantina@tudelft.nl).

Frédéric Boyer is with the Institute Mines Telecom Atlantique, 44307 Nantes, France (e-mail: frederic.boyer@imt-atlantique.fr).

Alessandro De Luca is with the Department of Computer, Control and Management Engineering, Sapienza University of Rome, 00185 Rome, Italy (e-mail: deluca@diag.uniroma1.it).

Federico Renda is with the Khalifa University Center for Autonomous Robotics System (KUCARS) and the Department of Mechanical and Nuclear Engineering, Khalifa University of Science and Technology, Abu Dhabi 127788, United Arab Emirates (e-mail: federico.renda@ku.ac.ae).

Digital Object Identifier 10.1109/TRO.2024.3370089

**Index Terms**—Control, dynamics, learning for soft robots, modeling, motion control, underactuated robots.

<table><thead><tr><th></th><th>NOMENCLATURE</th></tr></thead><tbody><tr><td>Symbol</td><td>Description.</td></tr><tr><td>R<sup>n</sup></td><td>Euclidean space of dimension n.</td></tr><tr><td>S<sup>n</sup></td><td>Unit sphere of dimension n.</td></tr><tr><td>R<sup>n×m</sup></td><td>Space of n × m matrices over R.</td></tr><tr><td>X</td><td>Manifold of interest with X = {M, N}.</td></tr><tr><td>R<sup>≥0</sup></td><td>Positive real numbers n.</td></tr><tr><td>T<sub>v</sub>X</td><td>Tangent space of manifold X at v ∈ X.</td></tr><tr><td>B(v)</td><td>Neighborhood of v ∈ X.</td></tr><tr><td>s<sub>e</sub>(3)</td><td>Special Euclidean algebra of dimension 3.</td></tr><tr><td>s<sub>o</sub>(3)</td><td>Special orthogonal algebra of dimension 3.</td></tr><tr><td>I<sub>n</sub> ∈ R<sup>n×n</sup></td><td>Identity matrix of dimension n.</td></tr><tr><td>O<sub>n×m</sub></td><td>Zero matrix of dimension n × m.</td></tr><tr><td>P</td><td>Symmetric positive-definite matrix.</td></tr><tr><td>S<sub>i</sub> ∈ R<sup>n</sup></td><td>Column i of matrix S ∈ R<sup>n×m</sup>.</td></tr><tr><td>S<sub>ij</sub> ∈ R</td><td>Element in row i and column j of S.</td></tr><tr><td>[v]<sub>i</sub> ∈ R<sup>i</sup></td><td>Vector containing the first i components of v ∈ R<sup>n</sup>, with i ≤ n.</td></tr><tr><td>||v||</td><td>Euclidean norm of v.</td></tr><tr><td>r&#773; ∈ s<sub>o</sub>(3)</td><td>Skew-symmetric matrix defined by r ∈ R<sup>3</sup>.</td></tr><tr><td>&alpha;&#770; = (</td><td>(</td><td>&beta;&#770;</td><td>&gamma;&#770;</td><td>) ∈ s<sub>e</sub>(3)</td><td>Tensor representation of &alpha; = (&beta;<sup>T</sup>&gamma;<sup>T</sup>)<sup>T</sup> with &beta;, &gamma; ∈ R<sup>3</sup>.</td></tr><tr><td>J<sub>f</sub>(x) =</td><td>(</td><td>&part;<sub>f</sub>/&part;x</td><td>) ∈ R<sup>l×h</sup></td><td>Jacobian of the vector function f(x) : R<sup>h</sup> → R<sup>l</sup>.</td></tr><tr><td>tanh(v)</td><td>Vector obtained by applying tanh(&middot;) componentwise to v.</td></tr></tbody></table>

## I. INTRODUCTION

ELECTRICAL, hydraulic, and mechanical systems, or their combinations, are Lagrangian systems that usually exhibit complex behavior. However, their physical nature displays special properties, such as symmetry and passivity, which have been
exploited to solve many control problems [1], [2], [3], [4], [5], otherwise difficult to address for generic nonlinear dynamics. To cope with their high nonlinearity and large number of degrees of freedom (DOF), representations with specific structures play a crucial role in simplifying analysis, as well as control design and synthesis. For example, coordinate transformations are often used to highlight some internal structures that simplify the derivation of feedback controllers for robotic systems and prove their stability [6], [7], [8], [9]. This article considers the input decoupling (ID) problem for input-affine Lagrangian systems. In particular, given a Lagrangian system whose inputs enter the second-order equations of motion through a configuration-dependent actuation matrix, we study under which conditions on the definition of generalized coordinates each input affects one and only one equation of motion. The coordinates solving the ID problem constitute a set of variables that simplify control design and synthesis because the input directly affects the equations of motion in a decentralized form. If the decoupling coordinates are regarded as system outputs, any control law implemented in these coordinates automatically decouples also the input-output channels. For fully actuated or overactuated Lagrangian systems, i.e., when the number of independent inputs is equal to or larger than the number of generalized coordinates, the ID problem can be solved by a configuration-dependent input transformation that inverts (or pseudoinverts) the actuation matrix. On the other hand, in the case of underactuated systems, a solution is available so far if the actuation matrix is constant and requires only a linear change of generalized coordinates [10].

For multi-input-multioutput plants described in state-space form, the input-output decoupling problem is solved by an inversion-based feedback controller, both in the linear [11] and nonlinear case [12], [13], [14]. For this, a feedback transformation is needed that requires a state and input transformation together with a feedback action. In particular, input-output decoupling is possible by a static state feedback if and only if the system has a well-defined vector relative degree, namely a nonsingular decoupling matrix [15], [16]. This result can be used to solve the ID problem for Lagrangian systems by taking the configuration variables (or a subset of them) as system outputs. In fact, in this way, the system has a vector relative degree, with all outputs having uniform relative degree two. It can be shown that, in the case of underactuation, the necessary and sufficient condition is the involutivity of the distribution spanned by the columns of the actuation matrix [17]. Unfortunately, involutivity is not easy to check, especially for high-dimensional dynamics, and has to be evaluated case by case, see [18], [19], [20], and [21]. Furthermore, finding the ID state variables in this case requires solving a system of nonlinear partial differential equations [16, Ch. 5], which could be impractical for control synthesis. Finally, in [22], Skogestad et al. have proposed input transformations for dynamic processes that achieve exact linearization and ID, under complete model knowledge. However, the analysis is restricted to systems with equal number of states, inputs, and outputs.

In this article, we focus on Lagrangian systems and we show that the choice of particular coordinates, called *actuation coordinates*, solves the ID problem without requiring a configuration-dependent transformation of the input nor a state

<figure><img src="image_3.png" /><figcaption>Fig. 1. Graphical representation of the proposed change of coordinates addressing the ID problem for a fully actuated Lagrangian system. In the *q* coordinates, the input *u* undergoes a nonlinear transformation through the actuation matrix *A*(*q*) when performing work on *q*, i.e., $\tau_q = A(q)u$. The proposed change of coordinates $\theta = h(q)$ bends the configuration space so that each component of *u* acts directly on one component of $\theta$, namely $\tau_\theta = u$. The existence of this transformation is possible because of the conservation of power $\dot{H}_q$, represented by the yellow area, under a change of coordinates, i.e., $\dot{H}_q = \dot{H}_\theta$.</figcaption></figure>

feedback. We derive necessary and sufficient conditions under which the actuation coordinates exist and constructively show that these coordinates solve the ID problem for fully actuated, overactuated, and underactuated Lagrangian systems. In particular, the unactuated coordinates remain arbitrary when the system is underactuated.

These results stem from power invariance under a change of coordinates, as graphically illustrated in Fig. 1 for the fully actuated case. Remarkably, similar considerations hold also when the dynamics is underactuated. We apply the results to several mechanical examples as archetypal Lagrangian systems. In addition, we prove that robotic systems driven by threadlike actuators, such as inelastic tendons or thin fluidic chambers, always admit actuation coordinates. Our results have relevant consequences on the control of continuum soft robots and other underactuated mechanical systems. Indeed, recent control laws for planar underactuated soft robots [10], [23], [24], [25], [26] generalize to all collocated mechanical systems with damping, such as soft robots moving in 3-D. Interestingly, the energy-based regulator in [27] uses in fact actuation coordinates to decouple the equations of motion of a 3-DOF underactuated soft robotic system. Similarly, in flexible link robots, the actuation coordinates are the clamped angles at the base of each beam, which have been extensively used in control design [28].

The rest of the article is organized as follows. In Section II, we introduce the notation and formalize the ID problem. Section III defines the class of Lagrangian systems, called *collocated*, for which the ID problem is solvable. Necessary and sufficient conditions for solving the ID problem are then derived for fully actuated or overactuated dynamics in Section IV and for underactuated systems in Section V. In Section VI, we prove that threadlike actuators yield collocated mechanical systems. Section VII extends two control strategies derived for underactuated mechanical systems with constant actuation matrix to
exploited to solve many control problems [1], [2], [3], [4], [5], otherwise difficult to address for generic nonlinear dynamics. To cope with their high nonlinearity and large number of degrees of freedom (DOF), representations with specific structures play a crucial role in simplifying analysis, as well as control design and synthesis. For example, coordinate transformations are often used to highlight some internal structures that simplify the derivation of feedback controllers for robotic systems and prove their stability [6], [7], [8], [9]. This article considers the input decoupling (ID) problem for input-affine Lagrangian systems. In particular, given a Lagrangian system whose inputs enter the second-order equations of motion through a configuration-dependent actuation matrix, we study under which conditions on the definition of generalized coordinates each input affects one and only one equation of motion. The coordinates solving the ID problem constitute a set of variables that simplify control design and synthesis because the input directly affects the equations of motion in a decentralized form. If the decoupling coordinates are regarded as system outputs, any control law implemented in these coordinates automatically decouples also the input-output channels. For fully actuated or overactuated Lagrangian systems, i.e., when the number of independent inputs is equal to or larger than the number of generalized coordinates, the ID problem can be solved by a configuration-dependent input transformation that inverts (or pseudoinverts) the actuation matrix. On the other hand, in the case of underactuated systems, a solution is available so far if the actuation matrix is constant and requires only a linear change of generalized coordinates [10].

For multi-input-multioutput plants described in state-space form, the input-output decoupling problem is solved by an inversion-based feedback controller, both in the linear [11] and nonlinear case [12], [13], [14]. For this, a feedback transformation is needed that requires a state and input transformation together with a feedback action. In particular, input-output decoupling is possible by a static state feedback if and only if the system has a well-defined vector relative degree, namely a nonsingular decoupling matrix [15], [16]. This result can be used to solve the ID problem for Lagrangian systems by taking the configuration variables (or a subset of them) as system outputs. In fact, in this way, the system has a vector relative degree, with all outputs having uniform relative degree two. It can be shown that, in the case of underactuation, the necessary and sufficient condition is the involutivity of the distribution spanned by the columns of the actuation matrix [17]. Unfortunately, involutivity is not easy to check, especially for high-dimensional dynamics, and has to be evaluated case by case, see [18], [19], [20], and [21]. Furthermore, finding the ID state variables in this case requires solving a system of nonlinear partial differential equations [16, Ch. 5], which could be impractical for control synthesis. Finally, in [22], Skogestad et al. have proposed input transformations for dynamic processes that achieve exact linearization and ID, under complete model knowledge. However, the analysis is restricted to systems with equal number of states, inputs, and outputs.

In this article, we focus on Lagrangian systems and we show that the choice of particular coordinates, called *actuation coordinates*, solves the ID problem without requiring a configuration-dependent transformation of the input nor a state

<figure><img src="image_3.png" /><figcaption>Fig. 1. Graphical representation of the proposed change of coordinates addressing the ID problem for a fully actuated Lagrangian system. In the *q* coordinates, the input *u* undergoes a nonlinear transformation through the actuation matrix *A*(*q*) when performing work on *q*, i.e., $\tau_q = A(q)u$. The proposed change of coordinates $\theta = h(q)$ bends the configuration space so that each component of *u* acts directly on one component of $\theta$, namely $\tau_\theta = u$. The existence of this transformation is possible because of the conservation of power $\dot{H}_q$, represented by the yellow area, under a change of coordinates, i.e., $\dot{H}_q = \dot{H}_\theta$.</figcaption></figure>

feedback. We derive necessary and sufficient conditions under which the actuation coordinates exist and constructively show that these coordinates solve the ID problem for fully actuated, overactuated, and underactuated Lagrangian systems. In particular, the unactuated coordinates remain arbitrary when the system is underactuated.

These results stem from power invariance under a change of coordinates, as graphically illustrated in Fig. 1 for the fully actuated case. Remarkably, similar considerations hold also when the dynamics is underactuated. We apply the results to several mechanical examples as archetypal Lagrangian systems. In addition, we prove that robotic systems driven by threadlike actuators, such as inelastic tendons or thin fluidic chambers, always admit actuation coordinates. Our results have relevant consequences on the control of continuum soft robots and other underactuated mechanical systems. Indeed, recent control laws for planar underactuated soft robots [10], [23], [24], [25], [26] generalize to all collocated mechanical systems with damping, such as soft robots moving in 3-D. Interestingly, the energy-based regulator in [27] uses in fact actuation coordinates to decouple the equations of motion of a 3-DOF underactuated soft robotic system. Similarly, in flexible link robots, the actuation coordinates are the clamped angles at the base of each beam, which have been extensively used in control design [28].

The rest of the article is organized as follows. In Section II, we introduce the notation and formalize the ID problem. Section III defines the class of Lagrangian systems, called *collocated*, for which the ID problem is solvable. Necessary and sufficient conditions for solving the ID problem are then derived for fully actuated or overactuated dynamics in Section IV and for underactuated systems in Section V. In Section VI, we prove that threadlike actuators yield collocated mechanical systems. Section VII extends two control strategies derived for underactuated mechanical systems with constant actuation matrix to
the collocated case, validating one of these controllers on a 3-D
tendon-driven underactuated soft robot. Finally, Section VIII
concludes this article.

II. PRELIMINARIES

A. Notation

We denote vectors and matrices with bold letters. Arguments
of the functions are omitted when clear from the context. The
notations adopted in the article are given in the Nomenclature.

B. *Dynamic Model*

Let $q \in \mathcal{M}$ be the generalized coordinates of a dynamical
system evolving on an $n$-dimensional smooth manifold $\mathcal{M}$ with
Lagrangian $\mathcal{L}_q(q, \dot{q})$. The system trajectories satisfy the Euler–
Lagrange equations of motions

$$
\left\{
\begin{array}{l}
\frac{d}{dt} \left( \frac{\partial \mathcal{L}_q(q, \dot{q})}{\partial \dot{q}} \right)^T - \left( \frac{\partial \mathcal{L}_q(q, \dot{q})}{\partial q} \right)^T = \tau_q(q, u) \\
\tau_q(q, u) = \mathbf{A}(q)u
\end{array}
\right.
\qquad (1)
$$

where $u \in \mathbb{R}^m$ are the available actuation inputs, $\mathbf{A}(q) \in \mathbb{R}^{n \times m}$ is the actuation matrix, and $\tau_q(q, u) \in \text{Im}(\mathbf{A}(q))$ collects the generalized forces performing work on $q$. For all $q \in \mathcal{M}$, we assume that $\mathbf{A}(q)$ is a full-rank matrix, i.e., $r = \text{rank}(\mathbf{A}(q)) = \min(m, n)$. When the dynamics is fully actuated ($m = n$) or underactuated ($m < n$), this is equivalent to asking that the actuation channels are all independent. On the other hand, if (1) is overactuated ($m > n$), we assume that there are exactly $n$ independent inputs. Note that, for the following derivations, when $r = m^* < \min(m, n)$ one can consistently discard $m - m^*$ linearly dependent columns of $\mathbf{A}(q)$ and consider the dynamics as underactuated.

In Appendix A, we recall two basic properties of Lagrangian
systems used in the following results.

C. Problem Statement

We look for a change of coordinates $\theta = h(q)$ from $\mathcal{B}(q) \subset \mathcal{M}$ to $\mathcal{N}$ where each of the first $r$ equations of motion in (1) is affected by one, and only one, independent actuator input, i.e., the right-hand side of the transformed equations of motion (see Property 2 in Appendix A) takes the form

$$
\tau_{\theta}(\theta, u) = \mathbf{A}_{\theta}(\theta)u = \begin{pmatrix} \mathbf{I}_r \\ \mathbf{O}_{n-r \times r} \end{pmatrix} u. \quad (2)
$$

We refer to such problem as the ID problem for the Lagrangian
dynamics (1). If a solution exists, then we say that (1) admits
a *collocated form*. Note that (2) covers only the fully- and
underactuated cases. When the dynamics is overactuated, i.e.,
*m* > *n*, it is impossible to obtain (2) because only *r* = *n* input
channels can be decoupled, and the remaining *m* − *n* inputs will
affect the dynamics through a configuration-dependent actuation
matrix.

III. COLLOCATED LAGRANGIAN SYSTEMS

In this section, we characterize a new class of Lagrangian
systems, which we call *collocated* because only such Lagrangian

dynamics admit a collocated form under a change of generalized
coordinates. In addition, a set of coordinates that solve the ID
problem comes for free without further system analysis.

To this end, we will exploit a concept that is known as the
passive output in the context of passivity-based control [29].
Consider the following vector function linear in the velocity:

$$
\dot{y} = \mathbf{A}^T(q) \dot{q}
$$

which is called the passive output because (1) is passive with
respect to the pair (u, dot y), with the storage function being the
system Hamiltonian (see Appendix A). We will assume that (1)
has *dot y* integrable as follows.

Integrability assumption: For all $q \in \mathcal{M}$, there exists a function $g(q) : \mathcal{M} \to \mathbb{R}^m$ such that

$$
\mathbf{J}_g(\mathbf{q}) = \frac{\partial g}{\partial \mathbf{q}} = \mathbf{A}^T(\mathbf{q}). \qquad (3)
$$

If the passive output is integrable, then we say that the
Lagrangian system (1) is collocated because it admits a col-
located form as defined in Section II-C. Furthermore, we define
$\mathbf{y} = \mathbf{g}(\mathbf{q})$ as actuation coordinates because, in such coordinates,
$\mathbf{u}$ acts directly on the equations of motion according to (2).

Remark: Each component of g(q) is defined up to a constant since any function $\bar{g}(q) = g(q) + k$, with $k \in \mathbb{R}^m$, satisfies the condition $\mathbf{J}_{\bar{g}} = \mathbf{J}_g$.

The integrability assumption requires each column of $\mathbf{A}(q)$
to be the gradient of a scalar function of the configuration vari-
ables. If $\mathbf{A}(q)$ is constant, i.e., $\mathbf{A}(q) = \mathbf{A}$, then $g(q) = \mathbf{A}^T q$.
More in general, when the column $\mathbf{A}_i(q)$ has continuous partial
derivatives, $\dot{y}_i = \mathbf{A}_i^T(q)\dot{q}$ is integrable [30, Ch. 2] if and only if

$$
\frac{\partial A_{ji}}{\partial q_k} = \frac{\partial A_{ki}}{\partial q_j} \quad \forall j, k \in \{1, \dots, n\}. \qquad (4)
$$

Note that this condition is equivalent to asking that, when $u_i$ is constant, the generalized work done by $u_i$ on $q$ does not depend on the system trajectories but only on the initial and final configurations $q_a$ and $q_b$, respectively, i.e.,

$$
W_{u_i}(q) := \int_{q_a}^{q_b} u_i \mathbf{A}_i^T(q) dq = u_i [g_i(q_b) - g_i(q_a)]
$$

where the last equality follows from the gradient theorem [30,
Prop. 1, Ch. 2]. In other words, $P_i(u_i, q) := u_i g_i(q)$ plays the
role of a potential energy for the dynamics. If the actuation
matrix is obtained using a differential formalism, such as the
virtual works principle, it is reasonable to expect—although
without any guarantee—that the integrability holds due to the
inherent differentiation involved. It is also worth observing the
following.

Remark: If $\mathbf{A}^T(q)$ is integrable according to (3), then the orthogonal complement to the codistribution spanned by $\mathbf{A}^T(q)$ satisfies the Frobenius theorem. However, in general, the inverse implication does not hold without also an input transformation.

Even when (4) is satisfied, it could be challenging to integrate
the passive output *ẏ* in the closed form. Nonetheless, it is always
possible to perform the numerical integration online based on the
measure of *q* and *dot q*. Assuming an exact knowledge of *A(q)* and
neglecting integration errors, there is formally no difference in
the collocated case, validating one of these controllers on a 3-D
tendon-driven underactuated soft robot. Finally, Section VIII
concludes this article.

II. PRELIMINARIES

A. Notation

We denote vectors and matrices with bold letters. Arguments
of the functions are omitted when clear from the context. The
notations adopted in the article are given in the Nomenclature.

B. Dynamic Model

Let $q \in \mathcal{M}$ be the generalized coordinates of a dynamical
system evolving on an $n$-dimensional smooth manifold $\mathcal{M}$ with
Lagrangian $\mathcal{L}_q(q, \dot{q})$. The system trajectories satisfy the Euler–
Lagrange equations of motions

$$
\left\{
\begin{array}{l}
\frac{d}{dt} \left( \frac{\partial \mathcal{L}_q(q, \dot{q})}{\partial \dot{q}} \right)^T - \left( \frac{\partial \mathcal{L}_q(q, \dot{q})}{\partial q} \right)^T = \tau_q(q, u) \\
\tau_q(q, u) = \mathbf{A}(q)u
\end{array}
\right.
\qquad (1)
$$

where $u \in \mathbb{R}^m$ are the available actuation inputs, $\mathbf{A}(q) \in \mathbb{R}^{n \times m}$ is the actuation matrix, and $\tau_q(q, u) \in \text{Im}(\mathbf{A}(q))$ collects the generalized forces performing work on $q$. For all $q \in \mathcal{M}$, we assume that $\mathbf{A}(q)$ is a full-rank matrix, i.e., $r = \text{rank}(\mathbf{A}(q)) = \min(m, n)$. When the dynamics is fully actuated ($m = n$) or underactuated ($m < n$), this is equivalent to asking that the actuation channels are all independent. On the other hand, if (1) is overactuated ($m > n$), we assume that there are exactly $n$ independent inputs. Note that, for the following derivations, when $r = m^* < \min(m, n)$ one can consistently discard $m - m^*$ linearly dependent columns of $\mathbf{A}(q)$ and consider the dynamics as underactuated.

In Appendix A, we recall two basic properties of Lagrangian
systems used in the following results.

C. Problem Statement

We look for a change of coordinates $\theta = h(q)$ from $\mathcal{B}(q) \subset \mathcal{M}$ to $\mathcal{N}$ where each of the first $r$ equations of motion in (1) is affected by one, and only one, independent actuator input, i.e., the right-hand side of the transformed equations of motion (see Property 2 in Appendix A) takes the form

$$
\tau_{\theta}(\theta, u) = \mathbf{A}_{\theta}(\theta)u = \begin{pmatrix} \mathbf{I}_r \\ \mathbf{O}_{n-r \times r} \end{pmatrix} u. \quad (2)
$$

We refer to such problem as the ID problem for the Lagrangian
dynamics (1). If a solution exists, then we say that (1) admits
a *collocated form*. Note that (2) covers only the fully- and
underactuated cases. When the dynamics is overactuated, i.e.,
*m* > *n*, it is impossible to obtain (2) because only *r* = *n* input
channels can be decoupled, and the remaining *m* − *n* inputs will
affect the dynamics through a configuration-dependent actuation
matrix.

III. COLLOCATED LAGRANGIAN SYSTEMS

In this section, we characterize a new class of Lagrangian
systems, which we call *collocated* because only such Lagrangian

dynamics admit a collocated form under a change of generalized
coordinates. In addition, a set of coordinates that solve the ID
problem comes for free without further system analysis.

To this end, we will exploit a concept that is known as the
passive output in the context of passivity-based control [29].
Consider the following vector function linear in the velocity:

$$
\dot{y} = A^T(q) \dot{q}
$$

which is called the passive output because (1) is passive with
respect to the pair (u, dot y), with the storage function being the
system Hamiltonian (see Appendix A). We will assume that (1)
has *dot y* integrable as follows.

Integrability assumption: For all $q \in \mathcal{M}$, there exists a function $g(q) : \mathcal{M} \to \mathbb{R}^m$ such that

$$
J_g(q) = \frac{\partial g}{\partial q} = A^T(q). \quad (3)
$$

If the passive output is integrable, then we say that the
Lagrangian system (1) is collocated because it admits a col-
located form as defined in Section II-C. Furthermore, we define
$\mathbf{y} = g(\mathbf{q})$ as actuation coordinates because, in such coordinates,
$\mathbf{u}$ acts directly on the equations of motion according to (2).

Remark: Each component of g(q) is defined up to a constant since any function $\bar{g}(q) = g(q) + k$, with $k \in \mathbb{R}^m$, satisfies the condition $J_{\bar{g}} = J_g$.

The integrability assumption requires each column of $\mathbf{A}(q)$ to be the gradient of a scalar function of the configuration variables. If $\mathbf{A}(q)$ is constant, i.e., $\mathbf{A}(q) = \mathbf{A}$, then $g(q) = \mathbf{A}^T q$. More in general, when the column $\mathbf{A}_i(q)$ has continuous partial derivatives, $\dot{y}_i = \mathbf{A}_i^T(q) \dot{q}$ is integrable [30, Ch. 2] if and only if

$$
\frac{\partial A_{ji}}{\partial q_k} = \frac{\partial A_{ki}}{\partial q_j} \quad \forall j, k \in \{1, \dots, n\}. \quad (4)
$$

Note that this condition is equivalent to asking that, when $u_i$ is constant, the generalized work done by $u_i$ on $q$ does not depend on the system trajectories but only on the initial and final configurations $q_a$ and $q_b$, respectively, i.e.,

$$
W_{u_i}(q) := \int_{q_a}^{q_b} u_i \mathbf{A}_i^T(q) dq = u_i [g_i(q_b) - g_i(q_a)]
$$

where the last equality follows from the gradient theorem [30, Prop. 1, Ch. 2]. In other words, $P_i(u_i, q) := u_i g_i(q)$ plays the role of a potential energy for the dynamics. If the actuation matrix is obtained using a differential formalism, such as the virtual works principle, it is reasonable to expect—although without any guarantee—that the integrability holds due to the inherent differentiation involved. It is also worth observing the following.

Remark: If $\mathbf{A}^T(q)$ is integrable according to (3), then the orthogonal complement to the codistribution spanned by $\mathbf{A}^T(q)$ satisfies the Frobenius theorem. However, in general, the inverse implication does not hold without also an input transformation.

Even when (4) is satisfied, it could be challenging to integrate
the passive output *ẏ* in the closed form. Nonetheless, it is always
possible to perform the numerical integration online based on the
measure of *q* and *dot q*. Assuming an exact knowledge of *A(q)* and
neglecting integration errors, there is formally no difference in
<figure><img src="image_4.png" /><figcaption>Fig. 2. Geostationary satellite actuated by a normal force $u_1$ and a tangential force $u_2$. The body configuration is described by the distance $q_1$ from the Earth center and the angle $q_2$ with respect to the horizontal axis. Only the normal force is collocated because it performs work directly and only on $q_1$.</figcaption></figure>

having $y$ in a closed form or computing it online. Furthermore, in many cases, the numerical integration should not be necessary. This is because the actuation coordinates are inherently related to the system inputs and should be easily measurable.

In some cases, the integrability may come directly as a consequence of the physical nature of the system, as illustrated in the following.

**Example 1 (Cartesian forces on a robot):** Inspired by the work in [31], consider a manipulator with n-DOF subject to m external forces applied to its structure. Assume that each force $f_i$ changes its magnitude over time but keeps the same direction in the global Cartesian frame so that

$$f_i = d_{f,i}u_i; \quad i \in \{1, \dots, m\}$$

with $\|d_{f,i}\| = 1$, and where $d_{f,i} \in \mathbb{R}^3$ represents the direction of $f_i$ and $u_i \in \mathbb{R}^{>0}$ its magnitude. If $p_i(q) \in \mathbb{R}^3$ denotes the point of application of $f_i$ in the global frame, then the effect of the force in the dynamics is

$$\tau_{q,i} = J_{p_i}^T(q) f_i = J_{p_i}^T(q) d_{f,i} u_i = A_i(q) u_i.$$

Thus, one can integrate $A_i(q)$ as $g_i(q) = d_{f,i}^T p_i(q)$. $\triangleleft$

On the other hand, the integrability conditions may not hold even for elementary dynamics.

**Example 2 (Geostationary satellite):** Consider a geostationary satellite orbiting around the Earth in a plane. The configuration $q = (q_1 q_2)^T$ describes its motion, being $q_1$ the distance from the Earth center and $q_2$ the angle with respect to a reference axis, as shown in Fig. 2. The satellite is controlled by a normal and tangential force, denoted as $u_1$ and $u_2$, respectively. The application of the Lagrangian formalism yields

$$\tau_q = A(q)u = \begin{pmatrix} 1 & 0 \\ 0 & q_1 \end{pmatrix} u.$$

The first column of $A(q)$ is already in the collocated form because $u_1$ performs work only on $q_1$. However, the second

column of $A(q)$ is nonintegrable because

$$\frac{\partial A_{12}}{\partial q_2} = 0 \neq 1 = \frac{\partial A_{22}}{\partial q_1}.$$

Indeed, the generalized work done on $q$ by the tangential force $u_2$ depends on the time evolution of $q_1$. $\triangleleft$

In the following, we consider the three different actuation scenarios, namely fully actuated, overactuated, and underactuated dynamics.

#### IV. CHANGE OF COORDINATES FOR FULLY ACTUATED SYSTEMS

In this section, we show that when the system is fully actuated, the existence of n actuator coordinates is necessary and sufficient to bring the dynamics to the collocated form. This case allows illustrating the results in the most simple scenario, where an input transformation is sufficient to solve the ID problem. However, the same arguments will be used also for underactuated systems.

**Theorem 1:** Suppose the system is fully actuated, i.e., $m = n$. There exists a change of coordinates $\theta = h(q) : \mathcal{B}(q) \to \mathcal{N}$ such that (1) takes the form

$$\frac{d}{dt} \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \dot{\theta}} \right)^T - \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \theta} \right)^T = u \quad (5)$$

i.e., $A_{\theta}(\theta) = I_n$, if and only if the integrability assumption (3) holds. Let $g(q)$ be the integral of $A^T(q)\dot{q}$. Then, a possible choice for $\theta$ is $\theta = g(q)$.

*Proof:* Under the integrability assumption, $\theta = g(q)$ defines a change of coordinates because its Jacobian $J_g(q) = A^T(q)$ has rank n at $q$.

Since the generalized power is coordinate invariant (Property 2 in Appendix A), it follows

$$\dot{\theta}^T \tau_{\theta} = \dot{q}^T \tau_{q}.$$

Noting that $\dot{\theta} = J_g(q)\dot{q}$ and using $\tau_q = A(q)u$, the previous equation rewrites as

$$\ddot{q}^T J_g^T(q) \tau_{\theta} = \dot{q}^T A(q) u$$

or, equivalently

$$\dot{q}^T (J_g^T(q) \tau_{\theta} - A(q)u) = 0. \quad (6)$$

Since (6) holds for all $\dot{q} \in T_q M$, it follows that

$$J_g^T(q) \tau_{\theta} - A(q)u = 0.$$

Furthermore, $J_g(q) = A^T(q)$ leads to

$$A(q) (\tau_{\theta} - u) = 0.$$

The previous equation defines a homogeneous linear system in the unknown $\tau_{\theta} - u$, which admits the unique solution $\tau_{\theta} = u$ since $A(q)$ is nonsingular, thus yielding the sufficiency of (5).

As for the necessity, suppose that a change of coordinates $\theta = h(q)$ exists such that (5) holds. Property 2 implies that, for all $\dot{q} \in T_q M$

$$\dot{q}^T (J_h^T(q) - A(q)) u = 0$$
<figure><img src="image_4.png" /><figcaption>Fig. 2. Geostationary satellite actuated by a normal force $u_1$ and a tangential force $u_2$. The body configuration is described by the distance $q_1$ from the Earth center and the angle $q_2$ with respect to the horizontal axis. Only the normal force is collocated because it performs work directly and only on $q_1$.</figcaption></figure>

having $y$ in a closed form or computing it online. Furthermore, in many cases, the numerical integration should not be necessary. This is because the actuation coordinates are inherently related to the system inputs and should be easily measurable.

In some cases, the integrability may come directly as a consequence of the physical nature of the system, as illustrated in the following.

**Example 1 (Cartesian forces on a robot):** Inspired by the work in [31], consider a manipulator with n-DOF subject to m external forces applied to its structure. Assume that each force $f_i$ changes its magnitude over time but keeps the same direction in the global Cartesian frame so that

$$f_i = d_{f,i}u_i; \quad i \in \{1, \dots, m\}$$

with $\|d_{f,i}\| = 1$, and where $d_{f,i} \in \mathbb{R}^3$ represents the direction of $f_i$ and $u_i \in \mathbb{R}^{>0}$ its magnitude. If $p_i(q) \in \mathbb{R}^3$ denotes the point of application of $f_i$ in the global frame, then the effect of the force in the dynamics is

$$\tau_{q,i} = J_{p_i}^T(q) f_i = J_{p_i}^T(q) d_{f,i} u_i = A_i(q) u_i.$$

Thus, one can integrate $A_i(q)$ as $g_i(q) = d_{f,i}^T p_i(q)$. $\triangleleft$

On the other hand, the integrability conditions may not hold even for elementary dynamics.

**Example 2 (Geostationary satellite):** Consider a geostationary satellite orbiting around the Earth in a plane. The configuration $q = (q_1 q_2)^T$ describes its motion, being $q_1$ the distance from the Earth center and $q_2$ the angle with respect to a reference axis, as shown in Fig. 2. The satellite is controlled by a normal and tangential force, denoted as $u_1$ and $u_2$, respectively. The application of the Lagrangian formalism yields

$$\tau_q = A(q)u = \begin{pmatrix} 1 & 0 \\ 0 & q_1 \end{pmatrix} u.$$

The first column of $A(q)$ is already in the collocated form because $u_1$ performs work only on $q_1$. However, the second

column of $A(q)$ is nonintegrable because

$$\frac{\partial A_{12}}{\partial q_2} = 0 \neq 1 = \frac{\partial A_{22}}{\partial q_1}.$$

Indeed, the generalized work done on $q$ by the tangential force $u_2$ depends on the time evolution of $q_1$. $\triangleleft$

In the following, we consider the three different actuation scenarios, namely fully actuated, overactuated, and underactuated dynamics.

#### IV. CHANGE OF COORDINATES FOR FULLY ACTUATED SYSTEMS

In this section, we show that when the system is fully actuated, the existence of $n$ actuator coordinates is necessary and sufficient to bring the dynamics to the collocated form. This case allows illustrating the results in the most simple scenario, where an input transformation is sufficient to solve the ID problem. However, the same arguments will be used also for underactuated systems.

**Theorem 1:** Suppose the system is fully actuated, i.e., $m = n$. There exists a change of coordinates $\theta = h(q) : \mathcal{B}(q) \to \mathcal{N}$ such that (1) takes the form

$$\frac{d}{dt} \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \dot{\theta}} \right)^T - \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \theta} \right)^T = u \quad (5)$$

i.e., $A_{\theta}(\theta) = I_n$, if and only if the integrability assumption (3) holds. Let $g(q)$ be the integral of $A^T(q)\dot{q}$. Then, a possible choice for $\theta$ is $\theta = g(q)$.

*Proof:* Under the integrability assumption, $\theta = g(q)$ defines a change of coordinates because its Jacobian $J_g(q) = A^T(q)$ has rank $n$ at $q$.

Since the generalized power is coordinate invariant (Property 2 in Appendix A), it follows

$$\dot{\theta}^T \tau_{\theta} = \dot{q}^T \tau_q.$$

Noting that $\dot{\theta} = J_g(q)\dot{q}$ and using $\tau_q = A(q)u$, the previous equation rewrites as

$$\ddot{q}^T J_g^T(q) \tau_{\theta} = \dot{q}^T A(q) u$$

or, equivalently

$$\dot{q}^T (J_g^T(q) \tau_{\theta} - A(q)u) = 0. \quad (6)$$

Since (6) holds for all $\dot{q} \in T_q M$, it follows that

$$J_g^T(q) \tau_{\theta} - A(q)u = 0.$$

Furthermore, $J_g(q) = A^T(q)$ leads to

$$A(q) (\tau_{\theta} - u) = 0.$$

The previous equation defines a homogeneous linear system in the unknown $\tau_{\theta} - u$, which admits the unique solution $\tau_{\theta} = u$ since $A(q)$ is nonsingular, thus yielding the sufficiency of (5).

As for the necessity, suppose that a change of coordinates $\theta = h(q)$ exists such that (5) holds. Property 2 implies that, for all $\dot{q} \in T_q M$

$$\dot{q}^T (J_h^T(q) - A(q)) u = 0$$
<figure><img src="image_5.png" /><figcaption>Fig. 3. Planar mechanism with two rotational joints having configuration $q \in \mathbb{S}^1 \times \mathbb{S}^1$. Positive (relative) rotations are counted counterclockwise. For $i=1,2$, the force $u_i$ actuates a cart coupled to the robot by a linear spring of stiffness $k_i$. The dynamics of the carts is negligible, and thus, the forces $u$ act instantaneously on the mechanism.</figcaption></figure>

leading to

$$ (\mathbf{J}_h^T(\mathbf{q}) - \mathbf{A}(\mathbf{q})) \mathbf{u} = \mathbf{0} \quad \forall \mathbf{u} \in \mathbb{R}^n. $$

If one chooses $\mathbf{u} = (\mathbf{I}_n)_i; i \in \{1, \dots, n\}$, then

$$ (\mathbf{J}_h^T(\mathbf{q}) - \mathbf{A}(\mathbf{q})) \mathbf{u} = (\mathbf{J}_h^T(\mathbf{q}) - \mathbf{A}(\mathbf{q}))_i = \mathbf{0}. $$

Thus, it holds $\mathbf{J}_h^T(\mathbf{q}) = \mathbf{A}(\mathbf{q})$ and $\mathbf{g}(\mathbf{q}) = \mathbf{h}(\mathbf{q})$. □

The following example illustrates the previous result.

**Example 3 (Spring actuated mechanism):** Consider a planar mechanism with two passive revolute joints, having angles $q_1$ and $q_2$ so that $q = (q_1 \ q_2)^T$. A spring with stiffness $k_i$ is attached to the distal end of each link, whose length is $l_i$, $i = 1, 2$. The springs are also connected to two carts moving on linear rails under the forces $u_1$ and $u_2$, with reference to Fig. 3. Assuming that the dynamics of the carts is negligible, it can be shown that the actuator inputs $\mathbf{u}$ directly affect the joint motion through the actuation matrix

$$ \mathbf{A}(\mathbf{q}) = \begin{pmatrix} -l_1 s_1 & l_1 c_1 + l_2 c_{12} \\ 0 & l_2 c_{12} \end{pmatrix} $$

where $s_i(c_i) = \sin(q_i)(\cos(q_i))$ and $s_{ij}(c_{ij}) = \sin(q_i + q_j)(\cos(q_i + q_j))$. Outside of singularities $q_1 \in \{0, \pi\}$ and $q_1 + q_2 \in \{\pm\pi/2\}$, $\mathrm{rank}(\mathbf{A}(\mathbf{q})) = 2$ and the integrability test (4) is successful because

$$ \frac{\partial A_{11}}{\partial q_2} = \frac{\partial (-l_1 s_1)}{\partial q_2} = 0 = \frac{\partial A_{21}}{\partial q_1} $$

$$ \frac{\partial A_{21}}{\partial q_2} = \frac{\partial (l_1 c_1 + l_2 c_{12})}{\partial q_2} = -l_2 s_{12} \\ = \frac{\partial (l_2 c_{12})}{\partial q_1} = \frac{\partial A_{22}}{\partial q_1}. $$

The passive output is integrable as

$$ \mathbf{y} = \begin{pmatrix} l_1 c_1 \\ l_1 s_1 + l_2 s_{12} \end{pmatrix}. $$

Note that $y_1$ is the x-coordinate of the position of the spring end attached to link 1. Similarly, $y_2$ is the y-coordinate of the spring attached to link 2. Indeed, the forces $u_1$ and $u_2$ perform work on the distal ends of the spring attached to the mechanism along these directions. ▷

### A. Overactuated Case

The previous result extends to overactuated systems, namely dynamics with more inputs than generalized coordinates. Thus, we have $r = n < m$. We partition $\mathbf{A}(\mathbf{q})$, which is a wide matrix, as

$$ \mathbf{A}(\mathbf{q}) = (\mathbf{A}_a(\mathbf{q}) \quad \mathbf{A}_o(\mathbf{q})) \quad (7) $$

where $\mathbf{A}_a(\mathbf{q}) \in \mathbb{R}^{n \times n}$ and $\mathbf{A}_o(\mathbf{q}) \in \mathbb{R}^{n \times (m-n)}$. Without loss of generality, we can have that $\mathrm{rank}(\mathbf{A}_a(\mathbf{q})) = n$ and the integrability condition holds for $\mathbf{A}_a(\mathbf{q})$.

**Corollary 1:** If the system is overactuated, i.e., $r = n < m$ and the same hypotheses of Theorem 1 hold for $\mathbf{A}_a(\mathbf{q})$, then there exists a change of coordinates $\theta = \mathbf{h}(\mathbf{q}) : \mathcal{B}(\mathbf{q}) \to \mathcal{N}$ such that (1) takes the form

$$ \frac{\mathrm{d}}{\mathrm{d}t} \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \dot{\theta}} \right)^{\mathrm{T}} - \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \theta} \right)^{\mathrm{T}} = (\mathbf{I}_n \quad \mathbf{A}_{o,\theta}(\theta)) \mathbf{u} \quad (8) $$

where

$$ \mathbf{A}_{o,\theta}(\theta) = \mathbf{A}_a^{-1} \mathbf{A}_o(\mathbf{q} = \mathbf{h}^{-1}(\theta)) \in \mathbb{R}^{n \times (m-n)}. $$

If $\mathbf{g}(\mathbf{q})$ is the integral of $\mathbf{A}_a^T(\mathbf{q})\dot{\mathbf{q}}$, then $\theta$ can be chosen as $\theta = \mathbf{g}(\mathbf{q})$.

**Proof:** Choosing again $\theta = \mathbf{g}(\mathbf{q})$ and following steps similar to those of the proof of Theorem 1, one obtains

$$ \mathbf{A}_a(\mathbf{q}) \tau_{\theta} = \mathbf{A}(\mathbf{q}) \mathbf{u}. $$

Expanding $\mathbf{A}(\mathbf{q})$ into (7) and left-multiplying the previous equation by $\mathbf{A}_a^{-1}(\mathbf{q})$ gives

$$ \tau_{\theta} = \mathbf{A}_a^{-1}(\mathbf{q}) (\mathbf{A}_a(\mathbf{q}) \quad \mathbf{A}_o(\mathbf{q})) \mathbf{u} \\ = \left( \mathbf{I}_n \quad \frac{\mathbf{A}_a^{-1}(\mathbf{q}) \mathbf{A}_o(\mathbf{q})}{\mathbf{A}_{o,\theta}(\theta)} \right) \mathbf{u} $$

with $\mathbf{q} = \mathbf{g}^{-1}(\theta)$.

Now, assume a change of coordinates $\theta = \mathbf{h}(\mathbf{q})$ exists such that (8) holds. After some computations, power invariance leads to the algebraic system

$$ (\mathbf{J}_h^T(\mathbf{q})(\mathbf{q}) - \mathbf{A}_a(\mathbf{q}) \quad \mathbf{J}_h^T(\mathbf{q})(\mathbf{A}_{o,\theta}(\mathbf{q}) - \mathbf{A}_o(\mathbf{q})) \mathbf{u} = \mathbf{0} $$

which must hold for all $\mathbf{u} \in \mathbb{R}^m$. By taking $\mathbf{u}_i = (\mathbf{I}_m)_i; i \in \{1, \dots, m\}$, it follows $\mathbf{J}_h^T(\mathbf{q}) = \mathbf{A}_a(\mathbf{q})$ and $\mathbf{A}_o(\mathbf{q}) = \mathbf{J}_h^T(\mathbf{q}) \mathbf{A}_{o,\theta}(\mathbf{q})$. Hence, at least $n$ passive outputs are integrable as $\mathbf{y} = \mathbf{h}(\mathbf{q})$. □

Note that it is not possible, in general, to simplify the expression of both terms in the actuation matrix because there are too many input variables to be decoupled.

**Example 4 (Tendon-driven joint):** Consider the tendon-driven finger of [32, Ch. 6.4] with 1-DOF $q$ and two actuator inputs
<figure><img src="image_5.png" /><figcaption>Fig. 3. Planar mechanism with two rotational joints having configuration $q \in \mathbb{S}^1 \times \mathbb{S}^1$. Positive (relative) rotations are counted counterclockwise. For $i = 1, 2$, the force $u_i$ actuates a cart coupled to the robot by a linear spring of stiffness $k_i$. The dynamics of the carts is negligible, and thus, the forces $u$ act instantaneously on the mechanism.</figcaption></figure>

leading to

$$ (\mathbf{J}_h^T(\mathbf{q}) - \mathbf{A}(\mathbf{q})) \mathbf{u} = \mathbf{0} \quad \forall \mathbf{u} \in \mathbb{R}^n. $$

If one chooses $\mathbf{u} = (\mathbf{I}_n)_i; i \in \{1, \dots, n\}$, then

$$ (\mathbf{J}_h^T(\mathbf{q}) - \mathbf{A}(\mathbf{q})) \mathbf{u} = (\mathbf{J}_h^T(\mathbf{q}) - \mathbf{A}(\mathbf{q}))_i = \mathbf{0}. $$

Thus, it holds $\mathbf{J}_h^T(\mathbf{q}) = \mathbf{A}(\mathbf{q})$ and $\mathbf{g}(\mathbf{q}) = \mathbf{h}(\mathbf{q})$. □

The following example illustrates the previous result.

**Example 3 (Spring actuated mechanism):** Consider a planar mechanism with two passive revolute joints, having angles $q_1$ and $q_2$ so that $q = (q_1 \ q_2)^T$. A spring with stiffness $k_i$ is attached to the distal end of each link, whose length is $l_i$, $i = 1, 2$. The springs are also connected to two carts moving on linear rails under the forces $u_1$ and $u_2$, with reference to Fig. 3. Assuming that the dynamics of the carts is negligible, it can be shown that the actuator inputs $\mathbf{u}$ directly affect the joint motion through the actuation matrix

$$ \mathbf{A}(\mathbf{q}) = \begin{pmatrix} -l_1 s_1 & l_1 c_1 + l_2 c_{12} \\ 0 & l_2 c_{12} \end{pmatrix} $$

where $s_i(c_i) = \sin(q_i)(\cos(q_i))$ and $s_{ij}(c_{ij}) = \sin(q_i + q_j)(\cos(q_i + q_j))$. Outside of singularities $q_1 \in \{0, \pi\}$ and $q_1 + q_2 \in \{\pm\pi/2\}$, $\mathrm{rank}(\mathbf{A}(\mathbf{q})) = 2$ and the integrability test (4) is successful because

$$ \frac{\partial A_{11}}{\partial q_2} = \frac{\partial (-l_1 s_1)}{\partial q_2} = 0 = \frac{\partial A_{21}}{\partial q_1} $$

$$ \frac{\partial A_{21}}{\partial q_2} = \frac{\partial (l_1 c_1 + l_2 c_{12})}{\partial q_2} = -l_2 s_{12} \\ = \frac{\partial (l_2 c_{12})}{\partial q_1} = \frac{\partial A_{22}}{\partial q_1}. $$

The passive output is integrable as

$$ \mathbf{y} = \begin{pmatrix} l_1 c_1 \\ l_1 s_1 + l_2 s_{12} \end{pmatrix}. $$

Note that $y_1$ is the x-coordinate of the position of the spring end attached to link 1. Similarly, $y_2$ is the y-coordinate of the spring attached to link 2. Indeed, the forces $u_1$ and $u_2$ perform work on the distal ends of the spring attached to the mechanism along these directions. ▷

### A. Overactuated Case

The previous result extends to overactuated systems, namely dynamics with more inputs than generalized coordinates. Thus, we have $r = n < m$. We partition $\mathbf{A}(\mathbf{q})$, which is a wide matrix, as

$$ \mathbf{A}(\mathbf{q}) = (\mathbf{A}_a(\mathbf{q}) \quad \mathbf{A}_o(\mathbf{q})) \quad (7) $$

where $\mathbf{A}_a(\mathbf{q}) \in \mathbb{R}^{n \times n}$ and $\mathbf{A}_o(\mathbf{q}) \in \mathbb{R}^{n \times (m-n)}$. Without loss of generality, we can have that $\mathrm{rank}(\mathbf{A}_a(\mathbf{q})) = n$ and the integrability condition holds for $\mathbf{A}_a(\mathbf{q})$.

**Corollary 1:** If the system is overactuated, i.e., $r = n < m$ and the same hypotheses of Theorem 1 hold for $\mathbf{A}_a(\mathbf{q})$, then there exists a change of coordinates $\theta = \mathbf{h}(\mathbf{q}) : \mathcal{B}(\mathbf{q}) \to \mathcal{N}$ such that (1) takes the form

$$ \frac{\mathrm{d}}{\mathrm{d}t} \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \dot{\theta}} \right)^T - \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \theta} \right)^T = (\mathbf{I}_n \quad \mathbf{A}_{o,\theta}(\theta)) \mathbf{u} \quad (8) $$

where

$$ \mathbf{A}_{o,\theta}(\theta) = \mathbf{A}_a^{-1} \mathbf{A}_o(\mathbf{q} = \mathbf{h}^{-1}(\theta)) \in \mathbb{R}^{n \times (m-n)}. $$

If $\mathbf{g}(\mathbf{q})$ is the integral of $\mathbf{A}_a^T(\mathbf{q})\dot{\mathbf{q}}$, then $\theta$ can be chosen as $\theta = \mathbf{g}(\mathbf{q})$.

**Proof:** Choosing again $\theta = \mathbf{g}(\mathbf{q})$ and following steps similar to those of the proof of Theorem 1, one obtains

$$ \mathbf{A}_a(\mathbf{q}) \tau_{\theta} = \mathbf{A}(\mathbf{q}) \mathbf{u}. $$

Expanding $\mathbf{A}(\mathbf{q})$ into (7) and left-multiplying the previous equation by $\mathbf{A}_a^{-1}(\mathbf{q})$ gives

$$ \begin{aligned} \tau_{\theta} &= \mathbf{A}_a^{-1}(\mathbf{q}) (\mathbf{A}_a(\mathbf{q}) \quad \mathbf{A}_o(\mathbf{q})) \mathbf{u} \\ &= \left( \mathbf{I}_n \quad \frac{\mathbf{A}_a^{-1}(\mathbf{q}) \mathbf{A}_o(\mathbf{q})}{\mathbf{A}_{o,\theta}(\theta)} \right) \mathbf{u} \end{aligned} $$

with $\mathbf{q} = \mathbf{g}^{-1}(\theta)$.

Now, assume a change of coordinates $\theta = \mathbf{h}(\mathbf{q})$ exists such that (8) holds. After some computations, power invariance leads to the algebraic system

$$ (\mathbf{J}_h^T(\mathbf{q})(\mathbf{q}) - \mathbf{A}_a(\mathbf{q}) \quad \mathbf{J}_h^T(\mathbf{q})(\mathbf{A}_{o,\theta}(\mathbf{q}) - \mathbf{A}_o(\mathbf{q})) \quad \mathbf{u} = \mathbf{0} $$

which must hold for all $\mathbf{u} \in \mathbb{R}^m$. By taking $\mathbf{u}_i = (\mathbf{I}_m)_i; i \in \{1, \dots, m\}$, it follows $\mathbf{J}_h^T(\mathbf{q}) = \mathbf{A}_a(\mathbf{q})$ and $\mathbf{A}_o(\mathbf{q}) = \mathbf{J}_h^T(\mathbf{q}) \mathbf{A}_{o,\theta}(\mathbf{q})$. Hence, at least $n$ passive outputs are integrable as $\mathbf{y} = \mathbf{h}(\mathbf{q})$. □

Note that it is not possible, in general, to simplify the expression of both terms in the actuation matrix because there are too many input variables to be decoupled.

**Example 4 (Tendon-driven joint):** Consider the tendon-driven finger of [32, Ch. 6.4] with 1-DOF $q$ and two actuator inputs
<figure><img src="image_6.png" /><figcaption>Fig. 4. Tendon-driven rotational joint having a single configuration variable $q \in S^1$. The cables tension $u_1$ and $u_2$ generate a torque at the joint.</figcaption></figure>

$u = (u_1 \ u_2)^T$, as sketched in Fig. 4. Assume for simplicity that the angle $q > 0$ (similar results hold for $q < 0$). The system input is $\tau_q = A(q)u$, where

$$A(q) = \left( \sqrt{a^2 + b^2} \sin \left( \tan^{-1} \left( \frac{a}{b} \right) + \frac{q}{2} \right) - R \right)$$

and $u$ collects the cable tensions. Thus, the passive output is

$$\dot{y} = A^T(q)\dot{q} = \begin{pmatrix} \sqrt{a^2 + b^2} \sin \left( \tan^{-1} \left( \frac{a}{b} \right) + \frac{q}{2} \right) \dot{q} \\ -R\dot{q} \end{pmatrix}.$$

Since the system is overactuated, there are two possible choices of the actuation coordinates out of singularities of $A(q)$. When $q = -2 \tan^{-1}(a/b)$, the first column of $A(q)$ becomes zero, and the first actuator does not have any effect on the motion. In this case, one can easily integrate $\dot{y}$ as

$$y = g(q) = \begin{pmatrix} -2\sqrt{a^2 + b^2} \cos \left( \tan^{-1} \left( \frac{a}{b} \right) + \frac{q}{2} \right) \\ -Rq \end{pmatrix}.$$

As expected, the choice of either the first or second component of $y$ yields dynamics in the form of (8). Note that $y$ is the tendon displacement with respect to the straight configuration. In fact, the system inputs perform work directly on the finger tendons. $\triangle$

## V. CHANGE OF COORDINATES FOR UNDERACTUATED SYSTEMS

We now focus on underactuated systems, specifically dynamics for which the number of independent actuator inputs is less than that of DOF. In the following, we show that also in this case, the actuation coordinates solve the ID problem. First, it is convenient to expand the actuation matrix as

$$A(q) = \begin{pmatrix} A_a(q) \\ A_u(q) \end{pmatrix} \qquad (9)$$

where $A_a(q) \in \mathbb{R}^{m \times m}$ is nonsingular and $A_u(q) \in \mathbb{R}^{(n-m) \times m}$. Note that such partition is always possible after reordering the linearly independent rows of $A(q)$. Furthermore, since $n > r = m$, $\dot{y}$ is an $m$-dimensional vector.

**Theorem 2:** Assume that (1) is underactuated, i.e., $m < n$. The integrability assumption is a necessary and sufficient condition for a change of coordinates $\theta = h(q) : \mathcal{B}(q) \to \mathcal{N}$ to exist such that (1) becomes

$$\frac{d}{dt} \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \dot{\theta}} \right)^T - \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \theta} \right)^T = \begin{pmatrix} u \\ \mathbf{0}_{n-m} \end{pmatrix} \quad (10)$$

i.e., $A_{\theta}(\theta)$ has the form

$$A_{\theta}(\theta) = \begin{pmatrix} I_n \\ O_{(n-m) \times m} \end{pmatrix}.$$

Let $g(q)$ be the integral of $A^T(q)\dot{q}$. Then, a possible choice of $\theta$ is

$$\theta = \begin{pmatrix} g(q) \\ \mathbf{0}_{n-m} \end{pmatrix} + \begin{pmatrix} O_{m \times m} & O_{m \times (n-m)} \\ O_{(n-m) \times m} & \mathbf{I}_{n-m} \end{pmatrix} \cdot q. \quad (11)$$

*Proof:* The choice of $\theta$ as given in (11) qualifies as a change of coordinates because its Jacobian

$$J_h(q) = \begin{pmatrix} A_a^T(q) & A_u^T(q) \\ O_{(n-m) \times m} & \mathbf{I}_{n-m} \end{pmatrix} \qquad (12)$$

is nonsingular at $q$. Power invariance and (12) imply

$$\dot{\theta}^T \tau_{\theta} = \dot{q}^T J_h^T(q) \tau_{\theta} = \dot{q}^T A(q) u = \dot{q}^T \tau_{q}.$$

Furthermore, being $\dot{q}$ arbitrary, it follows $J_h^T(q)\tau_{\theta} = A(q)u$, which can be rewritten as

$$\underbrace{\begin{pmatrix} A_a(q) & O_{m \times (n-m)} \\ A_u(q) & \mathbf{I}_{n-m} \end{pmatrix}}_{J_h^T(q)} \begin{pmatrix} \tau_{\theta_a} - u \\ \tau_{\theta_u} \end{pmatrix} = 0 \qquad (13)$$

where we expanded $\tau_{\theta}$ into the two vectors $\tau_{\theta_a} \in \mathbb{R}^m$ and $\tau_{\theta_u} \in \mathbb{R}^{n-m}$ performing work on $\theta_a$ and $\theta_u$, respectively. Equation (13) describes a homogeneous linear system with the unique solution

$$\begin{pmatrix} \tau_{\theta_a} - u \\ \tau_{\theta_u} \end{pmatrix} = 0 \qquad (14)$$

being that $J_h^T(q)$ is nonsingular.

To prove the necessary part of the statement, suppose there exists $\theta = h(q)$ such that (10) holds and partition

$$J_h(q) = \begin{pmatrix} J_{h_a}(q) \\ J_{h_u}(q) \end{pmatrix}$$

with $J_{h_a}(q) \in \mathbb{R}^{m \times n}$ and $J_{h_u}(q) \in \mathbb{R}^{(n-m) \times n}$. Exploiting once again power invariance, we obtain, after some computations

$$\dot{q}^T (J_{h_a}^T(q) - A(q)) u = 0$$

or, equivalently

$$J_{h_a}(q) = A^T(q).$$

Thus, the first $m$ components of $h(q)$ satisfy the integrability assumption. $\square$

In Appendix B, we report an alternative proof of the sufficient part of Theorem 2, which uses algebraic arguments instead of power invariance.

**Remark:** There is no constraint on choosing the unactuated variables, except that the corresponding Jacobian is nonsingular. Indeed, the factorization given in (13) holds independently of $\theta_u$. In other words, Theorem 2 does not rely on a specific choice of $\theta_u$, which could be used to further simplify the structure of the equations of motion.
<figure><img src="image_6.png" /><figcaption>Fig. 4. Tendon-driven rotational joint having a single configuration variable $q \in S^1$. The cables tension $u_1$ and $u_2$ generate a torque at the joint.</figcaption></figure>

$u = (u_1 \ u_2)^T$, as sketched in Fig. 4. Assume for simplicity that the angle $q > 0$ (similar results hold for $q < 0$). The system input is $\tau_q = A(q)u$, where

$$A(q) = \left( \sqrt{a^2 + b^2} \sin \left( \tan^{-1} \left( \frac{a}{b} \right) + \frac{q}{2} \right) - R \right)$$

and $u$ collects the cable tensions. Thus, the passive output is

$$\dot{y} = A^T(q)\dot{q} = \begin{pmatrix} \sqrt{a^2 + b^2} \sin \left( \tan^{-1} \left( \frac{a}{b} \right) + \frac{q}{2} \right) \dot{q} \\ -R\dot{q} \end{pmatrix}.$$

Since the system is overactuated, there are two possible choices of the actuation coordinates out of singularities of $A(q)$. When $q = -2 \tan^{-1}(a/b)$, the first column of $A(q)$ becomes zero, and the first actuator does not have any effect on the motion. In this case, one can easily integrate $\dot{y}$ as

$$y = g(q) = \begin{pmatrix} -2\sqrt{a^2 + b^2} \cos \left( \tan^{-1} \left( \frac{a}{b} \right) + \frac{q}{2} \right) \\ -Rq \end{pmatrix}.$$

As expected, the choice of either the first or second component of $y$ yields dynamics in the form of (8). Note that $y$ is the tendon displacement with respect to the straight configuration. In fact, the system inputs perform work directly on the finger tendons. $\triangle$

## V. CHANGE OF COORDINATES FOR UNDERACTUATED SYSTEMS

We now focus on underactuated systems, specifically dynamics for which the number of independent actuator inputs is less than that of DOF. In the following, we show that also in this case, the actuation coordinates solve the ID problem. First, it is convenient to expand the actuation matrix as

$$A(q) = \begin{pmatrix} A_a(q) \\ A_u(q) \end{pmatrix} \qquad (9)$$

where $A_a(q) \in \mathbb{R}^{m \times m}$ is nonsingular and $A_u(q) \in \mathbb{R}^{(n-m) \times m}$. Note that such partition is always possible after reordering the linearly independent rows of $A(q)$. Furthermore, since $n > r = m$, $\dot{y}$ is an $m$-dimensional vector.

**Theorem 2:** Assume that (1) is underactuated, i.e., $m < n$. The integrability assumption is a necessary and sufficient condition for a change of coordinates $\theta = h(q) : \mathcal{B}(q) \to \mathcal{N}$ to exist such that (1) becomes

$$\frac{d}{dt} \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \dot{\theta}} \right)^T - \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \theta} \right)^T = \begin{pmatrix} u \\ \mathbf{0}_{n-m} \end{pmatrix} \quad (10)$$

i.e., $A_\theta(\theta)$ has the form

$$A_{\theta}(\theta) = \begin{pmatrix} I_n \\ O_{(n-m) \times m} \end{pmatrix}.$$

Let $g(q)$ be the integral of $A^T(q)\dot{q}$. Then, a possible choice of $\theta$ is

$$\theta = \begin{pmatrix} g(q) \\ \mathbf{0}_{n-m} \end{pmatrix} + \begin{pmatrix} O_{m \times m} & O_{m \times (n-m)} \\ O_{(n-m) \times m} & \mathbf{I}_{n-m} \end{pmatrix} \cdot q. \quad (11)$$

*Proof:* The choice of $\theta$ as given in (11) qualifies as a change of coordinates because its Jacobian

$$J_h(q) = \begin{pmatrix} A_a^T(q) & A_u^T(q) \\ O_{(n-m) \times m} & \mathbf{I}_{n-m} \end{pmatrix} \qquad (12)$$

is nonsingular at $q$. Power invariance and (12) imply

$$\dot{\theta}^T \tau_{\theta} = \dot{q}^T J_h^T(q) \tau_{\theta} = \dot{q}^T A(q) u = \dot{q}^T \tau_q.$$

Furthermore, being $\dot{q}$ arbitrary, it follows $J_h^T(q)\tau_{\theta} = A(q)u$, which can be rewritten as

$$\underbrace{\begin{pmatrix} A_a(q) & O_{m \times (n-m)} \\ A_u(q) & \mathbf{I}_{n-m} \end{pmatrix}}_{J_h^T(q)} \begin{pmatrix} \tau_{\theta_a} - u \\ \tau_{\theta_u} \end{pmatrix} = 0 \qquad (13)$$

where we expanded $\tau_{\theta}$ into the two vectors $\tau_{\theta_a} \in \mathbb{R}^m$ and $\tau_{\theta_u} \in \mathbb{R}^{n-m}$ performing work on $\theta_a$ and $\theta_u$, respectively. Equation (13) describes a homogeneous linear system with the unique solution

$$\begin{pmatrix} \tau_{\theta_a} - u \\ \tau_{\theta_u} \end{pmatrix} = 0 \qquad (14)$$

being that $J_h^T(q)$ is nonsingular.

To prove the necessary part of the statement, suppose there exists $\theta = h(q)$ such that (10) holds and partition

$$J_h(q) = \begin{pmatrix} J_{h_a}(q) \\ J_{h_u}(q) \end{pmatrix}$$

with $J_{h_a}(q) \in \mathbb{R}^{m \times n}$ and $J_{h_u}(q) \in \mathbb{R}^{(n-m) \times n}$. Exploiting once again power invariance, we obtain, after some computations

$$\dot{q}^T (J_{h_a}^T(q) - A(q)) u = 0$$

or, equivalently

$$J_{h_a}(q) = A^T(q).$$

Thus, the first $m$ components of $h(q)$ satisfy the integrability assumption. $\square$

In Appendix B, we report an alternative proof of the sufficient part of Theorem 2, which uses algebraic arguments instead of power invariance.

**Remark:** There is no constraint on choosing the unactuated variables, except that the corresponding Jacobian is nonsingular. Indeed, the factorization given in (13) holds independently of $\theta_u$. In other words, Theorem 2 does not rely on a specific choice of $\theta_u$, which could be used to further simplify the structure of the equations of motion.
<figure><img src="image_7.png" /><figcaption>Fig. 5. Continuum soft robot discretized into two bodies. Under the PCC assumption, the configuration of each body is described by its curvature $\kappa_i$, bending angle $\phi_i$, and elongation $\delta L_i$ for $i = 1, 2$. Three tendons run from the base to the tip and spatial motion is obtained by applying a suitable cable tension $\mathbf{u} \in \mathbb{R}^3$.</figcaption></figure>

Note that the previous results can also be derived in a Hamiltonian formulation by considering the type 2 generating function [33] $G_2(\mathbf{q}, \tilde{\mathbf{p}}) = \mathbf{h}^T(\mathbf{q})\tilde{\mathbf{p}}$, where $\tilde{\mathbf{p}}$ is the momentum in the actuation coordinates.

We illustrate the application of Theorem 2 on a soft robotic arm.

**Example 5 (Continuum soft robot):** Consider a continuum soft robot discretized into two bodies, modeled under the piecewise constant curvature (PCC) hypothesis. Then, each body has three DOF, corresponding to its curvature $\kappa_i$, bending direction $\phi_i$, and elongation $\delta L_i$, $i = 1, 2$, so that

$$ \mathbf{q} = (\kappa_1 \quad \phi_1 \quad \delta L_1 \quad \kappa_2 \quad \phi_2 \quad \delta L_2)^T. $$

Three tendons that run from the base to the tip actuate the robot. Each actuator is located at a distance $d \in \mathbb{R}^+$ from the backbone and rotated from the previous by $120^\circ$, as illustrated in Fig. 5. By exploiting the principle of virtual works, one can show that $\tau_q = A(\mathbf{q})\mathbf{u}$, where

$$ A(\mathbf{q}) = \begin{pmatrix} -dc_2 & d\left(\frac{1}{2}c_2 - \frac{\sqrt{3}}{2}s_2\right) & d\left(\frac{1}{2}c_2 + \frac{\sqrt{3}}{2}s_2\right) \\ dq_1s_2 & -dq_1\left(\frac{1}{2}s_2 + \frac{\sqrt{3}}{2}c_2\right) & -dq_1\left(\frac{1}{2}s_2 - \frac{\sqrt{3}}{2}c_2\right) \\ 1 & 1 & 1 \\ -dc_5 & d\left(\frac{1}{2}c_5 - \frac{\sqrt{3}}{2}s_5\right) & d\left(\frac{1}{2}c_5 + \frac{\sqrt{3}}{2}s_5\right) \\ dq_4s_5 & -dq_4\left(\frac{1}{2}s_5 + \frac{\sqrt{3}}{2}c_5\right) & -dq_4\left(\frac{1}{2}s_5 - \frac{\sqrt{3}}{2}c_5\right) \\ 1 & 1 & 1 \end{pmatrix} \quad (15) $$

and $\mathbf{u} \in \mathbb{R}^3$ collects the cables tension. It is easy to verify that $\text{rank}(A(\mathbf{q})) = 3$ except when $q_1 = q_4 = 0$ and $q_2 = q_5 + k\pi$

with $k \in \mathbb{Z}$. In fact, control authority is lost when the arm is in the stretched configuration. However, this is an artifact due to the choice of the bending direction as a generalized coordinate [34]. Thus, the system is underactuated with $m = 3$ and $n = 6$.

It can be shown that (derivations are omitted for the sake of space) $\dot{\mathbf{y}} = A^T(\mathbf{q})\dot{\mathbf{q}} \in \mathbb{R}^3$ can be integrated as

$$ \mathbf{y} = \begin{pmatrix} q_3 + q_6 + d(q_1c_2 + q_4c_5) \\ q_3 + q_6 - \frac{d}{2}(q_1c_2 + q_4c_5) + \frac{\sqrt{3}}{2}d(q_1s_2 + q_4s_5) \\ q_3 + q_6 - \frac{d}{2}(q_1c_2 + q_4c_5) - \frac{\sqrt{3}}{2}d(q_1s_2 + q_4s_5) \end{pmatrix}. $$

According to Theorem 2, the ID problem is solvable through a coordinate change having the form

$$ \theta = \begin{pmatrix} \theta_a \\ \theta_u \end{pmatrix} = \begin{pmatrix} \mathbf{y} \\ \mathbf{\theta}_u \end{pmatrix} \quad (16) $$

where $\theta_u \in \mathbb{R}^3$ is any complement to $\mathbf{y}$. For example, a possible choice is

$$ \theta_u = \begin{pmatrix} q_3 + dq_1c_2 \\ q_3 - \frac{d}{2}q_1c_2 + \frac{\sqrt{3}}{2}dq_1s_2 \\ q_3 - \frac{d}{2}q_1c_2 - \frac{\sqrt{3}}{2}dq_1s_2 \end{pmatrix}. $$

Remarkably, the components of $\theta_a$ correspond to the change in tendon's length inside the arm. Similarly, $\theta_u$ collects the tendons elongation in the first body only.

The reader can verify that $J_h^{-T}$ takes the expression given by (17) shown at the bottom of the next page, which yields

$$ \tau_{\theta}(\theta, \mathbf{u}) = J_h^{-T}(\mathbf{q}) A(\mathbf{q}) \mathbf{u} = \begin{pmatrix} I_3 \\ O_{3 \times 3} \end{pmatrix} \mathbf{u} $$

△

## VI. INTEGRABILITY OF THREADLIKE ACTUATORS

This section shows that the conclusions drawn in Examples 4 and 5 hold for any mechanical system driven by threadlike actuators. This type of actuation is growing in popularity because it allows creating lightweight structures with high power density, and precise and distributed actuation [35]. We first prove the existence of the actuation coordinates for chains of rigid bodies. We then extend such results to continuum bodies described by reduced-order models.

Consider a mechanical system of rigid bodies with n-DOF actuated through m inelastic tendons. As described in [32, Ch. 6], one can always define m extension functions $g_i(\mathbf{q})$; $i \in \{1, \dots, m\}$, which measure the tendons displacement as a function of $\mathbf{q}$. The application of the principle of virtual works yields

$$ \tau_q = J_g^T(\mathbf{q})\mathbf{u} = A(\mathbf{q})\mathbf{u} $$

where $\mathbf{u} \in \mathbb{R}^m$ collects the tendons tension. It immediately follows that the passive output $\dot{\mathbf{y}}$ is integrable as $\mathbf{y} = g(\mathbf{q})$.

This result extends to mechanical systems with continuum bodies modeled under the geometric variable strain (GVS) technique, see [36] and [37] for a detailed presentation of all the quantities defined in the following. We denote the strain as $\xi \in \mathbb{R}^6$, where $X \in [0, L]$ is the curvilinear abscissa with $L$
<figure><img src="image_7.png" /><figcaption>Fig. 5. Continuum soft robot discretized into two bodies. Under the PCC assumption, the configuration of each body is described by its curvature $\kappa_i$, bending angle $\phi_i$, and elongation $\delta L_i$ for $i = 1, 2$. Three tendons run from the base to the tip and spatial motion is obtained by applying a suitable cable tension $\mathbf{u} \in \mathbb{R}^3$.</figcaption></figure>

Note that the previous results can also be derived in a Hamiltonian formulation by considering the type 2 generating function [33] $G_2(\mathbf{q}, \tilde{\mathbf{p}}) = \mathbf{h}^T(\mathbf{q})\tilde{\mathbf{p}}$, where $\tilde{\mathbf{p}}$ is the momentum in the actuation coordinates.

We illustrate the application of Theorem 2 on a soft robotic arm.

**Example 5 (Continuum soft robot):** Consider a continuum soft robot discretized into two bodies, modeled under the piecewise constant curvature (PCC) hypothesis. Then, each body has three DOF, corresponding to its curvature $\kappa_i$, bending direction $\phi_i$, and elongation $\delta L_i$, $i = 1, 2$, so that

$$ \mathbf{q} = (\kappa_1 \quad \phi_1 \quad \delta L_1 \quad \kappa_2 \quad \phi_2 \quad \delta L_2)^T. $$

Three tendons that run from the base to the tip actuate the robot. Each actuator is located at a distance $d \in \mathbb{R}^+$ from the backbone and rotated from the previous by $120^\circ$, as illustrated in Fig. 5. By exploiting the principle of virtual works, one can show that $\tau_q = A(\mathbf{q})\mathbf{u}$, where

$$ A(\mathbf{q}) = \begin{pmatrix} -dc_2 & d\left(\frac{1}{2}c_2 - \frac{\sqrt{3}}{2}s_2\right) & d\left(\frac{1}{2}c_2 + \frac{\sqrt{3}}{2}s_2\right) \\ dq_1s_2 & -dq_1\left(\frac{1}{2}s_2 + \frac{\sqrt{3}}{2}c_2\right) & -dq_1\left(\frac{1}{2}s_2 - \frac{\sqrt{3}}{2}c_2\right) \\ 1 & 1 & 1 \\ -dc_5 & d\left(\frac{1}{2}c_5 - \frac{\sqrt{3}}{2}s_5\right) & d\left(\frac{1}{2}c_5 + \frac{\sqrt{3}}{2}s_5\right) \\ dq_4s_5 & -dq_4\left(\frac{1}{2}s_5 + \frac{\sqrt{3}}{2}c_5\right) & -dq_4\left(\frac{1}{2}s_5 - \frac{\sqrt{3}}{2}c_5\right) \\ 1 & 1 & 1 \end{pmatrix} \quad (15) $$

and $\mathbf{u} \in \mathbb{R}^3$ collects the cables tension. It is easy to verify that $\text{rank}(A(\mathbf{q})) = 3$ except when $q_1 = q_4 = 0$ and $q_2 = q_5 + k\pi$

with $k \in \mathbb{Z}$. In fact, control authority is lost when the arm is in the stretched configuration. However, this is an artifact due to the choice of the bending direction as a generalized coordinate [34]. Thus, the system is underactuated with $m = 3$ and $n = 6$.

It can be shown that (derivations are omitted for the sake of space) $\dot{\mathbf{y}} = A^T(\mathbf{q})\dot{\mathbf{q}} \in \mathbb{R}^3$ can be integrated as

$$ \mathbf{y} = \begin{pmatrix} q_3 + q_6 + d(q_1c_2 + q_4c_5) \\ q_3 + q_6 - \frac{d}{2}(q_1c_2 + q_4c_5) + \frac{\sqrt{3}}{2}d(q_1s_2 + q_4s_5) \\ q_3 + q_6 - \frac{d}{2}(q_1c_2 + q_4c_5) - \frac{\sqrt{3}}{2}d(q_1s_2 + q_4s_5) \end{pmatrix}. $$

According to Theorem 2, the ID problem is solvable through a coordinate change having the form

$$ \theta = \begin{pmatrix} \theta_a \\ \theta_u \end{pmatrix} = \begin{pmatrix} \mathbf{y} \\ \mathbf{\theta}_u \end{pmatrix} \quad (16) $$

where $\theta_u \in \mathbb{R}^3$ is any complement to $\mathbf{y}$. For example, a possible choice is

$$ \theta_u = \begin{pmatrix} q_3 + dq_1c_2 \\ q_3 - \frac{d}{2}q_1c_2 + \frac{\sqrt{3}}{2}dq_1s_2 \\ q_3 - \frac{d}{2}q_1c_2 - \frac{\sqrt{3}}{2}dq_1s_2 \end{pmatrix}. $$

Remarkably, the components of $\theta_a$ correspond to the change in tendon's length inside the arm. Similarly, $\theta_u$ collects the tendons elongation in the first body only.

The reader can verify that $J_h^{-T}$ takes the expression given by (17) shown at the bottom of the next page, which yields

$$ \tau_{\theta}(\theta, \mathbf{u}) = J_h^{-T}(\mathbf{q})A(\mathbf{q})\mathbf{u} = \begin{pmatrix} I_3 \\ O_{3 \times 3} \end{pmatrix}\mathbf{u} $$

△

## VI. INTEGRABILITY OF THREADLIKE ACTUATORS

This section shows that the conclusions drawn in Examples 4 and 5 hold for any mechanical system driven by threadlike actuators. This type of actuation is growing in popularity because it allows creating lightweight structures with high power density, and precise and distributed actuation [35]. We first prove the existence of the actuation coordinates for chains of rigid bodies. We then extend such results to continuum bodies described by reduced-order models.

Consider a mechanical system of rigid bodies with n-DOF actuated through m inelastic tendons. As described in [32, Ch. 6], one can always define m extension functions $g_i(\mathbf{q})$; $i \in \{1, \dots, m\}$, which measure the tendons displacement as a function of $\mathbf{q}$. The application of the principle of virtual works yields

$$ \tau_q = J_g^T(\mathbf{q})\mathbf{u} = A(\mathbf{q})\mathbf{u} $$

where $\mathbf{u} \in \mathbb{R}^m$ collects the tendons tension. It immediately follows that the passive output $\dot{\mathbf{y}}$ is integrable as $\mathbf{y} = g(\mathbf{q})$.

This result extends to mechanical systems with continuum bodies modeled under the geometric variable strain (GVS) technique, see [36] and [37] for a detailed presentation of all the quantities defined in the following. We denote the strain as $\xi \in \mathbb{R}^6$, where $X \in [0, L]$ is the curvilinear abscissa with $L$
the body rest length. The GVS approach reduces the infinite-dimensional state of the system by assuming that $\xi$ admits a representation of the form

$$ \xi = \phi(X, q) \quad (18) $$

where $q \in \mathbb{R}^n$ is the configuration vector, parameterizing the strain. Under (18), the dynamic model of a continuum takes the form of (1) [38]. For threadlike actuators, such as tendons and thin fluidic chambers, the generalized actuation force is $\tau_q = A(q)u$ where

$$ A(q) = \int_{0}^{L} J_{\phi}^{T} \Phi_{a}(X, q) dX $$

and $\Phi_a(X, q) \in \mathbb{R}^{6 \times m}$ is the spatial actuation matrix, whose $i$th column

$$ (\Phi_a(X, q))_i = \begin{pmatrix} \tilde{d}_i(X) t_i(X, q) \\ t_i(X, q) \end{pmatrix} \in \mathbb{R}^6 \quad (19) $$

represents the distributed force of the $i$th actuator. In the previous expression, $d_i(X) \in \mathbb{R}^3$ and $t_i(X, q) \in \mathbb{R}^3$ are the actuator distance to the body backbone and its unit tangent vector [39], respectively. The latter can be computed as

$$ t_i(X, q) = \frac{[\hat{\xi}d_i + d'_i]_3}{\|\hat{\xi}d_i + d'_i\|} $$

where $d_i$ is expressed in homogeneous coordinates and $\cdot)' := \frac{\partial(\cdot)}{\partial X}$. Given $t_i(X, q)$, it is also possible to compute the length $L_{c_i}$ of the actuator as

$$ L_{c_i}(q) = \int_{0}^{L} t_i^T(X, q) [\hat{\xi} d_i(X) + d'_i(X)]_3 dX \quad (20) $$

or, after some manipulations

$$ L_{c_i}(q) = \int_{0}^{L} (\Phi_a(X, q))_i^T \left( \xi + \begin{pmatrix} 0_3 \\ d'_i \end{pmatrix} \right) dX. $$

The time derivative of (20) is

$$ \begin{aligned} \dot{L}_{c_i}(q) &= \int_0^L t_i^T [\hat{\xi} d_i]_3 dX = \int_0^L (\Phi_a(X, q))_i^T J_\phi dX \dot{q} \\ &= (A(q))_i^T \dot{q} \end{aligned} $$

which implies that the passive output $\dot{y} = A^T(q)\dot{q}$ is integrable as

$$ y = g(q) = \begin{pmatrix} L_{c_1}(q) \\ \vdots \\ L_{c_m}(q) \end{pmatrix} \in \mathbb{R}^m. $$

Thus, the actuation coordinates correspond to the length of the actuators, as for rigid systems. The earlier results are independent of the number of DOF and actuators. In other words, finite-dimensional models of mechanical systems actuated via tendons always admit a collocated form, independently of being fully-over-, or underactuated. Recalling that the actuator coordinates are defined up to a constant, it is also possible to consider the actuator elongation $\delta L_{c_i} := L_{c_i} - L_{c_i}^*$; $i \in \{1, \dots, m\}$, with respect to a reference length $L_{c_i}^* \in \mathbb{R}$, such as that in the stress-free configuration. This way, proprioceptive sensors such as encoders can easily measure the actuation coordinates. Consequently, the proposed change of coordinates is also helpful for control synthesis.

In Appendix C, we show that similar arguments apply to soft robots with volumetric actuators.

## VII. CONTROL OF COLLOCATED UNDERACTUATED MECHANICAL SYSTEMS

The aforementioned results prove a fact empirically observed in soft and continuum robot control. In particular, several works [40], [41], [42], [43] have shown that it is possible to obtain excellent closed-loop performance in shape and position tasks by controlling the actuator length. This is the case for both model-based and model-free approaches. However, to the best of the authors' knowledge, it has never been clarified why these coordinates represent a better choice than others, such as the curvature and bending direction. In the actuation coordinates, the dynamics is collocated, which is expected to simplify and robustify the closed loop, especially when the control law does not require significant system knowledge. When a controller is implemented in the actuation coordinates, explicit inversion of the actuation matrix is unnecessary because these coordinates inherently incorporate the inversion. It is also worth noting that any control problem formulated in the initial configuration space can be reformulated in the actuation coordinates. Remarkably, the previous considerations remain true also when the dynamics is underactuated. Furthermore, note that direct inversion of $A(q)$ is not possible in this case being $A(q)$ a tall matrix. These results

$$ J_h^{-T} = \begin{pmatrix} 0 & 0 & 0 & -\frac{2}{3d}c_5 & \frac{2}{3dq_4}s_5 & \frac{1}{3} \\ 0 & 0 & 0 & \frac{2}{3d}\left(\frac{1}{2}c_5 - \frac{\sqrt{3}}{2}s_5\right) & -\frac{2}{3dq_4}\left(\frac{1}{2}s_5 + \frac{\sqrt{3}}{2}c_5\right) & \frac{1}{3} \\ 0 & 0 & 0 & \frac{2}{3d}\left(\frac{1}{2}c_5 + \frac{\sqrt{3}}{2}s_5\right) & -\frac{2}{3dq_4}\left(\frac{1}{2}s_5 - \frac{\sqrt{3}}{2}c_5\right) & \frac{1}{3} \\ -\frac{2}{3d}c_2 & \frac{2}{3dq_1}s_2 & \frac{1}{3} & \frac{2}{3d}c_5 & -\frac{2}{3dq_4}s_5 & -\frac{1}{3} \\ \frac{2}{3d}\left(\frac{1}{2}c_2 - \frac{\sqrt{3}}{2}s_2\right) & -\frac{2}{3d}\left(\frac{1}{2}c_2 + \frac{\sqrt{3}}{2}s_2\right) & \frac{1}{3} & -\frac{2}{3d}\left(\frac{1}{2}c_5 - \frac{\sqrt{3}}{2}s_5\right) & \frac{2}{3dq_4}\left(\frac{1}{2}s_5 + \frac{\sqrt{3}}{2}c_5\right) & -\frac{1}{3} \\ \frac{2}{3d}\left(\frac{1}{2}c_2 + \frac{\sqrt{3}}{2}s_2\right) & -\frac{2}{3d}\left(\frac{1}{2}c_2 - \frac{\sqrt{3}}{2}s_2\right) & \frac{1}{3} & -\frac{2}{3d}\left(\frac{1}{2}c_5 + \frac{\sqrt{3}}{2}s_5\right) & \frac{2}{3dq_4}\left(\frac{1}{2}s_5 - \frac{\sqrt{3}}{2}c_5\right) & -\frac{1}{3} \end{pmatrix}. \quad (17) $$
the body rest length. The GVS approach reduces the infinite-dimensional state of the system by assuming that $\xi$ admits a representation of the form

$$ \xi = \phi(X, q) \qquad (18) $$

where $q \in \mathbb{R}^n$ is the configuration vector, parameterizing the strain. Under (18), the dynamic model of a continuum takes the form of (1) [38]. For threadlike actuators, such as tendons and thin fluidic chambers, the generalized actuation force is $\tau_q = A(q)u$ where

$$ A(q) = \int_{0}^{L} J_{\phi}^{T} \Phi_{a}(X, q) dX $$

and $\Phi_a(X, q) \in \mathbb{R}^{6 \times m}$ is the spatial actuation matrix, whose $i$th column

$$ (\Phi_a(X, q))_i = \begin{pmatrix} \tilde{d}_i(X) t_i(X, q) \\ t_i(X, q) \end{pmatrix} \in \mathbb{R}^6 \quad (19) $$

represents the distributed force of the $i$th actuator. In the previous expression, $d_i(X) \in \mathbb{R}^3$ and $t_i(X, q) \in \mathbb{R}^3$ are the actuator distance to the body backbone and its unit tangent vector [39], respectively. The latter can be computed as

$$ t_i(X, q) = \frac{[\hat{\xi}d_i + d'_i]_3}{\|\hat{\xi}d_i + d'_i\|} $$

where $d_i$ is expressed in homogeneous coordinates and $\cdot)' := \frac{\partial(\cdot)}{\partial X}$. Given $t_i(X, q)$, it is also possible to compute the length $L_{c_i}$ of the actuator as

$$ L_{c_i}(q) = \int_{0}^{L} t_i^T(X, q)[\hat{\xi}d_i(X) + d'_i(X)]_3 dX \quad (20) $$

or, after some manipulations

$$ L_{c_i}(q) = \int_{0}^{L} (\Phi_a(X, q))_i^T \left( \boldsymbol{\xi} + \begin{pmatrix} \mathbf{0}_3 \\ d'_i \end{pmatrix} \right) dX. $$

The time derivative of (20) is

$$ \begin{aligned} \dot{L}_{c_i}(q) &= \int_0^L t_i^T [\hat{\xi} d_i]_3 dX = \int_0^L (\Phi_a(X, q))_i^T J_\phi dX \dot{q} \\ &= (A(q))_i^T \dot{q} \end{aligned} $$

which implies that the passive output $\dot{y} = A^T(q)\dot{q}$ is integrable as

$$ y = g(q) = \begin{pmatrix} L_{c_1}(q) \\ \vdots \\ L_{c_m}(q) \end{pmatrix} \in \mathbb{R}^m. $$

Thus, the actuation coordinates correspond to the length of the actuators, as for rigid systems. The earlier results are independent of the number of DOF and actuators. In other words, finite-dimensional models of mechanical systems actuated via tendons always admit a collocated form, independently of being fully-over-, or underactuated. Recalling that the actuator coordinates are defined up to a constant, it is also possible to consider the actuator elongation $\delta L_{c_i} := L_{c_i} - L_{c_i}^*$; $i \in \{1, \dots, m\}$, with respect to a reference length $L_{c_i}^* \in \mathbb{R}$, such as that in the stress-free configuration. This way, proprioceptive sensors such as encoders can easily measure the actuation coordinates. Consequently, the proposed change of coordinates is also helpful for control synthesis.

In Appendix C, we show that similar arguments apply to soft robots with volumetric actuators.

## VII. CONTROL OF COLLOCATED UNDERACTUATED MECHANICAL SYSTEMS

The aforementioned results prove a fact empirically observed in soft and continuum robot control. In particular, several works [40], [41], [42], [43] have shown that it is possible to obtain excellent closed-loop performance in shape and position tasks by controlling the actuator length. This is the case for both model-based and model-free approaches. However, to the best of the authors' knowledge, it has never been clarified why these coordinates represent a better choice than others, such as the curvature and bending direction. In the actuation coordinates, the dynamics is collocated, which is expected to simplify and robustify the closed loop, especially when the control law does not require significant system knowledge. When a controller is implemented in the actuation coordinates, explicit inversion of the actuation matrix is unnecessary because these coordinates inherently incorporate the inversion. It is also worth noting that any control problem formulated in the initial configuration space can be reformulated in the actuation coordinates. Remarkably, the previous considerations remain true also when the dynamics is underactuated. Furthermore, note that direct inversion of $A(q)$ is not possible in this case being $A(q)$ a tall matrix. These results

$$ J_h^{-T} = \begin{pmatrix} 0 & 0 & 0 & -\frac{2}{3d}c_5 & \frac{2}{3dq_4}s_5 & \frac{1}{3} \\ 0 & 0 & 0 & \frac{2}{3d}\left(\frac{1}{2}c_5 - \frac{\sqrt{3}}{2}s_5\right) & -\frac{2}{3dq_4}\left(\frac{1}{2}s_5 + \frac{\sqrt{3}}{2}c_5\right) & \frac{1}{3} \\ 0 & 0 & 0 & \frac{2}{3d}\left(\frac{1}{2}c_5 + \frac{\sqrt{3}}{2}s_5\right) & -\frac{2}{3dq_4}\left(\frac{1}{2}s_5 - \frac{\sqrt{3}}{2}c_5\right) & \frac{1}{3} \\ -\frac{2}{3d}c_2 & \frac{2}{3dq_1}s_2 & \frac{1}{3} & \frac{2}{3d}c_5 & -\frac{2}{3dq_4}s_5 & -\frac{1}{3} \\ \frac{2}{3d}\left(\frac{1}{2}c_2 - \frac{\sqrt{3}}{2}s_2\right) & -\frac{2}{3d}\left(\frac{1}{2}c_2 + \frac{\sqrt{3}}{2}s_2\right) & \frac{1}{3} & -\frac{2}{3d}\left(\frac{1}{2}c_5 - \frac{\sqrt{3}}{2}s_5\right) & \frac{2}{3dq_4}\left(\frac{1}{2}s_5 + \frac{\sqrt{3}}{2}c_5\right) & -\frac{1}{3} \\ \frac{2}{3d}\left(\frac{1}{2}c_2 + \frac{\sqrt{3}}{2}s_2\right) & -\frac{2}{3d}\left(\frac{1}{2}c_2 - \frac{\sqrt{3}}{2}s_2\right) & \frac{1}{3} & -\frac{2}{3d}\left(\frac{1}{2}c_5 + \frac{\sqrt{3}}{2}s_5\right) & \frac{2}{3dq_4}\left(\frac{1}{2}s_5 - \frac{\sqrt{3}}{2}c_5\right) & -\frac{1}{3} \end{pmatrix}. \quad (17) $$
<figure><img src="image_8.png" /><figcaption>Fig. 6. Time evolutions of (a)–(c) actuator elongations and (d)–(f) configuration variables under the P-satI-D (24) in $\mathbf{y}$ space and the PD+ (25) in $\mathbf{q}$ space, respectively. The P-satI-D regulates the actuation coordinates to the desired target. On the other hand, the closed-loop system under the PD+ regulator is stable but has a steady-state error. (a) P-satI-D in $\mathbf{y}$ space; $t \in [0; 1]$ s. (b) P-satI-D in $\mathbf{y}$ space; $t \in [2; 3]$ s. (c) P-satI-D in $\mathbf{y}$ space; $t \in [4; 5]$ s. (d) PD+ in $\mathbf{q}$ space; $t \in [0; 1]$ s. (e) PD+ in $\mathbf{q}$ space; $t \in [2; 3]$ s. (f) PD+ in $\mathbf{q}$ space; $t \in [4; 5]$ s.</figcaption></figure>

allow extending the controllers in [10], [23], [24], [25], [26], [37] for planar underactuated mechanical systems with damping to those moving in 3-D. The following corollary formalizes such a statement for the regulators in [10] and [26].

**Corollary 2:** Consider an underactuated mechanical system satisfying the same hypotheses of Theorem 2. Suppose that there exists a dissipation function $\mathcal{F}(\mathbf{q}, \dot{\mathbf{q}})$ such that, for all $\dot{\mathbf{q}} \in T_q M$

$$ \frac{\partial \mathcal{F}(\mathbf{q}, \dot{\mathbf{q}})}{\partial \dot{\mathbf{q}}} \dot{\mathbf{q}} > 0 \qquad (21) $$

and, in the actuation coordinates

$$ \frac{\partial^2 \mathcal{L}_{\theta}(\theta, 0)}{\partial \theta_u^2} > 0. \qquad (22) $$

Let $K_P, K_D, K_I \in \mathbb{R}^{m \times m} > 0$, and $\gamma > 0$. There exist constants $\alpha > 0$ and $\bar{\gamma} > 0$ such that, if $K_P > \alpha I_m$ and $\gamma > \bar{\gamma}$, both the following regulators.

1) PD+ (with feedforward) [10]

$$ \mathbf{u} = - \left( \frac{\partial \mathcal{L}_{\theta}(\theta_{ad}, \theta_{ud}, 0, 0)}{\partial \theta_a} \right)^T + K_P(\theta_{ad} - \theta_a) - K_D \dot{\theta}_a. \qquad (23) $$

2) P-satI-D [26]

$$ \mathbf{u} = K_P (\theta_{ad} - \theta_a) - K_D \dot{\theta}_a + \frac{K_I}{\gamma} \int_0^t \tanh (\theta_{ad} - \theta_a) (z) dz \qquad (24) $$

will globally asymptotically stabilize the closed-loop system at $(\theta_a, \theta_u, \dot{\theta}_a, \dot{\theta}_u) = (\theta_{ad}, \theta_{ud}, 0, 0)$, where $\theta_{ad} \in \mathbb{R}^m$

<figure><img src="image_9.png" /><figcaption>Fig. 7. Time evolution of the control inputs under the P-satI-D (24) in $\mathbf{y}$ space and the PD+ (25) in $\mathbf{q}$ space, respectively. (a) P-satI-D in $\mathbf{y}$ space. (b) PD+ in $\mathbf{q}$ space.</figcaption></figure>

and $\theta_{ud} \in \mathbb{R}^{n-m}$ is the unique solution to

$$ \left( \frac{\partial \mathcal{L}_{\theta}(\theta_{ad}, \theta_u, 0, 0)}{\partial \theta_u} \right)^T = 0. $$

Under (22), for any value of the actuated coordinates $\theta_a$, there is a unique equilibrium of those unactuated, i.e., the system
<figure><img src="image_8.png" /><figcaption>Fig. 6. Time evolutions of (a)–(c) actuator elongations and (d)–(f) configuration variables under the P-satI-D (24) in $\mathbf{y}$ space and the PD+ (25) in $\mathbf{q}$ space, respectively. The P-satI-D regulates the actuation coordinates to the desired target. On the other hand, the closed-loop system under the PD+ regulator is stable but has a steady-state error. (a) P-satI-D in $\mathbf{y}$ space; $t \in [0; 1]$ s. (b) P-satI-D in $\mathbf{y}$ space; $t \in [2; 3]$ s. (c) P-satI-D in $\mathbf{y}$ space; $t \in [4; 5]$ s. (d) PD+ in $\mathbf{q}$ space; $t \in [0; 1]$ s. (e) PD+ in $\mathbf{q}$ space; $t \in [2; 3]$ s. (f) PD+ in $\mathbf{q}$ space; $t \in [4; 5]$ s.</figcaption></figure>

allow extending the controllers in [10], [23], [24], [25], [26], [37] for planar underactuated mechanical systems with damping to those moving in 3-D. The following corollary formalizes such a statement for the regulators in [10] and [26].

**Corollary 2:** Consider an underactuated mechanical system satisfying the same hypotheses of Theorem 2. Suppose that there exists a dissipation function $\mathcal{F}(\mathbf{q}, \dot{\mathbf{q}})$ such that, for all $\dot{\mathbf{q}} \in T_q M$

$$ \frac{\partial \mathcal{F}(\mathbf{q}, \dot{\mathbf{q}})}{\partial \dot{\mathbf{q}}} \dot{\mathbf{q}} > 0 \qquad (21) $$

and, in the actuation coordinates

$$ \frac{\partial^2 \mathcal{L}_{\theta}(\theta, 0)}{\partial \theta_u^2} > 0. \qquad (22) $$

Let $K_P, K_D, K_I \in \mathbb{R}^{m \times m} > 0$, and $\gamma > 0$. There exist constants $\alpha > 0$ and $\bar{\gamma} > 0$ such that, if $K_P > \alpha I_m$ and $\gamma > \bar{\gamma}$, both the following regulators.

1) PD+ (with feedforward) [10]

$$ \mathbf{u} = - \left( \frac{\partial \mathcal{L}_{\theta}(\theta_{ad}, \theta_{ud}, 0, 0)}{\partial \theta_a} \right)^T + \mathbf{K}_P(\theta_{ad} - \theta_a) - \mathbf{K}_D \dot{\theta}_a. \qquad (23) $$

2) P-satI-D [26]

$$ \mathbf{u} = \mathbf{K}_P (\theta_{ad} - \theta_a) - \mathbf{K}_D \dot{\theta}_a + \frac{\mathbf{K}_I}{\gamma} \int_0^t \tanh (\theta_{ad} - \theta_a) (z) dz \qquad (24) $$

will globally asymptotically stabilize the closed-loop system at $(\theta_a, \theta_u, \dot{\theta}_a, \dot{\theta}_u) = (\theta_{ad}, \theta_{ud}, 0, 0)$, where $\theta_{ad} \in \mathbb{R}^m$

<figure><img src="image_9.png" /><figcaption>Fig. 7. Time evolution of the control inputs under the P-satI-D (24) in $\mathbf{y}$ space and the PD+ (25) in $\mathbf{q}$ space, respectively. (a) P-satI-D in $\mathbf{y}$ space. (b) PD+ in $\mathbf{q}$ space.</figcaption></figure>

and $\theta_{ud} \in \mathbb{R}^{n-m}$ is the unique solution to

$$ \left( \frac{\partial \mathcal{L}_{\theta}(\theta_{ad}, \theta_u, 0, 0)}{\partial \theta_u} \right)^T = 0. $$

Under (22), for any value of the actuated coordinates $\theta_a$, there is a unique equilibrium of those unactuated, i.e., the system
<figure><img src="image_10.png" /><figcaption>Fig. 8. Frame sequences of robot motion, divided in three time windows. Panels (a)–(c) and (d)–(f) show the robot when controlled in the *y* or *q* coordinates using the laws (24) and (25), respectively. The initial and final configurations of each interval are shown in blue while light gray shapes represent intermediate configurations. A red star indicates the constant target position of the end-effector. (a) P-satI-D in *y* space; $t \in [0; 2]$ s. (b) P-satI-D in *y* space; $t \in [2; 4]$ s. (c) P-satI-D in *y* space; $t \in [4; 6]$ s. (d) PD+ in *q* space; $t \in [0; 2]$ s. (e) PD+ in *q* space; $t \in [2; 4]$ s. (f) PD+ in *q* space; $t \in [4; 6]$ s.</figcaption></figure>

equilibria are uniquely determined by $\theta_a$. Instead, (21) guarantees the internal stability of the closed-loop system. Note that the aforementioned controllers admit more general structures, see [10] and [26]. In addition, despite being developed for continuum soft robots, these apply to any underactuated mechanical system with damping on the unactuated variables.

We exploit the previous corollary to perform a shape regulation task for a continuum soft robot moving in 3-D. The robot has a rest length $L = 0.4$ m and cross-section radius $R \in [0.02; 0.008]$ m, which varies linearly from the base to the tip. The mass density is $\rho = 680$ kg/m$^3$. Furthermore, we consider a linear visco-elastic stress-strain curve with Young modulus $E = 8.88 \times 10^5$ N/m$^2$, Poisson ratio $P = 0.5$, and material damping $D = 1 \times 10^4$ N/m$^2$ s. Eight tendons actuate the robot. The first six have an oblique routing and are displaced 60° each. Their initial distance from the center line is 0.0016 m. Three of these run from the base to half of the robot while the remaining ones up to the tip. The last two tendons have a helical routing with pitch $0.4/2\pi$ m and are displaced 180°, with a distance from the backbone of 0.006 m. The strain is modeled as

$$ \xi(X, q) = \underbrace{(\Sigma^{-1}\Phi_a(X, q^*) \quad \Phi_u(X))}_{\Phi_{\xi}(X)} q + \xi^* $$

where $q \in \mathbb{R}^{15}$, $q^* = \mathbf{0}_{15}$ and $\xi^* = (0 \ 0 \ 1 \ 0 \ 0 \ 0)^T$ denote the stress-free configuration and strain, respectively, $\Sigma(X) \in \mathbb{R}^{6 \times 6}$ is the positive-definite body stiffness matrix and the columns of $\Phi_a \in \mathbb{R}^{6 \times 8}$ are defined as in (19). The strain basis

$\Sigma^{-1}\Phi_a(X, q^*)$ has proven to accurately describe the deformations due to the actuation forces [39]. Instead

$$ \Phi_u(X) = \begin{pmatrix} 1 & P_1(X) & P_2(X) & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & P_1(X) & P_2(X) & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix}_{0 \times 7} $$

with $P_1(X) := 2\frac{X}{L} - 1$ and $P_2(X) := 6(\frac{X}{L})^2 - 6\frac{X}{L} + 1$, encodes three Legendre polynomials modeling the angular deformations due to the gravitational field, not captured by $\Sigma^{-1}\Phi_a$. Since $n = 15$ and $r = m = 8$, the system is underactuated. Note that only shape regulation tasks can be achieved in general. We compare (24) with the PD+ regulator in $q$ space of [44]

$$ u = -A^{\dagger}(q_d) \left( \frac{\partial \mathcal{L}_q(q_d, 0)}{\partial q} \right)^T + A^T(q) [k_P(q_d - q) - k_D \dot{q}] \quad (25) $$

Due to $A(q)$, the aforementioned control law guarantees only local asymptotic stability [2], [44], and it requires information on the entire state of the robot to be implemented. The control gains of (24) and (25) are $K_P = k_P I_8$, $K_D = k_D I_8$, $K_I = k_I I_8$, and $\gamma = 1$, with $k_P = 2.5 \times 10^3$ N/m, $k_D = 10$ N \cdot s/m, and $k_I = 2 \times 10^3$ N/m \cdot s. Because of the underactuation, only the configurations satisfying the equilibrium equation

$$ \left( \frac{\partial \mathcal{L}_q(q_{\text{eq}}, 0)}{\partial q} \right)^T = A(q_{\text{eq}})u $$

with $u \in \mathbb{R}^8$, can be controlled. We command the three desired shapes given in (26) shown at the bottom of the next page,
<figure><img src="image_10.png" /><figcaption>Fig. 8. Frame sequences of robot motion, divided in three time windows. Panels (a)–(c) and (d)–(f) show the robot when controlled in the *y* or *q* coordinates using the laws (24) and (25), respectively. The initial and final configurations of each interval are shown in blue while light gray shapes represent intermediate configurations. A red star indicates the constant target position of the end-effector. (a) P-satI-D in *y* space; *t* ∈ [0; 2] s. (b) P-satI-D in *y* space; *t* ∈ [2; 4] s. (c) P-satI-D in *y* space; *t* ∈ [4; 6] s. (d) PD+ in *q* space; *t* ∈ [0; 2] s. (e) PD+ in *q* space; *t* ∈ [2; 4] s. (f) PD+ in *q* space; *t* ∈ [4; 6] s.</figcaption></figure>

equilibria are uniquely determined by $\theta_a$. Instead, (21) guarantees the internal stability of the closed-loop system. Note that the aforementioned controllers admit more general structures, see [10] and [26]. In addition, despite being developed for continuum soft robots, these apply to any underactuated mechanical system with damping on the unactuated variables.

We exploit the previous corollary to perform a shape regulation task for a continuum soft robot moving in 3-D. The robot has a rest length $L = 0.4$ m and cross-section radius $R \in [0.02; 0.008]$ m, which varies linearly from the base to the tip. The mass density is $\rho = 680$ kg/m$^3$. Furthermore, we consider a linear visco-elastic stress-strain curve with Young modulus $E = 8.88 \times 10^5$ N/m$^2$, Poisson ratio $P = 0.5$, and material damping $D = 1 \times 10^4$ N/m$^2$ s. Eight tendons actuate the robot. The first six have an oblique routing and are displaced 60° each. Their initial distance from the center line is 0.0016 m. Three of these run from the base to half of the robot while the remaining ones up to the tip. The last two tendons have a helical routing with pitch $0.4/2\pi$ m and are displaced 180°, with a distance from the backbone of 0.006 m. The strain is modeled as

$$ \xi(X, q) = \underbrace{(\Sigma^{-1}\Phi_a(X, q^*) \quad \Phi_u(X))}_{\Phi_{\xi}(X)} q + \xi^* $$

where $q \in \mathbb{R}^{15}$, $q^* = \mathbf{0}_{15}$ and $\xi^* = (0 \ 0 \ 1 \ 0 \ 0 \ 0)^T$ denote the stress-free configuration and strain, respectively, $\Sigma(X) \in \mathbb{R}^{6 \times 6}$ is the positive-definite body stiffness matrix and the columns of $\Phi_a \in \mathbb{R}^{6 \times 8}$ are defined as in (19). The strain basis

$\Sigma^{-1}\Phi_a(X, q^*)$ has proven to accurately describe the deformations due to the actuation forces [39]. Instead

$$ \Phi_u(X) = \begin{pmatrix} 1 & P_1(X) & P_2(X) & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & P_1(X) & P_2(X) & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix}_{\substack{0 \\ 3 \times 7}} $$

with $P_1(X) := 2\frac{X}{L} - 1$ and $P_2(X) := 6(\frac{X}{L})^2 - 6\frac{X}{L} + 1$, encodes three Legendre polynomials modeling the angular deformations due to the gravitational field, not captured by $\Sigma^{-1}\Phi_a$. Since $n = 15$ and $r = m = 8$, the system is underactuated. Note that only shape regulation tasks can be achieved in general. We compare (24) with the PD+ regulator in $q$ space of [44]

$$ u = -A^{\dagger}(q_d) \left( \frac{\partial \mathcal{L}_q(q_d, 0)}{\partial q} \right)^T + A^T(q) [k_P(q_d - q) - k_D \dot{q}] . \quad (25) $$

Due to $A(q)$, the aforementioned control law guarantees only local asymptotic stability [2], [44], and it requires information on the entire state of the robot to be implemented. The control gains of (24) and (25) are $K_P = k_P I_8$, $K_D = k_D I_8$, $K_I = k_I I_8$, and $\gamma = 1$, with $k_P = 2.5 \times 10^3$ N/m, $k_D = 10$ N \cdot s/m, and $k_I = 2 \times 10^3$ N/m \cdot s. Because of the underactuation, only the configurations satisfying the equilibrium equation

$$ \left( \frac{\partial \mathcal{L}_q(q_{\text{eq}}, 0)}{\partial q} \right)^T = A(q_{\text{eq}})u $$

with $u \in \mathbb{R}^8$, can be controlled. We command the three desired shapes given in (26) shown at the bottom of the next page,
as step references spaced in time by 2s. Furthermore, $q_{d,i}$ is converted into a desired tendon displacement $\theta_{ad,i}(t) = y_{d,i}(t)$; $i = 1, 2, 3$, for (24). The robot starts from the straight (stress-free) configuration at rest, and the simulation runs for 6s. In the following simulations, $\mathbf{y}$ has been computed through numerical integration because it was impossible to derive its closed-form expression. On the other hand, in an experimental setup equipped with motor encoders $\mathbf{y}$ could have been directly evaluated or obtained from the available measurements. Fig. 6(a)–(c) and (d)–(f) shows the evolution of the actuation coordinates and the configuration variables under (24) and (25), respectively, for three subintervals of length 1s. As expected, the P-satI-D regulates the actuation coordinates to the desired set point. On the other hand, the PD+ in $\mathbf{q}$ space fails this task, always showing a steady-state error. However, the closed-loop system remains stable. The control action for the two closed-loop systems is reported in Fig. 7. Note that the controller outputs are quite different. Finally, Fig. 8 presents a photo sequence of the two closed-loop systems. The end-effector reaches the correct position only under the P-satI-D. Indeed, the average norm of the steady-state Cartesian error is $5.8 \times 10^{-5}$ m for the P-satI-D and, respectively, $1.2 \times 10^{-2}$ m for the PD+ regulator.

VIII. CONCLUSION

This article has considered the ID problem for Lagrangian systems. We have shown that there exists a class of Lagrangian dynamics, called collocated, for which a coordinate transformation decouples actuator inputs entering the equations of motion through a configuration-dependent actuation matrix. These coordinates have a physical interpretation and can be easily computed. Under mild conditions on the differentiability of the actuation matrix, a simple test allows verifying if the dynamics is collocated or not. As a consequence of power invariance, the results equally apply to fully actuated, overactuated, and underactuated systems. In the case of underactuated dynamics, inputs are collocated with the actuation coordinates while some freedom in left in the definition of the unactuated coordinates. Since we consider only coordinate transformations, the proposed method differs from the standard differential geometric approach used for controlling input-affine nonlinear systems, which typically makes use of complete feedback transformations. As a byproduct of our approach, we have shown that all mechanical systems driven by threadlike actuators are collocated. Moreover, we were able to extend control laws recently developed for underactuated systems with constant actuation matrix to collocated mechanical systems with damping.

Future work will be devoted to the experimental validation of the proposed method and to special choices of the unactuated variables that further simplify the equations of motion, ease the

check of conditions for obtaining input-state or input-output exact linearization via feedback, or even reveal the existence of flat outputs for the system. Additionally, we will consider relaxing the integrability hypothesis at the cost of transforming the input.

APPENDIX

A. Properties of Lagrangian Systems

We recall two important properties of Lagrangian systems that play a key role in deriving the results in this article.

**Property 1:** Let $\mathcal{H}_q(\mathbf{q}, \dot{\mathbf{q}}) := \frac{\partial \mathcal{L}_q(\mathbf{q}, \dot{\mathbf{q}})}{\partial \dot{\mathbf{q}}} \dot{\mathbf{q}} - \mathcal{L}_q(\mathbf{q}, \dot{\mathbf{q}})$ be the system Hamiltonian. For all $\mathbf{q} \in \mathcal{M}, \dot{\mathbf{q}} \in T_q\mathcal{M}$, and $\mathbf{u} \in \mathbb{R}^m$, it holds

$$ \dot{\mathcal{H}}_q(\mathbf{q}, \dot{\mathbf{q}}) = \dot{\mathbf{q}}^T \tau_q(\mathbf{q}, \mathbf{u}). \quad (27) $$

Equation (27) states that the time rate of change of the Hamiltonian, i.e., the system total energy, equals the input power. It also follows from (1) and (27) that the dynamics is passive with respect to the pair $(\mathbf{u}, \dot{\mathbf{y}}) = (\mathbf{u}, A^T(\mathbf{q})\dot{\mathbf{q}})$ with storage function $\mathcal{H}_q(\mathbf{q}, \dot{\mathbf{q}})$.

**Property 2:** If $\theta: \mathcal{B}(\mathbf{q}) \to \mathcal{N} = h(\mathbf{q})$ is a (local) diffeomorphism, with $\mathbf{J}_h(\mathbf{q}) = \frac{\partial h}{\partial \dot{\mathbf{q}}}$, then

$$ \frac{d}{dt} \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \dot{\theta}} \right)^T - \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \theta} \right)^T = \tau_{\theta}(\theta, \mathbf{u}) \quad (28) $$

where

$$ \mathcal{L}_{\theta}(\theta, \dot{\theta}) = \mathcal{L}_{q}(q = h^{-1}(\theta), \dot{q} = \mathbf{J}_{h}^{-1}\dot{\theta}) $$

and

$$ \tau_{\theta}(\theta, \mathbf{u}) = \mathbf{J}_{h}^{-T} \tau_{q}(q = h^{-1}(\theta), \mathbf{u}). $$

This also implies that, for all $\theta \in \mathcal{N}, \dot{\theta} \in T_{\theta}\mathcal{N}, q \in \mathcal{M}, \dot{q} \in T_q\mathcal{M}$, and $\mathbf{u} \in \mathbb{R}^m$, it holds

$$ \dot{\mathcal{H}}_{\theta}(\theta, \dot{\theta}) = \dot{\theta}^{T} \tau_{\theta}(\theta, \mathbf{u}) = \dot{\mathbf{q}}^{T} \tau_{q}(\mathbf{q}, \mathbf{u}) = \dot{\mathcal{H}}_{q}(\mathbf{q}, \dot{\mathbf{q}}) \quad (29) $$

$$ \text{being } \mathcal{H}_{\theta}(\theta, \dot{\theta}) = \mathcal{H}_{q}(\mathbf{q} = \mathbf{h}^{-1}(\theta), \dot{\mathbf{q}} = \mathbf{J}_{h}^{-1}\dot{\theta}). $$

According to the previous property, the Euler–Lagrange equations and the power are invariant, i.e., they do not depend on the choice of coordinates representing the dynamics.

B. Alternative Proof of Theorem 2

We provide an alternative proof of the if part of Theorem 2. Similar considerations also hold for Theorem 1 and Corollary 1.

*Proof:* From Property 2, we have

$$ \tau_{\theta} = \mathbf{J}_{h}^{-T} A(\mathbf{q} = \mathbf{h}^{-1}(\theta)) \mathbf{u}. \quad (30) $$

$$ \begin{aligned} q_{d,1} &= (-8.85 \quad -4.70 \quad -1.39 \quad -26.30 \quad -23.41 \quad -26.33 \quad -1.74 \quad -1.79 \quad 0.19 \quad -0.08 \quad -0.09 \quad 0.70 \quad -0.68 \quad -0.05 \quad -1.55), \\ q_{d,2} &= (-13.21 \quad -11.74 \quad -9.91 \quad -15.40 \quad -14.20 \quad -14.25 \quad -0.26 \quad -0.72 \quad 0.09 \quad -0.35 \quad -0.11 \quad 0.20 \quad -2.08 \quad -0.75 \quad -1.75), \\ q_{d,3} &= (-12.99 \quad -11.52 \quad -20.49 \quad -22.94 \quad -19.77 \quad -20.40 \quad -0.86 \quad -0.87 \quad -0.55 \quad -0.50 \quad -0.43 \quad -2.01 \quad -1.33 \quad 0.10 \quad 1.51). \end{aligned} \quad (26) $$
as step references spaced in time by 2s. Furthermore, $q_{d,i}$ is converted into a desired tendon displacement $\theta_{ad,i}(t) = y_{d,i}(t)$; $i = 1, 2, 3$, for (24). The robot starts from the straight (stress-free) configuration at rest, and the simulation runs for 6s. In the following simulations, $\mathbf{y}$ has been computed through numerical integration because it was impossible to derive its closed-form expression. On the other hand, in an experimental setup equipped with motor encoders $\mathbf{y}$ could have been directly evaluated or obtained from the available measurements. Fig. 6(a)–(c) and (d)–(f) shows the evolution of the actuation coordinates and the configuration variables under (24) and (25), respectively, for three subintervals of length 1s. As expected, the P-satI-D regulates the actuation coordinates to the desired set point. On the other hand, the PD+ in $\mathbf{q}$ space fails this task, always showing a steady-state error. However, the closed-loop system remains stable. The control action for the two closed-loop systems is reported in Fig. 7. Note that the controller outputs are quite different. Finally, Fig. 8 presents a photo sequence of the two closed-loop systems. The end-effector reaches the correct position only under the P-satI-D. Indeed, the average norm of the steady-state Cartesian error is $5.8 \times 10^{-5}$ m for the P-satI-D and, respectively, $1.2 \times 10^{-2}$ m for the PD+ regulator.

VIII. CONCLUSION

This article has considered the ID problem for Lagrangian systems. We have shown that there exists a class of Lagrangian dynamics, called *collocated*, for which a coordinate transformation decouples actuator inputs entering the equations of motion through a configuration-dependent actuation matrix. These coordinates have a physical interpretation and can be easily computed. Under mild conditions on the differentiability of the actuation matrix, a simple test allows verifying if the dynamics is collocated or not. As a consequence of power invariance, the results equally apply to fully actuated, overactuated, and underactuated systems. In the case of underactuated dynamics, inputs are collocated with the actuation coordinates while some freedom in left in the definition of the unactuated coordinates. Since we consider only coordinate transformations, the proposed method differs from the standard differential geometric approach used for controlling input-affine nonlinear systems, which typically makes use of complete feedback transformations. As a byproduct of our approach, we have shown that all mechanical systems driven by threadlike actuators are collocated. Moreover, we were able to extend control laws recently developed for underactuated systems with constant actuation matrix to collocated mechanical systems with damping.

Future work will be devoted to the experimental validation of the proposed method and to special choices of the unactuated variables that further simplify the equations of motion, ease the

check of conditions for obtaining input-state or input-output exact linearization via feedback, or even reveal the existence of flat outputs for the system. Additionally, we will consider relaxing the integrability hypothesis at the cost of transforming the input.

APPENDIX

A. Properties of Lagrangian Systems

We recall two important properties of Lagrangian systems that play a key role in deriving the results in this article.

**Property 1:** Let $\mathcal{H}_q(\mathbf{q}, \dot{\mathbf{q}}) := \frac{\partial \mathcal{L}_q(\mathbf{q}, \dot{\mathbf{q}})}{\partial \dot{\mathbf{q}}} \dot{\mathbf{q}} - \mathcal{L}_q(\mathbf{q}, \dot{\mathbf{q}})$ be the system Hamiltonian. For all $\mathbf{q} \in \mathcal{M}, \dot{\mathbf{q}} \in T_q\mathcal{M}$, and $\mathbf{u} \in \mathbb{R}^m$, it holds

$$ \dot{\mathcal{H}}_q(\mathbf{q}, \dot{\mathbf{q}}) = \dot{\mathbf{q}}^T \tau_q(\mathbf{q}, \mathbf{u}). \quad (27) $$

Equation (27) states that the time rate of change of the Hamiltonian, i.e., the system total energy, equals the input power. It also follows from (1) and (27) that the dynamics is passive with respect to the pair $(\mathbf{u}, \dot{\mathbf{y}}) = (\mathbf{u}, A^T(\mathbf{q})\dot{\mathbf{q}})$ with storage function $\mathcal{H}_q(\mathbf{q}, \dot{\mathbf{q}})$.

**Property 2:** If $\theta: \mathcal{B}(\mathbf{q}) \to \mathcal{N} = h(\mathbf{q})$ is a (local) diffeomorphism, with $\mathbf{J}_h(\mathbf{q}) = \frac{\partial h}{\partial \dot{\mathbf{q}}}$, then

$$ \frac{d}{dt} \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \dot{\theta}} \right)^T - \left( \frac{\partial \mathcal{L}_{\theta}(\theta, \dot{\theta})}{\partial \theta} \right)^T = \tau_{\theta}(\theta, \mathbf{u}) \quad (28) $$

where

$$ \mathcal{L}_{\theta}(\theta, \dot{\theta}) = \mathcal{L}_{q}(q = h^{-1}(\theta), \dot{q} = \mathbf{J}_{h}^{-1}\dot{\theta}) $$

and

$$ \tau_{\theta}(\theta, \mathbf{u}) = \mathbf{J}_{h}^{-T} \tau_{q}(q = h^{-1}(\theta), \mathbf{u}). $$

This also implies that, for all $\theta \in \mathcal{N}, \dot{\theta} \in T_{\theta}\mathcal{N}, q \in \mathcal{M}, \dot{q} \in T_q\mathcal{M}$, and $\mathbf{u} \in \mathbb{R}^m$, it holds

$$ \dot{\mathcal{H}}_{\theta}(\theta, \dot{\theta}) = \dot{\theta}^{T} \tau_{\theta}(\theta, \mathbf{u}) = \dot{\mathbf{q}}^{T} \tau_{q}(\mathbf{q}, \mathbf{u}) = \dot{\mathcal{H}}_{q}(\mathbf{q}, \dot{\mathbf{q}}) \quad (29) $$

$$ \text{being } \mathcal{H}_{\theta}(\theta, \dot{\theta}) = \mathcal{H}_{q}(\mathbf{q} = \mathbf{h}^{-1}(\theta), \dot{\mathbf{q}} = \mathbf{J}_{h}^{-1}\dot{\theta}). $$

According to the previous property, the Euler–Lagrange equations and the power are invariant, i.e., they do not depend on the choice of coordinates representing the dynamics.

B. Alternative Proof of Theorem 2

We provide an alternative proof of the if part of Theorem 2. Similar considerations also hold for Theorem 1 and Corollary 1.

*Proof:* From Property 2, we have

$$ \tau_{\theta} = \mathbf{J}_{h}^{-T} A(\mathbf{q} = \mathbf{h}^{-1}(\theta)) \mathbf{u}. \quad (30) $$

$$ \begin{aligned} \mathbf{q}_{d,1} &= (-8.85 \quad -4.70 \quad -1.39 \quad -26.30 \quad -23.41 \quad -26.33 \quad -1.74 \quad -1.79 \quad 0.19 \quad -0.08 \quad -0.09 \quad 0.70 \quad -0.68 \quad -0.05 \quad -1.55), \\ \mathbf{q}_{d,2} &= (-13.21 \quad -11.74 \quad -9.91 \quad -15.40 \quad -14.20 \quad -14.25 \quad -0.26 \quad -0.72 \quad 0.09 \quad -0.35 \quad -0.11 \quad 0.20 \quad -2.08 \quad -0.75 \quad -1.75), \\ \mathbf{q}_{d,3} &= (-12.99 \quad -11.52 \quad -20.49 \quad -22.94 \quad -19.77 \quad -20.40 \quad -0.86 \quad -0.87 \quad -0.55 \quad -0.50 \quad -0.43 \quad -2.01 \quad -1.33 \quad 0.10 \quad 1.51). \end{aligned} \quad (26) $$
By exploiting the block triangular structure of $J_h(q)$, it follows

$$J_h^{-T}(q) = \begin{pmatrix} A_a^{-1}(q) & O_{m \times (n-m)} \\ -A_u A_a^{-1}(q) & I_{n-m} \end{pmatrix}$$

which yields

$$J_h^{-T} A = \begin{pmatrix} A_a^{-1}(q) & O_{m \times (n-m)} \\ -A_u A_a^{-1}(q) & I_{n-m} \end{pmatrix} \begin{pmatrix} A_a(q) \\ A_u(q) \end{pmatrix} = \begin{pmatrix} I_m \\ O_{(n-m) \times m} \end{pmatrix}.$$

<figure><img src="image_11.png" /></figure>

### C. Integrability of Volumetric Actuators

Following arguments similar to those in [45], it is possible to extend the results of Section VI to robotic systems with volumetric actuators [46].

Let $V_i^*$ be the volume of the ith actuator inside the robot when in the reference configuration. Similarly, denote with $V_i(q)$ the volume in the current deformed configuration. The work performed by the actuator on the robot is

$$W_{u_i} = (V_i(q) - V_i^*) u_i.$$

To determine the effect of $u_i$ on the generalized coordinates, it is possible to apply the principle of virtual works

$$\delta W_{(\tau_q)_i} = \delta \boldsymbol{q}^T (\tau_q)_i = \delta \boldsymbol{q}^T \left( \frac{\partial V_i(q)}{\partial \boldsymbol{q}} \right)^T u_i = \delta W_{u_i}$$

obtaining

$$A_i(q) = \left( \frac{\partial V_i(q)}{\partial q} \right)^T.$$

The previous equations imply that $\dot{\boldsymbol{y}}_i = A_i^T(q)\dot{\boldsymbol{q}}$ is integrable and the corresponding actuation coordinate can be chosen as $\delta V_i := V_i(q) - V_i^*$, which is the volume variation in the actuator chamber.

### REFERENCES

[1] G. Chen and F. L. Lewis, “Distributed adaptive tracking control for synchronization of unknown networked Lagrangian systems,” *IEEE Trans. Syst., Man, Cybern.*, vol. 41, no. 3, pp. 805–816, Jun. 2011.

[2] R. Ortega, J. A. L. Perez, P. J. Nicklasson, and H. J. Sira-Ramirez, *Passivity-Based Control of Euler-Lagrange Systems: Mechanical, Electrical and Electromechanical Applications*. Berlin, Germany: Springer, 2013.

[3] A. Loría, “Observers are unnecessary for output-feedback control of Lagrangian systems,” *IEEE Trans. Autom. Control*, vol. 61, no. 4, pp. 905–920, Apr. 2016.

[4] M.-F. Ge, Z.-W. Liu, G. Wen, X. Yu, and T. Huang, “Hierarchical controller-estimator for coordination of networked Euler-Lagrange systems,” *IEEE Trans. Cybern.*, vol. 50, no. 6, pp. 2450–2461, Jun. 2020.

[5] Y. Sun, D. Dong, H. Qin, and W. Wang, “Distributed tracking control for multiple Euler-Lagrange systems with communication delays and input saturation,” *ISA Trans.*, vol. 96, pp. 245–254, 2020.

[6] A. M. Giordano, C. Ott, and A. Albu-Schäffer, “Coordinated control of spacecraft’s attitude and end-effector for space robots,” *IEEE Robot. Autom. Lett.*, vol. 4, no. 2, pp. 2108–2115, Apr. 2019.

[7] B. Yi, R. Ortega, I. R. Manchester, and H. Sigueridjane, “Path following of a class of underactuated mechanical systems via immersion and invariance-based orbital stabilization,” *Int. J. Robust Nonlinear Control*, vol. 30, no. 18, pp. 8521–8544, 2020.

[8] R. Mengacci, F. Angelini, M. G. Catalano, G. Grioli, A. Bicchi, and M. Garabini, “On the motion/stiffness decoupling property of articulated soft robots with application to model-free torque iterative learning control,” *Int. J. Robot. Res.*, vol. 40, no. 1, pp. 348–374, 2021.

[9] M. Keppler, C. Ott, and A. Albu-Schäffer, “From underactuation to quasi-full actuation: Aiming at a unifying control framework for articulated soft robots,” *Int. J. Robust Nonlinear Control*, vol. 32, no. 9, pp. 5453–5484, 2022.

[10] P. Borja, A. Dabiri, and C. Della Santina, “Energy-based shape regulation of soft robots with unactuated dynamics dominated by elasticity,” in *Proc. IEEE 5th Int. Conf. Soft Robot.*, 2022, pp. 396–402.

[11] W. Zhang, Y. Wang, Y. Liu, and W. Zhang, “Multivariable disturbance observer-based $H_2$ analytical decoupling control design for multivariable systems,” *Int. J. Syst. Sci.*, vol. 47, no. 1, pp. 179–193, 2016.

[12] A. J. Krener, “On the equivalence of control systems and the linearization of nonlinear systems,” *SIAM J. Control*, vol. 11, no. 4, pp. 670–676, 1973.

[13] R. W. Brockett, “Feedback invariants for nonlinear systems,” *IFAC Proc.*, vol. 11, no. 1, pp. 1115–1120, 1978.

[14] A. Isidori, A. Krener, C. Gori-Giorgi, and S. Monaco, “Nonlinear decoupling via feedback: A differential geometric approach,” *IEEE Trans. Autom. Control*, vol. TAC-26, no. 2, pp. 331–345, Apr. 1981.

[15] H. Nijmeijer and A. Van der Schaft, *Nonlinear Dynamical Control Systems*. New York, NY, USA: Springer, 1990.

[16] A. Isidori, *Nonlinear Control Systems: An Introduction*. New York, NY, USA: Springer, 1992.

[17] R. Ortega, J. G. Romero, P. Borja, and A. Donaire, *PID Passivity-Based Control of Nonlinear Systems With Applications*. Hoboken, NJ, USA: Wiley, 2021.

[18] J. Wu and Y. Lu, “Adaptive backstepping sliding mode control for boost converter with constant power load,” *IEEE Access*, vol. 7, pp. 50797–50807, 2019.

[19] A. Lazrak and A. Abbou, “An improved control strategy for DFIG wind turbine to ride-through voltage dips,” in *Proc. IEEE 6th Int. Renewable Sustain. Energy Conf.*, 2018, pp. 1–6.

[20] R. N. Mishra and K. B. Mohanty, “Design and implementation of a feedback linearization controlled IM drive via simplified neuro-fuzzy approach,” *IETE J. Res.*, vol. 64, no. 2, pp. 209–230, 2018.

[21] X. Ding, J. Liang, S. Lu, F. Kong, and Y. Chen, “Robust back-stepping sliding mode control for LCL-type grid-connected inverters in weak grids,” *J. Power Electron.*, vol. 23, pp. 758–768, 2022.

[22] S. Skogestad, C. Zotică, and N. Alsop, “Transformed inputs for linearization, decoupling and feedforward control,” *J. Process Control*, vol. 122, pp. 113–133, 2023.

[23] E. Franco and A. Garriga-Casanovas, “Energy-shaping control of soft continuum manipulators with in-plane disturbances,” *Int. J. Robot. Res.*, vol. 40, no. 1, pp. 236–255, 2021.

[24] B. Caasenbrood, A. Pogromsky, and H. Nijmeijer, “Energy-based control for soft manipulators using Cosserat-beam models,” in *Proc. 18th Int. Conf. Informat. Control, Automat., Robot.*, 2021, pp. 311–319.

[25] P. Pustina, C. Della Santina, and A. De Luca, “Feedback regulation of elastically decoupled underactuated soft robots,” *IEEE Robot. Autom. Lett.*, vol. 7, no. 2, pp. 4512–4519, Apr. 2022.

[26] P. Pustina, P. Borja, C. Della Santina, and A. De Luca, “P-sat-D shape regulation of soft robots,” *IEEE Robot. Autom. Lett.*, vol. 8, no. 1, pp. 1–8, Jan. 2023.

[27] G. Soleti, J. Prechtl, P. R. Massenio, M. Baltes, and G. Rizzello, “Energy based control of a bi-stable and underactuated soft robotic system based on dielectric elastomer actuators,” in *Proc. 22th IFAC World Congr.*, 2023, pp. 7796–7801.

[28] A. De Luca, “Flexible robots,” in *Encyclopedia of Systems and Control*, 2nd ed., J. Baillieul and T. Samad, Eds. Berlin, Germany: Springer, 2021, pp. 814–822.

[29] A. Tsolakis and T. Keviczky, “Distributed IDA-PBC for a class of non-holonomic mechanical systems,” *IFAC-PapersOnLine*, vol. 54, no. 14, pp. 275–280, 2021.

[30] M. P. do Carmo, *Differential Forms and Applications*. New York, NY, USA: Springer, 1998.

[31] A. De Luca, R. Mattone, and G. Oriolo, “Steering a class of redundant mechanisms through end-effector generalized forces,” *IEEE Trans. Robot. Autom.*, vol. 14, no. 2, pp. 329–335, Apr. 1998.

[32] R. M. Murray, Z. Li, and S. S. Sastry, *A Mathematical Introduction to Robotic Manipulation*. Boca Raton, FL, USA: CRC, 1994.

[33] H. Goldstein, *Classical Mechanics*. Reading, MA, USA: Addison-Wesley, 1980.
By exploiting the block triangular structure of $J_h(q)$, it follows

$$J_h^{-T}(q) = \begin{pmatrix} A_a^{-1}(q) & O_{m \times (n-m)} \\ -A_u A_a^{-1}(q) & I_{n-m} \end{pmatrix}$$

which yields

$$J_h^{-T} A = \begin{pmatrix} A_a^{-1}(q) & O_{m \times (n-m)} \\ -A_u A_a^{-1}(q) & I_{n-m} \end{pmatrix} \begin{pmatrix} A_a(q) \\ A_u(q) \end{pmatrix} \\ = \begin{pmatrix} I_m \\ O_{(n-m) \times m} \end{pmatrix}.$$

<figure><img src="image_11.png" /></figure>

### C. Integrability of Volumetric Actuators

Following arguments similar to those in [45], it is possible to extend the results of Section VI to robotic systems with volumetric actuators [46].

Let $V_i^*$ be the volume of the ith actuator inside the robot when in the reference configuration. Similarly, denote with $V_i(q)$ the volume in the current deformed configuration. The work performed by the actuator on the robot is

$$W_{u_i} = (V_i(q) - V_i^*) u_i.$$

To determine the effect of $u_i$ on the generalized coordinates, it is possible to apply the principle of virtual works

$$\delta W_{(\tau_q)_i} = \delta \boldsymbol{q}^T (\tau_q)_i = \delta \boldsymbol{q}^T \left( \frac{\partial V_i(\boldsymbol{q})}{\partial \boldsymbol{q}} \right)^T u_i = \delta W_{u_i}$$

obtaining

$$A_i(\boldsymbol{q}) = \left( \frac{\partial V_i(\boldsymbol{q})}{\partial \boldsymbol{q}} \right)^T.$$

The previous equations imply that $\dot{\boldsymbol{y}}_i = A_i^T(q)\dot{\boldsymbol{q}}$ is integrable and the corresponding actuation coordinate can be chosen as $\delta V_i := V_i(q) - V_i^*$, which is the volume variation in the actuator chamber.

### REFERENCES

[1] G. Chen and F. L. Lewis, “Distributed adaptive tracking control for synchronization of unknown networked Lagrangian systems,” *IEEE Trans. Syst., Man, Cybern.*, vol. 41, no. 3, pp. 805–816, Jun. 2011.

[2] R. Ortega, J. A. L. Perez, P. J. Nicklasson, and H. J. Sira-Ramirez, *Passivity-Based Control of Euler-Lagrange Systems: Mechanical, Electrical and Electromechanical Applications*. Berlin, Germany: Springer, 2013.

[3] A. Loría, “Observers are unnecessary for output-feedback control of Lagrangian systems,” *IEEE Trans. Autom. Control*, vol. 61, no. 4, pp. 905–920, Apr. 2016.

[4] M.-F. Ge, Z.-W. Liu, G. Wen, X. Yu, and T. Huang, “Hierarchical controller-estimator for coordination of networked Euler-Lagrange systems,” *IEEE Trans. Cybern.*, vol. 50, no. 6, pp. 2450–2461, Jun. 2020.

[5] Y. Sun, D. Dong, H. Qin, and W. Wang, “Distributed tracking control for multiple Euler-Lagrange systems with communication delays and input saturation,” *ISA Trans.*, vol. 96, pp. 245–254, 2020.

[6] A. M. Giordano, C. Ott, and A. Albu-Schäffer, “Coordinated control of spacecraft’s attitude and end-effector for space robots,” *IEEE Robot. Autom. Lett.*, vol. 4, no. 2, pp. 2108–2115, Apr. 2019.

[7] B. Yi, R. Ortega, I. R. Manchester, and H. Sigueridjane, “Path following of a class of underactuated mechanical systems via immersion and invariance-based orbital stabilization,” *Int. J. Robust Nonlinear Control*, vol. 30, no. 18, pp. 8521–8544, 2020.

[8] R. Mengacci, F. Angelini, M. G. Catalano, G. Grioli, A. Bicchi, and M. Garabini, “On the motion/stiffness decoupling property of articulated soft robots with application to model-free torque iterative learning control,” *Int. J. Robot. Res.*, vol. 40, no. 1, pp. 348–374, 2021.

[9] M. Keppler, C. Ott, and A. Albu-Schäffer, “From underactuation to quasi-full actuation: Aiming at a unifying control framework for articulated soft robots,” *Int. J. Robust Nonlinear Control*, vol. 32, no. 9, pp. 5453–5484, 2022.

[10] P. Borja, A. Dabiri, and C. Della Santina, “Energy-based shape regulation of soft robots with unactuated dynamics dominated by elasticity,” in *Proc. IEEE 5th Int. Conf. Soft Robot.*, 2022, pp. 396–402.

[11] W. Zhang, Y. Wang, Y. Liu, and W. Zhang, “Multivariable disturbance observer-based $H_2$ analytical decoupling control design for multivariable systems,” *Int. J. Syst. Sci.*, vol. 47, no. 1, pp. 179–193, 2016.

[12] A. J. Krener, “On the equivalence of control systems and the linearization of nonlinear systems,” *SIAM J. Control*, vol. 11, no. 4, pp. 670–676, 1973.

[13] R. W. Brockett, “Feedback invariants for nonlinear systems,” *IFAC Proc.*, vol. 11, no. 1, pp. 1115–1120, 1978.

[14] A. Isidori, A. Krener, C. Gori-Giorgi, and S. Monaco, “Nonlinear decoupling via feedback: A differential geometric approach,” *IEEE Trans. Autom. Control*, vol. TAC-26, no. 2, pp. 331–345, Apr. 1981.

[15] H. Nijmeijer and A. Van der Schaft, *Nonlinear Dynamical Control Systems*. New York, NY, USA: Springer, 1990.

[16] A. Isidori, *Nonlinear Control Systems: An Introduction*. New York, NY, USA: Springer, 1992.

[17] R. Ortega, J. G. Romero, P. Borja, and A. Donaire, *PID Passivity-Based Control of Nonlinear Systems With Applications*. Hoboken, NJ, USA: Wiley, 2021.

[18] J. Wu and Y. Lu, “Adaptive backstepping sliding mode control for boost converter with constant power load,” *IEEE Access*, vol. 7, pp. 50797–50807, 2019.

[19] A. Lazrak and A. Abbou, “An improved control strategy for DFIG wind turbine to ride-through voltage dips,” in *Proc. IEEE 6th Int. Renewable Sustain. Energy Conf.*, 2018, pp. 1–6.

[20] R. N. Mishra and K. B. Mohanty, “Design and implementation of a feedback linearization controlled IM drive via simplified neuro-fuzzy approach,” *IETE J. Res.*, vol. 64, no. 2, pp. 209–230, 2018.

[21] X. Ding, J. Liang, S. Lu, F. Kong, and Y. Chen, “Robust back-stepping sliding mode control for LCL-type grid-connected inverters in weak grids,” *J. Power Electron.*, vol. 23, pp. 758–768, 2022.

[22] S. Skogestad, C. Zotică, and N. Alsop, “Transformed inputs for linearization, decoupling and feedforward control,” *J. Process Control*, vol. 122, pp. 113–133, 2023.

[23] E. Franco and A. Garriga-Casanovas, “Energy-shaping control of soft continuum manipulators with in-plane disturbances,” *Int. J. Robot. Res.*, vol. 40, no. 1, pp. 236–255, 2021.

[24] B. Caasenbrood, A. Pogromsky, and H. Nijmeijer, “Energy-based control for soft manipulators using Cosserat-beam models,” in *Proc. 18th Int. Conf. Informat. Control, Automat., Robot.*, 2021, pp. 311–319.

[25] P. Pustina, C. Della Santina, and A. De Luca, “Feedback regulation of elastically decoupled underactuated soft robots,” *IEEE Robot. Autom. Lett.*, vol. 7, no. 2, pp. 4512–4519, Apr. 2022.

[26] P. Pustina, P. Borja, C. Della Santina, and A. De Luca, “P-sat-D shape regulation of soft robots,” *IEEE Robot. Autom. Lett.*, vol. 8, no. 1, pp. 1–8, Jan. 2023.

[27] G. Soleti, J. Prechtl, P. R. Massenio, M. Baltes, and G. Rizzello, “Energy based control of a bi-stable and underactuated soft robotic system based on dielectric elastomer actuators,” in *Proc. 22th IFAC World Congr.*, 2023, pp. 7796–7801.

[28] A. De Luca, “Flexible robots,” in *Encyclopedia of Systems and Control*, 2nd ed., J. Baillieul and T. Samad, Eds. Berlin, Germany: Springer, 2021, pp. 814–822.

[29] A. Tsolakis and T. Keviczky, “Distributed IDA-PBC for a class of non-holonomic mechanical systems,” *IFAC-PapersOnLine*, vol. 54, no. 14, pp. 275–280, 2021.

[30] M. P. do Carmo, *Differential Forms and Applications*. New York, NY, USA: Springer, 1998.

[31] A. De Luca, R. Mattone, and G. Oriolo, “Steering a class of redundant mechanisms through end-effector generalized forces,” *IEEE Trans. Robot. Autom.*, vol. 14, no. 2, pp. 329–335, Apr. 1998.

[32] R. M. Murray, Z. Li, and S. S. Sastry, *A Mathematical Introduction to Robotic Manipulation*. Boca Raton, FL, USA: CRC, 1994.

[33] H. Goldstein, *Classical Mechanics*. Reading, MA, USA: Addison-Wesley, 1980.
[34] C. Della Santina, A. Bicchi, and D. Rus, "On an improved state parametrization for soft robots with piecewise constant curvature and its use in model based control," *IEEE Robot. Autom. Lett.*, vol. 5, no. 2, pp. 1001-1008, Apr. 2020.

[35] G. Palli, G. Borghesan, and C. Melchiorri, "Modeling, identification, and control of tendon-based actuation systems," *IEEE Trans. Robot.*, vol. 28, no. 2, pp. 277-290, Apr. 2012.

[36] F. Renda, C. Armanini, V. Lebastard, F. Candelier, and F. Boyer, "A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation," *IEEE Robot. Autom. Lett.*, vol. 5, no. 3, pp. 4006-4013, Jul. 2020.

[37] F. Boyer, V. Lebastard, F. Candelier, and F. Renda, "Dynamics of continuum and soft robots: A strain parameterization based approach," *IEEE Trans. Robot.*, vol. 37, no. 3, pp. 847-863, Jun. 2021.

[38] C. Armanini, F. Boyer, A. T. Mathew, C. Duriez, and F. Renda, "Soft robots modeling: A structured overview," *IEEE Trans. Robot.*, vol. 39, no. 3, pp. 1728-1748, Jun. 2023.

[39] F. Renda, C. Armanini, A. Mathew, and F. Boyer, "Geometrically-exact inverse kinematic control of soft manipulators with general threadlike actuators' routing," *IEEE Robot. Autom. Lett.*, vol. 7, no. 3, pp. 7311-7318, Jul. 2022.

[40] B. A. Jones and I. D. Walker, "Kinematics for multisection continuum robots," *IEEE Trans. Robot.*, vol. 22, no. 1, pp. 43-55, Feb. 2006.

[41] D. Braganza, D. M. Dawson, I. D. Walker, and N. Nath, "A neural network controller for continuum robots," *IEEE Trans. Robot.*, vol. 23, no. 6, pp. 1270-1277, Dec. 2007.

[42] V. Falkenhahn, A. Hildebrandt, R. Neumann, and O. Sawodny, "Dynamic control of the bionic handling assistant," *IEEE/ASME Trans. Mechatronics*, vol. 22, no. 1, pp. 6-17, Feb. 2017.

[43] M. Li, R. Kang, S. Geng, and E. Guglielmino, "Design and control of a tendon-driven continuum robot," *Trans. Inst. Meas. Control*, vol. 40, no. 11, pp. 3263-3272, 2018.

[44] C. Della Santina, C. Duriez, and D. Rus, "Model based control of soft robots: A survey of the state of the art and open challenges," *IEEE Control Syst. Mag.*, vol. 43, no. 3, pp. 30-65, Jun. 2023.

[45] M. Stölzle and C. Della Santina, "Piston-driven pneumatically-actuated soft robots: Modeling and backstepping control," *IEEE Control Syst. Lett.*, vol. 6, pp. 1837-1842, Dec. 2021.

[46] H. Su et al., "Pneumatic soft robots: Challenges and benefits," *Actuators*, vol. 11, no. 3, pp. 1-26, 2022.

<figure><img src="image_12.png" /></figure>

**Pietro Pustina** (Graduate Student Member, IEEE) received the M.Sc. degree in control engineering in 2021 from the Sapienza University of Rome, Rome, Italy, where he has been working toward the Ph.D. degree in automatic control since 2021.

His research interests include modeling and control of continuum soft robots.

<figure><img src="image_13.png" /></figure>

**Cosimo Della Santina** (Senior Member, IEEE) received the Ph.D. degree (*cum laude*) in robotics from the University of Pisa, Pisa, Italy, in 2019.

He is currently an Assistant Professor with TU Delft, Delft, The Netherlands, and a Research Scientist with German Aerospace Institute (DLR), Cologne, Germany. He was a visiting Ph.D. student and a Postdoc with Computer Science and Artificial Intelligence Laboratory, Massachusetts Institute of Technology, from 2017 to 2019. He was a Senior Postdoc and a Guest Lecturer with the Department

of Informatics, Technical University of Munich, in 2020 and 2021, respectively. His research interest is in providing motor intelligence to physical systems, focusing on elastic and soft robots.

Dr. Santina is a recipient of several awards, including the euRobotics Georges Giralt Ph.D. Award in 2020 and the IEEE RAS Early Academic Career Award in 2023. He is the PI for TU Delft of the European Projects Natural Intelligence and EMERGE, the co-Director of Delft AI Lab SELF, and involved in several Dutch projects.

<figure><img src="image_14.png" /></figure>

**Frédéric Boyer** (Member, IEEE) was born in France in 1967. He received the Diploma degree in mechanical engineering from the Institut Nationale Polytechnique de Grenoble, Grenoble, France, in 1991, the Master of Research degree in mechanics from the University of Grenoble, Grenoble, in 1991, and the Ph.D. degree in robotics from the University of Paris VI, Paris, France, in 1994.

He is currently a Professor with the Department of Automatic Control, Institut Mines Telecom Atlan-tique, Nantes, France, and a coleader of Bio-Robotics

Team, Laboratoire des Sciences du Numérique de Nantes. He has coordinated several research projects including one European FP7-FET project on a re-configurable eel-like robot able to navigate with electric sense. His current research interests include continuum and soft robotics, geometric mechanics, and biorobotics.

Dr. Boyer is the recipient of the Monpetit Prize from the Academy of Science of Paris in 2007 for his work in dynamics and the French "La Recherche Prize" in 2014 for his works on artificial electric sense.

<figure><img src="image_15.png" /></figure>

**Alessandro De Luca** (Life Fellow, IEEE) received the Ph.D. degree in systems engineering from the University of Rome, Rome, Italy, in 1987. Since 2000, he is currently a Professor of Robotics and Automation with the Sapienza University of Rome, Rome, Italy. His research interests include modeling, motion planning, and control of robotic systems (flexible manipulators, kinematically redundant arms, underactuated robots, and wheeled mobile robots), as well as physical human-robot interaction.

Dr. De Luca was the first Editor-in-Chief for the IEEE TRANSACTIONS ON ROBOTICS from 2004 to 2008, the RAS Vice-President for Publication Activities during 2012-2013, the General Chair of ICRA 2007, and the Program Chair of ICRA 2016. He is a recipient of three conference awards (Best Paper at ICRA 1998 and BioRob 2012, Best Application Paper at IROS 2008), the Helmholtz Humboldt Research Award in 2005, the IEEE-RAS Distinguished Service Award in 2009, and the IEEE George Saridis Leadership Award in Robotics and Automation in 2019. He was the Scientific Coordinator of the FP7 project SAPHARI—Safe and Autonomous Physical Human-Aware Robot Interaction from 2011 to 2015.

<figure><img src="image_16.png" /></figure>

**Federico Renda** (Member, IEEE) received the B.Sc. and M.Sc. degrees in biomedical engineering from the University of Pisa, Pisa, Italy, in 2007 and 2009, respectively, and the Ph.D. degree in biorobotics from Biorobotics Institute, Scuola Superiore Sant'Anna, Pisa, in 2014.

He is currently an Associate Professor with the Department of Mechanical and Nuclear Engineering, Khalifa University, Abu Dhabi, United Arab Emirates. His research interests include dynamic modeling and control of soft and continuum robots using

principles of geometric mechanics.

Open Access funding provided by Università degli Studi di Roma "La Sapienza" 2' within the CRUI CARE Agreement
[34] C. Della Santina, A. Bicchi, and D. Rus, "On an improved state parametrization for soft robots with piecewise constant curvature and its use in model based control," *IEEE Robot. Autom. Lett.*, vol. 5, no. 2, pp. 1001-1008, Apr. 2020.

[35] G. Palli, G. Borghesan, and C. Melchiorri, "Modeling, identification, and control of tendon-based actuation systems," *IEEE Trans. Robot.*, vol. 28, no. 2, pp. 277-290, Apr. 2012.

[36] F. Renda, C. Armanini, V. Lebastard, F. Candelier, and F. Boyer, "A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation," *IEEE Robot. Autom. Lett.*, vol. 5, no. 3, pp. 4006-4013, Jul. 2020.

[37] F. Boyer, V. Lebastard, F. Candelier, and F. Renda, "Dynamics of continuum and soft robots: A strain parameterization based approach," *IEEE Trans. Robot.*, vol. 37, no. 3, pp. 847-863, Jun. 2021.

[38] C. Armanini, F. Boyer, A. T. Mathew, C. Duriez, and F. Renda, "Soft robots modeling: A structured overview," *IEEE Trans. Robot.*, vol. 39, no. 3, pp. 1728-1748, Jun. 2023.

[39] F. Renda, C. Armanini, A. Mathew, and F. Boyer, "Geometrically-exact inverse kinematic control of soft manipulators with general threadlike actuators' routing," *IEEE Robot. Autom. Lett.*, vol. 7, no. 3, pp. 7311-7318, Jul. 2022.

[40] B. A. Jones and I. D. Walker, "Kinematics for multisection continuum robots," *IEEE Trans. Robot.*, vol. 22, no. 1, pp. 43-55, Feb. 2006.

[41] D. Braganza, D. M. Dawson, I. D. Walker, and N. Nath, "A neural network controller for continuum robots," *IEEE Trans. Robot.*, vol. 23, no. 6, pp. 1270-1277, Dec. 2007.

[42] V. Falkenhahn, A. Hildebrandt, R. Neumann, and O. Sawodny, "Dynamic control of the bionic handling assistant," *IEEE/ASME Trans. Mechatronics*, vol. 22, no. 1, pp. 6-17, Feb. 2017.

[43] M. Li, R. Kang, S. Geng, and E. Guglielmino, "Design and control of a tendon-driven continuum robot," *Trans. Inst. Meas. Control*, vol. 40, no. 11, pp. 3263-3272, 2018.

[44] C. Della Santina, C. Duriez, and D. Rus, "Model based control of soft robots: A survey of the state of the art and open challenges," *IEEE Control Syst. Mag.*, vol. 43, no. 3, pp. 30-65, Jun. 2023.

[45] M. Stölzle and C. Della Santina, "Piston-driven pneumatically-actuated soft robots: Modeling and backstepping control," *IEEE Control Syst. Lett.*, vol. 6, pp. 1837-1842, Dec. 2021.

[46] H. Su et al., "Pneumatic soft robots: Challenges and benefits," *Actuators*, vol. 11, no. 3, pp. 1-26, 2022.

<figure><img src="image_12.png" /></figure>

**Pietro Pustina** (Graduate Student Member, IEEE) received the M.Sc. degree in control engineering in 2021 from the Sapienza University of Rome, Rome, Italy, where he has been working toward the Ph.D. degree in automatic control since 2021.

His research interests include modeling and control of continuum soft robots.

<figure><img src="image_13.png" /></figure>

**Cosimo Della Santina** (Senior Member, IEEE) received the Ph.D. degree (*cum laude*) in robotics from the University of Pisa, Pisa, Italy, in 2019.

He is currently an Assistant Professor with TU Delft, Delft, The Netherlands, and a Research Scientist with German Aerospace Institute (DLR), Cologne, Germany. He was a visiting Ph.D. student and a Postdoc with Computer Science and Artificial Intelligence Laboratory, Massachusetts Institute of Technology, from 2017 to 2019. He was a Senior Postdoc and a Guest Lecturer with the Department

of Informatics, Technical University of Munich, in 2020 and 2021, respectively. His research interest is in providing motor intelligence to physical systems, focusing on elastic and soft robots.

Dr. Santina is a recipient of several awards, including the euRobotics Georges Giralt Ph.D. Award in 2020 and the IEEE RAS Early Academic Career Award in 2023. He is the PI for TU Delft of the European Projects Natural Intelligence and EMERGE, the co-Director of Delft AI Lab SELF, and involved in several Dutch projects.

<figure><img src="image_14.png" /></figure>

**Frédéric Boyer** (Member, IEEE) was born in France in 1967. He received the Diploma degree in mechanical engineering from the Institut Nationale Polytechnique de Grenoble, Grenoble, France, in 1991, the Master of Research degree in mechanics from the University of Grenoble, Grenoble, in 1991, and the Ph.D. degree in robotics from the University of Paris VI, Paris, France, in 1994.

He is currently a Professor with the Department of Automatic Control, Institut Mines Telecom Atlan-tique, Nantes, France, and a coleader of Bio-Robotics

Team, Laboratoire des Sciences du Numérique de Nantes. He has coordinated several research projects including one European FP7-FET project on a re-configurable eel-like robot able to navigate with electric sense. His current research interests include continuum and soft robotics, geometric mechanics, and biorobotics.

Dr. Boyer is the recipient of the Monpetit Prize from the Academy of Science of Paris in 2007 for his work in dynamics and the French "La Recherche Prize" in 2014 for his works on artificial electric sense.

<figure><img src="image_15.png" /></figure>

**Alessandro De Luca** (Life Fellow, IEEE) received the Ph.D. degree in systems engineering from the University of Rome, Rome, Italy, in 1987. Since 2000, he is currently a Professor of Robotics and Automation with the Sapienza University of Rome, Rome, Italy. His research interests include modeling, motion planning, and control of robotic systems (flexible manipulators, kinematically redundant arms, underactuated robots, and wheeled mobile robots), as well as physical human-robot interaction.

Dr. De Luca was the first Editor-in-Chief for the IEEE TRANSACTIONS ON ROBOTICS from 2004 to 2008, the RAS Vice-President for Publication Activities during 2012-2013, the General Chair of ICRA 2007, and the Program Chair of ICRA 2016. He is a recipient of three conference awards (Best Paper at ICRA 1998 and BioRob 2012, Best Application Paper at IROS 2008), the Helmholtz Humboldt Research Award in 2005, the IEEE-RAS Distinguished Service Award in 2009, and the IEEE George Saridis Leadership Award in Robotics and Automation in 2019. He was the Scientific Coordinator of the FP7 project SAPHARI—Safe and Autonomous Physical Human-Aware Robot Interaction from 2011 to 2015.

<figure><img src="image_16.png" /></figure>

**Federico Renda** (Member, IEEE) received the B.Sc. and M.Sc. degrees in biomedical engineering from the University of Pisa, Pisa, Italy, in 2007 and 2009, respectively, and the Ph.D. degree in biorobotics from Biorobotics Institute, Scuola Superiore Sant'Anna, Pisa, in 2014.

He is currently an Associate Professor with the Department of Mechanical and Nuclear Engineering, Khalifa University, Abu Dhabi, United Arab Emirates. His research interests include dynamic modeling and control of soft and continuum robots using

principles of geometric mechanics.

Open Access funding provided by Università degli Studi di Roma "La Sapienza" 2' within the CRUI CARE Agreement