Morales Bieze

Alexandre Kruszewski
email: alexandre.kruszewski@centralelille.fr

Bruno Carrez

Christian Duriez

Thor Morales Bieze

Design, implementation and control of a deformable manipulator robot based on a compliant spine

Keywords: Soft robots, Continuum robots, Deformable manipulators, Finite Element Method

come

Design, implementation and control of a deformable manipulator robot based on a compliant spine

Thor

Introduction

Industrial inspection tasks very often involve the transportation of instrumentation equipment to hard-to-reach areas and the manipulation of this equipment 1 in constraint spaces. The requirements of such tasks present a challenge to conventional rigid robots which are very often heavy and bulky due to their design specifications. A possible solution to such applications would be a deformable continuum manipulator robot with a high degree of dexterity which is able to negotiate its way around delicate equipment in narrow spaces. Some of the desired capabilities for such robot are:

• Large workspace: Ideally, the robot would be deployed away from the area to be inspected. It is then necessary for the robot to have a considerable reach and a large task space.

• Dexterity: Due to the constraints of the task, it is desired that the robot is able to reach a position inside its task space with more than one orientation of its end-effector.

• Ability to sneak between obstacles: The robot should have a high number of degrees of freedom in order to avoid obstacles while still being able to reach the desired position/orientation with its end-effector.

• Very compliant: In some cases, the contact with the environment cannot be avoided, therefore, in order to avoid or minimize the damage to the surroundings, it is necessary that the robot is compliant to external forces.

In order to design such robot, one can turn to nature to seek for inspiration, however, the lack of design methodologies for bio-inspired continuum robots pose certain issues: First, using the same principle of functionality in different environments can cause the robot to under-perform. For example, the idea of mimicking the morphology of an octopus tentacle may seem attractive at first glance, however, in this particular case, the medium (water) in which the source of inspiration normally operates provides a key support component (buoyancy) that enables the octopus to reach and grab objects with its boneless tentacles. However, outside of the water and without a backbone, a robotic tentacle may lack the structural support to perform the required task. Second, the lack of simulation tools to validate those principles of functionality may translate in a considerable amount of iterations in the process of design and implementation. And third, the conceptual stages of a continuum manipulator are often disjointed from its modeling and control methodologies, which may result in a device that is too complex to control. In the work presented in this paper, a simulation tool is at the very core of the robot implementation, that allows for a reduced number of design iterations (reduced cost of time and materials due to failed prototyping) and that also provides a suitable control strategy for the robot for both, direct piloting and path planning.

By a rapid shortcut, one might be tempted to equate rigid robots with vertebrates and soft robots with invertebrates. Yet the skeleton is structured around a spine that is mechanically designed to ensure flexibility. The spine is an important mechanical support and the vertebra repetition interspersed with soft inter-vertebral discs allows great freedom of movement, while absorbing shocks.

Additionally, the spine protects the transmission of information through the spinal cord.

All these elements have led us to consider the idea of drawing inspiration from the skeletal spine to design a compliant manipulator robot, capable of great freedom of movement, especially for inspection tasks in confined environments. We called this manipulator Echelon 3 (see Figure 1).To the best of our knowledge, this continuum robot is the very first to be designed, modeled and controlled in a single modeling framework: the finite element method. The methodology used is generic and can be used to design other robots.The simulation-assisted design, modeling and control strategies of the robot are explained in the following sections.

Figure 1: The Echelon 3 deformable manipulator.

Related Work

Due to their compliance to external forces and their dexterity (the variety of tasks that the robot can complete, and also how well it can perform those tasks ( [START_REF] Ma | On dexterity and dexterous manipulation In Advanced Robotics[END_REF]]), continuum manipulators have become a suitable alternative to rigid manipulators in applications that involve dynamically changing environments or the direct contact with humans. Such applications include skeletal trauma treatment ( [Wilkening et al(2011)], [START_REF] Alambeigi | A curved-drilling approach in core decompression of the femoral head osteonecrosis using a continuum manipulator[END_REF]]), endoscopy ( [START_REF] Cianchetti | STIFF-FLOP surgical manipulator: mechanical design and experimental characterization of the single module In Intelligent Robots and Systems[END_REF]], [START_REF] Fraś | New STIFF-FLOP module construction idea for improved actuation and sensing In Robotics and Automation[END_REF]]), minimally invasive surgery ( [START_REF] Orekhov | Analysis and validation of a teleoperated surgical parallel continuum manipulator[END_REF]], [START_REF] Mahoney | Reconfigurable parallel continuum robots for incisionless surgery In Intelligent Robots and Systems[END_REF]]), minimally invasive inspection ( [START_REF] Tonapi | generation rope-like robot for in-space inspection[END_REF]], [START_REF] Greer | A Soft, Steerable Continuum Robot That Grows via Tip Extension Soft robotics[END_REF]]) and search and rescue ( [START_REF] Bajo | Finding lost wrenches: Using continuum robots for contact detection and estimation of contact location In Robotics and Automation (ICRA)[END_REF]]).

Continuum manipulators can be broadly classified with respect to the actuation system implemented to apply forces and torques to the backbone. Intrinsic or locally actuated designs utilize pneumatic or hydraulic actuators made of elastometers, embedded in the structure of the robot to apply forces directly to the backbone. These actuators, refereed to as pneumatic artificial muscles (PAM's) reproduce the way longitudinal muscles work in an animal trunks or tentacles. Some examples of locally actuated continuum robots include the Oct-Arm ( [START_REF] Walker | Continuum robot arms inspired by cephalopods In Unmanned Ground Vehicle Technology[END_REF]], [START_REF] Mcmahan | Field trials and testing of the OctArm continuum manipulator In Robotics and Automation[END_REF]]), the "European Octopus" ( [START_REF] Laschi | Design of a biomimetic robotic octopus arm Bioinspiration & biomimetics[END_REF]], [START_REF] Calisti | ) Design and development of a soft robot with crawling and grasping capabilities In Robotics and Automation (ICRA)[END_REF]]) and the Bionic Handling Assistant ( [START_REF] Rolf | Constant curvature continuum kinematics as fast approximate model for the Bionic Handling Assistant In Intelligent Robots and Systems[END_REF]]). A disadvantage of using pneumatic actuators is that they tend to exhibit non-linear behaviors such as hysteresis, which makes the control of the robot a complex task.

Extrinsic actuation designs employ filaments or tendons routed along the structure of the manipulator. The tendons are pulled by motors positioned at the base of the robot to apply torques to the backbone to produce bending ( [START_REF] Cojocaru | Experiments with Tentacle Robots In ISR/ROBOTIK[END_REF]], [START_REF] Ataka | Towards safer obstacle avoidance for continuumstyle manipulator in dynamic environments[END_REF]]). Walker and Hanna presented a study on a soft manipulator inspired by the elephant trunk ( [START_REF] Walker | A novel'elephant's trunk'robot[END_REF]]), that allowed the identification of critical hardware necessities, in particular the need for elastic structures and backbone support. In the work presented by [START_REF] Gravagne | Large deflection dynamics and control for planar continuum robots[END_REF]]), an incompressible elastic rod is used as a backbone to guide the bending of the manipulator. This choice of backbone has been the most used in tendon based designs [START_REF] Simaan | A dexterous system for laryngeal surgery IEEE International Conference on Robotics and Automation[END_REF]], [START_REF] Bardou | Control of a multiple sections flexible endoscopic system IEEE/RSJ International Conference on Intelligent Robots and Systems[END_REF]], [START_REF] Zhao | Design and analysis of a kind of biomimetic continuum robot IEEE International Conference on Robotics and Biomimetics[END_REF]). When the application requires more compliance on the robot, a central spring can be used as a backbone ( [START_REF] Mehling | A minimally invasive tendril robot for in-space inspection IEEE[END_REF]], [START_REF] Yoon | Active bending endoscope robot system for navigation through sinus area IEEE/RSJ International Conference on Intelligent Robots and Systems[END_REF]]), however, these designs tend to be more difficult to control due to the spring absorbing a part of the actuation forces.

Modeling and Control of Soft Manipulators

The advantages of soft continuum robots come with the trade off of complexity in their modeling and control. In order to control a soft manipulator, a transformation between the low level degrees of freedom (actuator space) and the high level degrees of freedom which are the position and orientation of the tip of the robot (task space) needs to be formulated. A very common approach towards this model is to pass through an intermediate transformation involving the geometry of the robot. Using the assumption that, after actuation is applied, each section of the robot follows a constant curvature, the state of the tip of each section of robot can be related to its curvature ( [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review The International[END_REF]]), this approach have been implemented several times in the literature ( [START_REF] Rolf | Constant curvature continuum kinematics as fast approximate model for the Bionic Handling Assistant In Intelligent Robots and Systems[END_REF]], [START_REF] Escande | Geometric modelling of multisection bionic manipulator: Experimental validation on robotinoxt In Robotics and Biomimetics (ROBIO)[END_REF]], [START_REF] Escande | Kinematic calibration of a multisection bionic manipulator[END_REF]]). When the manipulator is composed of several serially connected constant curvature sections, the same approach can be used by concatenating the transformation matrices of each section, as presented in ( [START_REF] Jones | Kinematics for multisection continuum robots[END_REF]]), this approach is also compatible with closed-loop controllers where the state of the robot is acquired by external sensors or cameras [START_REF] Marchese | Design, kinematics, and control of a soft spatial fluidic elastomer manipulator[END_REF]]). When it is assumed a deformation with constant curvature between sections, the system can then be considered within a framework approaching that of a rigid articulated model. However, these models ignore the notion of internal forces equilibrium which is very important when studying deformations: material properties, self-collision and deformations caused by the environment could all have an influence on the kinematics of a soft manipulator.

Due to their ability to by-pass the complexity of modeling the highly nonlinear behavior of soft continuum manipulators, data driven approaches have gained popularity in the literature. These approaches are based on previous knowledge of input-output data derived from experimentation and can provide accurate fast approximations to the kinematic model ( [START_REF] Rolf | Goal babbling permits direct learning of inverse kinematics IEEE Transactions on Autonomous Mental Development[END_REF]], [Rolf and Steil(2014)]). Neural network approaches have been also proposed in the context of learning the kinematics of continuum manipulators using geometry models for data recollection ( [START_REF] Giorelli | Neural network and jacobian method for solving the inverse statics of a cable-driven soft arm with nonconstant curvature[END_REF]]), or as part of low level controllers to compensate for the dynamic uncertainties of the robots ( [START_REF] Braganza | A neural network controller for continuum robots IEEE transactions on robotics[END_REF]]). While data-driven models have shown great success in real world scenarios, the fact that the learning base changes when the operation conditions of the robot change make these models limited in their application.

The need of accurately modelling the physical behavior of continuum manipulator has led to the use of the theory of Cosserat rods ( [START_REF] Renda | Dynamic model of a multibending soft robot arm driven by cables[END_REF]). Based on such model, the control of the position of the remote end of these robots has been implemented in [START_REF] Torres | Motion planning for concentric tube robots using mechanics-based models[END_REF]]. Most proposed models only consider the deformation of the central line of the robot. However, given the structure chosen for the robot presented in this paper, it is important to take into account the deformations created by the lattice structure. Yet, the deformation of lattice structure has been widely modeled in FEM using a combination of beam elements ([Przemieniecki (1985)]).

Outline of the paper

The paper presents new methodological steps for the conception, model and control design of a continuum manipulator, implemented for a compliant spine robot, actuated by tendons.

The contributions of this paper, with respect to the current state-of-the-art, are as follows:

• In the design of the manipulator, we proved that it is possible to obtain large global deformations while limiting the amplitude of the local deformations in the structure of the robot. The robot is made of 3D printed vertebrae, however, it retains a high compliance and its motion is created by deformation.

• We proved that taking the design inspiration from vertebrate animals does not undermine the soft capabilities of the continuum manipulator.

• We showed how a simulation of the FEM model of the robot is used to unify the different stages of its implementation. This allows for a reduction of iterations on the design stage and also provides a suitable modeling and control strategies for the continuum manipulator.

• We demonstrated teleoperation with haptic feedback for a more intuitive manner of control, in the case where the robot is controlled by a human operator.

The paper is organized as follows: first, the mechanical design of the deformable continuum robot is presented in the next section, then the mechanical model with 3D simulation and the control design are explained. The control is further extended for tele-operation with haptic rendering and the experimental results are discussed before the conclusion.

Design

On the skeletal spine, ligaments make the connection between the vertebrae and stabilize them. Their deformation allows for a good overall range of motion but they have also a limiting movement role to protect the link. Their local deformation is relatively small. We have kept these principles in our design to create large amplitudes of global movement from relatively small local deformations which are constrained to remain small by the design.

Another source of inspiration in the design of the manipulator is the engineering truss (i.e. a framework consisting of beams and rods that is often used in civil engineering to support a roof, bridge, or other structures.)

Hardware description

The manipulator is composed by 3 sections with a total length of 78.3cm. Each section is actuated extrinsically by 3 tendons attached at the section end points (see Figure 2 and Extension 1). By pulling the tendons in different configurations, one can change the orientation and position of the tip of each section. To provide the tension in the tendons, a Pololu 3215 DC motor actuates a timing belt based system that pulls the tendon mounted at the back end of the robot. This choice of transmission is to avoid using a pulley system where the tendon can slip out of the spool when not in tension. The actuator setup is presented in Figure 3. The backbone of the robot is made out of corrugated high density Polyethylene with an elastic module of 0.8 GPa and a tensile strength of 15 MPa. The role of the backbone in the design is twofold: to give support to the structure of the robot while preventing the shearing between the vertebrae, and to provide housing for transmission and instrumentation elements (tendon guides and camera and sensor wires).

As illustrated in Figure 4, a transmission tube made of Teflon-PTFE is allocated inside the backbone. This tube connects the actuators at the back end of the robot to the middle and distal sections of the manipulator. The internal routing of the tendons reduces the required actuator stroke and decouples the deformation of each section.

The vertebrae that form the sections of the manipulator are 3D printed using regular Polylactic Acid (PLA) with an elastic module of 3.5GPa and a tensile strength of 50MPa and the tendons are made of regular Nylon fishing line rated for 36kg of tension. As mentioned before, the main role of the vertebrae is to guide the tendons along the backbone and to limit its local deformation by means of self-collision. In the base and middle sections, the dorsal part of the vertebrae is longer in comparison to the lateral parts (see figure 5). The objective of this feature is to increase the torque applied by the tendons at the termination of the sections to help the manipulator cope with gravitational forces.

The manipulator lacks extension since the backbone is passive. To compensate the lack of extension, the manipulator is mounted on top of a sliding platform which is actuated simultaneously by 2 Nema 17 stepper motors. This provides an additional degree of freedom along the longitudinal axis of the manipulator.

Modeling and simulation

In this section, we describe the mechanical model used to calculate the kinematic model of the robot. We use the Finite Element Method (FEM) with curved beam type elements. These beams are then assembled by applying the rigid transformations corresponding to the geometry of the vertebra structure. We reproduce the internal forces of the structure, the mass, the forces that apply to the passage of cables, as well as the external forces applied to the robot. From this mechanical model we derive the direct model as well as the inverse model.

Curved beam Model

The robot structure is based on repetitive pattern (the trellis). Each element is composed of tubes and curved rods that can be modeled in a generic way by FEM using curved beam elements. The main assumption is that the beam elements have a constant section along their main axis. The sections are circular (to model the central tube like in Figure 7 (d) and(e) ) or rectangle (to model the vertebras like in Figure 7 (b) and (c)). Each beam is parameterized by 2 frames X 1 and X 2 , i.e. 2 nodes that have 6 Degree Of Freedom (DOF): 3 translations, 3 rotations. Thus X 1 , X 2 ∈ SE(3) A precise interpolation along the main axis of the element is obtained by using cubic spline curve (see Figure 6).

Mechanically, the beam model integrates on the chosen geometry an elastic law of linear behavior with common parameters such as the Young's modu- lus and the Poisson ratio. After computing the area moment of inertia of the section geometry, a stiffness matrix K b can be obtained in the local frame of each beam (a full description of how this matrix is obtained can be found in [Przemieniecki (1985)]). This 12 × 12 matrix provides the internal forces and moments f b observed on the two nodes when it is multiplied by u loc the displacement of the nodes in the local frame of the beam.

f b = K b u loc (1)
It accounts for the energies of deformation in the 6 possible directions: 1 direction of stretching, 2 directions of bending, 1 direction of torsion, 2 directions of shearing.

In our formulation, the local frame corresponds to the frame interpolated in the middle of the beam X b ∈ SE(3). To compute this frame, we use the spline interpolation of the element as explained in Figure 6.

Once X b has been computed, it provides the transformation from the world coordinate to the local frame of the beam. The rotation T b defined by this transformation will be used to rotate the matrix of the beam K b in the global frame, to obtain a corotational approach, like in [Crisfield (1990)].

u loc ∈ R 12 is composed by two displacement twists u 1,loc ∈ se(3) and Figure 6: The beam position is defined by two frames X 1 and X 2 placed at the extremities. These frames are the Nodes of the beam. Based on the position of these two frame, we place the 4 control points of the spline p 0 , p 1 , p 2 , p 3 , using the initial length L of the beam. Vectors p 0 p 1 and p 2 p 3 are respectively aligned with the x axis (in red in the picture) of X 1 and X 2 . After construction of the spline, the frame of the beam X b is interpolated so that its x axis is aligned with the spline tangent. The torsion (rotation along x) between X 1 and X 2 is also interpolated (not shown on the picture) u 2,loc ∈ se(3). Their computation is based on the motion of frames X 1 and X 2 in the reference frame of the beam X b from the initial position to the current one.

The formulation of the beam allows different types of geometrical sections and not straight initial shape. To model the robot, we used both rectangular and tubular sections as shown in Figure 7.

Assembling and definition of the degrees of freedom

The mechanical structure of the robot is a truss, which corresponds to the assembly of beams. The geometrical assembling of the truss follows a series of parameters that have been described in the Design section: number of sections, number of vertebrae per sections, size of the vertebrae, placement of the cables, etc... It has been implemented in a generic manner, allowing for testing, in simulation, different designs.

One of the key component in the modeling is the mechanical assembly of beam models. The assembly is based on the fact that the structure is composed of beam elements connected by nodes. Each node of the beam model described bellow is defined by rigid frame. However, we have distinguished the concept of a frame node , i.e. extremity of a beam element and the concept of frame DOF, (for Degree Of Freedom), i.e. number of independent unknowns that define the configuration of the robot. Between node and DOF frames, we can use rigid transformations, traditionally used in rigid robotics, to perform the assembly.

We can model that a frame node is rigidly attached to a frame DOF during assembly. Consequently, there is a rigid transformation between the DOF and the node. We can compute the rigid position of the node in the DOF frame:

DOF X node ∈ SE(3). (We have DOF X node = ( node X DOF ) -1 )
At the beam level, the force f b ∈ R 12 defined in equation 1 is applied on two rigid frames. So it is a composition of two wrenches f 1,b , f 2,b ∈ se * (3). In f 1,b , the node of application of the wrench is 1 but the orientation is given by the beam frame X b .

If the node 1 is rigidly attached to a DOF i, we have to transport the point of application of the wrench and orientate the wrench in the world space. First, we rotate the wrench to obtain it the frame of node 1

f 1,1 = 1 R b f 1,b .
Then we can transpose the wrench in the node frame (both point of application and orientation):

f i,i = Ad T iX1 * f 1,1 .
Where Ad T is the transposed Adjoint of the transformation i X 1 between frame i and frame node 1. Finally, this wrench is rotated to obtain it the frame in the world coordinates f i,0 = R T i f i,i . We gather and add all the wrenches defined on DOFs in the world coordinates to obtain the vector of internal elastic forces F(q), where q = {X i } gathers all the positions of the DOFs. The same transformations are used to compute the mass matrix M(q).

Figure 8, illustrates the placement of the Frame DOFs in the structure and the rigid transformations we have implemented between beam nodes and DOFs. To model the robot, we have used 41 beams to model the backbone and 123 beams to model the vertebrae. In total, we used 102 frame DOFs, so the size of the model is 612. The global system of equation follows the second law of shows the discretization of the structure in beam elements (in yellow). Two beams are highlighted. The actuation is modeled using constraints: H T a provides the direction of the constraint and λ is the intensity of the pulling force. (b) illustrates the assembling of the two highlighted beams. The frame node X 1 is attached to a frame DOF X i . Same for a frame node X 2 rigidly attached to frame DOF X j .

Newton:

M(q)q = R -F(q) + H T a (q)λ (2)
where R contains the known external forces like the gravity and H T a (q)λ represents the contribution of the actuators through cables. Like in [START_REF] Largilliere | Real-time control of soft-robots using asynchronous finite element modeling[END_REF]], we use Lagrange multipliers to model the actuation. It allows to let as unknown λ the tension force applied on the cable. The distribution of this force on the DOFs depends on the geometrical configuration of the cable in the structure and the current position the robot q. The way this distribution is computed is explained in more details in [START_REF] Coevoet | Software toolkit for modeling, simulation, and control of soft robots[END_REF]]. We can gather this distribution in a matrix H T a (q). Figure 8 illustrate the construction of the matrix H a for the robot.

Forward simulation

The range of the displacements applied on the beam during the deformation of the manipulator justifies the use of a non-linear corotational beam model. The model provides good results for large displacement but small strain. At each time step of the simulation we update a linearisation of the internal forces:

F(q t+h ) = F(q t ) + ∂F ∂q (q t )dq (3)
We use implicit Euler integration with these linearized forces, the time step is h and the current state is (q t , v t ) (with v the velocities of the degrees of freedom). α and β are the two coefficients of Rayleigh (dynamic viscosity and static viscosity). We obtain the following system of equations:

( 1 h + α)M + (β + h) ∂F ∂q A(qt) dv = -f (q t ) -h ∂F ∂q v t + R b(qt,vt) +H T a (q t )λ (4)
In practice, in the simulation, we used α = 1 and β = 0.1. When solved, the value of the solution dv is used to update the Euler scheme: v t+h = v t + dv and q t+h = q t + hv t+h When solving the direct simulation problem, we want to impose reference lengths δ * a (t) to the cables and compute the new state of the robot as well as the new lengths of the cable. This is only a unilateral constraint, as we can only pull on the cable (not push). The effective length δ a (q t ) depends on the position of the degrees of freedom and can be written in the actuation space via a nonlinear mapping between motion space and actuation space. This nonlinear function can also be linearized at each time step using the expression of the matrix H a (q t ) at the current step: δ a (q t+h ) = δ a (q t ) + H a (q t )dq

(5)

Then we combine the equations of the deformable model and the cables:

       Adv -b -hH T a λ = 0 subject to : λ ≥ 0 λ ⊥ δ a + H a dq -δ * a ≤ 0 (6)
The symbol ⊥ states that each component of the vector λ that represents the forces in the cables is not null if and only if the cable is tight (δ a +H a dq-δ * a (t) = 0). On the other hand, when a cable is not tight (δ a + H a dq -δ * a (t) < 0 ), the cable cannot apply any forces λ = 0. To solve these equations, we first factorize the sparse matrix A using a LDL T decomposition. Then we decompose dv in two elements dv = dv free + dv a where dv free = A -1 b is the free motion (by fixing λ = 0 which represent a case where the cables are not pulling) and dv a = hA -1 H T a λ is the contribution of the actuator. We integrate the velocity to obtain a free displacement dq free .

Then, we project the mechanics in the constraint space using the Schur complement of the equation ( 6):

0 ≥ h 2 H a A -1 H T a Waa λ + δ a + H a dq free -δ * a ⊥ λ ≥ 0 (7)
We can compute the value of λ using a Linear Complementarity Problem (LCP) solver. Once λ is computed, we can find the final value of dv = dv free + hA -1 H T a λ and equation ( 6) is solved for the current time step. This approach, based on complementarity constraints, leads to realistic simulation of the cable/structure interaction. We could also used complementarity constraints to simulate dry friction but the level of friction was quite low on this robot and difficult to parameterize. Consequently, we suppose that we have frictionless interaction between the cable and the structure of the robot.

Inverse model

In our previous work [START_REF] Largilliere | Real-time control of soft-robots using asynchronous finite element modeling[END_REF]] and [START_REF] Coevoet | Software toolkit for modeling, simulation, and control of soft robots[END_REF]], we have proposed to use convex optimization to find the inverse model of a soft robot modeled and simulated with finite element method. For the first time, we apply the approach for deformable models with frame nodes (i.e. nodes which position are defined in SE(3), whereas previous approach were dealing only with node positions defined in R 3

For this robot, we have 10 actuators (9 cables and the translation of the base) and for the most simple case, we want to pilot only 5 degrees of freedom at the tip (3 translation and 2 rotation as the robot is not able to rotate by torsion along its main axis). For some inspection tasks, we may want to have more control on the posture of the robot, and add 2 additional directions on a node placed at the end of the second section, like explained in Figure 16. In both cases, the robot is redundant: there is more actuators than directions to be piloted in the task space. In some configurations, these redundancies decrease with cables at full-travel position and which can only act unilaterally.

Let X d i be a desired position of a DOF frame i of the model, defined in the task space. This frame will be called effector. From the current position of the same frame in the FEM simulation X i , we deduce a twist t e i ∈ se(3) that represents the gap to the desired position. This twist is defined in the world coordinate. We suppose that the gap is not large and can be corrected in one step. In case of large distance between X d i and X i , t e i could represent a fraction of the gap and the desired position will be reached in several steps.

By selecting the directions that we want to pilot in the twists (translations along X, Y , Z and/or rotations along Y , Z), and gathering the effector frames we want to pilot, we define a gap of the effectors δ e , that we want to minimize during the optimization. For instance for the configuration described in Figure 16, we have:

δ e =       I 3×3 [ 0 3×3 ] 0 1×3 0 γ 0 0 1×3 0 0 γ 0 1×3 [ 0 1×3 ] 0 1×3 [ 0 1×3 ]       t e 3 +       0 3×3 [ 0 3×3 ] 0 1×3 [ 0 1×3 ] 0 1×3 [ 0 1×3 ] 0 1×3 0 γ 0 0 1×3 0 0 γ       t e 2
where δ e ∈ R 7 represents the gap in the 7 directions that are defined in the task space, t e 3 and t e 2 stands for the twists defined the terminal points of section 3 and 2. The parameter γ is used to balance the errors in translation (defined in m) and in rotation (defined in rad).

The value of the twist t e i linearly depends on dq. So that at each step we can have the following linearized relationship:

δ e = H e dq + δ 0 e (q t ) = h 2 H e A -1 H a Wea λ + H e dq free + δ 0 e δ free e (8)
Where H e and H a are constraint projection matrices from DOF space to effector space and from actuator space to DOF space.

The optimisation goal is to find the value of λ that minimizes the norm 1 2 δ e δ T e . This optimisation is done by solving a QP (Quadratic Programming), instead of using the LCP presented in equation ( 7):

       min 1 2 hW ea λ + δ free e 2
subject to : λ ≥ 0 hW aa λ + H a dq free ≤ 0 (9)

The QP provides the value of λ and the following process is the same than for the direct case to obtain the final value of the motion dv. As the number of actuators is greater than the size of the task space (i.e. the number of effector directions), the matrix of the QP is not positive-definite and the solution could be non-unique. To make the QP convex, we introduce in the minimization expression an additional energy term E = h 2 λ T W aa λ which represents the mechanical work of the actuator forces. We then replace the equation ( 9) by:

min 1 2 δ e δ T e 2 + eE (10)
with e chosen sufficiently small so that the deformation energy does not disrupt the quality of the effectors positioning. Note that the QP problem is based on a local model of the FEM thus to be able to achieve large deformation, the linearization point (and the associated matrices) are updated at each simulation step. Moreover, since the model used by the QP is only valid locally, we bound the variation of the actuation (-δλ ≤ λ -λ prev ≤ δλ).

Control

Description

The 9 tendons of the robot are actuated by 9 DC motors controlled in position with a sampling rate of 10ms. These inner loop controllers ensure displacements with an accuracy under the millimeter range. A full stroke displacement can be achieved in approximately 1 second. The sliding platform motion is achieved using stepper motors controlled in position.

The orientation and the position of both the tip of the robot and the end of the second section are sensed using a magnetic-based sensor (Polhemus Liberty https://polhemus.com/motion-tracking/all-trackers/liberty). The resolution of this sensor is around 0.1mm and 0.01deg. This technology does not require a direct line of sight of the sensors (compared to vision based sensors) which allows the robot to inspect enclosed spaces as long as no magnetic disturbances are present in the working area (large iron objects, magnets, coils ...). These measurements can be updated with a rate up to 240Hz.

A desktop computer runs the main controller which is based on the FEM forward and inverse models. It communicates with the actuators and sensors via USB. The main controller updates the required tendon lengths and sliding position with a rate around 30Hz-40Hz depending on the simulation state.

Open loop control

A first way to control this robot is to consider the use of the inverse model as in [Duriez (2013)]. The references r i (t) are directly sent to the QP solver which computes the required actuation λ. During the QP solving step, we take into account the limitations of the actuators (both strokes [δ a , δ a ] and speed limits [dδ a , dδ a ]), the solution is guaranteed to stay in the work space of the robot and in the domain of validity of the QP problem by adding the following constraints:

δ a ≤ δ a (q t+h ) ≤ δ a dδ a ≤ δ a (q t+h ) -δ a (q t ) ≤ dδ a (11)
Where δ a (q t+h ) = hW aa λ + δ a + H a dq free is the new cable lengths to apply and δ a (q t ) is the previous cable lengths applied.

Some experiments are performed in order to evaluate the quality of the inverse model. To do that, the objective of the QP problem is set to minimize the error between the references and the position of the tip of the robot as well as its orientation. Two test are performed: The first test consists in following step trajectories along the Y axis and remain still in the others; The second test consists in following a circle trajectory in the YZ plane. Figure 9 depicts the evolution of the position of the tip of the robot during the Y step response. The settling time is around 1s and the maximal steady state error around 50mm on the Y axis. Some oscillations and asymmetry in the behavior can be noticed which illustrate the non linear behavior of the system. Deflection is present on the X and Z axis which shows that the inverse model is not perfect due to unavoidable modeling inaccuracies. Even in the nominal case, i.e. without environment interactions like contacts, a robust output feedback controller is required to track the references more accurately.

Closed loop control

In order to improve the accuracy of the positioning of the robot, a closed loop control strategy is implemented. This strategy relies on the measurement coming from the Polhemus sensor and the use of the inverse model. The main idea is to decouple the static input-output behavior of the robot by using the inverse model obtained in the previous section. Similar control schemes are used in [START_REF] Morales-Bieze | FEM-based kinematics and closed-loop control of soft, continuum manipulators soft Robotics[END_REF], Morales-Bieze (2017), [START_REF] Largilliere | Real-time control of soft-robots using asynchronous finite element modeling[END_REF]] in which a QP-based inverse model and pure integral controllers are considered to control soft robots. The results presented in [START_REF] Morales-Bieze | FEM-based kinematics and closed-loop control of soft, continuum manipulators soft Robotics[END_REF]] proved that, in the quasi-static case, this kind of control law have robustness properties towards the model inversion error. However these results are based on small speed assumptions which reduces the performances that can be obtained from this robot. Based on this knowledge and without any assumption on the speed magnitude during QP solving step, we implemented PI controllers in order to speedup the settling time while preserving some quality in the transient (limit the oscillations) and ensuring disturbance rejection.

PIs QP Sensors Robot u δ a r ε q y Figure 11: Control architecture. The block PIs represent a collection of PI controllers (one for each output), the QP block represent the inverse kinematic optimization under constraints. u is the vector of control signals, λ the actuation signals (tendon displacement), q the vector of coordinates of the robot, y are the measurements and r the references The control strategy is the described by Figure 11. The outputs of the robot q i (t) (generalized coordinates) are measured by the magnetic sensor. These measurements y i (t) are compared to the references and send to the PI controllers (one for each coordinate). The control signals u i (t) computed by the PI controllers are send to the QP solver which computes the required actuation δ a (tendon lengths and base position). Note that the number of actuators must be greater or equal than the number of outputs otherwise some PI controller may never converge and create unstable behavior.

u i (t) = K P i ε i (t) + K I i ε i (t) ε i (t) = r i (t) -y i (t) (12)
Thanks to the decoupling provided by the inverse model, each PI controller controls a single-input/single-output system. Thus they can be tuned empirically or by identification-based method around a given set-point. In this paper, the second method is employed by considering second order transfer functions with input delay during the identification stage. Also, all these PI controllers have clamping-based anti-windup scheme implemented to avoid unstable behavior around the limit of the working space of the robot.

Experiments

The efficiency of the closed loop is compared to the open loop strategy through 2 scenarios. These two cases consider the position and the orientation of the tip of the robot as the variables to be controlled. The first experiment consists in following a step reference signal on the Y axis (see also Extension 2). Figure 12 depicts the time response of the 3 controlled variables and shows how the control input u(t) provided to the QP solver compensates the external disturbances and modeling error. The results show a steady state error close to the sensor sensitivity (few mm). Figure 13 and Table 1 provide a comparison of the performances between the open loop control scheme and the closed loop. It appears that the closed loop improves the accuracy and the settling time. The closed loop also reduces the uncompensated coupling effect which is the most noticeable on the z axis (see Figure 13). The second case consists in following a circle trajectory in the YZ plane, as shown in Extension 3. The frequency of the trajectory is set to 0.05Hz. Figures 14 and 15 depicts the evolution of the position of the tip of the robot with respect to time and in the YZ plane respectively. Table 2 shows the performance indexes of the open loop and the closed loop in terms of mean error and standard deviation of the error signal. As in the step case, the closed loop improves the accuracy and is able to track the sinusoidal signal with less deflection than the open loop.

Axis

A final experiment is performed with a more practical goal. The robot task is to inspect through holes of a cube (see Extension 1). To do so, the position of the second section and the orientation of sections 2 and 3 are considered as the variables to be controlled. This setup is showed in Figure 16. This choice allows to reach a given position for the second section and orientate the camera placed at the tip of the robot inside the cube. This experiment shows that the controlled variables can be changed according to the robot task, only the PI control gains need to be re-tuned. The direct analogy that exist between the human arm and rigid manipulators provides an easy and obvious method of human operation for this type of robots.

In contrast, human operation of continuum robots remains a non-intuitive task due to the need to coordinate their (infinite) degrees of freedom and their coupled actuation [START_REF] Csencsits | User interfaces for continuum robot arms[END_REF]], [McMahan W and Walker (2009)]). In order to provide a method of tele-operation for the manipulator presented in this paper, a Touch Haptic device (3D Systems https://www.3dsystems.com/hapticsdevices/touch) is linked to the degrees of freedom of the end-effector of the robot in the simulation. By manipulating the stylus of the haptic device, the user can control the position and orientation of the end-effector of the robot. Addition- A basic solution to provide a haptic rendering is to apply a force on the stylus that is proportional to the difference between the reference position provided by the user and the measured position of the robot. However this strategy creates unwanted drag forces during fast movement due to the response time of the robot and its actuators. To overcome this problem, a filtering strategy is proposed based on the knowledge of the closed loop behavior of the robot. Instead of considering directly the reference position, the force feedback is computed from the predicted position if the robot was free to move. A gain G allows to tune the strength of the haptic feedback. To check the validity of this approach, we consider the following cases : slow and fast motions (resp. < 5mm/s and > 10mm/s) in which the user should feel nearly zero force, contact with a wall in which the user should feel when the contact occurs and finally moving an object (see Extension 4) in which the user should feel a friction-like resistance. The results of these experiments are depicted in figures 18 and 19. During the experiments we can notice that the operator can feel the wall contact along the Z axis at t = 92s. When the operator moves away from the wall there is a stick effect which is typical when a PI controller is considered (the integrator tries to reject the wall). When the operator tries to move an object (t = 138s), he can feel an opposing force as long as the object is still as the force applied is not enough to make the object move. When the force applied by the robot on the object is high enough, the object starts to move and the object friction forces suddenly reduce, which generate an accompanying haptic force due to inertia of the robot. Finally, during the free movement phases (t < 90s and 105s < t < 138s), the operator can feel some small parasitic forces which are filtered by our model based approach strategy.

Figure 16: The inspection task. The inspection entry points and the degrees of freedom selected as variables to be controlled. For this particular task, the position of the terminal point of sections 2 and the orientation of both section 2 and 3 are the controlled variables.

Discussion

Note that the proposed approach aims at being generic. The use of finite element methods permit the model of parallel structures and to include new structural parts incrementally without having to retrieve the new equations or parameters of the system. To illustrate this feature, consider the robotic arm to which we have added mechanical constrains (black trusses on figures 20) between some selected ribs to cause partial curvatures along a single section and asymmetric deformation. In order to include this modification using the methodology presented, it is enough to just add the new constraints into the model description (which is 3 lines of codes in the simulation framework).

In a framework like Coserat rods or constant curvature model, this would necessitate the addition of more sections and to take into account the coupling between non adjacent sections (the parallelism) and reconsider the controller design, all this probably by hand. Moreover, an important benefit of the use of the same model all along the design (modeling, simulation and control design) is the ability to define only once these new constraints as they will propagate all these modifications automatically to the control law design. This reduces the development time and risk of error by a non negligible factor.

To illustrate the importance of taking into account theses structural modifications we made another experiment. This experiment consists in trying to reach a given position using the inverse model in open loop. During this, we measured the error of the tip position between the model (which always reach the good position) and the real robot when the trusses are not taken into account (but present in the robot) and when they are. Note that this feature could be used in the case of mechanical failure. For example if a tendon breaks and we detect it, the simulation and the QP problem could be updated to take this new situation into account.

The current implementation of the controller does not include a correction of the simulated model (used to build the QP problem). But in one of our previous work [Morales-Bieze (2017)], we proved for a similar control that this kind of control scheme is robust to modeling error (expressed as an error on the actuator space to the effector space Jacobian matrix). This previous study showed that as long as the main entries of the Jacobians (of the model and of the real robot) have the right sign, there exist a controller that guarantee the convergence of the error to 0. Given the errors observed in open loop (10% of the length of the robot) we are always in this case. If the errors were greater, we have also presented a method to refresh the simulated state of the robot based on sensor information, thanks to a second inverse model. But that was not needed in this case.

Conclusion and future work

A novel design of an extrinsically actuated deformable continuum manipulator inspired by the skeletal spine and the engineering lattice is presented. The design makes use of a hollow backbone to prevent the shearing between the sections of the robot and the self-collision of the vertebrae to limit the local deformation. The simulation scheme used to unify the different stages in the implementation of the robot is showcased. Lastly the experimental validation of the modeling and control strategies is presented, as well as a proposed haptic feedback loop to enhance the tele-operation capabilities of the controller.

We are currently designing the next generation of the manipulator. In this iteration, some improvements over the current hardware are planned: We plan to add two additional sections to the manipulator to further enhance its workspace. In order to optimize the dimensions of the back end of the robot, we are currently working on the design of a single pulley system with a tension sensor to ensure that the tendon is always under a minimum tension to avoid slacking and to improve the control of the manipulator. The magnetic sensors currently employed to measure the position and orientation of the robot proved to be very sensitive to the noise induced by the endoscopic camera. We will replace these sensors with string potentiometers and flex sensors to reconstruct the state of the robot and use the estimated states to improve the inverse model used to control the manipulator.

We plan to improve the control of the robot by accounting for the deformations of the cables in the FEM model as well as implementing a dynamic controller [START_REF] Thieffry | [END_REF]]) to reduce vibrations and allow for the imposition of the acceleration during the trajectory of the manipulator. The implementation of a dynamic controller, when based on a reduced order model, will also decrease the computational requirements of this methodology. We plan also to include an automatic control gain tuning procedure in our framework. Finally, we plan to complement the model of the manipulator by including the self collision of the vertebrae and the use of a non-linear constitutive material law to account for plasticity phenomenon. 9 Appendix A: Index to Multimedia Extensions

Figure 2 :

 2 Figure 2: The manipulator assembly. The 3 sections of the manipulator: (1) base section, (2) middle section and (3) the distal section. The red vertebrae in the arm are the section end points where the tendons are fixed.

Figure 3 :

 3 Figure 3: Actuator setup. Each tendon is attached to a carriage mounted on top of a timing belt. The belt is moved by a Pololu 3215 DC motor.

Figure 4 :

 4 Figure 4: Longitudinal cut of the backbone and dorsal vertebrae. (a) Upper tendon of base section (b) Upper tendon of middle section (c) Backbone (d) Teflon-PTFE transmission tube allocated inside the backbone.

Figure 5 :

 5 Figure 5: Frontal view of the vertebra. (a) Vertebra lateral segments (b) Backbone attachment point (c) Vertebra dorsal segment.

Figure 7 :

 7 Figure 7: Mechanical model of the robot decomposed in beam elements. (a) element of the structure of the robot (b) geometrical parameters of the vertebrae (c) curved beam, left: 4 points of the splines, right: interpolation of the rectangular section between the two extremities of the beam (d) geometrical parameters of the central tube (e) beam model with the interpolation of the section of the tube between the two extremities of the beam.

Figure 8 :

 8 Figure 8: Model of a section with 3 vertebrae by assembling of beams. (a)shows the discretization of the structure in beam elements (in yellow). Two beams are highlighted. The actuation is modeled using constraints: H T a provides the direction of the constraint and λ is the intensity of the pulling force. (b) illustrates the assembling of the two highlighted beams. The frame node X 1 is attached to a frame DOF X i . Same for a frame node X 2 rigidly attached to frame DOF X j .

Figure 9 :

 9 Figure 9: Open loop step response of the tip of the robot y(t) and the associated QP inputs u(t).

Figure 10 :

 10 Figure 10: Open loop trajectory tracking of the tip of the robot y(t) and the associated QP inputs u(t).

Figure 13 :

 13 Figure 13: Evolution of the position of the tip of the robot in open loop y o and closed loop y c

Figure 14 :

 14 Figure 14: Evolution of the position of the tip of the robot in open loop y o and closed loop y c

Figure 15 :

 15 Figure 15: Evolution of the position of the tip of the robot in open loop y o and closed loop y c

Figure 17 :

 17 Figure 17: Haptic feedback architecture is based on the closed loop model

Figure 18 :Figure 19 :

 1819 Figure 18: Evolution of the haptic device force feedback with (red) and without using the closed loop model (blue) in various situations

Table 1 :

 1 Open loop and closed loop performances in terms of mean error, error standard deviation and settling time (10%)

		control	mean error std. dev. set. time
	X	Open loop	12.6	11.47
		Closed loop	1.61	4.82
	Y	Open loop	30.44	57.1	0.96s
		Closed loop	0.62	55.35	0.53s
	Z	Open loop	80.63	14.34
		Closed loop	0.001	7.22

Table 2 :

 2 Open loop and closed loop performances en terms of mean error and error standard deviation when the robot follows a circle trajectory in the YZ plane

		800			X
	position(mm)	750						r(t) y(t) u(t)
		700
		0	5	10	15	20	25	30	35
					time(s)
		200			Y
	position(mm)	0						r(t) y(t) u(t)
		-200
		0	5	10	15	20	25	30	35
					time(s)
		200			Z
	position(mm)	0 100						r(t) y(t) u(t)
		-100
		0	5	10	15	20	25	30	35
					time(s)
	Figure 12: Evolution of the position of the tip of the robot y(t) and the associ-
	ated QP inputs u(t) in closed loop
		Axis	control	mean error std. dev.
		X	Open loop	14.9		19.72
			Closed loop	0.24		8.90
		Y	Open loop	31.31		14.72
			Closed loop	0.35		11.70
		Z	Open loop	76.59		13.06
			Closed loop	0.042		9.04

Table 3 :

 3 Table 3 and figures 20 and 21 show the obtained modelling errors. It appears that the positioning error in open loop is 4 times more important in the case where the trusses are not modeled. Open loop modelling error with and without taking into account the truss

	Axis w.o. trusses with trusses
	X	-187mm	33mm
	Y	70mm	-44mm
	Z	-18mm	-6mm
	norm	200mm	55mm

  The authors would like to thanks Fabien Verbbrughe and Gilles Marguerite from Centrale Lille for their valuable advice.

	Extension Type	Description
	1	Video Description of the manipulator
			and inspection scenario
	2	Video	Point-to-point trajectory
			tracking experiment
	3	Video	Circle trajectory
			tracking experiment
	4	Video	Teleoperation with
			haptic feedback loop
	10 Appendix B: Table of variables