Frederick Largilliere

Mario Sanz Lopez

Eulalie Coevoet

Olivier Goury

Maxime Thieffry

Zhongkai Zhang

Bruno Carrez

Erwan Douaille

Marwa Eldiwiny

Piyush Jain

Damien Marchal

Sébastien Nelissen

Félix Vanneste

Sandra Hage Chehade

Contribution to Kinematic Modeling and Control of Soft Manipulators using Computational Mechanics

Keywords: Continuum manipulators, Finite Element Method, Continuum kinematics, Continuum manipulator control, Soft robots, Continuum Mechanics, Computational Mechanics

This work provides new methods for the kinematic modeling and control of soft, continuum manipulators based on the Finite Element Method. Contrary to the case of rigid manipulators, soft and continuum manipulators generate their motion by deformation, therefore, the proposed methodology accounts for the deformation mechanics to better describe the kinematics of these type of robots. This methodology does not produce analytic solutions, instead, a numerical approximation is provided by methods derived from Computational Mechanics. The methodology is applied to a continuum manipulator, namely, the Compact Bionic Handling Assistant (CBHA). A closed-loop control scheme based on control allocation is also presented. The models and controller are validated experimentally.

iii To my parents, Jessica and Jose Luis, for teaching me that a man, with love, can fly. v All models are wrong, but some are useful.

George 1.1 General Introduction 50 years after the construction of the first continuum manipulator, researchers and engineers still consider the modeling and control of this type of robot an open problem. The development of a generic modeling approach for continuum manipulators that can capture the main behavior characteristics regardless of its design seems like a far fetched idea. From our point of view, this issue is linked to some key factors. First, the bio-inspiration behind continuum manipulators: nature offers a myriad of shapes and motion mechanisms from which we can gather inspiration for the design of soft, continuum manipulators. However, most of these principles require a great insight of the biological object of inspiration. As we close the gap between the more traditional rigid designs of manipulators currently populating the factories and the organic bodies of snakes and tentacles, the need to trespass the knowledge boundaries of robotics as a field arises.

On the other hand, the tools used to study this type of robots usually aim at the over-simplification of the problem; as for the writing of this manuscript, the most popular approach towards the modeling of continuum manipulators is based on conventional rigid robotic methodologies that consider only the structural geometry of the robot as an important feature in the modeling. In reality, the geometrical description of the manipulator is only one part of 1.4. Chapter summary 3 freedom of continuum manipulators. The imposition of constraints in the form of Lagrange Multipliers allows for the model of sensors, actuators and endeffector. In this way, the forces applied to the manipulator by the actuators are related to the end-effector and sensors spaces. Forward and inverse kinematic models are derived from the FEM model of the manipulators. To account for non-modeled nonlinear behaviors of the robot, a closed-loop control strategy based on the implementation of the simulation as a state estimator is proposed.

Chapter summary

The manuscript is organized in 4 chapters as follows:

Chapter 2 introduces the definition of the kinematic model of continuum manipulators. To this end, the challenges of modelization are identified. The current state of the art in kinematic models and the methods used are given, along with the most recent efforts towards the formalization of a design methodology and control of continuum manipulators. The proposed methodology is then positioned in the context of the literature.

Chapter 3 formally introduces the methodology. First, a brief introduction to Continuum Mechanics is given and the fundamental steps of the Finite Element method are presented. Then, the method of modeling of continuum manipulators based on Computational Mechanics is explained. The concept of constraints, fundamental in the development and application of the method, is introduced. Also, in this chapter we talk about the simulation framework used for the implementation of the method and the general steps required to obtain a robot simulation. Then, the method is applied to a continuum manipulator, namely, the Compact Bionic Handling Assistant (CBHA). The kinematic relationships between 3 different spaces are explained and the experimental validation for forward and inverse kinematics is shown.

The objective of Chapter 4 is to present a closed-loop control strategy for continuum manipulators based on the methodology given in the previous chapters. The feed-forward control strategy based on real-time simulation is showcased and, based on the experimental results, we explain the need for more accurate controllers. The closed-loop control based on PI control law and control allocation based on Quadratic Programming (QP) formulation is proposed. The implementation of the controller is validated with experiments 4 Chapter 1. Introduction using the CBHA and the robustness analysis for the close-loop system is also performed.

Finally, Chapter 5 gives the conclusions for this work and discuss the perspectives of future work, in particular, the use of the methodology towards the design of soft, continuum manipulators based on simulation and possible improvements to the model of this type of robots.

List of Publications

The results obtained during the development of this work have made the topic of the following list of publications: Journal Publication • Bieze T., Largilliere F., Merzouki R., Duriez C. (Under review). FEMbased forward and inverse kinematics of soft, continuum manipulators. Soft Robotics Journal. Submitted on June 2017.

• Coevoet E., Bieze T., Largilliere F., Zhang Z., Thieffry M., Carrez B., Marchal D., Goury O., Dequidt J., Duriez • Duriez, C., Coevoet, E., Largilliere, F., Bieze T., Zhang, Z., Sanz-Lopez, M., Dequidt, J. (2016, December). Framework for online simulation of soft robots with optimization-based inverse model. In Simulation, Modeling, and Programming for Autonomous Robots (SIMPAR), IEEE International Conference on (pp. 111-118). IEEE.

• Lakhal, O., Melingui, A., Bieze, T. M., Escande, C., Conrard, B., Merzouki, R. (2014, December). On the kinematic modeling of a class of continuum manipulators. In Robotics and Biomimetics (ROBIO), 2014 IEEE International Conference on (pp. 368-373). IEEE.

List of Publications

• Bosman, J., Bieze, T. M., Lakhal, O., Sanz, M., Merzouki, R., Duriez, C. (2015, May). Domain decomposition approach for FEM quasistatic modeling and control of Continuum Robots with rigid vertebras. In Robotics and Automation (ICRA), 2015 IEEE International Conference on (pp. 4373-4378). IEEE.

Book Chapter

Introduction

In this chapter, we will, first, give a formal definition of a soft, continuum manipulator in the context of this work. The State of the Art related to the modeling, design and control of this type of robots is provided in order to position our work in the framework of the current literature.

Continuum manipulator definition

To grasp the concept of continuum manipulator it is always helpful to visualize a conventional manipulator composed, in general, by a set of joints, rotational or spherical, serially connected by rigid links. In this class of manipulators, the links provide a kinematic relationship between the joints; in other words, by knowing the lengths of the links, one can always know the relative position of the joints with respect to a common reference frame [Siciliano 2016]. Now, let us place additional joints along the rigid links of our imaginary manipulator.

Chapter 2. State of the Art Soon enough, our manipulator will have more active1 degrees of freedom than the ones required to position an object in space. In the literature, these type of manipulators are known as hyper-redundant [Hirose 2004]. If we keep repeating the process of adding joints until their number approaches infinity and the link lengths to zero, our robot will eventually approach to what is known as a continuum manipulator.

Figure 2.1: Left) A classic industrial manipulator, (Center) A hyper-redundant robot, the meeting ground between a rigid and a continuum manipulator, (Right) A continuum manipulator

The previous definition may imply that all the degrees of freedom in a continuum manipulator are controllable, which is in general, not the case. Often bio-inspired by the morphology and functionality of biological limbs and appendices like octopus tentacles [START_REF] Zheng | [END_REF]], elephant trunks [START_REF] Neppalli | OctArm-A soft robotic manipulator[END_REF]] [McMahan 2005[START_REF] Zhao | [END_REF]], plant tendrils [START_REF] Mehling | [END_REF]] [Yamada 2014] and other mammal tongues, continuum manipulators are composed by an elongated structure with no identifiable joints, which is continuously bending via elastic deformation [START_REF][END_REF]].

Soft continuum robot applications

The characteristics of soft continuum manipulators, like their natural compliance, high power to weight ratio, and reduced dimensions compared to their rigid counterparts, make them particularly suitable for applications in which the contact with humans is unavoidable or even desired. Applications such as skeletal trauma treatment [START_REF] Wilkening | [END_REF]] [Alambeigi 2017], endoscopy [Conrad 2013[START_REF] Cianchetti | STIFF-FLOP surgical manipulator: mechanical design and experimental characterization of the single module[END_REF][START_REF] Fraś | New STIFF-FLOP module construction idea for improved actuation and sensing[END_REF]] and minimally invasive surgery [START_REF] Qu | Motion control of a bio-inspired wire-driven multi-backbone continuum minimally invasive surgical manipulator[END_REF][START_REF] Orekhov | [END_REF]] [Mahoney 2016] have proved the great potential of application of continuum robots in the medical field, as extensively reviewed in [Burgner-Kahrs 2015].

As suggested by the very first prototypes of continuum manipulators [Anderson 1967]. The slender shape and high dexterity of continuum manipulators can be exploited in tasks such as minimally invasive inspection [START_REF] Mehling | [END_REF][START_REF] Tonapi | [END_REF]] and search and rescue [Bajo 2010] [Li 2017].

Continuum robots have been studied with the goal of exploiting their locomotive capabilities [START_REF] Godage | [END_REF][START_REF] Kang | [END_REF][START_REF] Arienti | [END_REF]], although the review presented in the following of this manuscript is concerned mainly on manipulation.

Bio-inspiration

If one sees for the first time a continuum manipulator, without any previous knowledge on the concept behind it, one can immediately identify the remarkable morphological similarities that this type of devices have with some soft-bodied animals, particularly with the muscular hydrostats. Muscular hydrostats, commonly found in elephant trunks, mammal tongues and octopus tentacles are soft muscular structures that can bend, extend and twist and provide the force required for movement and skeletal support to animals (or limbs) that lack a rigid skeleton, see Fig. 2.2.

Muscular hydrostats are typically composed by a fluid-filled cavity surrounded by a muscular wall reinforced with connective tissue fibres. The arrangement of the muscle fibres in a hydrostatic muscle include both circular and longitudinal muscle fibres. These two muscle fibres can antagonize one another to produce a variety of shape changes including elongation and bending [START_REF] William | Hydrostatic skeletons and muscular hydrostats[END_REF]]. The fluid inside the cavity of a hydrostatic limb is mainly a liquid which resists to volume change, thus, to create an elongation of the limb, the circular muscles contract to decrease the diameter while increasing the length to allow for a constant volume inside the cavity. The study of the biomechanics of hydrostatic structures have shown also additional fibres with more intricate configuration patterns that allow for more complex shape changes like twisting, present in octopus tentacles and mammals and reptile tongues [START_REF] William | Tongues, tentacles and trunks: the biomechanics of movement in muscularhydrostats[END_REF]]. The complete replication of hydrostats is very complex, but the study of their underlying function principles have given roboticists an interesting insight and a solid starting point in the design of soft, continuum manipulators.

Classification

Continuum manipulators can be broadly classified with respect to the type of backbone they possess. Single backbone manipulators (see Fig. Another classification of continuum manipulators is based on the actuation system implemented to apply forces and torques to the backbone. When the actuators are embedded in the structure and apply forces directly to the backbone, the actuation scheme is called intrinsic. Most continuum robot designs with multi-backbones use pneumatic cavities to conform the structure, essentially making a backbone composed of actuators that can be reshaped by applying pressure to the cavities. Extrinsic actuation systems place the actuators outside of the structure. These apply torques and forces at localized points via a mechanical link to the backbone, which is made by an incomprehensible elastic rod. A hybrid actuation scheme replaces the central elastic rod in an extrinsic design by an actively controlled actuator, which is in many designs, a pneumatic cavity [Chirikjian 2015[START_REF] Ataollahi | [END_REF]].

State of the art in soft, continuum manipulators

This section describes the research developed in the last 20 years, towards the design, modeling and control of soft, continuum manipulators. The review is divided into each of the subtopics mentioned in order to provide a clear overview of the field.

Design of continuum manipulators

As stated in [Walker 2013a], the first prototype of a continuum robot reported in the literature was the Tensor Arm [Anderson 1967], designed by [START_REF] Victor | Tensor arm manipulator design[END_REF] Conceived to be used under water, the prototype was able to achieve a wide range of shapes; however, the relationship between the shapes and inputs was highly complex and challenging for the computational resources of that time. Based on extrinsic actuation, the robot used nylon filaments routed along the structure through spacer discs that apply torques directly to the backbone to produce bending, see Fig. 2.6. As is usual in tendon-based designs, the termination points of the cables define the bending sections. This prototype has inspired since then, a significant number of designs based on the same principle.

Figure 2.6: The "Tensor Arm" by [START_REF] Victor | Tensor arm manipulator design[END_REF]].

Walker and Hanna presented the Elephant's Trunk Robot [Walker 1999], a continuum manipulator with extrinsic actuation composed of 4 main bending sections. The robot uses a pulley system outside the structure to pull cables Chapter 2. State of the Art and bend the sections. Spacer disks are used to give the robot its characteristic shape but they also provide a mechanical coupling; out of the 32 DoF only 8 are actively controllable, leaving the rest to be coupled by passive springs positioned between adjacent segments of the manipulator at a 90 degrees angle with respect to the tendons. Thanks to this work the identification of critical hardware necessities in the field was possible, in particular the need for elastic structures. Later work in design of continuum robots replaced the spring system by an incompressible elastic rod [START_REF] Gravagne | [END_REF]]. This choice of backbone has been, since then, the most used in tendon-based designs of continuum manipulators [Hemami 1984[START_REF] Simaan | [END_REF][START_REF] Bardou | [END_REF][START_REF] Zhao | [END_REF]]. Another choice of backbone is the use of a central spring [START_REF] Mehling | [END_REF]] that provides the desired compliance to the manipulator but makes the operation of this type of manipulators a difficult task, since the forces applied by the tendons tend to be absorbed by the spring. However, spring-backbone manipulators have been implemented as actuated endoscopes with great success [Yoon 2011].

At the same time, research in the area of artificial muscle design facilitated the conception of locally actuated designs, such as the Oct-Arm [START_REF] Walker | [END_REF][START_REF] Mcmahan | [END_REF]], the "European Octopus" [Calisti 2012] [Guglielmino 2010[START_REF] Laschi | [END_REF]], and the Bionic Handling Assistant [Rolf 2012], which is inspired by the trunk of an elephant, with possible appealing advantages like full-body manipulation. Intrinsic design manipulators proved to be very dexterous by reproducing the way longitudinal muscles work in an animal trunk or tentacles [START_REF] Mcmahan | [END_REF]]. Originally developed to provide actuation to a prosthesis device for severely paralyzed hand patients, the McKibben artificial muscle [START_REF] Vernon L Nickel | Development of Useful Function in the Severely Paralyzed Hand[END_REF]] has been used in continuum robot designs with great success [START_REF] Mcmahan | [END_REF][START_REF] Pritts | [END_REF]]. This type of pneumatic muscle is composed by a soft cavity constrained in its deformation by a braided sleeve. The muscle is considered contractor if the cavity expands radially when the pressure is increased in the cavity, or extensor if the cavity expands longitudinally. This operation is analogous to the longitudinal muscles. Experiments conducted on the McKibben muscle [Klute 1999] show that it provides a first order approximation to a biological muscle when contracting.

Another venue for researchers to design extrinsic actuation is the implementation of concentric tubes as the structure of the manipulator [START_REF] Sears | [END_REF]] [Gilbert 2016]. The body of the robot is composed by a telescopic array of concentric tubes that are free to rotate and translate with respect to each other. To provide bending to the robot, the tubes are often pre-curved, so that when a distal tube is translated outside of the preceding one, it will nat- urally bend [Furusho 2006]. Pre-curved concentric tubes designs have found a niche of applications in the medical field [Wu 2017], where their slim and clean design is well suited to act as actuated endoscopes in minimally invasive surgery [Furusho 2006] [Su 2012]. Given the advantages of intrinsic actuation in terms of dexterity, particularly the ability to extend the backbone longitudinally and in turn, improve its Chapter 2. State of the Art ability to bend [START_REF] Walker | [END_REF]], robot designers started to experiment with hybrid actuation schemes that incorporate both pneumatic cavities and tendons. An example of this type of design is the Air-Octor [Jones 2004a] [Jones 2004b] depicted in Fig. 2.9, a manipulator composed by a soft pneumatic hose that acts as the backbone, which is actuated by 3 cables to achieve bending. While the prototype proved to be capable of full-body manipulation, it suffered from sagging and kinking of the soft backbone cavity [McMahan 2005], a problem that is common to this type of designs. Another example of a realized hybrid actuation scheme is the KSI manipulator [START_REF] Immega | Tentacle-like manipulators with adjustable tension lines[END_REF]], developed by Kinetic Sciences Inc. in 1994 [Immega 1995].

Backbone variable stiffness

In order to exploit to the maximum the capabilities of continuum robots, the ability to variate the stiffness of its structure is needed, particularly in medical applications were excessive contact forces at the tip of the manipulator can cause severe trauma to the surrounding tissue. In the literature, one can find different approaches to the variable stiffness of the backbone. In [START_REF] Mahvash | [END_REF]], a stiffness control approach is presented. The method is based on the derivation of the kinematic model of the manipulator as the product of two transformations: the first transformation describes the non contact kinematics of the manipulator and it is specific to the robot itself, while the second transformation calculates the tip deflection due to the applied forces and is derived using the Cosserat rod formulation [START_REF] Mahvash | [END_REF]]. To implement a desired tip stiffness, the two transformations are used to solve for the actuator positions that deform the manipulator so as to generate the required tip force at the measured tip position. This method is applied to a concentric tube design intended for minimally invasive surgery. In [START_REF] Goldman | [END_REF]] a method for the stiffness control of intrinsic actuated manipulators is presented. The proposed algorithm relies on the measurement of actuation forces, in conjunction with the stiffness of the backbone, to move in a direction to minimize the environment interaction.

The combination of pneumatic backbones and tendons provide the means to increase the apparent stiffness of the structure by antagonizing the set of redundant actuators that compose the robot [Walker 2013b]. In [Maghooa 2015] the pneumatic backbone and tendons of a hybrid manipulator are controlled simultaneously to achieve a certain position of the tip of the arm while actively changing the stiffness of the robot applying the principle of antagonistic actuation. This functionality is also bio-inspired [START_REF] Shiva | [END_REF]], being present in a significant number of animal arms, trunks and even the human tongue.

Figure 2.9: The Air-Octor hybrid manipulator composed by a pneumatic backbone which is bent by tendons, as presented in [Jones 2004b] Some continuum manipulator prototypes utilize local rigidification mechanisms to simplify the stiffness control; such is the case of the STIFF-FLOP manipulator [START_REF] Cianchetti | STIFF-FLOP surgical manipulator: mechanical design and experimental characterization of the single module[END_REF][START_REF] Fraś | New STIFF-FLOP module construction idea for improved actuation and sensing[END_REF]], that utilizes granular jamming [Steltz 2010], in conjunction with pneumatic cavities to achieve a certain degree of variable stiffness of the structure. Each section of the robot is composed by an elastomeric cylinder that hosts three pneumatic cavities to provide bending and extension as well as a central chamber that holds the granular media. The section is surrounded by a braided sleeve to prevent the radial deformation of the elastomeric cavities. The granular media utilized is coffee powder, which jams when the central cylinder is under negative pressure causing drastic stiffness changes of the structure [START_REF] Jiang | [END_REF]]. The STIFF-FLOP is one of the most advanced continuum manipulators prototypes currently under development, having force [Noh 2014] and tactile [START_REF] Sareh | [END_REF]] capabilities provided by sensors embedded in the structure of the robot.

More advanced actuation systems have been studied with the intention of stiffness control. In [Sadeghi 2012], Sadeghi et al. demonstrated a method to control the stiffness of a worm-like soft robot utilizing electro-rheological materials that can change phase, from fluid to solid and vice versa, under the presence of magnetic fields.

Chapter 2. State of the Art

Design methodologies

The key challenge for creating soft machines that achieve their full potential is the development of controllable soft bodies using materials that integrate sensors, actuators and computation, and that together enable the body to deliver the desired behaviour [START_REF] Rus | Design, fabrication and control of soft robots[END_REF]]. The conception of a design methodology for continuum robots is not a trivial task due to the wide variety of possible backbone-actuation schemes, the applications to which the design is aimed at and more importantly, due to the interdisciplinary nature of the design tasks; very often, electronics, mechanics and computer science come into play during the integration of the system. Nevertheless, some design guidelines can be found in the literature. For example, a pragmatic approach to the design of continuum arms inspired by the octopus tentacle is presented in [START_REF] Guglielmino | [END_REF]]. In [START_REF] Bedell | [END_REF]] an algorithm for the optimization of concentric tube robot design is proposed; given the description of the task space and the number of sections in the robot, the algorithm solves for a design with the desired workspace.

The fabrication of artificial pneumatic muscles is also a daunting enterprise that requires, in most cases, multiple design iterations to achieve a desired behavior. A direct2 finite element simulation has been used in [START_REF] Connolly | [END_REF]] for the programming of desired deformations of pneumatic muscles with braided sleeves.

The work of Hauser et al. in [Hauser 2011] presents the concept of morphological computation. This concept aims at the relegation of difficult control tasks in soft bodied robots to the body itself, making use of inherit characteristics of the robot, like the natural compliance of its structure [Paul 2006]. While it is difficult to directly design smart bodied robots, this work provides useful design guidelines compatible with any kind of bio-inspired robot [START_REF] Pfeifer | [END_REF]].

Modeling of continuum robots

In order to control the movement of continuum robots, kinematic models, which relate the configuration (pose) of the backbone and task space (endeffector position) and actuator variables must be established. These models are fundamental for the development of control strategies and path planning for continuum manipulator robots. Equation (2.1) depicts the general form of a kinematic model for conventional manipulators

x = f (q) (2.1)
where x is the position of the end-effector inside the task space of the robot which is, more often than not, expressed in cartesian coordinates, q is the set of configuration variables related to the joint variables (which are, in general, observable and directly controlled) and f is the unknown function that relates task space variables to the configuration of the robot. For soft, continuum robots however, f depends on the configuration variables and also on the mechanics of the deformable material. We can then rewrite Eq. 2.1 for continuum manipulators as

x = f (q, ϑ) (2.2)
where ϑ represents the mechanics of the soft material. The modeling task deals with finding this relationship and the assumptions used in order to accomplish this goal dictates the type or flavor of the modeling approach.

In general, kinematic modeling approaches of continuum manipulators can be broadly classified as quantitative models, which deal with the description of the robots in a mathematical fashion making use of geometry and elasticity theories, and qualitative models, which are numerical abstractions of the more complex model and often make use of experimental data to find the closest solution to the kinematic problem.

Qualitative models Due to the complexity involved in the modeling of continuum robots, qualitative approaches, which aim at learning the kinematic problem by the use of learning algorithms, have a significant level of popularity due to their ability to by-pass the modeling task. These approaches are based on previous knowledge of input-output data derived from experimentation and can provide accurate fast approximations to the kinematic solutions. Rolf and Stein proposed in [START_REF][END_REF]] a control scheme for the Bionic Handling Assistant based on goal babbling [START_REF][END_REF]]. The control scheme is capable of dealing with highly non-linear behavior like hyperelasticity, plasticity and non-stationarities, which are very complex to address analytically.

Control schemes based on Neural Networks have been also investigated. [START_REF] Giorelli | [END_REF]] used a feedforward neural network to learn the in-Chapter 2. State of the Art verse kinematics of a tendon driven manipulator. In this approach, a geometrical model of the manipulator is used for data generation and in [START_REF] Giorelli | [END_REF]] a controller based on the static model Jacobian is used to control a nonconstant curvature manipulator. Another feedforward neural network controller was proposed by Braganza et al. [START_REF] Braganza | [END_REF][START_REF] Braganza | [END_REF]]. In this work, no model is required; instead, a low level controller uses the neural network component to compensate for the dynamic uncertainties of the manipulator. More recently, the work of Melingui et al. [START_REF] Melingui | [END_REF]] employed a neural network to solve the inverse kinematics based on the measurements of the end-effector position. In [START_REF] Melingui | [END_REF]] an adaptive algorithm is implemented to improve the performance of the controller by allowing rapid position convergence of the end-effector. While qualitative controllers have shown great success in real world scenarios, the fact that the learning base changes when the operation conditions of the robot change make data-driven models limited in their application. This is the main reason as to why quantitative approaches still dominate the literature.

Quantitative models Current quantitative models in the literature can be classified into 2 different main approaches: those that describe the behavior of the continuum backbone curve with respect to geometry parameters, which we will call Geometric Models, and those that use classical elasticity theories, such as strings or rods, to describe the manipulator with respect to the way its structure behaves under external loading, which we will call Elasticity Models. In this section, a review of these modeling approaches is given in order to highlight the fundamental differences between these models and the approach presented in this work and to position clearly the contribution of this work in the current literature.

Geometric models

We refer to the first type of quantitative models as Geometric models. These approaches compute the kinematic relationships making use of the geometry of the manipulators, in particular, the curvature that the robot exhibit under the actuation forces; however, the forces themselves are not accounted for in the model.

The first kinematic analysis of continuum robot backbones was presented by Chirikjian [START_REF] Chirikjian | Theory and applications of hyperredundant robotic mechanisms[END_REF]] [Chirikjian 1994b]. In this approach, also called Non-Constant curvature approach, the Cartesian position of backbone curve points x(s, t) is given by the integral representation:

x(s, t) = � s 0 l(σ, t)u(σ, t)dσ (2.3)
where u(s, t) = [0 1 0] T is the unit vector tangent to the curve at s and l(s, t) is a scaling factor that controls the length of the curve tangent according to local backbone extension or contraction. In [START_REF] Chirikjian | Theory and applications of hyperredundant robotic mechanisms[END_REF]], the follolwing backbone curve representation is given:

x(s, t) = [x 1 (s, t) x 2 (s, t) x 3 (s, t)] T =   � s 0 l(σ, t) sin(K(σ, t)) cos(T (σ, t))dσ � s 0 l(σ, t) cos(K(σ, t)) cos(T (σ, t))dσ � s 0 l(σ, t) sin(T (σ, t))   (2.4)
Figure 2.10: Backbone curve parametrization with K and T being the curve parameters as illustrated in Fig. 2.10. The backbone reference frame Φ(s, t) attached to the curve at s and which origin coincides with x(s, t) has its orientation described by

Φ(s, t) =   cos(K(s, t)) sin(K(s, t)) cos(T (s, t)) -sin(K(s, t)) sin(T (s, t)) -sin(K(s, t)) cos(K(s, t)) cos(T (s, t)) -cos(K(s, t)) sin(T (s, t)) 0 sin(T (s, t)) cos(T (s, t))   (2.5)
The inverse kinematics of the continuous curve, that is, the problem of finding the backbone position and orientation given a specific task constraint, Chapter 2. State of the Art is then solved by finding the values of the independent functions l(s, t), T (s, t) and K(s, t). A modal approach for the selection of these values is developed in [START_REF] Chirikjian | Theory and applications of hyperredundant robotic mechanisms[END_REF]]. The disadvantage of this approach is that the set of backbone shapes available in the model is restricted by the combination of modal functions, and further analysis of the hardware is required to tune the selection of the basis functions to be able to model arbitrary backbone shapes. In [START_REF] Mochiyama | [END_REF]], Mochiyama and Suzuki used the same approach to derive the kinematics and dynamics of a cable-like flexible manipulator.

In [START_REF] Hannan | [END_REF]], Hannan and Walker presented what would later become the most used approach towards the kinematic modeling of continuum manipulators. This work makes the assumption that, after actuation, the shape of the backbone of the manipulator can be assimilated as a piecewise constant curvature curve and its evolution, from base to end, is then described by a set of 3 discrete transformations for the planar curve. The first transformation is a rotation to point the tangent direction at the curve beginning to the end of the curve, followed by a translation along the new direction, from the curve starting point to the curve end point, and finally, a second rotation of same magnitude as the first to realign with the tangent at the curve end point. See Fig. 2.11

Figure 2.11: The geometry of a planar constant curvature curve segment Thus, a virtual rigid link manipulator with two rotational and one prismatic joint can be used to derive the kinematic model of the planar constant curvature backbone. The corresponding kinematic model is found using standard Denavit-Hartenberg notation [START_REF] Hartenberg | [END_REF]] for the virtual rigid manipulator. The Denavit-Hartenberg parameters are shown in the following table:

Link θ d a α 1 * 0 0 -90 • 2 0 * 0 90 • 3 * 0 0 0 •
where θ and d are the joint parameters for rotational and prismatic joints, respectively, α is the link twist and a is the length of the common normal vector between consecutive joints. Using the previous parameters, the homogeneous transformation matrix, which relates the Cartesian position of the end of the curve with respect to the base is formed as:

A 3 0 =     cos(θ 1 + θ 3 ) -sin(θ 1 + θ 3 ) 0 -d 2 sin(θ 1 ) sin(θ 1 + θ 3 ) cos(θ 1 + θ 3 ) 0 d 2 cos(θ 1 ) 0 0 1 0 0 0 0 1     (2.6)
Further manipulation of the derived homogeneous transformation matrix is required to relate the joint variables of the virtual robot to the configuration space variables of the continuous curve. Noting that

θ 1 = θ 3 = ks 2 (2.7)
and also

d 2 =|| x(s) ||= 2 sin θ k (2.8)
with k being the inverse of the radius of the arc and s being the arc length, the previous transformation matrix turns into

A 3 0 =     cos(sk) -sin(sk) 0 1 k (cos(sk) -1) sin(sk) cos(sk) 0 1 k sin(sk) 0 0 1 0 0 0 0 1     (2.9)
which relates the continuum curve arc length and curvature to the task space. It is important to note that the arc length s can be chosen arbitrarily to model the transformation from the curve shape to any point along the backbone curve.

Spatial continuum kinematics can also be modelled by adding an additional pair of rotational joints to the virtual manipulator that rotates the curvature plane as shown in Fig. 2.12

Link θ d a α 1 * 0 0 90 • 2 * 0 0 -90 • 3 0 * 0 90 • 4 * 0 0 -90 •
As reviewed in [Walker 2013a], the constant curvature approach produces the same kinematic relationships as the non-constant curvature approach. This approach has been used subsequently. Mahl et al. implemented a constant curvature continuum kinematics as fast approximate model for the Bionic Handling Assistant [Mahl 2014] and Escande et al. used the same approach to model the smaller version of the robot, the Compact Bionic Handling Assistant [START_REF] Escande | [END_REF][START_REF] Escande | [END_REF]].

Continuum manipulators composed by several serially connected sections can be modelled using the same approach by concatenating the transformation matrices of each section, as presented in [START_REF] Jones | [END_REF]] and [START_REF] Bardou | [END_REF]].

As is common when using geometric approaches, further modeling is required to find the relationships between the curve parameters and the actuator variables (displacement of tendons or pressure inside pneumatic muscles are the most common). This transformation is specific to each manipulator hardware and therefore, cannot be generalized. Nevertheless, geometric approaches provide a fast and accurate approximation to the kinematics of continuum manipulators in scenarios where gravity can be countered by a clever orientation of the manipulator or when the weight of the robot is small compared to its structural stiffness. More important, they provide close form kinematic equations which can be easily implemented in control algorithms.

Elasticity models

The second class of quantitative kinematic models are, what we will call, Elasticity based models. These models substitute the backbone of the robot by continuum mechanic objects, like rods and strings, to better capture the elastic behavior of the backbone. In contrast with geometric models, mechanic-based models provide a description of the deformation of the backbone based on the constitutive law of the material from which the backbone is made, by establishing static equilibrium between external loading, such as gravity, and internal forces.

The Cosserat rod theory [Antman 1973] has become a great resource in the research of kinematics for continuum manipulators. In Cosserat rod theory, a homogeneous transformation matrix g(s) is used to describe the rod with s ∈ [0 l] being the reference parameter. The evolution of g(s) along s is defined by

Ṙ(s) = R(s)û(s), ṗ(s) = R(s)v(s) (2.10)
where the dot denotes derivative with respect to s, R and p are the rotation matrix and position vector of g at s and R(s) and u(s) are the kinematic variables that represent the linear and angular rates of change of g(s). The ˆoperator is a linear mapping that converts � to the correspondent skew symmetric matrix as

û =   0 -u z u y u z 0 -u x -u y u x 0   (2.11)
Given an undeformed reference configuration of the rod g * (s), the reference kinematic variables v * and u * can be obtained by

[v * u * ] T = (g ( * -1) )(s) ġ * (s)) � (2.12)
where the operator � represents the inverse operation of ˆ. The internal force and moment vectors (in global coordinates) are denoted by n and m, the applied force distribution per unit of s is f , and the applied moment distribution per unit of s is l. The classical form of equilibrium differential equations for a Cosserat rod are

Chapter 2. State of the Art

ṅ(s) + f (s) = 0 (2.13) ṁ(s) + ṗ(s) × n(s) + l(s) = 0 (2.14)
Using the constitutive law of the rod to map the kinematic variables to the internal forces and moments we have

n(s) = R(s)K se (s)(v(s) -v * (s)) (2.15) m(s) = R(s)K bt (s)(u(s) -u * (s)) (2.16)
where K se is the stiffness matrix for shear and extension and K bt is the stiffness matrix for bending and torsion. Given the assumption that the stiffness matrices are constant with respect to s, a full set of explicit model equations that provide the values of ṗ(s), Ṙ(s), v(s) and u(s) can be derived. Usually, a linear constitutive law is used, but the approach can work under nonlinear constitutive laws as well.

The work of [START_REF] Davis | [END_REF][START_REF] Davis | A model for the embedded tendon control of a slender three-dimensional flexible robot link[END_REF]] on the modeling of flexible robot links with tendon control introduced this method in the context of robotics, although it was Trivedi et al. in 2008[Trivedi 2008a] who conducted the first work directly in the field of soft manipulators using this approach, as reviewed in [Burgner-Kahrs 2015].

One of the advantages of mechanic-based modeling methods is the great variety of elements from which one can conform the model. As an example, one of the special cases of Cosserat rod, the Kirchhoff rod, has been used to derive models for concentric tube robots [Rucker 2010b]. The Kirchhoff theory assumes inextensibility and neglects transverse shear strain, which are generally regarded as good assumptions for long thin rods like the tubes that make up active cannulas [Rucker 2010a]. Wenlong et al. used Timoshenko beam theory, which accounts for shear deformation and beam twisting, to map a driven load to the pose of a continuum manipulator [START_REF] Wenlong | [END_REF]]. Euler-Bernoulli beam theory, a special case of Timoshenko beam theory that only considers lateral external loading, has being used to simplify the computation of the mechanics of concentric tubes as well [START_REF] Webster | [END_REF]].

While elasticity models can be fast to compute (up to a certain number of elements), they can be quite complex to parametrize and implement, also, as explained before, since close form solutions for continuum mechanic objects exist only for a very limited number of cases, numerical approximations of the solutions of the partial differential equations that arise in their development must be performed.

Dynamics and control of continuum robots

The use of geometric models to describe the pose of a continuum manipulator is hardly sufficient due to the fact that internal and external loading in the manipulators play an important role in the configuration of the robot. Even in free space, gravity can cause significant deflections in the curvature of the robot even when no payload is being carried by the manipulator. Dynamic models, which aim at the description between the backbone shape and the forces that cause the deformation of it, have been the objective of intense research in the field since the very beginning [Chirikjian 1993] [Chirikjian 1994a].

As in the case with the modeling of kinematics, researchers have tried to apply methodologies commonly used for rigid manipulators to describe the dynamics of continuum ones. In [Mochiyama 2002], Mochiyama and Suzuki used the Lagrangian approach to describe the dynamics of a string-like manipulator. In this development, a discretization stage in the modeling is performed in which the backbone is assumed to be formed by circular slices located at σ along the backbone. Each slice has a mass. After calculating the kinetic and potential energy for each slice, the total energies are found by integration along the backbone and substituted in the Lagrangian L = K -P . Once the Lagrangian is calculated, the dynamic model is found as

∂ ∂t � ∂L ∂ θi (σ, t) � - ∂L ∂θ i (σ, t) = τ i (σ, t) (2.17)
where (θ i , τ i ) are the configuration variables and the forces that change them, respectively. In the initial study [Mochiyama 2002], non-extensible backbones were considered. Later work by Tatlicioglu [START_REF] Tatlicioglu | [END_REF]] considered extensible backbones for a planar manipulator; however, the complexity of the closed form models generated renders the application to non-planar backbones too computationally expensive to be implemented in the control of manipulators.

More computationally efficient formulations of the dynamics of continuum backbones have been developed. In [Kang 2011], the Newton-Euler formula-Chapter 2. State of the Art tion is applied to derive the dynamics of a manipulator for underwater applications. In [Giri 2011], Giri and Walker presented an approximation to the dynamics of a continuum arm using lumped parameters (mass, spring, and structural damping) and the Lagrangian formulation. This approach was implemented to compute the dynamics of an octopus inspired arm in [START_REF] Zheng | [END_REF]] and in [START_REF] Falkenhahn | [END_REF]] to describe the dynamics of the Bionic Handling Assistant manipulator.

The elasticity models discussed in section 2.3.2.2 provide a formulation that explicitly takes the internal and external forces into account. However, for more complex applications in which the external loading of the manipulator comes from the contact with its environment, or the velocity of the manipulator induces inertial effects in its behavior, these models are not sufficient. In [START_REF] Yu | [END_REF]] inverse dynamics and sliding mode control schemes for a continuum manipulator are proposed based on the Euler-Lagrange formulation and in [START_REF] Falkenhahn | [END_REF]] the dynamic control of the Bionic Handling Assistant based on a lumped mass model is presented.

Given the complexity involved in the derivation of the dynamics of continuum manipulators, model-less approaches have been proposed. Ivanescu et al [START_REF] Ivanescu | [END_REF]] avoid the difficulties of solving the complex PDE system that describes the dynamics of continuum manipulators by proposing controllers based on an energy formulation. A model-less controller based on the estimation of the Jacobian under spatial constraints is presented in [START_REF] Yip | [END_REF]] and [Yip 2016]. In this regard, controllers based on fuzzy-logic [START_REF] Qi | [END_REF][START_REF] Qi | [END_REF]] and adaptive algorithms [START_REF] Melingui | [END_REF]] have been presented.

In medical applications, the use of in-vivo feedback information gathered by radio-graphic and electromagnetic images [START_REF] Arai | [END_REF][START_REF] Bertocchi | [END_REF]] have been used to control continuum manipulators with great success. Camera vision systems [Chitrakaran 2007[START_REF] Boudjabi | [END_REF]] have also been employed to provide feedback information for the controller. However, all these techniques rely on external sensors to close the loop. When the environment in which the manipulator works makes impossible the use of external sensors, controllers based on embedded sensors in closed-loop [Penning 2011[START_REF] Penning | [END_REF][START_REF] Bajo | [END_REF]] can be used instead.

Work contextualization and contributions 29

Work contextualization and contributions

Given the classification of the approaches towards the modeling and control of continuum manipulators presented in section 2.3.2, it is now possible to position the methodology presented in this work, in the context of the current literature. Fig. 2.13 provides a helpful visualization of the position of our method with respect of other previously proposed methods and their classification.

Figure 2.13: The context of the proposed method in the current state of the art As it is illustrated in the figure, the proposed method is based on the mechanics of the materials, which means that no geometrical assumptions are done in order for the method to be applicable. Moreover, being a mechanicsbased approach, the description provided by the method explicitly accounts for the internal and external loading of the continuum backbone, which means that the deformation caused by gravity or the load of a payload will be reflected in the model. Nevertheless, there are a few important differences between the method here discussed and other mechanic-based methods: First, while most methods are based on Cosserat rod and other elasticity theories, the proposed approach is based on Finite Element Analysis which provides geometrically accurate descriptions of the backbone. Second, the modeling Chapter 2. State of the Art of actuators, sensors and end-effector is done implicitly by the formulation of the constraints, so no further modeling is required.

In the following, the contributions to the field derived from the development of this Ph.D. work are enlisted:

• FEM-based quasi-static modeling approach that accounts for complex structural shapes of single and multi-backbones.

• The model implementation of two types of actuation systems (i.e. tendons and pneumatic muscles) currently implemented in the majority of designs of soft manipulators.

• The integration of sensors in the model that allows for a description of the manipulator in the configuration space.

• A closed-loop control strategy for continuum manipulators based on the state estimation of the robot.

In the next chapter, the assumptions considered during modeling and the study framework of this methodology (namely Continuum Mechanics and Finite Element Method) are introduced.

Chapter 3

FEM-based model of Continuum Manipulators

Introduction

In contrast with rigid robots, soft manipulator kinematics not only depend on the geometry of the robot, but also on the mechanical properties of its material, in particular the stiffness, that dictates the net configuration of the robot under a set of external forces. This statement represents the core of the modeling principle considered in this manuscript. As an example, consider the manipulator in Fig. 3.1 (a).

A tendon routed through the structure is pulled to achieve a desired endeffector position. In Fig. 3.1 (b), a manipulator with the same geometry but with inhomogeneous material is actuated in the same manner by the same amount, however, the resulting configuration differs from (a). This is because of the influence of the material stiffness in the kinematics of the manipulator. While rigid manipulator kinematics can be used to solve positioning problems with the assumption of resistance/counter-actuation to gravity or load effects, soft manipulators easily comply to these forces and deform. To answer the same problems of positioning, it is then necessary to take into account the current deformation (i.e. change of geometry) induced by these forces to obtain a kinematic relation between the po-sition of the end-effector and the position of the actuators.

Another difference between soft continuum manipulators and rigid ones is the impact that external forces (particularly gravity) have on the final pose of the manipulator [Trivedi 2008b]. As mentioned before, in a continuum manipulator, there can be only a finite number of actuated degrees of freedom1 . The state of the remainder of the (infinite) degrees of freedom in a continuum robot backbone will be determined by both the constraints of the controlled degrees of freedom and internal and external forces [Walker 2013a].

In this work, this problem is addressed by a discretization of the degrees of freedom of the continuum manipulator, through numerical methods provided by Computational Mechanics. In the following, the Continuum Mechanics framework is briefly introduced to help us describe the motion of the continuum manipulator based on its deformation and to help the reader to understand the assumptions made in this modeling methodology. This work also explains the methodology of kinematics estimation on a continuum manipulator based on FEM.

Continuum mechanics framework

From a Continuum Mechanics perspective, we will consider the structure of a continuum manipulator at a macroscopic level, that is, there is no space between the particles that compose the manipulator, the body of the manipulator fills the entire space it occupies. In this section, we present the basic notions of Continuum Mechanics with the example of a continuum manipulator. The interested reader may refer to the Appendix A.1 for a more formal introduction.

Consider the manipulator in Fig. 3.2. with all its material particles occupying a 3-dimensional domain with a starting configuration k 0 . Given a set of external loads, the manipulator will deform, changing its geometry. The new deformed configuration of the manipulator is called k i .

A material particle of the manipulator will have a position P in k 0 and a position P in k i . Then, we can always find a function M that maps P from P. In the study of deformation, two important physical quantities arise: the stress and the strain.

Stress

The deformation of the manipulator will cause internal body reaction forces. These forces are called stress which is defined as the force per unit area.

Strain The displacement of the particles in a continuum manipulator can be considered, in the general case, as a result of a rigid body motion and a deformation. A measure of strain is the computation of the deformation part of the displacement.

Constitutive material law

As stated before, the internal forces in a continuum manipulator are caused by its deformation. Therefore, it is of interest to find a relation between the stress and the strain. The stress-strain relationship is given by the Constitutive law and it depends on the material from which the manipulator is build. We assume that the manipulator is built from linear elastic materials. For linear elastic materials, the constitutive law is represented by Hook's law.

In the experimental characterization of linear elastic materials, Hook's law is expressed in terms of the Young's modulus, which is a measure of the linear stiffness, the Poisson's ratio, which expresses a measure of the change in volume of the material, and the shear modulus. We will also make the assumption that the material of the manipulator is isotropic, which means that the material properties are independent from the direction they are measured. It means that we do not consider material made of fibers; however, this assumption does not prevent the use of different materials in the structure of the robot. In this case, only one Young's modulus, one Poisson's ratio, and one shear modulus characterize the material properties.

Forces in the continuum manipulator

In order to compute the response of the continuum manipulator to a set of known forces, the Principle of Conservation of Linear Momenta, commonly known as Newton' s second law of motion is applied. This principle states that the rate of linear momentum of a collection of particles is always proportional to the forces acting on the particles.

The forces acting on the particles are the sum of internal and external forces. Internal forces are of two types, volumetric and boundary forces. The boundary forces are related to the stress by a tensor called Cauchy's stress tensor and they need to be integrated along the entire manipulator. This integration produces a differential equation describing the motion of the continuum manipulator. In general, this equation cannot be solved analytically because it involves an infinite number of degrees of freedom. Instead, numerical methods provided by the field of Computational Mechanics are used to approximate the solution. In the following, we will introduce the Finite Element Method as means of discretization of continuum manipulators.

Chapter 3. FEM-based model of Continuum Manipulators

Finite Element Method

In order to model the kinematics of continuum manipulators, the equation of motion provided by Continuum Mechanics needs to be solved. This equation however, cannot be solved analytically, hence the need of numerical methods to approximate its solution. One of these numerical methods is the Finite Element Method (FEM). For a more formal introduction to FEM, the interested reader can refer to the Appendix A.2

The FEM can be conceptually divided in 3 main steps: discretization, discrete solution and assembly. The main idea behind the FEM is the consideration of the continuum body, in this case the continuum manipulator, as being comprised of a set of subdomains, called finite elements.

In this step, the theoretical infinite number of degrees of freedom in a continuum manipulator is reduced to a finite number, which makes the approximation of the solution attainable. The collection of finite elements that compose the structure under study is called mesh, as shown in Fig. 3.3. The generation of the mesh is an important process in FEM, since it defines the fidelity of the model. As a general rule, the model fidelity increases as we increase the number of elements in the mesh. In practice, the mesh is generated from a surface model of the manipulator, like those provided by Computer Assisted Design (CAD) software.

Surface model Volumetric mesh Many choices exist for selecting the type of element that will compose the mesh. One of the considerations to make in this regard is the dimensions of the problem. For example, for 1-dimensional problems, the typical choice is the two-node truss; for 2-dimensional problems plate or shell elements are often selected, and for 3-dimensional problems the hexahedron or tetrahedron element are a good example of volumetric elements (see Fig. 3.4).

The type of element selected for the mesh will determine the nodal interpolation, which is basically the expression of any physical quantity inside the mesh depending only on the values at the nodes. The accuracy of the interpolation depends on the amount of nodes an element has and on the choice of interpolation or shape functions. These functions can be linear or higher order polynomials. Once the mesh is generated, the solution to the differential equations is approximated over each element.

For the solution of the problem at the element level, each element is disassembled from the mesh and is referred to a reference element in which the computations are made easier. Because of this process, the nodes in each element will have a local numeration (at the element level) and a global numeration (at the mesh level). The final step in the FEM deals with the re-assembly of the elements using its global numeration.

To find the global equation system we must assemble all the element equations. In other words, we must combine local element equations for all the elements used in the discretization. Before the solution of the global system, boundary conditions should be imposed. Boundary conditions in the most simple case, occur at the fixation points of the manipulator (where the manipulator is attached to its base). In this case, the imposition of this boundary in the global system is simple; the equations for the nodes that are fixed to the base are eliminated from the global system. Another, more complex boundary condition is the one that represents the actuators, end-effector and sensors in our soft, continuum manipulator. In our work, these conditions are imposed as constraints in the model and allow for the kinematic descriptions of the manipulator. In the next part of this chapter, the Chapter 3. FEM-based model of Continuum Manipulators derivation of the kinematic models for continuum manipulators based on the aforementioned constraints will be explained.

FEM-based kinematics of soft manipulators

The computation of the kinematic model of continuum manipulators developed in this section depend on the simulation of the FEM model explained in 3.3. To arrive at the models of continuum manipulators, we shall start from the formulation given by Newton's second law that models the dynamic behavior of the robot as

M(x)ẍ = F ext -F(x, ẋ) + H T λ (3.1)
where x ∈ R n is the vector of generalized degrees of freedom (in this case, the position of the nodes), M(x) : R n � → M n×n is the inertia matrix of the manipulator, ẋ ∈ R n is the velocity vector, ẍ ∈ R n is the acceleration vector, F ext represents the external forces (like gravity), F(x, ẋ) represents the non-linear internal forces that depends on the generalized state of the manipulator and H T λ ∈ R n gathers the contribution of the constraints forces(actuators, sensors and end-effector) as Lagrange multipliers.

In this work, the study is limited to quasi-static behavior on purpose. Thus, in a first approach, the assumption is that the control of the robot is performed at low velocities. As such, the term M(x)ẍ in Eq. 3.1 that corresponds to the inertia effects is removed. On this basis, Eq. 3.1 for the quasi-static analysis of the continuum manipulator becomes

F ext -F(x) + H T λ = 0 (3.2)
The internal forces F(x), which are obtained from the FEM of the manipulator are non-linear. We use the Taylor series expansion to compute a unique linearization of F(x) per simulation step. During each step i of the simulation, the internal forces are updated as

F(x i ) ≈ F(x i-1 ) + K(x i-1 )dx (3.3)
where K(x) is the tangent stiffness matrix that depends on the actual positions of the nodes and dx is the difference between positions dx = x ix i-1 . This linearization is valid as long as the displacement of the nodes dx is small.

As mentioned before, the lines and columns that correspond to fixed nodes are removed from the system to get a full rank for matrix K. In F and K, the rows (and columns for K) contain the components of the internal forces (x, y, z) for the nodes, in the order corresponding to their global node numbering in the mesh.

The way H is obtained is explained in the following, but it is important to mention that its computation is performed with the values obtained from the previous simulation step. We then use the expression H(x i-1 ) and through the linearization explained in Eq. 3.3, we obtain the following formulation :

-K(x i-1 )dx = F ext + F(x i-1 ) + H(x i-1 ) T λ (3.4)
The variables dx and λ are both unknown and are found during the optimization process. It should also be noticed that the matrix K is highly sparse. In the implementation, a conjugate gradient solver is used and preconditioned by a sparse LDL T decomposition. For a mesh composed of about 1000 nodes and about 3000 tetrahedral elements, a refresh rate of 60Hz is obtained with the implementation available in the simulation framework.

Constraint for the end-effector

To set the Lagrange multiplier on the end-effector, a point or a set of points of the robot need to be considered as the end-effector; these are the points we seek to control in position. It could be any points mapped on the FE mesh. For each point, the constraint objective is to reduce the difference between the end-effector position and its desired position p des . Thus, a function δ e (x) : R 3n → R 3 with n being the number of nodes, evaluates this difference along x, y and z (see Fig. 3.5). If the end-effector corresponds to a node i of the mesh, the function is δ e (x) = x ip des , where x i is the position of node i. If the effector is set inside an element, we use

δ e (x) = n � i=0 φ i (p ef f )x i (3.5)
where p ef f is the position of the end-effector in the rest configuration of the FEM model and φ i is the shape (interpolation) function associated to node i. If several points are used for the end-effector position, the vector δ e (x) gathers the evaluation of the difference for all the points. The function is then R 3n → R 3m , where m is the number of end-effector points.

The matrix H used for the end-effectors corresponds to H e (x) = ∂δe(x) ∂x . The matrix H e is highly sparse: A row, that corresponds to a component of a point of the end-effector, will contain non null values on a very small number of columns. As the point is mapped on a single tetrahedral element, there is a maximum of 4 non-null values per row. Of course, the column should match with the components of the nodes, given the fact that the non-null values are gathered in 3x3 diagonal block matrices.

Finally, an important point is the effort value that is put on the Lagrange multiplier that corresponds to the terminal effector. The value of λ e will depend on the load applied on the end-effector. Two cases can be considered: I If the points defined as end-effector move freely in the space, there is no physical interaction, so the contribution of the constraint vanishes λ e = 0.

II If one or several points of the end-effector carries one object l which mass creates a load that could deform the structure. In such cases, the corresponding load should be set on λ e = m l g with m l being the mass of the object and g the gravity field.

In both cases, it is assumed that the value of the Lagrange multiplier λ e is known, whether by experimentation or identification.

Actuator constraint model

One can find different actuation schemes for soft, continuum manipulators in the literature. From piezo electric actuation [Su 2012], electro-rheological fluids [Sadeghi 2012], to electro-active polymers [START_REF] Chikhaoui | [END_REF]]; however, as implied in 2.3.1, the two most prominent actuation schemes for continuum robots remain the pneumatic actuators, and cables or tendons. As such, these two types of actuation are considered in this work.

The model of the actuators considers their physical characteristics. The contributions of these actuator constraints are unknown before the optimization process. However, given the type of actuation, the constraint is not set in the same way.

Cable actuator

In a first case (Fig 3.6), a cable is used to actuate the manipulator. The cable can simply be attached at one point of the structure, but it can also go through several other points or guides; in that case, frictionless guides are considered. In all cases, the unknown λ a is the stretching force inside the cable. It is evident that this force is unilateral (λ a ≥ 0) since we cannot push with the cable.

Let's suppose now that the points are numbered starting from the extremity where the cable is being pulled. The matrix H is computed this way: At each point p, representing one of the cable guides (Fig. 3.6), we take the direction of the cable before p and after p

d b = x p -x p-1 �x p -x p-1 � (3.6)
d a = x p+1 -x p �x p+1 -x p � (3.7)
To obtain the constraint direction that is applied to the point, we use

d p = d a -d b (3.8)
Note that the direction of the final point is equal to the direction "before" as d a does not exist. These constraint directions are mapped on the nodes using the interpolation

   . . . F n . . .    =    . . . φ n (α, β, γ) d p . . .    λ a = H T a λ a (3.9)
A function δ a (x) is defined to provide the length of the cable, given the position of the constrained nodes. This constraint formulation allows also to consider the physical limitations of the actuators, which in this case is the cable stroke. The actuator stroke can also be included by imposing

δ a (x) ∈ [δ min δ max ] (3.10)
Through the use of this function, we get

H a = ∂δ a (x) ∂x (3.11)

Pneumatic actuator

The formulation is compatible with pressure-based actuation of cavities that are placed in the structure. In that case, the effort λ a is the uniform pressure exerted on the wall of the cavity, while δ a represents the volume.

As the pressure is uniform inside the cavity, a single constraint can be set for each pneumatic actuator; however, all triangles of the cavity wall will be involved.

For each triangle t, the area a t and the normal direction n t are computed.

If this result is multiplied by the pressure, we obtain the force applied by the pneumatic actuator on the triangle t. We distribute this contribution to each of its nodes by dividing the resulting vector by 3. We sum up the results of each triangle in the corresponding column of H T a . This gives us the corresponding force acting on each of the nodes of all the triangles that compose the cavity

F i = � t∈S,i∈t a t 3 n t λ a = (H T a ) i λ a (3.12)
where F i is the pressure force assigned to the node i and S is the set of triangles composing the cavity. Usually, pneumatic actuators only provide positive pressure so λ a ≥ 0. However, in some cases, it is also possible to create both negative and positive pressure using vaccum/pressure actuation. In that case, there is no particular constraint on the unknown value of λ a , despite an eventual limit (max / min) of pressure that can be achieved by the actuator.

Sensor constraint model

In order to relate the end-effector position and the geometry of the manipulator, one needs sensors that can measure the geometrical state or shape of the robot. In the literature, this task is done either by using external sensors, such as stereo-vision systems and x-ray cameras [Otake 2014], or by sensors embedded in the body of the manipulator.

In this study, we assume that the robot has embedded sensors. Various studies have been conducted to evaluate the feasibility of different types of embedded sensors for shape reconstruction, including electromagnetic and piezoelectric polymers [Cianchetti 2012], cable length measurement [START_REF] Segreti | [END_REF][START_REF] Murphy | [END_REF]] and Fiber Bragg grating (FBG) [Farvardin 2016[START_REF] Searle | [END_REF]]. We can summarize the common aim of these techniques as the estimation of the configuration of the manipulator by measuring the relative position of a set of points along a string.

In this study, we implement the sensors in the simulation, using the same string concept. As in the case of the cable actuator, the string of the sensor is routed through several friction-less guides, at n points x n , as shown in Fig. 3.8. In the model, the measure of the lengths read by the sensor will be

n-1 � i=1 �x i+1 -x i � (3.13)
sensor string sensor guides which evaluates the distance between each sensor guide after the position of the nodes has been updated. A function δ s is defined to represent the current length of the sensor string given the position of the constrained nodes. The matrix H s that gathers the directions of the sensor constraint is obtained in the same way as for the cable actuator. Since the sensor does not apply any forces to the manipulator, λ s is equal to 0.

Reduced model in the constraint space

The classical resolution of a deformable body motion using FEM (like solving the static equilibrium of the structure described at equation 3.4) provides a forward model: it allows to compute the displacements of the structure, given the values of the efforts put on the actuator λ a . However, in the case of position control, the actuation λ a is the unknown. Consequently, for controlling the motion of the manipulator, an inverse model is needed, which is challenging to compute in real-time as the size of the system is in the range of several thousands degrees of freedom. In this work, another approach is used, based on the projection of the mechanics in the constraint space that drastically reduces the size of the optimization problem. This approach, initially developed in [Duriez 2013], is generalized. A new formulation of the inverse problem in the form of a quadratic programming (QP) optimization (developed in [Largilliere 2015]) is used.

Reduced compliance on the constraint space

As stated above, the optimization process relies on a projection of the mechanics in the constraint space. Each constraint has a direction that is set by a line of the matrices H e , H a and H s . This matrix is usually sparse, as the direction of the constraints is mapped on few nodes of the FE mesh. The values of the effort applied by the actuators λ a are not known at the beginning of the optimization process, whereas the value of λ e is supposed to be known and the value of λ s is 0.

The first step consists of obtaining a free configuration x free of the robot which is found by solving the equation 3.4 while considering that there is no actuation applied to the deformable structure. In practice, the known value of λ e is used and λ a = 0 is imposed.

The linear equation 3.4 is solved using a LDL T factorization of the matrix K. Given this new free position x free for all the nodes of the mesh, one can evaluate the values of From the FEM formulation of the problem that uses the global matrix K (Eq. A.55 ), a formulation that accounts for the directions of the constraints placed for actuators and end-effectors is derived. Using the Schur complement of matrix K in the Lagrange multiplier-augmented system [START_REF] Przemieniecki | [END_REF]], a small formulation of δ e is obtained. This variable expresses the difference between the desired position for the endeffector and its current position in terms of the actuators contributions λ a :

δ e = � H e K -1 H T a � � �� � Wea λ a + δ free e (3.16)
The Schur complement also provides similar formulations for the difference between a desired sensor or actuator position and its current position:

δ a = � H a K -1 H T a � � �� � Waa λ a + δ free a (3.17) δ s = � H s K -1 H T a � � �� � Wsa λ a + δ free s (3.18)
This step is central in the method. It consists of projecting the mechanics into the constraint space. As the constraints are the inputs (effector position shift and sensor length shift) and outputs (effort to apply on the actuators) of the inverse problem, the smallest possible projection space for the inverse problem is obtained. It allows for a projection that drastically reduces the size of the search space without loss of information. Indeed, the following subsection shows how the matrices W ea and W aa provide the mechanical coupling equations between actuators and effector point(s).

After this projection, the optimization is processed in the reduced constraint space to get the values of λ a . This part is described in 3.5.3.

The final configuration of the soft robot, at the end of the time step, is obtained as

x = x free + K -1 H T a λ a (3.19)
It should be emphasized that one of the main difficulties in this implementation is the computation W ea and W aa in a fast manner. No pre-computation is possible as their value changes at each iteration. However, this type of projection problem is frequent when solving friction contact on deformable objects. Thus, several strategies are already implemented in SOFA.

Coupled Kinematic Equations

Using the compliance operator W ea , one can get a measure of the mechanical coupling between effector and actuator, and with W aa , the coupling between actuators. For instance, the displacement δ i e created on the end-effector (along a direction stored on the line i of matrix H e ) by a unitary force λ j a applied by the actuator (which is stored at the line j of matrix H a ) is directly obtained by

Δδ i e = w ij ea λ j a + δ i,free e (3.20)
As the motion is created by deformation, the motion of the actuator j is influenced by actuator k. Through the same principle, actuator k also influences the displacement of the end-effector. To get a kinematic link between actuators and effector, the method needs to account for the mechanical coupling that can exist between actuators. It is captured by W aa that can be inverted if actuators are defined on independent degrees of freedom. Thus one can get a kinematic link by rewriting equation (3.17 ) is composed of a reduced number of linear equations that relate the displacement of the actuators to the displacement of the effector. Consequently, matrix W ea W -1 aa is equivalent to the Jacobian matrix of a rigid manipulator. This matrix is a local linearization provided by the FEM model on a given position and it needs to be updated for deformations with large displacements.

In conclusion, with basic manipulation of the matrices W ea and W aa , a condensed direct model of the soft manipulator can be obtained.

Inverse kinematic model solution by convex optimization

The goal of the optimization is to find how to actuate the structure so that the end-effector of the robot reaches a desired position. This was initially proposed in [Largilliere 2015]. It consists of reducing the norm of δ e which actually measures the shift between the end-effector and its desired position. Thus, computing min( 1 2 δ T e δ e ) leads to a Quadratic Programming (QP) problem

min � 1 2 λ a T W T ea W ea λ a + λ a T W T ea δ free e � (3.22)
subject to (course of actuators) : δ min ≤ δ a = W aa λ a + δ free a ≤ δ max and (case of unilateral effort actuation) :

λ a ≥ 0 (3.23)
The use of a minimization allows to find a solution even when the desired position is out of the workspace of the robot. In such a case, the algorithm will find the point that minimizes the distance to the desired position while staying in the limits introduced by the course of the actuators.

In practice, the QP solver available in the Computational Geometry Algorithms Library (CGAL) [Fabri 1998] is used. The matrix of the QP W T ea W ea is symmetric.

If the number of actuators is equal or inferior to the size of the endeffector space, the matrix W T ea W ea is also definite. In such a case, the solution of the minimization is unique. In the case when the number of actuators is greater than the degrees of freedom of the effector points, the matrix of the QP is only semi-definite. Consequently, the solution could be non-unique.

A new criterion for the minimization can be introduced, based on the deformation energy. Indeed, this energy E def is linked to the mechanical work of the forces exerted by the actuators. E def can be computed by evaluating the dot product between λ a and the displacements of the actuators

Δδ a = δ a -δ free a (3.24)
due to the actuator forces

E def = λ a T Δδ a = λ a T W aa λ a (3.25)
Matrix W aa is definite positive if the actuators are placed on different nodes of the FEM or with different directions (i.e. if there is no linear dependencies between lines of H a . Thus, one can add this energy in the minimization process by replacing (3.22) with:

min � 1 2 λ a T (W T ea W ea + νW aa )λ a + λ a T W T ea δ free e � (3.26)
with ν chosen sufficiently small so that the deformation energy does not disrupt the quality of the effector positioning. In practice,

ν = tr(W T ea W ea ) tr(W aa ) * 10 -3 (3.27)
is chosen so that the term νW aa does not alter the value of the QP matrix. Thanks to this modification, the QP matrix becomes definite positive and a unique solution of the problem can be found.

Method implementation

In this section, the inverse FEM method is implemented in order to obtain forward and inverse kinematic relationships for a geometrically com-Chapter 3. FEM-based model of Continuum Manipulators plex continuum manipulator. Two different descriptions are obtained for the inverse kinematics; one that relates the position of the end-effector to the actuator space, and another that relates the position of the endeffector to the geometric configuration of the robot. The kinematic models are used to pilot the robot in feedforward configuration. The forward kinematic model obtained from this approach is compared in terms of position accuracy to two different geometric models, developed for the same robot. Finally, the results are discussed.

Simulation framework

SOFA [START_REF] Faure | [END_REF]] is a simulation framework first released in 2007, originally intended for the interactive computation of biomechanical models for soft tissue. The framework is an open source C++ library with a modular architecture in which the simulations are built as scenes. The scenes are composed by various components that encapsulate one of the aspects of a simulation, such as the mechanical objects representing the degrees of freedom, the forces and constraints, the differential equations, the main loop algorithms, the linear solvers, the collision detection algorithms or the interaction devices.

At the present time, an auxiliary module, or plug-in for SOFA, dedicated entirely towards the modeling and control of soft robots and its environment is in development. This plug-in allows for a rapid generation of mechanic-based models of soft devices and their interaction with the environment. The algorithm proposed by Duriez in [Duriez 2013] is at the core of this plug-in and allows for the direct control of soft robots based on the simulation of the inverse FEM of the robots.

In order to solve non-linear deformations, while providing a computationally economic method that is able to simulate the robots in real time, a corotational formulation of volume FEM, like the one developed in [START_REF] Felippa | [END_REF]], [Felippa 2000] is implemented in SOFA.

Corotational FEM

In physical simulation of deformable bodies with small strains, the nonlinear part of the deformation of the body is in most cases neglected. However, when the elements undergo a rigid transformation, this nonlinear deformation is no longer small compared with the linear part and can no longer be neglected. Felippa and Haugen [START_REF] Felippa | [END_REF] proposed the computation of the rigid transformation by adding an additional frame of reference, called corotated configuration, attached to the elements that moves with it like a "shadow". Element deformations are measured with respect to the corotated configuration. The element stress field and tangent stiffness matrix are corrected using the element rotations. A significant feature of this implementation is that it allows for large rigid motions, as long as the strains are small. The implementation is element-independent and compatible with linear finite elements, which makes it particularly suitable for real-time FEM simulation.

Mesh generation

The mesh generation is a critical part in the simulation process. In our case, volumetric meshes (volume elements) are required. The Computational Geometry Algorithms Library (CGAL) is a software written in C++ that allows for the generation of surface and volume meshes. The input to the CGAL software are surface objects in stereolithography format (.stl), like those generated by a Computer Assisted Design Software (CAD). In practice, we use only one type of element for a given mesh; however, in some cases, different mesh subdomains are required, for example, to represent heterogeneous materials. CGAL is able to generate meshes composed of different subdomains while respecting the physical interfaces; elements do not cross interfaces.

Description of the Compact Bionic Handling Assistant

The CBHA is the bionic continuum manipulator component of the RobotinoXT, a didactic mobile platform designed by Festo Robotics [Festo 2012]. The system is shown in Fig. 3.9 (a). The bionic continuum manipulator is formed by 2 serially connected sections of pneumatic actuators, an axially rotating wrist and a compliant gripper. Without actuation, the manipulator has a length of 206mm, with each section having a length of 103mm. The width at the base of the manipulator is 100mm long and the top has 80mm of width. In our study, the end of the second section is considered as the end-effector.

Each section of the manipulator is composed of an array of pneumatic actuators or bellows, connected in parallel configuration, as shown in The pressure inside the cavities is regulated by a piezo proportional valve terminal, which includes its own pressure regulators. The air is compressed by two small membrane pumps that can attain a pressure of 2.5 to 2.7 bars. However, the nominal operation pressure for a single cavity in practice is in the range of 0.3 to 1 bar. The pose of the manipulator is obtained as the contribution of the poses of the 2 sections.

In order to sense the state of the robot, string potentiometers measure the lengths of the actuators. The control unit in the RobotinoXT is a EA09 control board that handles the sensor readings and drive units. An Application Programming Interface (API) embedded in the system provides communication with external applications via ethernet connection The manipulator is built by a process of Additive Manufacturing (AM), which is the name to describe the technologies that build 3D objects in a layer-upon-layer fashion. The material employed in the CBHA construction is Polyamide-Nylon, which has a Young's Modulus in the range of 1.6 to 2.3 GPa and a Poisson's ratio of 0.39.

Simulation of the CBHA

In order to simulate the CBHA and kinematically control it, a mesh, composed by 3528 tetrahedra is created from the surface representation of its structure (Fig. 3.11). To allow for a fast computation of such a big model, a model reduction technique, based, on a domain decomposition is performed (see Appendix B). We select the end-effector point to correspond to the tip of the second section and each pneumatic actuator is modeled following the procedure described in 3.4.2.

The implementation of the simulation that allows for the interactive control of the robot is handled by an external computer featuring an Intel Xeon CPU W3690 at 3.47GHz x6 with a performance, in terms of simulation steps, of 18.4 frames per second. The communication between the simulation and the robot is done via W-LAN with a latency of 3ms, and when actuation commands are sent, the pressure regulators in the CBHA have a response time of 10ms.

As implied in 3.4.1, one of the features that sets apart this modeling approach is the ability to predict yielding of the structure under external loading given an accurate parametrization of the material mechanics. This feature is of special importance when the robot is large in comparison to its thickness or when the orientation of the structure as a whole cannot be changed to counteract the effects of gravity. When the loads are known a priori, the magnitude of the force acting on the end-effector λ e is used in Eq. 3.4 to compute the configuration of the manipulator that accounts for said force.

As a preliminary experiment to assess the predictive capabilities of the model, we observe the deflection of the robot when a set of known loads are applied at the end-effector position. First, an initial configuration of the manipulator without loading is selected and the position of the endeffector is measured with respect to a reference frame; then, the load is applied and the new position of the end-effector is recovered. The loads are applied statically, so only the stationary state of the manipulator is observed. The same load is applied to the model of the manipulator using the same initial pose and the resulting end-effector position is also recovered. The model of the sensors presented in 3.4.3 is used to attain the initial configuration of the robot in the simulation. This configuration is measured by the string potentiometers.

A vector that connects the initial and final end-effector positions represents the deflection caused by the load. In order to assess the repeatability of the measurements, the loading sequence described is performed 40 times for each loading value and the average value is then used for the model validation. The results are presented in Fig. 3.12 Figure 3.12: Comparison between measured and predicted deflections caused by external loading on the CBHA manipulator

The results illustrate the compliance to loading profile of the manipulator, with a maximum deflection error of 4.107mm with an average error of 2.104mm between the measured values (blue line) and the prediction of the model (red line). The results also show that the CBHA presents strain hardening/necking stages of plastic behavior which corresponds to the compliance of the plastic material used to make the manipulator. This non-linear behavior is not considered in our model and therefore, the model prediction is accurate only for a small region of the profile. This behavior can be accounted for, in order to improve the predictive capabilities of the model, by implementing two constitutive laws with a switching condition, depending on the region of compliance in which the manipulator is operating. However, this would modify significantly the way the inverse FEM is formulated.

Kinematic models

In this section, the methodology developed in this chapter is applied to obtain the forward and inverse kinematic relationships of the CBHA based on the FEM simulation. An experiment is performed to assess the accuracy of the forward model using the real robot. The results are then compared to those obtained using 2 different modeling approaches. The inverse kinematic model, derived from the inverse FEM simulation of the robot, is also shown and validated experimentally.

Forward kinematic models

The forward kinematic model of a soft manipulator deals with the problem of finding the end-effector position, given a defined configuration of the manipulator. For a rigid manipulator, this configuration is simply the set of variables associated with the joints of the robot. In contrast with the rigid robots, the variables that express the configuration of a soft manipulator change with respect to the morphology of the robot and its type of actuation, and therefore, cannot be obtained in a straight forward manner. The FEM-based methodology presented provides a formulation to obtain the kinematic relation between the end-effector and the configuration of the manipulator.

Given the intrinsic nature of the CBHA, the configuration of the robot is represented by the lengths of the sensor strings, that correspond to a position of the end-effector, as illustrated in Fig. 3.13.

string lengths (model input)

e ector position (model output) The input to our model becomes the configuration of the robot, and the output will be the actuation required to achieve said configuration. When the actuation is computed and applied to the FEM model, the position of the end-effector that corresponds to the given configuration is retrieved. Eq. 3.18, which is reproduced here for clarification, is used to relate the configuration of the manipulator to the actuation space.

The lengths are recovered by the string potentiometers of the CBHA.

δ s = � H s K -1 H T a � � �� � Wsa λ a + δ free s (3.28)
In this approach, no geometrical assumptions are needed. Each part of the manipulator is modeled in detail as described in 3.6.5. Once the constraints have been incorporated in the model, the convex optimization finds each actuator contribution required to achieve the desired sensor lengths. The position of the end-effector is not represented explicitly in Eq. 3.18. Instead, the corresponding position is recovered after the position of the nodes is updated.

In the following, we briefly present 2 different forward models of the CBHA previously introduced in the literature, in order to compare the results obtained experimentally.

Constant curvature model

This approach follows the development shown in 2.3.2.1. It was tailored to the case of the CBHA in [START_REF] Escande | [END_REF]] and experimentally validated in [START_REF] Escande | [END_REF]]. In this approach, the evolution from end-to-end of a section i is described, in terms of backbone parameters, by 2 coupled rotations and one translation, given by the homogeneous transformation

i j T =     cos 2 φ i cos θ i + sin 2 φ i cos φ i sin φ i (cos θ i -1) cos φ i sin θ i x i cos φ i sin φ i (cos θ i -1) sin 2 φ i cos θ i + cos 2 φ i sin φ i sin θ i y i -cos φ i sin θ i -sin φ i sin θ i cos θ i z i 0 0 0 1     (3.29)
The end-effector cartesian coordinates of the bending section i are given by (x i , y i , z i ), where The backbone variables, namely the rotation angles φ i and θ i and the curvature radius r i , are illustrated in Fig. 3.14.

x i = r i cos φ i (1 -cos θ i ) (3.30) y i = r i sin φ i (1 -cos θ i ) (3.31) z i = r i sin φ i (3.32)
In order to have the kinematic model relating the configuration of the manipulator to the end-effector, the backbone variables need to be expressed in terms of the actuator lengths. These relationships are

φ i = tan -1 √ 3(l 3 -l 1 ) 2l 1 -l 2 -l 3 (3.33) i = D i 3d i (3.34) r i = (l 1 + l 2 + l 3 )d i D i (3.35)
where

D i = 2 � l 2 1 + l 2 2 + l 2 3 -l 1 l 2 -l 1 l 3 -l 2 l 3 (3.36)
The lengths of each actuator in section i are represented by l 1 , l 2 and l 3 . The parameter d i represents the diameter of section i. In this model, each section is considered to be a cylinder with constant radius. In the Hybrid approach to the modeling of the CBHA [START_REF] Lakhal | [END_REF]], the structure of the robot is considered as 17 vertebrae serially connected.

Between each pair of vertebrae, an inter-vertebra section is modeled as a 3UPS-1UP joint. The behavior of a sub-structure composed by 2 vertebrae and the inter-vertebra is represented by a parallel robot with 3 DoF, as depicted in Fig. 3.15.

Figure 3.15: Sub-structure of the CBHA modeled as a paralel robot

The position and orientation of the upper vertebra k + 1 with respect to the lower vertebra k is given by the transformation matrix

k k+1 T =     cos θ k sin θ k sin Ψ k sin θ k cos Ψ k 0 0 cos Ψ k -sin Ψ k 0 -sin θ k sin Ψ k cos θ k cos θ k cos Ψ k z k 0 0 0 1     (3.37)
where the angles θ k and Ψ k represent the pitch and roll angles of the vertebra, respectively. The jaw of vertebra k + 1 with respect to the vertebra k is constrained by a passive prismatic joint. The position of the center of vertebra k + 1 is given by the longitudinal translation z k .

The prismatic variable q n,k , shown in Fig. 3.15, represents the length of the inter-vertebra, which is a percentage of the total length of the Chapter 3. FEM-based model of Continuum Manipulators pneumatic actuator. The value of q n,k is computed by considering the minimum and maximum elongation of each inter-vertebra.

Forward kinematics experimental validation

In order to validate the position tracking capabilities of the model, a trajectory inside the task space of the robot is defined. The trajectory is composed of a set of via points for which the position can be measured in the real robot. The cartesian position of the end-effector is measured using a stereo-vision system composed by a set of 6 infrared cameras. For each position, the corresponding set of actuator lengths was recorded using the string potentiometers. The set of lengths was then used as an input for the forward kinematic model in the simulation.

After imposing the sensor constraints, the configuration of the robot is updated and the position of the end-effector is then recovered. In this experiment we assume that there is no payload being carried by the end-effector. Fig. 3.16 shows the error between the measured trajectory and the trajectory estimated by the FEM simulation. Fig. 3.17 and 3.18 show the estimated trajectory for the constant curvature and hybrid approaches, respectively. The euclidean error comparison between the 3 approaches is summarized in Table 3.1, where RMSD stands for the root-mean-square deviation, ME is the maximum error, CC is the constant curvature model, HA is the hybrid model and FEM is the finite element-based model. From Fig. 3.18, we can observe that the Hybrid approach is the least accurate. This can result from the fact that one cannot precisely measure the displacement of each vertebra (only 6 string potentiometers are available). Moreover, this model was initially developed to be able to inverse it, more than for the attainable precision of the forward kinematic model. The constant curvature model has a standard deviation of 12.87mm, which represents 6.43% of the total length of the manipulator; however, for longer manipulators with bigger variations in their diameter with respect to their length, the error increases considerably.

Moreover, as soon as a deformation in the manipulator is induced by an external force, both of these approaches lose relevance.

Inverse kinematic model

Inverse kinematics deal with the problem of finding the correct configuration of the robot, given a desired position of the end-effector. This configuration, as for the forward kinematic model, is represented by the set of sensor lengths that describe the state of the manipulator. Of course, the description of the robot can be given in the actuator space directly [START_REF] Bosman | [END_REF]], using in this case the relationship between the end-effector position and the actuators contribution. This provides a pressure-to-position model that requires a precise control over the actuation (pressure inside the pneumatic actuators). Due to the embedded control of the CBHA, the precise control cannot be attained however.

Instead, the relationship given by Eq. 3.16 is used to find the contribution of the actuators required to achieve the desired end-effector position. After these forces are computed and applied to the model in the simulation, the reading from the sensor model given in Eq. 3.13 will be taken as the solution of the inverse kinematic problem. Fig. 3.19 illustrates the input and output relationships in this model.

To validate the method, a set of 50 end-effector positions are selected inside the task space of the robot and the corresponding set of lengths for each position is recorded by the sensors of the robot. The same set of positions is used as inputs for the inverse model simulation, and the resultant length of each sensor is computed after the model is updated. This study is summarized in Table 3.2, where l 1 , . . . , l 6 represent the lengths of the sensors and their values are in mm, µ represents the mean The results show a mean error between 2.43mm and 4.08mm across all lengths, which represents between 1.21% and 2.04% of the total length of the manipulator. Again, the set of actuator contributions (in this case the pressures applied to the cavities) obtained from the optimization process can be used as input for the real robot to obtain a feed-forward control signal, as will be explained in the next chapter.

Conclusion of the chapter

The results of the experimental comparison show that the FEM-based forward kinematic model has greater accuracy in terms of position tracking, compared to the geometric models. The model of the sensors in the Figure 3.20: Comparison of the 6 measured and estimated lengths for a predefined set of end-effector positions for the CBHA manipulator allows for a representation of the kinematics in the configuration space that is less sensitive to model uncertainties, compared to a description in the actuator space. Nevertheless, given an accurate model of the constitutive law of the material, the method can provide a description in the actuator space that can be used to control the robot directly. Unlike other methods, the kinematic model involving the actuator space can be derived without any additional transformations.

The FEM model has a few limitations in its development, however. These limitations represent the main source of error in the experimental results: for the moment, the constitutive law used to model the material of the trunk is only an approximation, as shown in 3.6.5. Moreover, nonlinear effects like the plasticity of the material are not yet implemented in the model. When the trunk is bent at a maximum angle, the outer walls of the pneumatic cavities collide with each other, as shown in Fig. 3.21. The consideration of these collisions is not yet implemented in the simulation; however, the formulation of the contact forces is already developed and implemented in SOFA.

Additionally, the extension of the CBHA is limited by a steal cable that runs through its backbone. This constraint is not considered in the model and so, the extension of the robot in simulation is constrained only by the stiffness of the structure and load mass. Thus, this modeling uncertainty can be compensated by adding a closed-loop controller, which is introduced in the following chapter.

Chapter 4

Closed loop control of soft, continuum manipulators

Introduction

In the previous chapter, the kinematics of continuum manipulators based on the inverse FEM simulation were presented and the limitations of the model in terms of model uncertainties were discussed. The aim now, is to compensate for the limitations of the kinematic modeling given in the previous chapter in terms of accuracy of the position tracking and rapidity by using a kinematic controller.

Controllers for soft manipulators have been investigated in the past with the intention of rejecting non-linear behaviours and model uncertainties that result from the complex dynamics of the manipulators (see 2.3.3). Control based on energy formulations [START_REF] Ivanescu | [END_REF]], modelless approaches [START_REF] Yip | [END_REF]] and feedback controllers [Penning 2011[START_REF] Penning | [END_REF]] have been proposed before with the intention of achieving accurate positioning of the manipulators in presence of non-modeled dynamics. In this chapter we propose a feedback control scheme for soft, continuum manipulators. The control is based on the Jacobian estimation of the robot, obtained from the FEM model simulation and 68Chapter 4. Closed loop control of soft, continuum manipulators the implementation of a state estimator. We have tested the closed-loop controller on the CBHA, but the same kind of control can be used with any soft, continuum manipulator as long as it can be modelled using the FEM-based method explained in the previous chapter.

In the next section, the feed-forward control based on the inverse kinematic model of the CBHA is presented in order to highlight the need for closed-loop controllers.

Feed-forward control of continuum manipulators

Based on the inverse kinematic model of the CBHA presented in the previous chapter, it is possible to use a real-time simulation of the robot to directly control it. A feed-forward controller which, in the case of the CBHA, computes the pressures required to reach a desired position for the end-effector is implemented using the relationship given in Eq. 3.16 and 3.17 We call this type of relationship a pressure-to-position model, since it relates these two spaces. In order to describe accurately the implemen-4.2. Feed-forward control of continuum manipulators 69 tation of this type of controller, let us consider the discrete form of Eq. 3.16. Assuming the external forces to be constant, we have

δ e,k+1 = δ e,k + W ea (x k )Δλ a,k+1 (4.1)
δ a,k+1 = δ a,k + W aa (x k )Δλ a,k+1 (4.2)
where δ e,k+1 is the current end-effector position and Δλ a,k+1 is the increment of the actuators force contribution at time k + 1, x k is the position of the nodes at time k, also referred as the state of the robot and W ea is the coupling between end-effector and actuators. In order to describe the system in a language that is more familiar to robotic control theory, in the following, the matrix W ea will be recalled as J ea , which is the Jacobian matrix between the inputs of the system (the pressure in the cavities) and the output error. In this way, Eq. 4.1 can be re-written as

δ e,k+1 = δ e,k + J ea (x k )Δλ a,k+1 (4.3)
The error is defined as

e k+1 = δ e,k+1 -δ d e,k+1 (4.4)
where δ d e,k+1 is the desired position. The control signal sent to the robot (pressure commands) computed at each simulation time is

λ a,k+1 = λ a,k + Δλ a,k+1 (4.5)
The configuration of the robot at time k + 1 is updated by

x k+1 = x k + K -1 k H T a Δλ a,k+1 (4.6)
where the matrix H a results from the imposition of the model constraints as Lagrange multipliers, as explained in 3.4. In this controller implementation, the estimation of the Jacobian matrix for the real robot performed by the simulation is inaccurate. This inaccuracy is due to the fact that a model will never match perfectly the real system because of unmodelled phenomenons, parameter uncertainties and disturbances. This results in a significant disparity between the Jacobian of the robot and the Jacobian estimated by the simulation. The results show a maximum error of 23mm, which represents the 11% of the total length of the manipulator. The pressure-to-position model used in this controller is highly sensible to model disparities, since the pressure constraints are no longer modelled as Dirichlet boundary conditions [START_REF] Reddy | [END_REF]] in the FEM simulation. Moreover, in the case of the CBHA, because of the lack of sensors, we cannot be certain that the pressures computed by the simulation are the ones applied to the cavities.

One way to correct this type of model is the use of a closed-loop control strategy. To this end, some questions arise in the development of the controller: First, what type of controller could be used? and more importantly, Is the closed-loop system robustly stable?

Closed-loop control of continuum manipulators

In the development of the closed-loop controller for the CBHA, a fundamental change in the nature of the model used in simulation is performed. As explained previously, the pressure-to-position model used in the feed-forward control required an accurate control of the pressures applied to the cavities. Instead, the pressure-to-length model given in 3.7.2 is used. In this way, the inputs of our model become the lengths of the sensors that correspond to a desired position. Eq. 4.1 is replaced by

δ s,k+1 = δ s,k + J sa (x k )Δλ a,k+1 (4.7)
where J sa is the Jacobian matrix between sensors and actuators.

Instead of imposing a desired end-effector position δ d e , a reference computation is performed in order to obtain the sensor lengths δ d s required to achieve said position (Fig. 4.3). The reference computation provides the desired motion in terms of desired sensor lengths. The desired trajectory can also be set by the user by imposing the motion to the trunk and recording the sensor lengths. When these desired sensor lengths are provided, we can propose a closedloop approach, as shown in Fig. 4.4.

reference computation

In Fig. 4.4 the blue blocks represent the computations performed by simulation. The output of the reference computation (Fig. 4.3) is the input of the closed-loop. Two simulations executing simultaneously are implemented in the closed-loop system: One main simulation that computes the Inverse kinematic model and a second simulation that acts as a state estimator for the system. The state estimator is the Forward kinematic model simulation of the robot that computes an estimated In this way, we make sure that the configurations of both simulation model and the manipulator are the same before the estimation of the Jacobian is computed.

Closed-loop control law design

The tracking error e k in the closed-loop system is computed as

e k = δ s,k -δ d s,k (4.8)
with δ d s,k represents the desired lengths of the sensors and δ s,k represents the current lengths in the robot. We define the control vector v k as

v k = Ĵsa (x k )r k (4.9)
where Ĵsa (x k ) is the estimated Jacobian matrix between the sensors and actuators and r k = Δλ a,k+1 . Using Eq. 4.9, the kinematic model can be rewritten as

δ s,k+1 = δ s,k + v k (4.10)
The control law is based on proportional integrative strategy, therefore, the control vector v k is designed in the sensor space as

v k = -k p e k -k i h k (4.11)
with k p and k i being the proportional and integrative gains of the controller, respectively. The integrative term h at time k + 1 is computed as

h k+1 = h k + e k (4.12)
Then, the control allocation based on a Quadratic Programming (QP) formulation [START_REF] Johansen | [END_REF]] is employed to find a unique solution to

r k = Ĵ+ sa (x k )v k (4.13)
where Ĵ+ sa is the pseudo-inverse of the estimated Jacobian. In practice, as Ĵsa (x k ) may not be fully invertible, we introduce a variable O defined as

O = Ĵsa (x k )r k -v k (4.14)
Using O, the QP problem formulation (3.5.3) becomes

min v k (O T O) (4.15)
the resulting r k will be the best possible inversion of Eq. 4.9 in the least square sense. In addition, the QP formulation allows to define constraints like actuator saturation or positive direction of actuation. Using Eq. 4.11 in Eq. 4.13, r k is rewritten as

r k = -Ĵ+ sa (x k )(k p e k + k i h k ) (4.16)
Using Eq. 4.16 in Eq. 4.7, the closed-loop system is defined as

e k+1 = e k + J sa (x k )r k (4.17)
which in the ideal case in which Ĵsa is invertible, can be written as

e k+1 = e k + v k (4.18)
The system of Eq. 4.18 is a first order discrete model that can be controlled with any standard controller. We choose the control strategy to be based on proportional-integrative control law because we want to 74Chapter 4. Closed loop control of soft, continuum manipulators improve the convergence rate and remove any steady state error (in the sensor space at least).

The gains of the controller are tuned by extensive testing using a simulated version of the closed-loop; we use a direct FEM simulation of the manipulator that takes any external force as input and computes the displacement field based on these forces. No computation of the inverse of the stiffness matrix K -1 is required. This direct simulation is put into the place of the real robot in Fig. 4.4. The control input computed by the simulation is applied to the direct simulation and the performance of the robot is observed.

In order to communicate the different simulations, a couple of clientserver scripts are written in Python based on the sockets provided by ZeroMQ communication protocol [Hintjens 2013], that allows for a fast communication between the simulations with an average latency of 75µs.

After testing, the selected gain values are k p = 0.14 and k i = 0.0003 as a compromise between the rise time of the signal and its overshot. show that both, the simulation of the robot and the robot itself have the same settling time t s ≈ 400 simulation steps. We can also see that the curve that represents the measured value of the lengths in the robot jumps between two values. This behavior is a consequence of the poor resolution of the string potentiometers. While the precision of the sensors cannot be improved, the signal of the sensors can be filtered by using a numerical implementation of an Infinite Impulse Response filter [START_REF] Rabiner | [END_REF]]. However, after implementing the filter, a considerable delay in the measuring signal was noticed. This delay hampers considerably the performance of the closed-loop system and, therefore, we opted to not implement it in the final version of the controller. 5 also shows a different response in the transitory stage of the curve of the measured lengths. This response can be attributed to different factors; first, there is the time required to compute the configuration of the manipulator from the measured sensor lengths; second, there is a time delay for the desired pressure to be applied to the robot, and finally, the plasticity of the material from which the manipulator is built (polyamide-nylon), that introduces a dampening effect at large deformations, which is not represented in our FEM model. There is hysteresis in the behavior of the pneumatic cavities. On the other hand, the pneumatic valves that control the pressure inside the actuators have a small dead zone, so, when the manipulator starts its motion from a zeropressure condition, very small increments in the pressure do not produce any motion until this dead zone is surpassed, which is not considered in the FEM model.

Robustness analysis

Because of modeling uncertainties, the estimated Jacobian matrix Ĵsa (x k ) is, in general, different from the Jacobian of the robot J sa (x k ). We introduce the vector ω k that represents the disparities between the real Jacobian and the estimated Jacobian. We call this vector the inversion error and is defined as

ω k = [I -J sa (x k ) Ĵ+ sa (x k )]v k (4.19)
Then, the closed-loop system is re-written as

e k+1 = e k + v k + ω k = e k -k p e k -k i h k + ω k (4.20)
The disturbed closed-loop system is

� e k+1 h k+1 � = � I -k p I -k i I I I � � e k h k � + � I 0 � ω k (4.21)
It can be disturbing that we end up with such a simple linear system. We emphasize to the reader that the non-linearities are taken into account by the two simulation blocks (FKM and IKM in Fig. 4.4) in the closedloop control. In Eq. 4.21, we are writing the system in terms of e k and h k and if the model was perfect, the system would be trivial. However, we can have modeling errors, that is why, in the following, we will prove that the control is robust to these modeling uncertainties ω k .

To simplify the notation of the problem, we define the following vectors Using this notation, matrix ω k is written as

X k+1 = � e k+1 h k+1 � , X k = � e k h k � , D = � I 0 � and F = � k p k i � (4.22) Also � I -k p I -k i I I I � = A -BF (
ω k = [I -J sa (x k ) Ĵ+ sa (x k )]FX k (4.25)
78Chapter 4. Closed loop control of soft, continuum manipulators

If the Jacobian is bounded and invertible, then there always exist a parameter γ such that

ω T k ω k = X T k F T [I-J sa (x k ) Ĵ+ sa (x k )] T [I-J sa (x k ) Ĵ+ sa (x k )]FX k ≤ γ 2 X T k F T FX k (4.26) with [I -J sa (x k ) Ĵ+ sa (x k )] T [I -J sa (x k ) Ĵ+ sa (x k )] ≤ γ 2 I (4.27)
For the proof of stability, we use Lyapunov's second method of stability [START_REF] Lyapunov | The general problem of the stability of motion[END_REF]]. We define the Lyapunov candidate function as

V = X T k PX k (4.28)
where P is the Lyapunov matrix with the properties

P T = P > 0 (4.29)
From Eq. 4.28 and the notation given in Eq. 4.22, the variation of the Lyapunov function is defined as

ΔV = X T k+1 PX k+1 -X T k PX k (4.30)
Using Eq. 4.24, Eq. 4.30 is re-defined as

ΔV = ((A -BF)X k + Dω k ) T P((A -BF)X k + Dω k ) -X T k PX k (4.31)
By making

A -BF = C (4.32)
Eq. 4.31 is written as

ΔV = (CX k + Dω k ) T P(CX k + Dω k ) -X T k PX k = X T k C T PCX k + X T k C T PDω k + ω T k D T PCX k + ω T k D T PDω k -X T k PX k (4.33)
Reverting the notation in Eq. 4.24, Eq. 4.33 can be written in matrix form as

ΔV = � X k ω k � T � (A -BF) T P(A -BF) -P (A -BF) T PD D T P(A -BF) D T PD � � X k ω k � (4.34)
For the proof, we introduce an accessory parameter α ≥ 0 in Eq. 4.26, such that

Υ = αω T ω -αγ 2 X T k F T FX k < 0 (4.35)
From Eq. 4.35, the left hand side of the inequality is written in matrix form as

Υ = � X k ω k � T � -αγ 2 F T F 0 0 αI � � X k ω k � < 0 (4.36)
Adding and subtracting this term to Eq. 4.34 allows us to find a bounding for ΔV as

ΔV -Υ + Υ = � X k ω k � T Q � X k ω k � + Υ (4.37) with Q = � (A -BF) T P(A -BF) -P + αγ 2 F T F (A -BF) T PD D T P(A -BF) DPD T -αI � (4.38)
We know from Eq. 4.35 that Υ < 0. Therefore, if Q is definite negative, then ΔV < 0. To prove the closed-loop system to be stable, the values for matrix P > 0 and α ≥ 0 need to be found such as matrix Q is definite negative, given the predefined values of the boundary parameter γ and the tuned controller parameter k p and k i . To this end, a Linear Matrix Inequality [START_REF] Boyd | [END_REF]] Solver called SeDuMi [START_REF] Sturm | [END_REF]] is used in the software Matlab. In order to describe the LMI given by Eq. 4.32, Yalmip [Lofberg 2004], a toolbox for optimization that is compatible with Matlab is employed. Given a value of γ = 0.98 and the gain values k p = 0.14 and k i = 0.0003, the LMI was solved successfully. The resulting matrix P and the parameter α that make matrix Q negative definite are Using the LMI solver, we can also compute the maximum value of γ, which provides an insight on the robustness of the closed-loop system.

After some iterations we have max γ = 0.98685 < 1 (4.40)

One can see that if γ > 1 the sign in the control input could potentially change which would lead to instability, so γ = 0.98685 is close to the limit case. The proposed closed-loop system is robustly stable and can handle high Jacobian inversion errors in the change of control variables.

Conclusions of the chapter

In this chapter, feed-forward and closed-loop controllers based on the FEM model for the CBHA were studied. Given the fact that the estimation of the Jacobian matrix depends on the configuration of the robot, the feed-forward control was ineffective in terms of precise positioning.

Based on this observation, a closed-loop control strategy that relies on the forward kinematic model as a state estimator was proposed. The state estimator is able to apply the configuration of the robot to the simulation. The closed-loop was tuned using a direct FEM simulation of the CBHA and then implemented in the real robot. The performance of the controller was showed and the proof of stability was given.

As explained in previous chapters, the model employed in the controller works under quasi-static conditions, which means that it is assumed that once the robot reaches a desired position, a period of stabilization is required before moving to the next position. The natural next step in the control of continuum manipulators would be to extend the controller to the dynamic case in which the inertial effects of the manipulator are considered. This would allow for the control of the acceleration of the robot, required for more complex applications. In general however, in order to capture high frequency vibrations of the structure, a higher sampling rate is required.

Conclusions of the chapter 81

In 4.3, we mentioned that a reference computation was performed in order to change the spaces related by our models. In this case, the robot is actuated directly to achieve a certain cartesian position, then, a sensor reading is performed to obtain the lengths that correspond to said position. However, thanks to the implementation of the state estimator, an end-effector trajectory can be set manually by manipulating the trunk directly while recording the lengths with the estimator. This manner of pre-computing the reference is more intuitive and can be used when the control commands for a desired position are not known.

In this chapter, the closed-loop control strategy for continuum manipulators was presented and the proof of stability was given. However, in the robustness analysis of the controller, the only assumption imposed on the inversion uncertainty ω k was its norm boundedness. No structure was assumed on ω k . This simplification of the uncertainties may lead to extreme conservative results. For example, consider a system composed by two actuators with a proportional gain k p described by The LMI proposed before is unfeasible with these values but if we study this model it can be proved to be stable. The system described in Eq. 4.41 and 4.42 is not realistic in the context of soft continuum manipulators, since more often than not the actuators are coupled. Nevertheless, the example helps to highlight the need of structured system uncertainties [Doyle 1982]. For a structured system, the inversion uncertainty takes the form

ω k = EΔFX (4.44)
where E and F are called uncertainty weighting filters. The boundedness is now applied to the uncertainty structure Δ ≤ γI. Special procedures have been developed to deal with structured bounded uncertainties [START_REF] Morari | [END_REF]].

Chapter 5 The main objective of this work was to provide a methodology to study and model soft, continuum manipulators based on computational mechanics. The methodology was used in the kinematic modelization of a continuum manipulator, the CBHA. We provided new control methods based on such models. We have validated the methodology with experiments that involve the real-time simulation of the FEM model. In the following, the conclusions for each chapter are given.

Conclusion and Perspectives

Summary of conclusions

In chapter 2, a concise and complete view of the field as it is was given.

The main problem on the modelization of continuum manipulators was introduced and a brief study on the most popular approaches towards the modeling in the literature was given and a qualitative evaluation of these approaches was presented. In order to make room for our contribution in the field, the role that external forces play in the motion of continuum manipulators was also presented. Also, a taxonomy of the different types of designs was presented in order to highlight the main actuation differences. In our opinion, the most important part of the state of the art was the introduction of the main issue that most of the modeling approaches fail to capture: the deformation mechanics! Continuum manipulators produce their motion by deformation, and yet, very few approaches make emphasis on this feature.

Chapter 5. Conclusion and Perspectives

Chapter 3 started as an introduction on continuum mechanics in an attempt to show that the study of deformations seems like a more natural and compatible tool to model continuum manipulators. The Finite Element Method was presented as a method to solve differential equations with boundary conditions; a problem that is encountered when studying continuum media. Special emphasis was put in the concept of discretization, that allow us to reduce the infinite number of degrees of freedom of continuum manipulators. The methodology we use towards the modeling and control of soft, continuum manipulators was introduced and the model reduction based on the projection of the FEM model was explained. The model of sensors, actuators and end-effector as constraints in the FEM was explained. The mechanical coupling between the different constraints, provided by the compliance operators, was shown to be a necessity in the computation of the Jacobian matrix for continuum manipulators. The model of actuators and sensors based on Lagrange multipliers was used in the implementation of the kinematic models for the CBHA. The CBHA manipulator was introduced as a case of study to exemplify the implementation of the methodology. While our methodology does not produce analytic solutions for the kinematic model, the numerical approximation of the kinematic solution was shown superior to those obtained by 2 other geometric methods in terms of tracking position.

As a conclusion for chapter 3, some of the uncertainties in the model of the CBHA were discussed. This gave room for the implementation of the closed-loop control presented in chapter 4. The simulation of the FEM model was used to control the manipulator in real time. The discussion about the feed-forward control of the CBHA provided a good example of the need for closed-loop controllers. While there are a lot of non-linear behaviors that are currently not included in the model, the closed-loop control is still able to provide good performance in terms of the positioning of the manipulator. The proposed controller was proved to be robustly stable.

As a conclusion to this work, we have shown the feasibility of using FEM simulation as a generic methodology for modeling and control of soft, continuum manipulators, based on the premise that the study of deformations is fundamental if one seeks to characterize this type of robots.

One interesting venue for future research of the methodology presented in this manuscript is its application to the design of soft, continuum manipulators. Since the simulation of the backbone of continuum robots can be done in a fast manner. The simulation can be used to reduce the number of design iterations performed to achieve a certain behaviour. The fabrication of soft actuators, in particular fluidic actuators, is a complex task. Being able to observe the behavior that a certain design has, prior to its fabrication, is of great help when undertaking the design of a new prototype. This approach towards the design and integration of soft, continuum manipulators has been explored in the past, with the design and construction of the FeTCh manipulator.

The FeTCh manipulator

The prototype of the FeTCh manipulator was designed with the intention of using the principle of antagonistic actuation to achieve the rigidification of its structure [START_REF] Shiva | [END_REF]]. The backbone of the robot is composed by actuators in an intrinsic fashion, while including extrinsic tendon actuators. Each section of the robot is formed by three pneumatic actuators arranged in parallel configuration between two rigid platforms that guide the deformation of the actuators. This configuration gives the section 2 axis bending, and also extension. The arrangement of the actuators also helps to counter the lack of shearing stiffness that is characteristic to pneumatic actuators. The described section of the manipulator is shown in Fig. 5.1.

The extrinsic actuation of the manipulator is composed by tendons that are pulled by a servomotor-pulley system. These tendons are routed through the rigid platforms at 120 • around the rigid platform at the end of each section and apply a force that is regulated. The inclusion of the tendons provides more stability to the manipulator when it moves, but also provide a degree of selective rigidification of the structure.

The idea behind the actuator design was to create an extensor that can elongate when pressurized, as opposed to the more traditional contractor artificial muscles. The extension of the actuator is caused by the unfolding of its walls, rather than pure deformation caused by the pressure. The simulation of the actuator was used in order to observe The design allows for the positioning of the final rigid platform, but also for its orientation. Further testing on the rigidification effect based on the FEM model is still required but the preliminary results are promising.

The simulation of soft structures prior its fabrication allows also for the testing of the design robustness. Since the integration of soft robots in general is not a perfected process, defects in the soft parts are often encountered. For example, the most common flaw in soft pneumatic actuators are artificial heterogenous regions formed by air bubbles. In the simulation, the inclusion of these regions is easily done. This allows the validation of designs when fabrication flaws are expected. In addition, the range of actuation can also be determined, this allows for the selection of key components like the servos for the tendons and the type of pneumatic manifolds to be used in the back-end of the manipulators.

Once the final assembly for the manipulator has been established, the identification of the task space of the robot can also be performed.

Perspectives

The possible improvements in the current work regarding the modeling of continuum manipulators are numerous and it is difficult to mention them all here. However, based on the experience we have with the CBHA, some immediate modifications can be performed to improve the model without making it too complex for fast computation. For example, it is observed that at different temperatures, the manipulator behaves slightly differently; at higher temperature, the plastic material is more ductile. The inclusion of the model of thermal expansion in the element potential energy functional Π (Eq. A.52), from the FEM, can be done without changing its general formulation. By adding the term

h = � V B T Eε T dV (5.1)
where B is the relation between strains and nodal displacements and E is the elasticity matrix, we can represent the fictitious forces that correspond to the thermal expansion of the material. Of course, in small scale manipulators, this term can be completely neglected, but for applications that require bigger manipulators, the effect of temperature in their behavior can be considerable. It is also observed that material fatigue affects the motion of the manipulator. After prolonged use of the robot, the pneumatic actuators expand more easily. Creeping and relaxation stages of plastic materials should also be considered in the model.

In terms of design, there is a great need for novel actuation technologies to be developed with the purpose of providing continuum robots with bigger force generation capabilities. At the present time, it is clear that soft manipulators have a niche of applications in the medical field.

However, this narrow niche can be expanded provided that this type of robots can compete with rigid ones in terms of force generation and payload handling.

As a final perspective, much work needs to be done in terms of the formalization of the theory behind continuum manipulators. In order to expand and teach soft robotics as a field, a unified and comprehensive terminology and methodology of modeling and control needs to be firmly established. It is our hope that the work presented in this brief manuscript contributes towards that important objective. A material particle of B, with a position vector P 0 in the initial configuration κ 0 will have a new position vector P i in κ i after deformation given by

P i = M (P 0 , t) (A.1)
The function M is called a deformation mapping and describes a 1to-1 relation between particle positions in the initial configuration and the new position in the deformed configuration. The displacement of a particle in a continuum can be defined as

u(P i ) = P i -P 0 = M (P 0 , t) -P 0 (A.2)
This description of motion of continuum is referred as the Lagrangian description (Fig. A.1).

Now consider an infinitely small line segment in the body B given by dP 0 and dP i before and after deformation, respectively, as shown in Fig.

A.2. The relationship between these two lines is given by the gradient of deformation F, defined as

F = � ∂M ∂P i � � = ∇u + I (A.3)
where ∇ is the gradient operator with respect to P i and I is the identity matrix. The determinant of F, J = detF is referred to as the Jacobian of motion. The displacement of particles in a continuum media can be considered, in the general case, as a result of a rigid body motion (a translation and a rotation) and a deformation. A measure of strain is the computation of the deformation part of the displacement, independent of the rigid body motion.

In order to arrive at the definition of the strain tensor, we define the right Cauchy-Green tensor as

C = F • F � (A.4)
which gives us the square of local displacement due to the deformation. Similarly, the left Cauchy-Green tensor is defined as

B = F � • F (A.5)
Using Eq. A.4, we define the Green -Lagrange strain tensor in vector form as

E = 1 2 (F • F � -I) = 1 2 (C -I) (A.6)
with E = 0 when the displacement is a result of a rigid body motion. Eq. A.6 is expressed in Cartesian coordinates (X 1 , X 2 , X 3 ), as

E i,j = 1 2 � ∂u i ∂X j + ∂u j ∂X i + ∂u k ∂X i ∂u k ∂X j � , ∀ 1 ≤ (i, j) ≤ 3 (A.7)         σ 1 σ 2 σ 3 σ 4 σ 5 σ 6         =         C 11 C 12 C 13 C 14 C 15 C 16 C 21 C 22 C 23 C 24 C 25 C 26 C 31 C 32 C 33 C 34 C 35 C 36 C 41 C 42 C 43 C 44 C 45 C 46 C 51 C 52 C 53 C 54 C 55 C 56 C 61 C 62 C 63 C 64 C 65 C 66                 ε 1 ε 2 ε 3 2ε 4 2ε 5 2ε 6         (A.17) with C ij = C ji .
The components of strain are related to the components of stress by

ε i = W ij σ j (A.18)
where W ij are the coefficients of the compliance tensor W = C -1 . Eq. A.18 is written in matrix notation as

        ε 1 ε 2 ε 3 2ε 4 2ε 5 2ε 6         =         W 11 W 12 W 13 W 14 W 15 W 16 W 21 W 22 W 23 W 24 W 25 W 26 W 31 W 32 W 33 W 34 W 35 W 36 W 41 W 42 W 43 W 44 W 45 W 46 W 51 W 52 W 53 W 54 W 55 W 56 W 61 W 62 W 63 W 64 W 65 W 66                 σ 1 σ 2 σ 3 σ 4 σ 5 σ 6         (A.19)
Orthotropic materials Orthotropic materials have three mutually orthogonal planes of material symmetry, meaning that their material properties differ when measured from different directions. In this case, the elastic coefficients are reduced to 9. Most of the time, experimental characterization of material is performed using known loads; it is then convenient to express the relation of stress and strain in terms of the compliance. Eq. A.19 for orthotropic materials has the form

        ε 1 ε 2 ε 3 2ε 4 2ε 5 2ε 6         =         W 11 W 12 W 13 0 0 0 W 21 W 22 W 23 0 0 0 W 31 W 32 W 33 0 0 0 0 0 0 W 44 0 0 0 0 0 0 W 55 0 0 0 0 0 0 W 66                 σ 1 σ 2 σ 3 σ 4 σ 5 σ 6         (A.20)
In the experimental characterization of linear elastic materials, the coefficients of the compliance tensor are computed using the Young' s modulus, the shear modulus and the Poisson ratio [Wortman 1965]. After introducing these engineering constants, Eq. A.20 takes the form

        ε 1 ε 2 ε 3 2ε 4 2ε 5 2ε 6         =          1 E 1 -v 21 E 2 -v 31 E 3 0 0 0 -v 12 E 1 1 E 2 -v 32 E 3 0 0 0 -v 13 E 1 -v 23 E 2 1 E 3 0 0 0 0 0 0 1 G 23 0 0 0 0 0 0 1 G 13 0 0 0 0 0 0 1 G 12                  σ 1 σ 2 σ 3 σ 4 σ 5 σ 6         (A.21)
where E 1 , E 2 and E 3 are the Young's moduli in the three main material directions, v ij is Poisson's ratio, defined as the ratio of traverse strain in the jth direction to the axial strain in the ith direction when stressed in the i-direction, and G 23 , G 13 and G 12 are the shear moduli in the 2-3, 1-3 and 1-2 planes respectively.

When computed experimentally, the Young's modulus, which is a measure of linear stiffness, is represented by the slope of the linear part in the stress-strain relationship, while the Poisson's ratio is a measure of the change in volume of a material.

Isotropic materials If the material properties are independent from the direction, meaning that

E 1 = E 2 = E 3 = E , G 23 = G 13 = G 12 = G , v 12 = v 23 = v 13 = v (A.22)
the material is said to be isotropic. For isotropic materials, the stressstrain relationship is reduced to

σ ij = E 1 + v ε ij + vE (1 + v)(1 -2v) ε kk δ ij (A.23)
where δ ij is called the Kronecker Delta and it is used as a simple representation of the scalar product

δ ij = êi • êj (A.24)
that takes the values δ ij = 1 if i = j and δ ij = 0 otherwise. The vectors êi and êj are the orthonormal unit vectors in a right-handed basis system. Introducing the Lamé coefficients [Akamatsu 1991] .25) the constitutive equation for isotropic linearly elastic material takes the form

µ = E 2(1 + v) λ = Ev (1 + v)(1 -2v) (A
σ ij = 2µε ij + 3 � k=1 ε kk δ ij (A.26) A.1.1.

Equations of motion

The evolution of the motion in continuum bodies is derived from the Principle of Conservation of Linear Momenta commonly known as Newton's second law of motion that can be stated as the time rate of linear momentum of a collection of particles equals the net force exerted on the collection. This principle allow us to express the relationships between the primary unknown displacement field of the particles of a continuum body in terms of stresses and strains. The solution of these equations represent the response of the system to a set of input (known forces) data. In vector form, the principle is represented by

m dv dt = F (A.27)
where m is the constant total mass, v is the velocity, and F is the net force applied on the particles. For a body over the control volume Ω, Eq. A.27 takes the form In Cartesian rectangular coordinates, we have ∂σ ji ∂x j + ρf i = 0, ∀ 1 ≤ (i, j) ≤ 3 (A.36)

In the general case, Eq. A.34 is the partial differential equation that describes the motion of a deformable body. However, an analytical solution to Eq. A.34 is intractable for most cases. In turn, an approximation of the solution is computed using numerical methods.

A.2 Finite Element Method for linear elastic bodies

In order to reconstruct the sequence of configurations that an elastic body will take, given a set of forces acting on it, Eq. A.36 remains to be solved. This equation, however, cannot be solved analytically, since it involves an infinite number of degrees of freedom. In general, finding the solution for continuum elastic solids is done by making use of numerical methods provided by the field of Computational Mechanics, since analytical solutions exist only for trivial problems. One of these numerical methods is the Finite Element Method (FEM).

In the finite element method, a geometrically complex domain Ω is viewed as being comprised by a finite set of subdomains (elements) Ω e over which the solution of the governing equation is approximated using variational methods1 . The division of the domain in subdomains allows for an accurate representation of complex geometry.

We can summarize the Finite Element Method in three main steps:

1. The division of the domain of interest in smaller subdomains. We call this step the discretization step 2. The approximation of the solution at the element level as a linear combination of nodal values and shape functions 3. The assembly of the elements that allows to approximate the global solution over the whole original domain

In the following, these three steps are briefly explained for the formulation of the displacement-based FEM [Bathe 2006], which is used in the manuscript to derive the kinematics of soft manipulators. If more detailed explanation is required regarding the general formulation of the problem using variational methods, the reader is advised to refer to [START_REF] Reddy | [END_REF]].

A.2.1 Discretization of the domain The first step in FEM is the partition of the domain of interest Ω into smaller, geometrically independent domains called finite elements Ω e . The collection of finite elements that form the discrete domain is called mesh (Fig. A.3). In this step, the first source of error in a finite element analysis is introduced: In most cases, the discrete representation will differ geometrically from the original domain. It is evident that this error can be diminished by reducing the dimensions of the elements to better conform to Ω, but this results in an increased number of elements involved in the computation which, in turn, require more computational resources. A compromise between an accurate representation of the complex geometry and the number of elements is needed. In fact, mesh refinement and convergence is still an active research topic.

While there are some commercially available applications for automatic mesh generation, it is important to mention the general set of rules that a finite element mesh must follow. This set of rules are:

Finite elements in a mesh do not overlap in space. This property goes by the name of local support

Elements with high aspect ratios should be avoided. The aspect ratio of a three-dimensional element is the ratio between its largest and smallest dimensions A physical interface ( In the following, the treatment given to an element is explained.

A.2.2 Element solution

At the element level, the physical quantities of interest, which in our case is the displacement field u at any point inside the element, is approximated by a linear combination of u at the nodes of the element. This approximation is expressed as

u e ≈ n � i=1 u e i φ e i (x) (A.37)
where u e is the displacement inside the element, u e i is the displacement of node i and φ e i are the interpolation functions that depend on the element geometry. In the following, we introduce the linear tetrahedron element and derive its interpolation models. In this case φ e i are taken to be linear polynomials. for isotropic materials the relation in expanded form simplifies to

        σ xx σ yy σ zz σ xy σ yz σ zx         = E (1 + v)(1 -2v)         1 -v 0 0 0 0 0 0 1 -v 0 0 0 0 0 0 1 -v 0 0 0 0 0 0 1 2 -v 0 0 0 0 0 0 1 2 -v 0 0 0 0 0 0 1 2 -v         (A.47)
where E is the Youngs modulus and v is the Poisson ratio. Introducing Eq. A.43 and A.46 into the Total Potential Energy functional restricted to the element volume and assuming that the elasticity modulus do not vary over the element, we can compute the element stiffness matrix as K e = V B T CB (A.48)

A.2.3 Assembly

The purpose of the finite element solution of the elastic problem is to find the displacement field u that provides minimum to the total potential energy functional .49) with the second term corresponding to the volume forces, where f is the body force vector and the third term corresponding to the surface forces where m is the surface force vector. The displacement boundary conditions do not appear in the functional Π; boundary conditions are implemented after the assembly of the finite element equations.

Π Π = � V 1 2 ε T σdv - � V u T f dV - � S u T mdS (A
Then, we fill the mth row and mth column of matrix K with zeros, except for the main diagonal term, which is replaced by 1 K mj = 0, j = 1...n, K im = 0, i = 1...n,

K mm = 1 (A.62)
After boundary conditions are applied, matrix K has full rank.

A.2.5 Solution method

For practical applications, the system Eq. A.55 results in a large number of simultaneous linear algebraic equations. Fortunately, matrix K possesses some properties that enables an easy computation of the system of equations. For instance, K is symmetric, which means that it is enough to store only half of the matrix entries including those in the diagonal. Also, K is definite positive after setting the boundary conditions (Eq. A.62) and highly sparse, which makes storage and computations economic.

Solution methods for the simultaneous linear system of equations come in two different types: direct methods, such as those based on the LDL T decomposition or LDU method with profile matrix [Bunch 1971], which are usually used for problems of moderated size, and iterative methods that require less computations for large scale problems. In practice, we use the iterative method of the preconditioned conjugate gradient method or PCG, since we often deal with problems involving thousands of degrees of freedom. Let us consider the example system given by Ax = b (A.63)

In the PCG, the convergence of the conjugate gradient is improved by the preconditioning of the equation system

M -1 Ax = M -1 b (A.64)
where M -1 is the preconditioning matrix which in some sense approximates A -1 . The simplest preconditioning is diagonal preconditioning, in which M contains only diagonal entries of the matrix A. There are many existing solvers and libraries developed to solve large sparse linear matrix systems.

computed matrix K, the compliance that corresponds to the rigid frames can also be computed as

C = JK -1 J T (B.1)
where

J = � 0 6×n I 6X6 � (B.2)
with n being the number of independent DoF. The identity matrix in J corresponds to the constraints of each DoF in the upper frame. Given two consecutive frames i and j, the equivalent stiffness matrix K eq ij can be expressed as

K eq ij = R i � H ij C -1 ij H T ij H ij C -1 ij C -1 ij H T ij C -1 ij � R T i (B.3)
where R i is the rotational matrix of frame i. Matrix H ij describes the relationship between linear displacement on frame j and angular displacement on frame i and has the form

H ij = � H 3×3 0 3×3 (x j -x i ) ∧ H 3×3 � (B.4)
with x i and x j are the position of frames i and j respectively and the symbol ∧ represents the skew symmetric matrix representation of a vector.

Using this approach, the FE model of a continuum manipulator can be replaced by a chain of frames linked by equivalent springs.

Figure 2 . 2 :

 22 Figure 2.2: Diagram of the hydrostatic muscle in the arm of Octopus from [Trivedi 2008b]. AN, axial nerve cord; AR, artery; CM, circumferential muscle layer; CT, connective tissue; DCT, dermal connective tissue; EP, epidermis; IN, intramuscular nerve; LM, longitudinal muscle fibres; OME, external oblique muscle layer; OMI, internal oblique muscle layer; OMM, median oblique muscle layer; SU, sucker; TM, traverse muscle fibres; TR, trabeculae; and V, vein.

  2.3 (left)) have a central elongated structure that supports the passage of actuation/transmission elements along the body of the manipulator [Burgner-Kahrs 2015]. Many single-backbone designs utilize tendons routed along the structure which are spaced by discs attached to the backbone as a way of transmission. The termination points of the tendons in this type of design dictate the lengths of a bending section. Multibackbone (Fig. 2.3 (right)) continuum manipulators are typically composed by a parallel arrangement of elastic elements which are constrained with respect to each other in some way.

Figure 2 . 3 :

 23 Figure 2.3: Single and Multibackbone continuum robots:(left) A single backbone tendon-driven steerable cardiac catheter from [Camarillo 2008].(Right) The DDU multibackbone manipulator designed by Simaan et al. [Simaan 2009].

Figure 2

 2 Figure2.4: The KSI hybrid actuated manipulator[Immega 1995]

Figure 2

 2 Figure 2.7: (Left) Extensor and contractor pneumatic muscles, (Right) An example of a locally actuated design, the Oct-Arm manipulator composed by pneumatic muscles, from [McMahan 2006]

Figure 2

 2 Figure 2.8: A concentric tube robot or Active Cannula composed by 3 pre-curved tubes, from [Wu 2017]

Figure 2 .

 2 Figure 2.12: The geometry of a spatial constant curvature curve segment

  Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . 32 3.2 Continuum mechanics framework . . . . . . . . . 33 3.2.1 Constitutive material law . . . . . . . . . . . . . 34 3.2.2 Forces in the continuum manipulator . . . . . . . 35 3.3 Finite Element Method . . . . . . . . . . . . . . . 36 3.4 FEM-based kinematics of soft manipulators . . 38 3.4.1 Constraint for the end-effector . . . . . . . . . . 39 3.4.2 Actuator constraint model . . . . . . . . . . . . . 41 3.4.3 Sensor constraint model . . . . . . . . . . . . . . 44 3.5 Reduced model in the constraint space . . . . . 45 3.5.1 Reduced compliance on the constraint space . . . 45 3.5.2 Coupled Kinematic Equations . . . . . . . . . . . 47 3.5.3 Inverse kinematic model solution by convex optimization . . . . . . . . . . . . . . . . . . . . . . . 48 3.6 Method implementation . . . . . . . . . . . . . . 49 3.6.1 Simulation framework . . . . . . . . . . . . . . . 50 3.6.2 Corotational FEM . . . . . . . . . . . . . . . . . 50 3.6.3 Mesh generation . . . . . . . . . . . . . . . . . . 51 3.6.4 Description of the Compact Bionic Handling Assistant . . . . . . . . . . . . . . . . . . . . . . . . 51 3.6.5 Simulation of the CBHA . . . . . . . . . . . . . . 53 3.7 Kinematic models . . . . . . . . . . . . . . . . . . 55 3.7.1 Forward kinematic models . . . . . . . . . . . . . 56 3.7.2 Inverse kinematic model . . . . . . . . . . . . . . 62 3.8 Conclusion of the chapter . . . . . . . . . . . . . 63

Figure 3

 3 Figure 3.1: In this example a tendon is pulled to create the motion of an elastic soft robot. Starting with the same geometry, the material stiffness has an influence on the kinematics (output vs input displacements).

  Figure 3.2: Illustration of the initial and deformed configurations of the continuum manipulator.

Figure 3

 3 Figure 3.3: Illustration of a mesh. The volumetric mesh is created from a surface model of the manipulator.

  Figure 3.4: Illustration of an element. One of the many (tetrahedra) elements composing the mesh.

Figure 3 . 5 :

 35 Figure 3.5: Illustration of the end-effector constraint. δ e represents the difference between current and desired position.

Figure 3 . 6 :

 36 Figure 3.6: Cable actuation. d b and d a on the figure, represent the direction of the tendon before and after the cable guide, respectively, which are used to compute the normal forces at the guides.

Figure 3

 3 Figure 3.7: A pneumatic actuator. The cavity at rest (left) is then pressurized with air (right). The pressure applies a normal force to the walls of the cavity.

Figure 3 . 8 :

 38 Figure 3.8: Example of an embedded sensor. The sensor string is routed though n frictionless guides.

  shift between the effector point(s) position and the desired position introduced in section 3.4.1. One can also evaluateδ free a = δ a (x free ) (3.15)which represents the position of the actuated points without actuation effort.

  ) as: δ e = W ea W -1 aa (δ aδ free a ) + δ free e (3.21) Chapter 3. FEM-based model of Continuum Manipulators Equation (3.21

Figure 3

 3 Figure 3.9: The CBHA by Festo Robotics. (a) The anatomy of the Compact Bionic Handling Assistant. (b) A section of the manipulator, composed by 3 pneumatic actuators and their correspondent length sensor.

Fig. 3

 3 Fig. 3.9 (b). By applying different pressures to the bellows, each section can bend independently (Fig.3.10).

Figure 3 .

 3 Figure 3.10: Examples of configurations of the CBHA based on the pressure inside the pneumatic actuators.

Figure 3 .

 3 Figure 3.11: Visual surface model of the CBHA (left) and the FEM mesh used in simulation (right)

Figure 3 .

 3 Figure 3.13: Illustration of the forward kinematic model (inputs and output)

Figure 3 .

 3 Figure 3.14: Parametrization of the Constant curvature model of the CBHA

Figure 3 .

 3 Figure 3.16: Estimated trajectory of the FEM-based approach compared to the measured trajectory

Figure 3 .

 3 Figure 3.17: Estimated trajectory of the Constant curvature approach compared to the measured trajectory

Figure 3 .

 3 Figure 3.19: Illustration of the inverse kinematic model (input and output)

Figure 3 .

 3 Figure 3.21: Collision of the outer walls of the cavities when bent at maximum angles

  . . . . . . . . . . . . . . . . . . . . . 67 4.2 Feed-forward control of continuum manipulators 68 4.3 Closed-loop control of continuum manipulators 71 4.3.1 Closed-loop control law design . . . . . . . . . . 72 4.3.2 Robustness analysis . . . . . . . . . . . . . . . . 76 4.4 Conclusions of the chapter . . . . . . . . . . . . . 80

  Figure 4.1: Illustration of the feedforward controller based on the pressure-to-position model

70Chapter 4 .

 4 Closed loop control of soft, continuum manipulatorsTo show the performance of this type of control, a semicircular trajectory formed by via points inside the task space of the robot is computed. Each point of the trajectory is introduced in the simulation sequentially as a desired position for the end-effector and the actuation computed by the simulation is then used as set point for the pressure regulator of the robot. The results from this test are shown in Fig.4.2.

Figure 4 . 2 :

 42 Figure 4.2: Comparison between desired and measured trajectories using the simulation of the CBHA as a feed-forward controller

Figure 4 . 3 :

 43 Figure 4.3: Reference computation that is performed before the feedback system.

Figure 4

 4 Figure 4.4: Closed-loop control of the CBHA based on IKM and FKM simultaneous simulations and the controller

  Fig 4.5 shows the lengths of the simulated robot and the real robot, both in closed-loop, given a pre-computed set points corresponding to an end-effector position inside the task space of the robot. The position is chosen so that the actuators are far from their saturation points. The model simulation and the real robot have different initial conditions. After 1000 simulation steps, the set points are changed in both the simulation and the real robot.Some noticeable features inFig 4.5 can be observed. First, the results

Figure 4 . 5 :

 45 Figure 4.5: Comparison of real and estimated actuator lengths of the CBHA. A second set point is applied to the system after 1000 simulation steps in order to observe the performance of the controller. The time step is 0.1s for the experiment.

Fig

  Fig 4.5 also shows a different response in the transitory stage of the curve of the measured lengths. This response can be attributed to different factors; first, there is the time required to compute the configuration of the manipulator from the measured sensor lengths; second, there is a time delay for the desired pressure to be applied to the robot, and finally, the plasticity of the material from which the manipulator is built (polyamide-nylon), that introduces a dampening effect at large deformations, which is not represented in our FEM model. There is hysteresis in the behavior of the pneumatic cavities. On the other hand, the pneumatic valves that control the pressure inside the actuators have a small dead zone, so, when the manipulator starts its motion from a zeropressure condition, very small increments in the pressure do not produce any motion until this dead zone is surpassed, which is not considered in the FEM model.

76Chapter 4 .

 4 Closed loop control of soft, continuum manipulatorsA second experiment was performed with the real robot in the loop. In this experiment an external unknown force was applied to the manipulator in order to see the uncertainty rejection capabilities of the controller. Fig.4.6 shows the results of this experiment.

Figure 4 . 6 :

 46 Figure 4.6: Measured lengths of the CBHA in closed-loop. An external force is applied to the manipulator after 1050 time steps. The time step is 0.1s for the experiment.

Fig. 4 .

 4 Fig.4.6 shows the controller disturbance rejection capabilities. After the external force is applied, the manipulator is able to return to the desired configuration.

e

  1k+1 = (1kp)e 1k + ω 1k (4.41) e 2k+1 = (1kp)e 2k + ω 2k (4.42) given the values of uncertainties ω 1k = 0, ω 2k = 10 5 e 1k (4.43)

Contents 5 . 1

 51 Summary of conclusions . . . . . . . . . . . . . . 83 5.2 The FeTCh manipulator . . . . . . . . . . . . . . 85 5.3 Perspectives . . . . . . . . . . . . . . . . . . . . . . 88

Figure 5

 5 Figure 5.1: Single section of the FeTCh manipulator showing different configurations depending on the pressure inside the actuators

Figure 5 . 2 :

 52 Figure 5.2: Simulation of the soft actuator.After the validation of the design by the simulation, the actuator was built and tested. The actuator has a maximum extension of 81.25% when pressurized with 0.41bars. The final design of the manipulator is achieved by serially connecting 3 sections like the one shown in Fig.5.1.

Figure 5 . 3 :

 53 Figure 5.3: The FeTch manipulator.

  Figure A.1: Illustration of the Lagrangian description of motion between initial and deformed configuration for an elastic body

  Figure A.2: Illustration of the gradient of deformation in the Lagrangian description of motion

  ds is the vectorial representation of an area element, ρ is the mass density and Γ is the boundary surface of the body. Let f be the external body force per unit of mass, which acts on the distribution of mass inside the body. Considering an elemental volume dΩ inside Ω, the body force of dΩ is equal to ρdΩf . The total body force in the control volume is� Ω ρf dΩ (A.29)Let m be the surface force per unit area. The surface force acting on an elemental surface ds of the volume is mds. The total surface force acting on the closed surface of Ω is� Γ mds (A.30)rewriting m in terms of Cauchy's stress tensor (Eq. A.11), and using the divergence theorem, we have transport theorem[Cengel 2006], the global form of the equation of motion of the body occupying the region Ω with bounding surface Γ and acted upon the surface force m and body force f denotes the material time derivative. The local form is∇ • σ σ σ + ρf = ρ

Figure A. 3 :

 3 Figure A.3: Discretization of the domain. A continuum domain (a) is approximated by triangle elements (b) which produces a discretization error. The mesh is refined (c) by increasing the number of elements. This in turn reduces the error.

  Figure A.4: Illustration of the mesh rules. Elements should not cross physical interfaces formed by a change of material (a). Elements with high aspect ratios (b) should be avoided.

  material, the interpolation functions are simply the tetrahedron coordinates ζ i . The displacement field over the tetrahedron u e is defined by the three components u x , u y and u z . These are linearly interpolated from their node valuesu x2 u x3 u x4 u y1 u y2 u y3 u y4 u z1 u z2 u z3 u z4 A.40 and A.41, we obtain the isoparametric definition of the tetrahedron as a displacement model 2 x 3 x 4 y 1 y 2 y 3 y 4 z 1 z 2 z 3 z 4 u x1 u x2 u x3 u x4 u y1 u y2 u y3 u y4 u z1 u z2 u z3 u z4 the strains and the nodal displacements is given byε = Bu e (A.43) where matrix B is constant over the element, and has the form1 0 0 c 2 0 0 c 3 0 0 c 4 b 1 a 1 0 b 2 a 2 0 b 3 a 3 0 b 4 a 4 0 0 c 1 b 1 0 c 2 b 2 0 c 3 b 3 0 c 4 b 4 c 1 0 a 1 c 2 0 a 2 c 3 0 a 3 c4 0 a 4 the volume of the element. Similarly to Eq. A.20, the stress field is related to the strain field by σ = Cε (A.46)

  Finally, I want to thank my cat Brisby for keeping me company during the nights while writing this manuscript. General Introduction . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Framework and context . . . . . . . . . . . . . . . . . . . . . . General Introduction . . . . . . . . . . . . . . . . . . . 1 1.2 Framework and context . . . . . . . . . . . . . . . . . .

	ix Pseudo-inverse of the estimated Jacobian matrix between sensors and Acronyms AM -Additive Manufacturing API -Application Programming Interface CAD -Computer Assisted Design CBHA -Compact Bionic Handling Assistant CC -Constant Curvature (model) CGAL -Computational Geometry Algorithms Library CPU -Central Processing Unit DoF -Degrees of Freedom FBG -Fiber Bragg Grating FEM FKM HA -Hybrid Approach (model) IKM -Inverse Kinematic Model KSI -Kinetic Sciences Inc. LMI -Linear Matrix Inequality ME -Maximum Error PDE -Partial differential Equation PI -Proportional Integrative (controller) QP -Quadratic Programming RMSD -Root-Mean-Square Deviation SOFA -Simulation Open Framework Architecture W-LAN -Wireless Local Area Network g -Gravity field G -h -H -H a -H e -H s -J ea -Jacobian matrix between end-effector and actuators J sa -Ĵsa -Ĵ+ sa -k -k i -k p -K -K e -m -M -M -n -p des -p ef f -P -Position vector Ω -Three-dimensional domain End-effector position in free configuration ω -Vector of inversion error Desired position σ -Cauchy stress tensor Vector normal to a surface Φ -Total potential energy functional Inertia matrix φ i -Shape (interpolation) function associated to node i Mapping function λ s -Vector of sensor constraint contributions Constant total mass λ e -Vector of end-effector constraint contributions Element stiffness matrix λ a -Vector of actuator constraint contributions Tangent stiffness matrix λ -Vector of constraints contributions Proportional gain ζ i -Tetrahedron natural coordinates Integrative gain ε -Infinitesimal strain tensor Spatial configuration Δλ a -Increment of actuators contribution actuators δ f ree ΔV -Variation of the Lyapunov function s -Sensor string length in free configuration Estimated Jacobian matrix between sensors and actuators δ f ree of end-effector in free configuration e -Vector representing the difference between current and desired position Jacobian matrix between sensors and actuators δ f ree in free configuration a -Cable length (cable actuator) or cavity volume (pneumatic actuator) Matrix of sensors constraint directions δ s -Sensor string length Matrix of end-effector constraint directions of end-effector Matrix of actuators constraint directions δ e -Vector representing the difference between current and desired position Matrix of constraints directions δ a -Cable length (cable actuator) or cavity volume (pneumatic actuator) Integrative term Γ -Boundary surface Shear modulus γ -Bounding parameter -Forward Kinematic Model F e -Vector of element loads -Finite Element Method xi Glossary a t -Area of triangle C -Stiffness tensor d a -Direction after cable guide d b -xii v -Poisson's ratio v -Closed-loop system control vector V -Lyapunov candidate function W -Compliance matrix W aa -Matrix of compliance (mechanical coupling) between actuators W ea -Matrix of compliance (mechanical coupling) between end-effector and Contents Introduction 1 Chapter 1 Introduction 1.1 xiv actuators Direction before cable guide e -System error E -Young's modulus E def -Energy of deformation F -F ext -Vector of external forces x -Vector of estimated states Vector of internal forces W sa -Matrix of compliance (mechanical coupling) between sensors and ac-tuators x -Contents Vector of generalized states ẋ -Velocity vector 1.1
	S	-	Set of triangles composing a cavity
	u	-	Displacement field
	u e	-	Element displacement field

E. P. Box vii viii 2 1.3 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.4 Chapter summary . . . . . . . . . . . . . . . . . . . . . 3 1.5 List of Publications . . . . . . . . . . . . . . . . . . . . 4

•

  Duriez, C., Bieze, T. (2017). Soft Robot Modeling, Simulation and Control in Real-Time. In Soft Robotics: Trends, Applications and Challenges (pp. 103-109). Springer International Publishing. Hage S., Sanz-Lopez M., Duriez C., Project page at https://softroboticstoolkit.com/fetch • Bieze, T. (2015). On the Kinematic Modeling Methodology of Soft Manipulator Robots. Poster presentation. In Robotics and Automation (ICRA), 2015 IEEE International Conference on. IEEE. Introduction . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Continuum manipulator definition . . . . . . . . . . .

	Other contributions • Runner up Prize at the Soft Robotics Toolkit Design Competition (2015). FeTCH Mark I Manipulator. Contributors: Bieze T., State of the Art Contents Largilliere F., Chapter 2 2.1

7 2.2.1 Soft continuum robot applications . . . . . . . . . . . 8 2.2.2 Bio-inspiration . . . . . . . . . . . . . . . . . . . . . . 9 2.2.3 Classification . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 State of the art in soft, continuum manipulators . . . 13 2.3.1 Design of continuum manipulators . . . . . . . . . . . 13 2.3.2 Modeling of continuum robots . . . . . . . . . . . . . . 18 2.3.3 Dynamics and control of continuum robots . . . . . . 27 2.4 Work contextualization and contributions . . . . . . . 29

Table 3 .

 3

	1: Error metrics for the 3 modeling approaches for the
	CBHA in free load configuration
		CC model HA model FEM model
	RMSD(mm)	12.87	17.09	4.66
	ME(mm)	16.35	22.4	7.71

Table 3

 3

	.2: Error between measured and estimated lengths for
	the CBHA
	l(mm)	l 1	l 2	l 3	l 4	l 5	l 6
	µ(mm) 3.2 2.43 3.86 4.08 3.6 3.69
	σ(mm) 1.55 1.76 2.05 2.12 2.56 2.06

Thèse de Thor Morales Bieze, Lille 1, 2017 © 2017 Tous droits réservés. lilliad.univ-lille.fr Thèse de Thor Morales Bieze, Lille 1, 2017 © 2017 Tous droits réservés. lilliad.univ-lille.fr

Thèse de Thor Morales Bieze, Lille 1, 2017 © 2017 Tous droits réservés. lilliad.univ-lille.fr

https://team.inria.fr/defrost/

https://www.cristal.univ-lille.fr/

http://www.festo-didactic.com/int-en/

http://conacyt.gob.mx/ Thèse de Thor Morales Bieze, Lille 1, 2017 © 2017 Tous droits réservés. lilliad.univ-lille.fr

In robotics, the degrees of freedom can be of two types: passive degrees of freedom are the ones that cannot be directly controlled and conform passively to the movements of the robot, while active degrees of freedom are the ones that are controlled to provide the required forces to move the manipulator.Thèse de Thor Morales Bieze, Lille 1,

© 2017 Tous droits réservés. lilliad.univ-lille.fr

Thèse de Thor Morales Bieze, Lille 1, 2017 © 2017 Tous droits réservés. lilliad.univ-lille.fr2.2. Continuum manipulator definition

Thèse de Thor Morales Bieze, Lille 1, 2017 © 2017 Tous droits réservés. lilliad.univ-lille.fr 2.3. State of the art in soft, continuum manipulators

In this work, we refer to the direct finite element when the unknown nodal displacement is computed given the known external forces Thèse de Thor Morales Bieze, Lille 1, 2017 © 2017 Tous droits réservés. lilliad.univ-lille.fr

Thèse de Thor Morales Bieze, Lille 1, 2017 © 2017 Tous droits réservés. lilliad.univ-lille.fr3.2. Continuum mechanics framework

This statement can be extended to any physically realizable system.Thèse de Thor Morales Bieze, Lille 1,

© 2017 Tous droits réservés. lilliad.univ-lille.fr

Thèse de Thor Morales Bieze, Lille 1, 2017 © 2017 Tous droits réservés. lilliad.univ-lille.fr3.3. Finite Element Method

Thèse de Thor Morales Bieze, Lille 1, 2017 © 2017 Tous droits réservés. lilliad.univ-lille.fr3.4. FEM-based kinematics of soft manipulators

Thèse de Thor Morales Bieze, Lille 1, 2017 © 2017 Tous droits réservés. lilliad.univ-lille.fr 3.5. Reduced model in the constraint space

Thèse de Thor Morales Bieze, Lille 1, 2017 © 2017 Tous droits réservés. lilliad.univ-lille.fr 4.3. Closed-loop control of continuum manipulators

In fact, FEM can be viewed as a piecewise variational method Thèse de Thor Morales Bieze, Lille 1,

© 2017 Tous droits réservés. lilliad.univ-lille.fr

Acknowledgments

First, I would like to dedicate this manuscript to my parents, Jessica Bieze and Jose Luis Morales, who never forget to send me their love every day, and to my three sisters Luisa, Jessica and Melissa.

I would like to thank my two supervisors, Rochdi Merzouki and Christian Duriez for believing in me. Thank you for your guidance and for sharing your knowledge and wisdom with me. I could not have better mentors and for that, I will be forever grateful.

This Ph.D. thesis was developed within the research groups, Deformable Robotic This Ph.D. was funded by the Mexican National Council of Science and Technology (CONACYT) 4 .

Appendix A

Introduction to Continuum

Mechanics and FEM

A.1 Introduction

The objective of this appendix is to give the reader a brief introduction to continuum mechanics and FEM. The interested reader should refer to [START_REF] Constantinescu | [END_REF]], [START_REF] Reddy | [END_REF]] for a more detailed explanation of the concepts presented in the sequel.

A.1.1 Continuum Mechanics

When analyzing the effects of forces and energy on bodies from a mechanical point of view, we can choose to perform the analysis at different physical scales. Continuum Mechanics performs this study at the macroscopic scale, by assuming that the body is a continuous mass. By making this assumption, the molecular substructures can be disregarded and the body is then considered a continuum. In the following, the description of motion in continuum bodies is presented.

A.1.1.1 Concept of configuration and description of motion

To arrive at a description of motion in a continuum media, let us consider a body B with all its material particles occupying a three dimensional domain Ω ∈ R. For a given set of loads (external forces), the body will deform changing its geometry and this deformation will be accompanied by internal body reaction forces, called stress. The new region κ i occupied by the body after the i-th deformation is called the i-th configuration of the continuum. In this case, Ω = κ 0 is the initial or underformed configuration and it is used as a reference to express the deformation of the continuum body.

where u i are the cartesian components of the displacement, E 11 , E 22 and E 33 are called normal strains and E 12 , E 23 and E 13 are called shear strains. The tensor can also be written in terms of the gradient of the displacement

In the case of infinitesimal strains, the displacement gradients are small and the nonlinear terms in Eq. A.8 can be neglected. The infinitesimal strain tensor becomes

Also important in the study of continuum motion is the measure of stress. As mentioned before, when a body deforms, the deformation causes reaction forces within the body. These internal forces are called stress and are defined as the force per unit area.

Consider a small element of area a in a continuum media. The stress vector m acting on the surface of the element is defined as

where Δf (n) is the force acting on a small surface of area Δa and n is the unit vector normal to that area. We can also write m in terms of the Cauchy stress tensor σ as

The Cauchy stress tensor σ is defined as the current force per unit of deformed area. The first Piola-Kirchhoff stress tensor or Lagrangian stress tensor

is defined as the current force per unit of undeformed area. In a similar way, we introduce the second Piola-Kirchhoff stress tensor S, defined as the transformed current force per unit of undeformed area. The second Piola-Kirchhoff stress tensor is associated to the first Piola-Kirchhoff stress tensor by the equation

In practice, the second Piola-Kirchhoff stress tensor is used in the study of large deformations, useful for the simulation of elastic materials.

A.1.1.3 Constitutive equations

In order to reconstruct the deformed configuration of a body under a set of external forces, a relationship between stress and strain is necessary.

In our study, we will continue to consider the case of small strains.

For linear elastic materials, the relationship between the components of stress and strain tensors at each point of the body is given by Hooke's law, which is written in tensor form as

where the scalar coefficients C ijkl are called stiffness coefficients of the stiffness tensor C. In the general case, C is formed by 81 coefficients, but due to the symmetry of the stress and strain tensors, the coefficients are reduced to 21. Using single subscript notation for stresses and strains and two subscripts for the stiffness coefficients

Eq. A.14 takes the form

which is refered as the Voigt-Kelvin notation [Dellinger 1998]. In matrix notation, Eq. A.16 is written as

The linear Tetrahedron

The simplest solid element is the four-node tetrahedron, also known as the linear tetrahedron [Felippa 2001], depicted in fig A.5. The tetrahedron geometry is defined by the position of its nodes in Cartesian coordinates x i , y i and z i with i = 1, 2, 3, 4 or in the tetrahedral natural coordinates denoted by The value of ζ i is one at corner i and zero at the other 3 corners, including the entire opposite face and it varies linearly with distance as one traverses the distance from the corner to that face.

Any function f (x, y, z) linear in x, y, z that takes values f i (i = 1, 2, 3, 4) at the corners of the tetrahedron may be interpolated in terms of the coordinates of the tetrahedron by

The interpolation functions are expressed in terms of the tetrahedron coordinates because they depend on the geometry of the element, however, other quantities like stress, strain and displacements are expressed in Cartesian coordinates. Therefore, a transformation between these two coordinate systems is needed. This transformation is given by 104 Appendix A. Introduction to Continuum Mechanics and FEM

Using the relationships for stress and strain, derived in A.2.2.1, we can express Eq. A.49 in terms of the displacements of the nodes

where the superindex e denotes a value at the element level. Nodal displacement u e which corresponds to the minimum of the functional Π are determined by the condition ∂Π ∂u e = 0 (A.51)

Differentiating Π with respect to u e produces the following equilibrium equations for the finite element

Eq. A.52 is usually represented in the equivalent form .53) with

where matrix K e is the element stiffness matrix and F e is the vector of loads. The aim of the assembly is to form the global equation system

using the element equations of the ith finite element

Let us introduce the matrix K d and vectors u d and F d formed by the concatenation of Eq. A.56, such as

with the relations

where matrix K is the global stiffness matrix. Here, matrix A provides the transformation from global to local node enumeration. In practice, one node can be shared by several elements. The nodes have a local index (at the element level) and a global index (at the level of the mesh). Finally, we arrive at the following global equation system

A.2.4 Displacement boundary conditions

Displacement boundary conditions were not present in the total potential energy functional presented in Eq. A.49. They are usually applied to the global equation system, given by Eq. A.55, after its assembly for convenience. In the general case, before the application of boundary conditions, matrix K is singular, which means that its rows and columns are linearly dependent and, therefore, K is not a full rank matrix.

Let us consider the application of the displacement boundary condition

to the global system Eq. A.55. We will use explicit method for this purpose.

In the explicit method, we substitute the known value of the displacement u e m = d in the mth column of matrix K, of dimensions nxn, and move this column to the right-hand side of Eq. A.55

Domain decomposition of continuum manipulators

In some continuum manipulator designs, the structure of the robot is composed by a sequence of serially connected substructures that have a vertebra-inter-vertebra pattern. In such designs, the inter-vertebra section shows a deformation that is greater than the vertebra under actuation. In the context of FEM modelization, the model dimensions of such manipulators can be significantly decreased by means of a domain decomposition technique.

B.1 Domain decomposition of the CBHA

The body of the CBHA is composed by a set of substructures formed by the aforementioned vertebra-inter-vertebra pattern. The main idea of the domain decomposition is to discretize the deformation of the structure using a set of frames with 6 DoF, attached to the structure at the position of the vertebrae. In the simulation of the robot, the frames are nodes with 6 DoF which are kinematically linked to the intervertebra section. A non-linear stiffness matrix that correspond to the inter-vertebra section is pre-computed and condensed as springs with equivalent stiffness between the frames.

To pre-compute the equivalent stiffness, a substructure of the manipulator is isolated. The substructure is composed by 2 consecutive frames and the inter-vertebra between them. The nodes at the upper and lower boundary of the substructure are attached to the frames and have 6 DoF (pure rigid motion) while the nodes in-between the vertebrae have only 3 DoF (deformation). These 2 types of nodes need to be merged in the same mechanical system. The stiffness matrix K of the inter-vertebra can be computed using a conventional FEM approach that takes into account the nodes kinematically linked to the rigid frames. Given the