Matthias Tummers
email: matthias.tummers@univ-grenoble-alpes.fr.

Member, IEEE Vincent Lebastard

Frédéric Boyer

Fellow, IEEE Jocelyne Troccaz

Member, IEEE Benoît Rosa

Member, IEEE M Taha Chikhaoui

Cosserat Rod Modeling of Continuum Robots from Newtonian and Lagrangian Perspectives

Keywords: Modeling, Control, and Learning for Soft Robots, Flexible Robots, Tendon/Wire Mechanism, Cosserat Beams

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.

I. INTRODUCTION

C ONTINUUM robots are slender actuated flexible struc- tures 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 [START_REF] Burgner-Kahrs | Continuum robots for medical applications: A survey[END_REF], [START_REF] Dupont | Continuum robots for medical interventions[END_REF].

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 [START_REF] Burgner-Kahrs | Continuum robots for medical applications: A survey[END_REF]- [START_REF] Chikhaoui | Chapter 8 -modeling and control strategies for flexible devices[END_REF]. Two main modeling categories can be identified: geometry based models and mechanical models. Geometry based models [START_REF] Mahl | A variable curvature modeling approach for kinematic control of continuum manipulators[END_REF], [START_REF] Santina | Control oriented modeling of soft robots: The polynomial curvature case[END_REF], including the well known constant curvature approach [START_REF] Jones | Kinematics for multisection continuum robots[END_REF]- [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF], 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 [START_REF] Rao | How to model tendon-driven continuum robots and benchmark modelling performance[END_REF]. 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 [START_REF] Camarillo | Mechanics modeling of tendon-driven continuum manipulators[END_REF]- [START_REF] Liu | Effect of external and internal loads on tension loss of tendon-driven continuum manipulators[END_REF] 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 [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF]- [START_REF] Rucker | A geometrically exact model for externally loaded concentric-tube continuum robots[END_REF]. The approach was applied to various continuum manipulators as concentric tube continuum robots [START_REF] Ha | Elastic stability of concentric tube robots subject to external loads[END_REF], [START_REF] Ha | Modeling tube clearance and bounding the effect of friction in concentric tube robot kinematics[END_REF], tendon actuated continuum robots (TACRs) [START_REF] Rucker | Statics and dynamics of continuum robots with general tendon routing and external loading[END_REF]- [START_REF] Amanov | Tendon-driven continuum robots with extensible sections-a model-based evaluation of path-following motions[END_REF] and multi-backbone continuum robots [START_REF] Orekhov | Modeling parallel continuum robots with general intermediate constraints[END_REF]- [START_REF] Black | Parallel continuum robots: Modeling, analysis, and actuation-based force sensing[END_REF]. The method has been studied extensively for modeling the statics and dynamics [START_REF] Mahvash | Stiffness control of surgical continuum manipulators[END_REF]- [START_REF] Janabi-Sharifi | Cosserat rod-based dynamic modeling of tendon-driven continuum robots: A tutorial[END_REF] 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 [START_REF] Orekhov | Modeling parallel continuum robots with general intermediate constraints[END_REF], [START_REF] Wang | Eccentric tube robots as multiarmed steerable sheaths[END_REF]. 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 [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF]- [START_REF] Renda | A sliding-rod variablestrain model for concentric tube robots[END_REF]. 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 [START_REF] Goldstein | Classical Mechanics[END_REF]. 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 [START_REF] Renda | A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation[END_REF]. It allows the application of efficient numerical approximation 0000-0000/00$00.00 © 2023 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission See https://www.ieee.org/publications/rights/index.html for more information.

methods based on variational formulations, both in statics and in dynamics [START_REF] Zienkiewicz | The Finite Element Method: Its Basis and Fundamentals[END_REF]. 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 [START_REF] Rucker | Statics and dynamics of continuum robots with general tendon routing and external loading[END_REF], 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 [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF], 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 [START_REF] Trefethen | Spectral methods in MATLAB[END_REF] 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 [START_REF] Rucker | Statics and dynamics of continuum robots with general tendon routing and external loading[END_REF] 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 [START_REF] Antman | Nonlinear Problems of Elasticity[END_REF] and the existing literature on the modeling of continuum robots [START_REF] Rucker | A geometrically exact model for externally loaded concentric-tube continuum robots[END_REF]. On the other hand, the Lagrangian approach privileged notations from geometric mechanics [START_REF] Marsden | Introduction to mechanics and symmetry[END_REF] and its applications to rigid multi-body systems and continuous media, including fluid mechanics [START_REF] Boyer | Poincare Cosserat equations for the lighthill three-dimensional large amplitude elongated body theory: Application to robotics[END_REF]. This choice also eases the application of developments on manipulators to other contexts as e.g. eellike 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 R 3 whereas the Lagrangian approach makes extensive use of wrenches and twists, expressed in R 6 . In the literature on the Newtonian approach, these 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 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 R 3 can be expressed either in the inertial frame or in the crosssectional frames, which can sometimes lead to confusion1 . To avoid ambiguity, this paper adopts the convention of Juan Carlos Simo of noting a vector of 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 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) [START_REF] Marsden | Introduction to mechanics and symmetry[END_REF].

Newtonian approach

Lagrangian approach (+ this paper) Dimensions Units Reference length parameter

s ∈ [0, l] X ∈ [0, l] R m Position p(s) r(X) R 3 m Orientation R(s), (z-axis aligned) R(X), (x-axis aligned) R 3×3 Cross-sectional frame g(s) = (p, R)(s) g(X) = (R, r)(X) R 4×4 Linear rate of change (shear & extension) v(s) Γ(X) R 3 Angular rate of change (bending & torsion) u(s) K(X) R 3 m -1 Space-rate twist - ξ = (K T , Γ T ) T R 6 (m -1 , -) Initial value or prior to deformation * • * • 0 Total length of the robot l l R m Strain - ϵ = ξ -ξ 0 R 6 (m -1 , -) Related to the i th tendon • i • i Number of tendons n m R Length of the i th tendon l i l i R m Related to the j th segment • j • j Number of segments m χ R Length of the j th segment l j l j R m Cross-sectional frame position of the i th tendon r i (s) = [x i (s) y i (s) 0] T D i (X) = (0, D i,Y , D i,Z )(X) R 3 m Inertial frame tendon position p i (s) = p + Rr i r i (X) = r + RD i R 3 m Tangent to the path of the i th tendon ṗi /∥ ṗi ∥ t i = r ′ i /∥r ′ i ∥ R Tension in the i th tendon τ i τ i R N Vector of tensions - τ = (τ 1 .. τm) T R m N Stress (bending & torsion) m(s) C(X) R 3 Nm Stress (shear & extension) n(s) N (X) R 3 N Stress (wrench) - Λ(X) = C(X) T , N (X) T T R 6
(Nm, N) Angular and linear stiffness matrices

K bt (s), Kse(s) Ha(X), H l (X) R 3×3 N, Nm 2 Hookean stiffness matrix or H(X) = diag(Ha(X), H l (X)) R 6×6 (N, Nm 2 ) Generalized stiffness - Kϵϵ R k×k ** Total external distributed couples l(s) c(X), C(X) R 3 Nm/m Total external distributed forces f (s) n(X), N (X) R 3 N/m Total external distributed wrench - F (X) = C(X) T , N (X) T T R 6 (Nm/m, N/m) External tip couples - c + , C + R 3 Nm External tip forces - n + , N + R 3 N External tip wrench - F + = C T + , N T + T R 6 (Nm, N) Relative to external loads (not tendons) - •ext External distributed forces (not tendons) f e (s) next(X) R 3 N/m External distributed couples (not tendons) le(s) cext(X) R 3
Nm/m Relative to rod elastics -• rod Relative to actuation by the tendons -•act Distributed forces applied by the tendons

f t (s) nact(X) R 3 N/m Distributed couples applied by the tendons lt(s) cact(X) R 3 Nm/m Distributed forces applied to the i th tendon f i (s) ni (X) R 3 N/m Stress of actuation - Λact(X) R 6 (Nm, N) Generalized external forces - Qext R k ** Generalized restoring forces - Q rod R k ** Generalized forces of actuation - Qact R k ** Matrix of actuation - L R k×m ** Number of shape functions - k R Shape functions - Φ(X) = (Φ 1 .. Φ k ) T R 6×k ** Generalized strain coordinates - q = (q 1 .. q k ) T R k ** First and second derivative wrt. s or X •, • • ′ , • ′′ Residual vector - R
Another difference concerns the partial derivation with respect to the rod arc length variable in its reference (stressless) 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. n or 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 [START_REF] Goldstein | Classical Mechanics[END_REF]. 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 R 3 , then W denotes the skew-symmetric matrix of R 3×3 such that W =   0 -W 3 W 2 W 3 0 -W 1 -W 2 W 1 0   .
The reciprocal is denoted

• ∨ such that ( W ) ∨ = W . Similarly, if W = (W 1 , W 2 , W 3 ) T and V = (V 1 , V 2 , V 3 ) T are two vectors of R 3 , with ν = (W T , V T ) T a vector of R 6 , then ν denotes the matrix of R 4×4 as ν = W V 0 1×3 0 ,
where W is the skew-symmetric matrix associated to W ∈ 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 ν and homogeneous transformation g,

ν = W V and g = R r 0 1×3 1 , by the two matrices of R 6×6 ad ν = W 0 3×3 V W and Ad g = R 0 3×3 rR R .
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 ν will be used to model the inertial derivative of the cross-sectional basis as it moves along its centerline with a twist ν 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.

R(X)

x z y X X1 -n(X 1 ) -c(X 1 ) X2 n(X 2 ) c(X 2 ) c(X), n(X) inertial frame y z x r(X) g(X)
proximal side distal side 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 n(X) and couples c(X), and its internal restoring forces n(X) and couples c(X).

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 ∈ [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 crosssectional 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 ∈ [0, l] → r(X) ∈ R 3 , and another in the Lie group of rotations R : X ∈ [0, l] → R(X) ∈ SO(3). Together, r and R define a field of homogeneous transformation that describes the entire rod

g : X ∈ [0, l] → g(X) = R(X) r(X) 0 1×3 1 ∈ SE(3). (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) = 1 4×4 , i.e.:

C 1 = {g : X ∈ [0, l] → g(X) ∈ SE(3) , with: g(0) = 1 4×4 } .
(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 Γ and K of R 3 such that

r ′ = RΓ , R ′ = R K , (3)
where the • ′ symbol denotes differentiation with respect to X. When combined, the fields of linear and angular space-rate Γ and K, define a unique field of twist ξ such that

g ′ = g ξ (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 n and couple 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 ′ + n = 0 3×1 , c ′ +r ′ ×n+c = 0 3×1 , ∀X ∈]0, l[ . (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×1 , R(0) = 1 3×3 , n(l) = n + , c(l) = c + . (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 are counted positively when exerted by the the distal piece of rod (i.e. after arc-length X) onto the proximal piece (i.e. before arclength 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 • 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 Γ -Γ 0 and K -K 0 . Moreover, assuming linear material elasticity, the internal forces n and couples c are related to Γ and K through the constitutive laws

n = RH l (Γ -Γ 0 ) , c = RH a (K -K 0 ) , (7)
where, ∀X ∈ [0, l], H l (X) and 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Γ, n = RN , c = RC, n = R N , and c = R C, as follows:

(RN ) ′ + R N = 0 3×1 , (RC) ′ + (RΓ) × (RN ) + R C = 0 3×1 . (8)
Then, via composition rules of derivatives along with the relation

R ′ = R K RN ′ + R KN + R N = 0 3×1 , RC ′ + R KC + (RΓ) × (RN ) + R C = 0 3×1 , (9)
which, being true for any R and owing to the definition of the hat operator, provides

N ′ + K × N + N = 0 3×1 , C ′ + K × C + Γ × N + C = 0 3×1 . ( 10
)
Once supplemented with the BCs

r(0) = 0 3×1 , N (l) = R(l) T n + , R(0) = 1 3×3 , C(l) = R(l) T c + , (11)
which are deduced from their inertial frame version (6), equations [START_REF] Gonthina | Mechanics for tendon actuated multisection continuum arms[END_REF] stand for the ODEs of a Cosserat rod in the cross-sectional frame. Following the convention of [START_REF] Lynch | Modern Robotics[END_REF], all pairs of angular and linear vectors gather in R 6 vectors

C N ′ + K × C + Γ × N K × N + C N = 0 3×1 0 3×1 . (12)
Let us introduce the notations Λ = C T , N T T for the R 6 field of stress wrench, F = CT , N T T for the 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 [START_REF] Rao | How to model tendon-driven continuum robots and benchmark modelling performance[END_REF] in the more compact (geometric) form as a balance of wrenches expressed in the cross-sectional frame

Λ ′ -ad T ξ Λ + F = 0 6×1 . (13)
To which one needs to add the BCs

g(0) = 1 4×4 , Λ(l) = F + . (14)
Similarly to the twist form of ( 4), ( 7) can be rewritten as

Λ = Hϵ = H (ξ -ξ 0 ) , (15)
where

ϵ = ξ -ξ 0 = (K -K 0 ) T , (Γ -Γ 0 ) T T
is the R 6 field of strain and H = diag (H a , H l ) is the R 6×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) , h -(x(0)) = 0 , h + (x(l)) = 0 , ( 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 dim(h

-) + dim(h + ) = dim(x).
Based on this definition, it is straightforward to show that using the constitutive laws [START_REF] Jones | Kinematics for multisection continuum robots[END_REF] 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

    r R Γ K     ′ =          RΓ R K Γ ′ 0 -H -1 l KH l + H ′ l (Γ -Γ 0 ) + R T n K ′ 0 -H -1 a KH a + H ′ a (K -K 0 ) + ΓH l (Γ -Γ 0 ) + R T c          , r(0) = 0 3×1 , Γ(l) = Γ 0 + H -1 l R(l) T n + , R(0) = 1 3×3 , K(l) = K 0 + H -1 a R(l) T c + . ( 17
)
BC of (

(51)

(31) (32) (38)
(41)(42) (46) (46) (31)
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.

(17) defines a BVP in the explicit state form [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF], with x = (r, R, Γ, K) being the state variables. Note that using constitutive relations in alternate ways provides another equivalent BVP where strain variables (Γ, 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 [START_REF] Jones | Kinematics for multisection continuum robots[END_REF]), the considered robots have constant material properties along their length (i.e. H ′ = 0 6×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 [START_REF] Rucker | Mechanics-based modeling of bending and torsion in active cannulas[END_REF]. 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 [START_REF] Rucker | Mechanics-based modeling of bending and torsion in active cannulas[END_REF], the distributed loads are spelled out as

n = next + nact , c = cext + cact , (18)
where next and cext model the external loads applied by the environment such as gravity or contacts, while nact and cact 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 ni . Therefore, summing these contributions provides the distributed force of the m tendons on the backbone

nact = m i=1 -n i . (19)
As a tendon cannot transmit any couples, the distributed moment of the tendons exerted on the backbone cact is equal to only the sum of the cross products of each moment arm with each force

cact = m i=1 (r i -r) × (-n i ) = - m i=1 d i × ni , (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-crosssectional 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 ni , by applying the linear static balance (5) with its internal forces n i , to each tendon

individually ni + n ′ i = 0 3×1 . (21)
Again, as a perfectly flexible string cannot support internal couples or shear forces but only a constant tension τ i tangent to its path, the field of internal forces along each tendon i reads

n i = τ i t i , (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 (τ ′ i = 0 3×1 ), [START_REF] Ha | Modeling tube clearance and bounding the effect of friction in concentric tube robot kinematics[END_REF] provides the expression of the external forces of each tendon

ni = -n ′ i = -τ i t ′ i . (23)
B. Forward kineto-static BVP of a TACR Following the methodology of [START_REF] Rucker | Statics and dynamics of continuum robots with general tendon routing and external loading[END_REF], ( 23) is introduced into ( 19) and [START_REF] Ha | Elastic stability of concentric tube robots subject to external loads[END_REF], and the result into [START_REF] Dupont | Design and control of concentric-tube robots[END_REF] then in [START_REF] Mahl | A variable curvature modeling approach for kinematic control of continuum manipulators[END_REF]. These identifications provide the stress balance equations of the TACR actuated by a set of m tendons

n ′ + m i=1 τ i t ′ i + next = 0 3×1 , c ′ + r ′ × n + m i=1 d i × (τ i t ′ i ) + cext = 0 3×1 . (24)
To obtain an explicit form of the ODEs (24) compatible with [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF], 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

ni = -τ i t ′ i = -τ i r ′ i ∥r ′ i ∥ ′ = τ i r ′ i 2 ∥r ′ i ∥ 3 r ′′ i , (25)
referring to Appendix B of [START_REF] Rucker | Statics and dynamics of continuum robots with general tendon routing and external loading[END_REF] for further details. Next, one needs to expand the first and second derivatives of r i . Expressing these quantities in cross-sectional frame coordinates permits a more concise writing, as follows:

r ′ i = RΓ i = R KD i + D ′ i + Γ , r ′′ i = R Γ ′ i + KΓ i , (26)
where [START_REF] Orekhov | Modeling parallel continuum robots with general intermediate constraints[END_REF] in [START_REF] Amanov | Tendon-driven continuum robots with extensible sections-a model-based evaluation of path-following motions[END_REF] and the result in [START_REF] Rucker | A geometrically exact model for externally loaded concentric-tube continuum robots[END_REF] and [START_REF] Ha | Elastic stability of concentric tube robots subject to external loads[END_REF]

Γ i = R T r ′ i . Now, introducing
yields nact = R (a + AΓ ′ + GK ′ ) , cact = R (b + BΓ ′ + HK ′ ) , (27)
where A, B, G, H, a, and b depend on

Γ i = R T r ′ i = Γ + K × D i + D ′
i and are all proportional to τ i , as detailed in Appendix A. Finally, using [START_REF] Wang | Eccentric tube robots as multiarmed steerable sheaths[END_REF] in [START_REF] Dupont | Design and control of concentric-tube robots[END_REF] and the result in [START_REF] Rucker | Mechanics-based modeling of bending and torsion in active cannulas[END_REF], and rearranging the obtained expressions, yields the ODEs governing the statics of a TACR in the expected explicit form

    r R Γ K     ′ =     RΓ R K H l + A G B H a + H -1 c d     , (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 [START_REF] Rucker | Mechanics-based modeling of bending and torsion in active cannulas[END_REF]. 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

n + = m i=1 -n i (l) = - m i=1 τ i t i , c + = m i=1 d i × (-n i (l)) = - m i=1 d i × τ i t i . (29)
Finally, with these detailed BCs, (28) defines the same BVP as obtained in [START_REF] Rucker | Statics and dynamics of continuum robots with general tendon routing and external loading[END_REF]. 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 [START_REF] Renda | A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation[END_REF], [START_REF] Renda | A sliding-rod variablestrain model for concentric tube robots[END_REF], or by feeding a continuous Newton-Euler inverse algorithm with specific inputs [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF].

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 twostep 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 noncommutative 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×4 . Therefore the configuration space (2) of a clamped Cosserat rod can be alternatively defined by the functional space

C 2 = {ξ : X ∈ [0, l] → ξ(X) ∈ R 6 } . ( 30
)
C 2 is a (functional) linear space, since SE(3) in ( 2) is now replaced by the vector space R 6 . This alternate definition allows to apply usual procedures of linear reduction. According to the Ritz approach [START_REF] Macdonald | Successive approximations by the Rayleigh-Ritz variation method[END_REF], the strain field ϵ = ξ -ξ 0 is decomposed on a functional basis of "strain functions" as

ϵ = Φ(X)q , ( 31
)
where q is a R k vector of generalized strain coordinates and Φ is a R 6×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 (Φq + ξ 0 ) ∧ (32)
from X = 0 where g(0) = 1 4×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 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

δW ext + δW int = 0 , (33)
where δW ext and δ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 δg, with δg(0) = 0 4×4 . As for (4), this field of transformation perturbations is entirely described by a field of twist noted δζ and defined by

δg = gδ ζ , with: δζ(0) = 0 6×1 . (34)
The condition at X = 0 is imposed by the compatibility of variations δg with the geometric BC g(0) = 1 4×4 . When applied to the backbone, such virtual displacements generate variations of its strain field δϵ = δ(ξ -ξ 0 ) = δξ. These relative variations δξ are related to the absolute ones δζ through the commutation relation

δζ ′ = δξ -ad ξ δζ . (35)
To derive this relation, it suffices to remark that, since variations δ do not affect the labels X, the variation and Xderivation are exchangeable, i.e.

δ(g ′ ) = (δg) ′ . ( 36
)
Then, using the definitions of δζ and ξ in terms of g, namely δg = gδ ζ, and g ′ = g ξ in [START_REF] Zienkiewicz | The Finite Element Method: Its Basis and Fundamentals[END_REF], as well as the anti-hat operator, leads to (see Appendix B) the expected relation [START_REF] Goldstein | Classical Mechanics[END_REF]. Based on these preliminaries, one can now develop the balance of virtual works [START_REF] Renda | A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation[END_REF] as follows. The virtual work of internal elastic stress reads

δW int = - l 0 δϵ T Hϵ dX . (37)
As regards the work of external forces, it can be detailed as the sum

δW ext = δ Wext + δW act = l 0 δζ T F dX + δζ(l) T F + + m i=1 δl i τ i , ( 38
)
where δ Wext stands for the virtual work of external gravity and contact forces, while δW act is the virtual work of tensions exerted on the tendons, with δ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)

δl i = l 0
Jli δϵ dX,

with: Jli = 1 ∥Γ i ∥ Γ T i , Γ T i D T i 0 3×3 0 3×3 1 3×3 , (39)
where the Jacobian (continuum) operator Jli is introduced, with

Γ i = R T r ′ i = Γ + K × 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 [START_REF] Antman | Nonlinear Problems of Elasticity[END_REF], and exploiting the fact that τ ′ i = 0 3×1 , allows rewriting the virtual work of actuation in the form

δW act = m i=1 l 0 δϵ T JT li τ i dX = l 0 δϵ T m i=1
JT li τ i dX .

(40) This last equation shows that the effect of tendons naturally appears as stress (and not forces) since they work along strain variations δϵ (and not along configuration variations, i.e. virtual displacements δζ). Based on this remark, the field of actuation stress Λ act is defined as

δW act = l 0 δϵ T Λ act dX ⇒ Λ act = m i=1 τ i JT li = m i=1 1 ∥Γ i ∥ D i × Γ i Γ i τ i , (41)
which can be integrated to the constitutive law [START_REF] Liu | Effect of external and internal loads on tension loss of tendon-driven continuum manipulators[END_REF] yielding the more general active-passive law

Λ = Λ act + Hϵ . (42)
To detail the expression of the first term δ Wext of (38), the virtual displacements δζ must be expressed in terms of δϵ and subsequently δq. To this end, the first step is to introduce ad ξ = Ad -1 g Ad ′ g in the commutation relation [START_REF] Goldstein | Classical Mechanics[END_REF]. The second step is to integrate the result between 0 and X. This provides the expression of any virtual displacement field δζ that is compatible with the internal kinematics and the geometric BCs [START_REF] Kane | Dynamics: Theaory and Applications[END_REF] where J(X) defines the Jacobian matrix, which maps any variation δq to δζ(X). To obtain [START_REF] Kane | Dynamics: Theaory and Applications[END_REF], one needs the relation δζ(0) = 0 6×6 (consequence of g(0) = 1 4×4 ) as well as the variation of ( 31)

δζ(X) = Ad -1 g(X) X 0 Ad g(Y ) δξ(Y ) dY = Ad -1 g(X) X 0 Ad g(Y ) Φ dY δq = J(X)δq ,
δξ = δϵ = Φδq . ( 44
)
Using [START_REF] Kane | Dynamics: Theaory and Applications[END_REF] in δ Wext of [START_REF] Antman | Nonlinear Problems of Elasticity[END_REF], and [START_REF] Orekhov | Solving Cosserat rod models via collocation and the Magnus expansion[END_REF] in δW int of (37) as well as in δW act of (41), the balance of virtual works (33) reads

δq T Q ext + δq T K ϵϵ q = δq T Q act . (45)
(45), being true for any δq, provides the static balance of generalized forces

Q ext + K ϵϵ q = Q act . ( 46
)
This equation introduced the following three vectors of generalized forces. The R k vector of generalized restoring forces

Q rod = l 0 Φ T Hϵ dX = l 0 Φ T HΦ dX q = K ϵϵ q , (47)
with K ϵϵ the R k×k matrix of generalized stiffness. The R k vector of generalized forces of actuation

Q act = - l 0 Φ T Λ act dX = - m i=1 l 0 Φ T ∥Γ i ∥ D i × Γ i Γ i dX τ i = L(q)τ . ( 48
)
And the R k vector of generalized external forces

Q ext = - l 0 J T F dX -J(l) T F + , (49)
which can be alternatively written in the equivalent form (see Appendix D)

Q ext = - l 0 Φ T Λ dX , ( 50
)
where Λ is the solution of the inverse kineto-static BVP of a TACR, obtained by gathering ( 4) and [START_REF] Liu | Effect of external and internal loads on tension loss of tendon-driven continuum manipulators[END_REF],

g Λ ′ = g ξ ad T ξ Λ -F , g(0) = 1 4×4 , Λ(l) = F + . (51)
This new BVP is called the inverse BVP, because it allows to calculate force-outputs (here the stress Λ) from the knowledge of motion-inputs (here the strain ϵ = ξ -ξ 0 ). In other words, Λ is the field of stress that balances the imposed external wrenches F and F + when the backbone is in the configuration ξ.

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 [START_REF] Renda | A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation[END_REF], [START_REF] Renda | A sliding-rod variablestrain model for concentric tube robots[END_REF], 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 [START_REF] Kane | Dynamics: Theaory and Applications[END_REF]. On the other hand, if one uses (50) as in [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF], 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 [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF]), 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 [START_REF] Dupont | Design and control of concentric-tube robots[END_REF] and recall that, in the Newtonian approach (see Section IV), the action of the tendons is modeled as external forces nact and couples cact . Both are integrated in the wrench of external distributed loads F applied along the backbone as follows:

F = C N = Cext Next + Cact Nact = R T cext R T next + R T cact R T nact .
(52) Whereas in the Lagrangian approach (see Section V), the same actions are modeled by the wrench of stress across the backbone Λ act whose expression is defined by [START_REF] Lynch | Modern Robotics[END_REF], and which can be expressed in the inertial frame as

c act n act = RC act RN act = m i=1 d i × t i t i τ i , (53)
where remind that

RΓ 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 [START_REF] Rucker | Mechanics-based modeling of bending and torsion in active cannulas[END_REF] 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), [START_REF] Ha | Elastic stability of concentric tube robots subject to external loads[END_REF], and ( 23), but as the field of stress wrench across the backbone (53). This change of view point means that, in [START_REF] Rucker | Mechanics-based modeling of bending and torsion in active cannulas[END_REF], cact and nact are removed from the model of external loads [START_REF] Dupont | Design and control of concentric-tube robots[END_REF], 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 , n = n rod + n act . (54)
Note that once expressed in the cross-sectional frames, these two relations do define the active constitutive law [START_REF] Macdonald | Successive approximations by the Rayleigh-Ritz variation method[END_REF], which takes the detailed form

Λ = Λ act + Hϵ = R T c act R T n act + H a (K -K 0 ) H l (Γ -Γ 0 ) , (55)
where n act and c act are given by (53). As announced, substituting (53) and (54) in the stress balance (5), produces

n rod + m i=1 τ i t i ′ + next = 0 3×1 , c rod + m i=1 d i × (τ i t i ) ′ + r ′ × n rod + m i=1 τ i t i + cext = 0 3×1 . (56)
Rearranging (56) with the usual composition rules of derivatives and using d i = r i -r and τ ′ i = 0 3×1 , yields

n ′ rod + m i=1 τ i t ′ i + next = 0 3×1 , c ′ rod + r ′ × n rod + m i=1 d i × (τ i t ′ i ) + m i=1 (r ′ + d ′ i ) × (τ i t i ) + cext = 0 3×1 . (57)
Finally, using the relations r ′ + d ′ i = r ′ i = ∥r ′ i ∥ t i and t i × t i = 0 3×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 [START_REF] Black | Parallel continuum robots: Modeling, analysis, and actuation-based force sensing[END_REF]. 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 [START_REF] Black | Parallel continuum robots: Modeling, analysis, and actuation-based force sensing[END_REF] with the BCs of ( 17), with (29), of the Newtonian approach, is equivalent to the closed continuous formulation [START_REF] Featherstone | Rigid Body Dynamics Algorithms[END_REF], with ( 41), [START_REF] Macdonald | Successive approximations by the Rayleigh-Ritz variation method[END_REF], 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 [START_REF] Black | Parallel continuum robots: Modeling, analysis, and actuation-based force sensing[END_REF]. 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 [START_REF] Black | Parallel continuum robots: Modeling, analysis, and actuation-based force sensing[END_REF] with the shooting method consists in finding the unknown proximal BCs Γ(0) and K(0) such that the known distal BCs Γ(l) and K(l), defined by [START_REF] Rucker | Mechanics-based modeling of bending and torsion in active cannulas[END_REF], are fulfilled (see bottom left area of Fig. 2). This search is achieved iteratively by applying a root finding algorithm to the residual vector

R(Γ ⋆ (0), K ⋆ (0)) = Γ(l) -Γ ⋆ (l) K(l) -K ⋆ (l) , (58)
where Γ ⋆ (l) and K ⋆ (l) are the distal values of Γ and K obtained by forward integrating the ODEs of (28) from the known proximal BCs (r, R)(0) -supplemented with the guessed values (Γ ⋆ , K ⋆ )(0) of the missing ones -to X = l. This resolution proposes to iterate over the values of Γ(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 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) × R 3 or SE(3) can be used [START_REF] Renda | A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation[END_REF], [START_REF] Orekhov | Solving Cosserat rod models via collocation and the Magnus expansion[END_REF].

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 [START_REF] Bretl | Quasi-static manipulation of a Kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF].

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

t i j = lim ε→0 t i (X + j ) -t i (X - j ) , (59)
where X ± j = X j ± ε, 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 f 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

n j = τ i t i j , c j = d i × n j . ( 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 [START_REF] Black | Parallel continuum robots: Modeling, analysis, and actuation-based force sensing[END_REF], each time a tendon slope discontinuity is encountered. The slope discontinuity t i j can be obtained with the relation Γ i / ∥Γ i ∥ j = R T t i j . Hence, (60) is rewritten in its cross-sectional frame version

N j = τ i Γ i ∥Γ i ∥ j , C j = D i × N j , (61)
where Γ ± i is computed similarly to [START_REF] Orekhov | Modeling parallel continuum robots with general intermediate constraints[END_REF], with

Γ i (X ± j ) = Γ(X ± j ) + K(X ± j )D i (X j ) + D ′ i (X ± j ) . ( 62
)
In practice, this calculation is achieved by using the following jump relations on the space-rates:

Γ K + = Γ K - + H l + A G B H a + H -1 - c d j ,
(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×1 and D ′ i j ̸ = 0 3×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 [START_REF] Featherstone | Rigid Body Dynamics Algorithms[END_REF]. Note that this is different from [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF], where the resolution is performed by the explicit time integration of an overdamped equivalent system. Also, compared to [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF]- [START_REF] Renda | A sliding-rod variablestrain model for concentric tube robots[END_REF] 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 [START_REF] Featherstone | Rigid Body Dynamics Algorithms[END_REF], for any imposed vector of tendon tensions τ = (τ 1 .. τ 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 R that takes the generic form

R(q, τ ) = 0 k×1 , with: R(q, τ ) = Q ext +K ϵϵ q -Q act . ( 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)R(q, τ ) , (65)
where J(q) = ∂R ∂q denotes the R k×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 R and its Jacobian matrix J.

To numerically compute the residual vector R from any qvector, the expressions ( 47), [START_REF] Neumann | Considerations for follow-theleader motion of extensible tendon-driven continuum robots[END_REF], and (50), are added to the inverse BVP (51) through the state variable

y(X) = - l X Φ T (Λ act -Λ + HΦq) dY (66)
such that y(0) = R. The state variables (g, Λ, y) are then governed by the augmented inverse BVP, which is fed with ξ = ξ 0 + Φq, and reads

  g Λ y   ′ =   g (ξ 0 + Φq) ∧ ad T (ξ0+Φq) Λ - F Φ T (Λ act -Λ + HΦq)   , g(0) = 1 4×4 , Λ(l) = F + , y(l) = 0 k×1 , (67)
where Λ act is a function of q and τ , given by ( 41), while F depends on the context. For instance, if the robot is only subject to gravity, F = 0 1×3 , R T a T g Aρ T , with a g the acceleration gravity field expressed in the inertial frame, A the area of the robot's section, and ρ 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 [START_REF] Black | Parallel continuum robots: Modeling, analysis, and actuation-based force sensing[END_REF]. 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 Λ and y producing R(q) = y(0). This method has been proposed and interpreted as a Newton-Euler computed torque algorithm in [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF], similar to those developed for rigid multi-body systems [START_REF] Featherstone | Rigid Body Dynamics Algorithms[END_REF]. 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, Λ, 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 [START_REF] Trefethen | Spectral methods in MATLAB[END_REF].

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 ∆q as

  ∆ζ ∆Λ ∆y   ′ =   -ad (ξ0+Φq) ∆ζ + Φ∆q ad T Φ∆q Λ + ad T (ξ0+Φq) ∆Λ -∆ F Φ T (∆Λ act -∆Λ + HΦ∆q)   , ∆ζ(0) = 0 6×1 , ∆Λ(l) = ∆F + , ∆y(l) = 0 k×1 , (68)
where ∆ζ is the differential of g in se(3) defined by ∆ζ = g -1 g ∨ , while ∆Λ act , ∆ F , and ∆F + are the differentials of the model of actuation stress and external forces, respectively. Note here that, the variation ∆ being only generated by the variation of the configuration ∆q, ∆τ = 0 m×1 in ∆Λ act . Moreover, the dependence of Λ act on q is very weak and may be neglected in the calculation of the Jacobian. Thus, ∆Λ act = 0 6×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 ∆y(0) = ∆R = J(q)∆q, applying the unit input vector ∆q = δ α to (68), where δ α is a R k×k vector of zeros, except the entry α, which is equal to 1, (68) yields the column α of J. Hence, repeating the process for α = 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 [START_REF] Rucker | Mechanics-based modeling of bending and torsion in active cannulas[END_REF].

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 χ 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

χ j=1 δ Wext,j + δW act,j + δW int,j = 0 , (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 [START_REF] Antman | Nonlinear Problems of Elasticity[END_REF]. 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 (χ) ) 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 δq,

δq T (Q ext + K ϵϵ q -Q act ) = 0 k×1 , (70)
where [START_REF] Burgner-Kahrs | Continuum robots for medical applications: A survey[END_REF] .. Q ext,(χ) ) T , and Q act = (Q act,(1) .. Q act,(χ) ) T . In other words, the multi-segment TACR model is simply obtained by appending the matrices and vectors of each segment one after another.

K ϵϵ = diag(K ϵϵ,(1) .. K ϵϵ,(χ) ), Q ext = (Q ext,
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 χ BVPs and TBVPs, are now connected through their BCs as follows:

g j (1) = g j+1 (0) , Λ j (1) = Λ j+1 (0) , (71) ∆ζ j (1) = ∆ζ j+1 (0) , ∆Λ j (1) = ∆Λ j+1 (0) . (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×4 and Λ χ (1) = F + , respectively.

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 [START_REF] Rucker | Statics and dynamics of continuum robots with general tendon routing and external loading[END_REF], 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 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.

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 [START_REF] Oliver-Butler | Continuum robot stiffness under external loads and prescribed tendon displacements[END_REF], [START_REF] Laschi | Soft robot arm inspired by the octopus[END_REF]. 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, ±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 0, cos 2π l X + iπ , sin 2π l X + iπ 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 [START_REF] Rucker | Statics and dynamics of continuum robots with general tendon routing and external loading[END_REF], [START_REF] Starke | On the merits of helical tendon routing in continuum robots[END_REF]. 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 II. Full and dashed lines represent simulations obtained using the Newtonian approach and the Lagrangian approach, respectively. the backbone. Tendons (1,2,3) terminate in the first segment, tendons [START_REF] Chikhaoui | Chapter 8 -modeling and control strategies for flexible devices[END_REF][START_REF] Mahl | A variable curvature modeling approach for kinematic control of continuum manipulators[END_REF][START_REF] Santina | Control oriented modeling of soft robots: The polynomial curvature case[END_REF] in the second segment, and tendons [START_REF] Jones | Kinematics for multisection continuum robots[END_REF][START_REF]Practical kinematics for real-time implementation of continuum robots[END_REF][START_REF] Neppalli | Closedform inverse kinematics for continuum manipulators[END_REF] 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 [START_REF] Jones | Kinematics for multisection continuum robots[END_REF], [START_REF] Gonthina | Mechanics for tendon actuated multisection continuum arms[END_REF], [START_REF] Amanov | Tendon-driven continuum robots with extensible sections-a model-based evaluation of path-following motions[END_REF]. 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 φ j (X) (0, cos (θ j ) , sin (θ j )) T , with φ j (X) =      1 for X < (j -1) l j 1 - X-(j-1)lj lj for (j -1) l j ≤ X ≤ j l j 0 for X > j l j ,
where D i,j (X), stands for the routing path of tendon i = 1 .. 6 that ends in segment j = 1 .. 3 and θ j = (j -1) 2π/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 [START_REF] Oliver-Butler | Continuum robot stiffness under external loads and prescribed tendon displacements[END_REF], 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. II. Full and dashed lines represent simulations obtained using the Newtonian approach and the Lagrangian approach, respectively.

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 [START_REF] Rucker | Statics and dynamics of continuum robots with general tendon routing and external loading[END_REF]. 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 [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF]. 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 [START_REF] Black | Parallel continuum robots: Modeling, analysis, and actuation-based force sensing[END_REF], which is integrated from X = 0 to X = l, whereas in the Lagrangian approach, the distributed loads F are defined over X ∈ ]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 [START_REF] Starke | On the merits of helical tendon routing in continuum robots[END_REF].

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 [START_REF] Renda | A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation[END_REF], 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 convergence, 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 generally 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 [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF]- [START_REF] Renda | A sliding-rod variablestrain model for concentric tube robots[END_REF] 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 problems iteratively. And even more so, considering the fact that establishing the generalized stiffness matrix K ϵϵ , the most costly operation, is a constant function of the shape functions Φ 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 [START_REF] Rucker | Statics and dynamics of continuum robots with general tendon routing and external loading[END_REF], [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF]. 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 [START_REF] Rucker | Statics and dynamics of continuum robots with general tendon routing and external loading[END_REF], [START_REF] Starke | On the merits of helical tendon routing in continuum robots[END_REF], [START_REF] Neumann | Considerations for follow-theleader motion of extensible tendon-driven continuum robots[END_REF].

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. [START_REF] Wang | Eccentric tube robots as multiarmed steerable sheaths[END_REF] AND [START_REF] Black | Parallel continuum robots: Modeling, analysis, and actuation-based force sensing[END_REF] The ODEs ( 27) make use of the notations

APPENDIX A DETAILED EXPRESSIONS AND NOTATIONS USED IN

A i = -τ i Γ i 2 ∥Γ i ∥ 3 , A = m i=1 A i , B i = D i A i , B = m i=1 B i , G = - m i=1 A i D i = B T , H = - m i=1 B i D i , a i = A i KΓ i + KD ′ i + D ′′ i , a = m i=1 a i , b i = D i a i , b = m i=1 b i .
The BVP (28) makes use of the above notations, as well as

c = H l Γ ′ 0 -KH l (Γ -Γ 0 ) -R T next -a , d = H a K ′ 0 -KH a (K -K 0 ) -ΓH l (Γ -Γ 0 ) -R T cext -b .

APPENDIX B DERIVATION OF THE COMMUTATION RELATION

Let us recall the definitions δg = gδ ζ and g ′ = g ξ of ( 4) and [START_REF] Renda | A sliding-rod variablestrain model for concentric tube robots[END_REF], respectively. To derive [START_REF] Goldstein | Classical Mechanics[END_REF], these definitions are introduced in (36) as follows:

Figure 3 .

 3 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 Rt. 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.

Figure 4 .

 4 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 TableII. Full and dashed lines represent simulations obtained using the Newtonian approach and the Lagrangian approach, respectively.

Figure 5 .

 5 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 TableII. Full and dashed lines represent simulations obtained using the Newtonian approach and the Lagrangian approach, respectively.

  δ g ξ = gδ ζ ′ ⇔ δg ξ + gδ ξ = g ′ δ ζ + g δ ζ ′ . (73)Using the same definitions, once more, leads tog δ ζ ξ + δ ξ = g ξδ ζ + δ ζ ′ ,(74)which, holding for any g, yieldsδ ξ = δ ζ ′ + ξδ ζ -δ ζ ξ = δ ζ ′ + [ ξ, δ ζ ] ,(75)where [•, •] denotes the usual commutator of matrices. Finally, using the anti-hat (• ∨ ) operation, one can show that [ ξ, δ ζ ] ∨ = ad ξ δζ, and get the expected commutation relation in terms of twists[START_REF] Goldstein | Classical Mechanics[END_REF].

Table I COMPARISON

 I OF THE NOTATION CONVENTIONS OF THE LITERATURE.

Table II GEOMETRY

 II AND MATERIAL PARAMETERS OF THE SIMULATED ROBOT.

	Parameter Name	Parameter Symbol Value
	Backbone diameter	R b		0.4 mm
	Backbone length (A, B, C)	l		242 mm
	Backbone length (D, E)	l = l j × 3	300 mm
	Tendon offset	Rt		8 mm
	Gravity constant	∥ag∥	9.81 N/kg
	Robot equivalent density 2	ρ =	0.47 N/m ∥ag ∥π(R b ) 2	95 • 10 3 kg/m 3
	Young's modulus	E		210 GPa
	Poisson's ratio	ν		0.3125

Table III SYNTHESIS

 III OF RESULTS FOR THE FIVE SIMULATED SCENARIOS. * The letters of the scenarios correspond to the subsections in Section VIII. *** Differences between the results that are obtained with both approaches.

	Scenario**		A	B	C	D	E
	Robot characteristics	Type of routing	Parallel	Convergent	Helical	Parallel	Convergent
		Number of segments	1	1	1	3	3
		Robot length [mm]	242	242	242	300	300
	Tip position difference*** [mm] Median value	5.43 • 10 -4	9.40 • 10 -3	1.42 • 10 -1	3.42 • 10 -1	6.93 • 10 -2
		Lower bound 95% confidence interval 2.41 • 10 -9	8.77 • 10 -15	4.48 • 10 -3	1.00 • 10 -5	1.11 • 10 -2
		Upper bound 95% confidence interval	3.13 • 10 -3	2.32 • 10 -2	7.21 • 10 -1	1.29	2.76 • 10 -1
	Tip orientation difference*** [ • ]	Median value	3.26 • 10 -4	1.10 • 10 -3	1.99 • 10 -1	5.56 • 10 -2	1.32 • 10 -2
		Lower bound 95% confidence interval 1.03 • 10 -7	5.42 • 10 -8	8.09 • 10 -3	3.39 • 10 -7	2.78 • 10 -3
		Upper bound 95% confidence interval	1.65 • 10 -3	3.20 • 10 -3	5.74 • 10 -1	2.14 • 10 -1	4.74 • 10 -2

*

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'.

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[START_REF] Rucker | Statics and dynamics of continuum robots with general tendon routing and external loading[END_REF].

The algorithm parameters were chosen empirically during preliminary experiments such that the simulation results had converged for each approach.

This work was supported by grants ANR-11-LABX-0004-01, ANR-19-P3IA-0003, ANR-20-CE33-0001, ANR-10-IAHU-0002, and ANR-18-CE19-0012.

APPENDIX C

PROOF OF [START_REF] Marsden | Introduction to mechanics and symmetry[END_REF] The expression of the length of a tendon

leads to

In order to rewrite the last term of (77), recall that r ′ i = RΓ i , δr ′ i = δ(RΓ i ), and thus

Now remark that, since R ∈ SO(3), there always exists a vector δψ ∈ R 3 such that R T δR = δ ψ. As a result

Finally, introducing the expression Γ i = Γ + K × D i + D ′ i into (78) and the result of this substitution into (77) provides, after factorization of δϵ = (δK T , δΓ T ) T , the expected expression [START_REF] Marsden | Introduction to mechanics and symmetry[END_REF].

APPENDIX D PROOF OF (50)

The relation Ad ′ g = Ad g ad ξ , leads to

Then, introducing this relation in the balance [START_REF] Camarillo | Mechanics modeling of tendon-driven continuum manipulators[END_REF], one can integrate it, and obtain the balance of wrenches

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

which is no more than the expression of Q ext given by (49).