# Cosserat Rod Modeling of Continuum Robots from Newtonian and Lagrangian Perspectives

Matthias Tummers, Vincent Lebastard, Frédéric Boyer, Jocelyne Troccaz, Benoit Rosa, M. Taha Chikhaoui

► To cite this version:

Matthias Tummers, Vincent Lebastard, Frédéric Boyer, Jocelyne Troccaz, Benoit Rosa, et al. Cosserat Rod Modeling of Continuum Robots from Newtonian and Lagrangian Perspectives. IEEE Transactions on Robotics, 2023, 39 (3), pp.2360-2378. 10.1109/TRO.2023.3238171. hal-03935561

HAL Id: hal-03935561

https://hal.science/hal-03935561

Submitted on 18 Jan 2023

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_1.png" /></figure>
# Cosserat Rod Modeling of Continuum Robots from Newtonian and Lagrangian Perspectives

Matthias Tummers, Vincent Lebastard, Member, IEEE, Frédéric Boyer, Jocelyne Troccaz, Fellow, IEEE, Benoît Rosa, Member, IEEE, and M. Taha Chikhaoui, Member, IEEE

***Abstract***—Cosserat rod theory proved efficient modeling performances in robotics, especially in the context of continuum robots, in the past decade. The implementation of such theory is far from being unique and straightforward. We consider the illustrative example of multi-segment, general routing tendon actuated continuum robots in their nominal static operating regime. This paper details two main approaches based on Cosserat rod modeling, namely the Newtonian and Lagrangian approaches. We provide a walk-through guide regarding theoretical derivations and numerical implementation of both approaches, together with a proof of equivalence. This comparative study is supplemented with novel contributions and extensions of each approach and in-depth discussion of their performances and applicability, as well as highlighting their special features.

***Index Terms***—Modeling, Control, and Learning for Soft Robots, Flexible Robots, Tendon/Wire Mechanism, Cosserat Beams.

## I. INTRODUCTION

CONTINUUM robots are slender actuated flexible structures that can be compared to elephant trunks, octopus arms or squid tentacles. The main characteristics of continuum robots are their ability to conform to curvilinear paths and inherent compliance. In addition, their designs often offer ease of miniaturization and enhanced dexterity. These advantages enable continuum robots to have safe and soft interactions with unstructured environments, to evolve in confined spaces or avoid obstacles, and offer tolerance for geometric variations in grasped objects. This type of robots is suited for a myriad of applications, ranging from nuclear decontamination to medical applications [1], [2].

Modeling continuum robots is a challenging task due to the lack of discrete joints and their compliant nature, thus their deformation when in interaction with the environment. As such, numerous modeling approaches have been developed [1]–[4]. Two main modeling categories can be identified: geometry based models and mechanical models. Geometry based models [5], [6], including the well known constant curvature approach [7]–[11], have the advantage to lead to analytical formulations and are thus suitable for real-time

inverse kinematics and control. On the downside, these models are less accurate than their mechanical counterpart, are subject to error propagation, and cannot take external efforts into account [12]. Mechanical models, on the other hand, take the elastic behavior of the constitutive materials of the robots into account and solve the problem over boundary conditions (BCs) that can involve external loading of the structures.

Some mechanical models for continuum robots employ the Euler-Bernoulli [13]–[15] or Timoshenko beam theories that are geometrically linear and can model only small rotations of the beam’s cross-section. Thus, the community started exploiting the more general Cosserat rod theory [16]–[19]. The approach was applied to various continuum manipulators as concentric tube continuum robots [20], [21], tendon actuated continuum robots (TACRs) [22]–[25] and multi-backbone continuum robots [26]–[28]. The method has been studied extensively for modeling the statics and dynamics [29]–[31] and has become the standard over the past ten years in the robotics community. In a majority of cases, the authors propose to solve the boundary value problem (BVP) of the Cosserat rod with a shooting based approach. The limitations of the approach include the dependence of the convergence on the initial guess and the high dimension state vectors for complex robot structures [26], [27]. In this paper, this approach will be called the Newtonian approach, as it relies solely on Newton’s laws, and in particular on the action-reaction principle (Newton’s second law) to include actuation in the robot model.

Lately, another approach for solving the same mechanical problem was proposed based on a Lagrangian viewpoint [32]–[34]. In essence, the Lagrangian approach for mechanics is based on two key ingredients: first the concept of configuration space, second a variational principle that allows the static or dynamic balance on this space to be derived [35]. In contrast to Newtonian mechanics, this more abstract point of view allows to directly derive a minimal set of equations without the need for isolating subsystems and using the action-reaction principle. This advantage is particularly useful for complex continuum robots involving multiple rods in interaction with each other. In the case here considered, Lagrangian mechanics is used as a systematic reduction process of the Cosserat rod model. This reduction consists in parameterizing the configuration of a rod by the coefficients of its strain fields on a truncated functional basis according to the Ritz method [33]. It allows the application of efficient numerical approximation

<footnote>This work was supported by grants ANR-11-LABX-0004-01, ANR-19-P31A-0003, ANR-20-CE33-0001, ANR-10-IAHU-0002, and ANR-18-CE19-0012. M. Tummers, M.T. Chikhaoui, and J. Troccaz are with Univ. Grenoble Alpes, CNRS, UMR 5525, VetAgro Sup, Grenoble INP, TIMC, 38000 Grenoble, France matthias.tummers@univ-grenoble-alpes.fr. F. Boyer and V. Lebastard are with LS2N Laboratory, Institut Mines Telecom Atlantique, 44307 Nantes, France. B. Rosa is with ICube Laboratory, CNRS-Strasbourg University, 67400 Illkirch-Graffenstaden, France.
Digital Object Identifier 10.1109/TRO.2023.3238171</footnote>
methods based on variational formulations, both in statics and in dynamics [36]. On the other hand, since the residual vector is based on the weak form of the virtual works, its derivation is less straightforward than that of the Newtonian approach. Moreover, the dimension of the residual vector is generally larger, and increases with the number of segments.

This paper focuses on the modeling of continuum robots and chooses the use-case of TACRs in their nominal static operating regime (used in the representative medical applications of continuum robots), to illustrate theoretical models and numerical simulations. TACRs consist of a central elastic backbone, which is pulled upon by a set of tendons (i.e. cables) attached at its end. The tendons are routed through holes in disks, distributed along the length of the backbone. By routing tendons at different angles around the backbone and pulling on them, bending can be achieved in a 3D workspace.

The main motivation of this paper is to resolve the confusion of a certain number of concepts, shared or split between the Newtonian and the Lagrangian approaches, which arise due to the different viewpoints they are based on. In addition, the two approaches, coming from different communities, do not share a common set of notations or mathematical framework, which does not ease the comparison. Moreover, the continuum robotics community, being less familiar with Lagrangian mechanics, may face difficulties to make informed choices as to which approach should be used for each specific case. This paper provides an exhaustive derivation of both approaches in a common mathematical framework in order to clarify the links between them and objectively compare the strengths and weaknesses of each.

Beyond presenting, in a single structured picture, the two approaches that are currently dominant in the community, this article presents several specific contributions to each of these approaches. In the Newtonian approach, beyond [22], the case of a multi-segment TACR with slope discontinuities in the tendon routings between segments is addressed. For the Lagrangian approach, besides a more in-depth derivation of the static Lagrangian model through a canonical application of the principle of virtual work, new numerical perspectives are proposed. In [32], the numerical resolution is performed through explicit time integration of an overdamped equivalent system. In contrast, we here introduce a linearization of the nonlinear static balance equations, which enables to use Newton-Raphson's method. Another novelty of the present work is that efficient spectral methods [37] are used to calculate the residual vector and the Jacobian matrix thanks to a new BVP, called the inverse kineto-static BVP, and its tangent BVP.

Based on the comparison of the two methods, we provide a mathematical proof of the equivalence of the two approaches. This equivalence is backed by an extensive set of simulations, which compare the Lagrangian approach to the Newtonian approach. The latter was previously experimentally validated in [22] for TACRs. The simulations also put in light specifics of each approach owing to differences in numerical implementation, allowing the reader to make an advised choice. To ease prototyping of TACR applications, the associated code can be found at https://github.com/TIMClab-CAMI/Cosserat-

Rod-Modeling-of-Tendon-Actuated-Continuum-Robots.

Section II explains the reasons behind the differences of notations that can be found in the literature and the correspondence between them. In Section III, we establish the common mathematical framework as a building block for the derivation of both approaches in Sections IV and V. In these sections, the necessary keys for easily bridging with previously published work are given. Section VI provides the mathematical proof of equivalence between the two approaches. In Section VIII, both approaches are compared in terms of numerical performance using a common simulation framework. In Section IX, we discuss the implementation considerations and the simulation results, as well as general performance considerations of both approaches in different use-cases.

## II. NOTATIONS

When the Newtonian approach and the Lagrangian approach first occurred in the literature, two clear different systems of notations were recognizable. The symbols in the Newtonian approach were chosen to match as closely as possible both the notation proposed by S. S. Antman [38] and the existing literature on the modeling of continuum robots [19]. On the other hand, the Lagrangian approach privileged notations from geometric mechanics [39] and its applications to rigid multi-body systems and continuous media, including fluid mechanics [40]. This choice also eases the application of developments on manipulators to other contexts as e.g. eel-like robots. This paper purposely uses Lagrangian notations. This choice aims at better introducing the Lagrangian approach by familiarizing the reader with these notations throughout the developments of the Newtonian approach. A dictionary of the two systems of notations is presented in Table I.

Table I highlights some key differences that can be summarized as follows. In the Newtonian approach, vector fields are mostly expressed in $\mathbb{R}^3$ whereas the Lagrangian approach makes extensive use of wrenches and twists, expressed in $\mathbb{R}^6$. In the literature on the Newtonian approach, these $\mathbb{R}^3$ vectors are represented by bold symbols to differentiate them from scalars and matrices. As many developments of the Lagrangian approach are conducted in other (vector) spaces than $\mathbb{R}^3$, this practice is less applicable in this context, and will not be used in this article. In both approaches all vector fields in $\mathbb{R}^3$ can be expressed either in the inertial frame or in the cross-sectional frames, which can sometimes lead to confusion<fnref n="1" />. To avoid ambiguity, this paper adopts the convention of Juan Carlos Simo of noting a vector of $\mathbb{R}^3$ with a lowercase letter (preferably latin) when expressed in inertial frame, and uppercase, when expressed in cross-sectional frame. The only exception to this rule is the position vector field $r$, which is always expressed in the inertial frame, $R$ being used for the orientation. Vector fields in $\mathbb{R}^6$ are not concerned by this rule, since they represent vectors in the Lie algebra $se(3)$, denoted with Greek letters as it is usually the case in geometric mechanics on $SE(3)$ [39].

<footnote n="1">¹The term 'cross-sectional frame' will be used, which is equivalent to the terms 'body frame' and 'material frame'. Correspondingly, the reference frame is named the 'inertial frame' and is equivalent to the terms 'fixed frame' and 'spatial frame'.</footnote>
<table><caption>Table I
COMPARISON OF THE NOTATION CONVENTIONS OF THE LITERATURE.</caption><thead><tr><th></th><th>Newtonian approach</th><th>Lagrangian approach (+ this paper)</th><th>Dimensions</th><th>Units</th></tr></thead><tbody><tr><td>Reference length parameter</td><td>s &isin; [0, l]</td><td>X &isin; [0, l]</td><td>&#8477;</td><td>m</td></tr><tr><td>Position</td><td><b>p</b>(s)</td><td>r(X)</td><td>&#8477;<sup>3</sup></td><td>m</td></tr><tr><td>Orientation</td><td>R(s), (z-axis aligned)</td><td>R(X), (x-axis aligned)</td><td>&#8477;<sup>3</sup>&times;3</td><td></td></tr><tr><td>Cross-sectional frame</td><td>g(s) = (<b>p</b>, R)(s)</td><td>g(X) = (R, r)(X)</td><td>&#8477;<sup>4</sup>&times;4</td><td></td></tr><tr><td>Linear rate of change (shear &amp; extension)</td><td><b>v</b>(s)</td><td>&Gamma;(X)</td><td>&#8477;<sup>3</sup></td><td></td></tr><tr><td>Angular rate of change (bending &amp; torsion)</td><td><b>u</b>(s)</td><td>K(X)</td><td>&#8477;<sup>3</sup></td><td>m<sup>-1</sup></td></tr><tr><td>Space-rate twist</td><td>-</td><td>&xi; = (K<sup>T</sup>, &Gamma;<sup>T</sup>)<sup>T</sup></td><td>&#8477;<sup>6</sup></td><td>(m<sup>-1</sup>, -)</td></tr><tr><td>Initial value or prior to deformation*</td><td>&bull;*</td><td>&bull;0</td><td></td><td></td></tr><tr><td>Total length of the robot</td><td>l</td><td>l</td><td>&#8477;</td><td>m</td></tr><tr><td>Strain</td><td>-</td><td>&epsilon; = &xi; - &xi;<sub>0</sub></td><td>&#8477;<sup>6</sup></td><td>(m<sup>-1</sup>, -)</td></tr><tr><td>Related to the i<sup>th</sup> tendon</td><td>&bull;<sub>i</sub></td><td>&bull;<sub>i</sub></td><td></td><td></td></tr><tr><td>Number of tendons</td><td>n</td><td>m</td><td>&#8477;</td><td></td></tr><tr><td>Length of the i<sup>th</sup> tendon</td><td>l<sub>i</sub></td><td>l<sub>i</sub></td><td>&#8477;</td><td>m</td></tr><tr><td>Related to the j<sup>th</sup> segment</td><td>&bull;<sub>j</sub></td><td>&bull;<sub>j</sub></td><td></td><td></td></tr><tr><td>Number of segments</td><td>m</td><td>&chi;</td><td>&#8477;</td><td></td></tr><tr><td>Length of the j<sup>th</sup> segment</td><td>l<sub>j</sub></td><td>l<sub>j</sub></td><td>&#8477;</td><td>m</td></tr><tr><td>Cross-sectional frame position of the i<sup>th</sup> tendon</td><td><b>r</b><sub>i</sub>(s) = [x<sub>i</sub>(s) y<sub>i</sub>(s) 0]<sup>T</sup></td><td>D<sub>i</sub>(X) = (0, D<sub>i,Y</sub>, D<sub>i,Z</sub>)(X)</td><td>&#8477;<sup>3</sup></td><td>m</td></tr><tr><td>Inertial frame tendon position</td><td><b>p</b><sub>i</sub>(s) = <b>p</b> + R<b>r</b><sub>i</sub></td><td>r<sub>i</sub>(X) = r + RD<sub>i</sub></td><td>&#8477;<sup>3</sup></td><td>m</td></tr><tr><td>Tangent to the path of the i<sup>th</sup> tendon</td><td><b>p</b>&#775;<sub>i</sub>/||<b>p</b>&#775;<sub>i</sub>||</td><td>t<sub>i</sub> = r'<sub>i</sub>/||r'<sub>i</sub>||</td><td>&#8477;</td><td></td></tr><tr><td>Tension in the i<sup>th</sup> tendon</td><td>&tau;<sub>i</sub></td><td>&tau;<sub>i</sub></td><td>&#8477;</td><td>N</td></tr><tr><td>Vector of tensions</td><td>-</td><td>&tau; = (&tau;<sub>1</sub> .. &tau;<sub>m</sub>)<sup>T</sup></td><td>&#8477;<sup>m</sup></td><td>N</td></tr><tr><td>Stress (bending &amp; torsion)</td><td><b>m</b>(s)</td><td>C(X)</td><td>&#8477;<sup>3</sup></td><td>Nm</td></tr><tr><td>Stress (shear &amp; extension)</td><td><b>n</b>(s)</td><td>N(X)</td><td>&#8477;<sup>3</sup></td><td>N</td></tr><tr><td>Stress (wrench)</td><td>-</td><td>&Lambda;(X) = (C(X)<sup>T</sup>, N(X)<sup>T</sup>)<sup>T</sup></td><td>&#8477;<sup>6</sup></td><td>(Nm, N)</td></tr><tr><td>Angular and linear stiffness matrices</td><td>K<sub>bt</sub>(s), K<sub>se</sub>(s)</td><td>H<sub>a</sub>(X), H<sub>l</sub>(X)</td><td>&#8477;<sup>3</sup>&times;3</td><td>N, Nm<sup>2</sup></td></tr><tr><td>Hookean stiffness matrix</td><td></td><td>H(X) = diag(H<sub>a</sub>(X), H<sub>l</sub>(X))</td><td>&#8477;<sup>6</sup>&times;6</td><td>(N, Nm<sup>2</sup>)</td></tr><tr><td>Generalized stiffness</td><td>-</td><td>K<sub>&epsilon;&epsilon;</sub></td><td>&#8477;<sup>k&times;k</sup></td><td>**</td></tr><tr><td>Total external distributed couples</td><td>l(s)</td><td>&#8477;<sub>c</sub>(X), &#8477;<sub>c</sub>(X)</td><td>&#8477;<sup>3</sup></td><td>Nm/m</td></tr><tr><td>Total external distributed forces</td><td><b>f</b>(s)</td><td>&#8477;<sub>n</sub>(X), &#8477;<sub>n</sub>(X)</td><td>&#8477;<sup>3</sup></td><td>N/m</td></tr><tr><td>Total external distributed wrench</td><td>-</td><td>&#8477;<sub>F</sub>(X) = (&#8477;<sub>c</sub>(X)<sup>T</sup>, &#8477;<sub>n</sub>(X)<sup>T</sup>)<sup>T</sup></td><td>&#8477;<sup>6</sup></td><td>(Nm/m, N/m)</td></tr><tr><td>External tip couples</td><td>-</td><td>c<sub>+</sub>, C<sub>+</sub></td><td>&#8477;<sup>3</sup></td><td>Nm</td></tr><tr><td>External tip forces</td><td>-</td><td>n<sub>+</sub>, N<sub>+</sub></td><td>&#8477;<sup>3</sup></td><td>N</td></tr><tr><td>External tip wrench</td><td>-</td><td>F<sub>+</sub> = (C<sub>+</sub><sup>T</sup>, N<sub>+</sub><sup>T</sup>)<sup>T</sup></td><td>&#8477;<sup>6</sup></td><td>(Nm, N)</td></tr><tr><td>Relative to external loads (not tendons)</td><td>-</td><td>&bull;<sub>ext</sub></td><td></td><td></td></tr><tr><td>External distributed forces (not tendons)</td><td><b>f</b><sub>e</sub>(s)</td><td>&#8477;<sub>ext</sub>(X)</td><td>&#8477;<sup>3</sup></td><td>N/m</td></tr><tr><td>External distributed couples (not tendons)</td><td>l<sub>e</sub>(s)</td><td>&#8477;<sub>ext</sub>(X)</td><td>&#8477;<sup>3</sup></td><td>Nm/m</td></tr><tr><td>Relative to rod elastics</td><td>-</td><td>&bull;<sub>rod</sub></td><td></td><td></td></tr><tr><td>Relative to actuation by the tendons</td><td>-</td><td>&bull;<sub>act</sub></td><td></td><td></td></tr><tr><td>Distributed forces applied by the tendons</td><td><b>f</b><sub>t</sub>(s)</td><td>&#8477;<sub>act</sub>(X)</td><td>&#8477;<sup>3</sup></td><td>N/m</td></tr><tr><td>Distributed couples applied by the tendons</td><td>l<sub>t</sub>(s)</td><td>&#8477;<sub>act</sub>(X)</td><td>&#8477;<sup>3</sup></td><td>Nm/m</td></tr><tr><td>Distributed forces applied to the i<sup>th</sup> tendon</td><td><b>f</b><sub>i</sub>(s)</td><td>&#8477;<sub>i</sub>(X)</td><td>&#8477;<sup>3</sup></td><td>N/m</td></tr><tr><td>Stress of actuation</td><td>-</td><td>&Lambda;<sub>act</sub>(X)</td><td>&#8477;<sup>6</sup></td><td>(Nm, N)</td></tr><tr><td>Generalized external forces</td><td>-</td><td>Q<sub>ext</sub></td><td>&#8477;<sup>k</sup></td><td>**</td></tr><tr><td>Generalized restoring forces</td><td>-</td><td>Q<sub>rod</sub></td><td>&#8477;<sup>k</sup></td><td>**</td></tr><tr><td>Generalized forces of actuation</td><td>-</td><td>Q<sub>act</sub></td><td>&#8477;<sup>k</sup></td><td>**</td></tr><tr><td>Matrix of actuation</td><td>-</td><td>L</td><td>&#8477;<sup>k&times;m</sup></td><td>**</td></tr><tr><td>Number of shape functions</td><td>-</td><td>k</td><td>&#8477;</td><td></td></tr><tr><td>Shape functions</td><td>-</td><td>&Phi;(X) = (&Phi;<sub>1</sub> .. &Phi;<sub>k</sub>)<sup>T</sup></td><td>&#8477;<sup>6&times;k</sup></td><td>**</td></tr><tr><td>Generalized strain coordinates</td><td>-</td><td>q = (q<sub>1</sub> .. q<sub>k</sub>)<sup>T</sup></td><td>&#8477;<sup>k</sup></td><td>**</td></tr><tr><td>First and second derivative wrt. s or X</td><td>&bull;, &bull;</td><td>&bull;,&bull;''</td><td></td><td></td></tr><tr><td>Residual vector</td><td>-</td><td>&#8477;</td><td></td><td></td></tr></tbody></table>

<footnote>* The subscript ●<sub>o</sub> from the literature on the Lagrangian approach was, in this paper, replaced with a subscript ●<sub>0</sub> to prevent confusion.</footnote>

<footnote>** The units of these quantities depend on the choice of the shape functions.</footnote>
Another difference concerns the partial derivation with respect to the rod arc length variable in its reference (stress-less) configuration. It is denoted ● in the Newtonian approach and ●' in the Lagrangian approach. Related to this spatial differentiation are the densities of external forces, couples, and wrenches per unit arc length that are distinguished from the ordinary forces, couples, and wrenches by an overbar (e.g. $\bar{n}$ or $\bar{F}$ in Table I).

Using the symbols $q$ and $Q$, for the generalized strain coordinates, and the corresponding generalized forces, respectively, is a deliberate choice inherited from the standards of Lagrangian mechanics [35]. In the context of continuum robotics, this choice reflects the fact that the strain coordinates of a continuum robot are the distributed counterpart of the localized joint coordinates of a rigid robot.

A few notational conventions related to Lie group theory of SE(3) are presented below. In particular, the two approaches share common convention notations related to the ● operator, defined as follows. If $W = (W_1, W_2, W_3)^T$ is a vector of $\mathbb{R}^3$, then $\hat{W}$ denotes the skew-symmetric matrix of $\mathbb{R}^{3 \times 3}$ such that

$$ \hat{W} = \begin{pmatrix} 0 & -W_3 & W_2 \\ W_3 & 0 & -W_1 \\ -W_2 & W_1 & 0 \end{pmatrix}. $$

The reciprocal is denoted ●$^V$ such that $(\hat{W})^V = W$. Similarly, if $W = (W_1, W_2, W_3)^T$ and $V = (V_1, V_2, V_3)^T$ are two vectors of $\mathbb{R}^3$, with $\nu = (W^T, V^T)^T$ a vector of $\mathbb{R}^6$, then $\hat{\nu}$ denotes the matrix of $\mathbb{R}^{4 \times 4}$ as

$$ \hat{\nu} = \begin{pmatrix} \hat{W} & V \\ 0_{1 \times 3} & 0 \end{pmatrix}, $$

where $\hat{W}$ is the skew-symmetric matrix associated to $W \in \mathbb{R}^3$. Finally, the Lagrangian approach uses some of the standard notations of geometric mechanics on the Lie group SE(3), namely the adjoint maps $ad$ and $Ad$ defined for any field of twist $\nu$ and homogeneous transformation $g$,

$$ \nu = \begin{pmatrix} W \\ V \end{pmatrix} \quad \text{and} \quad g = \begin{pmatrix} R & r \\ 0_{1 \times 3} & 1 \end{pmatrix}, $$

by the two matrices of $\mathbb{R}^{6 \times 6}$

$$ ad_{\nu} = \begin{pmatrix} \hat{W} & 0_{3 \times 3} \\ \hat{V} & \hat{W} \end{pmatrix} \quad \text{and} \quad Ad_{g} = \begin{pmatrix} R & 0_{3 \times 3} \\ \hat{r}R & R \end{pmatrix}. $$

Throughout the paper, $Ad_g$ will be used to pass a twist from one frame to another, the two being separated by a pose $g$. And $ad_\nu$ will be used to model the inertial derivative of the cross-sectional basis as it moves along its centerline with a twist $\nu$ expressed in the moving frame.

### III. KINETO-STATIC MODEL OF COSSERAT RODS

Cosserat rod theory is a one-dimensional theory of continuous media used to describe slender bodies subject to finite deformations. In contrast to vibration theory or strength of materials, Cosserat rod theory is geometrically exact in the sense that it is not based on any approximation of small displacements or slopes. The equations of this theory are derived below as a basic element for the modeling of various continuum robots.

<figure><img src="image_2.png" /><figcaption>Figure 1. Free body diagram of a piece of rod. The yellow line represents an arbitrary finite piece of a rod from $X_1$ to $X_2$. The dash-dot line represents the centerline of the rod along which the reference length parameter $X$ evolves. Some rod sections are represented (in green) to each of which a cross-sectional frame is attached (in red). The transformation $g(X) = (R, r)(X)$ shows the link between the inertial frame (in blue) and the cross-sectional frame. The rod is at static equilibrium between the external forces $\bar{n}(X)$ and couples $\bar{c}(X)$, and its internal restoring forces $n(X)$ and couples $c(X)$.</figcaption></figure>

#### A. Kinematics

Let us consider a single rod in space, of length $l$. In Cosserat theory, the rod is modeled by a continuous stacking of rigid cross-sections labeled by a continuous index that is here chosen to be the arc length, noted $X \in [0, l]$, counted along the rod in a reference stress-less configuration. To each cross-section that is assumed to be circular, is attached a cross-sectional frame whose origin coincides with the cross-section center. The position and orientation of rigid cross-sections are described, respectively, by an $X$-parameterized curve in space $r: X \in [0, l] \mapsto r(X) \in \mathbb{R}^3$, and another in the Lie group of rotations $R: X \in [0, l] \mapsto R(X) \in \text{SO}(3)$. Together, $r$ and $R$ define a field of homogeneous transformation that describes the entire rod

$$ g: X \in [0, l] \mapsto g(X) = \begin{pmatrix} R(X) & r(X) \\ 0_{1 \times 3} & 1 \end{pmatrix} \in \text{SE}(3). \quad (1) $$

The matrix $g(X)$ represents the position and orientation of the $X$-cross-section of the Cosserat rod or the pose of the $X$-‘cross-sectional frame’ in the inertial frame (see Fig. 1). At this point, the configuration space of a rod clamped in a fixed basis is defined by the set of all possible fields of pose $g$, with $g(0) = \mathbf{1}_{4 \times 4}$, i.e.:

$$ C_1 = \{g: X \in [0, l] \mapsto g(X) \in \text{SE}(3), \text{ with: } g(0) = \mathbf{1}_{4 \times 4}\}. \quad (2) $$

The space-variations of position and orientation along the length are modeled by the linear and angular rates of change, respectively. In the cross-sectional frame, these two rates are defined by the two fields $\Gamma$ and $K$ of $\mathbb{R}^3$ such that

$$ r' = R\Gamma, \quad R' = R\hat{K}, \quad (3) $$

where the ●' symbol denotes differentiation with respect to $X$. When combined, the fields of linear and angular space-rate $\Gamma$ and $K$, define a unique field of twist $\xi$ such that

$$ g' = g\xi \quad (4) $$

is equivalent to the two relations of (3).
## B. Static Model of a Cosserat Rod

As shown in Fig. 1, the rod is at static equilibrium under the effect of some external force $\bar{n}$ and couple $\bar{c}$ density fields defined over $[0, l]$, and some tip force $n_+$ and couple $c_+$. Thus, the rod generates some internal restoring forces $n$ and couples $c$ over $[0, l]$ that fulfill Newton's balance of forces and couples in the inertial frame

$$n' + \bar{n} = 0_{3 \times 1}, \quad c' + r' \times n + \bar{c} = 0_{3 \times 1}, \quad \forall X \in [0, l]. \quad (5)$$

The rod is assumed to be clamped at $X = 0$, and these two sets of ordinary differential equations (ODEs) need to be supplemented by the BCs

$$r(0) = 0_{3 \times 1}, \quad R(0) = 1_{3 \times 3}, \quad n(l) = n_+, \quad c(l) = c_+. \quad (6)$$

Note here that $n(X)$ and $c(X)$ are transmitted across the X-cross-section. According to the conventions of continuum mechanics, if $X$ is oriented positively from the proximal to the distal end of the rod, $n(X)$ and $c(X)$ are counted positively when exerted by the distal piece of rod (i.e. after arc-length $X$) onto the proximal piece (i.e. before arc-length $X$). This first convention is in contrast to the convention of rigid multi-body dynamics where the inter-body wrenches are counted positively when applied from the basis to the tip. In what follows, the stress-less reference configuration of any kinematic variable is noted with index $\bullet_0$. As such, in the small strain regime, the two strain fields of the rod, related to shear and extension, and bending and torsion, respectively, are defined as $\Gamma - \Gamma_0$ and $K - K_0$. Moreover, assuming linear material elasticity, the internal forces $n$ and couples $c$ are related to $\Gamma$ and $K$ through the constitutive laws

$$n = R\mathcal{H}_l(\Gamma - \Gamma_0), \quad c = R\mathcal{H}_a(K - K_0), \quad (7)$$

where, $\forall X \in [0, l]$, $\mathcal{H}_l(X)$ and $\mathcal{H}_a(X)$ are the arc length parameterized stiffness matrices for shear and extension, and bending and torsion, respectively. (5) to (7) define a closed static formulation of a clamped-free Cosserat rod in the inertial frame. Alternatively, the same equations can be expressed in the cross-sectional frame. To derive this other version, all the vectors of the above formulation are shifted from the inertial frame to the mobile cross-sectional ones, using the notation conventions $r' = R\Gamma$, $n = RN$, $c = RC$, $\bar{n} = R\bar{N}$, and $\bar{c} = R\bar{C}$, as follows:

$$(RN)' + R\bar{N} = 0_{3 \times 1}, \quad (8) \\ (RC)' + (R\Gamma) \times (RN) + R\bar{C} = 0_{3 \times 1}.$$

Then, via composition rules of derivatives along with the relation $R' = R\hat{K}$

$$RN' + R\hat{K}N + R\bar{N} = 0_{3 \times 1}, \quad (9) \\ RC' + R\hat{K}C + (R\Gamma) \times (RN) + R\bar{C} = 0_{3 \times 1},$$

which, being true for any $R$ and owing to the definition of the hat operator, provides

$$N' + K \times N + \bar{N} = 0_{3 \times 1}, \quad (10) \\ C' + K \times C + \Gamma \times N + \bar{C} = 0_{3 \times 1}.$$

Once supplemented with the BCs

$$r(0) = 0_{3 \times 1}, \quad N(l) = R(l)^T n_+, \quad (11) \\ R(0) = 1_{3 \times 3}, \quad C(l) = R(l)^T c_+,$$

which are deduced from their inertial frame version (6), equations (10) stand for the ODEs of a Cosserat rod in the cross-sectional frame. Following the convention of [41], all pairs of angular and linear vectors gather in $\mathbb{R}^6$ vectors

$$\begin{pmatrix} C \\ N \end{pmatrix}' + \begin{pmatrix} K \times C + \Gamma \times N \\ K \times N \end{pmatrix} + \begin{pmatrix} \bar{C} \\ \bar{N} \end{pmatrix} = \begin{pmatrix} 0_{3 \times 1} \\ 0_{3 \times 1} \end{pmatrix}. \quad (12)$$

Let us introduce the notations $\Lambda = (C^T, N^T)^T$ for the $\mathbb{R}^6$ field of stress wrench, $\bar{F} = (\bar{C}^T, \bar{N}^T)^T$ for the $\mathbb{R}^6$ field of external loads, and $F_+ = (C_+^T, N_+^T)^T$, with $C_+ = R(l)^T c_+$ and $N_+ = R(l)^T n_+$ for the wrench of external tip loads. These notations enable to rewrite the two ODEs (12) in the more compact (geometric) form as a balance of wrenches expressed in the cross-sectional frame

$$\Lambda' - ad_{\xi}^{T} \Lambda + \bar{F} = 0_{6 \times 1}. \quad (13)$$

To which one needs to add the BCs

$$g(0) = 1_{4 \times 4}, \quad \Lambda(l) = F_+. \quad (14)$$

Similarly to the twist form of (4), (7) can be rewritten as

$$\Lambda = \mathcal{H}\epsilon = \mathcal{H}(\xi - \xi_0), \quad (15)$$

where $\epsilon = \xi - \xi_0 = ((K - K_0)^T, (\Gamma - \Gamma_0)^T)^T$ is the $\mathbb{R}^6$ field of strain and $\mathcal{H} = \text{diag}(\mathcal{H}_a, \mathcal{H}_l)$ is the $\mathbb{R}^{6 \times 6}$ Hookean stiffness matrix.

## C. Reshaping the Cosserat Model as a BVP

Whether in the Newtonian or in the Lagrangian approach, the notion of BVP plays an essential role in the modeling of continuum robots. In short, a BVP is a system of ODEs whose solutions must satisfy BCs that partially determine the state variables at the boundaries. In the present context, such a system can be set in state space form

$$x' = f(x), \quad h_-(x(0)) = 0, \quad h_+(x(l)) = 0, \quad (16)$$

where $x$ is the vector of state variables, while the two functions $h_+$ of the state vector $x$, at $X = 0$ and $l$, fix the BCs at the two ends of the rod, with $\text{dim}(h_-) + \text{dim}(h_+) = \text{dim}(x)$. Based on this definition, it is straightforward to show that using the constitutive laws (7) to remove the internal forces and couples in (5) and (6) and gathering the resulting equations with (3), yields a closed formulation describing the statics of an elastic rod

$$\begin{pmatrix} r \\ R \\ \Gamma \\ K \end{pmatrix}' = \begin{pmatrix} R\Gamma \\ R\hat{K} \\ \Gamma_0' - \mathcal{H}_l^{-1} \left[ \begin{pmatrix} \hat{K}\mathcal{H}_l + \mathcal{H}_l' \end{pmatrix} (\Gamma - \Gamma_0) + R^T \bar{n} \right] \\ K_0' - \mathcal{H}_a^{-1} \left[ \begin{pmatrix} \hat{K}\mathcal{H}_a + \mathcal{H}_a' \end{pmatrix} (K - K_0) + \hat{\Gamma}\mathcal{H}_l (\Gamma - \Gamma_0) + R^T \bar{c} \right] \end{pmatrix},$$

$$r(0) = 0_{3 \times 1}, \quad \Gamma(l) = \Gamma_0 + \mathcal{H}_l^{-1} R(l)^T n_+, \\ R(0) = 1_{3 \times 3}, \quad K(l) = K_0 + \mathcal{H}_a^{-1} R(l)^T c_+. \quad (17)$$
<figure><img src="image_3.png" /><figcaption>Figure 2. Synthesis of the steps in both approaches side by side. The corresponding steps meet horizontally across the columns. The blue and red areas are relative to the Newtonian and the Lagrangian approaches, respectively. The yellow area is common to both (up to Section III). The corresponding equation numbers are provided for each step. The abbreviations act., ext., and int. stand for actuation, external, and internal, respectively.</figcaption></figure>

(17) defines a BVP in the explicit state form (16), with $x = (r, R, \Gamma, K)$ being the state variables. Note that using constitutive relations in alternate ways provides another equivalent BVP where strain variables $(\Gamma, K)$ are replaced by stress ones $(n, c)$, while all the strain and stress vectors can be expressed in the cross-sectional or inertial frames.

### D. From Cosserat Rods to TACRs

The following two sections present two different approaches that use the above derived equations to model TACRs. The main steps of each of the approaches are synthesized in Fig. 2 and the accompanying video to help the reader follow through the developments and highlight the correspondence of the steps between approaches. Both approaches assume the same set of model simplifications. The tendons are inextensible, experience no friction, and are routed through continuous routing paths rather than through discrete disks. These are all reasonable assumptions for small curvatures and a large amount of routing disks, but may introduce significant discrepancies when this is not the case. Also, in addition to the assumption of material linear elasticity (see (7)), the considered robots have constant material properties along their length (i.e. $H' = 0_{6 \times 6}$).

## IV. NEWTONIAN MODEL OF TACRS

In the Newtonian approach, the model is deduced from Newton's laws, i.e. by isolating the backbone and each of the tendons separately, and by using the action-reaction principle to remove the interaction forces between them. To detail this process, let us first remark that, in a TACR, the backbone of the robot can be modeled by the BVP (17). Sections IV-A and IV-B describe how this BVP is complemented with a model of tendon actuation.

### A. Model of Tendon Actuation

To introduce the model of tendon actuation in (17), the distributed loads are spelled out as

$$ \bar{n} = \bar{n}_{ext} + \bar{n}_{act}, \quad \bar{c} = \bar{c}_{ext} + \bar{c}_{act}, \qquad (18) $$

where $\bar{n}_{ext}$ and $\bar{c}_{ext}$ model the external loads applied by the environment such as gravity or contacts, while $\bar{n}_{act}$ and $\bar{c}_{act}$ stand for those applied by the tendons. Invoking the action-reaction principle, the distributed force exerted by each tendon $i$ onto the backbone is equal to the opposite of the distributed force exerted by the backbone on this tendon, noted $\bar{n}_i$. Therefore, summing these contributions provides the distributed force of the $m$ tendons on the backbone

$$ \bar{n}_{act} = \sum_{i=1}^{m} -\bar{n}_i. \qquad (19) $$

As a tendon cannot transmit any couples, the distributed moment of the tendons exerted on the backbone $\bar{c}_{act}$ is equal to only the sum of the cross products of each moment arm with each force

$$ \bar{c}_{act} = \sum_{i=1}^{m} (r_i - r) \times (-\bar{n}_i) = - \sum_{i=1}^{m} d_i \times \bar{n}_i, \qquad (20) $$

where $r_i = r + RD_i$ is the position of the $i^{th}$ tendon expressed in the inertial frame, $D_i(X) = (0, D_{i,Y}, D_{i,Z})^T(X)$ the position of the $i^{th}$ tendon expressed in the X-cross-sectional frame, and $d_i = RD_i$. Modeling a tendon as an inextensible degenerate Cosserat rod, with no angular inertia, one can derive the expression for $\bar{n}_i$, by applying the linear static balance (5) with its internal forces $n_i$, to each tendon individually

$$ \bar{n}_i + n'_i = 0_{3 \times 1}. \qquad (21) $$

Again, as a perfectly flexible string cannot support internal couples or shear forces but only a constant tension $\tau_i$ tangent to its path, the field of internal forces along each tendon $i$ reads

$$ n_i = \tau_i t_i, \qquad (22) $$

with $t_i = r'_i/\|r'_i\|$ the unit tangent vector to the path of the tendon. Now, considering that the tension is constant ($\tau'_i = 0_{3 \times 1}$), (21) provides the expression of the external forces of each tendon

$$ \bar{n}_i = -n'_i = -\tau_i t'_i. \qquad (23) $$
## B. Forward kineto-static BVP of a TACR

Following the methodology of [22], (23) is introduced into (19) and (20), and the result into (18) then in (5). These identifications provide the stress balance equations of the TACR actuated by a set of m tendons

$$ 
\begin{aligned}
n' + \sum_{i=1}^{m} \tau_i t'_i + \bar{n}_{ext} &= 0_{3 \times 1}, \\
c' + r' \times n + \sum_{i=1}^{m} d_i \times (\tau_i t'_i) + \bar{c}_{ext} &= 0_{3 \times 1}.
\end{aligned} 
\quad (24) $$

To obtain an explicit form of the ODEs (24) compatible with (16), one needs to express the field of forces exerted by each tendon in terms of the backbone strains (Γ, K) and their space derivatives. To this end, let us first remark that we have

$$ \bar{n}_i = -\tau_i t'_i = -\tau_i \left( \frac{r'_i}{\|r'_i\|} \right)' = \tau_i \frac{\hat{r}_i^2}{\|r'_i\|^3} r''_i, \quad (25) $$

referring to Appendix B of [22] for further details. Next, one needs to expand the first and second derivatives of $\tau_i$. Expressing these quantities in cross-sectional frame coordinates permits a more concise writing, as follows:

$$ 
\begin{aligned}
r'_i &= R\Gamma_i = R(\hat{K}D_i + D'_i + \Gamma), \\
r''_i &= R(\Gamma'_i + \hat{K}\Gamma_i),
\end{aligned} 
\quad (26) $$

where $\Gamma_i = R^T r'_i$. Now, introducing (26) in (25) and the result in (19) and (20) yields

$$ 
\begin{aligned}
\bar{n}_{act} &= R(a + A\Gamma' + GK'), \\
\bar{c}_{act} &= R(b + B\Gamma' + HK'),
\end{aligned} 
\quad (27) $$

where A, B, G, H, a, and b depend on $\Gamma_i = R^T r'_i = \Gamma + K \times D_i + D'_i$ and are all proportional to $\tau_i$, as detailed in Appendix A. Finally, using (27) in (18) and the result in (17), and rearranging the obtained expressions, yields the ODEs governing the statics of a TACR in the expected explicit form

$$ \begin{pmatrix} r \\ R \\ \Gamma \\ K \end{pmatrix}' = \begin{pmatrix} R\Gamma \\ R\hat{K} \\ (\mathcal{H}_l + A) & G \\ B & \mathcal{H}_a + H \end{pmatrix}^{-1} \begin{pmatrix} c \\ d \end{pmatrix}, \quad (28) $$

where the expressions of c and d are detailed in Appendix A. This set of ODEs needs to be supplemented with the BCs of (17). Here the values of $n_+$ and $c_+$ are the forces and moments generated by the attachment of each tendon at the boundary, which, invoking the action-reaction principle, read

$$ 
\begin{aligned}
n_+ &= \sum_{i=1}^{m} -n_i(l) = -\sum_{i=1}^{m} \tau_i t_i, \\
c_+ &= \sum_{i=1}^{m} d_i \times (-n_i(l)) = -\sum_{i=1}^{m} d_i \times \tau_i t_i.
\end{aligned} 
\quad (29) $$

Finally, with these detailed BCs, (28) defines the same BVP as obtained in [22]. This BVP is the forward kineto-static BVP of a TACR, "forward" because its resolution provides outputs in motion (here strains), from inputs in force (here tendon tensions).

## V. LAGRANGIAN MODEL OF TACRS

In this section, the Lagrangian model of a TACR is derived from the principle of virtual work, starting from zero, i.e., calculating individually each contribution to the virtual work balance. In other words, to improve the intelligibility of the approach, we bring back TACR modeling to the canonical Lagrangian methodology. Note that this choice contrasts with recent publications on the topic, where the Lagrangian model is obtained either by projection of Cosserat partial differential equations with Jacobian matrices [33], [34], or by feeding a continuous Newton-Euler inverse algorithm with specific inputs [32].

Adopting this canonical viewpoint, Lagrangian modeling is achieved in two steps. In a first (kinematic) step, the Lagrangian approach consists of a reduction of the infinite dimensional configuration space (2) into a finite-dimensional configuration space. In a second step, this kinematic reduction is introduced in the principle of virtual work. As this principle holds in any definition of the configuration space, it allows to shift the static balance of Section III-B, from the configuration space (2), to the (finite dimensional) reduced one. This two-step process finally produces the reduced kineto-static model of a TACR.

### A. Strain Based Reduction

The purpose of the reduction is to shift the previously introduced continuous model to a finite dimensional model in terms of a vector of generalized strain coordinates q, similar to the joint coordinates of a rigid manipulator. g being in a non-commutative Lie group (i.e. a nonlinear manifold), the usual linear operations of interpolation or functional superimposition cannot directly be used to achieve the expected reduction. To circumvent this issue, the idea is to remark first that one can reconstruct any g by integrating (4) from $X = 0$ to $X = l$ starting from $g(0) = 1_{4 \times 4}$. Therefore the configuration space (2) of a clamped Cosserat rod can be alternatively defined by the functional space

$$ C_2 = \{\xi : X \in [0, l] \mapsto \xi(X) \in \mathbb{R}^6\}. \quad (30) $$

$C_2$ is a (functional) linear space, since $SE(3)$ in (2) is now replaced by the vector space $\mathbb{R}^6$. This alternate definition allows to apply usual procedures of linear reduction. According to the Ritz approach [42], the strain field $\epsilon = \xi - \xi_0$ is decomposed on a functional basis of "strain functions" as

$$ \epsilon = \Phi(X)q, \quad (31) $$

where q is a $\mathbb{R}^k$ vector of generalized strain coordinates and $\Phi$ is a $\mathbb{R}^{6 \times k}$ matrix of shape functions whose choice is fixed by the user. As a result, any configuration of the rod can be reconstructed in the inertial frame by integrating the ODE

$$ g' = g(\Phi q + \xi_0)^{\wedge} \quad (32) $$

from $X = 0$ where $g(0) = 1_{4 \times 4}$, to $X = l$. The purpose of the next section is to derive the static balance of a TACR in terms of the vector of generalized strain coordinates q of its backbone, i.e. on the finite dimensional configuration (vector)
space $\mathbb{R}^k$. This reduced model will take the usual matrix form of Lagrangian mechanics.

## B. Reduced Static Balance

Once the configuration space is defined, one needs to use one of the variational principles of Lagrangian mechanics to derive the model of the system in this space. In statics, the principle of virtual work applies. It considers the system in its static equilibrium where the external loads balance the internal restoring forces (stress), thus

$$ \delta W_{ext} + \delta W_{int} = 0, \quad (33) $$

where $\delta W_{ext}$ and $\delta W_{int}$ stand, respectively, for the work of external and internal forces along any virtual displacement (or variation) of the configuration, compatible with the geometric BCs and the parameterization used to describe the system. In the case of TACRs, the system considered is the backbone and the part of the tendons contained inside the robot. Since the tendons have negligible inertia and elasticity, the configuration of this system is that of the backbone alone. Thus any field of kinematically compatible virtual displacement along the robot can be defined by a small perturbation of g noted $\delta g$, with $\delta g(0) = 0_{4 \times 4}$. As for (4), this field of transformation perturbations is entirely described by a field of twist noted $\delta \zeta$ and defined by

$$ \delta g = g \delta \hat{\zeta}, \quad \text{with: } \delta \zeta(0) = 0_{6 \times 1}. \quad (34) $$

The condition at $X = 0$ is imposed by the compatibility of variations $\delta g$ with the geometric BC $g(0) = 1_{4 \times 4}$. When applied to the backbone, such virtual displacements generate variations of its strain field $\delta \epsilon = \delta(\xi - \xi_0) = \delta \xi$. These relative variations $\delta \xi$ are related to the absolute ones $\delta \zeta$ through the commutation relation

$$ \delta \zeta' = \delta \xi - a d_{\xi} \delta \zeta. \quad (35) $$

To derive this relation, it suffices to remark that, since variations $\delta$ do not affect the labels $X$, the variation and $X$-derivation are exchangeable, i.e.

$$ \delta(g') = (\delta g)' . \quad (36) $$

Then, using the definitions of $\delta \zeta$ and $\xi$ in terms of $g$, namely $\delta g = g \delta \hat{\zeta}$, and $g' = g \xi$ in (36), as well as the anti-hat operator, leads to (see Appendix B) the expected relation (35). Based on these preliminaries, one can now develop the balance of virtual works (33) as follows. The virtual work of internal elastic stress reads

$$ \delta W_{int} = - \int_{0}^{l} \delta \epsilon^{T} \mathcal{H} \epsilon \, dX. \quad (37) $$

As regards the work of external forces, it can be detailed as the sum

$$ \begin{aligned} \delta W_{ext} &= \delta \tilde{W}_{ext} + \delta W_{act} \\ &= \left( \int_{0}^{l} \delta \zeta^{T} \tilde{F} \, dX + \delta \zeta(l)^{T} \tilde{F}_{+} \right) + \left( \sum_{i=1}^{m} \delta l_{i} \tau_{i} \right), \quad (38) \end{aligned} $$

where $\delta \tilde{W}_{ext}$ stands for the virtual work of external gravity and contact forces, while $\delta W_{act}$ is the virtual work of tensions exerted on the tendons, with $\delta l_i$ the variation of length of the part of tendon $i$ contained in the robot. Now, assuming the tendons are inextensible, the kinematic relation between their length variations and the strain variations of the backbone is (see Appendix C)

$$ \begin{gathered} \delta l_i = \int_0^l \bar{J}_{li} \delta \epsilon \, dX, \\ \text{with: } \bar{J}_{li} = \frac{1}{\|\Gamma_i\|} (\Gamma_i^T, \Gamma_i^T) \begin{pmatrix} \hat{D}_i^T & 0_{3 \times 3} \\ 0_{3 \times 3} & 1_{3 \times 3} \end{pmatrix}, \quad (39) \end{gathered} $$

where the Jacobian (continuum) operator $\bar{J}_{li}$ is introduced, with $\Gamma_i = R^T r'_i = \Gamma + K \times D_i + D'_i$, and $D_i = R^T d_i$, the vector of the components of the radial position of the tendon $i$ in the cross-sectional frame of the backbone. Introducing this relation in the second term of (38), and exploiting the fact that $r'_i = 0_{3 \times 1}$, allows rewriting the virtual work of actuation in the form

$$ \delta W_{act} = \sum_{i=1}^{m} \int_{0}^{l} \delta \epsilon^{T} \bar{J}_{li}^{T} \tau_{i} \, dX = \int_{0}^{l} \delta \epsilon^{T} \left( \sum_{i=1}^{m} \bar{J}_{li}^{T} \tau_{i} \right) \, dX. \quad (40) $$

This last equation shows that the effect of tendons naturally appears as stress (and not forces) since they work along strain variations $\delta \epsilon$ (and not along configuration variations, i.e. virtual displacements $\delta \zeta$). Based on this remark, the field of actuation stress $\Lambda_{act}$ is defined as

$$ \begin{aligned} \delta W_{act} &= \int_{0}^{l} \delta \epsilon^{T} \Lambda_{act} \, dX \Rightarrow \\ \Lambda_{act} &= \left( \sum_{i=1}^{m} \tau_{i} \bar{J}_{li}^{T} \right) = \sum_{i=1}^{m} \frac{1}{\|\Gamma_{i}\|} \begin{pmatrix} D_{i} \times \Gamma_{i} \\ \Gamma_{i} \end{pmatrix} \tau_{i}, \quad (41) \end{aligned} $$

which can be integrated to the constitutive law (15) yielding the more general active-passive law

$$ \Lambda = \Lambda_{act} + \mathcal{H} \epsilon. \quad (42) $$

To detail the expression of the first term $\delta \tilde{W}_{ext}$ of (38), the virtual displacements $\delta \zeta$ must be expressed in terms of $\delta \epsilon$ and subsequently $\delta q$. To this end, the first step is to introduce $a d_{\xi} = Ad_{g}^{-1} Ad_{g}'$ in the commutation relation (35). The second step is to integrate the result between 0 and $X$. This provides the expression of any virtual displacement field $\delta \zeta$ that is compatible with the internal kinematics and the geometric BCs

$$ \begin{aligned} \delta \zeta(X) &= Ad_{g(X)}^{-1} \int_{0}^{X} Ad_{g(Y)} \delta \xi(Y) \, dY \\ &= \left[ Ad_{g(X)}^{-1} \int_{0}^{X} Ad_{g(Y)} \Phi \, dY \right] \delta q = J(X) \delta q, \quad (43) \end{aligned} $$

where $J(X)$ defines the Jacobian matrix, which maps any variation $\delta q$ to $\delta \zeta(X)$. To obtain (43), one needs the relation $\delta \zeta(0) = 0_{6 \times 6}$ (consequence of $g(0) = 1_{4 \times 4}$) as well as the variation of (31)

$$ \delta \xi = \delta \epsilon = \Phi \delta q. \quad (44) $$
Using (43) in $\delta \tilde{W}_{ext}$ of (38), and (44) in $\delta W_{int}$ of (37) as well as in $\delta W_{act}$ of (41), the balance of virtual works (33) reads

$$ \delta q^T Q_{ext} + \delta q^T K_{\epsilon\epsilon} q = \delta q^T Q_{act}. \quad (45) $$

(45), being true for any $\delta q$, provides the static balance of generalized forces

$$ Q_{ext} + K_{\epsilon\epsilon} q = Q_{act}. \quad (46) $$

This equation introduced the following three vectors of generalized forces. The $\mathbb{R}^k$ vector of generalized restoring forces

$$ Q_{rod} = \int_{0}^{l} \Phi^{T} \mathcal{H} \epsilon \, dX = \left( \int_{0}^{l} \Phi^{T} \mathcal{H} \Phi \, dX \right) q = K_{\epsilon\epsilon} q, \quad (47) $$

with $K_{\epsilon\epsilon}$ the $\mathbb{R}^{k \times k}$ matrix of generalized stiffness. The $\mathbb{R}^k$ vector of generalized forces of actuation

$$ \begin{aligned} Q_{act} &= -\int_{0}^{l} \Phi^{T} \Lambda_{act} \, dX \\ &= -\sum_{i=1}^{m} \left[ \int_{0}^{l} \frac{\Phi^{T}}{\|\Gamma_{i}\|} \left( \frac{D_{i} \times \Gamma_{i}}{\Gamma_{i}} \right) dX \right] \tau_{i} = L(q)\tau. \end{aligned} \quad (48) $$

And the $\mathbb{R}^k$ vector of generalized external forces

$$ Q_{ext} = - \int_{0}^{l} J^{T} \bar{F} \, dX - J(l)^{T} F_{+}, \quad (49) $$

which can be alternatively written in the equivalent form (see Appendix D)

$$ Q_{ext} = - \int_{0}^{l} \Phi^{T} \Lambda \, dX, \quad (50) $$

where $\Lambda$ is the solution of the inverse kineto-static BVP of a TACR, obtained by gathering (4) and (15),

$$ \begin{gathered} \begin{pmatrix} g \\ \Lambda \end{pmatrix}' = \begin{pmatrix} g \xi \\ ad_{\xi}^{T} \Lambda - \bar{F} \end{pmatrix}, \quad (51) \\ g(0) = 1_{4 \times 4}, \quad \Lambda(l) = F_{+}. \end{gathered} $$

This new BVP is called the inverse BVP, because it allows to calculate force-outputs (here the stress $\Lambda$) from the knowledge of motion-inputs (here the strain $\epsilon = \xi - \xi_0$). In other words, $\Lambda$ is the field of stress that balances the imposed external wrenches $\bar{F}$ and $F_+$ when the backbone is in the configuration $\xi$.

It should be noted here that, although expressions (49) and (50) are equivalent and dual to each other, choosing one or the other to calculate the generalized external forces, leads to fundamentally different algorithms. In details, if one uses (49), as in [33], [34], the external forces of (46) are calculated with Jacobian matrices according to a matrix projective process, often called Kane's method in the mechanical literature [43]. On the other hand, if one uses (50) as in [32], the same forces are calculated by solving the inverse BVP (51). In section VII.B, this second approach will be addressed yet with new numerical methods (compared to [32]), and applied to the simulation of a TACR in its quasi-static regime.

## VI. EQUIVALENCE OF THE TWO MODELS OF TENDON ACTUATION

Before Lagrangian reduction, the two approaches only differ by the model of tendon actuation (see Fig. 2). Hence, the purpose of this section is to demonstrate the equivalence of these two models. To this end, let us first reconsider (18) and recall that, in the Newtonian approach (see Section IV), the action of the tendons is modeled as external forces $\bar{n}_{act}$ and couples $\bar{c}_{act}$. Both are integrated in the wrench of external distributed loads $\bar{F}$ applied along the backbone as follows:

$$ \bar{F} = \begin{pmatrix} \bar{C} \\ \bar{N} \end{pmatrix} = \begin{pmatrix} \bar{C}_{ext} \\ \bar{N}_{ext} \end{pmatrix} + \begin{pmatrix} \bar{C}_{act} \\ \bar{N}_{act} \end{pmatrix} = \begin{pmatrix} R^T \bar{c}_{ext} \\ R^T \bar{n}_{ext} \end{pmatrix} + \begin{pmatrix} R^T \bar{c}_{act} \\ R^T \bar{n}_{act} \end{pmatrix}. \quad (52) $$

Whereas in the Lagrangian approach (see Section V), the same actions are modeled by the wrench of stress across the backbone $\Lambda_{act}$ whose expression is defined by (41), and which can be expressed in the inertial frame as

$$ \begin{pmatrix} c_{act} \\ n_{act} \end{pmatrix} = \begin{pmatrix} RC_{act} \\ RN_{act} \end{pmatrix} = \sum_{i=1}^{m} \begin{pmatrix} d_i \times t_i \\ t_i \end{pmatrix} \tau_i, \quad (53) $$

where remind that $R\Gamma_i = r'_i$, $t_i = r'_i/\|r'_i\|$ and $RD_i = d_i$. To demonstrate the equivalence of both approaches, the first step is to reconsider the BVP (17) on which the Newtonian approach is based. But, in contrast to Section IV-B, the action of tendons is no more modeled as a field of external forces and couples defined by (19), (20), and (23), but as the field of stress wrench across the backbone (53). This change of view point means that, in (17), $\bar{c}_{act}$ and $\bar{n}_{act}$ are removed from the model of external loads (18), to be replaced by internal forces and couples across the backbone $c_{act}$ and $n_{act}$. These internal forces superimpose to its usual elastic forces and couples, now noted $c_{rod}$ and $n_{rod}$, to form the full field of internal couples and forces along the robot considered as a stress-actuated backbone

$$ c = c_{rod} + c_{act}, \quad n = n_{rod} + n_{act}. \quad (54) $$

Note that once expressed in the cross-sectional frames, these two relations do define the active constitutive law (42), which takes the detailed form

$$ \Lambda = \Lambda_{act} + \mathcal{H} \epsilon = \begin{pmatrix} R^T c_{act} \\ R^T n_{act} \end{pmatrix} + \begin{pmatrix} \mathcal{H}_a (K - K_0) \\ \mathcal{H}_l (\Gamma - \Gamma_0) \end{pmatrix}, \quad (55) $$

where $n_{act}$ and $c_{act}$ are given by (53). As announced, substituting (53) and (54) in the stress balance (5), produces

$$ \begin{gathered} \left( n_{rod} + \sum_{i=1}^{m} \tau_i t_i \right)' + \bar{n}_{ext} = 0_{3 \times 1}, \\ \left( c_{rod} + \sum_{i=1}^{m} d_i \times (\tau_i t_i) \right)' + r' \times \left( n_{rod} + \sum_{i=1}^{m} \tau_i t_i \right) + \bar{c}_{ext} = 0_{3 \times 1}. \end{gathered} \quad (56) $$
Rearranging (56) with the usual composition rules of derivatives and using $d_i = r_i - r$ and $\tau'_i = 0_{3 \times 1}$, yields

$$ 
\begin{aligned}
n'_{rod} + \sum_{i=1}^{m} \tau_i t'_i + \bar{n}_{ext} &= 0_{3 \times 1}, \\
c'_{rod} + r' \times n_{rod} + \sum_{i=1}^{m} d_i \times (\tau_i t'_i) &+ \sum_{i=1}^{m} (r' + d'_i) \times (\tau_i t_i) + \bar{c}_{ext} = 0_{3 \times 1}.
\end{aligned} 
\quad (57) $$

Finally, using the relations $r' + d'_i = r'_i = \|r'_i\| t_i$ and $t_i \times t_i = 0_{3 \times 1}$ in (57), leads to the two stress balances of a TACR in the Newtonian approach (24) (with $n_{rod} = n, c_{rod} = c$), from which we deduced the BVP (28). Therefore, modeling the effect of tendons as external loads as in the Newtonian approach or as internal stress as in the Lagrangian approach leads to equivalent continuous models. In more details, the BVP (28) with the BCs of (17), with (29), of the Newtonian approach, is equivalent to the closed continuous formulation (46), with (41), (42), and (51), exploited by the Lagrangian approach.

## VII. NUMERICAL IMPLEMENTATION

Sections VII-A and VII-B deal with numerical implementation of the Newtonian and Lagrangian approaches, respectively. They both start with the case of a single-segment TACR and then extend the resolution to the multi-segment case.

### A. Numerical Implementation of the Newtonian Approach

In the Newtonian approach, the simulation of a TACR is achieved by solving the forward kineto-static BVP (28). Although this can be done by different numerical methods (finite differences, shooting, spectral methods ...), this paper adopts the shooting method, since it is largely dominant in the community.

1) **Newtonian Resolution of a Single-Segment TACR:** Solving the forward statics BVP (28) with the shooting method consists in finding the unknown proximal BCs $\Gamma(0)$ and $K(0)$ such that the known distal BCs $\Gamma(l)$ and $K(l)$, defined by (17), are fulfilled (see bottom left area of Fig. 2). This search is achieved iteratively by applying a root finding algorithm to the residual vector

$$ \mathcal{R}(\Gamma^*(0), K^*(0)) = \begin{pmatrix} \Gamma(l) - \Gamma^*(l) \\ K(l) - K^*(l) \end{pmatrix}, \quad (58) $$

where $\Gamma^*(l)$ and $K^*(l)$ are the distal values of $\Gamma$ and $K$ obtained by forward integrating the ODEs of (28) from the known proximal BCs $(r, R)(0)$ — supplemented with the guessed values $(\Gamma^*, K^*)(0)$ of the missing ones — to $X = l$.

This resolution proposes to iterate over the values of $\Gamma(0)$ and $K(0)$ but the relations (7) permit to equivalently iterate over $N(0)$ and $C(0)$ or any combination of these variables. Changing the state variables can be of interest, for instance, when the rod is very stiff for shear and extension, which may cause convergence issues. When computing the residual vector at each iteration, the kinematic ODE $R' = R\hat{K}$ is integrated

using quaternions, which ensures that $R$ remains in SO(3). Note also that other methods based on Magnus expansions in SO(3) $\times \mathbb{R}^3$ or SE(3) can be used [33], [44].

Choosing the initial guess is a crucial step that may determine the outcome of the implemented algorithm. When the initial guess is too far from the solution, algorithms may converge to a local minimum or fail to converge due to gradients that lead to infeasible points. Nonetheless, finding a good initial guess may be a tricky task. A solution to address this issue is to start from a vector of zeros. If the TACR is subject to small loads, the solution will be close to zero and the solver will likely converge. If, on the contrary, the TACR is subject to higher loads, the strategy consists in dividing the total load of the problem in a number of increasing loading steps such that the first step can be solved starting with the zero vector. Then, the subsequent steps are solved for, using the previous solution as input for the following guess. Another advantage of specifying the loading history is that the final solution can be controlled in the case of multiple possible solutions for a given load [45].

2) **Newtonian Resolution of a Multi-Segment TACR:** In a multi-segment robot, the poses as well as the internal forces and couples propagate from one segment to the next. Consequently, their values at the distal end of the first segment are used as inputs for the ODEs of the second segment and so on, until the most distal segment, where the BCs are evaluated. The implementation of multi-segment robots is thus straightforward based on the single-segment case and does not increase the size of the problem's state vector.

One particular case that, to the authors' knowledge, has not previously been reported in the literature, deserves a special focus: when a tendon undergoes a slope discontinuity. The tension applied to such a tendon will create a force, acting at the location of the discontinuity, which needs to be taken into account in the implementation. This happens for instance in a multi-segment TACR where a tendon is routed following a convergent routing path in the distal segment but runs parallel to the backbone in the proximal segments (see Fig. 3-Scenario C.). Practically, for a tendon $i$ undergoing a slope discontinuity between the segments $j$ and $j + 1$, the shift in routing orientation is equal to

$$ \llbracket t_i \rrbracket_j = \lim_{\epsilon \to 0} (t_i(X_j^+) - t_i(X_j^-)), \quad (59) $$

where $X_j^\pm = X_j \pm \epsilon$, $X_j$ is the arc length abscissa of the backbone cross-section to which is attached the disk connecting segments $j$ and $j + 1$, and $\llbracket f \rrbracket_j$ denotes the jump of any field $f$ when crossing the junction between segments $j$ and $j + 1$. Let us remind that $t_i = r'_i / \|r'_i\|$ denotes the unit tangent vector to the $i^{th}$ tendon. Applying Newton's laws to the tendon at the junction point, provides the jump of force and couple transmitted to the backbone by a tendon $i$ subject to a slope discontinuity at junction $j$

$$ \llbracket n \rrbracket_j = \tau_i \llbracket t_i \rrbracket_j, \quad \llbracket c \rrbracket_j = d_i \times \llbracket n \rrbracket_j. \quad (60) $$

Note that (60) is no more than a discrete version of (23). (60) will need to be added to the static balance (29) while forward integrating the BVP (28), each time a tendon slope
discontinuity is encountered. The slope discontinuity $[\|t_i\|]_j$ can be obtained with the relation $[\Gamma_i / \|\Gamma_i\|]_j = R^T [\|t_i\|]_j$. Hence, (60) is rewritten in its cross-sectional frame version

$$[N]_j = \tau_i \left[ \frac{\Gamma_i}{\|\Gamma_i\|} \right]_j, \quad [C]_j = D_i \times [N]_j, \quad (61)$$

where $\Gamma_i^{\pm}$ is computed similarly to (26), with

$$\Gamma_i(X_j^{\pm}) = \Gamma(X_j^{\pm}) + \hat{K}(X_j^{\pm}) D_i(X_j) + D_i'(X_j^{\pm}). \quad (62)$$

In practice, this calculation is achieved by using the following jump relations on the space-rates:

$$\begin{pmatrix} \Gamma \\ K \end{pmatrix}_{+} = \begin{pmatrix} \Gamma \\ K \end{pmatrix}_{-} + \begin{pmatrix} H_l + A & G \\ B & H_a + H \end{pmatrix}_{-}^{-1} \left[ \begin{bmatrix} c \\ d \end{bmatrix} \right]_j, \quad (63)$$

where the expressions of $A$, $G$, $B$, $H$, $c$, and $d$ are again given by Appendix A, except that the continuous $X$-derivatives $D'_i$ and $D''_i$ are replaced by their jumps $[D_i]_j = 0_{3 \times 1}$ and $[D'_i]_j \neq 0_{3 \times 1}$, respectively while, owing to the continuity of poses, all other quantities are evaluated at $X = X_j^-$. The contribution of this term can be observed in simulations in the accompanying video where the Newtonian approach is shown to fail without the slope discontinuity compensation term.

Finally, using (63) from segment to segment while piecewise integrating the ODEs of (28) from $X = 0$ to $l$, allows calculating the residual vector of (58) in the general case with possible tendon slope discontinuities between segments.

## B. Numerical Implementation of the Lagrangian Approach

Following the standard numerical techniques of the Lagrangian approach, the quasi-static simulation of a TACR is performed here by solving the implicit system of non-linear algebraic equations defined in (46). Note that this is different from [32], where the resolution is performed by the explicit time integration of an overdamped equivalent system. Also, compared to [32]–[34] regarding the Lagrangian approach, this paper further extends the simulations to Reissner rods (i.e. considering shear and extension).

### 1) Lagrangian Resolution of a Single-Segment TACR:

In the Lagrangian approach, the numerical resolution of the static problem consists in calculating the vector of generalized strain coordinates $q = (q_1 .. q_k)^T$ that fulfills the reduced static balance (46), for any imposed vector of tendon tensions $\tau = (\tau_1 .. \tau_m)^T$ (see bottom right area of Fig. 2). Since this balance defines a set of algebraic nonlinear equations, one can apply any root finding algorithm to solve for the vector of residuals $\mathcal{R}$ that takes the generic form

$$\mathcal{R}(q, \tau) = 0_{k \times 1}, \quad \text{with: } \mathcal{R}(q, \tau) = Q_{\text{ext}} + K_{\epsilon\epsilon}q - Q_{\text{act}}. \quad (64)$$

It is of common use in nonlinear structural statics (and dynamics) to use Newton-Raphson's method, which allows to update the $q$-vector at each step of a loop, iteratively reducing the residual, according to the linear relation

$$q^+ = q - J^{-1}(q)\mathcal{R}(q, \tau), \quad (65)$$

where $J(q) = \frac{\partial \mathcal{R}}{\partial q}$ denotes the $\mathbb{R}^{k \times k}$ Jacobian matrix of the residual vector, and $q^+$ is the updated value of $q$. Therefore,

applying this method requires to calculate the vector of residuals $\mathcal{R}$ and its Jacobian matrix $J$.

To numerically compute the residual vector $\mathcal{R}$ from any $q$-vector, the expressions (47), (48), and (50), are added to the inverse BVP (51) through the state variable

$$y(X) = - \int_{X}^{l} \Phi^{T} (\Lambda_{\text{act}} - \Lambda + \mathcal{H}\Phi q) dY \quad (66)$$

such that $y(0) = \mathcal{R}$. The state variables $(g, \Lambda, y)$ are then governed by the augmented inverse BVP, which is fed with $\xi = \xi_0 + \Phi q$, and reads

$$\begin{pmatrix} g \\ \Lambda \\ y \end{pmatrix}' = \begin{pmatrix} g(\xi_0 + \Phi q)^{\wedge} \\ ad_{(\xi_0+\Phi q)}^{T} \Lambda - \bar{F} \\ \Phi^{T}(\Lambda_{\text{act}} - \Lambda + \mathcal{H}\Phi q) \end{pmatrix}, \quad (67) \\ g(0) = 1_{4 \times 4}, \quad \Lambda(l) = F_{+}, \quad y(l) = 0_{k \times 1},$$

where $\Lambda_{\text{act}}$ is a function of $q$ and $\tau$, given by (41), while $\bar{F}$ depends on the context. For instance, if the robot is only subject to gravity, $\bar{F} = (0_{1 \times 3}, R^T a_g^T A \rho)^T$, with $a_g$ the acceleration gravity field expressed in the inertial frame, $A$ the area of the robot's section, and $\rho$ its density.

Note that, similarly to the inverse kineto-static BVP (51), (67) possesses interesting properties that makes its resolution straightforward compared to that of the forward BVP (28). In particular, it can be solved with two decoupled passes (i.e. without resorting to the shooting algorithm), with standard explicit space integrators. One can first integrate forward (from $X = 0$ to $l$), the ODE in $g$, and then integrate backward (from $X = l$ to 0) the two other ODEs in $\Lambda$ and $y$ producing $\mathcal{R}(q) = y(0)$. This method has been proposed and interpreted as a Newton-Euler computed torque algorithm in [32], similar to those developed for rigid multi-body systems [46]. In this article, we exploit a further property of the inverse BVP, and remark that since the $q$-vector is an imposed input, (67) is linear with respect to the state variables $(g, \Lambda, y)$. As a result, one can apply a spectral collocation method and replace the previous explicit ODE integrations, by the resolution of some linear algebraic systems with respect to the vector of the state variables on a Chebyshev grid [37].

In the continuum robotics community, Jacobian matrices are often calculated numerically. This paper proposes a more accurate calculation of $J$, based on the exact linearization of the residual vector. To this end, the BVP (67) is linearized with respect to its input $q$ by propagating the differential consequences of a variation $\Delta q$ as

$$\begin{pmatrix} \Delta \zeta \\ \Delta \Lambda \\ \Delta y \end{pmatrix}' = \begin{pmatrix} -ad_{(\xi_0+\Phi q)} \Delta \zeta + \Phi \Delta q \\ ad_{\Phi \Delta q}^{T} \Lambda + ad_{(\xi_0+\Phi q)}^{T} \Delta \Lambda - \Delta \bar{F} \\ \Phi^{T}(\Delta \Lambda_{\text{act}} - \Delta \Lambda + \mathcal{H}\Phi \Delta q) \end{pmatrix}, \quad (68) \\ \Delta \zeta(0) = 0_{6 \times 1}, \quad \Delta \Lambda(l) = \Delta F_{+}, \quad \Delta y(l) = 0_{k \times 1},$$

where $\Delta \zeta$ is the differential of $g$ in $\text{se}(3)$ defined by $\Delta \zeta = (g^{-1}g)^{\vee}$, while $\Delta \Lambda_{\text{act}}$, $\Delta \bar{F}$, and $\Delta F_{+}$ are the differentials of the model of actuation stress and external forces, respectively. Note here that, the variation $\Delta$ being only generated by the variation of the configuration $\Delta q$, $\Delta \tau = 0_{m \times 1}$ in $\Delta \Lambda_{\text{act}}$. Moreover, the dependence of $\Lambda_{\text{act}}$ on $q$ is very weak and
may be neglected in the calculation of the Jacobian. Thus, $\Delta\Lambda_{act} = 0_{6 \times 1}$ when solving (68).

The BVP (68) defines the augmented inverse tangent BVP (TBVP) of the robot. The full BVP obtained by gathering (67) and (68) is linear, and can, once again, be solved with a spectral method applied forward and backward, as for the BVP. Now, by virtue of the linear identity $\Delta y(0) = \Delta \mathcal{R} = J(q)\Delta q$, applying the unit input vector $\Delta q = \delta_{\alpha}$ to (68), where $\delta_{\alpha}$ is a $\mathbb{R}^{k \times k}$ vector of zeros, except the entry $\alpha$, which is equal to 1, (68) yields the column $\alpha$ of $J$. Hence, repeating the process for $\alpha = 1..k$ allows the Jacobian to be filled column by column. As a last remark, due to its nonlinearity, the spectral method cannot be directly applied to the forward statics BVP (17).

## 2) Lagrangian Resolution of a Multi-Segment TACR:

The Lagrangian model of a multi-segment TACR is simply obtained by applying the principle of virtual work to the system defined by the serial connection of the $\chi$ segments. All segments are subject to the external forces applied by the environment (e.g. gravity), those exerted by the tendons across the basis, as well as the restoring internal forces. Such a virtual work balance takes the generic form

$$ \sum_{j=1}^{\chi} \delta \tilde{W}_{ext,j} + \delta W_{act,j} + \delta W_{int,j} = 0, \quad (69) $$

with $j$ the index of segments and where each of the three contributions can be defined and detailed as in the single-segment case of (38). In particular, applying the same strain reduction segment by segment provides the vector of generalized strain coordinates of the entire robot $q = (q_{(1)}..q_{(\chi)})^T$. From now on, "•(j)" denotes a generalized vector or matrix, related to the segment $j$. Introducing these reduced virtual works in (69) allows to write, for any $\delta q$,

$$ \delta q^T (Q_{ext} + K_{\epsilon\epsilon} q - Q_{act}) = 0_{k \times 1}, \quad (70) $$

where $K_{\epsilon\epsilon} = \text{diag}(K_{\epsilon\epsilon,(1)}..K_{\epsilon\epsilon,(\chi)})$, $Q_{ext} = (Q_{ext,(1)}..Q_{ext,(\chi)})^T$, and $Q_{act} = (Q_{act,(1)}..Q_{act,(\chi)})^T$. In other words, the multi-segment TACR model is simply obtained by appending the matrices and vectors of each segment one after another.

To numerically solve this multi-segment model, the above single-segment method can be extended by considering a BVP of the form (67) for each of the segments, and its associated TBVP (68). Using normalized arc-length variables along each segment, these $\chi$ BVPs and TBVPs, are now connected through their BCs as follows:

$$ g_j(1) = g_{j+1}(0), \quad \Lambda_j(1) = \Lambda_{j+1}(0), \quad (71) $$

$$ \Delta\zeta_j(1) = \Delta\zeta_{j+1}(0), \quad \Delta\Lambda_j(1) = \Delta\Lambda_{j+1}(0). \quad (72) $$

Note that, as it is usually the case in Lagrangian mechanics, the above relations introduce no jumps on the stress. Indeed, the inter-segment forces of (60) are considered as internal forces of the multi-segment system that do not work in any virtual displacement field compatible with the inter-segment connections. Finally, the multi-segment solution is obtained by solving (65) applied to the whole structure. The Jacobian and the residual are computed, as in the single-segment case, by forward and backward integrations of (67) and (68), and starting from $g_1(0) = 1_{4 \times 4}$ and $\Lambda_{\chi}(1) = F_+$, respectively.

<table><caption>Table II
GEOMETRY AND MATERIAL PARAMETERS OF THE SIMULATED ROBOT.</caption><thead><tr><th>Parameter Name</th><th>Parameter Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Backbone diameter</td><td>R<sub>b</sub></td><td>0.4 mm</td></tr><tr><td>Backbone length (A, B, C)</td><td>l</td><td>242 mm</td></tr><tr><td>Backbone length (D, E)</td><td>l = l<sub>j</sub> × 3</td><td>300 mm</td></tr><tr><td>Tendon offset</td><td>R<sub>t</sub></td><td>8 mm</td></tr><tr><td>Gravity constant</td><td>||a<sub>g</sub>||</td><td>9.81 N/kg</td></tr><tr><td>Robot equivalent density<sup>2</sup></td><td>ρ = <sup>0.47 N/m</sup>⁄<sub>||a<sub>g</sub>||π(R<sub>b</sub>)<sup>2</sup></sub></td><td>95 · 10<sup>3</sup>kg/m<sup>3</sup></td></tr><tr><td>Young's modulus</td><td>E</td><td>210 GPa</td></tr><tr><td>Poisson's ratio</td><td>ν</td><td>0.3125</td></tr></tbody></table>

# VIII. PERFORMANCE COMPARISON

The objective of this section is to compare the behavior of both approaches when modeling various TACRs. This paper proceeds through the comparison of five example scenarios of increasing complexity and that are representative of the TACRs studied in the literature and numerous other configuration possibilities. The robots and their routings are represented in Fig. 3 and 3D views are shown in the accompanying video. To confront our simulations with results from the literature, we chose to model the TACR presented in [22], but with various tendon routings. Its geometry and material parameters are given in Table II. In all scenarios, the robot is oriented with its base pointing upwards, opposite to the action of gravity.

Both approaches also involve some numerical parameters. For the Newtonian approach, all segments' ODEs are integrated using the Runge-Kutta Dormand-Prince method with Matlab's `ode45` function. The 'ResTol' and 'AbsTol' parameters are both specified to $10^{-8}$. The global BVP (28) is solved using Matlab's `fsolve` function with the Levenberg-Marquardt algorithm with all parameters set to defaults. For the Lagrangian approach, the ODEs are integrated using a spectral method over a Chebyshev grid of 10 nodes. The strain distributions are projected on Legendre polynomials of orders 5 for torsion; 7 for bending in both directions; and 3 for extension and shearing in both directions. Summing all polynomial orders results in a total of 28 shape functions. Finally, the residual of the BVP (67) is brought to $< 10^{-8}$ with Newton-Raphson's method<fnref n="3" />.

When running simulations, particularly in the case of high loads, the algorithms may fail to converge if the actuation load is applied as a single step. Progressively increasing the tension in the tendons and iteratively building upon the previous solutions may solve this problem. This method will be used for some of the simulated cases and will be referred to as the number of needed loading steps.

As none of the codes for either approach were fully optimized, the computation times of the simulations are not provided. Information regarding which approach is more efficient for which case is discussed in Section IX.

<footnote n="2">²The equivalent density is the density applied to the backbone capturing the weight of the backbone itself, disks, and tendons. The value of the robot's self-weight is from [22].</footnote>

<footnote n="3">³The algorithm parameters were chosen empirically during preliminary experiments such that the simulation results had converged for each approach.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3. Left, front view (positive y-axis pointing upwards) of the different routings corresponding to the investigated scenarios. Right, top view of the robots’ routing disks, with the tendon offset $R_t$. The top views of the first three scenarios are superimposed. The numbers next to the tendons correspond to their index in each scenario. For details on each specific routing, please refer to the corresponding section.</figcaption></figure>

### A. Single-Segment Parallel Routing

Let us start with a simple case: a single-segment TACR with 3 tendons routed parallel at 120° around ($R_t = 8$ mm) the backbone (see Fig. 3-Scenario A.). 216 simulations are run with each approach corresponding to the possible combinations of the integer values of tension between 0 and 5 N for each tendon. The simulations are analyzed by pair (Newtonian and Lagrangian approaches) and the position and orientation of the tip are compared in Table III. This table also displays the main robot characteristics for each simulated scenario. Fig. 4a shows a sample of 20 pairs of the 432 simulations. One can see how both approaches render the same robot shape.

For this scenario, depending on the load combination in the tendons, all cases in both methods could be solved with only 1 loading step, except 3 cases (1.39%) with the Newtonian approach that required 2 loading steps.

### B. Single-Segment Convergent Routing

Let us now change the parallel routing from the TACR above to convergent routing. Examples of existing physical robots with this kind of routing can be found in [24], [47]. In this case, 2 tendons are routed at opposite positions of the robot, linearly convergent towards its end, starting at offset $R_t = 8$ mm at the proximal end of the robot and ending coincident with the backbone at its tip (see Fig. 3-Scenario B.). The cross-sectional frame position of the tendons now depends on the reference length parameter X as follows:

$$D_i(X) = (0, \pm R_t(1 - X/l), 0)^T .$$

Simulations are run for all combinations of tendon tensions over the integer values between 0 and 8 N. The maximum tension is higher than in Section VIII-A in order to operate the robot over a comparable workspace. With the two-tendon robot described, these combinations yield 162 simulations. Fig. 4b displays a sample of 20 pairs of these simulations. As in the previous case, one can see how both approaches render the same robot shapes. The values for the differences

between tip positions and orientations are similar to those for the previous scenario (see Table III). As for the previous scenario, a majority of cases of this scenario could be solved with 1 loading step. 6.17% of the cases with the Newtonian approach required 2 loading steps.

### C. Single-Segment Helical Routing

In this scenario, the routing is set to 2 helically routed tendons in a full turn around the length of the backbone, lying opposite to each other (see Fig. 3-Scenario C.). As such, the cross-sectional frame position of the tendon i is equal to

$$D_i(X) = R_t \left( 0, \cos\left(\frac{2\pi}{l}X + i\pi\right), \sin\left(\frac{2\pi}{l}X + i\pi\right) \right)^T,$$

with $R_t = 8$ mm.

The simulation results for tendon tensions varying over the integer values between 1 and 10 N for each tendon independently for both approaches are shown in Fig. 4c. The tension range is increased compared to the previous scenarios in order to cover a comparable workspace with the studied routing. Similar existing physical robots of the literature can be found in [22], [23]. The differences in position and orientation are still several orders of magnitudes lower than the robot’s length (see Table III). For this batch of simulations, all cases of both approaches were solved with a single loading step.

### D. Multi-Segment Parallel Routing

As the implementation of multi-segment TACRs considerably varies from one approach to the other (see Sections VII-A2 and VII-B2), it is interesting to compare simulation results with multi-segment TACRs as well. This section starts with a multi-segment TACR scenario involving no routing path discontinuities, thus not being subject to the matter. The next section follows up with a scenario that does involve routing path discontinuities (59).

The multi-segment TACR studied is a three-segment robot with 3 parallel routed tendons per segment at 120° around
<figure><img src="image_5.png" /><figcaption>Figure 4. A sample of 20 random pairs of the simulation results for single-segment TACRs with (a) parallel routing, (b) convergent routing, and (c) helical routing. The geometry and material parameters of the simulated robot are given in Table II. Full and dashed lines represent simulations obtained using the Newtonian approach and the Lagrangian approach, respectively.</figcaption></figure>

<table><caption>Table III
SYNTHESIS OF RESULTS FOR THE FIVE SIMULATED SCENARIOS.</caption><thead><tr><th>Scenario**</th><th></th><th>A</th><th>B</th><th>C</th><th>D</th><th>E</th></tr></thead><tbody><tr><td rowspan="3">Robot characteristics</td><td>Type of routing</td><td>Parallel</td><td>Convergent</td><td>Helical</td><td>Parallel</td><td>Convergent</td></tr><tr><td>Number of segments</td><td>1</td><td>1</td><td>1</td><td>3</td><td>3</td></tr><tr><td>Robot length [mm]</td><td>242</td><td>242</td><td>242</td><td>300</td><td>300</td></tr><tr><td rowspan="3">Tip position difference*** [mm]</td><td>Median value</td><td>5.43 · 10<sup>-4</sup></td><td>9.40 · 10<sup>-3</sup></td><td>1.42 · 10<sup>-1</sup></td><td>3.42 · 10<sup>-1</sup></td><td>6.93 · 10<sup>-2</sup></td></tr><tr><td>Lower bound 95% confidence interval</td><td>2.41 · 10<sup>-9</sup></td><td>8.77 · 10<sup>-15</sup></td><td>4.48 · 10<sup>-3</sup></td><td>1.00 · 10<sup>-5</sup></td><td>1.11 · 10<sup>-2</sup></td></tr><tr><td>Upper bound 95% confidence interval</td><td>3.13 · 10<sup>-3</sup></td><td>2.32 · 10<sup>-2</sup></td><td>7.21 · 10<sup>-1</sup></td><td>1.29</td><td>2.76 · 10<sup>-1</sup></td></tr><tr><td rowspan="3">Tip orientation difference*** [°]</td><td>Median value</td><td>3.26 · 10<sup>-4</sup></td><td>1.10 · 10<sup>-3</sup></td><td>1.99 · 10<sup>-1</sup></td><td>5.56 · 10<sup>-2</sup></td><td>1.32 · 10<sup>-2</sup></td></tr><tr><td>Lower bound 95% confidence interval</td><td>1.03 · 10<sup>-7</sup></td><td>5.42 · 10<sup>-8</sup></td><td>8.09 · 10<sup>-3</sup></td><td>3.39 · 10<sup>-7</sup></td><td>2.78 · 10<sup>-3</sup></td></tr><tr><td>Upper bound 95% confidence interval</td><td>1.65 · 10<sup>-3</sup></td><td>3.20 · 10<sup>-3</sup></td><td>5.74 · 10<sup>-1</sup></td><td>2.14 · 10<sup>-1</sup></td><td>4.74 · 10<sup>-2</sup></td></tr></tbody></table>

the backbone. Tendons (1,2,3) terminate in the first segment, tendons (4,5,6) in the second segment, and tendons (7,8,9) in the third segment (see Fig. 3-Scenario D.). The robot geometry and material parameters are the same as for the previous robots, except for the length of the segments, shortening them to $l_j = 100$ mm for a total robot length of $l = 3l_j = 300$ mm. Indeed, a too slender robot is subject to instabilities and for some tendon load combinations, multiple solutions would be possible. This simulation scenario can be compared to the physical multi-segment TACR examples of the literature [7], [10], [25].

505 tendon tension combinations that cover the workspace are simulated with integer values of tendon tensions varying from 0 to 5 N. Fig. 5a displays a sample of 20 pairs of these 1010 simulations. The robot shapes with both approaches superimpose and the differences in position and orientation of the tip remain very small (0.43% of the robot's length) (see Table III). In this scenario, 1.58% of the cases required 2 loading steps with the Newtonian approach. All other cases were solved with 1 loading step.

### E. Multi-Segment Convergent Routing

Here, the same robot as in the previous scenario is simulated, but with 2 convergent tendons per segment. The tendons are only routed convergent in their terminating segment and are routed parallel to the backbone elsewhere, following:

$$D_{i,j}(X) = (-1)^i R_t \varphi_j(X) (0, \cos(\theta_j), \sin(\theta_j))^T,$$

with

$$\varphi_j(X) = \begin{cases} 1 & \text{for } X < (j-1)l_j \\ 1 - \frac{X-(j-1)l_j}{l_j} & \text{for } (j-1)l_j \le X \le j l_j \\ 0 & \text{for } X > j l_j, \end{cases}$$

where $D_{i,j}(X)$, stands for the routing path of tendon $i = 1..6$ that ends in segment $j = 1..3$ and $\theta_j = (j-1)2\pi/3$. As such, the pair of tendons of each segment run on opposite sides of the backbone. The tendons of the first segment are aligned with the $xy$-plane. For the following segments, the tendons are shifted 120° counterclockwise with respect to the previous segment (see Fig. 3-Scenario E.).

Fig. 5b displays a sample of the 686 simulations, corresponding to 343 actuation combinations covering the workspace with integer values of tendon tensions ranging from 0 to 4 N. The differences in position and orientation of the tip are reported in the last column of Table III. For this batch of simulations, all cases of both approaches were solved with a single loading step. Building upon a design featuring convergent routings, as presented in [24], this scenario, with multiple such segments, involves tendon slope discontinuities between the segments. Still, the values are similar to those of other scenarios. The accompanying video shows how the Newtonian approach fails without the slope discontinuity compensation term. In some cases, the discrepancies can go up to 50% of the robots' length.

<footnote>** The letters of the scenarios correspond to the subsections in Section VIII.</footnote>

<footnote>*** Differences between the results that are obtained with both approaches.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 5. A Sample of 20 random pairs of the simulation results for the three-segment TACRs with (a) parallel routing and (b) convergent routing. The geometry and material parameters of the simulated robot are given in Table II. Full and dashed lines represent simulations obtained using the Newtonian approach and the Lagrangian approach, respectively.</figcaption></figure>

IX. DISCUSSION

In terms of mathematical modeling, the Newtonian approach makes use of Newton’s action-reaction principle and gathers ad hoc equations to form the global BVP of the problem. The main challenge consists in reshaping the BVP, to find an explicit form suited for solving with a shooting method. This challenge is addressed through complex algebraic manipulations [22]. On the other hand, the philosophy of Lagrangian mechanics is to consider the system as a whole, removing the need for isolating subsystems and making them interact. The difficulty here resides in reducing the infinite dimensional problem to a finite dimensional problem in both the spaces of kinematics and statics, the latter involving evolved calculus [32]. The advantage of this approach is that the reduced form of the problem resembles that of a classical rigid series manipulator. As such, the generalized strain coordinates are analogous to the joint coordinates, which is particularly useful for control purposes. Another advantage of the Lagrangian approach is that it can be applied systematically to a wide variety of systems (closed loops, lumped joints, etc.) without having to redefine the vector of residuals for each specific case. Redefining residuals may be delicate as soon as one deviates from the simple case of the multi-segment robots studied here.

Numerically, the results from Section VIII show that both approaches produce almost identical results. The differences can partially be explained by the way the external distributed loads are taken into account in each approach. In the Newtonian approach, these loads are included in (28), which is integrated from $X = 0$ to $X = l$, whereas in the Lagrangian approach, the distributed loads $\bar{F}$ are defined over $X \in ]0, l[$. Still, all tip position differences are comprised between 0 and 1.29 mm which is 8.00% of the robots’ diameters and stays within 0.43% of the robots’ lengths. These discrepancies are negligible compared to the modeling errors that are reported in studies that involve experimental setups, which go over 7% of the robots’ lengths [23].

The Lagrangian approach proposes a reduction of the problem that can be adapted for each particular system. Shape functions that determine the size of the residual must be specified to represent the strains in the various segments of

the modeled robot. A great variety of choices are available, and it may be challenging to find the strain basis that is best suited to a particular problem while keeping the computational cost low. Including the routing schemes of the tendons as shape functions was proven to be an efficient choice [33], but complex interactions between the robot and its environment will however require including higher dimension generic shape functions. Moreover, when modeling multi-segment robots, a number of shape functions must be added for every new segment. Thus, the combination of high dimension shape functions and multiple segments may lead to excessively long residual vectors. The benefit of this reduction is that the variational principles are transferred to the numerical implementation, which are proven to be numerically robust. As such, for virtually any given number of degrees of freedom, the problem remains well conditioned and can be solved with Newton-Raphson’s method, even for very high actuation and external loads. In our simulations, the benefit of preserving the variational principles through the reduction is demonstrated by the fact that all simulations could be solved with a single loading step.

The Newtonian approach is the easiest to implement, particularly for those not familiar with Lagrangian mechanics. Its infinite configuration space removes the need for determining shape functions. The consequence is that the robustness of the mechanical principles is lost at the implementation. Therefore, the approach is particularly dependent on the initial guess that is fed to the solver. This approach is thus well suited for problems where an approximation of the solution is known (which is seldom the case) or for problems whose solution is close to zero. In other words, the method is efficient for solving non-complex cases with small deformations, but increasing the problem complexity or modeling high strains requires proceeding in steps, which is computationally costly. In the Newtonian approach, multi-segment robots require more complex action-reaction considerations (see Section VII-A). Moreover, increasing the number of segments entangles the problem without providing more degrees of freedom to the state vector or the residual, and therefore makes the convergence less certain for this kind of robots. On the other hand,
the constant residual size goes together with an approximately
linear relation between the number of segments and the
computational cost. Besides, in this approach, computing the
Jacobian matrix of the residual vector is complicated and using
finite differences to estimate the gradient is less efficient and
may lead to inaccuracies.

Let us note that the above discussed reduction, influencing
the size of the residual of the problem in the Lagrangian
approach, must not be confused with the spatial discretization.
Both approaches require to numerically integrate the ODEs
along the backbones, requiring to discretize the backbones
over a number of nodes or integration points. Increasing the
number of these integration points contributes to the accuracy
of the simulations but also increases the computational cost,
meaning an equilibrium must be found. On this point, note that
the spectral integration benefits from an exponential conver-
gence, compared to the polynomial convergence of the usual
Runge-Kutta integrators. This enhanced convergence permits
to integrate over a reduced number of integration points with
the Lagrangian approach (typically 10 to 20).

We now draw conclusions from the whole processes used
to model and simulate the presented scenarios and more gen-
erally from our experience with these processes. Apart from
the numerical implementations presented in this paper, other
numerical methods can be used. For instance, the BVP of the
Newtonian approach can be solved with finite differences or
spectral methods, while [32]–[34] present different numerical
methods for the Lagrangian approach. In what follows, we
compare the derivation of the Newtonian approach solved
with a shooting method and the derivation of the Lagrangian
approach solved with Newton-Raphson’s method. Overall,
the Newtonian approach is more efficient for two types of
problems: (i) problems that include one or only a few segments
and involving low tendon tensions, and (ii) problems solved
iteratively with a known close solution as can be the case in
robot control. Indeed, starting from a close solution favors the
convergence, especially for complex structures. These more
complex structures with many segments and tendons, and
problems involving large deformations for which the solution
is completely unknown have a better chance on convergence
with the Lagrangian approach. As for the Newtonian approach,
the Lagrangian approach is more efficient when solving prob-
lems iteratively. And even more so, considering the fact that
establishing the generalized stiffness matrix $K_{ee}$, the most
costly operation, is a constant function of the shape functions
$\Phi$ and, therefore, is computed only once.

X. CONCLUSION

In this paper, we show how the Newtonian approach and
the Lagrangian approach derive from the same kineto-static
model of a Cosserat rod. We have unified the derivation to
highlight the similarities between both approaches, which was
previously a strenuous task due to their different frameworks
and community backgrounds [22], [32]. Substantial additional
material is formulated in the derivation of the Lagrangian
approach presented in this paper to improve its intelligibility.
Also, we address the implementation of tendon slope discon-

tinuities in the Newtonian approach, which was missing from
the literature [22], [23], [48].

It is, for the first time, proved mathematically that both
approaches are equivalent. The simulation results show that
either of the approaches can be used interchangeably to model
TACRs and obtain identical solutions (all tip positions differ
within of 0.43% of the simulated robots' lengths). Still, some
considerations summarized in Section IX, are to be taken into
account when choosing one or the other approach.

Finally, the developments and findings presented in this
paper apply to TACRs but can be extrapolated to model other
types of continuum robots. Our simulations did not take into
account any other external forces than gravity, but future work
including interactions with the environment can incorporate
them easily as described in the developments of Sections IV,
V, and VII. Also, building on the present paper that compared
the two approaches for robots in their quasi-static regime,
future work will focus on the dynamic case. We hope the
availability of the code will serve as a building block for the
continuum robotics community.

APPENDIX A

DETAILED EXPRESSIONS AND NOTATIONS USED IN (27)
AND (28)

The ODEs (27) make use of the notations

$$
A_i = -\tau_i \frac{\hat{\Gamma}_i^2}{\|\hat{\Gamma}_i\|^3}, \qquad A = \sum_{i=1}^{m} A_i,
$$

$$
B_i = \hat{D}_i A_i, \qquad B = \sum_{i=1}^{m} B_i,
$$

$$
G = - \sum_{i=1}^{m} A_i \hat{D}_i = B^T, \qquad H = - \sum_{i=1}^{m} B_i \hat{D}_i,
$$

$$
a_i = A_i (\hat{K}\Gamma_i + \hat{K}D_i' + D_i''), \quad a = \sum_{i=1}^{m} a_i,
$$

$$
b_i = \hat{D}_i a_i, \qquad b = \sum_{i=1}^{m} b_i.
$$

The BVP (28) makes use of the above notations, as well as

$$
c = \mathcal{H}_l \Gamma'_0 - \hat{K} \mathcal{H}_l (\Gamma - \Gamma_0) - R^T \bar{n}_{ext} - a,
$$

$$
d = \mathcal{H}_a K'_0 - \hat{K} \mathcal{H}_a (K - K_0) \\
\qquad - \hat{\Gamma} \mathcal{H}_l (\Gamma - \Gamma_0) - R^T \bar{c}_{ext} - b.
$$

APPENDIX B

DERIVATION OF THE COMMUTATION RELATION

Let us recall the definitions $\delta g = g\delta\hat{\zeta}$ and $g' = g\hat{\zeta}$ of (4)
and (34), respectively. To derive (35), these definitions are
introduced in (36) as follows:

$$
\delta (g \hat{\xi}) = (g \delta \hat{\zeta})' \Leftrightarrow \delta g \hat{\xi} + g \delta \hat{\zeta} = g' \delta \hat{\zeta} + g (\delta \hat{\zeta})'. \quad (73)
$$

Using the same definitions, once more, leads to

$$
g(\delta\hat{\zeta}\hat{\xi} + \delta\hat{\xi}) = g(\hat{\xi}\delta\hat{\zeta} + (\delta\hat{\zeta})'), \quad (74)
$$
which, holding for any $g$, yields

$$ \delta\hat{\xi} = (\delta\hat{\zeta})' + \hat{\xi}\delta\hat{\zeta} - \delta\hat{\zeta}\hat{\xi} = (\delta\hat{\zeta})' + [\hat{\xi}, \delta\hat{\zeta}], \quad (75) $$

where [•, •] denotes the usual commutator of matrices. Finally, using the anti-hat (●$^V$) operation, one can show that $[\hat{\xi}, \delta\hat{\zeta}]^V = ad_\xi\delta\zeta$, and get the expected commutation relation in terms of twists (35).

## APPENDIX C

### PROOF OF (39)

The expression of the length of a tendon

$$ l_i = \int_0^l \|r'_i\| \, dX = \int_0^l (r_i'^T r_i')^{1/2} \, dX \quad (76) $$

leads to

$$ \delta l_i = \int_0^l \delta \|r'_i\| \, dX = \int_0^l \frac{\delta r_i'^T r_i'}{\|r_i'\|} \, dX. \quad (77) $$

In order to rewrite the last term of (77), recall that $r'_i = R\Gamma_i$, $\delta r'_i = \delta(R\Gamma_i)$, and thus

$$ \delta r_i'^T r_i' = (\delta R \Gamma_i + R \delta \Gamma_i)^T (R \Gamma_i) = (R^T \delta R \Gamma_i)^T \Gamma_i + \delta \Gamma_i^T \Gamma_i. $$

Now remark that, since $R \in SO(3)$, there always exists a vector $\delta\psi \in \mathbb{R}^3$ such that $R^T \delta R = \delta\hat{\psi}$. As a result

$$ \frac{\delta r_i'^T r_i'}{\|r_i'\|} = \frac{(\delta\psi \times \Gamma_i)^T \Gamma_i + \delta\Gamma_i^T \Gamma_i}{\|\Gamma_i\|} = \frac{\delta\Gamma_i^T \Gamma_i}{\|\Gamma_i\|}. \quad (78) $$

Finally, introducing the expression $\Gamma_i = \Gamma + K \times D_i + D_i'$ into (78) and the result of this substitution into (77) provides, after factorization of $\delta\epsilon = (\delta K^T, \delta\Gamma^T)^T$, the expected expression (39).

## APPENDIX D

### PROOF OF (50)

The relation $Ad'_g = Ad_g ad_\xi$, leads to

$$ ad_\xi^T = -Ad_g^T Ad_g^{-T'}. \quad (79) $$

Then, introducing this relation in the balance (13), one can integrate it, and obtain the balance of wrenches

$$ \Lambda(X) = Ad_{g(X)}^T \int_X Ad_{g(Y)}^{-T} \bar{F}(Y) \, dY + Ad_{g(X)}^T Ad_{g(l)}^{-T} F_+ . \quad (80) $$

Projecting this balance on the strain basis (pre-multiplying by $\Phi^T$ and integrating the products over $[0, l]$), and using by-part integrations, yields the reduced balance of stress

$$ -\int_{0}^{l} \Phi^{T} \Lambda dX = -\int_{0}^{l} J^{T} \bar{F} dX - J(l)^{T} F_{+}, \quad (81) $$

which is no more than the expression of $Q_{ext}$ given by (49).

## REFERENCES

[1] J. Burgner-Kahrs, D. C. Rucker, and H. Choset, “Continuum robots for medical applications: A survey,” *IEEE Trans. Robot.*, vol. 31, no. 6, pp. 1261–1280, dec 2015.

[2] P. E. Dupont, N. Simaan, H. Choset, and D. C. Rucker, “Continuum robots for medical interventions,” *Proc. IEEE*, vol. 110, no. 7, pp. 847–870, jul 2022.

[3] S. M. H. Sadati, S. E. Naghibi, A. Shiva, I. D. Walker, K. Althoefer, and T. Nanayakkara, “Mechanics of continuum manipulators, a comparative study of five methods with experiments,” in *Proc. Annu. Conf. Towards Auton. Robot. Syst.*, jul 2017, pp. 686–702.

[4] M. T. Chikhaoui and B. Rosa, “Chapter 8 - modeling and control strategies for flexible devices,” in *Endorobotics*. New York, NY, USA: Academic, 2022, pp. 187–213.

[5] T. Mahl, A. E. Mayer, A. Hildebrandt, and O. Sawodny, “A variable curvature modeling approach for kinematic control of continuum manipulators,” in *Proc. IEEE Amer. Contr. Conf.*, jun 2013, pp. 4945–4950.

[6] C. D. Santina and D. Rus, “Control oriented modeling of soft robots: The polynomial curvature case,” *IEEE Robot. Automat. Lett.*, vol. 5, no. 2, pp. 290–298, apr 2020.

[7] B. A. Jones and I. D. Walker, “Kinematics for multisection continuum robots,” *IEEE Trans. Robot.*, vol. 22, no. 1, pp. 43–55, feb 2006.

[8] ——, “Practical kinematics for real-time implementation of continuum robots,” *IEEE Trans. Robot.*, vol. 22, no. 6, pp. 1087–1099, dec 2006.

[9] S. Neppalli, M. A. Csencsits, B. A. Jones, and I. D. Walker, “Closed-form inverse kinematics for continuum manipulators,” *Adv. Robot.*, vol. 23, no. 15, pp. 2077–2091, jan 2009.

[10] P. Gonthina, M. Wooten, I. S. Godage, and I. Walker, “Mechanics for tendon actuated multisection continuum arms,” in *Proc. IEEE Int. Conf. Robot. Automat.*, apr 2020, pp. 3896–3902.

[11] R. J. Webster and B. A. Jones, “Design and kinematic modeling of constant curvature continuum robots: A review,” *Int. J. Robot. Res.*, vol. 29, no. 13, pp. 1661–1683, jun 2010.

[12] P. Rao, Q. Peyron, S. Lilge, and J. Burgner-Kahrs, “How to model tendon-driven continuum robots and benchmark modelling performance,” *Front. Robot. AI*, vol. 7, no. 2, pp. 223–243, feb 2021.

[13] D. B. Camarillo, C. F. Milne, C. R. Carlson, M. R. Zinn, and J. K. Salisbury, “Mechanics modeling of tendon-driven continuum manipulators,” *IEEE Trans. Robot.*, vol. 24, no. 6, pp. 1262–1273, dec 2008.

[14] R. J. Webster, J. M. Romano, and N. J. Cowan, “Mechanics of precurved-tube continuum robots,” *IEEE Trans. Robot.*, vol. 25, no. 1, pp. 67–78, feb 2009.

[15] Y. Liu and F. Alambeigi, “Effect of external and internal loads on tension loss of tendon-driven continuum manipulators,” *IEEE Robotics and Automation Letters*, vol. 6, no. 2, pp. 1606–1613, apr 2021.

[16] D. Trivedi, A. Lotfi, and C. D. Rahn, “Geometrically exact models for soft robotic manipulators,” *IEEE Trans. Robot.*, vol. 24, no. 4, pp. 773–780, aug 2008.

[17] D. C. Rucker and R. J. Webster, “Mechanics-based modeling of bending and torsion in active cannulas,” in *Proc. IEEE Int. Conf. Biomed. Robot. Biomechatronics*, Oct 2008, pp. 704–709.

[18] P. E. Dupont, J. Lock, B. Itkowitz, and E. Butler, “Design and control of concentric-tube robots,” *IEEE Trans. Robot.*, vol. 26, no. 2, pp. 209–225, apr 2010.

[19] D. C. Rucker, B. A. Jones, and R. J. Webster, “A geometrically exact model for externally loaded concentric-tube continuum robots,” *IEEE Trans. Robot.*, vol. 26, no. 5, pp. 769–780, oct 2010.

[20] J. Ha, F. C. Park, and P. E. Dupont, “Elastic stability of concentric tube robots subject to external loads,” *IEEE Trans. Biomed. Eng.*, vol. 63, no. 6, pp. 1116–1128, jun 2016.

[21] J. Ha, G. Fagogenis, and P. E. Dupont, “Modeling tube clearance and bounding the effect of friction in concentric tube robot kinematics,” *IEEE Trans. Robot.*, vol. 35, no. 2, pp. 353–370, apr 2019.

[22] D. C. Rucker and R. J. Webster, “ statics and dynamics of continuum robots with general tendon routing and external loading,” *IEEE Trans. Robot.*, vol. 27, no. 6, pp. 1033–1044, dec 2011.

[23] J. Starke, E. Amanov, M. T. Chikhaoui, and J. Burgner-Kahrs, “On the merits of helical tendon routing in continuum robots,” in *Proc. IEEE/RSJ Int. Conf. Intell. Robots Syst.*, sep 2017, pp. 6470–6476.

[24] K. Oliver-Butler, J. Till, and D. C. Rucker, “Continuum robot stiffness under external loads and prescribed tendon displacements,” *IEEE Trans. Robot.*, vol. 35, no. 2, pp. 403–419, apr 2019.

[25] E. Amanov, T.-D. Nguyen, and J. Burgner-Kahrs, “Tendon-driven continuum robots with extensible sections—a model-based evaluation of path-following motions,” *Int. J. Robot. Res.*, vol. 40, no. 1, pp. 7–23, jan 2021.
[26] A. L. Orekhov, V. A. Aloi, and D. C. Rucker, "Modeling parallel continuum robots with general intermediate constraints," in *Proc. IEEE Int. Conf. Robot. Automat.*, may 2017, pp. 6142-6149.

[27] J. Wang, J. Peine, and P. E. Dupont, "Eccentric tube robots as multiarmed steerable sheaths," *IEEE Trans. Robot.*, 2021.

[28] C. B. Black, J. Till, and D. C. Rucker, "Parallel continuum robots: Modeling, analysis, and actuation-based force sensing," *IEEE Trans. Robot.*, vol. 34, no. 1, pp. 29-47, feb 2018.

[29] M. Mahvash and P. E. Dupont, "Stiffness control of surgical continuum manipulators," *IEEE Trans. Robot.*, vol. 27, no. 2, pp. 334-345, apr 2011.

[30] J. Till, V. Aloi, and D. C. Rucker, "Real-time dynamics of soft and continuum robots based on Cosserat rod models," *Int. J. Robot. Res.*, vol. 38, no. 6, pp. 723-746, may 2019.

[31] F. Janabi-Sharifi, A. Jalali, and I. D. Walker, "Cosserat rod-based dynamic modeling of tendon-driven continuum robots: A tutorial," *IEEE Access*, vol. 9, pp. 68703-68719, may 2021.

[32] 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.

[33] 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. Automat. Lett.*, vol. 5, no. 3, pp. 4006-4013, jul 2020.

[34] F. Renda, C. Messer, D. C. Rucker, and F. Boyer, "A sliding-rod variable-strain model for concentric tube robots," *IEEE Robot. Automat. Lett.*, vol. 6, no. 2, pp. 3451-3458, apr 2021.

[35] H. Goldstein, *Classical Mechanics*, 3rd ed. Reading, MA, USA: Addison Wesley, 2001.

[36] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu, *The Finite Element Method: Its Basis and Fundamentals*, 6th ed. Oxford, UK: Butterworth-Heinemann, 2005.

[37] L. N. Trefethen, *Spectral methods in MATLAB*. Philadelphia, PA, USA: SIAM, 2000.

[38] S. S. Antman, *Nonlinear Problems of Elasticity*, 2nd ed. New York, NY, USA: Springer-Verlag, 2005.

[39] J. E. Marsden and T. S. Ratiu, *Introduction to mechanics and symmetry*, 2nd ed. New York, NY, USA: Springer-Verlag, 1999.

[40] F. Boyer, M. Porez, and A. Leroyer, "Poincare Cosserat equations for the lighthill three-dimensional large amplitude elongated body theory: Application to robotics," *J. of Nonlinear Sci.*, vol. 20, no. 1, pp. 47-79, feb 2010.

[41] K. M. Lynch and F. C. Park, *Modern Robotics*. Cambridge, U.K.: Cambridge University Press, 2017.

[42] J. K. L. MacDonald, "Successive approximations by the Rayleigh-Ritz variation method," *Phys. Rev.*, vol. 43, no. 10, pp. 830-833, may 1933.

[43] T. R. Kane and D. A. Levinson, *Dynamics: Theory and Applications*. New York, NY, USA: McGraw-Hill (Inc.), 1985.

[44] A. L. Orekhov and N. Simaan, "Solving Cosserat rod models via collocation and the Magnus expansion," in *Proc. IEEE/RSJ Int. Conf. Intell. Robots Syst.*, oct 2020, pp. 8653-8660.

[45] T. Bretl and Z. McCarthy, "Quasi-static manipulation of a Kirchhoff elastic rod based on a geometric analysis of equilibrium configurations," *Int. J. Robot. Res.*, vol. 33, no. 1, pp. 48-68, jan 2014.

[46] R. Featherstone, *Rigid Body Dynamics Algorithms*. New York, NY, USA: Springer-Verlag, 2007.

[47] C. Laschi, M. Cianchetti, B. Mazzolai, L. Margheri, M. Follador, and P. Dario, "Soft robot arm inspired by the octopus," *Adv. Robot.*, vol. 26, no. 7, pp. 709-727, apr 2012.

[48] M. Neumann and J. Burgner-Kahrs, "Considerations for follow-the-leader motion of extensible tendon-driven continuum robots," in *Proc. IEEE Int. Conf. Robot. Automat.*, may 2016, pp. 917-923.

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

**Vincent Lebastard** joined the École Normale Supérieure de Cachan, France, in 2001, and received an aggregation from the ministry of the Education in Electrical Engineering in 2002. He received a Ph.D degree from the University of Nantes, France, in 2007. He is currently a Assistant Professor with the Department of Automatic Control, IMT-Atlantique, Nantes, France, where he is a member of the Robotics Team, Laboratoire des Sciences du Numérique de Nantes (LS2N). Dr. lebastard received the French "La Recherche Prize" in 2014, for his

works on artificial electric sense. His research interests include biorobotics and dynamic modeling.

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

**Frédéric Boyer** was born in France in 1967. He received the Diploma 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 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, IMT-Atlantique, Nantes, France, where he is a member of the Robotics Team, Laboratoire des Sciences du Numérique de Nantes (LS2N). His current research interests include structural dynamics, geometric mechanics, and biorobotics (locomotion dynamics and underwater electric sensing). Dr. Boyer received 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. He has coordinated several national projects and one European FP7-FET project on a reconfigurable eel-like robot able to navigate with electric sense.

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

**Jocelyne Troccaz** Research Director at CNRS, working in TIMC Laboratory, Grenoble, France. Graduated in Computer Science. PhD in robotics in 1986, Institut National Polytechnique de Grenoble. CNRS Research fellow from 1988. Specialized in Medical Robotics and Computer-Assisted Medical Interventions. Her main interests are in the development of new robotic paradigms and devices and in image registration. Active in several clinical areas (urology, radiotherapy, cardiac surgery, orthopedics, etc.) in collaboration with Grenoble University Hospital and La Pitié Salpétrière Paris Hospital. Thanks to industrial transfer, hundreds of thousands of patients, worldwide, benefited from technology and systems she developed. IEEE Fellow, MICCAI Fellow. Dr. Troccaz has been an associate editor of the IEEE Transactions on Robotics and Automation and of the IEEE Transactions on Robotics. She is currently member of the steering committee of IEEE Transactions on Medical Robotics and Bionics and editorial board member of Medical Image Analysis. Member of the French Academy of Sciences. Member of the French Academy of Surgery. Recipient of several awards and medals: Award from the French Academy of Surgery in 2014; Silver Medal from CNRS in 2015; Chevalier de la Légion d'Honneur in 2016, Enduring Impact Award MICCAI society in 2022.

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

**Matthias Tummers** was born in Belgium, in 1993. He received the University degree in biomedical engineering from the UCLouvain, Louvain-la-Neuve, Belgium, in 2016. From 2017 to 2019, he was a Research Assistant and Research Engineer with the Institute of Mechanics, Materials and Civil Engineering, Louvain-la-Neuve, Belgium and the TIMC Laboratory, Grenoble, France, respectively. He is currently working toward the Ph.D. degree within the Department of Computer Aided Medical Interventions team at the TIMC Laboratory. His research

interests include medical instrumentation, the design and control of miniature medical robots and mechanical modeling of continuum robots.

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

**Benoît Rosa** received the Ph.D. degree in robotics from Université Pierre et Marie Curie, Paris, France, in 2013, for which he received the best PhD award from the CNRS Research group in Robotics. From 2013 to 2015 he was a postdoctoral researcher at KU Leuven in Belgium, and from 2015 to 2017 he was a postdoctoral fellow at Boston Children's Hospital in the USA. Since 2017 he is a permanent CNRS Researcher at the ICube lab in Strasbourg. His research interests cover continuum robotics for minimally invasive surgery and image-guided surgery.
<figure><img src="image_12.png" /></figure>

**M. Taha Chikhaoui** received his Ph. D. in Automatic Control from the University of Franche-Comté in 2016. His research was performed at FEMTO-ST Institute in Besançon, France, where he was affiliated with the Automatic Control and Micro-Mechatronic Systems (AS2M) Department. He is recipient of the Prix de Thèse 2016. He received an Electromechanical Engineering Diploma in 2013 from the High Private School of Engineering and Technologies (ESPRIT) in Tunis, Tunisia, where he has ranked first of his specialty. He was postdoctoral researcher at the Laboratory for Continuum Robotics, Leibniz Universität Hannover in Hanover, Germany between 2017 and 2018. Since 2019, he is a Research Scientist at CNRS (The French National Centre for Scientific Research), appointed to TIMC Laboratory in Grenoble, France with the Computer-Assisted Medical Interventions (GMCAO) team. The main focus of his research is continuum robotics for (bio-)medical applications. His interests are based on an interdisciplinary approach including robotic design, modeling, innovative mechatronics, and control strategies based on perception, among others.