Instrumented bio-inspired cable-driven compliant continuum robot: static modeling and experimental evaluation

Abstract

Energy efficiency is inherent for autonomous robotic device. Snakes are well known for their ability to low energy consumption when swimming. However, the swimming know-how is poorly understood. Designing a snake robot inspired by snakes as a tool to find out the swimming energy efficiency crucial point will lead to the development of hyper efficient undulating locomotors. In this article, we introduce a four tendons driven continuum robot made of bio-inspired compliant vertebrae to assess the energy consumption of a planar and a spatial snake motion. The tendon-driven continuum robot constitutes the head–neck part of a locomotor snake robot. A static modeling coupled with an optimization method was implemented to generate bio-inspired motions recorded on snake swimming head. A friction model describing the friction between cables and the disks is investigated and compared to a frictionless model. The proposed prototype is equipped with exteroceptive sensors to record motion and proprioceptive sensors to measure cable forces applied at the tip of the robot. Hence, the work of the forces, thus the energy required to execute a trajectory are computed and analyzed. The energy is introduced as a key criterion to assess the swimming motion of a locomotor snake robot.

Nomenclature

$m$

Total number of tendons

$c$

Designated tendon among $m$ tendons

$v$

Total number of vertebrae. Each vertebra consists of a disc and a universal joint

$i$

Designated vertebra among $v$ vertebrae

$n$

Length of a trajectory

$j$

Designated position among $n$ positions of a trajectory

$\theta$

Bending angle

$\varphi$

Angle defining the direction of curvature

$\mathbf{d}(\mathbf{s}_{\mathbf{i}})$

Translation vector

$r_{i}$

Radius of curvature

$l_{i}$

Length of the vertebra $i$

$k_{i}$

Subsegment curvature

$\Gamma _{i}$

Curvature of the backbone along $x_{i-1}$ axis

$\beta _{i}$

Curvature of the backbone along $y_{i-1}$ axis

${ }^{i-1}{\mathbf{T}}{}_{i}$

Homogeneous transformation matrix from frame $i$ to $i-1$

${ }^{i-1}{\mathbf{T}}{}_{0}$

Homogeneous transformation matrix from frame $0$ to $i-1$

$\mathbf{R}_{z}(\varphi _{i})$

Rotation matrix around $z_{i-1}$ axis

$\mathbf{R}_{y}(\theta _{i})$

Rotation matrix around $y_{i-1}$ axis

$\mathbf{R}_{z}(-\varphi _{i})$

Rotation matrix around $z_{i-1}$ axis

${ }^{i}{\mathbf{p}_{\boldsymbol{t} }}{}_{i,c }$

Tendon distribution array on disk $i$

${ }^{i}{O}{}_{i}$

Center of disk $i$ in the frame $i$

${ }^{i-1}{O}{}_{{ucj_{i}}}$

Center of gravity of the universal compliant joint between two disks in the frame $i-1$

${ }^{i}{p}{}_{i,c}$

Tendon position on disk $i$

$r_{c}$

Radius from center of disk $i$

${{ }^{0}{\mathbf{F}}{}_{\mathbf{ucj}}}_{g,i}$

Gravity vector of the universal compliant joint $i$ in the frame $0$

${{ }^{0}{\mathbf{F}}{}_{\mathbf{d}}}_{g,i}$

Gravity vector of the disk $i$ in the frame $0$

${m_{f}}_{i}$

Mass of the flexible element $i$

${m_{d}}_{i}$

Mass of the disk $i$

${ }^{i-1}{\mathbf{F}_{\mathbf{d}}}{}_{{g_{i}}}$

Force from gravity of disk $i$ in the frame $i-1$

${ }^{i-1}{\mathbf{M}_{\mathbf{d}}}{}_{{g_{i}}}$

Moment from the gravity of the disk $i$ in the frame $i-1$

${ }^{i-1}{\mathbf{F}_{\mathbf{ucj}}}{}_{{g_{i}}}$

Force from gravity of the universal compliant joint $i$ in the frame $i-1$

${ }^{i-1}{\mathbf{M}_{\mathbf{ucj}}}{}_{{g_{i}}}$

Moment from the gravity of the universal compliant joint $i$ in the frame $i-1$

${ }^{v-1}{\mathbf{F}_{\mathbf{ext}}}{}$

External forces applied at the tip of the robot in the frame $i-1$

${ }^{v-1}{\mathbf{M}_{\mathbf{ext}}}{}$

Moment induced by external forces in the frame $i-1$

${ }^{0}{\mathbf{F}}{}_{e}$

External force applied at the tip of the robot in the global frame

${ }^{i-1}{\mathbf{F}_{i,c}}{}$

Driving tendons tensions applied by tendon $c$ on the disk $i$ in the frame $i-1$

$\mathbf{F}_{i,c}$

Tension magnitude along the tendon $c$ between the frame $i-1$ and $i$

$\mathbf{F}_{i+1,c}$

Tension magnitude along the tendon $c$ between the frame $i$ and $i+1$

${ }^{i-1}{\mathbf{M}_{i,c}}{}$

Moment induced by tendon tension in the frame $i-1$

${ }^{i-1}{\mathbf{F}}{}_{{O_{i}}}$

Resultant forces in the frame $i-1$

${ }^{i-1}{\mathbf{M}}{}_{{O_{i}}}$

Resultant moment in the frame $i-1$

$\mu$

Friction coefficient

${\unicode[Arial]{x03C3}}$

Angle between the tendon and the vertebra

${F_{f}}_{\sigma }$

Frictional force in function of angle $\sigma$

${F_{N}}_{\sigma }$

Normal force in function of angle $\sigma$

$M_{1}$

Fixed mass to determiner friction coefficient

$M_{2}$

Adjusted mass to determine friction coefficient

${ }^{i-1}{\mathbf{U}}{}_{i-1,c}$

Unit vector between the frames $i$ and $i-1$ in the frame $i-1$

${ }^{i-1}{\mathbf{U}}{}_{i+1,c}$

Unit vector between the frames $i$ and $i+1$ in the frame $i-1$

${ }^{i-1}{\mathbf{N}_{i,c}}{}$

Normal force in the frame $i-1$

${ }^{i-1}{{ }^{ }{\mathbf{n}}{}}{}-1_{i}$

Unit normal vector in the frame $i-1$

$\mathbf{K}_{\mathbf{x}}$

Lateral stiffness

$\mathbf{K}_{\mathbf{y}}$

Dorsal/ventral stiffness

${ }^{i-1}{\mathbf{M}}{}_{\mathbf{ben}{\mathbf{d}_{i}}}$

Bending moment of the universal compliant joint $i$

${ }^{i-1}{\tau }{}_{{x_{i}}}$

Torsion torque around y axis in the frame $i-1$

${ }^{i-1}{\tau }{}_{y}$

Torsion torque around x axis in the frame $i-1$

${ }^{i-1}{\tau }{}_{{z_{i}}}$

Twist moment

$G$

Shear modulus

$I_{0}$

Polar moment of inertia

$\varepsilon _{i}$

Twisting angle of vertebrae $i$

$\rho _{w}$

Water density

$V_{i}$

Volume of vertebrae $i$

${ }^{i-1}{\mathbf{f}}{}_{\mathbf{a},0}$

Archimed force applied to the center of gravity of vertebra $i$ in the frame $0$

${ }^{i-1}{\mathbf{F}}{}_{\mathbf{a},i}$

Archimed force in the frame $i-1$

${ }^{i-1}{\mathbf{M}}{}_{\mathbf{a},i}$

Moment induced by Archimed force in the frame $i-1$

${ }^{i-1}{\mathbf{M}}{}_{\mathbf{l}{O_{i}}}$

Moment induced by resultant forces in the frame $i-1$

$\alpha _{{lim_{i}}}$

Maximum bending angle of vertebra $i$

$\mathbf{F}_{c,v,j}$

Optimal tendon force of cable $c$ applied to vertebra $v$ for the position $j$

$\boldsymbol{V}$

Design vector for the complete trajectory

$\mathbf{va}\mathbf{r}_{\mathbf{n}}$

Design vector for each position of the trajectory

$E_{{x_{0}}}$

Relative errors between desired position and optimized position along $\overrightarrow {x_{0}}$

$E_{{y_{0}}}$

Relative errors between desired position and optimized position along $\overrightarrow {y_{0}}$

$e(\mathbf{va}\mathbf{r}_{\mathbf{j}})$

Optimization criteria

$\alpha _{f}$

Weighing factor

$x_{{{d_{{O_{v}}}}_{j}}},y_{{{d_{{O_{v}}}}_{j}}}$

Desired position of point $O_{v}$ derived from measured trajectories on biological snakes

$x_{{{c_{{O_{v}}}}_{j}}},y_{{{c_{{O_{v}}}}_{j}}}$

Computed positions taken at the tip of the continuum robot

$f$

Objective function

$C_{1},C_{2}$

Penalty function

$cst_{1},cst_{2}$

Scalars used to penalize in the objective function

${l_{\textit{tendons}}}_{c,j}$

Tendons lengths for a given position $j$ of the defined trajectory

$\gamma _{(1,2),j}$

Pulley angles of servomotors 1 and 2 for a given position $j$ of the defined trajectory

$W$

Work of a force in a plane

$\Delta d$

Displacement between two positions

$F$

Force applied in a plane to compute the work

$\theta _{w}$

Angle between the force and the displacement to compute the work

$\mathbf{F}_{\mathbf{e}{\mathbf{q}_{j}}}$

Equivalent force at position $j$ of the defined trajectory

$\mathbf{W}_{\mathbf{e}{\mathbf{q}_{j}}}$

Equivalent fork at each point of the trajectory

$\Delta \mathbf{d}_{(j,j+1)}$

Displacement between position $j$ and position $j+1$

${F_{eq}}_{{{x_{0}}_{j}}}$

Equivalent force applied by tendons on vertebra $v$ along $x_{0}$ axis

${F_{eq}}_{{{y_{0}}_{j}}}$

Equivalent force applied by tendons on vertebra $v$ along $y_{0}$ axis

${F_{eq}}_{{{z_{0}}_{j}}}$

Equivalent force applied by tendons on vertebra $v$ along $z_{0}$ axis

$\Delta d_{{{x_{0}}_{j}}}$

Displacement along $x_{0}$ axis

$\Delta d_{{{y_{0}}_{j}}}$

Displacement along $y_{0}$ axis

$\Delta d_{{{y_{0}}_{j}}}$

Displacement along $y_{0}$ axis

1. Introduction

Biomimetic underwater robots are mostly based on continuum robots’ solution and aim to replicate the locomotion of marine creatures, such as snakes. The purpose is to attain the same level of efficiency in swimming as observed in a natural environment. Most continuum robots are based on the use of tendons, which is a commonly employed actuation principle. These robots are called tendon-driven continuum robots.

1.1. Bio-inspired robot locomotors

Lead by biology, bio-inspired locomotors work toward reproducing specificities of biological animals such as locomotion. Undulatory locomotion of snakes is widely investigated. Crawling modes (concertina, sidewinding [1]) are mainly studied and reproduced by snake robots [2] on various substrates. Anguilliform swimming mode was also slightly studied and reproduced on snake robots since Hirose [3]. A growing interest on the capability of snake robots to evolve on complex and chaotic fields led to a wide variety of design of snake robots but still adopt the same assembly. The Amphibot II [4], the mamba robot [5], the HUMRS snake robot [6], and the ACM-R5 [7] are anguilliform underwater locomotors made of a series of identical rigid modules either linked by pivot joints for planar motions or assembled to form universal active joints. This assembly mode offers the possibility to control the position of the whole body by actuating each joint. However, musculoskeletal behavior is strongly degraded due to a lack of elasticity in active joints. In the case of compliant body adaptation, snake robots embed proprioceptive sensors such as force sensors [8, 9]. The body shape control strategy allows the snake robot to adapt to external environment. This type of design do not take into account the biological snake features such as the passive compliancy, the continuum aspect of snakes while swimming, the body shape and the weight. Also, the snake robot actuation is far from biological snakes [10]. An alternative design to avoid broken lines snake robot are the continuum snake robots. Well known as snake robotic arms, they are capable to adapt to external environments [11] and produce fluid motions in addition of being lightweight. Their actuation methods constraints robots being manipulators. Actuators are located at the base for tendon-driven continuum robots [12, 13] or required additional components such as pump for fluid-like actuators [14]. Other actuation methods were employed such as shape memory alloy Nitinol cables [15] used in underwater arm. Sensors complete the robot for assessment of external environment where [16] embed force sensors at the tip. In the case of embedded actuators, cable-driven continuum robots are embedded into locomotors robots [17–19] to play the role of a manipulator onto a locomotor. Continuum robots are also designed to become a locomotor as the soft snake robot [17]. Yet, to the best knowledge of authors, none of cable-driven continuum locomotor robot such as snake robot was investigated.

1.2. Robot modeling

Based on our robot design, the modeling method was oriented on continuum robots. The various modeling methods of cable-driven compliant continuum robots are depicted in literature to predict robot deformation and position. The kinematic model is a simplistic model that considers the robot deformations based of piecewise constant curvature theory [19, 20]. In fact, the robot deformation is approximated by a series of tangent curvature arcs even though curvatures are not fully circular. Although no interactions with external environment are considered. Contrary to kinematic modeling, static modeling predicts the robot deformation according to tendons forces, tendon routine, gravity, external forces [21–23], and any other forces according to environment [15]. The Cosserat static modeling using the variable curvature representation is widely used in continuum robots [24] but poorly investigated for locomotors [25]. The constant curvature static model assumes a constant arc curvature of each subsegment [23]. The assembly of subsegment constitutes the continuum robot. The robot tip position is obtained throughout various formalisms such as arc geometry and D-H parameters. Other modeling methods are described in literature such as pseudo-rigid body representation (PRB) [26, 27] and modal approaches [23].

1.3. Energy measurements

As the robot is a locomotor, the energy consumption is the crux, and robots design are led to be energy efficient. Energy efficiency was mainly explored by improving locomotion methods [28]. Locomotion methods were compared, and simulations demonstrate the anguilliform locomotion is the most efficient one for given tasks such as inspection. Gait patterns are optimized based on the minimization of power consumption and the maximization of forward velocity to minimize the energy consumption [29]. Also, the criterion was investigated to improve joint elasticity of such robots [30] based on a multi-rigid segment design (poly articulated robots). Meanwhile various investigations were conducted on gait pattern to minimize energy consumption, mechanical design is poorly investigated for this. Hence, we investigate the work for several trajectories to characterize the consumed energy using force gauge.

1.4. Current study

In this article, we investigate a tendon-driven continuum robot in the air and underwater using the constant curvature static modeling specifically adapted to the robot design. An optimization method based on static modeling was described for the robot path planning exported from biological snakes’ measurements. The robot is instrumented to allow comparison between numerical and experimental trajectories and tendon forces. The work of the robot is assessed for several trajectories and discussed as an optimization criterion to reduce the energy consumption of tendon-driven locomotors.

The main contributions of this paper are listed below:

• • A novel locomotor snake robot architecture based on compliant continuum robot design and bio-inspired by snakes.

• • A prototype equipped with exteroceptive as well as proprioceptive sensors to assess the energy consumption.

• • A static modeling based on the robot architecture coupled with an optimization approach to generate bio-inspired trajectories recorded on biological snakes swimming.

• • A new criterion introduced for robot designing: the work of forces-oriented criterion. This criterion assesses the energy consumption of the robot for any trajectories over within its workspace.

The paper is structured as follows. Section 2 depicts the static modeling associated with the snake robot design. We introduce the modeling of the torsion beam, the friction effects, and the Archimedes thrust as the robot is immersed. Section 3 deals with the trajectory characterization and generation. Snake motions are measured and duplicated through an optimization approach. Numerical and experimental evaluation is featured in section 4 and conclusion sum-up the presented work and the main contributions.

2. Robot modeling

2.1. Robot design

The proposed continuum compliant snake robot is made of an assembly of 45 compliant vertebrae all different from one to another in terms of length, weight, stiffness, and amplitudes. All the vertebrae are composed of a compliant universal joint and a disk where tendons can go through. A rigid module defined as the actuator (acting on tendons) is embodied in the snake robot. Thus, according to tendon routine, one or multiple segments of the robot can be controlled actively and passively. In the present article, the head–neck part of the robot is assumed to be active and actuated by tendons. As a first study, this snake part will be modeled, and experimentations will be performed. This model will be extended to the entire snake body to control the whole body and reproduce anguilliform motion.

The present robot is a part of the complete swimming snake robot as shown in Fig. 1(a) [31]. It is made of a continuum compliant part (Fig. 1(b)) and a rigid module (Figure 1(b)). The rigid part actuates two pairs of orthogonal and antagonistic cables using a Dynamixel^® 2XC-430W-250T servomotor, which represents the equivalent of muscles. The braided coated steel cables, with a 0.5 mm diameter, are used as tendons that link the muscles and vertebrae. The continuum module (Fig. 1(b)) is as monolithic 3D compliant part printed in nylon^® (PA12) with the selective laser sntering (SLS) method. It is assimilated as an assembly of $v=8$ vertebrae where the articulation between two disks reproduces the natural motion through a sized compliant universal joint. The joint compliancy is comparable to that of a biological snake, exhibiting a similar level of stiffness. Each vertebra (compliant universal joint and disk) is different in terms of disk diameter, length, stiffness, and limits (or mechanical stops) from one to another such as those on the biological snakes.

Figure 1. $(a)$ Swimming snake robot skeleton. $(b)$ Continuum robot (neck-head). $(c)$ Force sensors. $(d)$ experimental force sensors position.

We designed, manufactured, and assembled our proprioceptive force sensors. These latter are embedded in the robot enabling the assessment of the energy consumption of the robot for a given gait pattern. Each force sensor measures tendon tension through an aluminum deformable test piece on which a Wheatstone gauge (Vishay^® S5020) is bonded and allowing deformation measure. Knowing the displacement and the force for a given gait pattern, the corresponding mechanical energy consumption of the robot can be determined by computing the work of the force. To this end, we integrated force sensors at the tip of the robot (Fig. 1(c)). Each sensor is made of a bending test piece equipped with a Wheatstone bridge. Four force sensors are attached to the four tendons on the last disk that represents the snake head (Fig. 1(b)). Force sensors are calibrated using a set of mass varying from 30 to 250 g and a linear interpolation linking the measured tension to the force is defined.

The continuum robot is assimilated into a part of the complete snake robot. In this study, the robot represents the head–neck part. Assuming the swimming snake robot is entirely actuated through cables, the adapted solution for the continuum robot (neck-head) can be easily extended to the entire part of the robot from the head to the cloaca for actuation.

2.2. Static modeling

2.2.1 Kinematic modeling: constant curvature assumption

The robot vertebrae are assumed to locally bend as a constant curvature [19] meanwhile the curvature is varying from a vertebra to another as they are all different in terms of size and stiffness.

The constant curvature assumption was applied to each subsegment (bending between two vertebrae). Arc geometry formulation was used for vertebra kinematic modeling as described by [23]. Position of $\mathbf{d}(\mathbf{s}_{\mathbf{i}})$ at any point in the arc in function of the base in the $i-1$ frame is given by

(1) \begin{align} \mathbf{d}\left(\mathbf{s}_{\mathbf{i}}\right)=\left[\begin{array}{c} x_{i}\\[3pt] y_{i}\\[3pt] z_{i} \end{array}\right]=\left[\begin{array}{c} r_{i}\cdot \cos \left(\varphi _{i}\right)\cdot \left(1-\cos \left(\theta \left(s_{i}\right)\right)\right)\\[3pt] r_{i}\cdot \sin \left(\varphi _{i}\right)\cdot \left(1-\cos \left(\theta \left(s_{i}\right)\right)\right)\\[3pt] r_{i}\cdot \sin \left(\theta \left(s_{i}\right)\right) \end{array}\right] \end{align}

where $r_{i}=\frac{1}{k_{i}}$ is the radius of curvature. The subsegment curvature $k_{i}$ is defined by Eq. (2), where $, \beta _{i}, \Gamma _{i}$ are, respectively, the curvature of the backbone along $\overrightarrow {x_{i-1}}$ and $\overrightarrow {y_{i-1}}$ . $\varphi _{i}$ is defined by Eq. (3).

(2) \begin{align} k_{i}=\sqrt{\beta _{i}^{2}+\Gamma _{i}^{2}} \end{align}

(3) \begin{align} \varphi _{i}=atan2\left(\Gamma _{i},\beta _{i}\right) \end{align}

Given $\mathbf{d}_{i}$ for $s_{i}=l_{i}$ , the homogeneous transformation matrix from $i-1$ to $i$ is described by [23] as follows:

(4) \begin{align} { }^{i-1}{\mathbf{T}}{}_{i}=\left[\begin{array}{cc} \mathbf{R}_{z}\left(\varphi _{i}\right)\mathbf{R}_{y}\left(\theta _{i}\right)\mathbf{R}_{z}\left(-\varphi _{i}\right) & \ \mathbf{d}_{i}\\ 0 & \ 1 \end{array}\right] \end{align}

The direct geometrical model links the actuator space, giving the pulley angles, to the configuration space and to the cartesian space at the robot tip [19].

2.2.2 Forces acting on the continuum robot

Given the kinematic modeling, the constant curvature static modeling predicts the robot bending and the robot tip position. The constant curvature static modeling is described by [21–23]. The four tendons, $m=4$ , are distributed around the backbone at a distance $r_{c}$ from the center of each disk. This distance, $r_{c}$ , can vary as the snake shape (and disks diameters) is varying all along the body [31]. In this study, $r_{c}$ remain constant. The distribution of tendons on a disk $i$ is described by [22] as follows:

(5) \begin{align} { }^{i}{\mathbf{p}_{\mathbf{t}}}{}_{i,c}=r_{c,i}\cdot \left[\cos \left(\alpha _{m}\right) \sin \left(\alpha _{m}\right)\ 0\right], with\ \alpha\ \epsilon \left[0,\frac{\pi }{2}, \pi, \frac{3\pi }{2}\right] \end{align}

As defined by [21], the acting forces in the continuum robot are coming from the gravity, friction, and external forces. The gravity vector defined in the global frame for the mass of the universal compliant joint $(\text{flexible element})i$ , $m_{uc{j_{i}}},$ and the mass of the disk $i$ , $m_{{d_{i}}}$ are given by

(6) \begin{align} { }^{0}{\mathbf{F}_{\mathbf{ucj}}}_{i}=\left[0\ 0\,{m_{ucj}}_{i}\cdot g\ 0\right]^{\mathrm{T}} \end{align}

(7) \begin{align} { }^{0}\mathbf{F}_{\mathbf{d}_{i}}=\left[0\ 0\,{m_{d}}_{i}\cdot g\ 0\right]^{\mathrm{T}} \end{align}

The force $\mathbf{F}_{{\mathbf{d}_{i}}}$ and moments $\mathbf{M}_{{\mathbf{d}}_{i}}$ due to the gravity of the disk are given by

(8) \begin{align} { }^{i-1}\mathbf{F}_{\mathbf{d}_{i}}={ }^{i-1}{\mathbf{T}}{}_{0}\cdot { }^{0}{\mathbf{F}}{}_{i} \end{align}

(9) \begin{align} { }^{i-1}\mathbf{M}_{{\mathbf{d}}_{i}}=\overrightarrow {{{ }^{i-1}{O}{}_{i-1}}{ }^{i-1}{O}{}_{i}}\times { }^{i-1}\mathbf{F}_{\mathbf{d}_{i}} \end{align}

Note that ${ }^{i-1}{\mathbf{T}}{}_{0}$ is the homogeneous transformation matrix. The force ${\mathbf{F}_{\mathbf{ucj}}}_{i}$ and moments ${\mathbf{M}_{\mathbf{ucj}}}_{g}$ due to gravity of the universal compliant joint (flexible element) are given by

(10) \begin{align} { }^{i-1}{\mathbf{F}_{\mathbf{ucj}}}_{i}={ }^{i-1}{\mathbf{T}}{}_{0}\cdot { }^{0}{\mathbf{F}_{\mathbf{ucj}}}_{i} \end{align}

(11) \begin{align} { }^{i-1}{\mathbf{M}_{\mathbf{ucj}}}_{i}=\overrightarrow {{ }^{i-1}{\mathbf{O}}{}_{i-1}{ }^{i-1}{\mathbf{O}_{g}}_{i}}\times { }^{i-1}{\mathbf{F}_{\mathbf{ucj}}}_{i} \end{align}

With ${ }^{i-1}{\mathbf{O}_{\mathbf{ucj}}}_{i}$ the gravity center of the compliant universal joint between two disks. In case of external forces applied at the robot tip, the forces $\mathbf{F}_{\mathbf{ext}}$ and the moment $\mathbf{M}_{\mathbf{ext}}$ induced by the force $\mathbf{F}_{\mathbf{ext}}$ are described by

(12) \begin{align} { }^{v-1}{\mathbf{F}_{\mathbf{ext}}}{}={ }^{v}{\mathbf{T}}{}_{0}\cdot { }^{0}{\mathbf{F}}{}_{e} \end{align}

(13) \begin{align} { }^{v-1}{\mathbf{M}_{\mathbf{ext}}}{}=\overrightarrow {{ }^{v-1}{O}{}_{v-1}{ }^{v-1}{O}{}_{v}}\times { }^{v-1}{\mathbf{F}_{\mathbf{ext}}}{} \end{align}

where ${ }^{0}{\mathbf{F}}{}_{e}$ is the contact force applied to the robot body tip in the frame $R_{0}$ . The driving tendon tensions ${ }^{i-1}{\mathbf{F}}{}_{i,c}$ , triggered by the tendon $c$ at disk $i$ , are the sum of the two driving force vectors $F_{i,c}$ and $F_{i+1,c}$ (Fig. 2(a)) lead by their respective unit vector described as follows and given by [21] and [23]:

Figure 2. $(a)$ Continuum robot pattern. $(b)$ Kinematic diagram of the continuum robot. $(c)$ Subsegment kinematic diagram.

(14) \begin{align} { }^{i-1}{\mathbf{F}_{i,c}}{}=F_{i,c}\cdot \frac{\overrightarrow {{}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i-1,c}}{}}}{\overrightarrow {\left\| {}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i-1,c}}{}\right\| }}+F_{i+1,c}\cdot \frac{\overrightarrow {{}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i+1,c}}{}}}{\overrightarrow {\left\| {}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i+1,c}}{}\right\| }} \end{align}

(15) \begin{align} { }^{i-1}{\mathbf{M}_{i,c}}{}={ }^{i-1}{\mathbf{p}_{\mathbf{t}_{i,c}}}{}\times { }^{i-1}{\mathbf{F}_{i,c}}{} \end{align}

Note that when $i=v$ , the force reduces to:

(16) \begin{align} { }^{i-1}{\mathbf{F}_{i,c}}{}=F_{{i_{c}}}\cdot \frac{\overrightarrow {{}^{v-1}{p_{v,c}}{}{}^{v-1}{p_{v-1,c}}{}}}{\overrightarrow {\left\| {}^{v-1}{p_{v,c}}{}{}^{v-1}{p_{v-1,c}}{}\right\| }} \end{align}

The equilibrium equations are defined at the $i-1$ local frame for each subsegment. Then, the resultant moments and forces at vertebra $i=v$ and at point $O_{v-1}$ in the frame $v-1$ are described as follow:

(17) \begin{align} { }^{v-1}{\mathbf{F}}{}_{{O_{v}}}=\sum _{c}^{m}{ }^{v-1}{\mathbf{F}_{v,c}}{}+{ }^{v-1}{\mathbf{F}_{\mathbf{cu}}}{}_{{g_{v}}}+{ }^{v-1}{\mathbf{F}_{\mathbf{d}}}{}_{{g_{v}}} \end{align}

(18) \begin{align} { }^{v-1}{\mathbf{M}}{}_{{O_{v}}}=\sum _{c}^{m}{ }^{v-1}{\mathbf{F}_{v,c}}{}+{ }^{v-1}{\mathbf{M}_{\mathbf{cu}}}{}_{{g_{v}}}+{ }^{v-1}{\mathbf{M}_{\mathbf{d}}}{}_{{g_{v}}} \end{align}

If $i=v$ , the external force $\mathbf{F}_{\mathbf{ext}}$ and moment $\mathbf{M}_{\mathbf{ext}}$ should be included to the equilibrium equations.

2.2.3 Tendons frictions

When the continuum robot is achieving a trajectory, the vertebrae bend according to their stiffness. In the case of discrete interactions (Fig. 2(a)), forces from tendons are locally applied to each disk which led to friction between the tendons and the associated disk. This reflects that the friction coefficient will vary over the continuum robot from the motor to the end effector and during the trajectory execution. To extract the relation between friction coefficient and bending angle, a test bench was specially designed and manufactured as shown in Fig. 3(a).

Figure 3. $(a)$ Friction test bench. $(b)$ Friction coefficient $\mu$ according to deflection angle.

A braided steel-coated cable of 0.5 mm diameter is used for the experiment. The cable goes through a hole of 0.75 mm in the disk. In addition, the cable is guided by two ball bearings to limit the friction and to vary the $\sigma$ angle (Fig. 3(a)). A set of 10 fixed mass $M_{1}$ is set for a range of 12 angles varying from 30 degrees to 85 degrees. The variable mass $M_{2}$ is adjusted for cable balance according to the angle. A linear regression is applied for each angles according to the mass. Thus, a friction coefficient is extracted for each angle (Fig. 3(b)). The friction is described by Coulomb law and defined the coefficient $\mu$ following the Coulomb friction formula as follows:

(19) \begin{align} \mu \left(\sigma \right)=\frac{{F_{f}}_{\sigma }}{{F_{N}}_{\sigma }} \end{align}

With $F_{f}$ the frictional force defined by

(20) \begin{align} {F_{f}}_{\sigma }=\left| M_{1}-M_{2}\right| \end{align}

And the normal force $F_{N}$ defined by

(21) \begin{align} {F_{N}}_{\theta }=M_{1}+M_{2} \end{align}

$M_{1}$ and $M_{2}$ are the fixed mass and the adjusted mass (Fig. 3(a)), respectively. Exponential regression Eq. (22) was the most adapted to describe the friction coefficient $\boldsymbol{\mu }_{i,c}$ in regard of the bending angle ${\unicode[Arial]{x03C3}} .$ Applied to the continuum robot, the coefficient is depicted as following:

(22) \begin{align} \boldsymbol{\mu }_{i,c}=0.689e^{-0.027{\boldsymbol{\sigma }_{i,c}}} \end{align}

$\boldsymbol{\sigma }_{i,c}$ is the local bending angle (Fig. 2(a)) of the vertebrae determined using unit vectors ${ }^{i-1}{\mathbf{U}}{}_{i-1,c}$ and ${ }^{i-1}{\mathbf{U}}{}_{i+1,c}$ :

(23) \begin{align} { }^{i-1}{\mathbf{U}}{}_{i-1,c}=\frac{\overrightarrow {{}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i-1,c}}{}}}{\overrightarrow {\left\| {}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i-1,c}}{}\right\| }} \end{align}

(24) \begin{align} { }^{i-1}{\mathbf{U}}{}_{i+1,c}=\frac{\overrightarrow {{}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i+1,c}}{}}}{\overrightarrow {\left\| {}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i+1,c}}{}\right\| }} \end{align}

(25) \begin{align} \boldsymbol{\sigma }_{i,c}=\cos ^{-1} \left(\frac{{}^{i-1}{\mathbf{U}}{}_{i-1,c}\cdot {}^{i-1}{\mathbf{U}}{}_{i+1,c}}{\left\| {}^{i-1}{\mathbf{U}}{}_{i-1,c}\right\| \left\| {}^{i-1}{\mathbf{U}}{}_{i+1,c}\right\| }\right) \end{align}

The tension magnitude $F_{i+1,c}$ is updated according to the friction between the disks and the tendons. As depicted in ref. [21], the resulting tension magnitude can be calculated as

(26) \begin{align} F_{i+1,c}=F_{i,c}-\boldsymbol{\mu }_{i,c }\left| { }^{i-1}{\mathbf{N}_{i,c}}{}\right| \end{align}

where ${ }^{i-1}{\mathbf{N}_{i,c}}{}$ is the normal force applied by tendon tension on the considered disk $i$ defined on the plane $(O_{i},\overrightarrow {x_{i}},\overrightarrow {y_{i}})$ and direction is $(\overrightarrow {p_{ic},O_{i}})$ :

(27) \begin{align} { }^{i-1}{\mathbf{N}_{i,c}}{}={ }^{i-1}{\mathbf{F}_{i,c}}{}-\left({ }^{i-1}{\mathbf{F}_{i,c}}{}\cdot { }^{i-1}{\mathbf{n}}_{i}\right) \end{align}

where ${ }^{i-1}{\mathbf{n}}_{i}$ is the unit normal vector of plane $(O_{i},\overrightarrow {x_{i}},\overrightarrow {y_{i}})$ [21]. Thus, the net force ${ }^{i-1}{\mathbf{F}_{i,c}}{}$ is updated with the computed tension magnitude $F_{i+1,c}$ . Simulations were performed with the friction coefficient determined experimentally. By applying a range of forces similar to experimental tests ( $F_{0}\in [0,0.7]\,(N)$ ), we show that the friction force is very small and can be neglected during the numerical gait pattern generation (Fig. 4). One note that the friction coefficient is defined in the air without any lubrication. Thus, it is assumed that when the robot is immersed, the coefficient decreases even more due to water lubrication.

Figure 4. Numerical simulation with implementation of friction coefficient. $(a)$ Deformed continuum robot. $(b)$ Comparison of bending angles along the 8 vertebrae (7 angles). $(c)$ Friction coefficient distribution.

2.2.4 Deformation of the compliant vertebra: torsion beam modeling

The continuum compliant robot is designed to reproduce the snake swimming undulations by adopting an elastic behavior similar to that of biological snakes. Thus, the torsion beam diameters are optimally sized in orthogonal directions (both on dorsal/ventral side and lateral side) according to stiffness repartition along biological snakes as described by Gautreau et al. [32, 31]. Note that $K_{{y_{i}}}$ is associated to dorsal/ventral stiffness while $K_{{x_{i}}}$ is associated to the lateral stiffness. Thus, the robot bending will vary heterogeneously according to the applied force. This behavior requires static modeling to take into account the tendons tensions and external forces.

(28) \begin{align} \mathbf{K}_{\mathbf{x}}=\left[K_{{x_{1}}},\ldots, K_{{x_{v}}}\right]^{T} \end{align}

(29) \begin{align} \mathbf{K}_{\mathbf{y}}=\left[K_{{y_{1}}},\ldots, K_{{y_{v}}}\right]^{T} \end{align}

The design of each vertebra is then similar to a compliant universal joint where a comprehensive description of the system is provided in ref. [32]. The beams compliances are provided by the material flexibility (PA12). The beams are deformed in torsion with small deformations so that each vertebra locally bend to make the continuum robot bending under tendon tensions (Fig, 2(c)). The resulting large deformation of the continuum robot (bending) is allowed by small deformations of each beam in torsions. Small deformations enable to adopt a linear elastic behavior for beams torsion as follows:

(30) \begin{align} { }^{i-1}{\mathbf{M}}{}_{\mathbf{ben}{\mathbf{d}_{i}}}=\left[\begin{array}{c} { }^{i-1}{\tau }{}_{{x_{i}}}\\ { }^{i-1}{\tau }{}_{{y_{i}}}\\ { }^{i-1}{\tau }{}_{{z_{i}}} \end{array}\right] \end{align}

where ${ }^{i-1}{\mathbf{M}}{}_{\mathbf{ben}{\mathbf{d}_{i}}}$ is the torsion moment induced by the torsional stiffness of the torsion beam. ${ }^{i-1}{\tau }{}_{{x_{i}}}$ and ${ }^{i-1}{\tau }{}_{{y_{i}}}$ are the torque obtained from the vertebrae length $l_{i}$ and $\Gamma _{i}$ and $\beta _{i}$ , respectively, the curvature around $y$ axis and $x$ axis. Note that the curvatures $\Gamma _{i}$ and $\beta _{i}$ are the optimized parameters:

(31) \begin{align} { }^{i-1}{\tau }{}_{{x_{i}}}=K_{{x_{i}}}\cdot \theta _{{x_{i}}}=K_{{x_{i}}}\cdot l_{i}\cdot \Gamma _{i} \end{align}

(32) \begin{align} { }^{i-1}{\tau }{}_{{y_{i}}}=K_{{y_{i}}}\cdot \theta _{{y_{i}}}=K_{{y_{i}}}\cdot l_{i}\cdot \beta _{i} \end{align}

(33) \begin{align} { }^{i-1}{\tau }{}_{{z_{i}}}={ }^{i-1}{T}{}_{i}\cdot \left[0\ 0\ G_{i}\cdot I_{{0_{i}}}\cdot \frac{\varepsilon _{i}}{l_{i}}\ 0\right]^{T} \end{align}

where $l_{i}$ is defined by $l_{i}=L_{{11_{i}}}+L_{{12_{i}}}$ . $\varepsilon _{i}$ is the twisting angle of vertebra $i$ . ${ }^{i-1}{\tau }{}_{{z_{i}}}$ is the twist moment induced around z axis described by [21] and [23]. In the case of the present robot, the twisting moment has not been taken into account in the design phase. Thus, ${ }^{i-1}{\tau }{}_{{z_{i}}}$ is neglected to avoid any twisting moment in numerical simulations.

2.2.5 Archimedes force modeling

The robot is intended to be a locomotor swimming on the $(O_{0},\overrightarrow {x_{0}},\overrightarrow {y_{0}})$ plane. Assuming the skeleton is immersed, the Archimedes force is opposed to gravity and acts on the robot bending angle. For a configuration in which the robot is horizontally setup (Fig. 5), the Archimedes thrust (force and moment) are defined according to the volume of the vertebrae $V_{i}$ , as follows:

Figure 5. Experimental setup of the robot subject to gravity and Archimedes force.

(34) \begin{align} { }^{i-1}{\mathbf{f}}{}_{\mathbf{a},0}=\left(\rho _{w}V_{i}\right)\cdot \left[\begin{array}{cccc} 0 & -1 & 0 & 0 \end{array}\right]^{T} \end{align}

(35) \begin{align} { }^{i-1}{\mathbf{F}}{}_{\mathbf{a},i}={ }^{i-1}{\mathbf{T}_{0}}{}\cdot { }^{i-1}{\mathbf{f}}{}_{\mathbf{a},0} \end{align}

(36) \begin{align} { }^{i-1}{\mathbf{M}}{}_{\mathbf{a},i}=\overrightarrow {{ }^{i-1}{O}{}_{i-1}{ }^{i-1}{O_{g}}_{i}}\cdot { }^{i-1}{\mathbf{F}}{}_{\mathbf{a},i} \end{align}

Thus, the net force and the net moment applied at vertebra $i=v$ are set to:

(37) \begin{align} { }^{v-1}{\mathbf{F}}{}_{{O_{v}}}=\sum _{c}^{m}{ }^{v-1}{\mathbf{F}_{v,c}}{}+{ }^{v-1}{\mathbf{F}_{\mathbf{cu}}}{}_{{g_{v}}}+{ }^{v-1}{\mathbf{F}_{\mathbf{d}}}{}_{{g_{v}}}+{ }^{v-1}{\mathbf{F}}{}_{\mathbf{a},v} \end{align}

(38) \begin{align} { }^{v-1}{\mathbf{M}}{}_{{O_{v}}}=\sum _{c}^{m}{ }^{v-1}{\mathbf{F}_{v,c}}{}+{ }^{v-1}{\mathbf{M}_{\mathbf{cu}}}{}_{{g_{v}}}+{ }^{v-1}{\mathbf{M}_{\mathbf{d}}}{}_{{g_{v}}}+{ }^{v-1}{\mathbf{M}}{}_{\mathbf{a},v} \end{align}

Note that in the case of horizontal static modeling, the gravity is along y axis $.$ A comparison of end effector position between numerical optimization and experiment for a given trajectory will be given in section 4.

2.2.6 Equilibrium equations

Given the resultant forces and net moments at vertebra $i=v$ Eqs. (17)–(18), Eqs. (37)–(38), the equilibrium equations at vertebrae $i$ for $i\in \{0,..,v\}$ are described by [21] as follow:

(39) \begin{align} { }^{i-1}{\mathbf{F}}{}_{{O_{i-1}}}=\sum _{c}^{m}{ }^{i-1}{\mathbf{F}_{i,c}}{}+{ }^{i-1}{\mathbf{F}_{\mathbf{cu}}}{}_{{g_{i}}}+{ }^{i-1}{\mathbf{F}_{\mathbf{d}}}{}_{{g_{i}}}+{ }^{i-1}{\mathbf{F}}{}_{\mathbf{a},i}+{ }^{i-1}{\mathbf{F}}{}_{{O_{i}}} \end{align}

(40) \begin{align} { }^{i-1}{\mathbf{M}}{}_{{O_{i-1}}}=\sum _{c}^{m}{ }^{i-1}{\mathbf{F}_{i,c}}{}+{ }^{i-1}{\mathbf{M}_{\mathbf{cu}}}{}_{{g_{i}}}+{ }^{i-1}{\mathbf{M}_{\mathbf{d}}}{}_{{g_{i}}}+{ }^{v-1}{\mathbf{M}}{}_{\mathbf{a},i}+{ }^{i-1}{\mathbf{M}}{}_{\mathbf{l}{O_{i}}}+{ }^{i-1}{\mathbf{M}}{}_{{O_{i}}} \end{align}

The discretized forces and moments ${ }^{i-1}{\mathbf{F}}{}_{{O_{i}}}$ and ${ }^{i-1}{\mathbf{M}}{}_{{O_{i}}}$ are the lumped forces and moments at point $O_{i}$ in the frame $i-1$ , and defined as follows:

(41) \begin{align} { }^{i-1}{\mathbf{F}}{}_{{O_{i}}}={ }^{i-1}{\mathbf{T}}{}_{i}\cdot { }^{i}{\mathbf{F}}{}_{{O_{i}}} \end{align}

(42) \begin{align} { }^{i-1}{\mathbf{M}}{}_{{O_{i}}}={ }^{i-1}{\mathbf{T}}{}_{i}\cdot { }^{i}{\mathbf{M}}{}_{{O_{i}}} \end{align}

And ${ }^{i-1}{\mathbf{F}}{}_{{O_{i}}}$ induce a moment ${ }^{i-1}{\mathbf{M}}{}_{\mathbf{l}{O_{i}}}$ in the frame $i-1$ according to the vertebra length $i$ :

(43) \begin{align} { }^{i-1}{\mathbf{M}}{}_{\mathbf{l}{O_{i}}}=\overrightarrow {{ }^{i-1}{O}{}_{i-1}{ }^{i-1}{O}{}_{i}}\cdot { }^{i-1}{\mathbf{F}}{}_{{O_{i}}} \end{align}

As described in the previous subsection, the vertebrae have a known compliance allowing the continuous robot to bend in two planes. This results in the following equilibrium equation as depicted by [21]:

(44) \begin{align} { }^{i-1}{\mathbf{M}}{}_{{O_{i-1}}}={ }^{i-1}{\mathbf{M}}{}_{\mathbf{ben}{\mathbf{d}_{i}}} \end{align}

The static model of the robot was solved using the numeric method trust-region-recursive algorithm (fsove in MATLAB) [23].

Remark: The inertia effects are not taken into account in this study and have been neglected, thus a slow motion has been considered. Static modeling is sufficient and allows for the modeling of anguiform snake swimming. This approach is in line with the conditions of our experiments.

2.3. Workspace with and without limits (mechanical stops)

Inspired from nature, mechanical stops imposed by the vertebrae structure (pre-zygapophyses and post-zygapophyses) limit spine motions around $\overrightarrow {x_{i}}$ , $\overrightarrow {y_{i}}$ , and $\overrightarrow {z_{i}}$ [31]. This structure implies a workspace limitation between two vertebrae and for a range of $v$ vertebrae. The maximum vertebrae angles $\alpha _{{lim_{i}}}$ ( $i\in [1,v]$ , vertebrae) vary along the snake body. It was measured on two Python regius and an Helicops angulatus through x-rays at the National Museum of Natural History [31]. Values are applied to the static model both on dorsal/ventral side and lateral side. By varying the tendon forces alternatively from 0 to 1.2 N to scan the workspace, we show the robot workspace is wider in the lateral side than dorsal/ventral side (Fig. 6), which limit the motions while swimming. It is assumed the lateral side undulations implies greater propulsion than dorsal/ventral undulations side.

Figure 6. Robot workspace $(a)$ Without limits, $(b)$ With mechanical stops.

We assume that the mechanical stops play a role in energy efficiency and reduce the energy consumption while swimming. These mechanical limits prevent large motions while swimming and allows concerned muscles to contract in an optimized range of forces. Note that these limits are barely never reached while swimming. In the case of swimming locomotion optimization, to reduce energy consumption, these limits should artificially be shortened to reduce large deflection of the snake body.

3. Trajectory gait pattern characterization and generation

3.1. Trajectories measurements of biological snakes (planar and volume swimming locomotion)

With the aim of reproducing the trajectories of a part of a snake body achieved by biological snakes when swimming, a motion capture analysis was performed. The motion capture methodology is described in ref. [33]. As the continuum robot design simulate the neck (first quarter of a biological snake), we focused measurements on the head and the neck.

Two clusters of markers where stick on snake dorsal side. The racer terrestrial snake (Hierophis viridiflavus) is 1 m length and weight 200 g. The experiment is performed in a swimming test bench [33]. The motion was recorded using six Qualisys infrared cameras. It is important to mention that the snake was captured opportunistically and released few days after capture (permit issued to XB, DBEC 004/2022).

Regarding the robot configuration (Fig. 1(b)), the actuator is set in the neck, and the tip of the continuum robot is assumed to be the head. Thus, the neck is fixed, and the head is mobile. Consequently, for a given swimming sequence, the neck is considered as the reference fixed frame, and the head positions is computed in the neck local frame. The head trajectory is computed in the $(O_{0},\overrightarrow {x_{0}},\overrightarrow {y_{0}})$ plane to be planar (Fig. 7(a)). Two maximum positions were sufficient for trajectory optimization as the robot deformation imposed the trajectory transition due to the design. When the head trajectory is computed in the volume (Fig. 7(b)), all the positions that compose the trajectory pattern are optimized.

Figure 7. Input trajectories for the optimisation method. $(a)$ Head planar trajectory using LANA software where the two extrema extracted [33]. (b) Head trajectory in a volume for a given swimming sequence extracted using MOCAP.

An optimization formulation is depicted in the following subsection to make the continuum robot moving accordingly to the extracted trajectories. Consider that owing to the design of the robot, optimization occurs within the plane defined by the origin $(O_{0})$ , the x-axis $(\overrightarrow {x_{0}})$ , and the y-axis $(\overrightarrow {y_{0}})$ . The transition between two positions takes place in a volume due to the specific design of the continuum robot. Moving from one position to another involves the robot tip following a trajectory reminiscent of a hyperbolic path.

3.2. Optimization problem formulation

The aim of this section is to detail the optimization methodology to identify the optimal tendon forces $\mathbf{F}_{c,v,j}$ ( $c\in [1,4]$ tendons and $j\in [1,n]$ goal positions) of the continuum robot to go through a defined position and thus achieve a desired trajectory. The heuristic-based method genetic algorithm was used to solve the optimization problem through a global approach covering a wide range of solutions. The design vector $\mathbf{V}$ can be defined for each $j$ positions and for the four tendons of the snake as follow:

(45) \begin{align} \mathbf{V}=\left[\mathbf{va}\mathbf{r}_{\mathbf{1}},\ldots, \mathbf{va}\mathbf{r}_{\mathbf{n}}\right]^{\mathbf{T}};\;\mathbf{va}\mathbf{r}_{\mathbf{j}}=\left[F_{1,v,j},F_{2,v,j},F_{3,v,j},F_{4,v,j}\right]^{\mathbf{T}} \end{align}

Where the interval existence of each tendon force is $\mathbf{F}_{c,v,j}\in [0,2.5]\,N$ . The optimization problem is based on the functions $E_{{x_{0}}}$ and $E_{{y_{0}}}$ defined as a quadratic error. Thus, the continuum robot reaches the target position. The optimization criteria of the problem are defined as follow:

(46) \begin{align} e\left(\mathbf{va}\mathbf{r}_{\mathbf{j}}\right)=\alpha _{f}\cdot E_{{x_{0}}}\left(\mathbf{va}\mathbf{r}_{\mathbf{j}}\right)+\left(1-\alpha _{f}\right)\cdot E_{{y_{0}}}\left(\mathbf{va}\mathbf{r}_{\mathbf{j}}\right) \end{align}

with,

(47) \begin{align} E_{{x_{0}}}\left(\mathbf{va}\mathbf{r}_{\mathbf{j}}\right)=\frac{1}{2}\left(x_{{{d_{{O_{v}}}}_{j}}}-x_{{{c_{{O_{v}}}}_{j}}}\right)^{2}, \mathrm{j}\in \left[1;\;\mathrm{n}\right]\ \textit{positions} \end{align}

(48) \begin{align} E_{{y_{0}}}\left(\mathbf{va}\mathbf{r}_{\mathbf{j}}\right)=\frac{1}{2}\left(y_{{{d_{{O_{v}}}}_{j}}}-y_{{{c_{{O_{v}}}}_{j}}}\right)^{2}, \mathrm{j}\in \left[1;\;\mathrm{n}\right]\ \textit{positions} \end{align}

$\alpha _{f}$ is a weighing factor where $\alpha _{f}\in [0,1]$ . Note that $x_{{{d_{{O_{v}}}}_{j}}}$ and $y_{{{d_{{O_{v}}}}_{j}}}$ are the desired positions derived from measured trajectories on biological snakes. $x_{{{c_{{O_{v}}}}_{j}}}$ and $y_{{{c_{{O_{v}}}}_{j}}}$ are the computed positions in a plan at the tip of the continuum robot.

The optimization problem is subject to constraints in order to respect the operating properties as two pairs of antagonistic cables are implanted orthogonally (Fig. 1(b)). For the proper functioning of cable, the following constraints are defined:

(49) \begin{align} if\ F_{1,v,j}\ or\ F_{4,v,j}\gt 0, \textit{respectively}\ F_{4,v,j}\ or\ F_{1,v,j}=0 \end{align}

(50) \begin{equation} \textit{if}\ F_{2,v,j}\ or\ F_{3,v,j}\gt 0, \textit{respectively}\ F_{3,v,j}\ or\ F_{2,v,j}=0\end{equation}

In fact, when a tendon is tensioned, the antagonist tendon is released and cannot exert any force. The optimization problem was formulated to identify the design vector $\mathbf{V}$ corresponding to achieve a desired trajectory. A penalty function was used to handle the constraints and to ensure that the fitness of any feasible solution is better than the non-feasible ones. Hence, the objective function $f$ to be minimized is defined as follows:

(51) \begin{equation}\min\ :\ (f(\mathbf{V})=e\left(\mathbf{V}\right)+C_{1}+C_{2})\end{equation}

Subject to

(52) \begin{align} C_{1} \begin{cases} cst_{1} & if\ F_{1,v,j}\gt\ 0\ and\ F_{4,v,j}\gt\ 0 \\ & \quad 0 \quad \textit{otherwise} \end{cases} \end{align}

(53) \begin{align} C_{2}\begin{cases} cst_{2} & if\ F_{2,v,j}\gt\ 0\ and\ F_{3,v,j}\gt\ 0 \\ & \quad 0 \quad \textit{otherwise} \end{cases} \end{align}

with,

(54) \begin{align} & u_{k}\in \left[u_{k\min }, u_{k\max }\right] ;u_{k}\in \left\{F_{1,v,j},F_{2,v,j},F_{3,v,j},F_{4,v,j}\right\}\nonumber\\ & \mathbf{V}=\left[\mathbf{va}\mathbf{r}_{\mathbf{1}},\ldots, \mathbf{va}\mathbf{r}_{\mathbf{n}}\right]^{\mathbf{T}};\;\mathbf{va}\mathbf{r}_{\mathbf{j}}=\left[F_{1,v,j},F_{2,v,j},F_{3,v,j},F_{4,v,j}\right]^{\mathbf{T}}\\& u_{k\min }=0N;\ u_{k\max }=25N\nonumber \end{align}

where $cst_{1}$ and $cst_{2}$ are high scalars used to penalize the objective function.

Once the trajectory is defined, the genetic algorithm method is used to solve the optimization problem and identify the optimal design vector. The optimization problem was solved using MATLAB 2022R, running on a Windows 10 with $10^{th}$ Gen Intel^® Core ^TM i5-10500 (3.1 Ghz) processor and 16 GB Ram. The average computation time to obtain a solution and solve the optimization problem is a function of the number of points describing a given trajectory (Table. I). The obtained results are compared with experimental ones in the following section. Simulating the robot with the optimal design vector enables computing the tendons lengths ${l_{\textit{tendons}}}_{c,j}$ and therefore the two pulley angles $\gamma _{(1,2),j}$ as depicted in ref. [19].

Table I. Optimization problems: Computation time vs trajectory size.

Remark: The computing time varies significantly among the three optimization problems due to the trajectory length (number of points that compose the trajectory). Simulation N°1 consists of 2 positions, whereas simulation N°2 and N°3 consist of 241 positions. Additionally, Archimedean forces are taken into account for simulation N°3. Adding more positions to optimize is time consuming as is adding additional forces.

4. Numerical and experimental evaluation

4.1. Material and method

4.1.1 Experimental test bench in air

Two experimental test benches were used for the experiments. For a numerical/experimental coupling, the robot was tested vertically outside of water. The actuator (“neck”) (Fig. 8) was fixed on the test bench, and the tip of the robot (“head”) was mobile and equipped with four force sensors. The motion was recorded through motion capture. The head and the neck were equipped with a cluster of three passive reflective markers (Fig. 8), whereas disks of vertebrae were equipped with only one reflective marker. Velocity and orientation are computed from MOCAP using MATLAB. Tip of the robot is equipped with four force sensors for each of the four actuated tendons. Forces were obtained using Phidget^® control panel at a recording frequency of 128 milliseconds.

Figure 8. Experimental test bench and prototype in air with the schematic set-up.

The work of the combined forces is computed at each position of the imposed trajectory. The general definition of the work $W$ of a force in a plane with an angle $\theta _{w}$ at each point $j$ is described as

(55) \begin{align} W=F\cdot \Delta d\cdot \cos \left(\theta _{w}\right) \end{align}

In the case of the continuum robot, the equivalent force is the combination of the 4 forces measured at the tip of the 4 tendons as shown in Fig. 9.

Figure 9. Equivalent force $\boldsymbol{F}_{\boldsymbol{eq}}$ applied by tendons at the tip of the robot.

The equivalent force computed at the tip of the last disk $O_{v}$ due to tendon force is defined as

(56) \begin{align} \mathbf{F}_{\mathbf{e}{\mathbf{q}_{j}}}=\left[\begin{array}{c} {F_{eq}}_{{{x_{0}}_{j}}}\\[4pt] {F_{eq}}_{{{y_{0}}_{j}}}\\[4pt] {F_{eq}}_{{{z_{0}}_{j}}} \end{array}\right]=\left[\begin{array}{c} \sum _{c=1}^{4}{F_{{x_{0}}}}_{c,j}\\[4pt] \sum _{c=1}^{4}{F_{{y_{0}}}}_{c,j}\\[4pt] \sum _{c=1}^{4}{F_{{z_{0}}}}_{c,j} \end{array}\right] \end{align}

With $c\in [1,4]$ tendons and $j\in [1,n]$ positions on the trajectory. The total works along each axis, according to the equivalent force $\mathbf{F}_{\mathbf{e}{\mathbf{q}_{j}}}$ and the displacement $\Delta \mathbf{d}_{(j,j+1)}$ , result in an equivalent work at each point of the trajectory $\mathbf{W}_{\mathbf{e}{\mathbf{q}_{j}}}$ :

(57) \begin{align} \mathbf{W}_{\mathbf{e}{\mathbf{q}_{j}}}=\left[\begin{array}{c} {F_{eq}}_{{{x_{0}}_{j}}}\\[2pt] {F_{eq}}_{{{y_{0}}_{j}}}\\[2pt] {F_{eq}}_{{{z_{0}}_{j}}} \end{array}\right]\cdot \left[\begin{array}{c} \Delta d_{{{x_{0}}_{j}}}\\[2pt] \Delta d_{{{y_{0}}_{j}}}\\[2pt] \Delta {d_{{z_{0}}}}_{j} \end{array}\right]={F_{eq}}_{{{x_{0}}_{j}}}\cdot \Delta d_{{{x_{0}}_{j}}}+{F_{eq}}_{{{y_{0}}_{j}}}\cdot \Delta d_{{{y_{0}}_{j}}}+{F_{eq}}_{{{z_{0}}_{j}}}\cdot \Delta d_{{{z_{0}}_{j}}} \end{align}

The robot is manually set vertically to begin the imposed trajectory. The starting position is identical to numerical simulation (tip of the robot is at coordinate [0,0], respectively, along $x$ axis and $z$ axis). Each trajectory is repeated three times, and the starting position is not modified between trials. The trajectories N°1 and N°2 were performed on the test bench.

4.1.2 Experimental test bench in water

In the second case, the robot is vertically immersed underwater. The motion was recorded through GoPro^® hero 9 video cameras with the linear mode to avoid image distortion and placed in front of the robot outside of water (Fig. 10). Motion of the robot is extracted using the markerless DeepLabCut software. The software was trained using the resnet50 neuronal network with three labeled videos. In order to allow the deep learning software to track the positions of specific points, we labeled several videos from our experiments. We specified which points to follow for each video. Consequently, any video taken under the same experimental conditions as the labeled ones can be analyzed to track these specific points. Each video is labeled on 40 frames. Velocity and orientation are computed using MATLAB. The force sensors were coated using epoxy to protect strain gauges from water. Forces were obtained using Phidget^® control panel at a recording frequency of 128 milliseconds. The actuator remains outside of water as it is not waterproof. The work was computed based on measured force and displacements $\Delta d_{x},\Delta d_{y}$ , and $\Delta d_{z}$ .

Figure 10. Experimental test bench and prototype in water with the schematic set-up.

In this study, the displacement can only be measured on a plane $(O_{0},\overrightarrow {x_{0}},\overrightarrow {z_{0}})$ . The equivalent work is computed using Eq. (54). All the continuum flexible modules were immersed. Similarly, to the previous experimental test bench and the experimentations, the starting position of the robot is manually set vertically, and the robot position is not modified from a played trajectory to another one. The trajectories N°1 and N°2 were carried out on the test bench.

4.1.3 Experimental test bench in water and horizontally

The last experiments focus on the continuum robot horizontally oriented as the snake robot will swim on the $(O_{0},\overrightarrow {x_{0}},\overrightarrow {z_{0}})$ plane. The Archimedes force is considered for the continuum flexible module.

The actuator is partially immersed, and a waterproof case was design to protect the Dynamixel^® 2XC-430W-250T servomotor. Two GoPro Hero 9 cameras were positioned orthogonally to record the motion. The camera 1 is positioned to record lateral undulations, and camera 2 is positioned to record the dorsal/ventral undulations, as shown in Fig. 11. Two virtual markers were set on the labeled videos (Fig. 11) using DeepLabCut software to extract the 3D movement of the tip of the robot based on the same local frame, by combining the positions of each video. The trajectory N°3 was performed for the corresponding setup.

Figure 11. Prototype set-up dedicated for horizontal experiments.

4.2. Experimental generated trajectories

4.2.1 Planar trajectory

The obtained trajectory represents the location of the snake head undulating in function of the neck on a $(O_{0},\overrightarrow {x_{0}},\overrightarrow {y_{0}})$ plane, corresponding to the lateral side oscillations. Experimental positions of the robot end-effector along $\overrightarrow {x}$ , $\overrightarrow {y}$ , and $\overrightarrow {z}$ are similar to numerical ones as depicted in Fig. 12 and Fig. 13 which confirm the truthfulness of the proposed static modeling for the given continuum robot design.

Figure 12. Position comparison between numerical result, experimental in air, and experimental in water results.

Figure 13. (a) Position errors between numerical and experimental results in air. (b) Position errors between numerical and experimental results in water.

Figure 14. $(a)$ Relative errors between experimental forces in air and numeric forces. $(b)$ Relative errors between experimental forces in water and numeric forces.

Figure 15. Planar trajectory: in the air $(a)$ Tendons forces $(b)$ Displacements $(c)$ Work; Underwater $(d)$ Tendons forces $(e)$ Displacements $(f)$ Work.

Figure 16. Vertical 3D trajectory of the robot tip positions along $x$ , $y$ and $z$ axis.

Figure 17. (a) Position errors between numerical and experimental results in air. (b) Position errors between numerical and experimental results in water.

Figure 18. Computed work for vertical 3D trajectory. $(a)$ Work on the air. $(b)$ Work on water.

Figure 19. Horizontal 3D trajectory end effector positions along $x$ , $y$ and $z$ axis.

Figure 20. Position errors between numerical and experimental results in water.

Remark 1: The error positioning can come from several reasons, the open loop control, the manual initial positioning, as well as the MOCAP analysis, whereas still in ratio to the biological study and mimicry. Snake swimming undulations are approximative and do not require high accuracy (positioning error less than 10 mm is acceptable in this study).

Measured experimental forces present a behavior closely similar to those obtained during numeric optimization. Relative errors in the air between experimental forces and numeric forces are smaller than those obtained underwater (Fig. 14). Errors are slightly high but remains acceptable according to experimental conditions.

In fact, the magnitude force order is quite similar between air, water, and numeric environment results. A comparison of forces (air, water, and numeric environment) was performed to assess the force accuracy. This latter shows that the force error is similar between experimentations in the air and the numeric and the experimental underwater vs the numerical one. In both cases, one notes that during the experimentation, force sensors mounted on nonactuated cables detect a small force variation, which can be explained by the internal friction between cables and disks. These forces are not considered in the static model and as well minimized underwater as expected due to the water effect acting as lubricant.

Remark 2: Force sensors are slightly sensitive to temperature. In experimental conditions, temperature can vary between 20°C and 35°C, which alter the force measurement. Tendons were manually pretensioned, which can lead to a variation of measured forces. Note that these results are preliminary and promising to open the way to new developments.

Given the force and the displacement for a gait pattern, the computed work shows the energy consumption in function of position (Fig. 15(c) (f)). It appears that the work is higher underwater than in air. One assumes water resistance plays a significant role, and water force should be considered in the next experimentations. Work is high in transition phases from an extremum to another. The work is null at the vertical position as it is a resting position. Work shows where and when the robot is consuming energy. This could be a major criterion for robot design to adjust the robot stiffness and tendon routine in the objective of energy consumption minimization.

4.2.2 Vertical 3D trajectory

The 3D trajectory represents the head of the snake undulating in function of the neck on a volume. The positions of the robot tip (end-effector) shown in Fig. 16 follow faithfully the optimized trajectory, and the deviations (Fig. 17) are still acceptable for biomimicry motion (Fig. 7).

The obtained results demonstrate the robot capability to reproduce an imposed gait pattern in a volume. Experimental forces were measured in the air and in water to compute equivalent work (Fig. 18). The work order of magnitude appears higher for underwater motions than for air. The water forces are additional forces that should be considered for fluid–structure interaction.

4.2.3 Horizontal 3D trajectory

This section represents the head of the snake undulating in function of the neck on a volume and horizontally. This configuration was close to biology as the prototype was ¾ immersed. The neck was fixed, and the head was moving. The recorded 3D movement shows the robot achieve faithfully the imposed motion as shown in Fig. 19 and Fig. 20.

5. Conclusion

In this article, the tendon-driven continuum robot modeling adapted to the bioinspired robot was introduced. The design of the continuum robot neck-head and its static modeling were presented to allow the robot to achieve a specific task. The trajectory gait pattern was assessed on biological snakes, and an optimization gait pattern was set to predict the robot moving under given forces. Experimental test benches and robot prototypes were presented. Experimental results were assessed and discussed in the last section. Experiments shown the static modeling feasibility of predicting the robot’s movement in the air and underwater.

Future work will be focused on the improvement of the snake robot’s efficiency by minimizing its energy consumption. For an imposed path trajectory problem, we suggest that energy minimization should serve as the objective function in an optimization problem. This optimization should include additional variables such as tendon routing and body stiffness. By mixing these variables, we aim to achieve an optimal solution that balances the imposed trajectory with the system’s efficiency and performances. Lastly, future work will be focused on the swimming undulations of the snake in a completely immersed scenario as well as the coating with an artificial skin development.