# Magnetic concentric tube robots: introduction and analysis

## Journal Article

### Author(s):

Peyron, Quentin; Boehler, Quentin ID; Rougeot, Patrick; Roux, Pierre; Nelson, Bradley ID; Andreff, Nicolas; Rabenorosoa, Kanty; Renaud, Pierre

### Publication date:

2022-04

### Permanent link:

https://doi.org/10.3929/ethz-b-000543633

### Rights / license:

Creative Commons Attribution 4.0 International

### Originally published in:

The International Journal of Robotics Research 41(4), https://doi.org/10.1177/02783649211071113

### Funding acknowledgement:

743217 - Soft Micro Robotics (EC)

185039 - Arbeitstitel "Soft Magnetic Robots: Modeling, Design and Control of Magnetically Guided Continuum Manipulators" (SNF)
# Magnetic concentric tube robots: introduction and analysis

Quentin Peyron<fnref n="1" />,<fnref n="2" /> ID, Quentin Boehler<fnref n="3" /> ID, Patrick Rougeot<fnref n="2" />, Pierre Roux<fnref n="2" />, Bradley J. Nelson<fnref n="3" />, Nicolas Andreff<fnref n="2" />, Kanty Rabenorosoa<fnref n="2" /> and Pierre Renaud<fnref n="1" />

The International Journal of
Robotics Research
2022, Vol. 0(0) 1–23
© The Author(s) 2022
Article reuse guidelines:
sagepub.com/journals-permissions
DOI: 10.1177/02783649211071113
journals.sagepub.com/home/ijr

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

SAGE

## Abstract

In this paper, we propose a new type of continuum robot, referred to as a magnetic concentric tube robot (M-CTR), for performing minimally invasive surgery in narrow and difficult-to-access areas. The robot combines concentric tubes and magnetic actuation to benefit from the 'follow the leader' behaviour, the dexterity and stability of existing robots, while targeting millimetre-sized external diameters. These three kinematic properties are assessed through numerical and experimental studies performed on a prototype of a M-CTR. They are performed with general forward and inverse kineto-static models of the robot, continuation and bifurcation analysis, and a specific experimental setup. The prototype presents unique capabilities in terms of deployment and active stability management, while its dexterity in terms of tip orientability is also among the best reported for other robots at its scale.

## Keywords

Continuum robot, concentric tube robot, magnetic continuum robot, follow-the-leader deployment, orientability, stability analysis, medical robot

## 1. Introduction

### 1.1. Continuum robots for minimally invasive surgery

Continuum robots are frequently considered for minimally invasive surgery (MIS) due to their kinematic advantages (Burgner-Kahrs et al., 2015). They usually consist of a slender elastic backbone that is deformed by actuation torques or forces, in order to control their shape and tip pose (see Figure 1). In particular, it is possible to deploy continuum robots with a 'follow the leader' (FTL) behaviour. The robot backbone can be actuated so that its shape corresponds to the path followed by the tip without relying on contacts with its environment (Garriga-Casanovas and Rodriguez y Baena, 2018; Palmer et al., 2014). As a result, continuum robots are well-suited for deployment of medical tools through a patient's anatomy while avoiding sensitive areas, followed by precise manipulation at the operation site. They have successfully been considered in a number of interventions including vascular surgery, cardiac surgery (Gosline et al., 2012), neurology (Burgner et al., 2013), otolaryngology (Simaan et al., 2009), foetal (Dwyer et al., 2017) and abdominal interventions (Burdette et al., 2010).

MIS is constantly evolving towards interventions in smaller and more difficult to access operation sites. Two examples of interventions showing this trend are olfactory cells inspection (Girerd et al., 2018), and middle-ear

interventions (Dahroug et al., 2018; Fichera et al., 2017). They involve manipulating medical tools in millimetre-scale cavities along several tens of millimetre, inside the nose and the ears, for tissue manipulation and imaging. In particular, such tools must be deployed through tortuous and narrow orifices, while ensuring the integrity of sensitive wall tissues. For example, the olfactory cleft is covered with fragile neural receptors (Kavoi and Jameela, 2011) which would be damaged in case of contact with the tool, resulting in a degraded olfactory sense. In the case of middle-ear surgery, the tympanic membrane and facial nerves can be damaged, resulting in hearing losses and facial distortion (Olszewska et al., 2004). The tools must be manipulated at the operation site while complying with the constraints imposed by the patient's anatomy, which requires dexterity. Being able to rotate tools with a large angular displacement around the operation site is an important feature in tissue manipulation (Simaan et al., 2009) and imaging (Goldman

### Corresponding author:

Quentin Peyron, ICube/AVR and FEMTO-ST Institute/AS2M, Besançon 25044, France.

Email: quentin.peyron@inria.fr

<footnote n="1">¹Cube, University of Strasbourg, CNRS, INSA Strasbourg, France</footnote>

<footnote n="2">²FEMTO-ST Institute/AS2M, Univ. Bourgogne Franche-Comte/CNRS, Besançon, France</footnote>

<footnote n="3">³ETH Zurich/Multi-Scale Robotics Lab, Zürich, Switzerland</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic view introducing M-CTRs. The blue and orange parts are composed, respectively, of Nitinol and soft elastomer materials. Red parts designate magnetic elements. Proximal actuation includes rotation and translation of each tube. *B* represents the external magnetic field.</figcaption></figure>

et al., 2013). Limited dexterity may imply several complications such as partial resection and thus regrowth of unhealthy tissues, such as cholesteatoma (Olszewska et al., 2004), and wrong diagnostics due to partial imaging (Mowatt et al., 2011).

The evolution of MIS involves four kinds of requirements for the design of continuum robots. First, there is a need for robots with millimetre diameters, in order to navigate through narrow environments. Second, FTL behaviour is required to access the operation site and to manage the safety throughout the tool's deployment. Third, high dexterity is required in term of tip orientation capability, which has been referred as ‘orientability’ in the literature (Li et al., 2017; Wu et al., 2017). Fourth, performing FTL deployments and tip orientation control requires the robot to be stable during its motion, in particular in the vicinity of sensitive tissues. Many other requirements can be relevant for the design of continuum robots for minimally invasive surgery. Load bearing capabilities and stiffness properties are of importance in applications involving tissue resection and manipulation for example. Medical devices must also comply to requirements in terms of sterilization, usability and integration to the surgical workflow. In this paper, we focus on the four requirements detailed above, demonstrating an interesting concept of continuum robot for minimally invasive surgery rather than providing a ready-to-use and specific medical tool.

## 1.2. Magnetic concentric tube robot

Continuum robots have been classified (Burgner-Kahrs et al., 2015) according to the actuation strategy used to deform the backbone. Each actuation strategy has its own advantages and drawbacks, and it confers to the robot different performances and features as shown in Table 1. Three criteria are used according to the requirements we introduced. The first criterion is the capability of achieving

FTL behaviour during deployment. The behaviour can be perfect (P); that is, the robot backbone follows the tip path exactly, or approximated (A), depending on the intrinsic kinematic properties of the robot. The second criterion is the highest angular displacement that can be achieved by the robot tip, which is directly related to orientability (Li et al., 2017). The third criterion is the lowest external diameter achieved among the existing prototypes.

Amongst the intrinsic and extrinsic actuation strategies, only the concentric tube (CT) and tendon driven (TD) actuation have been considered to perform FTL deployments, the last one not being capable of perfect deployments. The angular displacements seem limited for the most part to ±120°, except for robots considering TD and magnetic (M) actuation. External diameters below 1 mm were only achieved with robots using magnetic and electro-active polymer (EAP) actuation. Consequently, there is no design of continuum robot that can achieve simultaneously perfect FTL deployments, over 120° angular displacement and have a millimetre diameter with a single actuation strategy. Considering hybrid strategies seems promising to obtain new robot design with improved performances. In particular, the hybrid design that combines tendon-driven actuation and concentric tubes seems to provide interesting compromises. However, the prototypes presented in Amanov et al. (2017) and Swaney et al. (2016) cannot satisfy the three requirements simultaneously. In particular, performing perfect FTL behaviour requires to assemble two robots concentrically while integrating tendons, making its miniaturization challenging. In this work, we introduce and investigate the use of concentric tubes with magnetic actuation in order to benefit from the combined advantages of both concentric tube robots (CTR) and magnetic continuum robots (M-CR).

A CTR consists of a telescopic assembly of pre-curved elastic tubes (Webster et al., 2006; Dupont et al., 2010) that are usually made out of Nitinol (NiTi). When assembled, the tube interactions create forces and torques along the backbone, which deform the robot. The tubes are rotated and translated at their base using actuators to control the robot shape and tip pose. CTRs are interesting because they can be deployed in a perfect FTL manner for specific shapes and configurations of the tubes (Gilbert et al., 2015; Garriga-Casanovas and Rodriguez y Baena, 2018) and because they can be fabricated with millimetre external diameters due to the simple assembly of the tubes. Thanks to these two features, they have been proven to be useful for deployment through narrow and sensitive areas such as the olfactory cleft (Girerd et al., 2018) and the brain (Comber et al., 2016). However, CTRs also exhibit complex kinematic behaviour that can make them cumbersome to control. Most notably, they suffer from elastic instabilities due to the torsion of the tubes (Dupont et al., 2010; Gilbert et al., 2016; Ha et al., 2016). For given rotation and translation of the tubes, the robot can have multiple configurations and may snap from one configuration to another. These unstable phenomena are usually avoided in order to ensure patient
<table><caption>**Table 1.** Performances of existing classes of continuum robots in terms of FTL deployment, angular displacement and external diameter.</caption><thead><tr><th>Actuation strategies</th><th>FTL deployment<sup>1</sup></th><th>Angular displacement (degree)<sup>2</sup></th><th>External diameter (mm)<sup>3</sup></th><th>References</th></tr></thead><tbody><tr><td rowspan="3">Extrinsic</td><td>Tendon-driven (TD)</td><td>A</td><td>180</td><td>1.6</td><td>Amanov et al. (2019)<sup>1</sup>Fichera et al. (2017)<sup>2,3</sup></td></tr><tr><td>Multi-backbone (MB)</td><td>—</td><td>± 90</td><td>4.2</td><td>Simaan et al. (2009)</td></tr><tr><td>Concentric tube (CT)</td><td>P</td><td>—</td><td><strong>1.14</strong></td><td>Comber et al. (2016)</td></tr><tr><td rowspan="4">Intrinsic</td><td>Fluidic (F)</td><td>—</td><td>± 120</td><td>4.9</td><td>Bailly and Amirat (2005)</td></tr><tr><td>Shape memory alloy (SMA)</td><td>—</td><td>± 90</td><td>1.67</td><td>Couture and Szewczyk (2018)<sup>2</sup>Jayender et al. (2009)<sup>3</sup></td></tr><tr><td>Electro-active polymer (EAP)</td><td>—</td><td>± 93</td><td>0.95</td><td>Chikhaoui et al. (2018)<sup>2</sup>Farajollahi et al. (2016)<sup>3</sup></td></tr><tr><td>Magnetic (M)</td><td>—</td><td>± 150</td><td><strong>0.8</strong></td><td>Chautems et al. (2018)<sup>2</sup>Charreyron et al. (2019)<sup>3</sup></td></tr><tr><td rowspan="4">Hybrid</td><td>TD-CT</td><td>P</td><td>± 164</td><td>2</td><td>Amanov et al. (2017)<sup>1</sup>Swaney et al. (2016)<sup>2,3</sup></td></tr><tr><td>MB-CT</td><td>P</td><td>± 90</td><td>30</td><td>Kang et al. (2016)</td></tr><tr><td>EAP-CT</td><td>P</td><td>± 20</td><td>1</td><td>Chikhaoui et al. (2018)</td></tr><tr><td>F-TD</td><td>A</td><td>± 90</td><td>7.5</td><td>Kundrat et al. (2016)<sup>1,3</sup>Maghooa et al. (2015)<sup>2</sup></td></tr></tbody></table>

The denomination of the classes is inspired from Burgner-Kahrs et al. (2015). A and P indicates approximate and perfect FTL deployment, respectively. The values correspond to the best results reported in the literature and may come as a result from different papers. Indices 1, 2 and 3 are used to indicate which paper the values are from.

safety. It is realized through careful tube design, which implies limiting the tube pre-curvature (Kim et al., 2014; Bergeles et al., 2015; Ha et al., 2017) and thus the tip angular displacement. These phenomena have recently been considered as useful for tasks requiring a minimal amount of force such as suturing (Riojas et al., 2018). The high dynamics induced by the snapping behaviour are then used to generate tip forces that could not be produced in a quasi-static fashion due to the robot's flexibility.

M-CRs, also called magnetic catheters in the literature, consist in a flexible backbone along which magnetic elements are fixed (Edelmann et al., 2017; Tunay, 2004). An external magnetic field is generated using an electromagnetic navigation system (eMNS), in order to generate torques and forces at specific locations along the backbone. Controlling the magnitude, the orientation, and the gradient of the magnetic field alters the robot's shape and the tip's pose. Due to low bending stiffness and the application of torques directly on the backbone, M-CRs can achieve high angular displacements during magnetic actuation (Chautems et al., 2018). M-CRs can also be fabricated at sub-millimetre scales when passive elements are used (Charreyron et al., 2019). Achieving FTL behaviour with M-CRs has never been considered to the best of our knowledge.

The continuum robot with hybrid actuation we introduce and analyze in this paper is designated as a magnetic concentric tube robot (M-CTR). It is designed to combine the properties of CTRs and M-CRs. As shown on Figure 1, M-CTRs consist of a telescopic assembly of tubes on which one or several passive magnetic elements are fixed.

The tubes can be pre-curved or initially straight and are composed of materials with different stiffness properties, such as NiTi or silicone. M-CTRs are actuated by rotating and translating the tubes and by varying the magnetic field applied on the robot. This concept aims at obtaining robots with millimetre diameter and with FTL deployment capabilities and orientability that could not be achieved using CTR and M-CR technologies alone. In addition, the M-CTR may possess unique abilities in its interaction with the environment. Indeed, magnetic actuation has an effect on concentric tube elastic stability. Applying external forces and torques on a CTR can stabilize or destabilize the robot as shown in Ha et al. (2016). As a result, active stability management may be possible, to ensure patient safety while generating higher tip forces.

## 1.3. Contributions

In this paper, the concept of M-CTRs was investigated throughout three studies, each leading to different contributions. First of all, FTL deployment capability of the hybrid continuum robot was assessed. As it requires the use of magnetic actuation to achieve FTL capability for the first time, a deployment strategy is proposed and validated numerically. Second, the robot orientability was evaluated. We especially consider a scenario inspired from tissue manipulation and imaging, where the robot tip must be rotated at the operation site about axes fixed with respect to the target. This leads us to evaluate orientability numerically while considering physical phenomena such as tube torsion and deformation due to gravity and magnetic actuation for
<table><caption>**Table 2.** Abbreviations and notations.</caption><thead><tr><th></th><th>Abbreviation</th></tr></thead><tbody><tr><td></td><td><strong>Variables</strong></td></tr><tr><td>MIS</td><td>Minimally invasive surgery</td></tr><tr><td>FTL</td><td>Follow the leader</td></tr><tr><td>CTR</td><td>Concentric tube robot</td></tr><tr><td>M-CR</td><td>Magnetic continuum robot</td></tr><tr><td>M-CTR</td><td>Magnetic concentric tube robot</td></tr><tr><td>eMNS</td><td>Electromagnetic navigation system</td></tr><tr><td>NiTi</td><td>Nitinol</td></tr><tr><td>d.o.f.</td><td>Degree of freedom</td></tr><tr><td></td><td><strong>Variables</strong></td></tr><tr><td>i</td><td>Index of tube</td></tr><tr><td>j</td><td>Index of section</td></tr><tr><td>n</td><td>Number of tubes</td></tr><tr><td>n<sub>j</sub></td><td>Number of tubes along section j</td></tr><tr><td>s</td><td>Arc-length of the backbone</td></tr><tr><td>L</td><td>Total length of the robot</td></tr><tr><td>ΔL<sub>jk</sub></td><td>Length of the subsection k of section j</td></tr><tr><td>α<sub>i</sub></td><td>Rotation angle of tube i</td></tr><tr><td><b>B</b></td><td>Magnetic field</td></tr><tr><td><b>G</b></td><td>Vector built with the independent components of the magnetic field gradient</td></tr><tr><td>R<sub>0</sub>, R<sub>b</sub>, R<sub>i</sub></td><td>World, backbone and tube i frames</td></tr><tr><td><sup>0</sup>R<sub>b</sub></td><td>Rotation matrix from R<sub>b</sub> to R<sub>0</sub></td></tr><tr><td><b>p</b></td><td>Cartesian position of the backbone</td></tr><tr><td><b>Q</b></td><td> Quaternion vector representing the backbone orientation</td></tr><tr><td><b>u</b></td><td>Darboux vector representing the backbone curvature</td></tr><tr><td>û<sub>i</sub></td><td>Darboux vector representing the pre-curvature of tube i</td></tr><tr><td>θ<sub>i</sub></td><td>Torsion angle of tube i</td></tr><tr><td>k<sub>bi</sub>, k<sub>ti</sub></td><td>Bending and torsional stiffness of tube i</td></tr><tr><td><b>τ</b><sub>m</sub><b>f</b><sub>m</sub></td><td>Magnetic torques and forces</td></tr><tr><td><b>m</b></td><td>Dipole moment</td></tr><tr><td><b>f</b><sub>g</sub></td><td>Gravity forces</td></tr><tr><td>ρ<sub>i</sub></td><td>Linear mass of tube i</td></tr><tr><td>g</td><td>Gravity acceleration</td></tr><tr><td>λ<sub>p</sub>λ<sub>q</sub></td><td>Lagrange multipliers</td></tr><tr><td><b>x</b></td><td>State space vector as a function of s</td></tr><tr><td><b>q</b></td><td>Actuation inputs vector</td></tr><tr><td><i>N</i></td><td>Number of nodes used to discretize the robot backbone</td></tr><tr><td><i>k</i></td><td>index of node</td></tr><tr><td><b>x</b><sub>k</sub></td><td>State space vector evaluated at node k</td></tr><tr><td><b>p</b><sub>N</sub>, <b>Q</b><sub>N</sub></td><td>Position and orientation of the robot tip</td></tr><tr><td><b>e</b><sub>3</sub></td><td>Unitary vector e<sub>3</sub> = [0 0 1]<sup>T</sup></td></tr><tr><td>φ</td><td>Relative angle between the magnet direction and the magnetic field in the (y<sub>0</sub>, z<sub>0</sub>) plane</td></tr><tr><td>φ<sub>M</sub></td><td>Maximum angle φ for which the robot stays stable</td></tr><tr><td></td><td><strong>Operators</strong></td></tr><tr><td>a<sub>x</sub>, a<sub>y</sub>, a<sub>z</sub></td><td>Components of vector <b>a</b>, <b>a</b> = [a<sub>x</sub> a<sub>y</sub> a<sub>z</sub>]<sup>T</sup></td></tr><tr><td><b>a</b> / {a<sub>1</sub>, a<sub>2</sub>}</td><td>Vector <b>a</b> without the components a<sub>1</sub> and a<sub>2</sub></td></tr><tr><td>|<b>a</b>|</td><td>Euclidean norm of vector <b>a</b></td></tr><tr><td><b>a</b> × <b>b</b></td><td>Cross product of <b>a</b> and <b>b</b></td></tr><tr><td><b>ab</b></td><td>Matrix product of <b>a</b> and <b>b</b></td></tr><tr><td>∇(<b>a</b>)</td><td>Gradient of <b>a</b></td></tr></tbody></table>

the first time. Open-loop control was also proposed for rotating the tip around the target which was validated experimentally. Third, the M-CTR stability was evaluated

numerically and experimentally, as well as active control of the robot's stability. A strategy for identifying adequate magnetic fields is provided and validated.
The paper is organized as follows. The methods and material used to realize the three studies are described in section 2. The results concerning the FTL deployment capabilities, orientability and stability of the M-CTR are described in sections 3, 4 and 5, respectively. According to these results, the interest of the proposed concept of continuum robot is discussed in section 6.

## 1.4. Notations

The abbreviations and notations used in the rest of the paper are gathered in Table 2.

# 2. Methods and materials

FTL deployment capabilities, orientability and stability are assessed using a combination of numerical and experimental approaches. The numerical evaluations are based on a general kineto-static model of M-CTR.

## 2.1. Kineto-static model

The kineto-static model relates the M-CTR configuration to its actuation inputs and external forces. It is elaborated for arbitrary M-CTR designs following the steps introduced in Peyron et al. (2018).

### 2.1.1. Description and assumptions.
M-CTR kineto-statics rely on the interaction of tubes with different stiffness and magnetic properties. We introduce a specific parametrization of the robot which can take into account these different properties. The M-CTR is composed of *n* tubes which are indexed from the outermost to the innermost. When assembled, the tubes form the backbone of the robot, which is then composed of *n* sections with a different number of tubes in interaction. In addition, due to the magnetic elements along the tubes, the sections can be decomposed into several subsections with different magnetic and mechanical properties. This robot segmentation is represented for the design under investigation on Figure 2. The 3 tubes form 3 sections, and section 3 is composed of two subsections due to the presence of the magnet at the tip of the inner tube. The

<figure><img src="image_3.png" /><figcaption>Figure 2. Representation of a M-CTR composed of two pre-curved tubes, an inner straight tube and a magnet at the tip with the associated parametrization.</figcaption></figure>

sections are indexed from the proximal to the distal end. The subsections are indexed from the proximal to the distal end of the corresponding section, the subsection *l* of section *j* being labelled *jl*. The length of subsection *jl* is denoted $\Delta L_{jl}$. The first subsection length $\Delta L_{j1}$ is determined by the relative translation of the tubes, while the others are fixed and depend on the length of the magnetic elements and their location along the tubes.

The actuation inputs are the proximal rotation and translation of each tube, and the magnetic field generated on the magnetic elements. The proximal rotation of tube *i* is denoted $\alpha_i$. The tube proximal translations are parametrized by $\Delta L_{j1}, j = 1...n$. We consider that the magnetic field can be controlled at one location in the workspace, which is generally the case for eMNSs, and that it is not homogeneous. As a result, it varies across the workspace, including along the robot backbone. The magnetic field generated on the magnetic elements situated at an arc-length *s* from the robot base is denoted $\mathbf{B}(s)$. These actuation inputs are related to the robot configuration under the following assumptions:

* The compression and shear deformations are negligible, which is true for thin and slender continuum robots (Rucker et al., 2010; Kratchman et al., 2017).

* The materials composing the tubes have a linear and isotropic behaviour (Rucker et al., 2010).

* The tubes are guided in the transmission section so that they remain straight, which is generally the case for CTRs (Dupont et al., 2010).

* The robot evolves in free-space. Gravity is the only external force applied on the robot backbone.

### 2.1.2. M-CTR configuration.
The M-CTR configuration comprises the robot backbone position, orientation and curvature, and the torsion angle of each tube. These geometrical entities vary along the backbone, which is parametrized by the arc-length *s*. Dependency in *s* is not mentioned in the following for sake of compactness. The backbone position is denoted $\mathbf{p}$ and is expressed in the eMNS reference frame denoted $\mathcal{R}_0(O, \mathbf{x}_0, \mathbf{y}_0, \mathbf{z}_0)$. The robot orientation is parameterized with a Bishop frame $\mathcal{R}_b(\mathbf{p}, \mathbf{x}_b, \mathbf{y}_b, \mathbf{z}_b)$ where $\mathbf{z}_b$ is tangent to the backbone. The rotation between $\mathcal{R}_0$ and $\mathcal{R}_b$ is represented by the unit quaternion vector $\mathbf{Q}$ which satisfies

$$C_{q1} = \mathbf{Q}^T \mathbf{Q} - 1 = 0 \qquad (1)$$

Compared to the Euler angles used in Peyron et al. (2018), this representation is singularity-free, induces quadratic equations more efficiently solved with numerical tools and is classically considered for M-CRs and elastic rods (Tunay, 2013; Edelmann et al., 2017). The backbone curvature is represented by a Darboux vector expressed in $\mathcal{R}_b$ and denoted as $^b\mathbf{u}$. Given the construction of $\mathcal{R}_b$, the coordinate along $\mathbf{z}_b$ of $^b\mathbf{u}$ always equals 0, so that $^b\mathbf{u} = [u_x \ u_y \ 0]$. The torsion angle $\theta_i$ of tube *i* is defined by attaching a frame $\mathcal{R}_i(\mathbf{p}, \mathbf{x}_i, \mathbf{y}_i, \mathbf{z}_i)$ at each cross section of the tube so that $\mathbf{z}_i$ is tangent to the tube centre-line. The pre-
curvature of tube i expressed in this frame is denoted by $\hat{u}_i$. Frame $\mathcal{R}_i$ is then obtained by applying a rotation $(\theta_i, z_b)$ on $\mathcal{R}_b$. In the following, we use left superscript 0, b and i to indicate in which frames vectors are expressed.

Position, orientation and curvature vectors constitute a redundant representation of the robot configuration. As demonstrated in Nikravesh et al. (1985), orientation and curvature are linked by the following relation

$$C_{q2} = {}^b\mathbf{u} - 2L_q Q' = 0 \quad (2)$$

where prime denotes the derivative with respect to s and $L_q$ is a matrix depending on $Q$. Its expression is given in Appendix B.

The backbone position and orientation are also linked. Since shear and extension are neglected, the axis $z_b$ remains tangent to the backbone. The relation between $p$ and $Q$ is then

$$C_p = p' - {}^0\mathcal{R}_b e_3 = 0 \quad (3)$$

where ${}^0\mathcal{R}_b$ is the rotation matrix from $\mathcal{R}_b$ to $\mathcal{R}_0$ and $e_3 = [0 \ 0 \ 1]^T$.

### 2.1.3. Mechanical equilibrium.
The equilibrium conditions relate the robot stiffness to the forces induced by the magnetic field and gravity. We first express them along subsection $jl$, where $n_j$ tubes interact. The robot stiffness depends on the stiffness of each tube present in the subsection, which is expressed as a stiffness matrix $K_i$ for tube i

$${}^i\mathbf{K}_i = \begin{bmatrix} k_{bi} & 0 & 0 \\ 0 & k_{bi} & 0 \\ 0 & 0 & k_{ti} \end{bmatrix} \quad (4)$$

where $k_{bi}$ and $k_{ti}$ are, respectively, the bending and torsional stiffness of tube i. The magnetic field and the gravity exert magnetic forces $f_m$, torques $\tau_m$ and gravity forces $f_g$ which are distributed along the robot's sections. With $m$ as the dipole moment of the magnetic element present on subsection $jl$, $\rho_i$ the linear mass of tube i and g the gravity acceleration, the forces and torques are written

$$\begin{aligned} {}^0\tau_m &= ({}^0\mathcal{R}_b {}^b m) \times {}^0\mathcal{B} \\ {}^0f_m &= F({}^0\mathcal{R}_b {}^b m) G^0 \end{aligned} \quad (5)$$

$${}^0f_g = -\sum_{i=1}^{n_j} \rho_i g \quad (6)$$

where $n_j$ is the number of tubes interacting along subsection $j$, $F({}^0\mathcal{R}_b {}^b m)$ is a 3 × 5 matrix formed with the components of ${}^0\mathcal{R}_b {}^b m$ and $G = [\partial B_x / \partial x \ \partial B_x / \partial y \ \partial B_x / \partial z \ \partial B_y / \partial y \ \partial B_y / \partial z]^T$ contains the 5 independent components of the magnetic field gradient as described in Petruska and Nelson (2015).

The equations describing mechanical equilibrium are obtained by first computing the total potential energy due to the robot deformation, and the forces and torques applied on

the backbone. The Euler-Lagrange formula is then applied to this energy formulation considering the geometrical constraints (1), (2) and (3). The equations consist in a set of differential equations which are written as

$$\left\{ \begin{aligned} & k_{ti} \theta_i'' - \hat{u}_i^T \mathbf{e}_3 + k_{ti}' (\theta_i' - \hat{u}_i^T \mathbf{e}_3) \quad \dots \\ & - {}^b\mathbf{u}^T \frac{\partial^b \mathbf{R}_i}{\partial \theta_i} {}^i\mathbf{K}_i \hat{u}_i + \tau_{mz} = 0 \\ & - \sum_{i=1}^{n_j} k_{bi} {}^b\mathbf{u} + \sum_{i=1}^{n_j} {}^b\mathbf{R}_i {}^i\mathbf{K}_i \hat{u}_i - \lambda_q = 0 \\ & \lambda_p' + {}^0f_g + {}^0f_m = 0 \end{aligned} \right. \quad (7)$$

with $S_q = -2[2\partial L_q^T / \partial s \ \mathbf{Q}]$ and $S_p = \partial^0 \mathcal{R}_b e_3 / \partial \mathbf{Q}$. The first equation describes the equilibrium in torsion of tube i. Since all tubes can rotate freely with respect to each other, it is written for each tube present on the subsection. The torsional magnetic torque $\tau_{mz}$ applies on tube i if a magnet is fixed on this tube, and equals $( {}^b\mathcal{R}_0 {}^0\tau_m )^T e_3$. For the other tubes, $\tau_{mz} = 0$. The other equations describe the equilibrium in bending of the robot backbone. The two matrices $\lambda_q$ and $\lambda_p$ contain the Lagrange multipliers due to the geometrical constraints. Multipliers $\lambda_q$ are related to constraints $C_q1, C_q2$ and $\lambda_p$ to constraint $C_p$. They correspond, respectively, to the internal moments and forces due to the gravity and the magnetic field. The intermediate steps of the model derivation are provided in Appendix B.

The resulting equilibrium equations are valid on the subsection $jl$. Following the approach in Ha et al. (2016) and Peyron et al. (2019a), we extend them to the whole length of the robot by introducing virtual tubes when needed. This means that each tube is considered to cover the entire length of the robot, but with different mechanical and magnetic properties according to the presence of the tubes and magnetic elements. When a tube i is not present along section j, its bending stiffness and linear density are considered as zero and its torsional stiffness as constant. This implies that the stiffness, linear density and dipole moment are piece-wise constant functions of the arc-length s. Equations (7) are then valid for every section of the M-CTR.

### 2.1.4. Boundary conditions.
The differential equation (7) and the geometrical constraints (2) and (3) are subject to proximal (at $s = 0$) and distal (at $s = L$) boundary conditions. First, the backbone position and orientation at $s = 0$, denoted by $p_0$ and $Q_0$, respectively, appear in the following proximal boundary condition

$$\begin{aligned} \mathbf{p}(0) &= p_0 \\ \mathbf{Q}(0) &= Q_0 \end{aligned} \quad (8)$$

The rotation of tube i at its proximal extremity also induces a proximal condition on the tube torsion angle $\theta_i$. In particular, the tube accumulates torsional deformations in the transmission section, resulting in a torsion angle at $s = 0$
that is lower than the actuation angle $\alpha_i$. Since the tubes are considered as straight along the transmission section, their torsional curvature is constant along this section (Dupont et al., 2010). The torsion angle of tube $i$ at $s = 0$ is then written

$$ \theta_i(0) = \alpha_i + \beta_i \theta_i'(0) \quad (9) $$

where the transmission length $\beta_i$ depends on the tube and section length as follows

$$ \beta_i = \sum_{i=1}^{i} \sum_{i}^{n_{ij}} \Delta L_{ij} - L_i \quad (10) $$

Since no force or torque is applied at the robot tip, the tubes are not deformed in torsion at their distal extremity, leading to the distal boundary condition

$$ \theta_i'(L) = \hat{u}_{iz} \quad (11) $$

This assumption also implies a distal condition for the internal moments and forces $\lambda_q$ and $\lambda_p$

$$ \begin{aligned} \lambda_q(L) &= 0_{4 \times 1} \\ \lambda_p(L) &= 0_{3 \times 1} \end{aligned} \quad (12) $$

The equations describing the robot geometry (1–3) and mechanical equilibrium (7), and the boundary conditions (8–12), constitute the kineto-static model of M-CTR. This model allows for computing the 16 + $n$ robot states, denoted by $\mathbf{x} = [\theta_1 \quad \dots \quad \theta_n \quad u_x \quad u_y \quad \lambda_q^T \quad \lambda_q^T \quad \mathbf{Q}^T \quad \mathbf{p}^T]$ as a function of the $2n + 8$ actuation inputs $\mathbf{q} = [\alpha_1 \quad \dots \quad \alpha_n \quad \Delta L_{11} \quad \dots \quad \Delta L_{n1} \quad \mathbf{B}^T \quad \mathbf{G}^T]$.

## 2.2. Numerical framework

CTRs and M-CRs are known to present multiple configurations for a single set of actuation inputs. Loss of stability can also be encountered during robot deployments. Therefore, numerical analysis is performed with the framework for cardinality and stability analysis of continuum robots developed in Peyron et al. (2019a). It has been used successfully to evaluate the FTL deployment, orientability and stability of CTRs (Peyron et al., 2019a, 2019b) and to study the magneto-elastic phenomena of M-CRs (Peyron et al., 2018).

### 2.2.1. Discretization.
The kineto-static model of a M-CTR is a set of non-linear differential equations with proximal and distal boundary conditions. In order to solve it, the model is discretized using finite differences. First, the robot is discretized as a finite number $N$ of nodes along its backbone. In particular, each section is discretized with 20 nodes, which proved to result in a reasonable accuracy in previous work. Each node is indexed from the proximal to the distal extremities of the M-CTR. The states composing the robot configuration are then evaluated at each node. As a result, the robot configuration is represented with a

finite number of states that we gather in a matrix $\mathbf{X} = [\mathbf{x}_1 \quad \dots \quad \mathbf{x}_N]$, where $\mathbf{x}_i$ is the state vector evaluated at node 1. Finally, the differential equations are evaluated at each node by replacing the derivatives by finite differences. A central second order finite difference scheme is used for the second order derivative of $\theta_i$, as proposed in Peyron et al. (2019a). Backward and forward first order finite difference schemes are used for $(\mathbf{p}, \mathbf{Q})$ and $(\lambda_q, \lambda_p)$, respectively, similarly to Peyron et al. (2018), in order to take into account respective proximal and distal boundary conditions.

After discretization, the kineto-static model consists in a set of non-linear equations of the general form

$$ \mathbf{G}(\mathbf{Y}, \mathbf{Q}_N, \mathbf{p}_N, \mathbf{q}) = 0 \quad (13) $$

where $\mathbf{Y} = \mathbf{X}/\{\mathbf{Q}_N, \mathbf{p}_N\}$ and $(\mathbf{Q}_N, \mathbf{p}_N)$ is the pose of the robot distal tip. We then obtain an implicit kineto-static model, which is used to solve forward and inverse kinematics with the same set of equations. Forward kinematics are obtained by considering $\mathbf{q}$ as the input and the robot configuration $\mathbf{X}$ as the output. Inverse kinematics are obtained by considering the tip pose $(\mathbf{Q}_N, \mathbf{p}_N)$ as the input and $\mathbf{Y}$ and $\mathbf{q}$ as outputs.

### 2.2.2. Kinematic and stability analysis.
FTL deployment, orientability and stability analysis were performed on Matlab R2015a (Mathwork Inc.) using the Matcont Toolbox, which implements a continuation method, a step size control algorithm and bifurcation analysis.

FTL deployment capabilities were assessed by simulating deployments using the continuation method. During the simulation, the section lengths were increased sequentially and the forward kinematics were solved at each step with a prediction and correction process. As a result, robot configurations are computed while being robust to non-linear behaviours as demonstrated in the case of CTR in Peyron et al. (2019a).

Orientability was evaluated by using the same continuation process as for the inverse kinematics. The tip position $\mathbf{p}_N$ was then fixed, and the tip orientation $\mathbf{Q}_N$ was varied while the corresponding robot configurations were computed. During the simulation, the orientation variation applied at each step was automatically controlled using a dedicated algorithm, in order to ensure convergence of the model solution. In particular, the variation decreased automatically and ultimately vanished when approaching orientation limits, beyond which the inverse kinematics do not have solutions. As a result, achievable tip orientations were obtained.

Stability changes due to magnetic actuation were evaluated by detecting the appearance of multiple robot configurations for the same actuation inputs. As with CTRs, unstable phenomena observed in M-CRs are characterized by cardinality changes as shown in Peyron et al. (2018). These changes appear as bifurcations during the continuation process. Following Peyron et al. (2018; 2019a), actuation inputs were varied, the forward kinematics were
<table><caption>**Table 3.** Parameters of the m-CTR prototype.</caption><thead><tr><th></th><th colspan="4">Tubes properties</th></tr><tr><th></th><th>Tube 1</th><th>Tube 2</th><th>Tube 3</th><th>Magnet</th></tr></thead><tbody><tr><td>L<sub>i</sub> (m)</td><td>64.00 × 10<sup>-3</sup></td><td>200.50 × 10<sup>-3</sup></td><td>213.50 × 10<sup>-3</sup></td><td>3.00 × 10<sup>-3</sup></td></tr><tr><td>OD (m)</td><td>1.07 × 10<sup>-3</sup></td><td>0.65 × 10<sup>-3</sup></td><td>0.25 × 10<sup>-3</sup></td><td>4.00 × 10<sup>-3</sup></td></tr><tr><td>ID (m)</td><td>0.77 × 10<sup>-3</sup></td><td>0.42 × 10<sup>-3</sup></td><td>0</td><td>2.00 × 10<sup>-3</sup></td></tr><tr><td>u&#770;<sub>tx</sub> (m<sup>-1</sup>)</td><td>14.40</td><td>11.02</td><td>0</td><td>0</td></tr><tr><td>k<sub>bi</sub> (N.m<sup>2</sup>)</td><td>3.89 × 10<sup>-4</sup></td><td>2.00 × 10<sup>-3</sup></td><td>1.20 × 10<sup>-6</sup></td><td>4.88 × 10<sup>-1</sup></td></tr><tr><td>k<sub>fl</sub> (N.m<sup>2</sup>)</td><td>2.88 × 10<sup>-4</sup></td><td>1.50 × 10<sup>-4</sup></td><td>8.94 × 10<sup>-7</sup></td><td>3.84 × 10<sup>-1</sup></td></tr><tr><td>ρ<sub>i</sub> (kg.m<sup>-3</sup>)</td><td>6.45 × 10<sup>3</sup></td><td>6.45 × 10<sup>3</sup></td><td>1.14 × 10<sup>3</sup></td><td>7.46 × 10<sup>3</sup></td></tr><tr><td>m<sub>z</sub> (N.m<sup>-2</sup>.T<sup>-1</sup>)</td><td>0</td><td>0</td><td>0</td><td>1.17 × 10<sup>4</sup></td></tr><tr><th colspan="5">Maximum section lengths</th></tr><tr><td>Section jl</td><td>11</td><td>21</td><td>31</td><td>32</td></tr><tr><td>ΔL<sub>jl</sub> (m)</td><td>50.00 × 10<sup>-3</sup></td><td>50.00 × 10<sup>-3</sup></td><td>13.00 × 10<sup>-3</sup></td><td>3.00 × 10<sup>-3</sup></td></tr><tr><th colspan="5">Initial position and orientation</th></tr><tr><td>p<sub>0</sub> (m)</td><td colspan="4">[ 0.001 -0.0360 0 ]</td></tr><tr><td>Q<sub>0</sub></td><td colspan="4">[ 1 0 0 0 ]</td></tr></tbody></table>

solved and the appearance or vanishing of bifurcations was assessed with bifurcation analysis.

## 2.3. Design of M-CTRs

In order to simplify the analysis, a M-CTR with 3 tubes and controlled with a homogeneous magnetic field was considered. The general structure of the considered robot is depicted in Figure 2. It is composed of two pre-curved tubes and one straight inner tube with a magnetic element fixed at the tip. This results in a concentric assembly of a 2-tube CTR and a M-CR. The two pre-curved tubes were actuated in rotation and translation, and the inner tube in translation only. Since the magnetic field is considered as homogeneous, its gradient is null and does not generate forces on the magnetic element. It can generate torques to bend the robot with 2 degrees of freedom (d.o.f.). As a consequence, the M-CTR has a total of 7 d.o.f. The robot motion due to these d.o.f. is presented in the video of Extension 1 (see Appendix A, Table 6).

Two objectives were considered for the determination of the tube geometrical, mechanical and magnetic properties. The first objective is to allow the investigation of FTL deployment capabilities, orientability and stability. The second objective is to keep consistency with the targeted application, which requires a robot with millimetre diameter and capable to perform in-vivo imaging. The initial shape of the pre-curved tubes was chosen in order to deploy the robot with a perfect FTL behaviour, following the conditions given in Gilbert et al. (2015). The tubes pre-curvature was then optimized in order to observe stability changes under actuation, so that the influence of the magnetic field on the robot stability could be observed. The optimization was conducted while limiting the tube deformation to ε = 0.5% for which the Nitinol is linear (Iasnii and Junga, 2018). The

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Prototype of M-CTR and markers used for the vision-based pose estimation. Colour beads (C1) were used for position measurements across the workspace. An April tag (C3) was used to measure the complete tip pose during evaluation of tip orientability. Another tag (C2) was used to obtain the pose of the actuation unit in the eMNS.</figcaption></figure>

expressions of tube deformation and elastic stability limit in Dupont et al. (2010) were used as objective functions. A Nelder-Mead numerical method, implemented in the fsolve function in Matlab, was used for the optimization. Finally, the inner tube stiffness was chosen so that large deformations of the distal subsection 31 were obtained under magnetic actuation, maximizing tip angular displacement of the M-CTR and therefore its orientability.

The resulting robot is shown in Figure 5. The two pre-curved tubes consist of NiTi tubes with planar and constant pre-curvatures. They were obtained from straight tubes
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Overview of the experimental setup, comprising the eMNS CardioMag eMNS (A1), the M-CTR and its actuation unit (A2), a stereo-vision camera system (A3) and a front camera (A4). Cameras are represented by schematics to improve readability of the figure as they are placed 1.7 m away from the eMNS.</figcaption></figure>

(Euroflex Inc.) which were pre-shaped with a heating process. The tube pre-curvatures after fabrication were measured with a FARO arm and the values are reported in Table 3. The external diameter of the outer tube is 1.07 mm. An optical fibre of 250 µm diameter (SMF28, Corning Inc.) and a permanent ring magnet were assembled on the M-CTR. The use of optical fibres is consistent with in-vivo imaging scenarios (Girerd et al., 2018), and the desired level of flexibility and linearity of material behaviour (Antunes et al., 2008). The magnet is chosen with a conventional axial magnetization, which allows to demonstrate the targeted properties while making the results simple to interpret. A 4 mm diameter magnet is used, which allowed to generate enough magnetic torques considering the specifications of the eMNS. All design parameters are gathered in Table 3.

## 2.4. Experimental setup

The experimental analysis was conducted using the setup presented in Figure 3. It is composed of a M-CTR prototype and its actuation unit (A2), the CardioMag eMNS (Edelmann et al., 2017) for magnetic field generation (A1) and a vision-based measurement system composed of a stereo-vision camera system (A3) and a front camera (A4). The M-CTR actuation unit for tube rotation and translation was specifically designed to be compatible with the CardioMag. The different components are interfaced with a master computer using the Robot Operating System (ROS). They are presented in more detail in the following.

### 2.4.1. Actuation unit.
The actuation unit is inspired from Hendrick et al. (2015b) but is designed to be compatible with the eMNS. The actuation unit is presented in Figure 4. The tubes are mounted on two translation stages (B2)

<figure><img src="image_6.png" /><figcaption>**Figure 5.** CAD representation of the actuation unit. It is composed of two stages (B2) supporting the pre-curved tubes (B3), and of stepper motors (B1) to actuate their rotation and translation. Translation is calibrated using contact switches (B4), and rotation using an optical sensor (B5), an optical marker on gear (B6) and mechanical coding pins (B7,B8).</figcaption></figure>

actuated in rotation and translation using electro-mechanical motors (B1). The entirety of the structure is composed of a non-magnetic material. The structural parts were 3D printed using an ABS polymer. Glass and polymer bearings were used to guide the different axes. The screws used for translation were made out of non-magnetic steel (Inox 1.4301), and the axes for rotation out of aluminium. The tubes were fixed to the translation stages using a 3-finger clamp (B3) machined from a block of ABS. The CTR actuation consists of four stepper motors (17HM15-0904S, Gotronic) located at the unit base, which are not sensitive to the magnetic fields used in this work. The motors are controlled using an Arduino Mega 2560 board.

The actuation unit is equipped with sensors in order to calibrate the tube translation and rotation. Mechanical switches (B4) are used to calibrate the tube translations. In order to calibrate the tube rotations, a mechanical coding system comprising two pins was developed. A first planar pin (B7) is glued to the tubes on a surface plate, so that the pin plane corresponds to the tube plane. A second cylindrical pin (B8) allows for fixing the orientation of the clamp, and thus of the tube, with respect to the gear (B6). An optical sensor (B5) is used to detect a coloured marker placed on one face of the gear (not represented here), which initializes the tube rotation.

### 2.4.2. Electromagnetic navigation system.
The CardioMag is composed of eight fixed coils supplied with independently controlled currents. It is able to generate any magnetic field or field gradient at one point of its workspace, with a maximum magnetic field amplitude of 80 mT (Edelmann et al., 2017). The magnetic field varies from that point across the workspace and was modelled in (Petruska et al., 2017). This model was used to compute the current reference for each coil, according to the desired magnetic field and application point. It was also used to predict the magnetic field along the robot B(s), and its gradient for given currents in the electromagnets.

### 2.4.3. Vision-based measurement system.
A vision-based strategy is adopted to measure the M-CTR pose without disturbing the magnetic field. A stereo tracking system was
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Simulation of the FTL deployment of the M-CTR for two values of the relative orientation $\phi$. The paths followed during the deployments are represented with dashed lines.</figcaption></figure>

implemented using two Basler A602fc cameras (656 × 490 pixels, 15 Hz) fixed on the ceiling (A3). Coloured hollow beads of 6 mm diameter were fixed to the robot (C1). The centres of the beads are detected using standard image processing functions, providing a position accuracy of 1 mm. The positions of the beads are expressed in frame $\mathcal{R}_0$, represented on Figure 5. The transformation between the cameras and $\mathcal{R}_0$ is measured using an April tag (C2), and the AprilTag library (Olson, 2011).

The tip orientation is measured when needed using an April tag fixed on the magnet (C3). Orientation estimation of the marker using the stereo system is not sufficient. A third Basler A602fc camera, focused on the robot tip (A4), provided the missing orientation measurement. The camera was positioned manually in order to maximize the April tag detection during its motion. The camera pose with respect to the robot frame is not known. It is thus used to measure tip position and angular displacements instead of absolute tip position and orientation. The measurement errors were determined to be less than 1 mm and 1°, respectively. They correspond to less than 2% of the robot length and 5% of the tip angular displacement, respectively, which is reasonable for assessing experimentally the orientability and stability of M-CTR.

## 3. FTL deployment

In this section, a strategy for deploying M-CTR in a FTL manner using magnetic actuation is developed, which is analyzed and validated numerically.

### 3.1. Perfect FTL deployment of a M-CTR

When a path is composed of sections with constant curvatures, FTL deployment with a continuum robot is possible. The robot must be composed of sections with constant curvature and variable length. CTRs can fulfil these two requirements by considering tubes with constant pre-curvature, planar or helical, and by placing them so that their curvatures are aligned or opposed (Gilbert et al., 2015). However, since the tube curvature is fixed, the robot can

follow a limited number of paths following the tube alignments.

Magnetic actuation can be used to overcome this limitation with the considered design of M-CTR. Aligning or opposing the pre-curved tubes implies that the robot is contained in a plane. Without loss of generality, we here choose to work in the $(y_0, z_0)$ plane as shown on Figure 6(a). Applying a magnetic field in this plane generates a constant bending moment along the backbone (Chautems et al., 2018). As a result, neglecting the effect of gravity, the curvature of the robot backbone stays constant along the sections, but varies according to the magnetic field. Considering that the robot is in the $(y_0, z_0)$ plane and the permanent magnet has an axial magnetization, the expression of the backbone curvature in equation (7) reduces to

$$ 
\begin{align} 
{}^b\mathbf{u} &= [u_x \quad 0 \quad 0]^T \nonumber \\ 
u_x &= \begin{cases} \frac{k_{b2}\hat{u}_{1x} + k_{b2}\hat{u}_{1x}}{k_{b1} + k_{b2}} + \frac{\lambda_{qx}}{k_{b1} + k_{b2}} & s \in [0, \Delta L_{11}] \\ \hat{u}_{2x} + \frac{\lambda_{qx}}{k_{b2}} & s \in [\Delta L_{11}, \Delta L_{11} + \Delta L_{21}] \\ \frac{\lambda_{qx}}{k_{b3}} & s \in [\Delta L_{11} + \Delta L_{21}, \Delta L_{11} + \Delta L_{21} + \Delta L_{31}] \end{cases} \nonumber \\ 
\lambda_{qx} &= \tau_{mx} \qquad \tau_{mx} = |m_z\mathbf{e}_3 \times ({}^b\mathbf{R}_0(\mathbf{Q}_N)^0\mathbf{B})|^ \tag{14} 
\end{align} 
$$

Controlling the magnetic field so that the magnetic torque $\tau_{mx}$ stays constant during the deployment leads to a constant value of ${}^b\mathbf{u}$, and the condition for perfect FTL behaviour is fulfilled. The torque depends on the magnetic field magnitude and orientation with respect to the magnet. We chose to keep the magnetic field magnitude constant during the deployment. Therefore, the magnetic field direction is computed so that its relative orientation with the magnet $\phi$ (see Figure 6(a)) stays constant with the relation

$$ {}^0\mathbf{B} = |\mathbf{B}|\mathbf{R}(\mathbf{Q}_x(\phi))\mathbf{Q}_N)e_z \tag{15} $$

where $\mathbf{Q}_x(\phi)$ is the quaternion representation of the rotation about $\mathbf{x}_0$ of angle $\phi$, and $\mathbf{R}(\mathbf{Q})$ is the rotation matrix
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Simulation of the influence of $\phi$ on the path followed by the M-CTR during FTL deployment.</figcaption></figure>

associated to the quaternion $Q$. The two quaternion vectors are multiplied using the Hamilton product operator.

The effectiveness of this strategy was analyzed by simulating the deployment for several values of $\phi$. The pre-curved tubes were placed so that their curvatures were aligned, that is, $\alpha_1 = \alpha_2 = 0$. The magnetic field was assumed to be controlled at the magnet centre and its magnitude was fixed at $|\mathbf{B}| = 80$ mT. The gravity was set to 0. The section lengths were increased sequentially from the proximal to the distal section until their maximum value was reached (see Table 3). The deployment consists of 3 steps during which $\Delta L_{11}$, $\Delta L_{21}$ and $\Delta L_{31}$ are varied, respectively. The section length variations and the resulting robot shape are computed using the continuation method. The error of the FTL deployment is evaluated as the mean distance between the desired backbone node position along the path and the actual node position for each value of section length.

The results of the simulations are presented on Figure 6, where robot configurations during the deployment and the followed paths are represented. A video of the complete deployment is provided in Extension 2 (see Appendix A, Table 6). The M-CTR exhibits a perfect FTL behaviour as expected. It is able to follow the desired paths with a mean FTL error during the deployments below 0.54 mm, which represents 0.47% of the robot total length. Such an error is negligible with respect to the values reported in the literature (Garriga-Casanovas and Rodriguez y Baena, 2018). Moreover, the M-CTR distal subsection 31 can be deployed along a curved path although the inner tube is initially straight, and different curvatures of this path can be obtained for different values of $\phi$. This is an improvement over the deployment capabilities of CTRs, for which this subsection would remain straight.

## 3.2. Influence of the magnetic field orientation

The range of paths that can be followed with perfect FTL behaviour depends on the achievable range of $\phi$. Following

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Evolution of $\phi_M$ according to the inner tube bending stiffness.</figcaption></figure>

equation (14), the magnetic bending moment is maximum when the magnetic field is orthogonal to the magnet axis, that is, when $|\phi| = \pi/2$. In practice however, the maximum value of $\phi$ is limited by the magneto-elastic instabilities inherent to the behaviour of M-CR (Peyron et al., 2018; Tunay, 2004). When the magnetic field is oriented in the robot plane, there is a critical magnetic field orientation above which the robot becomes unstable and can have multiple configurations for the same inputs. This critical orientation then determines the maximum value of $\phi$, denoted $\phi_M$ that can be achieved while conserving stability of the M-CTR.

The value of $\phi_M$ is computed by considering the worst case scenario in terms of magneto-elastic stability. Tunay (2004) proved that the longer the robot, the higher the risks of becoming unstable. The M-CTR is then considered fully deployed for the evaluation of $\phi_M$. We use then the same process as described in Peyron et al. (2018). The angle $\phi$ is varied with the continuation process until a stability change is detected. The value of $\phi_M$ corresponds to the value of $\phi$ at the stability change, which equals $\phi_M = 6.31^\circ$ with the parameters of the M-CTR prototype. The corresponding range of paths is computed with equations (2), (3) and (14) and is presented in Figure 7. The paths exhibit large variations in curvature for the distal section and small variations for the proximal ones. This is due to the low bending stiffness of the inner tube with respect to the NiTi tubes and the magnetic field magnitude. For low values of $\phi$, the bending moment induced by the magnetic field is too small to affect the shape of the NiTi tubes, but is high enough to deform the inner tube significantly.

As a consequence, the curvature of the robot distal section can be varied continuously by changing the relative orientation between the magnetic field and the magnet. Moreover, significant variations of curvature were observed. The minimum radius of curvature was 5.4 mm, which is interesting considering the scale and the geometry of environments in which such type of robot could be deployed. This implies that tip angular displacements of $\pm 160^\circ$ can be
achieved, which is comparable to the best displacements presented in the literature according to Table 1.

### 3.3. Influence of the inner tube stiffness

The inner tube stiffness is an important quantity for the deployment capabilities of the M-CTR. A low inner tube stiffness with respect to the NiTi tube sets the global shape of FTL paths to have constant curvatures along subsections 11 and 21, and a variable curvature for subsection 31. The stiffness also has a significant impact on the magneto-elastic phenomena and then on the value of $\phi_M$. The evolution of $\phi_M$ according to $k_{b3}$ is presented in Figure 8. The values of stiffness were chosen in the range $[1 \times 10^{-8}, 1 \times 10^{-4}]$ N.m$^2$, with a logarithmic distribution, to consider materials such as Nylon, NiTi, and stainless steel. For the three highest values of $k_{b3}$, the robot does not experience any stability changes under magnetic field rotations. Therefore, we report on the graph the value of $\phi$ which produces the maximum magnetic torque, which is equal to 90°, as explained previously. As the stiffness decreases, $\phi_M$ decreases rapidly and stays below 3° for $k_{b3} < 5 \times 10^{-7}$ N.m$^2$.

This variation of $\phi_M$ implies a variation of the bending moment applied on the robot, which impacts the deformation of the pre-curved tubes due to magnetic actuation.

The paths which can be followed in a FTL manner are computed for each value of stiffness and for the extreme values of $\phi$ denoted $\phi_M$ and $-\phi_M$. A subset of these paths is presented in Figure 9. Large values of $k_{b3}$ imply large values of $\phi_M$ and consequently large values of the magnetic bending moment. As a result, the NiTi tubes are deformed significantly by the magnetic field. Such values also imply that variations of subsection 11 curvature are smaller. Changing the value of $\phi$ globally modifies the shape of the path, while the curvature experiences small variations from one section to

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Evolution of the paths that can be followed according to the inner tube stiffness. For each value of stiffness, only the paths obtained for $\phi = \phi_M$ and $\phi = -\phi_M$ are represented.</figcaption></figure>

another. As the stiffness decreases, so does $\phi_M$ and the deformation of the NiTi tubes. For $k_{b3} \le 5.62 \times 10^{-7}$ N.m$^2$, their deformation is negligible and the curvature of section 2 is high, leading to paths with similar shapes than the paths obtained with the M-CTR prototype. Note that in practice this curvature is also limited by self-collisions between the robot tip and its backbone, which are not accounted for here.

Changing the inner tube stiffness allows for deployment of the robot with a perfect FTL behaviour along paths with different shapes. In particular, increasing $k_{b3}$ allows for the robot to withstand a higher magnetic torque while staying stable, high enough to deform the NiTi tubes. This results in a global modification of the path, with significant variations of the path position according to the magnetic field orientation.

## 4. Orientability

In this section, the orientability of M-CTR is analyzed numerically and demonstrated experimentally in a tip orientation scenario using open-loop control.

### 4.1. Orientability definition and investigation strategy

In the literature, the orientability of a robot at one point of its workspace $p_d$ is defined as the number of orientations achievable by the robot tip while maintaining its position at $p_d$ (Wu et al., 2017). In most related work, continuum robots are considered to have 6 d.o.f., which allows for full control of their tip pose (Chikhaoui et al., 2016; Li et al., 2017; Wu et al., 2017). In these works, the evolution of orientability across the robot workspace and the impact of actuation redundancy were studied. Orientability evaluation is based on solving simple kinematic models of continuum robot based on the constant curvature assumption (Chikhaoui et al., 2016; Li et al., 2017; Wu et al., 2017). As a result, phenomena such as torsion of the backbone and deformation due to external forces are neglected. The tube torsion has been considered in a preliminary work in Peyron et al. (2019b) for CTRs. Orientability evaluation is performed here for the first time accounting for the effect of external forces produced by the magnetic field and the gravity.

The study of M-CTR orientability is performed in three phases. First, an orientability map across the robot workspace is generated. A simplified forward static model is used that neglects torsion deformations and gravity effects. Second, orientability is evaluated numerically at several locations in the workspace by solving the inverse kinematic problem in (13) with our numerical framework, accounting for the first time for torsion deformations and gravity effects. Third, these results are used to achieve tip orientation control of the M-CTR prototype and to demonstrate its orientability experimentally.

In order to simplify the study of M-CTR orientability, a minimal set of actuation inputs among the 7 available variables was used. Since the tubes can be rotated at their
base, the control of the axial tip orientation is trivial. Consequently, we are not interested in this rotation and only 5 d.o.f. are considered. Two actuation inputs are considered as constant during the tip rotation, which are $\alpha_3$ and $\Delta L_{31}$. In particular, subsection 31 is completely deployed in order to maximize the angular displacement under magnetic actuation and $\alpha_3$ is set to 0. The minimal set of actuation inputs is then [$\alpha_2$, $\Delta L_{11}$, $\Delta L_{21}$, $\mathbf{B}$].

## 4.2. Orientability evaluation

As an initial step, the orientability map is generated following the method proposed in Wu et al. (2017). The forward kinematics of the M-CTR are derived assuming constant curvature, and neglecting gravity and deformations due to tube torsion. We also assume that the magnetic field only deforms subsection 31. Subsection 31 has a constant curvature $u_{mag}$ and it deforms in a plane of osculating angle $\phi_{mag}$. In that case, the curvature along the backbone is

$$ 
\begin{align*}
{}^b\mathbf{u} &= [u_x \quad u_y \quad 0]^T \\
u_x &= \begin{cases} \frac{k_{b1}\hat{u}_{1x} \cos(\alpha_1) + k_{b2}\hat{u}_{2x} \cos(\alpha_2)}{k_{b1} + k_{b2}}, & j=4 \\ \hat{u}_{2x} \cos(\alpha_2), & j=3 \\ u_{mag} \cos(\phi_{mag}), & j=2 \\ 0, & j=1 \end{cases} \\
u_y &= \begin{cases} \frac{k_{b1}\hat{u}_{1x} \sin(\alpha_1) + k_{b2}\hat{u}_{2x} \sin(\alpha_2)}{k_{b1} + k_{b2}}, & j=4 \\ \hat{u}_{2x} \sin(\alpha_2), & j=3 \\ u_{mag} \sin(\phi_{mag}), & j=2 \\ 0, & j=1 \end{cases} \tag{16}
\end{align*}
$$

The robot configuration is obtained by integrating (2, 3) with this expression of curvature.

The actuation space was discretized into $40 \times 10^6$ sets of actuation inputs using the Monte Carlo method, and accounting for the ranges of variation summarized in Table 4. The model was solved for each set of actuation inputs, leading to a set of robot tip positions and orientations. The workspace was then discretized into $4 \times 4 \times 4$ mm cubic regions, of which the centres correspond to the desired position $\mathbf{p}_d$. The robot configurations which tip positions lie in a same cubic region were then collected, producing a set of tip orientations.

The set of tip orientations was finally used to evaluate the orientability at the corresponding $\mathbf{p}_d$. The tip orientations are represented on a service sphere of area $A_S$ located at the tip position. This sphere is discretized into a finite number of patches, representing the different possible tip orientations. The area $A_R$ of the service region is determined from the number of patches intersected by the backbone tangent at the tip. The orientability $D(\mathbf{p}_d)$ is then computed as

$$ D(\mathbf{p}_d) = A_R/A_S \tag{17} $$

The best robot orientability is obtained when $D(\mathbf{p}_d)$ gets close to 1.

The M-CTR orientability was evaluated for each cubic region, leading to the heat map presented on Figure 10. The workspace of the M-CTR is symmetric with respect to the $(\mathbf{y}_0, \mathbf{z}_0)$ plane, where the orientability is maximum. Therefore, we represent orientability in the $(\mathbf{y}_0, \mathbf{z}_0)$ plane only. The maximum orientability is equal to 0.98. It indicates that almost any tip orientation can be achieved with the robot tip. Three locations, labelled 1, 2 and 3, were chosen to evaluate M-CTR tip orientability. Location 1 was chosen at the centre of the workspace to facilitate experimental evaluation while maximizing the orientability. Location 2 and 3 were obtained by setting $(\alpha_2, \alpha_3) = (0, 0)$ to limit stability issues and by moving progressively towards the workspace boundaries. The service spheres at the 3

<table><caption>**Table 4.** Ranges of actuation inputs considered during the actuation space discretization.</caption><thead><tr><th>a<sub>2</sub></th><th>&Delta;L<sub>21</sub></th><th>&Delta;L<sub>11</sub></th><th>u<sub>2</sub></th><th>&phi;<sub>2</sub></th></tr></thead><tbody><tr><td>[0, 2&pi;]</td><td>[0, 0.05]</td><td>[0, 0.05]</td><td>[-u<sub>2M</sub>, u<sub>2M</sub>]</td><td>[-&pi;, &pi;]</td></tr><tr><td colspan="4">u<sub>2M</sub> = m<sub>z</sub>|<b>B</b>|/k<sub>b2</sub></td><td></td></tr></tbody></table>

Subsection lengths and curvature are, respectively, in m and m⁻¹, respectively.

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Orientability of the m-CTR across its workspace. The three locations considered for orientability assessment are represented with red squares. Configurations leading to these locations according to the constant curvature model are also represented.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Service spheres obtained with the constant curvature model for three locations in the workspace. The service regions are depicted in green.</figcaption></figure>

locations are represented on Figure 11, and the orientability equals 0.93, 0.83 and 0.38, respectively.

Beyond the generation of the orientability map, the continuation method is used with the developed kinematic static model to evaluate M-CTR orientability, by considering 1D rotations of the tip starting from an initial configuration of the robot. The initial configurations of each location are selected from the set computed previously, so that the robot tip is in the corresponding cubic region and the robot is in the $(y_0, z_0)$ plane. They are represented on Figure 10. We consider then the tip to be oriented around $x_0$, $y_0$ and $z_0$. The corresponding rotations are represented with quaternions $Q_x$, $Q_y$ and $Q_z$, respectively, so that

$$ 
\begin{aligned}
Q_x &= [\cos(a_x/2) \quad \sin(a_x/2) \quad 0 \quad 0]^T \\
Q_y &= [\cos(a_y/2) \quad 0 \quad \sin(a_y/2) \quad 0]^T \\
Q_z &= [\cos(a_z/2) \quad 0 \quad 0 \quad \sin(a_z/2)]^T
\end{aligned} 
\quad (18) $$

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Initial and extreme M-CTR configurations obtained during tip orientation. (1,2,3): orientation around $x_0$. (1,4,5): orientation around $y_0$. (1,6,7): orientation around $z_0$.</figcaption></figure>

<table><caption>**Table 5.** Angular displacements obtained during M-CTR tip rotation at the selected locations in the workspace.</caption><thead><tr><th rowspan="2">Location</th><th colspan="3">Angular displacement</th></tr><tr><th>x<sub>0</sub></th><th>y<sub>0</sub></th><th>z<sub>0</sub></th></tr></thead><tbody><tr><td>1</td><td>28.4°</td><td>177.3°</td><td>188.3°</td></tr><tr><td>2</td><td>57.5°</td><td>347.2°</td><td>360.0°</td></tr><tr><td>3</td><td>111.0°</td><td>30.0°</td><td>83.6°</td></tr></tbody></table>

where $a_x$, $a_y$ and $a_z$ are the rotation angles. Moreover, the target around which the tip is rotated is considered at a distance $|d| = 3$ mm from the robot tip, which corresponds to the thickness of the April tag used to measure the orientation. To perform the 1D rotation, the inverse kinematic model is solved from (13) by considering the desired tip orientation $Q_N$ and position $p_N$

$$ 
\begin{aligned}
Q_N &= Q_{xyz}Q_{N,0}Q_0 \\
p_N &= p_d - |d|^0R_B(Q_N)e_z
\end{aligned} 
\quad (19) $$

where $Q_{x|z}$ is to be replaced with the desired rotation quaternion in (18), $Q_{N,0}$ is the initial tip orientation and $Q_0 = [\cos(\theta_{1,N}/2) \quad 0 \quad 0 \quad \sin(\theta_{1,N}/2)]^T$ represents the tip axial rotation of angle $\theta_{1,N}$. The initial robot configuration computed with the constant curvature model is used as an initial guess for the inverse kinematics solution. The numerical framework is finally used to vary the rotation angles ($a_x$, $a_y$, $a_z$) and to compute the corresponding robot configurations.

Robot configurations obtained at location 1 are presented in Figure 12. The initial M-CTR configuration is labelled 1 on the figure. Configurations (2,3), (4,5) and (6,7) are obtained when rotating the tip around $x_0$, $y_0$ and $z_0$, respectively. We obtain robot configurations that lead to the same tip position with different tip orientations. Since the initial tip tangent is almost aligned with $y_0$, rotating about
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Block diagram of the open-loop control used for the experimental demonstration of M-CTR orientability.</figcaption></figure>

this axis leads only to small backbone motions. The angular displacements obtained for the 3 locations are presented in Table 5. The results show significant differences of angular displacement between the axes of $R_0$. The M-CTR tip orientability is not isotropic in the workspace. Interestingly, relatively small angular displacement is obtained about $\mathbf{x}_0$ at the location 1, where the orientability is expected to be maximum. This is probably due to the impact of gravity and tube torsion, which could not be anticipated with simpler constant curvature forward kinematics.

The major part of the obtained values is in the range [90°, 360°], which is higher than the angular displacements reported in the literature (Goldman et al., 2013; Simaan et al., 2009). According to these simulations, the M-CTR can be used to rotate a tool fixed at a given position about a specific axis with large amplitudes.

### 4.3. Orientation control

A demonstration of the M-CTR tip rotation capabilities was performed experimentally at location 1 in the workspace. We focus especially on tip rotation around $\mathbf{x}_0$ and $\mathbf{z}_0$, since variations of the robot shape are the most important. We perform the tip rotations by following a pre-computed path in actuation space with quasi-static open-loop control. The block diagram of the control strategy is presented on Figure 13. The path was first planned in the task space, that is, for desired tip position $\mathbf{p}_d$ and tip orientation angle $a_x$, $a_z$. Starting from the initial configuration, $a_z$ (resp. $a_x$) was increased with a constant step size towards its maximum value, moving from configuration 1 to 6 (resp. 1 to 3) on Figure 12. The maximum value of $a_x$ is chosen as 25°, in order to leave some margin with the maximum angular displacement evaluated in simulation. The maximum value of $a_z$ is limited to 90° to ensure the April tag detection. Then angle $a_z$ (resp. $a_x$) was then decreased towards its minimum value (configuration 6 to 7, minimum values of -90°, resp. 0°) and increased again to return to the initial tip angle (configuration 7 to 1). The step size was computed so that the entire rotation cycle is achieved in 40 steps. The path in task space was then expressed in the actuation space using the M-CTR inverse kinematic model. The actuation inputs computed during the simulation are interpolated according to the desired path. The desired magnetic field $\mathbf{B}_d$ was then generated at $\mathbf{p}_d$ with the CardioMag, and the actuation inputs $\alpha_2$, $\alpha_3$, $\Delta L_{11}$ and $\Delta L_{21}$ were sent to the M-CTR actuation unit. The steps were performed with a period of 4s, so that the robot reaches its steady state at the end of each

step. Each cycle was repeated 5 times. Videos of the robot motion during the rotation cycles are presented in Extension 3 (see Appendix, Table 6). In the videos, the visual markers are removed and the open-loop trajectory is played continuously instead of step by step. During the rotation cycle, the robot tip pose was estimated using the camera A5 (see Figure 3). Details concerning the computation of the tip pose are provided in appendix C.

The mean values of the measured tip orientations over the 5 cycles are projected on service spheres shown in Figure 14. We observe that the points are globally contained in a plane, of which the normal vector corresponds to the tip rotation axis. It indicates that the tip is rotated around a fixed axis, as expected. The experimental angular displacement and the rotation errors were computed following the equations provided in appendix. We obtained angular displacements of $\Delta a = 22.24 \pm 0.01^\circ$ and $\Delta a = 136.36 \pm 1.2^\circ$ for the rotations around $\mathbf{x}_0$ and $\mathbf{z}_0$, respectively. They represent 11.04% and 24.24% of the desired angular displacement around the two axes. Rotation errors for the two axes were $\zeta = 0.18^\circ$ and $\zeta = 7.76^\circ$, respectively, which represent 0.81% and 4.97% of the angular displacement. Higher angular displacement and rotation errors for the rotation around $\mathbf{z}_0$ can be attributed to friction which is higher since the middle tube is rotated.

The evolution of the tip displacement error $\epsilon$ is presented as a function of the desired tip rotation angle in Figure 15. Each point corresponds to the mean error value obtained for the corresponding tip angle over the 5 cycles. The mean deviation from these mean values is 0.22 mm. The mean position error for tip rotations about $\mathbf{x}_0$ and $\mathbf{z}_0$ are 0.8 mm and 6.3 mm, respectively, which represents 0.93% and 7.33% of the robot length. Considering that the standard errors of kineto-static models of continuum robots are up to 5% of the robot length, the errors obtained here seem acceptable for a first prototype.

As a result, the proposed open-loop control can be used to perform M-CTR tip rotation about the axes of $R_0$ while staying at the same position in the workspace. As predicted by the model, angular displacements surpassed 90° for the rotation around $\mathbf{z}_0$, showing the key advantage of M-CTRs in terms of orientability.

## 5. Stability

In this section, the stability modulation capabilities of M-CTR are studied by developing a modulation strategy and validating it both numerically and experimentally.
<figure><img src="image_15.png" /><figcaption>Figure 14. Experimental tip orientation around $x_0$ (left) and $z_0$ (right).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15. Experimental tip position errors during tip orientation around $x_0$ (top) and $z_0$ (bottom).</figcaption></figure>

## 5.1. Magnetic field for stability modulation

Unstable phenomena of CTRs are similar to beam buckling, as demonstrated in Gilbert et al. (2016). When the tubes are placed in opposition, increasing the tube interaction length is analogous to increasing compression forces on a straight beam. This increase eventually leads to a critical interaction length (resp. a critical compression force) for which cardinality and stability changes are observed. The configuration where the tubes are placed in opposition is particularly important because it is the worst case scenario from a stability point of view. When tubes are rotated with respect to one another and their interaction length is gradually increased, this is the first configuration to become unstable, as shown in Hendrick et al. (2015a).

From the definition of stability, the M-CTR is stable when the tubes are placed in opposition, if it returns to the same configuration when a perturbation is applied. The robot can be ensured to keep this configuration by applying forces and

torques on the backbone, which has been shown for a buckled beam with tendons (Berlin and Sussman, 1994), shape memory alloy (Baz et al., 1992) and piezoelectric actuators (Schaeffner et al., 2016). In the case of the M-CTR, the magnetic field is used to generate these restoring torques. For the analysis, the inner tube is considered as fully retracted, so that the magnetic torques apply directly on the pre-curved tubes. For the analysis, the inner tube is considered as fully retracted, i.e. its tip aligns with the tip of the middle tube and $\Delta L_{31} = 0$mm. As a consequence, the magnetic torques apply directly on the pre-curved tubes. As the NiTi tubes are significantly stiffer than the inner tube, the robot does not experience the magneto-elastic instabilities mentioned in section 3. The robot is considered as contained in the ($y_0, z_0$) plane when the pre-curved tubes are in opposition. Our strategy to generate the restoring torques is to consider a magnetic field in this plane as well, of which the relative angle with the magnet axis is $\phi$, as introduced in section 3, Figure 6. The buckling of the pre-curved tubes implies out-of-plane motion of the robot and of the magnet. This motion is represented with a rotation of angle $\gamma$ about $y_b$. Assuming that the magnetic field is homogeneous and following equation (5), the magnetic torque applied on the magnet is written in $R_b$

$$ \tau_m = [-m_z|\mathbf{B}|c_\gamma s_\phi - m_z|\mathbf{B}|s_\gamma c_\phi \quad m_z|\mathbf{B}|s_\gamma s_\phi]^{\mathrm{T}} \quad (20) $$

where $c_\gamma$ and $s_\gamma$ (resp. $c_\phi$ and $s_\phi$) denote the cosine and sine of angle $\gamma$ (resp. $\phi$). We are especially interested in the component about $y_b$, which induces out-of-plane motion. When $\phi \in [-\pi/2, \pi/2]$, $\tau_{nr}y_b$ is of opposite sign to $\gamma$. The magnetic torque induces tip motion opposed to the out-of-plane perturbation and thus restores the initial configuration of the M-CTR. On the contrary, for $\phi \in [\pi/2, 3\pi/2]$, $\tau_{nr}y_b$ and $\gamma$ have the same sign. The magnetic torque favours out-of-plane motion. Choosing a magnetic field aligned with the magnet when the tubes are in opposition ($\phi = 0$) leads therefore to the largest restoring torque and improves the M-CTR's stability. Choosing a magnetic field opposed to the magnet ($\phi = \pi$) deteriorates the robot's stability. These two magnetic fields are denoted $B_1$ and $B_2$ in the following.

## 5.2. Stability analysis

The impact of the fields $B_1$ and $B_2$ was determined by performing a stability analysis with the numerical method. In the literature, the stability of CTRs was studied by detecting stability and cardinality changes during tube rotation, or during tube deployment when they are assembled in opposition. Only the first scenario has been used to study stability experimentally. The second scenario consists in placing the robot in its buckling configuration, which was proven to be sensitive to small defects in experimental setups for M-CR (Singh et al., 2013). Therefore, we consider the first scenario.

First, the kineto-static model of M-CTR is solved with the tubes being placed in the buckling configuration to
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Bifurcation diagrams of the M-CTR during inner tube rotation. The LP bifurcation at $α_2 = 180°$ belongs to the graph obtained with no magnetic field.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 17.** Tip coordinate along $x_0$ during inner tube rotation and for the different magnetic fields. Arrows indicate the evolution of the tip pose during clockwise and counter-clockwise rotations.</figcaption></figure>

obtain the magnet orientation. The tubes are placed in opposition (($α_2, α_3$) = (180, 0)°) and deployed until the critical interaction length is reached. In that way, the pre-curved tubes are at a stability limit, and the impact of the magnetic field on this limit can be directly observed. This critical length is computed with the criterion in Hendrick et al. (2015a) and equals 43.2 mm. The tube 2 is also fully deployed. The magnetic fields **B**₁ and **B**₂ are computed considering the magnet orientation and a magnitude $|**B**| = 80$ mT with the relation (15). They are then applied and considered as constant during the tube rotation. Starting

<figure><img src="image_19.png" /><figcaption>**Figure 18.** Vectors and points used to define the angular displacement during tip rotation and the rotation error.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 19.** Points on the service sphere during tip rotation around $z_0$ after projection.</figcaption></figure>

from the configuration where ($α_2, α_3$) = (0, 0)°, the tube rotation $α_2$ is varied with continuation while $α_3$ is held constant.

The resulting diagrams for the case without magnetic fields and with the fields **B**₁ and **B**₂ are presented in Figure 16 in terms of the tip coordinate along $x_0$. The case without magnetic fields reveals a S-shaped diagram with a vertical tangent at $α_2 = 180°$. As expected, the robot is at the stability limit for this tube rotation. In presence of **B**₁, the diagram is similar but with a lower tangent steep angle at $α_2 = 180°$, indicating that the robot is stable during the full tube rotation. In presence of **B**₂, two limit points (LP) appear during the tube rotation, the configurations between these points being unstable.

We can draw two conclusions from these simulations. First, the magnetic field has a significant impact on the robot stability, even if it does not deform the NiTi tubes significantly. This is due to the fact that in the buckling configuration, the robot behaviour is sensitive to small perturbations, produced here by the magnetic torque. Second, the strategy for selecting magnetic fields is efficient. The fields **B**₁ and **B**₂ stabilize and destabilize the robot, respectively.

## 5.3. Experimental analysis

We reproduced the scenario with the experimental setup. The tubes were initially aligned and deployed until the critical interaction length was reached. The critical interaction length of the prototype differs from the predicted
value. It was evaluated using a trial and error process, by slowly translating the tubes and performing full rotation of tube 2, until unstable phenomena were observed. We determined a value of 15 mm. The significant difference can be attributed to the friction between the tubes, which is not taken into account in the criterion of Hendrick et al. (2015a). The magnetic fields $B_1$ and $B_2$ were computed with this value of critical interaction length. They were generated with the CardioMag at the location corresponding to the tip position when the tubes are in opposition. The tube 2 was then rotated from 0° to 360° and in the inverse direction with 80 equally-spaced angular steps. Since the magnetic field is constant during the tube rotation and the stepper motors have a higher bandwidth than the CardioMag, the robot converges faster towards its steady state. Therefore, the steps were sent to the actuation unit with a period of 2 s. The rotation cycle was repeated 5 times. During the tube rotation, the tip position was measured using the stereoscopic vision system.

The resulting experimental diagrams are represented in Figure 17, and videos of the robot motion during tube rotation are provided in Extension 4 (see Appendix A, Table 6). In the videos, the tube rotation was performed continuously. Each point corresponds to the mean value of the tip coordinate over the 5 cycles. The mean deviation is 0.35 mm. When no magnetic field is applied, the robot becomes unstable near $α_2 = 180°$. When applying $B_1$, the robot stays stable during the full tube rotation in the two directions. Hysteresis is visible on the corresponding diagram, which is due to friction between the tubes. When applying $B_2$, the robot snaps with higher amplitudes and dynamics.

In conclusion, the numerical results were confirmed experimentally. The magnetic field has a significant impact on the robot stability, and it can be used to provoke or to prevent snapping. The strategy for computing the magnetic fields also showed its benefits experimentally, despite the presence of unmodelled phenomena such as friction.

## 6. Conclusion and perspectives

In this paper, M-CTRs are introduced and analyzed as candidates for performing MIS at millimetre-scales and in difficult-to-access areas. The interest of the M-CTR is assessed on a representative design of robot through numerical and experimental studies. These studies required the establishment of a generic kineto-static model, the application of a numerical framework and the implementation of a specific experimental setup. They led to a number of results concerning three important performance criteria in the targeted medical application, namely the FTL deployment capabilities, the orientability and the stability.

The robot can be deployed along planar paths with a perfect FTL behaviour by adjusting the magnetic field orientation during the deployment. A continuous set of paths with different curvatures can be followed by using magnetic actuation. In particular, the M-CTR can achieve

distal radii of curvature down to 5.4 mm and angular displacements of ±160° all with a backbone which external diameter does not exceed 1 mm. This angular displacement is on par with the best performing continuum robots at that size. In addition, the global shape of the paths can be modulated by changing the inner tube stiffness.

The M-CTR can be used to orient a tool at its tip around a target. Numerical studies show that the tip orientability depends on the tip position in the workspace, and the axis around which the tip is rotated. For three locations in the workspace, angular displacements in the range [90°, 360°] are obtained which are comparable to and even greater than the angular variations reported in the literature. M-CTR orientability was also demonstrated experimentally by implementing open-loop control of the tip pose. Angular displacements of 136.4° were obtained, which exceed the state of the art. Moreover, open-loop control is promising as a control strategy for M-CTRs, given its simplicity.

Magnetic fields can be used to alter the stability of the M-CTR. In particular, the robot can be stabilized or destabilized by proper selection of the magnetic fields, which was demonstrated numerically and experimentally.

This work opens several perspectives that we intend to develop in the future. The M-CTR will be designed for assessment in medical applications such as olfactory cells inspection and middle-ear surgery. We will in particular focus on the optimization of the number of tubes, the tube properties and the development of accurate control for teleoperated or autonomous procedures. The robot will be composed of a navigation section, with tubes dedicated to the FTL deployment, and an exploration section dedicated to tip pose control as proposed for CTR in Mitros et al. (2018). The number, pre-curvatures and stiffness of the tubes in the navigation section will be optimized to conform with the specific geometries of nasal cavities and ear canals. The exploration section will be designed to fit the dexterity requirements. Also, the dimensions of the magnet will be optimized in order to fit the required 1~mm external diameter. This will be a challenge as the magnet volume should be large enough to produce a sufficient magnetic torque on the concentric tubes, which in turn depends on the tube properties and the application requirements. Potential solutions to this problem are to include the constraint on the magnet size in the tube optimization, and to consider other magnetic elements such as tubes composed of flexible ferro-magnetic material (kim et al., 2019). Concerning control, the proposed open-loop control shows promising results but its accuracy could still be improved. Closed-loop control will be developed by considering sensor integration and in-vivo constraints. Also, we will consider the use of inhomogeneous magnetic field and the field gradient as an actuation input, in order to increase the number of degrees of freedom and the achievable robot shapes. Finally, the fundamental properties of M-CTRs will be further investigated. With the current robot structure, the FTL deployment strategy will be validated and analyzed experimentally. The magnet at the tip will be replaced by magnets with different magnetization directions, in order to expand the deployment, orientation and stability modulation capabilities. On the long term, other architectures of
M-CTRs will be studied. They will be generated by varying the number of tubes, the number and the location of the magnets along the backbone, and the mechanical, geometrical and magnetic properties of each element.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Agence Nationale de la Recherche within the Biomedical Innovation program (NEMRO ANR-14-CE17-0013), the Investissements d'Avenir (Robotex ANR-10-EQPX-44, Labex CAMI ANR-11-LABX-0004) and EUR EIPHI (Contract No. ANR-17-EURE-0002), the Swiss National Science Foundation through grant number 200020B_185039 and the ERC Advanced Grant 743217 Soft Micro Robotics (SOMBOT).

Supplemental Material

Supplemental material for this article is available online.

ORCID iDs

Quentin Peyron https://orcid.org/0000-0003-2628-1464

Quentin Boehler https://orcid.org/0000-0001-7573-3135

References

Amanov E, Granna J and Burgner-Kahrs J (2017) Toward improving path following motion: Hybrid continuum robot design. In: 2017 IEEE international conference on robotics and automation (ICRA), ■■■, ■■■, ■■■. pp. 4666–4672. DOI: 10.1109/ICRA.2017.7989542

Amanov E, Nguyen TD and Burgner-Kahrs J (2019) Tendon-driven continuum robots with extensible sections—A model-based evaluation of path-following motions. In: The International Journal of Robotics Research. SAGE Publications Ltd STM. DOI: 10.1177/0278364919886047

Antunes P, Lima H, Monteiro J, et al. (2008) Elastic constant measurement for standard and photosensitive single mode optical fibres. *Microwave and Optical Technology Letters* 50(9): 2467–2469. DOI: 10.1002/mop.23660. https://onlinelibrary.wiley.com/doi/pdf/10.1002/mop.23660

Bailly Y and Amirat Y (2005) Modeling and control of a hybrid continuum active catheter for aortic aneurysm treatment. In: Proceedings of the 2005 IEEE international conference on robotics and automation, ■■■, ■■■, ■■■. pp. 924–929. DOI: 10.1109/ROBOT.2005.1570235

Baz A, Poh S, Ro J, et al. (1992) Active control of nitinol-reinforced composite beam. In: HS Tzou and GL Anderson (eds), *Intelligent Structural Systems, Solid Mechanics and Its Applications*. Dordrecht: Springer Netherlands, pp. 169–212. ISBN. DOI: 10.1007/978-94-017-1903-2_4.

Bergeles C, Gosline AH, Vasilyev NV, et al. (2015) Concentric tube robot design and optimization based on task and anatomical constraints. *IEEE Transactions on Robotics* 31(1): 67–84. DOI: 10.1109/TRO.2014.2378431

Berlin AA and Sussman GJ (1994) *Towards Intelligent Structures: Active Control of Buckling*. Cambridge, USA: PhD Thesis, Massachusetts Institute of Technology.

Burdette EC, Rucker DC, Prakash P, et al. (2010) The ACUSITT ultrasonic ablator: the first steerable needle with an integrated interventional tool. In: *Medical Imaging 2010: Ultrasonic Imaging, Tomography, and Therapy*. International Society for Optics and Photonics, Vol. 7629, p. 76290V. DOI: 10.1117/12.845972

Burgner J, Swaney PJ, Lathrop RA, et al. (2013) Debulking from within: a robotic steerable cannula for intracerebral hemorrhage evacuation. *IEEE Transactions on Biomedical Engineering* 60(9): 2567–2575. DOI: 10.1109/TBME.2013.2260860

Burgner-Kahrs J, Rucker DC and Choset H (2015) Continuum robots for medical applications: a survey. *IEEE Transactions on Robotics* 31(6): 1261–1280. DOI: 10.1109/TRO.2015.2489500.

Charreyron SL, Gabbi E, Boehler Q, et al. (2019) A magnetically steered endolaser probe for automated panretinal photocoagulation. *IEEE Robotics and Automation Letters* 4(2): xvii–xxiii. DOI: 10.1109/LRA.2018.2888894.

Chautems C, Lyttle S, Boehler Q, et al. (2018) Design and evaluation of a steerable magnetic sheath for cardiac ablations. *IEEE Robotics and Automation Letters* 3(3): 2123–2128. DOI: 10.1109/LRA.2018.2809546

Chikhaoui MT, Benouhiba A, Rougeot P, et al. (2018) Developments and control of biocompatible conducting polymer for intracorporeal continuum robots. *Annals of Biomedical Engineering* 46(10): 1511–1521. DOI: 10.1007/s10439-018-2038-2

Chikhaoui MT, Rabenorosoa K and Andreff N (2016) Kinematics and performance analysis of a novel concentric tube robotic structure with embedded soft micro-actuation. *Mechanism and Machine Theory* 104: 234–254. DOI: 10.1016/j.mechmachtheory.2016.06.005

Comber DB, Pitt EB, Gilbert HB, et al. (2016) Optimization of curvilinear needle trajectories for transforamenal hippocampotomy. In: *Operative Neurosurgery*. DOI: 10.1227/NEU.0000000000001361

Couture T and Szewczyk J (2018) Design and experimental validation of an active catheter for endovascular navigation. *Journal of Medical Devices: American Society of Mechanical Engineers Digital Collection* 12(1): 011003. DOI: 10.1115/1.4038334.Publisher

Dahroug B, Tamadazte B, Weber S, et al. (2018) Review on otological robotic systems: toward microrobot-assisted cholesteatoma surgery. *IEEE Reviews in Biomedical Engineering* 11: 125–142. DOI: 10.1109/RBME.2018.2810605

Dupont PE, Lock J, Itkowitz B, et al. (2010) Design and control of concentric-tube robots. *IEEE Transactions on Robotics* 26(2): 209–225. DOI: 10.1109/TRO.2009.2035740
Dwyer G, Chadebecq F, Amo MT, et al. (2017) A continuum robot and control interface for surgical assist in fetoscopic interventions name: IEEE robotics and automation letters. *IEEE Robotics and Automation Letters* 2(3): 1656–1663. DOI: 10.1109/LRA.2017.2679902.Conference.

Edelmann J, Petruska AJ and Nelson BJ (2017) Magnetic control of continuum devices. *The International Journal of Robotics Research* 36(1): 68–85. DOI: 10.1177/0278364916683443

Farajollahi M, Woehling V, Plesse C, et al. (2016) Self-contained tubular bending actuator driven by conducting polymers. *Sensors and Actuators A: Physical* 249: 45–56. DOI: 10.1016/j.sna.2016.08.006

Fichera L, Dillon NP, Zhang D, et al. (2017) Through the eustachian tube and beyond: a new miniature robotic endoscope to see into the middle earconference name: IEEE robotics and automation letters. *IEEE Robotics and Automation Letters* 2(3): 1488–1494. DOI: 10.1109/LRA.2017.2668468

Garriga-Casanovas A and Rodriguez y Baena F (2018) Complete follow-the-leader kinematics using concentric tube robots. *The International Journal of Robotics Research* 37(1): 197–222. DOI: 10.1177/0278364917746222

Gilbert HB, Hendrick RJ and Webster III RJ (2016) Elastic stability of concentric tube robots: a stability measure and design test. *IEEE Transactions on Robotics* 32(1): 20–35. DOI: 10.1109/TRO.2015.2500422

Gilbert HB, Neimat J and Webster RJ (2015) Concentric tube robots as steerable needles: achieving follow-the-leader deployment. *IEEE Transactions on Robotics* 31(2): 246–258. DOI: 10.1109/TRO.2015.2394331.

Girerd C, Lihoreau T, Rabenorosoa K, et al. (2018) In vivo inspection of the olfactory epithelium: feasibility of robotized optical biopsy. *Annals of Biomedical Engineering* 46: 1951–1961. DOI: 10.1007/s10439-018-2076-9.

Goldman RE, Bajo A, MacLachlan LS, et al. (2013) Design and performance evaluation of a minimally invasive telerobotic platform for transurethral surveillance and intervention. *IEEE Transactions on Biomedical Engineering* 60(4): 918–925. DOI: 10.1109/TBME.2012.2226031

Gosline AH, Vasilyev NV, Butler EJ, et al. (2012) Percutaneous intracardiac beating-heart surgery using metal MEMS tissue approximation tools. *The International Journal of Robotics Research* 31(9): 1081–1093. DOI: 10.1177/0278364912443718

Ha J, Park FC and Dupont PE (2016) Elastic stability of concentric tube robots subject to external loads. *IEEE Transactions on Biomedical Engineering* 63(6): 1116–1128. DOI: 10.1109/TBME.2015.2483560

Ha J, Park FC and Dupont PE (2017) Optimizing tube precurvature to enhance the elastic stability of concentric tube robots. *IEEE Transactions on Robotics* 33(1): 22–37. DOI: 10.1109/TRO.2016.2622278

Hendrick RJ, Gilbert HB and Webster RJ (2015a) Designing snap-free concentric tube robots: A local bifurcation approach, 2015 IEEE international conference on robotics and automation (ICRA), Seattle, ■■■, ■■■, pp. 2256–2263. DOI: 10.1109/ICRA.2015.7139498

Hendrick RJ, Mitchell CR, Herrell SD, et al. (2015b) Hand-held transendoscopic robotic manipulators: a transurethral laser prostate surgery case study. *The International Journal of Robotics Research* 34(13): 1559–1572. DOI: 10.1177/0278364915585397

Iasnii VP and Junga R (2018) Phase transformations and mechanical properties of the nitinol alloy with shape memory. *Materials Science* 54(3): 406–411. DOI: 10.1007/s11003-018-0199-7

Jayender J, Patel RV and Nikumb S (2009) Robot-assisted active catheter insertion: algorithms and experiments. *The International Journal of Robotics Research* 28(9): 1101–1117. DOI: 10.1177/0278364909103785

Kang B, Kojcev R and Sinibaldi E (2016) The first interlaced continuum robot, devised to intrinsically follow the leader. *Plos One* 11(2): e0150278. DOI: 10.1371/journal.pone.0150278

Kavoi BM and Jameela H (2011) Comparative morphometry of the olfactory bulb, tract and stria in the human, dog and goat. *International Journal of Morphology* 29: 939–946.

Kim J-S, Lee D-Y, Kim K, et al. (2014) Toward a solution to the snapping problem in a concentric-tube continuum robot: grooved tubes with anisotropy. In: IEEE international conference on robotics and automation (ICRA), ■■■, ■■■, ■■■, pp. 5871–5876. DOI: 10.1109/ICRA.2014.6907723

Kim Y, Parada GA, Liu S, et al. (2019) Ferromagnetic soft continuum robots. *Science Robotics* 4(33): eaax7329. DOI: 10.1126/scirobotics.aax7329

Kratchman LB, Bruns TL, Abbott JJ, et al. (2017) Guiding elastic rods with a robot-manipulated magnet for medical applications. *IEEE Transactions on Robotics* 33(1): 227–233. DOI: 10.1109/TRO.2016.2623339.

Kundrat D, Schoob A, Kahrs LA, et al. (2016) First study on a monolithically manufactured variable-length continuum robot with hybrid actuation. *CURAC. Bern*: 87–92.

Lazarus A, Miller JT and Reis PM (2013) Continuation of equilibrium and stability of slender elastic rods using an asymptotic numerical method. *Journal of the Mechanics and Physics of Solids* 61(8): 1712–1736. DOI: 10.1016/j.jmps.2013.04.002

Li Z, Wu L, Ren H, et al. (2017) Kinematic comparison of surgical tendon-driven manipulators and concentric tube manipulators. *Mechanism and Machine Theory* 107: 148–165. DOI: 10.1016/j.mechmachtheory.2016.09.018.

Lock J, Laing G, Mahvash M, et al. (2010) Quasistatic modeling of concentric tube robots with external loads. In: IEEE/RSJ international conference on intelligent robots and systems (IROS), ■■■, ■■■, ■■■, pp. 2325–2332. DOI: 10.1109/IROS.2010.5651240.

Maghooa F, Stilli A, Noh Y, et al. (2015) Tendon and pressure actuation for a bio-inspired manipulator based on an antagonistic principle. In: IEEE international conference on robotics and automation (ICRA), ■■■, ■■■, ■■■, pp. 2556–2561. DOI: 10.1109/ICRA.2015.7139542

Mitros Z, Khadem M, Seneci C, et al. (2018) Towards modelling multi-arm robots: eccentric arrangement of concentric tubes. In:
2018 7th IEEE international conference on biomedical robotics and biomechatronics (Biorob), ■■■, ■■■, ■■■, pp. 43–48. DOI: 10.1109/BIOROB.2018.8488091. ISSN: 2155-1782

Mowatt G, N'Dow J, Vale L, et al. (2011) Photodynamic diagnosis of bladder cancer compared with white light cystoscopy: Systematic review and meta-analysis. *International Journal of Technology Assessment in Health Care* 27(1): 3–10. DOI: 10.1017/S0266462310001364.

Nikravesh PE, Wehage RA and Kwon OK (1985) Euler parameters in computational kinematics and dynamics. Part 1. *Journal of Mechanisms, Transmissions, and Automation in Design* 10(3): 358–365. DOI: 10.1115/1.3260722

Olson E. (2011) April tag: a robust and flexible visual fiducial system. In: 2011 IEEE international conference on robotics and automation, ■■■, ■■■, ■■■, pp. 3400–3407. DOI: 10.1109/ICRA.2011.5979561

Olszewska E, Wagner M, Bernal-Sprekelsen M, et al. (2004) Etiopathogenesis of cholesteatoma. *European Archives of Oto-Rhino-Laryngology* 261(1): 6–24. DOI: 10.1007/s00405-003-0623-x.

Palmer D, Cobos-Guzman S and Axinte D (2014) Real-time method for tip following navigation of continuum snake arm robots. *Robotics and Autonomous Systems* 62(10): 1478–1485. DOI: 10.1016/j.robot.2014.05.013.

Petruska AJ, Edelmann J. and Nelson B. J. (2017) Model-Based Calibration for Magnetic Manipulation. *IEEE Transactions on Magnetics* 53(7): 1–6. DOI: 10.1109/TMAG.2017.2653080.

Petruska AJ and Nelson BJ (2015) Minimum bounds on the number of electromagnets required for remote magnetic manipulation. *IEEE Transactions on Robotics* 31(3): 714–722. DOI: 10.1109/TRO.2015.2424051.

Peyron Q, Boehler Q, Rabenorosoa K, et al. (2018) Kinematic analysis of magnetic continuum robots using continuation method and bifurcation analysis. *IEEE Robotics and Automation Letters* 3(4): 3646–3653. DOI: 10.1109/LRA.2018.2855803.

Peyron Q, Rabenorosoa K, Andreff N, et al. (2019a) A numerical framework for the stability and cardinality analysis of concentric tube robots: Introduction and application to the follow-the-leader deployment. *Mechanism and Machine Theory* 132: 176–192. DOI: 10.1016/j.mechmachtheory.2018.10.020

Peyron Q, Rabenorosoa K, Andreff N, et al. (2019b) *Orientability Evaluation of Concentric Tube Robots Deployed in Natural Orifices*. Rennes: SURGETICA Conference, pp. 99–101.

Riojas KE, Hendrick RJ and Webster RJ (2018) Can elastic instability be beneficial in concentric tube robots?. *IEEE Robotics and Automation Letters* 3(3): 1624–1630. DOI: 10.1109/LRA.2018.2800779.

Rucker DC, Jones BA and Webster III RJ (2010) A geometrically exact model for externally loaded concentric-tube continuum robots. *IEEE Transactions on Robotics* 26(5): 769–780. DOI: 10.1109/TRO.2010.2062570

Schaeffner M, Götz B and Platz R (2016) Active buckling control of a beam-column with circular cross-section using piezo-elastic supports and integral LQR control. *Smart Materials and Structures* 25(6): 065008. DOI: 10.1088/0964-1726/25/6/065008

Simaan N, Kai Xu K, Wei Wei A, et al. (2009) Design and integration of a telerobotic system for minimally invasive surgery of the throat. *The International Journal of Robotics Research* 28(9): 1134–1153. DOI: 10.1177/0278364908104278

Singh K, Tipton CR, Han E, et al. (2013) Magneto-elastic buckling of an Euler beam. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences* 469(2155): 20130111. DOI: 10.1098/rspa.2013.0111

Swaney PJ, York PA, Gilbert HB, et al. (2016) Design, fabrication, and testing of a needle-sized wrist for surgical instruments. *Journal of Medical Devices* 11(1): 014501–014509. DOI: 10.1115/1.4034575

Tunay I (2004) Modeling magnetic catheters in external fields. In: The 26th annual international conference of the IEEE engineering in medicine and biology society, San Francisco, CA, USA, 1–5 September 2004, p. 2006. DOI: 10.1109/IEMBS.2004.1403591

Tunay I (2013) Spatial continuum models of rods undergoing large deformation and inflation. *IEEE Transactions on Robotics* 29(2): 297–307. DOI: 10.1109/TRO.2012.2232532

Webster RJ, Okamura AM and Cowan NJ (2006) Toward active cannulas: miniature snake-like surgical robots. In: 2006 IEEE/RSJ international conference on intelligent robots and systems, ■■■, ■■■, ■■■, pp. 2857–2863. DOI: 10.1109/IROS.2006.282073

Wu L, Crawford R and Roberts J (2017) Dexterity analysis of three 6-DOF continuum robots combining concentric tube mechanisms and cable-driven mechanisms. *IEEE Robotics and Automation Letters* 2(2): 514–521. DOI: 10.1109/LRA.2016.2645519
# Appendix

## A Appendix: Index to multimedia extensions

The following multimedia content is associated with this paper:

<table><caption>**Table 6.** Multimedia content.</caption><thead><tr><th>Extension</th><th>Media type</th><th>Description</th></tr></thead><tbody><tr><td>1</td><td>Video</td><td>Demonstration of the M-CTR degrees of freedom</td></tr><tr><td>2</td><td>Video</td><td>Simulation of follow-the-leader deployments</td></tr><tr><td>3</td><td>Video</td><td>Demonstration of orientability through orientation control</td></tr><tr><td>4</td><td>Video</td><td>Demonstration of active stability management</td></tr></tbody></table>

## B Appendix: Development of the kineto-static model

B.1 Expression of the matrix $L_q$. The matrix $L_q$ is obtained from the standard relation between the backbone curvature and its orientation

$$ ^0\mathbf{R}_b' = ^0\mathbf{R}_b[\mathbf{b}\mathbf{u}]_\times \quad (21) $$

where $\mathbf{b}$ denotes the derivative with respect to the arc-length $s$ and $^0\mathbf{R}_b$ is the rotation matrix from $\mathcal{R}_b$ to $\mathcal{R}_0$. The backbone orientation is represented in the kineto-static model with a quaternion $\mathcal{Q}$, which can be decomposed in a scalar $\mathcal{Q}_r$ and a vector $\mathcal{Q}_u$ such that $\mathcal{Q} = [\mathcal{Q}_r \mathcal{Q}_u^T]^T$. The relation between $\mathcal{Q}$ and $^0\mathbf{R}_b$ is given by the Rodrigues formula

$$ ^0\mathbf{R}_b = (2\mathcal{Q}_r^2 - 1)\mathcal{I} + 2(\mathcal{Q}_u\mathcal{Q}_u^T + \mathcal{Q}_r[\mathcal{Q}_u]_\times) \quad (22) $$

Injecting (22) in (21) leads to the geometrical constraint (2) and the expression of $L_q$

$$ \begin{aligned} C_{q2} &= u_B - 2L_q Q' = 0 \\ L_q &= [-Q_u Q_r \mathcal{I} - [Q_u]_\times] \end{aligned} \quad (23) $$

B.2 Mechanical equilibrium. The equilibrium equations of the M-CTR are obtained by computing the total potential energy due to the robot deformation and the forces and torques applied on the backbone.

The potential energy due to the robot deformation and stored in subsection *jl*, where $n_j$ tubes interact, is the sum of the energy due to the deformation of each tube. Considering that the material composing the tubes is linear and isotropic, the energy depends on the stiffness matrix of each tube, which is expressed in $\mathcal{R}_i$ for tube $i$

$$ \mathbf{K}_i = \begin{bmatrix} k_{bi} & 0 & 0 \\ 0 & k_{bi} & 0 \\ 0 & 0 & k_{ti} \end{bmatrix} \quad (24) $$

Considering that subsection *jl* is bounded by arc-lengths $s_j$ and $s_{j+1}$ along the backbone, the potential energy for the subsection is written

$$ \begin{aligned} E_e &= \int_{s=s_j}^{s_{j+1}} w_e ds \\ w_e &= \sum_{i=1}^{n_j} \left( \mathbf{R}_b^b \mathbf{u} - \hat{\mathbf{u}} \right)^T \mathbf{K}_i \left( \mathbf{R}_b^b \mathbf{u} - \hat{\mathbf{u}} \right) \end{aligned} \quad (25) $$

The potential energy due to gravity forces $\mathbf{f}_g$, magnetic forces $\mathbf{f}_m$ and torques $\boldsymbol{\tau}_m$ is derived by distributing these forces and torques along the robot sections. According to Lazarus et al. (2013), the potential energy induced by these forces and torques is written

$$ \begin{aligned} E_f &= - \int_{s=s_j}^{s_{j+1}} w_f ds \\ w_f &= \boldsymbol{\tau}_q^T \mathbf{Q} + \left( ^0\mathbf{f}_g + ^0\mathbf{f}_m \right)^T \mathbf{p} \end{aligned} \quad (26) $$

where $\boldsymbol{\tau}_q$ is the projection of $\boldsymbol{\tau}_m$ in the quaternion space. Finally, the total potential energy stored in the subsection *jl* is the sum of (25) and (26) and is written:

$$ \begin{aligned} E &= \int_{s=s_j}^{s_{j+1}} w ds \\ w &= w_e - w_f \end{aligned} \quad (27) $$

The Euler-Lagrange formula is then applied to the potential energy $E$ considering the kinematic constraints (1-3). It consists in expressing the following differential equation for each robot state coordinate $x_k$ in $\mathbf{x} = [\theta_1 \quad \dots \quad \theta_n u_x u_y \mathbf{Q}^T \mathbf{p}^T]$

$$ \begin{aligned} & \frac{\partial}{\partial s} \left( \frac{\partial w}{\partial x'_k} - \frac{\partial w}{\partial x_k} \right) + \frac{\partial}{\partial s} \left\langle \frac{\partial \mathbf{C}_q}{\partial x'_k}, \boldsymbol{\lambda}_q \right\rangle - \left\langle \frac{\partial \mathbf{C}_q}{\partial x_k}, \boldsymbol{\lambda}_q \right\rangle \\ & + \frac{\partial}{\partial s} \left\langle \frac{\partial \mathbf{C}_p}{\partial x'_k}, \boldsymbol{\lambda}_p \right\rangle - \left\langle \frac{\partial \mathbf{C}_p}{\partial x_k}, \boldsymbol{\lambda}_p \right\rangle = 0 \end{aligned} \quad (28) $$
where $C_q = [C_{q1} \ C_{q2}]$ and $\langle a,b \rangle$ denotes the scalar product between vectors $a$ and $b$. This leads to the following equilibrium equation

$$ egin{cases} k_{ii} \theta_i' - \hat{u}_i^{i'T} e_3 + k_{ii}' (\theta_i' - i \hat{R}_i^T e_3) \quad \dots \\ -{}^b \mathbf{u} \frac{\tau \dot{\mathbf{R}}_i^b}{\partial \theta_i} \mathbf{K}_i^i \hat{u}_i + \tau_{mz} = 0 \\ -\sum_{i=1}^{n_j} k_{bi} {}^b \mathbf{u} + \sum_{i=1}^{n_j} {}^i \mathbf{R}_b {}^i \mathbf{K}_i \hat{u}_i - \lambda_q = 0 \\ -2 \mathbf{L}_q^T \lambda_q' + \tau_q + \mathbf{S}_q \lambda_q + \mathbf{S}_p \lambda_p = 0 \\ \lambda_p' + {}^0 \mathbf{f}_g + {}^0 \mathbf{f}_m = 0 \end{cases} \quad (29) $$

The relation between $\tau_q$ and $\tau_m$ is obtained by comparing the equilibrium equations with the static model of CTR developed in Lock et al. (2010). In that work, the special Cosserat rod equations are used, leading to the following expression of the backbone curvature

$$ egin{aligned} & -\sum_{i=1}^{n_j} k_{bi} {}^b \mathbf{u} + \sum_{i=1}^{n_j} {}^b \mathbf{R}_i {}^i \mathbf{K}_i \hat{u}_i + {}^b \mathbf{m} = 0 \\ & {}^b \mathbf{m}' + {}^b \tau + [{}^b \mathbf{u}]_\times {}^b \mathbf{m} + [\mathbf{v}]_\times {}^b \mathbf{R}_0 {}^0 \mathbf{n} = 0 \\ & {}^0 \mathbf{n}' + {}^0 \mathbf{f} = 0 \end{aligned} \quad (30) $$

where $\mathbf{m}$ and $\mathbf{n}$ are, respectively, the internal torques and forces due to distributed torques ${}^b \tau$ and forces $\mathbf{0}f$, and $\mathbf{v} = [0 \ 0 \ 1]^T$. We conclude from the first and last line of equation (30) that the Lagrange multipliers $\lambda_q$ and $\lambda_p$ correspond, respectively, to the internal moment $-{}^b \mathbf{u}$ and to the internal force $\mathbf{0}v$. As a consequence, multiplying both side of the second line of equation (30) by $2\mathbf{L}_q^T$, and considering distributed magnetic torques ${}^b \tau = {}^b \mathbf{R}_0 {}^0 \mathbf{\tau}_m$, leads to

$$ -2\mathbf{L}_q^T \lambda_q' + 2\mathbf{L}_q^T {}^b \mathbf{R}_0 {}^0 \mathbf{\tau}_m + \dots + 2\mathbf{L}_q^T \left( -[{}^b \mathbf{u}]_\times \lambda_q + [\mathbf{v}]_\times {}^b \mathbf{R}_0 {}^b \lambda_p \right) = 0 \quad (31) $$

By comparing equations (31) and (29), we deduce that $\tau_q$ relates to the magnetic torque $\mathbf{0}\tau_m$ by the relation

$$ \tau_q = 2\mathbf{L}_q^T {}^b \mathbf{R}_0 {}^0 \mathbf{\tau}_m \quad (32) $$

As a side note, by injecting equations (21) and (22) into equation (31), we can also verify the expressions of $S_q$ and $S_p$

$$ egin{aligned} S_q &= -2\mathbf{L}_q^T [{}^b \mathbf{u}]_\times = -2 \left[ 2 \frac{\partial \mathbf{L}_q^T}{\partial s} \mathbf{Q} \right] \\ S_p &= 2\mathbf{L}_q^T [\mathbf{v}]_\times {}^b \mathbf{R}_0 = \frac{\partial {}^0 \mathbf{R}_q e_3}{\partial \mathbf{Q}} \end{aligned} \quad (33) $$

# C Appendix: Experimental measure of tip orientation ad angular displacement

## C.1 Computation of the tip pose.
The front camera A5 (see Figure 3) gives a relative estimation of the robot tip pose

with respect to its initial configuration. In order to compute the tip pose, the measured position and angular displacements are thus reported on the initial tip pose predicted by the kineto-static model. As a consequence, the position around which the tip is rotated and the axis of rotation differs from the measurements by a constant and unknown translation and rotation. We assume in the following that the mean position of the tip during the experiment, denoted $\bar{p}$ corresponds to $p_d$. To enforce this, a constant translation $t$ is applied at each measured tip position

$$ t = p_d - \bar{p} \quad (34) $$

We also assume that the tip is rotated about the desired axis. The rotation axis measured by the front camera is evaluated by constructing the service sphere, finding the orientation plane which best fits the experimental points, and computing the normal vector $z_p$ to this plane. The plane is computed using the least-squares method. A constant rotation can then be defined between the desired rotation axis, $x_0$ or $z_0$, and $z_p$. This constant rotation is applied to each measured orientation.

The camera pose is chosen in order to maximize the visibility of April tags during the rotation cycle. It was not possible to find a pose for which the April tag was always visible during the rotation about $z_0$, due to the large angular displacement of the M-CTR tip. As a consequence, the camera pose was chosen to measure half of the rotation cycle, that is, for $a_z \ge 0$. The robot showed a symmetric behaviour during the other half of the cycle. The symmetric results for $a_z \le 0$ are constructed, then presented in the following figure.

## C.2 Angular displacement and rotation error.
The angular displacement achieved during tip orientation and the rotation error are defined according to the points and vectors introduced in Figure 18. The angular displacement is computed by projecting the experimental points $p$ in the orientation plane. The projected points are denoted $p^*$. The points marking orientation limits, denoted $p_1^*$ and $p_2^*$, are then extracted as presented on Figure 19. The angular displacement $\Delta a$ is finally computed using the relation

$$ \Delta a = \arcsin \left( \frac{\| \mathbf{Op}_1^* \|^2 + \| \mathbf{Op}_2^* \|^2 - \| p_1^* p_2^* \|^2}{2 \| \mathbf{Op}_1^* \| \| \mathbf{Op}_2^* \|} \right) \quad (35) $$

The rotation error, denoted $\xi$, is defined as the angle between the points on the service sphere before and after projection. It is defined by

$$ \xi = \arcsin \left( \frac{\| \mathbf{Op}^* \|}{\| \mathbf{Op} \|} \right) \quad (36) $$