Actuation force vector (Input) dimension n u " 3

End-effector (Output) dimension q " " q T 1 ¨¨¨q T N ‰ T P R n Generalized strain vector P R 3n Displacement vector τ " rτ 1 ¨¨¨τ nτ s T P R nτ Soft robot actuators force vector.

K " K P R nˆn Generalized stiffness matrix " Kpqq P R 3nˆ3n Generalized stiffness matrix Hpqq P R nˆnτ Generalized actuation matrix P R 3nˆnτ Generalized actuation matrix F pqq P R n Generalized external force vector P R 3n Generalized external force vector δτ Actuation vector with respect to prior configuration. In an etymological sense, and according to [START_REF] Zunt | classique de la conception des robots mous est associée à l'incertitude d'atteindre de tels objectifs et aux dépenses économiques substantielles nécessaires pour les essais et les erreurs. Dans le but d'atteindre les objectifs souhaités, il est donc logique[END_REF], the word robotics originated from the word robot, which is derived from the Czech noun robota meaning labor. Its original appearance was in the RUR play of the first fictional automatons called robots published in 1920 [START_REF] Zunt | classique de la conception des robots mous est associée à l'incertitude d'atteindre de tels objectifs et aux dépenses économiques substantielles nécessaires pour les essais et les erreurs. Dans le but d'atteindre les objectifs souhaités, il est donc logique[END_REF]. Scientifically, robotics is a multidisciplinary field that aims at the design, modeling, control, and applications of robots for the purpose of increasing productivity and assisting human beings [START_REF] Cherubini | Collaborative manufacturing with physical human-robot interaction[END_REF].

In accordance with [START_REF] George Thuruthel | Control strategies for soft robotic manipulators: A survey[END_REF][START_REF] Brian | Lie groups, Lie algebras, and representations: an elementary introduction[END_REF][START_REF] Trivedi | Soft robotics: Biological inspiration, state of the art, and future research[END_REF], robots can be classified based on their fabrication material and their degrees of freedom. Consistently, there are two main classes of robots (as shown by Fig. 1.1): rigid and hyper-redundant (it is also important to state that some rigid robots can also be hyper-redundant, e.g., serpentine robots).

Rigid Robotics and the Rise of a New Robot Class

Rigid robots are a class of robots made of rigid materials with a Young's modulus greater than 10 9 P a [START_REF] Majidi | Soft robotics: a perspective-current trends and prospects for the future[END_REF] The first appearance of a rigid robot can be traced back to the third century B.C. and earlier by the mechanical engineer Yan Shi that created the first human-shaped figure of his mechanical handiwork [START_REF] Needham | Science and civilisation in china: Volume 2, history of scientific thought[END_REF]. The first mechanical humanoid robot Mechanical Knight was designed by Leonardo da Vinci (see A in Fig. 1.2). Mechanical robots have then continued to emerge until 1738 where a mechanical duck called Digesting Duck that was able to eat and flap its wings was created by the inventor and artist Jacques de Vaucanson (see B in Fig. 1.2). Following the birth of electricity, robotics have evolved first via the Teleautomaton device that was created in 1898 by Nicola Tesla, which was a radio-controlled vessel (see C in Fig. 1.2), followed by the first electronic humanoid robot Elektro in 1930 and created by the Westinghouse Electric Corporation (see D in Fig. 1.2). The first appearance of rigid robots imitating biological behavior is traced back to 1948 by William Grey Walter [START_REF] Sabbatini | Imitation of life: A history of the first robots[END_REF] (see E in Fig. 1.2). Next, the first commercial robot Unimate was created in 1956 by the Unimation company [START_REF] Waurzyniak | Masters of manufacturing: Richard morley[END_REF], the same company that introduced the first installed industrial robot in 1961 (see F in Fig. 1.2). After that, the first full-scale humanoid robot WABOT-1 [START_REF] Kuo | Robotics and Mechatronics: Proceedings of the 6th IFToMM International Symposium on Robotics and Mechatronics (ISRM 2019)[END_REF] was created by the Waseda University (see G in Fig. 1.2). Then, the famous six electro-mechanically driven axis robot Famulus [START_REF] Singh | Evolution of industrial robots and their applications[END_REF] (see H in Fig. 1.2) was developed by the KUKA Robot Group in 1973. Next, in 1975, a programmable universal manipulation arm PUMA was created by Victor Scheinman (see I in Fig. 1.2). In 1983, the creation of the first parallel programming language used for the robot's control was performed by Stevo Bozinovski and Mihail Sestakov. [START_REF] Sabbatini | Imitation of life: A history of the first robots[END_REF]. F-Unimate [START_REF] Waurzyniak | Masters of manufacturing: Richard morley[END_REF]. G-WABOT-1 [START_REF] Kuo | Robotics and Mechatronics: Proceedings of the 6th IFToMM International Symposium on Robotics and Mechatronics (ISRM 2019)[END_REF]. H-Famulus [START_REF] Singh | Evolution of industrial robots and their applications[END_REF]. I-PUMA. J-Delta and ABB Flex-Pitcher robots [START_REF] Douglas R Ewing | Robots in the operating room-the history[END_REF][START_REF] Gasparetto | A brief history of industrial robotics in the 20th century[END_REF][START_REF] Hägele | Industrial robotics[END_REF].

Since then, researchers made a huge progress on rigid robotics, and owing to their stiff material and high-torque joints which make them exceptionally stable and accurate, rigid robots have become widely and successfully deployed especially in the industrial environment [START_REF] Sylla | Ergonomic contribution of able exoskeleton in automotive industry[END_REF] and have been proved to increase productivity (e.g., Delta and ABB Flex-Pitcher robots in J of Fig. 1.2) [START_REF] Douglas R Ewing | Robots in the operating room-the history[END_REF][START_REF] Gasparetto | A brief history of industrial robotics in the 20th century[END_REF][START_REF] Hägele | Industrial robotics[END_REF].

However, as new and more complex applications arise, traditional rigid robots have shown many drawbacks due to their high stiffness, lack of compliance in conventional actuation mechanisms, and limited degrees of freedom [START_REF] Fazeli | Fundamental limitations in performance and interpretability of common planar rigid-body contact models[END_REF]. Moreover, rigid robots have been considered less practical when performing operations in dynamic environments [START_REF] Cetinkunt | Performance limitations of joint variable-feedback controllers due to manipulator structural flexibility[END_REF] and not safe when interacting in fragile environments and also with humans [START_REF] Cherubini | Collaborative manufacturing with physical human-robot interaction[END_REF][START_REF] Holly | A taxonomy for human-robot interaction[END_REF].

In the last decades, as a consequence of these drawbacks, many researchers have been motivated to seek a novel design of robots that can be flexible and safe in order to cope with new applications, leading to the rise of soft robotics.

Soft Robotics

Hyper-redundant robots consist of a large number of DoFs [START_REF] Trivedi | Soft robotics: Biological inspiration, state of the art, and future research[END_REF]. Continuum robots are capable of continuous deformations, but not all continuum robots are necessarily soft [START_REF] Trivedi | Soft robotics: Biological inspiration, state of the art, and future research[END_REF]. Soft robots are a sub-class of continuum robots that are hyper-redundant robots [START_REF] Morales | Contribution to the kinematic modeling and control of soft manipulators using computational mechanics[END_REF].

In a general sense, soft robotics is a branch of robotics that consists of the design, modeling, and control of robots that are fabricated from flexible and compliant materials equivalent to those found in living organisms such as silicone rubber and other materials that can be readily found in nature [START_REF] Kim | Soft robotics: a bioinspired evolution in robotics[END_REF][START_REF] Laschi | Soft robotics: Technologies and systems pushing the boundaries of robot abilities[END_REF][START_REF] Rus | Design, fabrication and control of soft robots[END_REF][START_REF] Trivedi | Soft robotics: Biological inspiration, state of the art, and future research[END_REF].

The history of soft robotics is provided in detail in [START_REF] Hughes | Soft manipulators and grippers: a review[END_REF][START_REF] Laschi | Soft robotics: Technologies and systems pushing the boundaries of robot abilities[END_REF][START_REF] Surya G Nurzaman | Active sensing system with in situ adjustable sensor morphology[END_REF][START_REF] Rus | Design, fabrication and control of soft robots[END_REF][START_REF] Wang | Soft-material robotics[END_REF]. In the following, we present only the main points of the history of soft robotics. The first appearance of soft robots is traced back to the late 1970s where robot grippers based on granular materials were first published in [START_REF] Cardaun | Stand der greiferentwicklung, fordern und heben[END_REF][START_REF] Ap Perovskii | Universal grippers for industrial robots[END_REF][START_REF] Schmidt | Flexible moulding jaws for grippers[END_REF]. After, a continuously-deforming body using elastomers was controlled by a pneumatic actuation in 1984 [START_REF] James F Wilson | Robotic mechanics and animal morphology[END_REF]. Next, robot grippers using electrorheological fluid were proposed in 1989 [START_REF] Gary | Electrorheological fluid-based robotic fingers with tactile sensing[END_REF]. In 1991, [START_REF] Suzumori | Development of flexible microactuator and its applications to robotic mechanisms[END_REF] developed a flexible micro-actuator and later extended its applications to soft robotic grasping and also mobile mechanisms [START_REF] Suzumori | Elastic materials producing compliant robots[END_REF][START_REF] Suzumori | Applying a flexible microactuator to robotic mechanisms[END_REF] (see A in Fig. 1.3).

Ever since the field of soft robotics continued to make outstanding progress and has gained impressive attention at the beginning of the 21st century. Consequently, a huge portion of new soft material has been synthesized and made commercially accessible [START_REF] Bauer | 25th anniversary article: a soft future: from robots and sensor skin to energy harvesters[END_REF][START_REF] Hughes | Soft manipulators and grippers: a review[END_REF][START_REF] Kim | Soft robotics: a bioinspired evolution in robotics[END_REF][START_REF] Laschi | Soft robotics: new perspectives for robot bodyware and control[END_REF][START_REF] Laschi | Soft robotics: Technologies and systems pushing the boundaries of robot abilities[END_REF][START_REF] Surya G Nurzaman | Active sensing system with in situ adjustable sensor morphology[END_REF][START_REF] Rus | Design, fabrication and control of soft robots[END_REF][START_REF] Trivedi | Soft robotics: Biological inspiration, state of the art, and future research[END_REF], diverse fabrication techniques for soft robots have been developed and made available [START_REF] Bauer | 25th anniversary article: a soft future: from robots and sensor skin to energy harvesters[END_REF][START_REF] Hughes | Soft manipulators and grippers: a review[END_REF][START_REF] Kim | Soft robotics: a bioinspired evolution in robotics[END_REF][START_REF] Laschi | Soft robotics: new perspectives for robot bodyware and control[END_REF][START_REF] Laschi | Soft robotics: Technologies and systems pushing the boundaries of robot abilities[END_REF][START_REF] Surya G Nurzaman | Active sensing system with in situ adjustable sensor morphology[END_REF][START_REF] Rus | Design, fabrication and control of soft robots[END_REF][START_REF] Trivedi | Soft robotics: Biological inspiration, state of the art, and future research[END_REF], a growing number of works showing the use of soft robotics has been published in high-profile scientific journals [START_REF] Hughes | Soft manipulators and grippers: a review[END_REF][START_REF] Kim | Soft robotics: a bioinspired evolution in robotics[END_REF][START_REF] Laschi | Soft robotics: new perspectives for robot bodyware and control[END_REF][START_REF] Laschi | Soft robotics: Technologies and systems pushing the boundaries of robot abilities[END_REF][START_REF] Rus | Design, fabrication and control of soft robots[END_REF][START_REF] Trivedi | Soft robotics: Biological inspiration, state of the art, and future research[END_REF], and researchers generally agree that soft robots should be used in robotic applications in the future as they are safer and more adaptive in dynamic environments compared to the conventional rigid robots.

Being made of soft and deformable material, soft robots provide many benefits such as high adaptability and dexterity, and safe collisions [START_REF] Delettre | Robust control of a planar manipulator for flexible and contactless handling[END_REF], which are very useful for various applications, especially the manipulation of fragile objects [START_REF] Michael | Kinematics and the implementation of an elephant's trunk manipulator and other continuum style robots[END_REF][START_REF] Hughes | Soft manipulators and grippers: a review[END_REF][START_REF] Laschi | Soft robot arm inspired by the octopus[END_REF][START_REF] Mcmahan | Field trials and testing of the octarm continuum manipulator[END_REF] (see B, C, and D in Fig. 1.3), environment exploration [START_REF] Aracri | Soft robots for ocean exploration and offshore operations: A perspective[END_REF][START_REF] Della | Soft robots. Encyclopedia of Robotics[END_REF][START_REF] Giorgio-Serchi | Hybrid parameter identification of a multi-modal underwater soft robot[END_REF][START_REF] Mintchev | A soft robot for random exploration of terrestrial environments[END_REF][START_REF] Shan | Untethered soft robots for future planetary explorations? Advanced Intelligent Systems[END_REF][START_REF] Zhu | Tuna robotics: A high-frequency experimental platform exploring the performance space of swimming fishes[END_REF] (see E, F, G, and H in Fig. 1.3), and medical operations [START_REF] Burgner | A telerobotic system for transnasal surgery[END_REF][START_REF] Burgner-Kahrs | Continuum robots for medical applications: A survey[END_REF][START_REF] Cianchetti | Biomedical applications of soft robotics[END_REF][START_REF] Cianchetti | Soft robotics technologies to address shortcomings in today's minimally invasive surgery: the stiff-flop approach[END_REF][START_REF] Haouas | 4-dof spherical parallel wrist with embedded grasping capability for minimally invasive surgery[END_REF][START_REF] Mattheis | Flex robotic system in transoral robotic surgery: the first 40 patients[END_REF][START_REF] Mitros | Design and modelling of a continuum robot for distal lung sampling in mechanically ventilated patients in critical care[END_REF][START_REF] Simaan | A dexterous system for laryngeal surgery[END_REF][START_REF] Tsz | Soft robotics in medical applications[END_REF] (see I, J, K in Fig. 1.3). Figure 1.3: Soft Robots Applications: A-grasping and mobile mechanisms [START_REF] Suzumori | Elastic materials producing compliant robots[END_REF][START_REF] Suzumori | Applying a flexible microactuator to robotic mechanisms[END_REF]. B-octArm [START_REF] Mcmahan | Field trials and testing of the octarm continuum manipulator[END_REF]. C-elephant's trunk [START_REF] Michael | Kinematics and the implementation of an elephant's trunk manipulator and other continuum style robots[END_REF]. D-octopus arm [START_REF] Laschi | Soft robot arm inspired by the octopus[END_REF]. E-terrestrial exploration [START_REF] Mintchev | A soft robot for random exploration of terrestrial environments[END_REF]. F-tuna robotics [START_REF] Zhu | Tuna robotics: A high-frequency experimental platform exploring the performance space of swimming fishes[END_REF]. G-multi-modal underwater [START_REF] Giorgio-Serchi | Hybrid parameter identification of a multi-modal underwater soft robot[END_REF]. H-planetary explorations [START_REF] Shan | Untethered soft robots for future planetary explorations? Advanced Intelligent Systems[END_REF]. I-transoral robotic surgery [START_REF] Mattheis | Flex robotic system in transoral robotic surgery: the first 40 patients[END_REF]. J-laryngeal surgery [START_REF] Simaan | A dexterous system for laryngeal surgery[END_REF]. K-distal lung sampling [START_REF] Mitros | Design and modelling of a continuum robot for distal lung sampling in mechanically ventilated patients in critical care[END_REF].

Challenges and Motivations

Even though soft robots have made an interesting progress during the last few years, the actual process of designing a soft robot is still mainly inspired by a range of biological systems [START_REF] Rus | Design, fabrication and control of soft robots[END_REF] including the structure and behavior of animal species such as the trunk of elephants [START_REF] Michael | Kinematics and the implementation of an elephant's trunk manipulator and other continuum style robots[END_REF], and the arms of an octopus [START_REF] Laschi | Soft robot arm inspired by the octopus[END_REF]. Such a design procedure can be useful for initial tests and experimentation of the designed soft robot. However, when faced with specific performance objectives such as path planning and control [START_REF] Delettre | Robust control of a planar manipulator for flexible and contactless handling[END_REF][START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF], it is possible that such a soft robot's design may not be able to achieve its objectives, especially because its range of reachability (i.e., its workspace) may be restricted.

Therefore, it is of benefit to the soft robotics community to propose generic methodologies in order to evaluate the workspace of soft robots, and also to assist and guide the systematic design of soft robots for the purpose of optimizing specific performance objectives.

The workspace estimation in soft robotics is still an open subject, and its importance is due to the numerous benefits that it can provide in order to solve different soft robotics' scientific challenges, mainly related to the robot's mechanical design and controller synthesis. On the one hand, the result of workspace evaluation can provide information about the accessibility of the object to be controlled [START_REF] Delettre | Robust control of a planar manipulator for flexible and contactless handling[END_REF][START_REF] Della | Dynamic control of soft robots interacting with the environment[END_REF][START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF], i.e. by identifying if the object's position belongs to the soft robot's workspace, sparing the controller's design time that can be wasted in trying to control an unreachable object outside a robot's workspace. On the other hand, workspace information is crucial as well when planning feasible trajectory to be tracked by soft robots since the start and finish positions of a feasible trajectory should belong to a robot's workspace [START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF]. Besides, workspace analysis is also a necessary step to assist an effective, purpose-driven, and optimal design of soft robots [START_REF] Hiller | Automatic design and manufacture of soft robots[END_REF], improving their stability, and increasing the range of their manipulability and reachability.

To achieve this aim, the actual thesis adopts two different mathematical models to describe the deformation of soft robots both with a slender shape and general form. Based on those adopted mathematical models, several approaches to estimate the workspace of soft robots are proposed. Finally, an effective method to optimize the design of soft robots for the purpose of achieving specific performance objectives is presented.

State of the Art

This work is concerned with the workspace evaluation and design optimization of soft robots based on their adopted mathematical models. Accordingly, a concise state of the art targeting the modeling, workspace evaluation, and design optimization of soft robots will be presented.

Mathematical Modeling of Soft Robots

Owing to their natural conformity and compliance, soft robots consists of a large number of DoFs and are characterized by different mechanical laws [START_REF] Morton E Gurtin | An introduction to continuum mechanics[END_REF] compared to that of rigid robots, which makes their kinematic and dynamic modeling highly nonlinear and more complex.

In the literature, several modeling methods were proposed to describe the deformation of soft robots, and the most important ones are the Piece-wise Constant Curvature (PCC) [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF], the Finite Element Method (FEM) [START_REF] Duriez | Control of elastic soft robots based on real-time finite element method[END_REF], and the Piece-wise Constant Strain (PCS) based Discrete Cosserat method [START_REF] Renda | Discrete cosserat approach for multisection soft manipulator dynamics[END_REF].

In the following, we will briefly investigate each modeling method.

Piece-wise Constant Curvature

The PCC approach [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF] was initially established for kinematics studies and later extended to the dynamics of soft robots. It describes the investigated soft robot via a fixed number of arcs characterized by the curvature's radius, the arc's angle, and the bending plan. The robot is kinematically described through n s successive transformations T j i (as shown in Fig. 1.4), with n s being the total number of segments that constitutes the investigated robot.

The PCC has been widely and successfully applied to several continuum robots controlled by various methods of actuation [START_REF] Della | On an improved state parametrization for soft robots with piecewise constant curvature and its use in model based control[END_REF][START_REF] Robert K Katzschmann | Dynamic motion control of multi-segment soft robots using piecewise constant curvature matched with an augmented rigid body model[END_REF][START_REF] Runge | A framework for the kinematic modeling of soft material robots combining finite element analysis and piecewise constant curvature kinematics[END_REF][START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF] . However, when the studied soft robot is subject to external loads, the assumption of a constant curvature does not always hold. Figure 1.4: PCC modeling approach [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF].

Finite Element Method

In [START_REF] Duriez | Control of elastic soft robots based on real-time finite element method[END_REF], a modeling approach based on FEM was developed to describe the deformation of soft robots with general geometries (as shown by Fig. 1.5). In this method, the deformable domain of the structure is discretized into smaller finite elements using a specific meshing geometry, composed of a finite number of DoFs, which is then used to interpolate the behavior of the deformable domain by measuring the variation of the associated nodal displacements.

This modeling technique was successfully applied to different geometries of soft robots controlled by different means of actuation [START_REF] Morales Bieze | Finite element method-based kinematics and closed-loop control of soft, continuum manipulators[END_REF][START_REF] Duriez | Control of elastic soft robots based on real-time finite element method[END_REF][START_REF] Goury | Fast, generic, and reliable control and simulation of soft robots using model order reduction[END_REF][START_REF] Rodríguez | Real-time simulation of hydraulic components for interactive control of soft robots[END_REF]. However, to obtain a good modeling precision, this approach requires a fine mesh when discretizing the robot's spatial space, which in turn increases the number of DoFs and results in high computation complexity.

To solve this problem, model-order reduction (MOR) techniques have been introduced into the FEM approach [START_REF] Goury | Fast, generic, and reliable control and simulation of soft robots using model order reduction[END_REF] in order to reduce the number of DoFs and eventually decrease the computation complexity, but this improvement slightly affects the modeling precision. Figure 1.5: FEM modeling approach [START_REF] Duriez | Control of elastic soft robots based on real-time finite element method[END_REF][START_REF] Goury | Fast, generic, and reliable control and simulation of soft robots using model order reduction[END_REF].

Piece-wise Constant Strain

The Discrete Cosserat method [START_REF] Renda | Discrete cosserat approach for multisection soft manipulator dynamics[END_REF][START_REF] Renda | Discrete cosserat approach for soft robot dynamics: A new piece-wise constant strain model with torsion and shears[END_REF][START_REF] Renda | A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation[END_REF][START_REF] Renda | Discrete cosserat approach for multisection soft manipulator dynamics[END_REF][START_REF] Renda | Discrete cosserat approach for soft robot dynamics: A new piece-wise constant strain model with torsion and shears[END_REF][START_REF] Renda | Dynamic model of a multibending soft robot arm driven by cables[END_REF] is a Piece-wise Constant Strain (PCS) method that can be viewed as an extension of the PCC technique.

This approach (see Fig. 1.6) describes the designed soft robot with a finite set of stain vectors (linear and angular) and provides the benefit of achieving a higher modeling precision with a lower number of DoFs. Besides, this modeling technique considers the material properties and geometric nonlinearity, thus taking into account large deformations. However, like the PCC, this approach can only be applied to slender-shaped soft robots. Moreover, the PCS assumes that all cross-sections (which are perpendicular to the centroidal line of slender soft robots) are rigid and undeformable. Such an assumption is imposed for the sake of modeling simplicity, but inevitably it may affect the modeling precision. To solve this problem, a generalization of the PCS method was recently investigated in [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF].

Figure 1.6: PCS modeling approach [START_REF] Renda | Discrete cosserat approach for multisection soft manipulator dynamics[END_REF][START_REF] Renda | Discrete cosserat approach for soft robot dynamics: A new piece-wise constant strain model with torsion and shears[END_REF].

Modeling Methods Comparison

The PCC modeling approach provides an effective and simple method for the modeling of soft robots. This method was largely deployed in the community of soft robotics, however, the adopted assumptions are not always valid which is very important for modeling the dynamics of soft robots.

Alternatively, the PCS method takes into account both the material properties of the studied robot and geometric nonlinearity, but with the limitation that it can only be applied to slendershaped soft robots. Accordingly, the PCS modeling method is adopted by the present thesis in order to achieve the mathematical model of soft robots with slender forms.

On the other hand, and although it is computationally expensive, the FEM modeling method is adopted by this thesis in order to provide the mathematical model of soft robots with general form.

Workspace Evaluation

By definition, the workspace (also known as the accessible output set) of a robot is an area/volume that contains the whole of its achievable equilibrium configurations.

Rigid Robotics Workspace Methodologies

The subject of workspace determination has been widely investigated by the community of rigid robotics. In the literature, we can find many contributions to this topic using different methods to represent the workspace of rigid serial/parallel manipulators.

As to analytical Jacobian-based methods, the method proposed by [START_REF] Abdel-Malek | On the determination of boundaries to manipulator workspaces[END_REF][START_REF] Abdel | Analytical boundary of the workspace for general3dof mechanisms[END_REF] (see A in Fig. 1.7) uses an explicit input-output function and consists of analytically determining singularities that make the sub-Jacobian row-rank deficient in order to map the boundaries of rigid manipulators. This method was indeed successful in determining the singularities for a wide variety of rigid mechanisms [START_REF] Khairallah | Serial 5dof manipulators: workspace, void, and volume determination[END_REF]. However, this method has its limitations [START_REF] Abdel-Malek | On the determination of boundaries to manipulator workspaces[END_REF]: it requires that the output is an explicit function of the input (i.e., x " Ωpqq), thus making it impossible to be used for soft robots. Moreover, a higher dimension of inputs increases the complexity of the analytical formulation of the sub-Jacobian which as a result makes it near impossible to analytically determine all singularities of such mechanisms.

In [START_REF] Edward | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF][START_REF] Ej Haug | Dextrous workspaces of manipulators. i. analytical criteria[END_REF] a Jacobian-based numerical method was presented for mapping the boundaries of rigid manipulator's workspaces using a continuation method based on the theory of bifurcation (see B in Fig. 1.7). To demonstrate the broadness of the applicability of this method, it was used to determine the workspace of planar and spatial Stewart platforms [START_REF] Masory | Workspace evaluation of stewart platforms[END_REF]. However, while seeking an initial boundary point, starting from a feasible initial point inside the workspace, this method may fail to identify enclosure curves (namely voids/holes) inside the workspace due to the discretization [START_REF] Abdel-Malek | On the determination of boundaries to manipulator workspaces[END_REF]. Moreover, it may also not be possible to distinguish an internal and an external boundary to the workspace when the ray emanating from an initial configuration passes through a bifurcation point [START_REF] Abdel-Malek | On the determination of boundaries to manipulator workspaces[END_REF]. These problems were properly handled in [START_REF] Bohigas | A complete method for workspace boundary determination on general structure manipulators[END_REF] where the methods presented in [START_REF] Edward | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF] were extended to treat more complex workspaces. Moreover, continuation methods were also used in [START_REF] Boehler | Definition and computation of tensegrity mechanism workspace[END_REF] for the computation of tensegrity mechanism workspace, and a general high-order continuation method [START_REF] Gauthier Hentz | Higher-order continuation for the determination of robot workspace boundaries[END_REF] with automatic differentiation was applied to a planar RRR mechanism and a three-dimensional Orthoglide parallel mechanism.

Regarding numerical methods, an optimization approach to determine the boundaries of rigid manipulator's workspaces has been presented in [START_REF] Ja Snyman | An optimization approach to the determination of the boundaries of manipulator workspaces[END_REF] (see C in Fig. 1.7). This numerical method consists of finding a suitable radiating point in the output coordinate space and then determining the points of intersection of a representative pencil of rays, emanating from the radiating point, with the boundary of the accessible set. This method is based on an optimization approach in which the accessibility of the workspace boundary is formulated in terms of a constrained optimization problem. However, the difficulties of this method lie in defining an appropriate objective function, non-convexity issues, and the inability to identify interior knowledge of the workspace.

Also, the wrench-feasible workspace of parallel rigid robots driven by cables has been analyzed in [START_REF] Gouttefarde | Interval-analysis-based determination of the wrench-feasible workspace of parallel cable-driven robots[END_REF] by applying an interval analysis bisection method. Such a method has been then extended to analyze the workspace of cable-driven parallel mechanisms in [START_REF] Berti | Workspace analysis of redundant cable-suspended parallel robots[END_REF][START_REF] Merlet | On the workspace of suspended cable-driven parallel robots[END_REF] (see D in Fig. 1.7). The advantage of this approach lies in the fact that it explores all feasible configurations in the workspace. However, this approach is in a way exhaustive since it consists of bisecting the whole space in order to find the feasible configurations which represent the workspace.

Finally, a numerical probabilistic method, Monte Carlo, was presented in [START_REF] Alciatore | Determining manipulator workspace boundaries using the monte carlo method and least squares segmentation[END_REF][START_REF] Cao | Accurate numerical methods for computing 2d and 3d robot workspace[END_REF] to find the 2D and 3D workspace of a 3-DOF serial manipulator (see E in Fig. 1.7). This numerical method yields a robust result since it explores all actuators combinations, and has the advantage that it does not rely on the inverse problem of the robot. However, this technique is computationally expensive for high DoFs mechanisms. Figure 1.7: Rigid workspace methods: A-Jacobian-based analytical methods [START_REF] Abdel | Analytical boundary of the workspace for general3dof mechanisms[END_REF]. B-Jacobianbased numerical methods [START_REF] Edward | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF]. C-optimization methods [START_REF] Ja Snyman | An optimization approach to the determination of the boundaries of manipulator workspaces[END_REF]. D-interval analysis methods [START_REF] Gouttefarde | Interval-analysis-based determination of the wrench-feasible workspace of parallel cable-driven robots[END_REF]. E-Monte-Carlo methods [START_REF] Cao | Accurate numerical methods for computing 2d and 3d robot workspace[END_REF].

Although the above methods have been successfully implemented in determining the workspace of several types of rigid robots. However, they cannot be directly extended to analyze soft robots because of the modeling complexity for soft robots. In fact, the input-output relation for open-loop rigid chains (also called serial-link rigid robots) can normally be explicitly defined through its direct kinematic model (u " f pqq). However, such a relation for soft robots is generally defined in an implicit way (f pu, τ q " 0) through a complex static/kinematic model (which will be explained in Chapter 2) for the purpose of modeling soft material deformations. On the other hand, closed-loop rigid chains (also called parallel rigid robots) necessitate the use of numerical schemes in order to solve the direct kinematics of the robot [START_REF] Vigen Arakelian | Increase of singularity-free zones in the workspace of parallel manipulators using mechanisms of variable structure[END_REF][START_REF] Briot | Are parallel robots more accurate than serial robots[END_REF][START_REF] Briot | Dynamics of parallel robots[END_REF][START_REF] Merlet | Parallel robots[END_REF], which is the same problem for soft robots. Hence, the workspace methodologies explored for the different kinds of rigid mechanisms need to be adapted in order to treat the same problem for soft robots since they consist of a large number of DoFs due to the compliance and the flexibility of the used material, and also because the modeling of soft robots is established through different mechanical laws, i.e., principles of continuum mechanics [START_REF] Morton E Gurtin | An introduction to continuum mechanics[END_REF] (which will be explained in Section 2.5.2 for the PCS model, and Section 2.5.3 for the FEM model).

Continuum Robotics Workspace Approaches

On the other side, in the community of continuum robotics, the workspace subject has been rarely investigated. To our best knowledge, three main contributions were provided and all of them consist of discretizing the inputs (actuators) space in different fashions.

The paper of [START_REF] Burgner-Kahrs | Workspace characterization for concentric tube continuum robots[END_REF] (see A in Fig. 1.8) presents the workspace characterization for concentric tube continuum robots by extending Monte-Carlo algorithms developed in [START_REF] Cao | Accurate numerical methods for computing 2d and 3d robot workspace[END_REF]. Also, in [START_REF] Yuan | Workspace analysis of cable-driven continuum manipulators based on static model[END_REF] (see B in Fig. 1.8), the workspace analysis of cable-driven continuum manipulators based on the static model was investigated using an approach that consists of discretizing the input domain and finding the set of all reachable equilibrium positions of the manipulator's end-effector.

The limitations of the above-stated approaches are: firstly, they consist of discretizing the inputs space, which varies depending on the configuration of the investigated continuum robot; secondly, the computational complexity (which is equal to n nτ d ) of these approaches can explode exponentially when increasing the dimension of the actuators n τ or when we want to get higher precision, i.e., increasing n d ; thirdly, valuable information about the robot's singularities and interior/exterior boundaries cannot be identified using these classical approaches since they only focus on solving the direct model of the investigated continuum robots in order to find the set of poses that constitutes the workspace.

In [START_REF] Cao | Workspace analysis of tendon-driven continuum robots based on mechanical interference identification[END_REF] (see C in Fig. 1.8), an approximation method was used to determine the workspace boundaries of continuum robots by sweeping the robot from the neutral configuration to the fully curved configuration for the purpose of deducing the workspace surface from the robot's tip trajectory. The proposed approximation method consists of estimating the workspace boundary by searching only the points distributed on the surface of the workspace. However, this method still consists of exploring the inputs space in a different manner by varying one configuration from its minimum to its maximum value while keeping the other configurations at a fixed neutral value. This strategy can map the workspace boundaries of the investigated robot and at the same time reduce the computational time of the pre-proposed scanning method from n nτ d combinations to 2 nτ combinations. Nevertheless, for robots controlled by a higher number of actuators, the approximation method will still be computationally expensive. Another limitation of this approach lies in failing to identify some workspace boundaries because singularities result not only from actuators saturation, but also due to physical limitations [START_REF] Bohigas | A complete method for workspace boundary determination on general structure manipulators[END_REF][START_REF] Edward | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF] such as the mechanism's length, self-collision, or actuation instabilities, because this approach doesn't apply the idea that the workspace interior/exterior boundaries of the investigated continuum robots are extracted from the set of their output singularities ignoring the row-rank deficiency of the Jacobian of the generalized mechanical constraints which doesn't comply with what was studied in [START_REF] Abdel | Analytical boundary of the workspace for general3dof mechanisms[END_REF][START_REF] Bohigas | A complete method for workspace boundary determination on general structure manipulators[END_REF][START_REF] Edward | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF].

Soft Robotics Workspace Methods

For a given soft robot's configuration, till the starting date (1 October 2018) of the present thesis, no work has yet been investigated about how to comprehensively evaluate its workspace.

Workspace Methods Summary

As stated above, different methodologies were investigated for the workspace evaluation of rigid robots, few for continuum robots, and almost none for soft robots.

The methods investigated for continuum robots are not efficient since they are based on the discretization of the inputs space, whose dimension depends on the studied configuration of the robot, but such methods can still be used as a comparative reference since they exhaustively explore all possible achievable configurations of the investigated robot's workspace.

On the other hand, it is more appropriate to extend the applicability of the main workspace approaches developed in the rigid robotics community in order to treat the soft robotics case. However, this cannot be done in a direct way since soft robots are modeled in a more complex manner, yielding an implicit input-output relation. Figure 1.8: Continuum workspace methods: A-Monte-Carlo algorithms [START_REF] Burgner-Kahrs | Workspace characterization for concentric tube continuum robots[END_REF]. B-direct method [START_REF] Yuan | Workspace analysis of cable-driven continuum manipulators based on static model[END_REF]. C-approximation method [START_REF] Cao | Workspace analysis of tendon-driven continuum robots based on mechanical interference identification[END_REF].

Design Optimization of Soft Robots

Several design optimization methods were proposed for specific soft robotics applications, which can be classified into three main categories: heuristic methods, model-based, and hybrid.

Heuristic Methods

Genetic algorithms were explored in [START_REF] Chapelle | Evaluation functions synthesis for optimal design of hyper-redundant robotic systems[END_REF][START_REF] Hiller | Automatic design and manufacture of soft robots[END_REF] for the design optimization of soft robots (see A in Fig. 1.9). These methods are particularly useful for treating highly nonlinear problems since they don't require the calculation of costly gradients. However, defining a proper fitness function is difficult.

Evolutionary algorithms were also used in [START_REF] Daniel M Bodily | Multi-objective design optimization of a soft, pneumatic robot[END_REF] (see B in Fig. 1.9) where it was proposed to optimize a small number of size parameters of pneumatic-actuated inflatable manipulators, and in [START_REF] Runge | Design optimization of soft pneumatic actuators using genetic algorithms[END_REF] (see C in Fig. 1.9) for the design optimization of soft pneumatic actuators. In [START_REF] Ansari | Multiobjective optimization for stiffness and position control in a soft robot arm module[END_REF] (see D in Fig. 1.9), the central concept was to develop an assistive manipulator that can automate the bathing task for elderly citizens where the development of a novel algorithm based on cooperative multi-agent reinforcement learning that simultaneously optimizes stiffness and position was proposed. Figure 1.9: Heuristic Methods: A-Automatic design of soft robots [START_REF] Hiller | Automatic design and manufacture of soft robots[END_REF]. B-Multi-objective design optimization of a soft pneumatic robot [START_REF] Daniel M Bodily | Multi-objective design optimization of a soft, pneumatic robot[END_REF]. C-Design optimization of soft pneumatic actuators [START_REF] Runge | Design optimization of soft pneumatic actuators using genetic algorithms[END_REF]. D-Multi-objective optimization for stiffness and position control in a soft robot arm [START_REF] Ansari | Multiobjective optimization for stiffness and position control in a soft robot arm module[END_REF].

Model-Based Approaches

As to model-based approaches, in [START_REF] Zhang | Topology optimized design, fabrication and evaluation of a multimaterial soft gripper[END_REF] (see A in Fig. 1.10), a density-based topology optimization method was employed for the design of the exterior shell of a soft pneumatic gripper with a given shape and position in order to grasp a wide range of objects using the FEM model. Next, improving the mechanical performance of soft pneumatic robots through parameterized layout and shape optimization of air chambers was studied in [START_REF] Guo | Chamber layout design optimization of soft pneumatic robots[END_REF] (see B in Fig. 1.10). In [START_REF] Nalbach | Modeling and design optimization of a rotational soft robotic system driven by double cone dielectric elastomer actuators[END_REF] (see C in Fig. 1.10), a model-based optimization for the design of a rotational soft robotic system driven by double cone dielectric elastomer actuators was investigated. The work of [START_REF] Tommaso Carella | Design optimization of a passive soft robotic device for neurorehabilitation[END_REF] (see D in Fig. 1.10) studied the design optimization of a soft robotic device for neurorehabilitation by computing a set of optimal parameters which, in turn, could be used to customize the device according to the patient's need. In [START_REF] Elsayed | Finite element analysis and design optimization of a pneumatically actuating silicone module for robotic surgery applications[END_REF] (see E in Fig. 1.10), the design of a pneumatically actuated silicone module with three pneumatic chambers is considered and optimized with the aim of using it in a soft robot arm for robotic surgery applications. In [START_REF] Md Ghazaly | Design optimization & analysis of a soft crawling robot[END_REF] (see F in Fig. 1.10), the design optimization of soft crawling robots (inchworm and quadrupedal) for maximum bending motion with respect to input pressure was evaluated using FEM. The work of [START_REF] Steven E Mchugh | Dynamics-based design of a soft robot[END_REF] (see G in Fig. 1.10) developed a methodology for converting caterpillar-like soft-bodied robot from a lumped-parameter form into a Computer-Aided Design (CAD) model that could be easily manufactured using the results of coarse dynamic simulations. Also, the subject of optimal actuator location was investigated in [START_REF] Morris | Linear-quadratic optimal actuator location[END_REF], where it is proposed to minimize a linear-quadratic (LQ) function to find the optimal actuators location. This method was then extended in [START_REF] Wu | Optimal actuator location for electro-active polymer actuated endoscope[END_REF] (see H in Fig. 1.10) to achieve the optimal actuator location for a medical endoscope controlled by electro-active polymers. Figure 1.10: Model-based approaches: A-topology optimized design of a multimaterial soft gripper [START_REF] Zhang | Topology optimized design, fabrication and evaluation of a multimaterial soft gripper[END_REF]. B-soft pneumatic robots chamber layout design optimization [START_REF] Guo | Chamber layout design optimization of soft pneumatic robots[END_REF]. C-rotational soft robotic system design [START_REF] Nalbach | Modeling and design optimization of a rotational soft robotic system driven by double cone dielectric elastomer actuators[END_REF]. D-soft robotic device design optimization [START_REF] Tommaso Carella | Design optimization of a passive soft robotic device for neurorehabilitation[END_REF]. E-design optimization of a pneumatically actuating silicone module [START_REF] Elsayed | Finite element analysis and design optimization of a pneumatically actuating silicone module for robotic surgery applications[END_REF]. F-design optimization of a soft crawling robot [START_REF] Md Ghazaly | Design optimization & analysis of a soft crawling robot[END_REF]. G-dynamics-based design of soft robots [START_REF] Steven E Mchugh | Dynamics-based design of a soft robot[END_REF]. H-optimal actuator location for electro-active polymer actuated endoscope [START_REF] Wu | Optimal actuator location for electro-active polymer actuated endoscope[END_REF].

Hybrid Approaches

In [START_REF] Raeisinezhad | Design optimization of a pneumatic soft robotic actuator using model-based optimization and deep reinforcement learning[END_REF] (see Fig. 1.11), two frameworks, a model-based optimization, and deep reinforcement learning were proposed for the design optimization of a multi-chamber pneumatic-driven soft actuator in order to optimize its mechanical performance. The design goal was to achieve maximal horizontal motion of the top surface of the actuator with a minimum effect on its vertical motion.

Conclusion

As stated above, different strategies (heuristic, model-based, and hybrid) were fully investigated for the design optimization of different soft robots. Those methods have helped achieve the optimal design of the investigated soft robots which were successfully deployed in real-life scenarios.

However, despite the significant progress in the conceptual design of soft robots, the mentioned methods are in a sense specific to a particular type of robot with a predefined concept. On the other hand, generic methods focusing on the systematic design optimization of soft robots that are controlled via a generic number of actuators for the purpose of optimizing specific performance objectives have rarely been addressed in the literature. This is due to the difficulty and complexity that lies in establishing an appropriate cost function combining multiple desired design objectives, particularly when these objectives are competing with each other.

To our best knowledge, no former work has been investigated on the optimization of the design of soft robots for the purpose of optimizing their workspace attainability and reachability which is of high importance and benefit to the design optimization of soft robots.

Contributions of the Thesis 1.3.1 Contributions and Organization of the Manuscript

For soft robotic applications such as design optimization, path planning, and pick-place, it is necessary and useful to evaluate the workspace of the studied soft robot. In this thesis, we investigate such a problem from two views, a direct and an inverse one.

On the one hand, the direct problem of the workspace can be formulated by the following question: given the designed configuration of the investigated soft robot, how should we efficiently estimate its workspace? (which will be studied in Chapter 3, Chapter 4, and Chapter 5). Estimating the workspace of a given soft robot's prototype is a crucial step for assessing the reachability and accessibility of particular locations which relates to the controller synthesis of the robot. Also, solving this problem helps in determining the robot's singularities and impediments of motion which relates to the design of the robot.

On the other hand, the following question formulates the inverse problem: how should we optimize the design of a soft robot in order to achieve specific performance objectives? (which will be investigated in Chapter 6). Given specific performance objectives (e.g., reaching specific locations), it is valuable for many reasons (both economic and scientific) to determine the soft robot's optimal design such that its workspace contains those desired locations in a virtual environment before proceeding to the physical design of the final prototype of the robot.

Accordingly, to answer these two questions, the manuscript was organized into 7 chapters as follows (See Fig. 1.12): Figure 1.12: Manuscript organization.

Chapter 2 introduces the dynamic model of the adopted PCS and FEM modeling methods. Also, in this chapter, we establish the workspace definition of soft robots. Finally, we present a classic approach (named the forward approach) that consists of discretizing the actuators space in order to provide the workspace of soft robots. This approach was then validated via a trunk-like soft robot for both the PCS and FEM cases.

In Chapter 3, we propose an optimization-based approach [START_REF] Amehri | Fem-based exterior workspace boundary estimation for soft robots via optimization[END_REF][START_REF] Walid | Discrete cosserat method for soft manipulators workspace estimation: An optimization-based approach[END_REF] that consists of estimating the workspace of soft robots by mapping their workspace's exterior boundaries. This approach was applied to both the PCS [START_REF] Walid | Discrete cosserat method for soft manipulators workspace estimation: An optimization-based approach[END_REF] and FEM [START_REF] Amehri | Fem-based exterior workspace boundary estimation for soft robots via optimization[END_REF] modeling methods.

Next, an interval analysis-based method [START_REF] Amehri | Fem based workspace estimation for soft robots: a forward-backward interval analysis approach[END_REF] (named the forward-backward approach) is outlined in Chapter 4 in order to achieve the workspace of soft robots by exploring all feasible configurations. The proposed approach was applied to both the PCS and FEM modeling methods.

Chapter 5 presents a continuation method [START_REF] Amehri | Workspace boundary estimation for soft manipulators using a continuation approach[END_REF] based on the theory of bifurcation in order to map the interior and exterior boundaries of soft robots. This approach was only applied to the PCS modeling method.

After, an optimization approach based on the adopted PCS model was established in Chapter 6 in order to optimize the design of soft robots for the purpose of achieving specific performance objectives.

Finally, Chapter 7 provides the conclusions of the present thesis and presents the perspective of future works, especially, the orientation aspect of the workspace, the generalization of the continuation approach, and the possibility of its extension to treat the FEM case, and also the design optimization of soft robots based on the FEM method.

Publications

The contributions of the present thesis are outlined in the following: International Conference:

1. Amehri, Walid, Gang Zheng, and Alexandre Kruszewski. "Fem based workspace estimation for soft robots: A forward-backward interval analysis approach." In 2020 3rd IEEE International Conference on Soft Robotics (RoboSoft), pp. 170-175. IEEE, 2020. Summary: In this paper, the preliminary result of an interval analysis-based approach in order to estimate the workspace of soft robots using the FEM modeling method was submitted.

Journal Articles: Summary: In this journal, the full version of the proposed interval analysis approach was presented which contains new results for FEM model with external forces, self-contained explanation and results on interval analysis, novel uniform spatial grid discretization strategy, a detailed algorithm with an improved stop condition, and comprehensive configuration simulation scenarios to emphasize the effectiveness of the proposed interval analysis based approach.

Submitted Articles:

1. Amehri, Walid, Gang Zheng, and Alexandre Kruszewski. ASME JMR, Position-Access Workspace of Soft Manipulators using an Interval Analysis Method. Summary: In this article, we extended the applicability of the interval analysis-based approach by applying it to the PCS modeling method in order to estimate the workspace of soft robots with slender shape.

2. Amehri, Walid, Gang Zheng, and Alexandre Kruszewski. IJRR, Cosserat-based Optimization Design for Slender Soft Manipulators. Summary: In this article, we proposed an optimization approach based on the adopted PCS model in order to optimize the design of soft robots with slender shape for the purpose of achieving specific performance objectives.

Chapter 2

Modeling and Forward Methods

Introduction

As stated in Chapter 1, the PCS method provides a good modeling precision with a lower number of DoFs which can only be applied to slender-shaped soft robots. On the other hand, the FEM method provides a precise model of soft robots with the general form, but it is computationally expensive.

Accordingly, the first step of this chapter is to present the mathematical models of the PCS and FEM methods for the purpose of modeling soft robots with a slender shape and general form, respectively. Next, we establish the definition of a soft robot's workspace. Finally, we propose a forward approach to estimate the workspace of soft robots. This approach was then validated using a trunk-like soft robot [START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF] for both the PCS and FEM cases.

In the following, we do not consider the material's nonlinearity, instead, we consider that the material is isotropic, and the deformation is entirely elastic.

Piece-wise Constant Strain (PCS)

The dynamic model based on the Piece-wise Constant Strain (PCS) approach was developed in [START_REF] Renda | A geometric and unified approach for modeling soft-rigid multi-body systems with lumped and distributed degrees of freedom[END_REF] for a soft-rigid multibody system. In the following, all the quantities are expressed in the local (body) coordinate frame if not specified. The superscript 1 represents the partial derivative with respect to the space variable.

In the following, we illustrate the mathematical model of the PCS model using the example of a tendon-driven soft robot. The method could be as well applied to other manners of actuation [START_REF] Renda | A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation[END_REF].

Kinematics

The configuration of a soft body i with respect to its predecessor in the chain (as illustrated by Fig. 2.1) is defined as a curve g i p¨q : X P r0, L i s Þ Ñ g i pXq P SEp3q with:

g i pXq " ˆRi u i 0 T 1 ˙(2.1a)
where R i is the orientation matrix, u i is the end-effector position, L i is the length, ξi is the strain matrix, of a body i, and X P r0, Ls is the abscissa along the robot arm.

The continuous models of the position, velocity and acceleration of a soft body can be derived from the Cosserat rod theory, which gives [START_REF] Boyer | Poincare's equations for cosserat media: Application to shells[END_REF]:

g 1 i " g i ξi (2.2a) η 1 i " 9 ξ i ´ad ξi η i (2.2b) 9 η 1 i " : ξ i ´ad 9 ξi η i ´ad ξi 9 η i (2.2c)
where ad p¨q P R 6ˆ6 is the adjoint operator of the Lie algebra, and ξi pXq "

ˆk i p i 0 T 0 ˙P sep3q ξ i pXq " `kT i , p i T ˘T P R 6 (2.3)
defines the strain state with ki pXq P sop3q, k i pXq P R 3 , and p i pXq P R 3 respectively the angular and linear strain (when compared to the reference values k i and p i );

ηi pXq "

ˆw i v i 0 T 0 ˙P sep3q η i pXq " `wT i , v T i ˘T P R 6 . (2.4)
is the cross-section velocity twist with wi pXq P sop3q, w i pXq P R 3 , and v i pXq P R 3 respectively the angular and linear velocity. To model constrained rod, such as the Kirchhoff-Love case with angular strain only, the strain field is specified as:

ξ i " B qi q i `ξi (2.5)
where B qi P R 6ˆni forms a basis for the allowed motion subspace, q i P R ni contains the values of the allowed strains and, ξ i P R 6 is the initial twist modeling the initial shape.

Assuming piece-wise constant strains [START_REF] Renda | Discrete cosserat approach for soft robot dynamics: A new piece-wise constant strain model with torsion and shears[END_REF], (2.2) can be analytically integrated using the matrix exponential method, leading to:

g i pXq "e X ξi
(2.6a)

η i pXq "Ad ´1 gi η h `Ad ´1 gi T gi B qi 9 q i (2.6b) 9 η i pXq "Ad ´1 gi 9 η h `Ad ´1 gi ż X 0
Ad gipsq ad ηipsq dsB qi 9 q i `Ad ´1 gi T gi B qi : q i (2.6c)

where Ad gi pXq P R 6ˆ6 is the Adjoint operator of SEp3q, and T gi pXq is the tangent operator of the exponential map. Successive applications of the kinematics (2.6) for all the bodies of the system yields the definition of the geometric Jacobian J i pq, Xq P R 6ˆn and its derivative 9 J i pq, 9 q, Xq P R 6ˆn (n being the total number of DOFs), which links the generalized strain vector q " " q T 1 q T 2 ¨¨¨q T N ‰ T P R n (N being the total number of bodies) and the velocity twist η i pXq, for each soft body i, as shown below:

η i pXq " i ÿ h"0 Ad ´1 g h ¨¨¨g i Tg h Bq h 9 q h " i ÿ h"0 i S h 9 q h " J i pq, Xq 9 q
(2.7a)

9 η i pXq " i ÿ h"0 i S h : q h `Ad ´1 g h ¨¨¨g i ż X 0 Ad g h psq ad η h psq dsBq h 9 q h " i ÿ h"0 i S h : q h `i 9 S h 9
q h " J i pq, Xq: q `9 J i pq, 9 q, Xq 9 q

(2.7b) ηi pXq "g ´1 i 9 g i (2.7c)
where the block elements of the i th Jacobian i S p¨q P R 6ˆn p¨q and its derivative i 9 S p¨q P R 6ˆn p¨q have been defined in [START_REF] Renda | Discrete cosserat approach for multisection soft manipulator dynamics[END_REF].

Continuous Dynamics

Once a Jacobian is found, the generalized dynamics of the system can be obtained by projecting the free dynamics of each soft body by virtue of the D'Alembert's principle. The free dynamic equation of a soft body, with its boundary conditions, is given by [START_REF] Renda | Discrete cosserat approach for multisection soft manipulator dynamics[END_REF]:

M i 9 η i `ad ηi M i η i " `Fi i ´Fa i ˘1 `ad ξi `Fi i ´Fa i ˘`F e i `Fi i ´Fa i ˘p0q " ´FJ i `Fi i ´Fa i ˘pL i q " ´Ad gij F J j ; (2.8)
where M i pXq " diagpJ xi , J yi , J zi , A i , A i , A i qρ i P R 6ˆ6 is the screw inertia matrix of the crosssection (J ¨i pXq being the second moment of the area about the axis ¨and A i pXq being the area of the cross-section); Fei pXq P R 6 is the distributed external load; F ai pXq P R 6 is the internal wrench due to the distributed actuation [START_REF] Renda | Screw-based modeling of soft manipulators with tendon and fluidic actuation[END_REF]; F ii pXq P R 6 is the internal wrench due to the elasticity of the soft body; F J p¨q P R 6 is the wrench transmitted across joint p¨q and ad p¨q (respectively Ad p¨q ) P R 6ˆ6 is the co-adjoint (respectively co-Adjoint) map of the Lie algebra (respectively Lie group).

Elasticity and Actuation Forces

Regarding the internal elastic force, a linear viscoelastic constitutive model is usually chosen:

F ii pXq " Σ i pξ i ´ξ˚q `Υi 9 ξ i " Σ i B qi q i `Υi B qi 9 q i (2.9)
where

Σ i pXq " diagpG i J xi , E i J yi , E i J zi , E i A i , G i A i , G i A i q Υ i pXq " diagpJ xi , 3J yi , 3J zi , 3A i , A i , A i qν i P R 6ˆ6
are the screw stiffness and viscosity matrix (E i being the Young's modulus, G i the shear modulus and ν i the shear viscosity).

With regards to the actuation load (as depicted in Fig 2 .2), by computing the force and moment exerted by the internal tendon on the mid-line of the rod [START_REF] Renda | Screw-based modeling of soft manipulators with tendon and fluidic actuation[END_REF]. This is given by:

F a pXq " nτ ÿ i"1 " di pXqλ i pXq λ i pXq  τ i " H τ pq, Xqτ (2.10)
where d i pXq P R 3 represents the distance from the mid-line to the internal i th actuator, di denotes d i in sop3q, λ i pXq P R 3 is the unit vector tangent to the actuator path, and τ P R nτ is the vector of magnitude of the actuators force given by the negative of the tendons tension (n τ being the total number of actuators). While the distance d i pXq is fixed by design, the tangent vector λ i pXq depends, in general, on the deformation ξ i pXq of the soft robot body [START_REF] Rucker | Statics and dynamics of continuum robots with general tendon routing and external loading[END_REF], [START_REF] Renda | Dynamic model of a multibending soft robot arm driven by cables[END_REF]. Going further into details, the unit tangent vector λ i pXq can be obtained by spatial differentiation of the position vector of the actuator, thus: .11) where d i pXq is expressed in homogeneous coordinates, r¨s 3 extracts the first three rows of a homogeneous vector and ¨ takes the Euclidean norm.

λ i pXq " " g ´1 i pg i d i q 1 ‰ 3 g ´1 i pg i d i q 1 " " ξi d i pXq `d1 i pXq ı 3 ξi d i pXq `d1 i pXq . ( 2
Remark 1. Considering d i as a small parameter, a useful approximation can be obtained, which allows for a unit tangent vector independent of the deformation with a slight loss of accuracy [START_REF] Renda | Discrete cosserat approach for multisection soft manipulator dynamics[END_REF]. In this case, one can consider λ i " r1 0 0s T . As a result, the generalized actuation matrix Hpqq may be approximated by Hpqq « H.

Discrete Dynamics

By Jacobian projection of the free dynamics (2.8), we obtain the generalized dynamics in its classical form:

M pqq : q `pC pq, 9 qq `Dq 9 q `Kq " H pqq τ `F pqq (2.12)

where M P R nˆn is the generalized mass matrix, C P R nˆn is the generalized Coriolis matrix, D P R nˆn is the block-diagonal generalized damping matrix, K P R nˆn is the block-diagonal generalized stiffness matrix, H P R nˆnτ is the generalized actuation matrix, F P R n is the vector of generalized position-dependent external force and τ P R nτ is the vector of applied actuators force. Going further into details, the coefficient matrices take the form:

M pqq " N ÿ i"1 ż Li 0 J T i M i J i dX (2.13a) C pq, 9 qq " N ÿ i"1 ż Li 0 J T i ´ad Ji 9 q M i J i `Mi 9 J i ¯dX (2.13b) D " diag ˜BT q1 ż L1 0 Υ 1 dXB q1 , ¨¨¨, B T q N ż L N 0 Υ N dXB q N ¸(2.13c) K " diag ˜BT q1 ż L1 0 Σ 1 dXB q1 , ¨¨¨, B T q N ż L N 0 Σ N dXB q N ¸(2.13d) H pqq " » -˜BT q1 ż L1 0 H τ dX ¸T , ¨¨¨, ˜BT q N ż L N 0 H τ dX ¸T fi fl T (2.13e) F pqq " N ÿ i"1 ż Li 0 J T i Fei (2.13f) 2.3 Finite Element Method (FEM) 2.3.

FEM Introduction

The dynamic model of soft robots can be obtained based on Galerkin method of FEM and the Euler Lagrange equation [START_REF] Singiresu | The finite element method in engineering[END_REF]. FEM consists firstly of spatially discretizing the deformable domain of the structure into smaller finite elements N through the creation of a corresponding mesh which contains a finite number of DoFs n to interpolate the behavior of the deformable domain by measuring the variation of the associated nodal displacements. Then, the Euler-Lagrange method is applied to deduce the dynamics of soft robots.

In the literature, various geometries are proposed to discretize the spatial domain of a given structure, such as prisms, beams, hexahedrons and tetrahedrons [96, p. 732]. For the purpose of generality and convenience for element matrix integration, we select the linear tetrahedron element (shown by Fig. 2.3) as the appropriate meshing structure in this thesis. Based on the discretized mesh model, the deduction of its dynamics can be realized element by element. In the following, we apply the modeling approach on a single linear tetrahedron element to get its motion equation, and then we perform an assembly process of all elements to obtain the motion equation for the whole structure.

Displacement Model

After having created the mesh, we have to define the primary variables of a single element, which corresponds to the displacement of the nodes. According to [96, p. 93], the displacement vector can be defined by the following formula: δ e px, y, z, tq " N e q e (2. [START_REF] Berti | Workspace analysis of redundant cable-suspended parallel robots[END_REF] where the index e denotes the investigated element, N e is the basis functions matrix, and q e represents the nodal (1, 2, 3, 4) displacement vector for a single element, with:

N e " » - N 1 0 0 ¨¨¨N 4 0 0 0 N 1 0 ¨¨¨0 N 4 0 0 0 N 1 ¨¨¨0 0 N 4 fi fl e (2.

15)

And q e " rq 1x q 1y q 1z q 2x q 2y q 2z q 3x q 3y q 3z q 4x q 4y q 4z s T e (2.16)

The basis functions N e represents the mathematical model of the displacement over space, which needs to satisfy some requirements such as the continuity over a specific order and the Kronecker property where each basis function has a value of one at its own node and zero at the other nodes in the corresponding element [95, p. 82].

In the literature, we can find different polynomial families to define the basis functions. In the following, we present the Lagrange family [START_REF] Reddy | An introduction to the finite element method[END_REF] because it is accurate and straightforward. It has been stated in [71, p. 252] that the order p of the polynomial is defined via the number of nodes of the element n e by the following equation:

n e " 1 6
pp `1qpp `2qpp `3q " 4 (2.17)

Using the linear tetrahedron element which contains four nodes, the polynomial order for the basis functions should be equal to 1. According to this, the basis functions can be written in the following polynomial form:

N i px, y, zq " " α p1,iq `αp2,iq x `αp3,iq y `αp4,iq z ‰ e (2.

18)

where α i for 1 ď i ď 4 are the coefficients to be determined. For this, we define V equ as the volume of the element in an equilibrium state, then the basis functions of the linear tetrahedron element can be formulated as:

N i px, y, zq " " 1 6V equ pa i `bi x `ci y `di zq  e (2.19) with α p1,iq " a i 6V equ , α p2,iq " b i 6V equ , α p3,iq " c i 6V equ , α p4,iq " d i 6V equ
, and

V equ " 1 6 1 x 1 y 1 z 1 1 x 2 y 2 z 2 1 x 3 y 3 z 3 1 x 4 y 4 z 4 e
where px i , y i , z i q are the coordinates of the node i in the corresponding element. The variables a i , b i , c i and d i in (2.19) are constant coefficients which are computed when the element is in equilibrium point, with the following form:

a i " p´1q m`1 x j y j z j x k y k z k x l y l z l e , b i " p´1q m 1 y j z j 1 y k z k 1 y l z l e c i " p´1q m x j 1 z j x k 1 z k x l 1 z l e , d i " p´1q m x j y j 1 x k y k 1 x l y l 1 e
where m " j `k `l, and i, j, k, l P t1, 2, 3, 4u.

Strain-Displacement Model

The linear strain-displacement relation is defined by the following form [START_REF] Singiresu | The finite element method in engineering[END_REF]:

ε e " B e q e (2.20)
where ε e is the strain vector and B e is the strain matrix calculated by:

B e " » - - - - - - - - - - - - - - - - - - - - - - - - - - BN 1 Bx 0 0 ¨¨¨B N 4 Bx 0 0 0 BN
B e " 1 6V equ » - - - - - - - b 1 0 0 ¨¨¨b 4 0 0 0 c 1 0 ¨¨¨0 c 4 0 0 0 d 1 ¨¨¨0 0 d 4 c 1 b 1 0 ¨¨¨c 4 b 4 0 0 d 1 c 1 ¨¨¨0 d 4 c 4 d 1 0 b 1 ¨¨¨d 4 0 b 4 fi ffi ffi ffi ffi ffi ffi fl e (2.22)

Stress-Displacement Model

The relation between the stress and the nodal displacement vector is defined using the constitutive relation by the following representation [START_REF] Singiresu | The finite element method in engineering[END_REF]:

σ e " CB e q e
(2. [START_REF] Burgner | A telerobotic system for transnasal surgery[END_REF] where σ e represents the stress vector.

For an isotropic material with given Young's Modulus E and Poisson's ratio ν, the elasticity matrix is defined by: With the deduced displacement relation, strain-displacement relation and stress-displacement relation, we are ready to apply the Euler-Lagrange method in order to derive the motion equation of the corresponding element.

C " » - - - - - - -

Euler-Lagrange Equation

The Euler-Lagrange method is based on energy formulation, and it allows us to find the motion equation of the corresponding element of the soft robot's mesh. In general, its form can be written as follows [START_REF] Singiresu | The finite element method in engineering[END_REF]:

d dt " BL B 9 q e  `BE D B 9 q e " BL Bq e (2.25)
where E D is the dissipation energy, and L is the Lagrangian defined by [START_REF] Singiresu | The finite element method in engineering[END_REF] Lpq e , 9 q e , tq " E K ´EP (2.26) with E K and E P representing the kinetic energy and the potential energy, respectively.

Potential Energy

The potential energy is defined by the difference between the strain energy and the potential work [START_REF] Singiresu | The finite element method in engineering[END_REF]:

E P " E S ´WP (2.27)

Strain Energy

It is known that the strain energy is the energy stored in the body under loading, and it can be modeled as follows [START_REF] Singiresu | The finite element method in engineering[END_REF]:

E S " 1 2
q T e K e pq e qq e (2.28)

where K e represents the stiffness matrix for the analyzed element.

The integral form of the strain energy has the following form [START_REF] Singiresu | The finite element method in engineering[END_REF]: where f b e and f s e are the unit body and surface forces, H e pq e q is the actuation matrix, τ is the applied actuators force, and F e pq e q is the external applied forces, for the corresponding element.

E S " 1

Kinetic Energy

The kinetic energy is the energy of the body resulting from its motion, it can be modeled as follows [START_REF] Singiresu | The finite element method in engineering[END_REF]:

E K " 1 2 9
q e T M e pq e q 9 q e (2.34)

where M e represents the mass matrix for the related element.

The integral form of this kinetic energy has the following form [START_REF] Singiresu | The finite element method in engineering[END_REF]: where ρ e is the volumetric mass density.

E K " 1

Dissipation Energy

The dissipation energy is the energy due to the phenomena of damping, it can be modeled as follows [START_REF] Singiresu | The finite element method in engineering[END_REF]:

E D " 1 2 9 q e T D e pqq 9 q e (2.38)
where D e represents the damping matrix for the analyzed element.

The integral form of this damping energy has the following form [START_REF] Singiresu | The finite element method in engineering[END_REF]: Lpq e , 9 q e , tq " 1 2 9 q e T M e pq e q 9 q e ´1 2 q T e K e pq e qq e `rH e pq e qτ `Fe pq e qs q e (2.42)

E D " 1
Since K e is a symmetric matrix, thus the right term of (2.25) becomes:

BL Bq e " ´Ke pq e qq e `He pq e qτ `Fe pq e q (2.43)

The terms on the left of (2.25) depends only on 9 q, since the mass matrix M e and the damping matrix D e are symmetric, therefore we obtain:

d dt " BL B 9 q e
 " M e pq e q : q e BE D B 9 q e " D e pq e q 9 q e (2.44)

Finally, the motion equation of the analyzed element can be written as follows:

M e pq e q : q e `De pq e q 9 q e `Ke pqqq " H e pq e qτ `Fe pq e q (2.45)

Similarly, we apply the same procedure for each element of the mesh model, then we can assemble all elements using the connectivity matrix [START_REF] Reddy | An introduction to the finite element method[END_REF] to obtain the motion equation for the whole structure:

M pqq: q `Dpqq 9 q `Kpqqq " Hpqqτ `F pqq (2.46) where M pqq is the global mass matrix, Kpqq is the global stiffness matrix and Dpqq is the global damping matrix, and formulated below:

M pqq " N ÿ e"

Workspace Definition

The workspace of a soft robot contains all its achievable equilibrium configurations where the end-effector can reach and stay there. In this thesis, we focus on the position aspect of the workspace, and we discuss the orientation aspect of the workspace in Chapter 7, Section 7.2.1.

For the purpose of defining the workspace of soft robots, let us note u P R nu as the position of the end-effector in the inertial frame, and q P R n (n being the total number of DoFs) has a different meaning according to the appropriate modeling approach (it denotes the generalized strain vector in PCS, and the displacement vector in FEM). Then, we can define its workspace as follows:

Definition 1. Consider a soft robot dynamically modeled by (2.12) for the PCS case or (2.46) for the FEM case, controlled by bounded actuators τ P T , the workspace W E of the end-effector u, is a subspace of R nu , defined below: W E " tu " ϕpqq | Dpq P R n , τ P T q s.t. Φpτ, qq " 0u where ϕ represents the geometric model of the studied soft robot, Φpτ, qq " Kq ´Hpqqτ ´F pqq for PCS, and Φpτ, qq " Kpqqq ´Hpqqτ ´F pqq for FEM.

Remark 2. Using the above definition, It is important to state that the workspace contains all equilibrium positions where the soft robot's end-effector can reach and stay there. Moreover, the stability of a calculated configuration in the workspace can be asserted by verifying the eigenvalues of the stiffness matrix [START_REF] Zheng | Controllability pre-verification of silicone soft robots based on finite-element method[END_REF]. Finally, the elastic limit can be asserted by verifying that the corresponding constraint Φpτ, qq " 0 is valid for the corresponding τ P T .

In this thesis, we adopt two strategies to estimate the workspace of soft robots;

• Discretizing the inputs (actuators) space T [START_REF] Alciatore | Determining manipulator workspace boundaries using the monte carlo method and least squares segmentation[END_REF][START_REF] Burgner-Kahrs | Workspace characterization for concentric tube continuum robots[END_REF][START_REF] Cao | Workspace analysis of tendon-driven continuum robots based on mechanical interference identification[END_REF][START_REF] Cao | Accurate numerical methods for computing 2d and 3d robot workspace[END_REF][START_REF] Yuan | Workspace analysis of cable-driven continuum manipulators based on static model[END_REF].

• Discretizing the outputs (end-effector) space [START_REF] Abdel | Analytical boundary of the workspace for general3dof mechanisms[END_REF][START_REF] Gouttefarde | Interval-analysis-based determination of the wrench-feasible workspace of parallel cable-driven robots[END_REF][START_REF] Edward | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF][START_REF] Merlet | Parallel robots[END_REF][START_REF] Ja Snyman | An optimization approach to the determination of the boundaries of manipulator workspaces[END_REF].

The first strategy is manifested through the so-called forward approach [START_REF] Alciatore | Determining manipulator workspace boundaries using the monte carlo method and least squares segmentation[END_REF][START_REF] Burgner-Kahrs | Workspace characterization for concentric tube continuum robots[END_REF][START_REF] Cao | Accurate numerical methods for computing 2d and 3d robot workspace[END_REF][START_REF] Yuan | Workspace analysis of cable-driven continuum manipulators based on static model[END_REF] (see Section 2.5), and the second strategy is manifested through three different approaches; optimization [START_REF] Walid | Discrete cosserat method for soft manipulators workspace estimation: An optimization-based approach[END_REF] (see Chapter 3), interval analysis (see Chapter 4) [START_REF] Amehri | Fem based workspace estimation for soft robots: a forward-backward interval analysis approach[END_REF], and continuation [START_REF] Amehri | Workspace boundary estimation for soft manipulators using a continuation approach[END_REF] (see Chapter 5), which will be established in what follows.

Forward Approach 2.5.1 Forward Approach Introduction

The most basic approach to determine the workspace of soft robots is to iteratively approximate the input-output map φ : T Ñ W E through the discretization of the actuators bound T with a prescribed precision n d , and then determine for each actuator value its corresponding end-effector position, as illustrated by Fig. 2.5.

Figure 2.5: Forward Approach

In the following, we will explain the details of the forward approach for both the PCS and the FEM modeling methods.

Forward Approach for Workspace Estimation: PCS case

According to Definition 1, the workspace of a soft robot contains all equilibrium points. Applying this property to equation (2.12), the workspace of the studied soft robot is determined by the following equations:

K q " Hpqqτ `F pqq (2.48a) upqq " E 1 gpqq E T 2 (2.48b) with gpqq " e L1 ξ1 ¨¨¨e L N ξN (2.49a) ξi " z B qi q i ; i P 1 ¨¨¨N (2.49b)
The purpose of E 1 and E 2 in (2.48b) is to extract the end-effector position u from the configuration matrix g, where E 2 is a constant elementary vector defined as E 2 " r0 1ˆ3 1s, and E 1 is a constant elementary matrix defined by the following: when n u " 2, it may take the form of E 1 " rI 2 , 0 2 s, or E 1 " rr1, 0, 0, 0s T , r0, 0, 1, 0s T s T (depending on the actuators location); when n u " 3, E 1 " rI 3 0 3ˆ1 s (with I and 0 are respectively the identity and zero matrices).

Using the PCS framework, the problem of estimating W E is equivalent to, for all admissible τ P T , measure all the possible configurations of its end-effector u defined in (2.48b). However, this is not trivial, mainly due to the reason that the mapping φ : T Ñ W E is implicit and nonlinear, thus no analytic φ can be found such that u " φpτ q.

In practice, for the purpose of estimating W E from the bounded actuators τ P T , we might use numerical methods to approximate the mapping φ. In the following, we may use the gradient of the generalized strain vector q with respect to the actuators force vector τ , denoted as ∇ τ pqq, which will be analyzed hereafter.

By definition, the gradient ∇ τ pqq is calculated by:

∇ τ pqq " " Bq Bτ ı T
Thus, by calculating the partial derivative of (2.48a) with respect to the actuators vector τ , we have:

K Bq Bτ "
BrHpqqτ s Bτ `BrF pqqs Bτ Using the principle of variable separation, the above equation is equivalent to:

K Bq Bτ " BrHpqqs Bq pI n b τ q Bq Bτ `Hpqq `BrF pqqs Bq Bq Bτ
From this we get:

˜K ´BrHpqqs Bq pI n b τ q ´BrF pqqs Bq ¸Bq Bτ " Bpqq
where the symbol b is the Kronecker tensor product. Let us define the matrix B as follows for the sake of simplicity: Bpτ, qq " K ´Hq pqqpI n b τ q ´Fq pqq (2.50)

where H q pqq " BrHpqqs Bq , F q pqq " BrF pqqs Bq

Obviously, the calculation of ∇ τ pqq depends on the invertibility of the matrix B in (2.50), which in turns is indirectly determined by the generalized actuation matrix H in (2.10), and the vector of generalized position-dependent external force F in (2.13f).

Therefore, the following assumption has been imposed in this thesis.

Assumption 1. For the studied soft robot, it is supposed that the actuators are installed in such a way ensuring that the matrix B defined in (2.50) is invertible for all q P W E .

Remark 3. Physically, the above assumption means that the actuators of the investigated soft robots are arranged in a way that a small actuators variation for a feasible configuration q P W E might locally and uniquely determine the relative variation of position. In other words, it implies that the mounted actuators are installed such that they are not mounted on each other. If such an assumption is violated, i.e., several actuators are mounted in a redundant way, then it is necessary to keep only 1 actuator and remove other ones.

Hence, if Assumption 1 is satisfied, we can then obtain:

Bq Bτ " ˜K ´Hq pqqpI n b τ q ´BrF pqqs Bq ¸´1 Hpqq which means ∇ τ pqq " H T pqq " B ´1pτ, qq ‰ T .
As stated earlier in the introduction of this chapter (Section 2.5.1), the straightforward approach to estimate the workspace via (2.48) is to firstly discretize the actuators bounds T (as indicated in Fig. 2.5), and then approximate the mapping φ in an iterative manner.

Precisely, after discretizing the actuators force bounds T , noted as T d , then for each τ P T d , with the knowledge of the end-effector position upq pj´1q q and the generalized strain vector q pj´1q (for the case j " 1, upq p0q q, q p0q represents the end-effector's initial position and the initial generalized strain vector of the studied soft robot respectively), the following proposition enables us to calculate the next generalized strain vector q pjq . Proposition 1. If Assumption 1 is satisfied, then for a given τ P T d the following dynamics: q pjq " q pj´1q ´B´1 pτ, q pj´1q q Φpτ, q pj´1q q (2.51)

exponentially converges to the corresponding solution of q in (2.48a), with B defined in (2.50) and Φpτ, q pj´1q q " Kq pj´1q ´Hpq pj´1q qτ ´F pq pj´1q q.

Remark 4. It is important to state that convergence towards the closest solution with reference to the initial solution can only be guaranteed if the conditions resulting from Kantorovich's theorem [START_REF] James | The newton-kantorovich theorem[END_REF] are respected.

Proof. Solving (2.48a) numerically can be achieved by minimizing the following quadratic cost function Spτ, q ˚q: Spτ, q ˚q " Φ T pτ, q ˚qΦpτ, q ˚q, where Φpτ, q ˚q " Kq ˚´H pq ˚qτ ´F pq ˚q First, we compute the time derivative of the cost function: Spτ, q ˚q " 2Φ T pτ, q ˚q BrΦpτ, q ˚qs Bq ˚9 q " 2Φ T pτ, q ˚qBpτ, q ˚q 9 q ẘith B defined in (2.50). In order to minimize the cost function Spτ, q ˚q, we need to define an adaptive law 9 q ˚such that 9

Spτ, q ˚q is negative. For this, and in accordance with Assumption 1, if the following adaptive law is chosen:

9 q ˚" ´γB ´1pτ, q ˚qΦpτ, q ˚q (2.52)
where γ is a positive parameter, then we can obtain:

9
Spτ, q ˚q " ´2γSpτ, q ˚q which implies that Spτ, q ˚q exponentially converges to 0. Equivalently, it means that q exponentially converges to q with the proposed continuous-time adaptive law (2.52). Finally, the discrete-time version of (2.52) can be written as: q pjq " q pj´1q ´∆tγB ´1pτ, q pj´1q q Φpτ, q pj´1q q and we prove Proposition 1 by choosing γ " 1 ∆t .

After, with q pjq calculated via (2.51), we compute the strain twist ξpjq i related to the strain vector q pjq i of each body i as follows:

ξpjq i " { B qi q pjq i ; i P 1 ¨¨¨N
with which the configuration matrix gpq pjq q associated to the generalized strain vector q pjq can be calculated as:

gpq pjq q " e L1 ξpjq 1 ¨¨¨e L N ξpjq N Then, we can derive the corresponding end-effector position as follows:

upq pjq q " E 1 gpq pjq q E T 2
Finally the estimation of W E can be successfully accomplished by iterating the above procedure for all feasible inputs τ P T d from the initial state j " 1 to the last iteration j " n d .

Forward Approach for Workspace Estimation: FEM case

In accordance with the workspace definition (see Definition 1), the workspace of a soft robot can finally be defined by the following set of equations:

Kpqq q " Hpqqτ `F pqq (2.53a) q E " Cq (2.53b) u " q E `up0q (2.53c)
where q E is the displacement vector of the end-effector with respect to its known initial position vector u p0q , u represents the position vector of the end-effector, and C P R 3ˆ3n is a selection matrix associated with the end-effector node coordinates.

To achieve the straightforward approach, given the information of the end-effector's position vector u pj´1q and the vector of nodal displacements q pj´1q (for the case j " 1, u p0q and q p0q are respectively the initial position of the end-effector and the initial nodal displacement vectors of the investigated soft robot), we apply Newton-Raphson method to obtain the next nodal displacement q pjq corresponding to its actuators value τ P T d , with T d being the discretized set of T , which can be realized as follows:

q pjq " q pj´1q ´Φq pu, q pj´1q q Φpτ, q pj´1q q, @τ P T d where Φ q pτ, qq is the gradient of Φpτ, qq with respect to q, and Φpτ, q pj´1q q " Kpq pj´1q qq pj´1q Hpq pj´1q qτ ´F pq pj´1q q. Then, the corresponding end-effector position vector can be deduced by: u pjq " Cq pjq `up0q

Finally, as depicted by Fig. 2.5, the estimation of the reachable workspace W E can be achieved by iterating the above process for all τ P T d .

Experimental Validation

The goal of this section is to validate the workspace estimation obtained by applying the forward approach on the PCS and FEM models with the real workspace obtained from the experimental setup using a trunk-like soft robot [START_REF] Wu | Optimal actuator location for electro-active polymer actuated endoscope[END_REF] (as depicted by Fig. 2.6) for both the PCS and FEM models.

The exact geometric parameters of the investigated trunk-like soft robot are illustrated in Fig. 2.6. The studied soft robot is controlled by 4 independent symmetric tendons (2 upward and 2 lateral) mounted through it from the base to the tip, and the tendons are actuated by 4 independent stepper motors. To obtain the value of the end-effector position, a sensor position was positioned on the tip of the robot, and for this purpose, a long and uniform hole was made along the whole length of the robot. Moreover, 3D-printed rigid rings (blue rings in Fig. 2.6) are mounted along the soft robot to minimize friction between the tendons and the soft rubber. Finally, the casting material of the investigated trunk-like robot is an isotropic silicone rubber with its Young's modulus E " 6.5 ˆ10 5 P a and a Poisson's ration ν " 0.45. Using the forward approach, the real workspace can be achieved by reading the value of the position of the end-effector for each value of τ . The actuators bounds were within a magnitude T " r0, 5s 1 ˆ¨¨¨ˆr0, 5s 4 N , and the obtained experimentation result of the workspace is depicted by gray-colored points in Fig. 2.7.

On the other hand, in a simulation environment, we apply both the PCS and FEM methods to model the investigated trunk-like soft robot. The studied soft robot was modeled by the PCS method using N " 15 cylindrical bodies with equal lengths and decreasing radius from the base to the tip of the robot. Using the FEM method with N " 5024 tetrahedral elements and n " 1484 nodes, the investigated soft robot was modeled. The tendons in the simulation environment were approximately positioned using their real coordinates. Based on the adopted mathematical models, we apply the forward approach to get the corresponding workspace estimation. The obtained simulation results of the workspace for both the PCS and FEM cases are depicted by red-colored points in Figs. 2.7a and 2.7b, respectively. The result of the workspace obtained from the real experimentation took approximately 360 seconds, while the results of the workspace obtained via the simulations took approximately 50 seconds. Visually, we can observe from Fig . 2.7, the workspace estimation results estimated by applying the forward approach on the PCS and FEM case approximately matches the workspace result obtained from the real experimentation. Numerically, we compare the volume of the set of points obtained from the real experiment V r with the set of points obtained via the PCS (noted for the accuracy of the FEM result. Based on this, and in the following chapters, the workspace obtained via the forward approach will be considered as a reference. Moreover, the workspace estimation achieved by the approaches that will be presented in the next chapters (Chapter 3, Chapter 4, and Chapter 5) will be superimposed and compared with the workspace estimation obtained from the forward approach in order to highlight their validity and efficiency.

V P CS

Conclusion

In this chapter, we have presented the mathematical models of two adopted modeling methods, PCS and FEM. Next, we have presented the definition of the workspace of a soft robot. Finally, we have proposed a forward approach to estimate the workspace of soft robots based on PCS and FEM, and we have used a trunk-like soft robot to validate the workspace estimation obtained via the forward approach for the PCS and FEM with the real workspace result.

However, the main disadvantage of the forward approach lies in its computation complexity which can be very high.

Given n d subsets (n d P N), since τ P T Ă R nτ (n τ P N), for each entry of τ we will obtain n nτ d possible combinations of input for τ . Based on this perspective, the computation complexity of the forward approach depends on two parameters, n d and n τ . To obtain a correct estimation of W E , the discretization precision n d needs to be sufficiently big, which will result in a significant increase in the computation complexity. Moreover, it will be exponentially exploded when the number of inputs n τ becomes larger. Note that the base of the computation complexity for the forward approach equals n τ (the inputs' dimension), which varies according to different configurations of soft robots. Whereas, the dimension of the end-effector space is constant (n u ď 3; for the position case) because it is independent of any configuration of the investigated soft robot. Based on this observation, to reduce the computation complexity when estimating the workspace of soft robots, it is therefore, reasonable to propose novel methods that are based on discretizing the output space, i.e. the end-effector space u, which is invariant in dimension, instead of discretizing the input space, i.e., the actuators space T with variant dimension.

As stated in Section 2.5.4, the forward approach will be used as a reference to validate the workspace methodologies presented in the following chapters.

Chapter 3

Optimization-Based Approach for Workspace Estimation

Introduction

Due to the high computation complexity of the forward approach when estimating the workspace of soft robots, it is logical to propose an approach that can discretize the end-effector space since its dimension is smaller and constant compared to the actuators space for hyper-redundant robots, such as soft robots. In addition, it is more efficient to only map the exterior boundary of the workspace and avoid the heavy computation of its interior points.

Such a method is called the optimization-based approach [START_REF] Ja Snyman | An optimization approach to the determination of the boundaries of manipulator workspaces[END_REF] in this thesis and it consists of discretizing the end-effector space in order to map the exterior boundary of W E , noted as BW E , and deducing that the area/volume enclosed by BW E is the reachable workspace W E .

Generally speaking, this method consists of numerical algorithms used to map the surface enveloping soft robots' workspaces [START_REF] Amehri | Fem-based exterior workspace boundary estimation for soft robots via optimization[END_REF][START_REF] Walid | Discrete cosserat method for soft manipulators workspace estimation: An optimization-based approach[END_REF], which will be realized by the following steps:

1. In the first step, we establish the definition of a soft robot's workspace boundary, denoted as BW E .

2. Next, consistent with the workspace boundary definition, we seek to find a boundary point u b on BW E . The strategy proposed consists of selecting a radiating point v, then from this point, a ray is emanated along a certain direction until BW E is met (as illustrated by 3. Finally, a map of BW E can be achieved by solving the procedure that consists of finding u b for multiple successive rays with respective directions emanating at angular intervals.

Unlike the forward method (Section 2.5), the computation complexity of the proposed optimization-based approach depends on two factors: the dimension of the space to be discretized, and the discretization precision. By managing these two parameters, the optimization-based approach can reduce the complexity of the workspace estimation.

1. It can reduce the dimension of the space to be discretized in the forward approach. In the forward approach, we discretize the actuators space whose dimension n τ is varying according to the soft robot structure, but the space to be discretized in the optimizationbased approach is always less than or equal to 3 for the position case;

2. When discretizing the space with a prescribed precision (with n d sub-spaces), the forward method cannot reduce the value of n d , while the optimization-based approach, which consists only in finding then mapping the boundary points u b , results in reducing the total number of iterations.

Briefly, the proposed optimization-based approach can be summarized into the following procedure:

[Step 1] Defining the workspace boundary BW E . [Step 2] Finding a point on BW E . [Step 3] Mapping the workspace boundary BW E .
In the following, we present the detailed procedure of the optimization-based approach for both the PCS [START_REF] Walid | Discrete cosserat method for soft manipulators workspace estimation: An optimization-based approach[END_REF] and the FEM [START_REF] Amehri | Fem-based exterior workspace boundary estimation for soft robots via optimization[END_REF] modeling methods.

Optimization Approach for the PCS case

Workspace Boundary BW E Definition

Given a configuration of a soft robot, with its workspace W E defined in Definition 1, and τ P T " rτ , τ s, @u P W E , we can then denote the neighborhood of u with a specified positive radius ε by:

B ε puq " t u ε P R nu | u ´uε ă ε u (3.1)
based on which the workspace boundary BW E according to [START_REF] Ja Snyman | An optimization approach to the determination of the boundaries of manipulator workspaces[END_REF] can be defined below.

Definition 2. The workspace boundary BW E is a subset of its workspace W E , and defined by the following:

BW E " tu P W E | @ε ą 0, Du ε P B ε puq, s.t. u ε R W E u Remark 5.
To address the mapping of BW E , the strategy proposed is based on solving a constrained minimization problem defined in Section 3.2.2 of the optimization-based approach procedure. For this, we first need to calculate the gradient of the end-effector position upqq with respect to the actuators vector τ , denoted as ∇ τ " upqq ‰ , and established in the following theorem.

Theorem 1. For a given controllable configuration of the soft robot, with ∇ τ pqq defined in Assumption 1, the geometric Jacobian Jpqq defined in (2.7a), and the orientation matrix Rpqq computed via the definition of the configuration matrix gpqq through the following:

Rpqq

" E 1 gpqq E T 1 (3.2)
then the gradient of the end-effector position upqq with respect to the actuators vector τ , denoted as ∇ τ " upqq ‰ , is established by:

∇ τ " upqq ‰ " ∇ τ pqqJ T pqqr0 3 Rpqqs T (3.3)
with ∇ τ pqq " H T pqq " B ´1pτ, qq ‰ T being defined in Assumption 1, and Bpτ, q, L, αq " KpL, αq Hq pq, L, αqpI n b τ q ´Fq pq, L, αq.

Proof. Using the principle of variable separation, the gradient of the end-effector position ∇ τ " upqq ‰ can be written as follows:

∇ τ " upqq ı " " Bupqq Bτ ı T " " u q pqq " ∇ τ pqq ‰ T ı T
Based on the velocity twist vector ηpqq defined in (2.7c), we can derive the linear velocity vpqq via the following: vpqq " 9

upqq " E 1 ηpqq E T 2 , with ηpqq " g ´1pqq 9 gpqq

Next, we develop the term ηpqqE T 2 :

ηpqq E T 2 " g ´1pqq 9 gpqq E T 2 " g ´1pqq Bgpqq Bq pI 4 b 9 qq E T 2
Using the definition of the configuration matrix gpqq, we develop the above equation, and we obtain:

ηpqq E T 2 " " R ´1pqq R ´1pqqupqq 0 1  » - BRpqq Bq Bupqq Bq 0 0 fi fl " 0 9 q  " " R ´1pqq Bupqq Bq 9 q 0 ı T
Substituting the above equation of ηpqqE T 2 in the linear velocity 9 upqq, we achieve the following expression of 9 upqq:

9 upqq " E 1 ηpqq E T 2 " R ´1pqq Bupqq Bq 9 q (3.4)
Using the relation between the velocity vector ηpqq and the geometric Jacobian Jpqq defined by (2.7a), we derive another form of the linear velocity 9 upqq, which can be formulated as follows:

ηpqq " " wpqq 9 upqq ‰ T " Jpqq 9 q which gives: 9 upqq " r0 3 I 3 s Jpqq 9 q (3.5)

Using the expressions (3.4) and (3.5), we deduce the following equation: R ´1pqq Bupqq Bq 9 q " r0 3 I 3 s Jpqq 9 q

Then, we get: u q pqq " Bupqq Bq " Rpqq r0 3 I 3 s Jpqq " r0 3 Rpqqs Jpqq Finally, by substituting u q pqq in ∇ τ " upqq ‰ , we deduce the following: Assume that a radiating point v is selected, and that v is exterior to the workspace W E . It is now proposed that, consistent with Definition 2, a boundary point u b in a particular direction emanating from the radiating point v (as depicted by Fig. 3.2), can be determined by solving the following constrained optimization problem [START_REF] Ja Snyman | An optimization approach to the determination of the boundaries of manipulator workspaces[END_REF]:

∇ τ " upqq ‰ " ∇ τ pqq " u q pqq ‰ T " ∇ τ pqqJ T pqqr0 3 Rpqqs T
τ b " arg min τ upqq ´v 2 2 " f pqq s.t. τ P T K q " Hpqqτ `F pqq upqq " E 1 gpqq E T 2 (3.6)
where τ b is the optimal and feasible actuators vector (since it satisfies the bounded constraint) for the purpose of minimizing the nonlinear cost function f pqq " upqq ´v 2 2 , which represents the distance between v and upqq.

Based on the solution of τ b , and in accordance with Proposition 1, we calculate the corresponding generalized strain vector q b through (2.51). Then, we compute the strain twist ξb i related to the strain vector q b i of each segment i using (2.49b). Next, we derive the configuration matrix gpq b q associated to the generalized strain vector q b from (2.49a). Finally, using (2.48b), we deduce the boundary point position u b .

Regarding the implementation of the methodology, we need to compute the gradient of the objective function f pqq with respect to τ , denoted by ∇ τ `f pqq ˘, which can be formulated as follows:

∇ τ `f pqq ˘" 2∇ τ " upqq ‰ `upqq ´v˘( 3.7)
with ∇ τ " upqq ‰ defined in (3.3) . Another additional question linked to the implementation of the proposed methodology is how a radiating point v may be generated. As depicted by Fig. 3.2, the radiating point v must be exterior to the workspace. In practice, since the length L of a soft robot is bounded (L ď L), we only need to choose v such that || Ov|| ą L, where O " r0, 0, 0s T is the origin of the inertial frame.

Mapping the Workspace Boundary BW E

In this section, we present a numerical methodology used for mapping the workspace boundary BW E . Solving this problem is equivalent to finding a boundary point u b for multiple successive directions (as indicated in Fig. 3.3). In the following, we represent a simplified illustration of mapping the workspace boundary (see Fig. 3.3). The workspace boundary BW E may be numerically mapped by solving the optimization problem (3.8) for n d successive rays, with respective directions v pjq , j " 0, 1, 2, . . . , n d emanating at angular intervals of the angles δ z , δ x and δ y (with δ z,x,y " 2π{n d ) from a radiating point v 0 [START_REF] Ja Snyman | An optimization approach to the determination of the boundaries of manipulator workspaces[END_REF]:

τ bj " min τ }upqq ´vpjq } 2 2 s.t. τ P T K q " Bpqqτ `F pqq upqq " E 1 gpqq E T 2 (3.8)
where τ bj is the optimal and feasible actuators force vector in order to minimize the distance between v pjq and upqq.

From the solution τ bj , and according to Proposition 1, we calculate the corresponding generalized strain vector q bj (2.51). Then, we compute the strain twist ξbj i related to the strain vector q bj i of each segment i using (2.49b), and based on which we derive the configuration matrix gpq bj q associated to the generalized strain vector q bj from (2.49a). Finally, we deduce the boundary point position u bj using (2.48b).

The direction vector v pjq is expressed by the following:

v pjq " E pjq δ v 0
where the initial radiating point v 0 is exterior to the workspace and can be calculated using the same technique suggested in Section 3.2.2. The rotation matrix E pjq δ is defined as:

E pjq δ " E δ pjq z E δ pjq x E δ pjq y
, and E δ pjq z,x,y represents the basic rotation matrices about the z, x and y axis respectively [START_REF] Brian | Lie groups, Lie algebras, and representations: an elementary introduction[END_REF].

Up to this point, the procedure of the optimization-based approach has been completely described in each step with its related specifics.

Convexity Issues during the Mapping of BW E

Going a step further, difficulties may now arise in relation to some details of the methodology implementation, and it concerns the strategy to be adopted if non-convexity of the workspace W E interferes with the mapping of the workspace boundary BW E as illustrated in Fig. 3.4, where as a result of the non-convexity problem, the workspace may be over-estimated. Eventually, in this particular situation, mapping the boundary BW E may result in an overestimation of the workspace since all cyan colored boundary points in Fig. 3.4 between the successfully determined boundary points u bn 1 and u bn 2 , may not be identified, due to fact that, for all admissible direction vectors v pjq between v pn1q and v pn2q , the only possible solutions that can be achieved by the minimization problem (3.8) will be either u bn 1 or u bn 2 , i.e., @v pjq P rv pn1q v pn2q s,

(3.8) ùñ u bn 1 _ u bn 2
To solve this specific problem, we need first to understand the phenomenon leading to the occurrence of this situation. Actually, this particular behavior is mainly due to a sudden leap between two successive iterations, either gradual or declining, of at least one actuators force τ bj k P τ bj of the optimization problem (3.8).

In this thesis, and consistent with the above observation, the strategy proposed to solve this particular problem is firstly based on discretizing the space between the boundary points u bn 1 and u bn 2 , based on which we obtain a list of radiating points v l . In the next step, for each radiating point v l P ru bn 1 u bn 2 s we proceed by solving the minimization problem (3.6), which allows us to obtain the feasible actuators force vector τ b corresponding to each radiating point v l . Next, by solving (2.48), we can identify the boundary points u b on the non-convex part of the workspace (See the cyan colored boundary points in Fig. 3.4). Finally, by mapping all boundary points, we obtain the correct estimation of the workspace (as shown in Fig. 3.5).

Briefly, the proposed strategy to solve the non-convexity problem can be summarized into the following procedure: [Step. A] Calculate v l by discretizing the space between the boundary points u bn 1 and u bn 2 .

[Step. B] For each discretized point v l , find its related feasible actuators vector by solving (3.6), then deduce its corresponding boundary point from (2.48).

In the following illustration (see Fig. In the following, a brief algorithm is presented to describe the main steps of the optimizationbased approach along with the adopted strategy to solve the non-convexity problem (See Algorithm 1).

Validation and Analysis

In this section, we will illustrate the presented results by evaluating the workspace for different configurations of tendon-actuated soft robots. We consider various scenarios, and we implement both the forward and optimization-based approaches to get W E and BW E respectively.

Algorithm 1 Calculate BW E Require: v 0 , T , Fei , F p0q , E 1 , E 2 , K, E i , G i , A i , L i , n d , n s , r δ Ð 2π{n d ; BW E Ð ∅; Ź initialization Y Ð rT , E 1 , E 2 , K, E i , G i , A i , L i s; F Ð F p0q for j Ð 1 to n d do v pjq Ð E pjq δ v 0 Ź direction vector τ pjq b Ð solve p3.8q pv pjq , Y q Ź optimal solution ru pjq b , q pjq b s Ð solve (2.48) pτ pjq b , Y, F q Ź boundary point F Ð solve (2.13f) pq pjq b , Fei q BW E Ð BW E ' u pjq b Ź append if }u pj´1q b ´upjq b } 2 ą r then Ź non-convexity v l Ð dispu pjq b , u pj´1q b , n s q Ź discretize for s Ð 1 to n s do τ psq b Ð solve (3.6) pv psq l , Y q Ź optimal solution ru psq b , q psq b s Ð solve (2.48) pτ psq b , Y, F q F Ð solve (2.13f) pq psq b , Fei q BW E Ð BW E ' u psq b Ź append end for end if end for
In the following table (see Table 3.1), we give the numerical values of the Young's modulus E i , shear modulus G i , the Poisson's ratio ν, the inertial length L i , the radius r i , the second moments of area J px,y,zqi and the cross-section area A i for each segment i P t1, 2, 3u. For the optimization-based approach we propose to discretize the angles with a discretization step size of 0.05 Radian (unit).

Soft robot segments

In the following scenarios and also for the scenarios of the following workspace estimation chapters (Chapter 4 and Chapter 5), the obtained results of the workspace estimation via the proposed approaches are superimposed with the forward method results, and the 3D results are presented via 2D-views (Oxz, Oyx, and Oyz). Also, since in this thesis we only focus on the position aspect of the workspace (see Definition 1), the obtained simulation results represent only the feasible positions of the workspace but not the orientation aspect. Besides, please note that one operation is a complete iteration that contains all the necessary steps to achieve one single feasible point in the workspace. Also, the presented simulation scenarios were implemented using MATLAB (In the forward approach, we used the function "fsolve" to solve the forward statics, while the function "fmincon" was used to solve the optimization problems formulated in this thesis) on an Intel Xeon(R) with a 16-GB RAM and a 3.50 GHz processor.

Finally, for the forward approach, we propose to discretize the actuators force vector with a discretization step size of 1 Newton (unit). In addition, since we know from simulation the cost of a single operation of the forward approach in order to find a feasible configuration in the workspace (0.00165 second for a single operation), then, we only gave an approximation to the computational time for the scenarios where the investigated soft robot is actuated by more than 4 actuators since the workspace estimation of these scenarios via the forward approach requires a long time.

Scenario 1: planar 2-segments 2-tendons soft robot

In the first scenario, we consider a soft robot composed of two segments and actuated by two tendons (See Fig. 3.6a). The two tendons are installed on the Ozx-plan, and parallel to each other. The first tendon is fixed at the position p0, 0, ´r1 {2q and extends along the first segment length L 1 . The second tendon is fixed at the position p0, 0, ´r2 {2q and extends to pL 1 `L2 , 0, ´r2 {2q. Such a chosen structure allows us to obtain a 2D workspace. We choose a tension magnitude within T 1,2 " r0, 200sN that will allow us to apply a couple magnitude within C T 1 " r0, 6sN ¨m and C T 2 " r0, 4sN ¨m.

First, we use the forward approach to estimate the workspace W E of this soft robot. We choose n d " 200, and we obtain the following estimation of W E (See Fig. 3.6b). The computation cost of the forward approach for the two segments -two tendons soft robot is equal to 200 2 operations, and the workspace estimation took 66.35 seconds to obtain the full estimation of its workspace W E .

Next, we apply the optimization-based approach with n d " 120, and based on (3.8) we calculate all boundary points on the surface enveloping the workspace (See Fig. 3.6b). Finally, as illustrated by Fig. 3.6b, we achieve BW E of this particular soft robot configuration. The computation cost of the optimization-based approach for the studied soft robot is equal to 168 operations, and it required nearly 1.82 seconds to fully estimate its BW E .

Scenario 2: planar 2-segments 4-tendons soft robot

In the second scenario, we consider the soft robot of the first scenario and we add two symmetric tendons (See Fig. 3.7a). We choose a tension magnitude within T 1,¨¨¨,4 " r0, 100sN that will allow us to apply a couple magnitude within C T 1,3 " r0, 3sN ¨m and C T 2,4 " r0, 2sN ¨m.

First, we use the forward approach to estimate the workspace W E of this soft robot. We choose n d " 100, and we obtain the following estimation of W E (See Fig. 3.7b). The computation cost of the forward approach for the two segments -four tendons soft robot is equal to 100 4 operations, and the workspace estimation will take 1.65ˆ10 5 seconds to obtain the full estimation of its workspace W E .

Next, we apply the optimization-based approach with n d " 120, and based on (3.8) we calculate all boundary points on the surface enveloping the workspace (See Fig. 3.7b). Finally, as illustrated by Fig. 3.7b, we achieve BW E of this particular soft robot configuration. The computation cost of the optimization-based approach for the studied soft robot is equal to 168 operations, and it required nearly 3.92 seconds to fully estimate its BW E .

Scenario 3: planar 3-segments 3-tendons soft robot

In the third scenario, we consider adding another segment and a tendon to the first scenario (See Fig. 3.8a). The additional tendon is mounted in a parallel manner to the tendons in the first scenario: it is fixed at the position p0, 0, ´r3 {2q and extends along the soft robot length L 1 `L2 `L3 . Also, this configuration enables us to achieve a larger 2D workspace estimation compared to the first scenario. We choose a tension magnitude within T 1,2,3 " r0, 100sN that will allow us to apply a couple magnitude within C T 1 " r0, 3sN ¨m, C T 2 " r0, 2sN ¨m, and C T 1 " r0, 1sN ¨m.

First, we use the forward approach to estimate the workspace W E of this soft robot. We choose n d " 100, and we obtain the following estimation of W E (See Fig. 3.8b). The computation cost of the forward approach for the three segments -three tendons soft robot is equal to 100 3 operations, and the workspace estimation took 2962 seconds to obtain the full estimation of its workspace W E .

Next, we apply the optimization-based approach with n d " 120, and based on (3.8) we calculate all boundary points on the surface enveloping the workspace (See Fig. 3.8b). Finally, as illustrated by Fig. 3.8b, we achieve BW E of this particular soft robot configuration. The computation cost of the optimization-based approach for the studied soft robot is equal to 220 operations, and it required nearly 3.52 seconds to fully estimate its BW E .

Scenario 4: planar 3-segments 6-tendons soft robot

In the fourth scenario, three additional tendons are routed in a fashion where they are symmetric to the tendon configuration in the third scenario (See Fig. 3.9a). We choose a tension magnitude within T 1,¨¨¨,6 " r0, 100sN that will allow us to apply a couple magnitude within C T 1,4

" r0, 3sN ¨m, C T 2,5

" r0, 2sN ¨m, and C T 3,6 " r0, 1sN ¨m.

First, we use the forward approach to estimate the workspace W E of this soft robot. We choose n d " 100, and we obtain the following estimation of W E (See Fig. 3.9b). The computation cost of the forward approach for the three segments -six tendons soft robot is equal to 100 6 operations, and the workspace estimation will take 1.65ˆ10 9 seconds to obtain the full estimation of its workspace W E .

Next, we apply the optimization-based approach with n d " 120 , and based on (3.8) we calculate all boundary points on the surface enveloping the workspace (See Fig. 3.9b). Finally, as illustrated by Fig. 3.9b, we achieve BW E of this particular soft robot configuration. The computation cost of the optimization-based approach for the studied soft robot is equal to 220 operations, and it required nearly 13.58 seconds to fully estimate its BW E .

Scenario 5: spatial 3-segments 9-tendons soft robot

In the final scenario, three additional tendons are added to the configurations of the above scenario in order to obtain a 3D workspace. The three additional tendons are mounted in the following way; the seventh, eighth and ninth tendons are routed in a manner where the seventh tendon is fixed at the position p0, ´r1 {2, 0q and extends along the first segment length L 1 , the eighth tendon is fixed at the position p0, ´r2 {2, 0q and extends to pL 1 `L2 , ´r2 {2, 0q, and the ninth tendon is fixed at the position p0, ´r3 {2, 0q and extends along the soft robot length L 1 `L2 `L3 (See Fig. 3.10a). With the tendons being mounted in this fashion, the actual configuration allows us to have a 3D workspace. We choose a tension magnitude within T 1,...,9 " r0, 50sN that will allow us to apply a couple magnitude within C T 1,4,7 " r0, 1.5sN ¨m, C T 2,5,8

" r0, 1sN ¨m, and C T 3,6,9 " r0, 0.5sN ¨m.

First, we use the forward approach to estimate the workspace W E of this soft robot. We choose n d " 100, and we obtain the following estimation of W E (See Figs. 3.10b, 3.10c, 3.10d). The computation cost of the forward approach for the three segments -nine tendons soft robot is equal to 100 9 operations, and the workspace estimation will take 1.65 ˆ10 12 seconds to obtain the full estimation of its workspace W E .

Next, we apply the optimization-based approach with n dy " 120 and n dz " 72, and based on (3.8) we calculate all boundary points on the surface enveloping the workspace (See Figs.

3.10b, 3.10c, 3.10d), and we finally achieve BW E of this particular soft robot configuration. The computation cost of the optimization-based approach for the studied soft robot is equal to 8640 operations, and it required nearly 510 seconds to fully estimate its BW E .

Notes on scenarios results

Based on the achieved results, we can clearly observe that the workspace boundary estimation result achieved via the optimization-based method encloses the estimated workspace obtained via the forward approach.

In addition, the proposed optimization-based approach can yield equivalent workspace estimation precision compared to the forward method, and simultaneously reduce both the number of operations and the computational time required for the workspace estimation for each scenario.

Notes on computational efficiency

For each particular configuration of a tendon-actuated soft robot, we showed the estimation of its workspace using both the forward and optimization-based approach. From the simulation results, we observe that the proposed optimization-based approach consisting of mapping the exterior boundaries of soft robots' workspaces has significantly reduced the computation cost of the workspace estimation for each soft robot configuration.

In the following table (see Table . 3.2), we summarize the operations complexity and the computational time for each scenario. Visually, the results of the above table are depicted in Figs. 3.11a and 3.11b, and we can clearly observe that, when the size of the actuators, (i.e., the number of tendons) increases, the computation complexity of the forward approach explodes exponentially, whereas the computation complexity of the optimization-based approach remains almost linearly stable.

Optimization Approach for the FEM case 3.3.1 Implementation of the Methodology

The implementation of the optimization-based approach on the FEM model [START_REF] Amehri | Fem-based exterior workspace boundary estimation for soft robots via optimization[END_REF] follows the same steps described for the PCS model and detailed in Section 3.2.

In the first step, the definition of the workspace boundary for the FEM model is similar to that of the PCS method (see Definition 2 in Section 3.2.1).

In the second step, the goal is the same, which is to find a point on BW E . However, the optimal problem to find a boundary point is formulated under different constraints for the FEM model, and established as follows:

rτ b , q b s " arg min τ,q }upqq ´v} 2 2 s.t. τ P T Kpqqq " Hpqqτ `F pqq u " Cq `up0q (3.9)
where the obtained actuators vector τ b with its associated displacement vector q b represents the optimal and feasible solution of (3.9). Finally, we deduce the related end-effector boundary position using the following:

u b " Cq b `up0q
Similar to the PCS case, the third step consists of mapping the workspace boundary BW E . However, the optimal problem needs to be adapted for the FEM model, and formulated as follows:

rτ bj , q bj s " arg min

τ,q }upqq ´vpjq } 2 2 s.t. τ P T Kpqqq " Hpqqτ `F pqq u " Cq `up0q (3.10)
where the achieved actuators vector u bj with its associated displacement vector q bj is the optimal and feasible solution to minimize the distance between v pjq and u . Accordingly, the associated end-effector boundary position can be deduced as:

u bj " Cq bj `up0q
Finally, the solution to non-convexity problems that may arise using the FEM model can be solved following the same procedure described in Section 3.2.4 for the PCS case. However, we need to use the corresponding optimization equations (3.9) and the related FEM static model established in (2.53) during Step. B of the proposed strategy.

The algorithm of the optimization-based approach for the FEM case is presented in the following (See Algorithm 2):

Validation and Analysis

In this section, we illustrate the presented results by evaluating the reachable workspace of a trunk-like soft robot with a Young's modulus E " 1.8 ˆ10 7 P a, and a Poisson's ratio ν " 0.45.

We consider three different scenarios where the trunk-like soft robot is actuated by different cable routing configurations, then we implement both the forward and optimization-based approach in order to estimate the reachable workspace W E and the workspace boundary BW E , respectively, for each scenario.

For the optimization-based approach we propose to discretize the angles with a discretization step size of 0.08 Radian (unit).

Require: v 0 , C, T , n d , n s , r, u p0q δ Ð 2π{n d ; BW E Ð ∅; Ź Initialization for j Ð 1 to n d do v pjq Ð R
pjq δ v 0 Ź direction vector rτ bj , q bj s Ð solve (3.10) pv pjq , C, T , u p0q q Ź optimal solution u bj Ð Cq bj `up0q Ź boundary point

BW E Ð BW E ' u bj Ź append if }u bj´1 ´ubj } 2 ą r then Ź non-convexity v l Ð dispu bj´1 , u bj , n s q Ź discretize for s Ð 1 to n s do
rτ bs , q bs s Ð solve (3.9) pv psq l , C, T , u p0q q u bs Ð Cq bs `up0q Ź boundary point BW E Ð BW E ' u bs Ź append end for end if end for 3.3.2.1 Scenario 1: trunk-like soft robot actuated by two tendons.

In the first scenario, we consider a cable routing configuration where the trunk-like soft robot is actuated by two tendons (as shown by Fig. 3.12a). We apply a force whose magnitude is within a specific range T p0q " T " r0 100s ˆr0 100s.

Following the procedure presented in Section 2.5.3, we firstly apply the forward approach with a discretization precision n d " 100, and we obtain the corresponding estimation of W E (illustrated by blue colored points in Figs. 3.12b, 3.12c and 3.12d) for this soft robot's configuration. In term of operations' computation complexity, the forward approach requires n nτ d " 100 2 operations, and in term of time's computation, the workspace estimation took 1463 seconds to obtain the full estimation of this scenario's workspace.

Next, we use (3.10) (See Fig. 3.3) of the optimization-based approach with n dz " 72 and n dx " 12 to map the boundary points of this particular configuration.

For this scenario, the optimization-based approach required a total of 864 operations with a time computation of 984 seconds. The obtained boundary points are depicted by green colored points in Figs. 3.12b, 3.12c and 3.12d, which represent the workspace boundary BW E of this soft robot's configuration. In this scenario, we consider a cable routing configuration where the trunk-like soft robot is actuated by three tendons (as shown by Fig. 3.13a). We apply a force whose magnitude is within a specific range T p0q " T " r0 100s ˆr0 100s ˆr0 100s.

First, we apply the forward approach with a discretization number equal n d " 100, and to achieve the workspace W E estimation (illustrated by blue colored points in Figs. 3.13b, 3.13c and 3.13d) of this particular configuration. In term of operations' computation complexity, the forward approach requires n nτ d " 100 3 operations, and in term of time's computation, the workspace estimation will take 1.4 ˆ10 5 seconds to obtain the full estimation of this scenario's workspace.

Next, using (3.10) of the optimization-based approach with n dz " 72 and n dx " 12, we map the boundary points of this particular configuration.

For this scenario, the proposed optimization-based approach required a total of 864 operations with a time computation of 1034 seconds. The achieved boundary points are illustrated by green colored points in Figs. 3.13b, 3.13c and 3.13d, which represent the workspace boundary BW E of this particular configuration. In the final scenario, we consider a cable routing configuration where the trunk-like soft robot is actuated by four symmetric tendons (as shown by Fig. 3.14a). We apply a force whose magnitude is within a specific range T p0q " T " r0 100s 1 ˆ¨¨¨ˆr0 100s 4 .

First, we apply the forward approach with a discretization precision equal n d " 100, and we obtain the corresponding estimation of W E (illustrated by blue colored points in Figs. 3.14b, 3.14c and 3.14d) for this soft robot's configuration. In term of operations' computation complexity, the forward approach requires n nτ d " 100 4 operations, and in term of time's computation, the workspace estimation will take 1.4 ˆ10 7 seconds to achieve the full estimation of this scenario's workspace.

Next, we use (3.10) of the optimization-based approach with n dz " 72 and n dx " 12 to map the boundary points of this particular configuration.

For this scenario, the optimization-based approach required a total of 864 operations with a time computation of 1039 seconds. The obtained boundary points are depicted by green colored points in Figs. 3.14b, 3.14c and 3.14d, which represent the workspace boundary BW E of the studied soft robot's configuration.

Notes on computational efficiency

In the following table (see Table 3.3), and Figures. 3.15a and 3.15b, we summarize the operations complexity and time computation for each scenario.

Conclusion

In this chapter, we have proposed an optimization-based approach that consists of mapping the exterior boundaries of the workspace of soft robots.

The proposed approach was successfully implemented to both the PCS and the FEM models, where we have shown its efficiency in reducing the computation complexity and computational time necessary to estimate the workspace of soft robots in contrast to the forward approach which explodes exponentially when increasing the dimension of the actuators.

On the other hand, since this approach consists only in mapping the exterior boundary of the workspace by determining the boundary points lying on the surface enclosing the workspace, it surely will not be able to identify inside information of the workspace such as interior boundaries.

65

Chapter 4

Interval Analysis Approach for Workspace Estimation

Due to the limitations encountered by the optimization approach that cannot provide knowledge on the interior configurations of the workspace, it is therefore necessary to provide a method that can overcome this limitation, but is also based on the second strategy, i.e., of discretizing the end-effector's space.

Such a method is called the forward-backward approach [START_REF] Amehri | Fem based workspace estimation for soft robots: a forward-backward interval analysis approach[END_REF] that is based on interval analysis techniques [START_REF] Gouttefarde | Interval-analysis-based determination of the wrench-feasible workspace of parallel cable-driven robots[END_REF][START_REF] Jaulin | Interval analysis[END_REF][START_REF] Merlet | Determination of 6d workspaces of gough-type parallel manipulator and comparison between different geometries[END_REF][START_REF] Merlet | Parallel robots[END_REF][START_REF] Merlet | On the workspace of suspended cable-driven parallel robots[END_REF], and it consists of discretizing the end-effector's space by starting from an initial reachable configuration then exploring the whole possible reachable space in order to finally estimate the workspace of a soft robot.

In the following, we will first provide a brief introduction about interval analysis techniques used by this approach, and then an introduction of the approach itself followed by its realization for both the PCS and FEM models.

Interval Analysis Introduction

The proposed forward-backward approach uses interval analysis techniques in order to estimate the workspace of soft robots. Thus, a brief introduction to some interval analysis techniques used in this thesis is presented in the following. More details can be found in [START_REF] Jaulin | Interval analysis[END_REF].

Consider two real intervals X 1 " rX 1 , X 1 s and X 2 " rX 2 , X 2 s, then according to [START_REF] Jaulin | Interval analysis[END_REF], the product of X 1 and X 2 can be calculated by the following:

X 1 X 2 " " min `CpX 1 , X 2 q ˘, max `CpX 1 , X 2 q ˘‰ with CpX 1 , X 2 q " tX 1 X 2 , X 1 X 2 , X 1 X 2 , X 1 X 2 u.
The product of a real matrix A P R by a real interval X " rX , X s with appropriate dimension is calculated as follows:

AX " # " AX , AX ‰ , A ą 0 " AX , AX ‰ , A ă 0
However, due to the wrapping effect ( AX Ą tAx | x P X u ) [START_REF] Jaulin | Interval analysis[END_REF], this technique yields an over-estimation of the resulted space which can be explained by considering the following example: Therefore, it is necessary to propose a technique to eliminate the resulted over-estimation in order to achieve the correct and accurate space, and this will be implemented in the proposed forward-backward approach and explained in the following.

A " ˆ1 1 0 1 ˙, X " ˆr´1, 0s r1

Forward-Backward Approach Introduction

In this section, we present a brief overview of how the proposed forward-backward approach operates.

Firstly, we need to discretize the end-effector space. To realize this, we need to have an idea about its bounds. In this approach, we approximate the relation between q and τ at the pj ´1qth iteration as follows:

Kδq pjq " Hpq pj´1q qδτ pj´1q `F pq pj´1q q, PCS and Kpq pj´1q qδq pjq " Hpq pj´1q qδτ pj´1q `F pq pj´1q q, FEM with δq pjq " q ´qpj´1q and δτ pj´1q " τ ´τ pj´1q , where δq pjq and δτ pj´1q represent respectively the generalized strain and the actuators vectors with respect to the prior configuration (for the initial state j " 1, upq p0q q and q p0q are the end-effector's initial position and the initial generalized strain vector respectively, and τ p0q is the initial actuators force vector). Based on the above result, a forward estimation of the end-effector workspace related to δq pjq can be obtained (see Step 1 in Fig. 4.2), noted as W pjq E , by using interval analysis techniques [START_REF] Jaulin | Interval analysis[END_REF] (which will be detailed in Section 4.3.1 for the PCS case and Section 4.4 for the FEM case). Once it is estimated, we can discretize W pjq , we calculate a feasible small neighborhood, noted as W pjq Es (subscript "s" stands for small neighborhood), around the current end-effector position u pj´1q (see Step 2 in Figs. [START_REF] Alciatore | Determining manipulator workspace boundaries using the monte carlo method and least squares segmentation[END_REF].2 and 4.3). After that, we proceed backwardly to the actuators domain T in order to compute the corresponding actuators τ pjq (see Step 3 in Figs. 4.2 and 4.4), with which the investigated soft robot may reach W pjq Es . This step, i.e. Step 3, is necessary to obtain the accurate corresponding end-effector position u pjq which is approximately equal to the associated point in the set W pjq

Es . In addition, this step is also necessary for the computation of the new actuators bounds T pjq in Step 4. After, by utilizing the actuators vector τ pjq inferred from Step 3, we calculate the new actuators bound T pjq (shown as Step 4 in Figs. 4.2 and 4.5), allowing us to forwardly construct the new end-effector bound W pj`1q E which will be required for the next iteration of the proposed approach. However, since interval analysis techniques lead to an over-estimation of W pjq E (due to the wrapping effect [START_REF] Jaulin | Interval analysis[END_REF]), this thesis proposes a simple strategy to accurately reduce the dimension of W pjq E . The strategy proposed to reduce the overestimation lies within two steps of the forward-backward approach. In part, Step 2 indeed allows us to reduce the over-estimation by discarding certain unfeasible small neighborhoods because the feasible small neighborhoods W pjq

Es must belong to the discretized space W pjq 2). On the other hand, the over-estimation correction occurs mainly at the level of Step 5 of the approach where we dismiss pre-explored feasible small neighborhoods, and therefore efficiently and precisely reduce the over-estimated workspace (which is illustrated by Fig. 4.6, and will be explained in detail in

Step 5, Section 4.3.5).

Eventually, to achieve the whole estimation of W E using the forward-backward method, we iterate the above procedure starting from the initial state j " 1 until the stop condition (which will be detailed in Compared to the optimization-based approach, the forward-backward approach is exhaustive since it seeks to explore the whole workspace providing its interior and exterior points.

Forward-Backward Approach for the PCS case

In this section, we present the detailed realization of the forward-backward approach applied to the PCS model following the portrayed scheme in Fig. 4.2.

4.3.1

Step 1: Forward estimation of W pjq E In the first step, we use interval analysis techniques to forwardly estimate the workspace of the end-effector.

Precisely, suppose that at the pj ´1qth iteration we have the information of the end-effector position u pj´1q , the generalized strain vector q pj´1q and the corresponding actuators vector τ pj´1q P T pj´1q (for the case j " 1, it corresponds to the initial configuration of the studied soft robot). Then, it is now proposed that (2.48) can be approximated by the following:

Kδq pjq " Hpq pj´1q qδτ pj´1q `F pq pj´1q q (4.1a)

u pjq " E 1 gpδq pjq q E T 2 (4.1b)
with δq pjq " q ´qpj´1q and δτ pj´1q " τ ´τ pj´1q .

Let us note the matrix A as:

Apq pj´1q q " K ´1H pq pj´1q q Define now A `" maxt0, Au
where the maxt¨u operation is applied element-wise [START_REF] Zheng | Design of interval observer for a class of uncertain unobservable nonlinear systems[END_REF], and

A
´" A ´At hen using the following theorem, we achieve an over-estimation of the end-effector workspace W pjq E corresponding to the given actuators bound T pj´1q .

Theorem 2. For a given controllable configuration of the soft robot, with the actuators bound T pj´1q defined as T pj´1q " rδτ pj´1q , δτ pj´1q s where δτ pj´1q and δτ pj´1q represent the lower and upper bounds of the actuators vector with respect to the prior configuration respectively. An over-estimated end-effector workspace W pjq E is given by

W pjq E " ! u pjq P R nu |u pjq P " u pjq , u pjq ı) (4.2)
with u pjq " E 1 gpδq pjq qE T 2 , and u pjq " E 1 gpδq pjq qE T 2 the configuration matrix bounds gpδq pjq q P rgpδq pjq q, gpδq pjq qs , and the generalized strain vector bounds δq pjq P rδq pjq , δq `A´δ τ pj´1q `K´1 F pq pj´1q q Proof. From (4.1a) we can get the generalized strain vector as follows:

δq pjq " K ´1H pq pj´1q qδτ pj´1q `K´1 F pq pj´1q q " Apq pj´1q qδτ pj´1q `K´1 F pq pj´1q q Using Lemma 3 in [START_REF] Zheng | Design of interval observer for a class of uncertain unobservable nonlinear systems[END_REF], we get the following:

A `δτ pj´1q `A´δ τ pj´1q ď Aδτ pj´1q ď A `δτ pj´1q `A´δ τ pj´1q
and deduce the generalized strain vector bounds as: `A´δ τ pj´1q `K´1 F pq pj´1q q

δq pjq ď
Next, applying the same lemma to (2.49b), and since B qi ě 0, we obtain the strain twist bounds as:

ξpjq i ď ξpjq i ď ξpjq i
After, by implementing the above result in (2.49a), and according to [START_REF] Jaulin | Interval analysis[END_REF], we get the configuration matrix bounds: gpδq pjq q ď gpδq pjq q ď gpδq pjq q Finally, we deduce the bounds of the end-effector position: The proposed strategy is to discretize W pjq E using a globally uniform spatial grid with a global small pre-defined precision s E d in each direction x, y and z, i.e.,

u pjq ď
W pjq E d " ! u pjq P W pjq E |u pjq " u pj´1q `sE d Γ, for Γ P Z 3 ) (4.3)
Next, an η-neighborhood of the configuration u pj´1q can be defined as:

S η pu pj´1q q " ! ũ P R 3 |}u pjq ´upj´1q } 8 " η ) with η " s E d . Consequently, with the discretized space W pjq E d
being established, we can then define the feasible small neighborhood of the end-effector position u pj´1q as follows:

W pjq Es " W pjq E d X S η pu pj´1q q (4.4)
The relationship between the estimated workspace W pjq E , its discretization space W pjq

Step 3: Backward Estimation of the Corresponding Configuration

The third step's purpose is to choose the discretized points in the feasible small neighborhood space W pjq Es , except the point u pj´1q , and then backwardly calculate the feasible actuators vector τ pjq which can drive the soft robot's end-effector to reach those points. The corresponding scheme is depicted by Fig.

Precisely, since any point ûpjq P W pjq

Es zu pj´1q is located inside the feasible small region of u pj´1q , then we can use this value to compute the associated configuration by solving the following nonlinear constrained optimization problem:

τ pjq " arg min τ upqq ´û pjq 2 2 s.t. τ P T K q " Hpqqτ `F pqq upqq " E 1 gpδq pjq q E T 2 (4.5)
where τ pjq is the optimal and feasible solution (since it satisfies the bounded and nonlinear constraints) for the purpose of minimizing the nonlinear cost function upqq ´û pjq 2 2 , which represents the distance between ûpjq and upqq.

With the calculated τ pjq , and in accordance with Proposition 1, we can calculate the corresponding generalized strain vector q pjq through (2.51). Then, we compute its related strain twist ξpjq i using (2.49b). Next, we derive the associated configuration matrix gpq pjq q via (2.49a). Finally, using (2.48b), we may deduce the corresponding position u pjq .

Regarding the implementation of the methodology, the gradient of f pqq with respect to τ , denoted by ∇ τ `f pqq ˘can be analytically formulated as:

∇ τ `f pqq ˘" 2∇ τ " upqq ‰ `upqq ´û pjq with ∇ τ " upqq ‰ " ∇ τ pqqJ T pqqr0 3
Rpqqs T being already established in Theorem 1, and ∇ τ pqq " H T pqq " B ´1pτ, qq ‰ T being defined in Assumption 1, with Bpτ, q, L, αq " KpL, αq´H q pq, L, αqpI n b τ q ´Fq pq, L, αq.

Step 4: Calculation of New Actuators Bound

In this step, we use the value of τ pjq to compute the associated new bound, noted as T pjq , from T pj´1q . For this, denote τ pjq k for 1 ď k ď n τ as the kth value of τ pjq , and diagpτ pjq q as the related n τ ˆnτ diagonal matrix whose entries are the n τ elements of the vector τ pjq , i.e.,

diagpτ pjq q " diag ! τ pjq 1 , ¨¨¨, τ pjq nτ )
Then, the new bound T pjq , illustrated by Fig. 4.5, can be determined using the value of τ pjq from Step 3 and the value of T pj´1q as follows:

T pjq " T pj´1q ´diagpτ pjq q ˆ1nτ ˆ2 (4.6)

where 1 nτ ˆ2 represents the n τ ˆ2 matrix with all entries equal to 1.

Step 5: Iteration and Stop Condition

For a particular configuration u pj´1q with its corresponding over-estimated workspace W pjq ). An important property of Step 2 allows us to not take into account certain non-feasible neighborhoods around u pj´1q , which results in diminishing the number of feasible small neighborhoods W pjq Es to be explored for the next iterations of the forward-backward approach, as illustrated in Fig. 4.3, by verifying the condition (4.4).

At the same time, given the information of W pj´1q

Es , and the newly calculated feasible small neighborhood W pjq Es from (4.4), an additional characteristics allowing us to further shrink the space to be explored is that some of the points belonging to W pjq

Es have already been explored by the former iterations, noted as W pj´1q EY , and defined by recurrence as follows:

W pj´1q EY " Y j´1 k"1 W pkq Es
The above elimination procedure is executed in Step 5, as depicted by Fig. 4.6, by verifying the following condition: W pjq Es Ð S η pu pj´1q qz `Wpj´1q EY X S η pu pj´1q q ˘(4.7) Combining these two conditions (Eq. (4.4) in Step 2 and Eq. (4.7) in Step 5) enables us to shrink the space to be explored efficiently and accurately, resulting in a significant decrease of the computation complexity of the proposed forward-backward approach (as it will be shown in the simulation results in Section 4.3.6).

Consequently, the stop condition of the proposed forward-backward method can be formulated as follows: Dj P N, such that ũ P W pj´1q EY , @ũ P W pjq Es or equivalently Dj P N, such that W pjq EY " W pj´1q EY implying that when backwardly estimating the related actuators vector τ pj´1q of each point ũ P W pjq Es , and then, if the associated end-effector position u pjq belongs to the union of all feasible small neighborhood set of the last iterations W pj´1q EY (as shown in Fig. 4.7), we can then conclude that these configurations have already been explored. Finally, we deduce that all feasible small neighborhoods were explored and no more reduction of the over-estimated workspace is possible. At this level, the estimated workspace is:

W E " W pj´1q EY .
In the following, a brief algorithm is presented in the table below to describe the main steps of the forward-backward approach.

Algorithm 3 Calculate tu, q, τ, T , W Es , W EY u pjq , W E Require: Γ, s E d , tu, q, τ, T , W Es , W EY u pj´1q while W pjq Es R W pj´1q EY do Ź Stop Condition. W pjq E Ð f 1 ptu, q, T u pj´1q q Ź Step 1: Forward estimation of W pjq E , eq. (4.2) W pjq E d Ð f 2 pu pj´1q , Γ, s E d q Ź Step 2: Discretized space W pjq E d
, eq. ( 4.3)

W pjq Es Ð W pjq E d
X S η pu pj´1q q Ź Step 2: Feasible neighborhood W pjq Es , eq. (4.4) ûpjq Ð W pjq Es zu pj´1q Ź Exclusion. tu, q, τ u pjq Ð f 3 pû pjq , q pj´1q q Ź Step 3: Backward estimation of u pjq , eq. (4.5) T pjq Ð f 4 pτ pjq , T pj´1q q Ź Step 4: New actuators bound T pjq , eq. (4.6)

W pjq

Es Ð S η pu pj´1q qzpW pj´1q EY X S η pu pj´1q qq Ź Step 5: Dismiss explored area, eq. (4.7)

W pjq EY Ð W pj´1q EY Y W pjq Es Ź Append. tu, q, τ, T , W EY u pj´1q Ð tu, q, τ, T , W EY u pjq Ź Update Values. j Ð j `1 Ź Next iteration. end while W E Ð W pj´1q EY Ź The Estimated Reachable Workspace.

Validation and Analysis

In this section, the same soft robots' configurations presented in the scenarios of Section 3.2.5 in Chapter 3 are considered, where we implement both the forward and forward-backward approaches to deduce the workspace estimation W E . Table 3.1 presents the numerical values of the investigated soft robot.

In the following scenarios, the workspace obtained from the forward approach is depicted by black-colored points. We suggest discretizing the end-effector space with a discretization step size of s E d " 0.02 Meter (unit).

Scenario 1: planar 2-segments 2-tendons soft robot

In the first scenario, we consider a soft robot composed of two segments and actuated by two tendons (See Fig. 4.8a). The two tendons are installed as the first scenario of Section 3.2.5.1 in Chapter 3. We choose a tension magnitude within T 1,2 " r0, 200sN that will allow us to apply a couple magnitude within C T 1 " r0, 6sN ¨m and C T 2 " r0, 4sN ¨m. Next, we apply the proposed forward-backward interval analysis approach by following the procedure presented in Section 4.3. Using T p0q we can forwardly estimate W p1q For this scenario, the forward-backward approach took 36 complete iterations and found a total of 1820 feasible points in the workspace with a time computation of 17.5 seconds. Those obtained feasible points are illustrated by red-colored points in Fig. 4.8b, and the union of those feasible points represents the workspace W E estimation, which is depicted by the gray zone with green contour. Moreover, for a complete 36 iterations it is expected to have W E d " 6 36 « 10 28

Scenario 2: planar 2-segments 4-tendons soft robot

In the second scenario, we consider the soft robot of the first scenario and we add two symmetric tendons (See Fig. 4.9a). We choose a tension magnitude within T 1,2 " r0, 150sN that will allow us to apply a couple magnitude within C T 1,3 " r0, 4.5sN ¨m and C T 2,4 " r0, 3sN ¨m.

Next, we apply the proposed forward-backward interval analysis approach by following the procedure presented in Section 4.3. The forward-backward approach in this scenario took 47 complete iterations and computed a total of 3640 feasible points in 187.63 seconds. Those obtained feasible points are illustrated by red-colored points in Fig. 4.9b, and their union represents the workspace W E estimation, which is illustrated by the gray zone with green contour. In addition, we actually explored a total of W E d " 13464 points, and a final total of W Es " 3640 feasible points, instead of an expected W E d " 6 47 « 3 ˆ10 36 points to explore.

Scenario 3: planar 3-segments 3-tendons soft robot

In the third scenario (See Fig. 4.10a), we consider the same tendons routing of the third scenario in Section 3.2.5.3 from Chapter 3. We choose a tension magnitude within T 1,2,3 " r0, 100sN that will allow us to apply a couple magnitude within C T 1 " r0, 3sN ¨m, C T 2 " r0, 2sN ¨m, and C T 1 " r0, 1sN ¨m.

Next, following the procedure presented in Section 4.3, we apply the proposed interval analysis approach. For this scenario, the forward-backward approach took 32 complete iterations and computed a total of 2275 feasible points in a duration of 72.84 seconds. Those obtained feasible points are depicted by red-colored points in Fig. 4.10b, and their union is illustrated by the gray zone with green contour, which gives the workspace W E estimation. Besides, we actually explored a total of W E d " 8415 points, and a final total of W Es " 2275 feasible points, instead of an expected W E d " 6 32 « 7 ˆ10 24 points to explore.

Scenario 4: planar 3-segments 6-tendons soft robot

In this scenario, three additional tendons are routed in a fashion where they are symmetric to the tendons configuration in the third scenario (See Fig. 3.9a). We choose a tension magnitude within T 1,¨¨¨,6 " r0, 100sN that will allow us to apply a couple magnitude within C T 1,4 " r0, 3sN ¨m, C T 2,5

" r0, 2sN ¨m, and C T 3,6 " r0, 1sN ¨m.

Next, the proposed interval analysis approach is applied by following the procedure presented in Section 4.3. The forward-backward approach of this scenario took 55 complete iterations and computed a total of 4550 feasible points in the workspace within 2436 seconds. The red-colored points in Fig. 4.11b depict those obtained feasible points, where their union represents the workspace W E estimation, which is illustrated by the gray zone with green contour. Furthermore, for a complete 55 iterations it is expected to have W E d " 6 55 « 6 ˆ10 42 points to explore. However, we actually explored a total of W E d " 16830 points, and a final total of W Es " 4550 feasible points. In the final scenario (See Fig. 4.12a), we consider the same tendons routing of the last scenario in Section 3.2.5.5 from Chapter 3, but in this scenario, the three additional tendons are mounted upward (i.e., on the positive level of the Oy axis). We choose a tension magnitude within T 1,...,9 " r0, 50sN that will allow us to apply a couple magnitude within C T 1,4,7 " r0, 1.5sN ¨m, C T 2,5,8

" r0, 1sN ¨m, and C T 3,6,9 " r0, 0.5sN ¨m.

Next, we apply the proposed forward-backward interval analysis approach by following the procedure presented in Section 4.3. This approach took 35 complete iterations and computed a total of 13350 feasible points with a computational time of 4679 seconds. Those obtained feasible points are depicted by red-colored points in Figs.4.12b, 4.12c and 4.12d, and their union of represents the workspace W E estimation, which is depicted by the gray area with green contour. Furthermore, we actually explored is a total of W E d " 50490 points, and a final total of W Es " 13350 feasible points, instead of the expected W E d " 6 35 « 10 27 points to explore. .12: Scenario 5 -W E estimation via forward approach (black points) and forwardbackward approach (gray area).

Notes on scenarios results

The proposed forward-backward approach yields comparable workspace estimation precision compared to the forward method, and at the same time reduces the complexity required to estimate the workspace of soft robots.

In fact, this significant reduction introduced by the proposed forward-backward approach was partially due to the elimination process of inappropriate configurations from Step. 2 (See Fig. 4.3), but mainly due to the dismissal of all feasible configurations that have already been explored during Step. 5 (See Fig. 4.6).

Notes on computational efficiency

In Table 4.1, and Figures. 4. 4.4 Forward-Backward Approach for the FEM case

Implementation of the Methodology

The implementation of the forward-backward approach for the FEM model follows the same steps described for the PCS model (which was detailed in Section 4.3, and illustrated by Fig.

4.2). However, Step 1 and

Step 3 of the forward-backward needs to be adapted for the FEM case.

For the first step of the this approach, we firstly need to linearly approximate the workspace defined by (2.53) with the knowledge of the position vector of the end-effector u pj´1q , the displacement vector q pj´1q , and the corresponding actuators vector τ pj´1q P T pj´1q , as follows: Kpq pj´1q qδq pjq " Hpq pj´1q qδτ pj´1q `F pq pj´1q q u pjq " Cδq pjq `upj´1q (4.8) with δq pjq " q ´qpj´1q and δu pj´1q " u ´upj´1q . Next, suppose that the actuators are well installed in a manner that the investigated soft robot is always controllable [START_REF] Zheng | Controllability pre-verification of silicone soft robots based on finite-element method[END_REF], then for any given value of q pj´1q the matrix Kpq pj´1q q is always invertible. Let us introduce the matrix A as: Apq pj´1q q " CK ´1pq pj´1q qHpq pj´1q q and define the matrix A `as before by A `" maxt0, Au and A ´" A ´AẀ ith the above necessary elements, for the related displacement δq pjq , the over-estimated workspace W pjq E corresponding to the actuators bound T pj´1q can be calculated via the following theorem:

Theorem 3. Given a soft robot's controllable configuration, with the actuators bound T pj´1q defined as T pj´1q " rδτ pj´1q , δτ pj´1q s with δτ pj´1q and δτ pj´1q represents the lower and upper bounds of the actuators vector with respect to the prior configuration respectively. Then the over-estimated end-effector workspace W pjq E is defined as

W pjq E " ! u pjq P R nu |u pjq P " u pjq , u pjq ı) (4.9)
with u pjq " A `δτ pj´1q `A´δ τ pj´1q `CK ´1pq pj´1q qF pq pj´1q q `upj´1q and u pjq " A `δτ pj´1q `A´δ τ pj´1q `CK ´1pq pj´1q qF pq pj´1q q `upj´1q .

Proof. Using the linear approximation defined in (4.8), the nodal displacement vector can be directly deduced as:

δq pjq " K ´1pq pj´1q q " Hpq pj´1q qδτ pj´1q `F pq pj´1q q ‰

With the defined matrix A, we can then use the above equation to write the end-effector displacement in function of the matrix A, which yields: u pjq " Cδq pjq `upj´1q " CK ´1pq pj´1q q " Hpq pj´1q qδτ pj´1q `F pq pj´1q q ‰ `upj´1q " CK ´1pq pj´1q qHpq pj´1q qδτ pj´1q `CK ´1pq pj´1q qF pq pj´1q q `upj´1q " Apq pj´1q qδτ pj´1q `CK ´1pq pj´1q qF pq pj´1q q `upj´1q Next, according to Lemma 3 in [START_REF] Zheng | Design of interval observer for a class of uncertain unobservable nonlinear systems[END_REF], the following over-estimation is obtained:

A `δτ pj´1q `A´δ τ pj´1q ď Aδτ ď A `δτ pj´1q `A´δ τ pj´1q
By simply defining u pjq " A `δτ pj´1q `A´δ τ pj´1q `CK ´1pq pj´1q qF pq pj´1q q `upj´1q and u pjq " A `δτ pj´1q `A´δ τ pj´1q `CK ´1pq pj´1q qF pq pj´1q q `upj´1q that leads to: u pjq ď u pjq ď u pjq Finally, we obtain the over-estimated workspace W pjq E defined by (4.9).

The third step of the forward-backward approach needs also to be adapted in order to treat the FEM case. Therefore, as explained in Step 3 of Section 4.3.3 for the PCS model, we follow the same procedure for the FEM model, and we determine the associated displacement vector of the end-effector as: qpjq E " ûpjq ´upj´1q

Consequently, we can calculate the equivalent actuators and displacement vector by solving the optimization problem formulated below:

rτ pjq , q pjq s " arg min

τ,q }q E ´q pjq E } 2 2 s.t. τ P T Kpqqq " Hpqqτ `F pqq q E " Cq (4.10)
Solving the backward mechanism (4.10) yields the feasible actuators vector τ pjq with its related feasible displacement vector q pjq , which will eventually drive the soft robot's end-effector to reach ûpjq . Finally, we deduce the end-effector position associated with the optimal displacement vector q pjq as u pjq " Cq pjq `upj´1q .

The algorithm of the forward-backward approach for the FEM case consists of the same main steps that was described for the PCS model in Table 3. However, the functions of Steps 1 and 3 need to be adapted for the FEM case as explained above.

Validation and Analysis

This section presents the comparative results between the forward and forward-backward approaches to the same soft robots' configurations presented in the scenarios of Section 3.3.2 in Chapter 3.

In the following scenarios, the workspace obtained from the forward approach is depicted by blue-colored points. We propose to discretize the end-effector space with a discretization step size of s E d " 0.005 Meter (unit).

Scenario 1: trunk-like soft robot actuated by two cables.

In the first scenario (as shown by Fig. 4.14a), we consider the same configuration of the first scenario in Section 3.3.2.1 from Chapter 3. The applied force is of the following range: T p0q " T " r0 100s ˆr0 100s.

For this scenario, the forward-backward approach took 37 complete iterations and found a total of 635 feasible points in the workspace with a time computation complexity of 757 seconds. Those obtained feasible points are depicted by red colored points in Figs. 4.14b, 4.14c and 4.14d), and the union of all small neighborhood of those feasible points gives the estimation of workspace W E , which is represented by the gray zone with green contour. It is clear to see from the three-view drawing that the workspace estimation result obtained via forward-backward method coincides with that obtained via forward approach.

Moreover, for a complete 37 iterations it is expected to have W E d " 6 37 " 6.18 ˆ10 28 points to explore. However, what we actually explored is a total of W E d " 1680 points, and a final total of W Es " 635 feasible points. This significant reduction was partially due to the elimination process of inappropriate configurations from Step. 2 (See Fig. 4.3), but mainly due to the elimination of all feasible configurations that were already explored during Step. 5 (See Fig. In this scenario (as shown by Fig. 4.15a), we consider the same configuration of the second scenario in Section 3.3.2.2 from Chapter 3. The applied force is of the following range: T p0q " T " r0 100s ˆr0 100s ˆr0 150s. Next, we apply the forward-backward interval analysis approach using the equations presented in Section 4.4. For this scenario, the forward-backward approach took 44 complete iterations and calculated a total of 1941 feasible points in the workspace with a time computation complexity of 2063 seconds. Those achieved feasible points are illustrated with red colored points in Figs. 4.15b, 4.15c and 4.15d), and the union of all small neighborhood of those feasible points yield the estimated workspace W E , which is represented by the gray zone within the green contour. Using the three-view illustration, we can clearly observe that the estimated workspace calculated via the proposed forward-backward approach conforms with that obtained through the forward approach.

Furthermore, for a complete 44 iterations it is expected to get W E d " 6 44 " 1.17 ˆ10 34 points to investigate. Yet, the total of points that we actually explored is W E d " 5466, with a final total of W Es " 1941 feasible points. This substantial reduction was slightly due to Step. 2 (refer to Fig. 4.3) elimination process of irrelevant configurations, but considerably due to the exclusion of already explored feasible configurations throughout Step. 5 (See Fig. In this scenario (as shown by Fig. 4.16a), we consider the same configuration of the last scenario in Section 3.3.2.3 from Chapter 3. The applied force is of the following range: T p0q " T " r0 100s ˆr0 100s ˆr0 150s ˆr0 150s.

Next, we apply the proposed forward-backward interval analysis approach by following the procedure presented in Section 4.4. For this scenario, the forward-backward approach took 42 complete iterations and computed a total 3667 feasible points in the workspace with a time computation complexity of 3733 seconds. Those obtained feasible points are depicted by red colored points in Figs. 4.16b, 4.16c and 4.16d, and the union of all small neighborhood of those feasible points gives the estimation of the workspace W E , which is represented by the gray zone with green contour. It is clear to see from the three-view drawing that the workspace estimation result obtained via forward-backward method coincides with that obtained via forward approach.

In addition, for a complete 42 iterations it is expected to have W E d " 6 42 " 4.81 ˆ10 32 points to explore. However, what we actually explored is a total of W E d " 10626 points, and a final total of W Es " 3667 feasible points. This significant reduction was partially due to the elimination process of inappropriate configurations from Step. 2 (See Fig. 4.3), but mainly due to the elimination of all feasible configurations that were already explored during Step. 5 (See Fig.

Notes on computational efficiency

In the following table (see Table 4.2) and illustrations (see Fig. 4.17a and Fig. 4.17b), we summarize the operations and time computation complexity of each scenario. .17: Operations complexity and computational time of the forward and forward-backward approaches for the workspace estimation of the investigated scenarios using the FEM model.

Conclusion

In this chapter, we have proposed an interval analysis approach, named the forward-backward approach after its procedure, which consists of exploring the whole feasible configurations of the workspace by discretizing the end-effector's domain.

The proposed approach was implemented for both the PCS and the FEM models, where we have shown its effectiveness in reducing the computation complexity and computational time necessary to estimate the workspace of soft robots in contrast to the forward approach which explodes exponentially when increasing the dimension of the actuators.

On the other hand, since this approach consists of determining all the possible configurations that are feasible, it is useful in identifying interior and exterior information about the workspace. However, this approach is exhaustive in the sense that it explores all feasible configurations of the workspace instead of mapping only its interior and exterior boundaries.

Chapter 5

Continuation Approach for Workspace Estimation: PCS case

Introduction

Due to the limitations encountered by the optimization approach in Chapter 3 that cannot provide knowledge of the interior configurations and boundaries of the workspace of soft robots, and also due to the limitations of the interval analysis based forward-backward approach in Chapter 4 is exhaustive by exploring the whole workspace, it is logical to provide an alternative methodology allowing us to determine only but all boundaries of soft robots workspaces, both interior and exterior.

Such a method is called the continuation approach [START_REF] Amehri | Workspace boundary estimation for soft manipulators using a continuation approach[END_REF] that is based on the theory of bifurcation. This method is founded on the mathematical theory of bifurcation that was originally introduced by the French mathematician Henri Poincaré [START_REF] Poincaré | Sur l'équilibre d'une masse fluide animée d'un mouvement de rotation[END_REF], and it consists of studying the solutions of nonlinear differential equations, in addition to the study of possible variations in the topological structure of a vector field family. Subsequently, Henri Poincaré classified various types of bifurcation points [START_REF] Poincaré | New methods of celestial mechanics[END_REF]. Accordingly, this approach takes into account possible bifurcation behavior in continuation calculation while mapping the interior and exterior boundaries of soft robots workspaces. The continuation method [START_REF] Edward | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF] applies the idea that the workspace boundary (noted as BW E ) of a soft robot can be extracted from the set of its output singularities. Starting from an initial point on the boundary (i.e., point A as shown in Fig. 5.1), this method firstly classifies the status of the boundary point, based on which it calculates the associated tangent vectors, and then numerically predicts its evolution tendency which gives point B. Note that this predicted point B might not be still on BW E , thus a correction procedure is applied to seek a point of BW E by minimizing the distance to B which yields point C P BW E . Such a procedure is for the purpose of mapping the point B to one-dimensional solution curves.

Briefly speaking, the proposed continuation approach can be summarized in the following procedure:

1. Formulating the analytical criteria of the workspace boundary.

Finding an initial point on the workspace boundary.

3. Calculating tangent vectors to continuation curves at boundary points.

Mapping one-dimensional solution curves.

In the following, we will outline the realization of the continuation approach for the PCS model and later discuss the possibility of its extension to the FEM model.

Implementation of the Methodology

In this section, we present the specific details of the continuation approach applied to the PCS model.

Formulating the Analytical Criteria of the Workspace Boundary

Intuitively, a workspace boundary might be caused by the physical constraint (such as bounded length) of the soft robots, or by the boundedness of actuators, i.e., τ k P rτ k , τ k s where τ k and τ k represent respectively the lower and upper bound of the kth actuator. In order to take into account actuator boundedness, we introduce two intermediate variables s 2k´1 P R and s 2k P R satisfying τ k ´τk `1 2 s 2 2k´1 " 0 (equivalent to τ k ě τ k ), and ´τ k `τk `1 2 s 2 2k " 0 (equivalent to τ k ď τ k ). The purpose of the introduction of the variable s P R 2nτ is to treat the boundedness of the actuators in the same way as the mathematical model (e.g., (2.48)) via algebraic equations.

Following the above thoughts, in order to formulate the analytical criteria of the workspace boundary by taking into account the saturation of actuators, let us introduce the following generalized coordinates vector:

Q " " u T , τ T , s T ‰ T " " u T , Z T ‰ T (5.1)
where the actuators vector τ P R nτ and the intermediate variable vector s P R ns are combined in the vector Z "

" τ T , s T ‰ T P R n Z , Q P R n Q with n Q " n u `nZ , n Z " n τ `
ns , and n s " 2n τ . With the introduced intermediate variable vector s, the bounded inequality constraints of each actuator can be then transformed into equality constraints, and this enables us to achieve a generalized formulation of the studied soft robot with the following mechanical constraints ΨpQq:

ΨpQq " Ψpu, Zq " » - - - - - - - - - - - - - - - - - - - - - - - - - u ´E1 g `qpτ q ˘ET 2 τ 1 ´τ1 `1 2 s 2 1 ´τ 1 `τ1 `1 2 s 2 2 . . . τ k ´τk `1 2 s 2 2k´1 ´τ k `τk `1 2 s 2 2k
. . . where Ψpu, Zq P R m is the generalized mechanical constraints, with m " n u `ns , and the generalized strain vector q corresponding to the actuators vector τ is calculated via (2.51).

τ
Remark 6. Note that for a given actuator configuration τ " rτ 1 , ¨¨¨, τ nτ s T with τ k " rτ k , τ k s for 1 ď k ď n τ , its associated intermediate variables s 2k´1 and s 2k can be uniquely determined by:

s 2k´1 " a 2|τ k ´τ k |, s 2k " a 2|τ k ´τk | (5.3)
i.e., Z " " τ T , s T ‰ T can be uniquely determined by τ . Furthermore, by using numerical approach to solve Ψpu, Zq " 0, like Newton method, we can then obtain the associated value of u.

Considering the above formulation of the generalized mechanical constraints (5.2), the generalized form of a soft robot W E can therefore be explicitly reformulated as follows:

W E " t u P R nu | Ψpu, Zq " 0, for some Z u (5.4)
Now, consider a soft robot configuration in which its end-effector position lies on the workspace boundary of (5.4). Physically, this situation is equivalent to a singular behavior in which the soft robot end-effector cannot move in certain directions. Mathematically, we can deduce that the soft robot's workspace boundary is a subset of its workspace at which the sub-Jacobian matrix of the generalized mechanical constraints Ψpu, Zq with respect to the actuators vector τ and the intermediate variable s, denoted as Ψ Z pu, Zq P R mˆn Z , is row-rank deficient [START_REF] Ej Haug | Dextrous workspaces of manipulators. i. analytical criteria[END_REF], and it can be defined by the following: Definition 3. The workspace boundary BW E of a soft robot is a subset of its workspace W E , at which the sub-Jacobian Ψ Z pu, Zq P R mˆn Z is row-rank deficient, i.e.,

BW E " t u P W E | rank Ψ Z pu, Zq ă m, for Z P R n Z with Ψpu, Zq " 0 u
Equivalently, the row-rank deficient of Ψ Z pu, Zq implies that there exists an orthonormal vector γ P R m forming the nullspace of Ψ T Z , i.e., Ψ T Z pu, Zqγ " 0. Thus, an analytical condition [START_REF] Jo | Workspace analysis of multibody mechanical systems using continuation methods[END_REF] for a soft robot's workspace boundary can be reformulated as:

BW E " t u P W E | Ψ T Z pu, Zqγ " 0, γ T γ " 1, Ψpu, Zq " 0u (5.5)
In accordance with (5.5), a composite vector is introduced:

Y " " u T , Z T , γ T ‰ T " " Q T , γ T ‰ T (5.6)
where

Y P R n Y , n Y " n u `nZ `m " n Q `m.
Since boundary points must be characterized by the generalized vector Q, that satisfy (5.5), then, stated in terms of the composite vector Y , points on BW E can be characterized by the following composite set of equations:

GpY q " » - - Ψpu, Zq Ψ T Z pu, Zqγ 1 2 pγ T γ ´1q fi ffi fl " 0 (5.7)
Clearly, the above equation is a system of m `nZ `1 equations with n u `nZ `m variables. Thus the dimension of its solutions, equivalently the dimension of BW E , is equal to n u ´1.

Accordingly, when the end-effector position is two-dimensional u P R 2 , then BW E is a one-dimensional curve. On the other hand, in the case where soft robots are used in spatial applications, the output vector is three-dimensional u P R 3 , implying that BW E is two-dimensional.

Remark 7. In spatial applications, this thesis proposes an approach to calculating sets that are of dimension two by mapping the achievable one-dimensional BW E curves on the Oxy, Oyz, and Oxz plans of the workspace boundary solution set (5.7).

Moreover, a practical approach [START_REF] Jo | Workspace analysis of multibody mechanical systems using continuation methods[END_REF] to map two-dimensional workspace boundary sets has been proposed, where the purpose is to define a linear relation within the output coordinates by introducing a plan intersecting with the workspace boundary in a one-dimensional curve.

Finding an Initial Point on the Workspace Boundary

In order to map the workspace boundary, we need first to figure out how to determine an initial boundary point u b0 . For concrete applications, it is a nontrivial task to find a point on the workspace boundary BW E (as shown in Fig. 5.2). Assume that a radiating point u 0 is selected, and that u 0 is exterior to W E . It is now proposed that, consistent with Definition 3, a boundary point u b0 in a particular direction emanating from the radiating point u 0 (as depicted by Fig. 5.2), can be determined by solving the following constrained optimization problem:

τ b0 " min τ upqq ´u0 2 2 s.t. τ P T K q " Hpqqτ `F pqq upqq " E 1 g `qpτ q ˘ET 2 (5.8)
where τ b0 is the optimal and feasible actuators vector (since it satisfies the bounded constraint) for the purpose of minimizing the nonlinear cost function f pqq " upqq ´u0 2 2 , which represents the distance between u 0 and upqq.

From the solution τ b0 , and according to Proposition 1, we calculate the corresponding generalized strain vector q b0 (2.51). Then, we compute the strain twist ξb0

i related to the strain vector q bj i of each body i using (2.49b), and based on which we derive the configuration matrix gpq b0 q associated to the generalized strain vector q b0 from (2.49a). Finally, we deduce the boundary point position u b0 using (2.48b).

Using the values of u b0 , τ b0 , and s b0 , we establish the corresponding generalized vector Q b0 using (5.2). Since the obtained solution verifies the workspace boundary definition (see Definition 3), what remains now is to find a vector γ b0 that satisfies (5.5), which can be done by calculating the nullspace of Ψ T Z pQ b0 q. Finally, the obtained vector γ b0 with the generalized vector Q b0 yields the desired initial composite vector Y b0 (5.6) from which the algorithm that consists of mapping the workspace boundary solution curves (5.7) will start.

Regarding the implementation of the methodology, we need to calculate the sub-Jacobian Ψ Z pu, Zq, which can be formulated as follows:

Ψ Z pu, Zq " " ´∇τ " upqq ‰ 0 nuˆns ∆ ˘1 ∆ s  (5.9)
with ∇ τ " upqq ‰ " ∇ τ pqqJ T pqqr0 3 Rpqqs T being already established in Theorem 1, and ∇ τ pqq " H T pqq " B ´1pτ, qq ‰ T being defined in Assumption 1, with Bpτ, q, L, αq " KpL, αq´H q pq, L, αqpI n b τ q ´Fq pq, L, αq.

The matrices ∆ ˘1 and ∆ s are equal to the following: The initial radiating point u 0 in 5.2, may be generated using the same approach explained in Section 3.2.2 of Chapter 3.

∆ ˘1 " » - - - - - - ´1 ¨¨¨0 `1 ¨¨¨0 . . . . . . . . . 0 ¨¨¨´1 0 ¨¨¨`1 fi ffi ffi ffi ffi ffi fl , ∆ s " » - - s 1 . . .
Alternatively, we could also directly select a limit point, i.e., a point for which the related actuators vector is equal to a combination set of upper and lower values of each actuator, as an initial boundary point. Such a point is always located on BW E , according to Definition 3, since its associated sub-Jacobian is row-rank deficient due to the actuators saturation (because one or multiple variables s will be zero when some actuators reach their upper or lower bounds, and thus the matrix ∆ s in (5.10) will cause the sub-Jacobian Ψ Z to be row-rank deficient).

Calculating Tangent Vectors to Continuation Curves at Boundary Points

As pointed out by Fig. 5.1, mapping BW E consists of finding tangent vectors hpY q to solution curves at each boundary point met during the continuation process. But before tangent calculation, we need to identify the status of each boundary point met along the mapping of the workspace boundaries.

Actually, there exist two different types of boundary points [START_REF] Edward | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF], a regular boundary point, and a bifurcation point. A boundary point can be classified as regular [START_REF] Edward | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF] if the Jacobian matrix of the composite set that characterizes candidate points on the workspace boundary GpY q defined in (5.7), denoted as G Y pY q P R n Y ´1ˆn Y , has full row-rank, and Ψ Z pu, Zq is row-rank deficient 1.

On the other hand, a bifurcation point is where at least two solution curves can cross, and it can be classified into two different categories, simple and multiple [START_REF] Edward | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF]. A bifurcation point may be classified as simple [START_REF] Wang | Dextrous workspaces of manipulators, part 2: Computational methods[END_REF] if both G Y pY q and Ψ Z pu, Zq are row-rank deficient 1. Finally, a multiple bifurcation point [START_REF] Charles Qiu | Dextrous workspaces of manipulators, part iii: Calculation of continuation curves at bifurcation points[END_REF] can be identified by two conditions, the first is that G Y pY q is row-rank deficient of degree 1, and the second is that Ψ Z pu, Zq is row-rank deficient of degree 2.

These possible cases are summarized in Table 5.1.

Boundary point

Row-Rank Ψ Z pu, Zq G Y pY q Regular m ´1 n Y ´1 Bifurcation Simple m ´1 n Y ´2 Multiple m ´2 n Y ´2
Table 5.1: Classification of the possible categories of boundary points.

In the following, we will explain in detail how, in each case, tangent vectors to continuation curves may be calculated.

Regular boundary points

Within the solution curve (5.7), a regular point u br characterized by the composite vector Y br (as shown by Fig. 5.3), where at most a single solution curve is passing through, may be met. At such a point, the row-rank deficiency of Ψ Z is 1, and G Y pY br q has full row-rank. At such a point, a tangent vector hpY br q (as illustrated by Fig. 5.3) to the set (5.7) may be calculated based on the following theorem [START_REF] Eugene | Introduction to numerical continuation methods[END_REF][START_REF] Werner C Rheinboldt | Numerical analysis of parametrized nonlinear equations[END_REF]:

Theorem 4. Let G: S Ă R n Y Þ Ñ R n Y ´1, be of class C r , r ě 1, on an open set S Ă R n Y . Thus, at a point Y br , if G Y
has full row-rank, then a unit tangent vector hpY br q P R n Y is uniquely defined by: G Y pY br qhpY br q " 0 h T pY br qhpY br q " 1 G Y pY br q h T pY br q ą 0 (5.11)

The tangent vector hpY br q determined by (5.11) can be calculated at each point along the workspace boundary solution set (5.7) as long as G Y pY br q has full row-rank.

However, difficulties may arise when multiple solution curves cross at bifurcation points of (5.7). According to Table 5.1, G Y is row-rank deficient at such configurations, and (5.11) is unable to determine a unique tangent vector. Therefore, new techniques must be adopted [START_REF] Charles Qiu | Dextrous workspaces of manipulators, part iii: Calculation of continuation curves at bifurcation points[END_REF][START_REF] Wang | Dextrous workspaces of manipulators, part 2: Computational methods[END_REF] in order to compute tangent vectors hpY q to continuation curves at both simple and multiple bifurcation points.

Bifurcation points -Simple bifurcation point

While proceeding along the workspace boundary solution set defined by (5.7), a point u bs characterized by its composite vector Y bs (as depicted by Fig. 5.4) may be found at which the Jacobian G Y pY bs q becomes row-rank deficient of degree 1, and the rank deficiency of Ψ Z is 1. Such a point is classified as a simple bifurcation point [START_REF] Wang | Dextrous workspaces of manipulators, part 2: Computational methods[END_REF], where a pair of solution curves (See Fig. 5.4) characterized by a pair of tangent vectors hpY bs q satisfying (5.11) are traversing it.

For the purpose of identifying the nature of a singular behavior at such a point Y bs , we perform an elementary row decomposition, e.g. singular value decomposition (SVD) or QR decomposition, to simplify G Y pY bs q into the following form:

EG Y pY bs q " " F Y pY bs q 0  (5.12)
where E is an elementary non-singular matrix, and F Y pY bs q is an pn Y ´2q ˆnY matrix with full row-rank. Applying elementary row operation to the workspace boundary solution set (5.7) using E, we get the following equivalent equation: q GpY bs q " " FpY bs q fpY bs q  " EGpY bs q " 0 (5.13) Using (5.12) and (5.13), we observe that F Y pY bs q is the Jacobian matrix of FpY bs q, and the Jacobian of the function fpY bs q is f Y pY bs q " 0. Since F Y pY bs q is full row-rank, i.e., rank F Y pY bs q " n Y ´2, then the nullspace of F Y pY bs q (equivalently, the nullspace of G Y pY bs q) is of dimension 2, which generates two different non-zero solutions of the following equation:

F Y pY bs qζ i " 0, i " 1, 2 (5.14)
that can be orthonormalized;

ζ i ζ j " δ ij , i, j " 1, 2
, where δ is the Dirac delta function.

The purpose now is to determine the pair of tangent vectors hpY bs q to continuation curves at Y bs by writing it in the form of hpY bs q " ř j α j ζ j . With ζ being calculated via (5.14), the remaining question now is about how to compute the vector α.

Using Taylor's expansion, the tangent vectors of the solution curve can be calculated. Taylor's expansion for functions of n Y variables is stated in [START_REF] Barnett | Calculus of several variables and differentiable manifolds[END_REF] as follows:

Theorem 5. According to [START_REF] Barnett | Calculus of several variables and differentiable manifolds[END_REF]

, let G: S Ă R n Y Þ Ñ R n Y ´1, n Y ě 2, be of class C r , r ě 1, on the open set S Ă R n Y .
Then, the quadratic Taylor approximation is defined by: q GpY bs `∆Y q « q GpY bs q `«

F Y pY bs qα i ζ i 1 2 pα j ζ j qf Y Y pY bs qpα k ζ l q ff (5.15)
where the Hessian matrix of fpY q, evaluated at Y bs , is the matrix:

f Y Y pY bs q " « BfpY q BY i Y j ff Y "Y bs (5.16)
Since Y bs `∆Y and Y bs must satisfy the boundary workspace solution set (5.7), then we have: q GpY bs `∆Y q " q GpY bs q " 0 Taking that into consideration, (5.15) becomes:

« F Y pY bs qα i ζ i 1 2 pα j ζ j qf Y Y pY bs qpα k ζ k q ff " 0 (5.17)
Noting that, for all α i , the upper sub-vector term F Y pY bs qα i ζ i " 0, since ζ i represents the nullspace of F Y pY bs q as stated in (5.14). Finally, (5.17) is reduced to the following:

α T " ζ i f Y Y pY bs qζ j ‰ α " 0 (5.18)
where α " rα 1 , α 2 s T . Using elementary row decomposition, e.g. eigenvalue decomposition, the above quadratic equation may be solved. Another way is to transform (5.18) to a quadratic equation in terms of α 1 and α 2 and solving it using the quadratic formula. No matter which way is applied, solving (5.18) yields two normalized solution vectors α 1 and α 2 . Finally, with ζ 1 and ζ 2 being calculated via (5.14), we obtain a pair of tangent vectors hpY bs q " rh 1 pY bs q T h 2 pY bs q T s T depicted by Fig. 5.4, and formulated by the following:

h 1 pY bs q " α 1 1 ζ 1 `α1 2 ζ 2 h 2 pY bs q " α 2 1 ζ 1 `α2 2 ζ 2
(5.19)

Figure 5.4: Simple bifurcation point -tangent vector.

Bifurcation points -Multiple bifurcation point

During the continuation calculation, a point u bm characterized by the composite vector Y bm (as depicted by Fig. 5.5) may be encountered at which the Jacobian Ψ Z becomes row-rank deficient of degree 2, and the rank deficiency of G Y pY bm q is 1. Such a point is classified as a multiple bifurcation point [START_REF] Charles Qiu | Dextrous workspaces of manipulators, part iii: Calculation of continuation curves at bifurcation points[END_REF], where numerous branches may occur (See Fig. 5.5) depending on the design of the soft robot. This thesis uses a numerical method [START_REF] Charles Qiu | Dextrous workspaces of manipulators, part iii: Calculation of continuation curves at bifurcation points[END_REF] that consists of finding, for all branches intersecting at Q bm , all vectors γ such that G Y pQ bm , γq for each branch is row-rank deficient.

The adopted strategy can be separated into two major steps:

[

Step 1] Calculating all solutions for γ at a multiple bifurcation point.

[

Step 2] Calculating all tangent vectors hpQ bm , γq to continuation curves. Step 1: Calculating all solutions for γ at a multiple bifurcation point In the first step, vector γ for which G Y pQ bm , γq is row-rank deficient will be determined. Since the row-rank deficiency of Ψ Z is 2, thus, we have 2 orthonormal solutions γ 1 and γ 2 such that:

Ψ T Z γ i " 0 γ T i γ j " δ ij , i, j " 1, 2 (5.20)
where γ can be represented as a linear combination of γ 1 and γ 2 as:

γ " γpβq " tanpβqγ 1 `γ2 (5.21)
with β P r0 πs. Hence, the purpose of this step is to solve the following problem:

Problem 1. Calculate all γ for which each corresponding G Y pQ bm , γq is row-rank deficient, where:

G Y pQ bm , γq " » - - Ψ u Ψ τ Ψ s 0 0 " Ψ T Z γ ı τ " Ψ T Z γ ı s Ψ T Z 0 0 0 γ T fi ffi fl (5.22)
Since the number of γ is finite, it can be assumed that neither γ 1 nor γ 2 are solutions. Therefore, the solution of (5.21) must verify the condition sinpβq cospβq ‰ 0.

Lemma 1. G Y pQ bm , γq is row-rank deficient if and only if there exists a non-zero vector c " rc T 1 , c T 2 , c 3 s T , where c 1 P R m , c 2 P R n Z , and c 1 P R:

c T G Y pQ bm , γq " 0 (5.23)
By developing (5.23) through (5.22), we obtain the following set of equations:

c T 1 Ψ τ `cT 2 " Ψ T Z γ ı τ " 0 (5.24a) c T 1 Ψ s `cT 2 " Ψ T Z γ ı s " 0 (5.24b) c T 2 Ψ T Z `c3 γ T " 0 (5.24c)
Introduce now µ as: µ " µpβq " tanpβq, then, the problem may be equivalently reformulated as follows:

Problem 2. Assume that neither G Y pQ bm , γ 1 q nor G Y pQ bm , γ 2 q is row-rank deficient, the goal now would be to find all µpβq such that G Y pQ bm , γpβqq is row-rank deficient.

A necessary condition for G Y pQ bm , γpβqq has already been outlined by the set of equations (5.24a), (5.24b), and (5.24c). Multiplying (5.24c) by γ, we obtain:

c T 2 Ψ T Z γ `c3 γ T γ " 0 Using (5.

20), we deduce:

Ψ T Z γ " 0 ùñ c 3 " 0 Taking that into consideration, (5.24c) finally becomes:

Ψ Z c 2 " 0
And since the row-rank deficiency of Ψ Z : R n Z Þ Ñ R m is 2, then there exists an orthonormal basis σ of the nullspace of Ψ T Z such that:

c 2 " d 1 σ 1 `d2 σ 2 `¨¨¨`d nσ σ nσ " σd (5.25)
where σ " tσ 1 , σ 2 , ¨¨¨, σ nσ u, d " td 1 , d 2 , ¨¨¨, d nσ u, and n σ " n Z ´m `2. From (5.24a), and (5.24b), we obtain the following:

c T 1 " Ψ τ , Ψ s ı " ´cT 2 « " Ψ T Z γ ı τ , " Ψ T Z γ ı s ff
For the purpose of formulating c 1 in function of c 2 , and since " Ψ τ , Ψ s ‰ is right invertible, we can write the above equation in the following form:

c T 1 " c T 2 Mpγq (5.26)
with:

Mpγq " ´«" Ψ T Z γ ı τ , " Ψ T Z γ ı s ff " Ψ τ , Ψ s ı `(5.27)
After, by substituting c T 1 from (5.26) into (5.24b), we get:

c T 2 ˜MpγqΨ s `"Ψ T Z γ ı s ¸" 0
Next, we replace (5.25) into the above equation, we obtain:

d T Hpγpβqq " 0 with Hpγpβqq " σ T ˜MpγqΨ s `"Ψ T Z γ ı s ¸(5.28)
As stated in Lemma 2.3 of [START_REF] Charles Qiu | Dextrous workspaces of manipulators, part iii: Calculation of continuation curves at bifurcation points[END_REF], G Y pQ bm , γq is row-rank deficient if and only if Hpγpβqq is row-rank deficient. By substituting (5.21) into (5.28), we achieve:

Hpγpβqq " tanpβqA `B (5.29)
where:

A " σ T ˜Mpγ 1 qΨ s `"Ψ T Z γ 1 ı s ¸(5.30a) B " σ T ˜Mpγ 2 qΨ s `"Ψ T Z γ 2 ı s ¸(5.30b)
Based on Lemma 2.4 in [START_REF] Charles Qiu | Dextrous workspaces of manipulators, part iii: Calculation of continuation curves at bifurcation points[END_REF], if both G Y pQ bm , γ 1 q and G Y pQ bm , γ 2 q have full row-rank, then A is right invertible. Equivalently, according to Theorem 2.5 in [START_REF] Charles Qiu | Dextrous workspaces of manipulators, part iii: Calculation of continuation curves at bifurcation points[END_REF], G Y pQ bm , γq is row-rank deficient if only and if µ is a real eigenvalue of ´BA `.

Step 2: Calculating all tangent vectors hpQ bm , γq to continuation curves. In this step, we will describe how to calculate the tangent vectors hpQ bm , γq to all solution curves passing through a multiple bifurcation point.

Consistent with Lemma 2.4 in [START_REF] Charles Qiu | Dextrous workspaces of manipulators, part iii: Calculation of continuation curves at bifurcation points[END_REF], denote χ as the eigenvector corresponding to the eigenvalue µ of the matrix ´BA `. By definition, we obtain:

´BA `χ " µχ (5.31)
Let us define ρ 2 P R ns as the following [START_REF] Charles Qiu | Dextrous workspaces of manipulators, part iii: Calculation of continuation curves at bifurcation points[END_REF]:

ρ 2 " A `χ (5.32)
From the above equation, we can easily deduce that χ " Aρ 2 . Using this, and substituting (5.32) in (5.31), we obtain:

´Bρ 2 " µAρ 2
which can be reformulated as: pµA `Bqρ 2 " 0 Since µ " tanpβq, we get: ptanpβqA `Bqρ 2 " 0 Using (5.29), then (5.28), the above equation can finally be rewritten as follows:

σ T S " 0 (5.33) with S " ˜MpγqΨ s `"Ψ T Z γ ı s ¸ρ2 .
Since σ is an orthonormal basis of the nullspace of Ψ T Z , and from (5.33) S is orthogonal to σ, then S represents the column space of Ψ T Z [START_REF] Evar | Linear algebra and matrix theory[END_REF]. Therefore, there exists a unique vector ρ 3 P R m in the column space of Ψ T Z , for which the following is deduced:

S " ´ΨT Z ρ 3 (5.34)
Given that ρ 3 and γ are respectively the column space and the nullspace of Ψ T Z . Therefore, using the first formula of S from (5.33), its second formula from (5.34), we obtain the following set of equations:

˜MpγqΨ s `"Ψ T Z γ ı s ¸ρ2 " ´ΨT Z ρ 3 (5.35a) 0 " γ T ρ 3 (5.35b)
which finally allows us to obtain the final form of ρ 3 P R m as follows:

ρ 3 " " ´ΨT Z γ T  `«MpγqΨ s `"Ψ T Z γ ı s 0 ff ρ 2 (5.36)
Next, we define ρ 1 P R nu`nτ by the following:

ρ 1 " ´"Ψ u , Ψ τ ı `Ψs ρ 2 (5.37)
Substituting both Mpγq and ρ 1 respectively from (5.26) and (5.37) into (5.35a), and using (5.35b), we can deduce the following system of equations:

" Ψ u , Ψ τ ı ρ 1 `Ψs ρ 2 " 0 (5.38a) « 0, " Ψ T Z γ ı τ ff ρ 1 `"Ψ T Z γ ı s ρ 2 `ΨT Z ρ 3 " 0 (5.38b) γ T ρ 3 " 0 (5.38c)
Eventually, the above set can be reformulated into the following form:

0 " » - - - - - " Ψ u , Ψ τ ı Ψ s 0 « 0, " Ψ T Z γ ı τ ff " Ψ T Z γ ı s Ψ T Z 0 0 γ T fi ffi ffi ffi ffi fl » - ρ 1 ρ 2 ρ 3 fi fl 0 " G Y pQ bm , γq hpQ bm , γq (5.39)
Finally, it is clear that:

hpQ bm , γq " rρ T 1 , ρ T 2 , ρ T 3 s T (5.40)
is the tangent vector of the continuation curve corresponding to its associated eigenvalue µ and eigenvector χ.

Mapping One-Dimensional Solution Curves

With tangent vectors being calculated for each possible boundary point situation (as shown in Fig. 5.3, Fig. 5.4, and Fig. 5.5), the final step of the continuation method is to map the solution curves of (5.7). Given a boundary point u b0 characterized by its composite vector Y b0 , where Y b0 is a solution of the workspace boundary set (5.7), i.e. GpY b0 q " 0, and its tangent vector hpY b0 q being calculated depending on its category (regular boundary point, simple or multiple bifurcation point), a ray is emanated from Y b0 , with an infinitesimally small discretization precision ε, on which a prediction Y In this thesis, we propose to calculate Y b (as depicted in Fig. 5.6) using a constrained optimization problem, subject to nonlinear constraints, and formulated by the following:

Y b " min Y Y ´Y p s.t. GpY q " 0 (5.42)
where GpY q defined in (5.7) to characterize the workspace boundary set. With Y b being computed, we proceed by repeating the same process for each boundary point met during the continuation calculation until all solution curves of the workspace boundary are mapped.

An equivalent technique to (5.42) presented in [START_REF] Eugene | Introduction to numerical continuation methods[END_REF] has also been proven to be an efficient approach to find the next possible solution on the continuation curve.

The stop condition can be evaluated by verifying whether all limit points have been explored. When it is the case, then the stop condition is true, and the continuation process is terminated by yielding the inner and outer boundaries of the soft robot. If not, there might exist some isolated boundaries which cannot be reached from the chosen initial point. In this situation, we can choose another combination set of upper and lower values of each actuator, which has not yet been explored by the former process, as a new initial point and restart the same process.

In following, we present a brief algorithm (see Algorithm 4) explaining the major steps of the proposed continuation approach in order to map the interior and exterior boundaries of soft robots.

Algorithm 4 Calculate BW E Require: u 0 , T , ε BW E Ð ∅, stop Ð 0 Ź Initialization.
rτ b0 , u b0 s Ð solve p5.8q pu 0 , T q Ź Initial boundary point. Q b0 Ð solve (5.1) pu 0 , T q Ź Initial boundary configuration.

γ b0 Ð nullpΨ T Z pQ b0 qq Y b0 Ð rQ b0 T , γ b0 T
s T Ź Initial composite vector. while !stop do f lag Ð T able 5.1 pY b0 q Ź Boundary point class. switch f lag do Ź Tangent vector. case (Regular boundary point) hpY b0 q Ð solve (5.11) pQ b0 , Y b0 q case (Simple bifurcation point) hpY b0 q Ð solve (5.19) pQ b0 , Y b0 q case (Multiple bifurcation point)

hpY b0 q Ð solve (5.40) pQ b0 , Y b0 q Y p Ð Y b0 ˘εhpY b0 q Ź Predict. Y b Ð solve (5.42) pY b0 , Y p q Ź Correct. BW E Ð BW E ' Y b Ź Append. Y b0 Ð Y b Ź Update initial composite vector. stop Ð updatepBW E , Y b0 q Ź Update stop condition. end while

Validation and Analysis

We implement both the forward and the continuation approaches on planar and spatial soft robots to deduce the workspace estimation W E , and the workspace boundary BW E . Table 3.1 presents the numerical values of the investigated soft robot.

In the following scenarios, the workspace estimation obtained from the forward approach is depicted by gray-colored points. For the continuation approach, we suggest to proceed with a step ε " 0.01.

Scenario 1: planar 2-segments 2-tendons soft robot

In the first scenario, we consider the same configuration of the first scenario in Section 3.2.5.1 of Chapter 3. We choose a tension magnitude within T 1,2 " r0, 200sN that will allow us to apply a couple magnitude within C T 1 " r0, 6sN ¨m and C T 2 " r0, 4sN ¨m. Starting from an initial boundary configuration calculated via (5.8), which is a regular boundary point in this case (blue point on Fig. 5.7b), we use (5.11) to calculate its corresponding tangent vectors, which yields one tangent vector. Based on this, we calculate its associated prediction via (5.41), then deduce its corresponding solution using (5.42). Proceeding along the solution curve (5.7), we meet four bifurcation points (green points on Fig. 5.7b), between which regular boundary points (red points on Fig. 5.7b) are met. Finally, workspace boundary curves are mapped as shown in Fig. 5.7b.

For this scenario, the continuation approach computed a total of 87 boundary points and took only 9.54 seconds to map the workspace boundaries for this scenario. Figure 5.7: Scenario 1 -W E estimation via forward approach (gray points) and workspace boundaries BW E computed via the continuation approach (red curves and green bifurcation points).

Scenario 2: planar 2-segments 4-tendons soft robot

In the second scenario, we consider the soft robot of the first scenario and we add two symmetric tendons (See Fig. 5.8a). We choose a tension magnitude within T 1,2 " r0, 150sN that will allow us to apply a couple magnitude within C T 1,3 " r0, 4.5sN ¨m and C T 2,4 " r0, 3sN ¨m.

The numerical continuation algorithm starts from the blue point on Fig. 5.8b, which is a limit configuration in this case. Next, using (5.40) we obtain 4 distinct tangent vectors, each one corresponding to its associated solution curve since we have 4 that crosses through this point. After, we calculate the prediction corresponding to each tangent vector (5.41), and we obtain for each curve its associated solution (5.42). Regular boundary points (red points in Fig. 5.8b) are computed between each of the sixteen total bifurcation points, including the initial point (green points in Fig. 5.8b) found for this particular soft robot structure. Finally, interior and exterior boundaries to the planar two segments -four tendons soft robot's workspace are mapped and illustrated by Fig. 5.8b.

For this scenario, the continuation approach computed a total of 687 boundary points and took only 96.35 seconds to map the workspace boundaries for this scenario. Figure 5.8: Scenario 2 -W E estimation via forward approach (gray points) and workspace boundaries BW E computed via the continuation approach (red curves and green bifurcation points).

Scenario 3: planar 3-segments 3-tendons soft robot

In the third scenario (See Fig. 4.10a), we consider the same tendons routing of the third scenario in Section 3.2.5.3 from Chapter 3. We choose a tension magnitude within T 1,2,3 " r0, 100sN that will allow us to apply a couple magnitude within C T 1 " r0, 3sN ¨m, C T 2 " r0, 2sN ¨m, and C T 1 " r0, 1sN ¨m.

Interior and exterior workspace boundaries of this particular soft robot structure are mapped (as shown by Fig. 5.9b), using (5.8) we find an initial boundary configuration (blue point on Fig. 5.9b) that is a regular boundary point, based on which, we calculate its unique corresponding tangent vector through (5.11). Next, through (5.41) we predict the possible solution, and we determine its associated solution (5.42), proceeding all the way along solution curves we identify 8 bifurcation points (green points in Fig. 5.9b), with regular boundary points (magenta points in Fig. 5.9b) lying between them.

For this scenario, the continuation approach computed a total of 327 boundary points and took only 40.62 seconds to map the workspace boundaries for this scenario. Figure 5.9: Scenario 3: W E estimation via forward approach (gray points) and workspace boundaries BW E computed via the continuation approach (red curves and green bifurcation points).

Scenario 4: planar 3-segments 6-tendons soft robot

In this scenario, three additional tendons are routed in a fashion where they are symmetric to the tendons configuration in the third scenario (See Fig. 3.9a). We choose a tension magnitude within T 1,¨¨¨,6 " r0, 100sN that will allow us to apply a couple magnitude within C T 1,4

" r0, 3sN ¨m, C T 2,5 " r0, 2sN ¨m, and C T 3,6 " r0, 1sN ¨m.

Solution curves of the workspace boundary set (5.7) are mapped beginning from an initial boundary point corresponding to a limit point calculated via a random combination set of lower and upper values of each actuators vector bounds and is located in the interior boundary. From this initial point, relevant tangent vectors are calculated (5.40), yielding 6 possible directions to explore. Next, we predict the solution corresponding to each boundary curve from (5.41), and deduce their appropriate boundary configuration (5.42) (see Fig. 5.10b). Finally, solution sets are calculated and exterior and interior workspace boundaries are mapped (as shown by Fig.

5.10b).

For this scenario, the continuation approach computed a total of 4126 boundary points and took only 1096 seconds to map the workspace boundaries for this scenario. In this scenario, we consider a soft robot composed of two segments and actuated by four tendons (See Fig. 5.11a). The two lateral tendons are installed as in the first scenario. The two upward tendons are installed on the Oxy plan, and parallel to each other. The first tendon is fixed at the position p0, r 1 {2, 0q and extends along the first segment length L 1 . The second tendon is fixed at the position p0, r 2 {2, 0q and extends to the end-effector position. We choose the bounded tension magnitude T 1,2,3,4 " r0, 150sN enabling us to apply a couple magnitude C T 1,3 " r0, 4.5s, and C T 2,4 " r0, 3sN ¨m.

Interior and exterior workspace boundaries of this particular structure are mapped starting from an initial limit point. Next, using (5.40) we compute all possible tangent vectors for each plan Oxy, Oyz and Oxz. Based on this, we calculate a prediction corresponding to each tangent vector and we obtain each associated solution (5.42). The obtained result is depicted by red curves in a three-view drawing (Oxz, Oyx, and Oyz plans, in Figs. 5.11b, 5.11c, and 5.11d respectively).

For this scenario, the continuation approach computed a total of 688 boundary points and took only 107 seconds to map the workspace boundaries for this scenario. In this scenario, we consider a soft robot composed of two segments and actuated by six tendons (See Fig. 5.12a). There are four lateral symmetric tendons (as in Fig. 5.8a), and two upward tendons as in the above scenario. We choose the bounded tension magnitude T 1,¨¨¨,6 " r0, 100sN enabling us to apply a couple magnitude C T 1,2,3 " r0, 3s, and C T 4,5,6 " r0, 2sN ¨m.

The obtained result is depicted by red curves in a three-view drawing (Oxz, Oyx, and Oyz plans, in Figs. 5.12b, 5.12c, and 5.12d respectively). For this scenario, the continuation approach computed a total of 4107 boundary points and took only 980 seconds to map the workspace boundaries for this scenario. In this scenario, we consider a soft robot composed of three segments and actuated by four tendons (See Fig. 5.13a). The three lateral tendons are installed as in the third scenario (see Fig. 5.9a). The three upward tendons are installed on the Oxy plan, and are parallel to each other. The first tendon is fixed at the position p0, r 1 {2, 0q and extends along the first segment length L 1 . The second tendon is fixed at the position p0, r 2 {2, 0q and extends to pL 1 `L2 , r 2 {2, 0q, and the third tendon is fixed at the position p0, r 3 {2, 0q and extends all along to the end-effector position. We choose a tension magnitude within T 1,¨¨¨,6 " r0, 100sN that will allow us to apply a couple magnitude within C T 1,4 " r0, 3sN ¨m, C T 2,5 " r0, 2sN ¨m, and C T 3,6 " r0, 1sN ¨m.

The obtained result is depicted by red curves in a three-view drawing (Oxz, Oyx, and Oyz plans, in Figs. 5.13b, 5.13c, and 5.13d respectively). For this scenario, the continuation approach computed a total of 4122 boundary points and took only 1392 seconds to map the workspace boundaries for this scenario. In the final scenario (See Fig. 4.12a), we consider the same tendons routing of the last scenario in Section 4.3.6.5 from Chapter 4. We choose a tension magnitude within T 1,...,9 " r0, 50sN that will allow us to apply a couple magnitude within C T 1,4,7 " r0, 1.5sN ¨m, C T 2,5,8 " r0, 1sN ¨m, and C T 3,6,9

" r0, 0.5sN ¨m.

The obtained result is depicted by red curves in a three-view drawing (Oxz, Oyx, and Oyz plans, in Figs. 5.14b, 5.14c, and 5.14d respectively). For this scenario, the continuation approach computed a total of 7621 boundary points and took only 2674 seconds to map the workspace boundaries for this scenario. boundaries BW E computed via the continuation approach (red curves and green bifurcation points).

Notes on computational efficiency

In Table 5.2, and Figures 5.15a and 5.15b, we summarize the operations' complexity and computational time of each scenario. In this thesis, starting from one boundary point, the proposed continuation approach enables us to map the exterior and interior workspace boundaries whether such boundaries result from actuators saturation, or physical limitations such as the mechanism's length, self-collision, or elastic instabilities. However, as stated at the end of Section 5.2.4, isolated boundaries (see Fig. 5.16a) and voids (see Fig. 5.16b) might arise in complex soft robot configurations. In either case, both BW p1q E and BW p2q E verify the workspace boundary Definition 3. Suppose that in the first instance, the initial boundary point from which the continuation algorithm starts is on BW p1q E as shown in Fig. 5.16. In this case, using this point as a starting configuration of the continuation algorithm, we can surely map BW p1q E in the first process. Then, to map BW p2q E , we need to verify if all limit points have been explored by the first process, and since we have a second workspace boundary, it is evident that there exists at least one limit configuration that was not explored by the former process. Then, using this configuration as a new initial point, we can run a second continuation algorithm process allowing us to map BW p2q E . Another important issue concerns the characterization of different types of boundaries: exterior boundaries, interior boundaries, and voids. Firstly, we would like to emphasize that the proposed continuation approach contains necessary information to calculate the normal vector of each point on the obtained boundaries. In fact, for any boundary point Y b " rQ b , γ b s P BW E , its normal vector can be calculated as [START_REF] Edward | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF]: η b " ηpY b q " Ψ T u pQ b qγ b . To explain the main idea, we take Fig. E , with its normal vector η b E , therefore in the normal direction we can find a point (noted as u p E ) which will not belong to the workspace. In summary, this judgment can be formulated as the following condition:

Du p E " u b E `εη b E , s.t. u p E R W E
where ε P R is a pre-defined small value. Conversely, for a boundary point Y b I lying on the interior curve, with its normal vector η b I as shown by Fig. 5.16b, an interior boundary may be sorted by verifying: @u p I " u b I `εη b I , s.t. u p I P W E Finally, for Y b V lying on BW p2q E , with its normal vector η b V as illustrated by Fig. 5.16b, a void may be distinguished by verifying the following condition:

Du p V " u b V `εη b V , s.t. u p V R W E
Note that the above condition is identical to that of the exterior boundary. However, the fact that BW p2q E is enveloped by BW p1q E enables us to assert that BW p2q E is the void and BW p1q E is the exterior boundary.

Conclusion

The continuation approach allows the mapping of interior and exterior boundaries of soft robots with slender shapes using the PCS model.

However, it is important to state that this approach cannot provide information on all points in the workspace, but only on boundary points. Moreover, this approach was only applied to slender-shaped soft robots via the PCS model since this modeling method provides the analytical form of all necessary parameters. Thus, it is of interest to extend the applicability of this approach in order to also treat soft robots with general shape which can be modeled by FEM. However, since the FEM model is established differently from the PCS approach, some steps of the proposed continuation approach will need to be adapted in order to consider the FEM case (which will be outlined in Chapter 7).

113 Chapter 6 Design Optimization: PCS case

Introduction

In Chapter 3, Chapter 4, and Chapter 5, we have proposed different approaches in order to answer the direct problem of soft robots, i.e., to estimate the workspace of a given designed configuration of a soft robot. On the other hand, this chapter treats the inverse problem, i.e., to optimize the design of a soft robot in order to guarantee that certain points belong to its workspace. Solving such a problem is useful for robotic applications where we don't have an idea about the exact geometric parameters of the soft robot, but we do know what task is the soft robot trying to achieve, such as Pick-and-Place, and trajectory planning.

The actual process of a soft robot's design still follows intuition and the procedure of trial-and-error. Considering specific performance objectives, the classical design approach of soft robots is coupled with the uncertainty of achieving such objectives and also with the substantial economical expenditures necessary for the trial-and-error endeavors. For the purpose of accomplishing the desired objectives, it is therefore logical, for both economic and scientific reasons, to optimize the design of soft robots in a virtual environment before proceeding to its final physical conception. This thesis follows the above thought and seeks to answer the following question: Given a soft robot (illustrated by Fig. 6.1a), composed of a finite number of segments with bounded length L P L, and driven by mounted actuators (i.e., tendons) with bounded actuators length (i.e., tendons' length) α P A, and bounded actuators magnitude (i.e., tendons' tension) τ P T . How should we optimize the design of a soft robot in order to reach the desired position (e.g., blue-colored points u d in Fig. 6.1a). In other words, we seek to determine the optimal length parameter L of each segment, and the optimal actuators length parameter α such that the objective points (e.g., blue-colored points u d in Figs. [START_REF] Eugene | Introduction to numerical continuation methods[END_REF].1a) must be located inside the workspace of the optimally designed soft robot, as depicted by Fig. 6.1b.

In the present work, we assume that the segments radius r and the distance of the actuators to the mid-line d are both pre-determined and fixed parameters. However, the proposed method can be generalized to treat the case where r and d are varying (which will be explained in Chapter 7).

The reachable workspace of a soft robot has been defined in Definition 1 for the case that L and α are constant parameters. However, since we now seek to achieve the optimal values of L and α, the definition of the reachable workspace of a soft robot needs to be adapted, which should be a function of L and α. Definition 4. The workspace W E of the end-effector u of a soft robot modeled by Ψpτ, q, L, αq " 0 (with q P R n being the generalized strain vector), composed of bounded length of segments L P L, and driven via bounded actuators τ P T , with bounded actuators length α P A, is a subspace of R 3 , defined by the following: W E pL, αq " tu " ψpq, L, αq | Dpq P R n , τ P T q, s.t. Ψpτ, q, L, αq " 0u where ψ represents the geometric model of the soft robot which depends not only on q, but also on the parameters to be optimized L and α (more details will be discussed in 6.2.1), and Ψpτ, q, L, αq " KpL, αq q ´Hpq, L, αqτ ´F pq, L, αq.

To accomplish the vision depicted by Fig 6 .1b, we propose an optimization approach, which will be implemented on the PCS model and explained in the following.

Implementation of the Methodology

Influence of the Parameters to be Optimized on the PCS Static Model

Before we step into the details of the proposed optimization approach, we need to clarify the influence of the parameters to be optimized on the adopted mathematical model of soft robots.

In this work, we seek to optimize two parameters; the length L i of each segment i P 1 ¨¨¨N , and the actuators location, particularly the length parameter α k of each actuator k P 1 ¨¨¨n τ .

To be consistent with Definition 4, the workspace of the studied soft robot is now governed by: KpL, αq q " Hpq, L, αqτ `F pq, L, αq (6.1a) upq, L, αq " E 1 gpq, L, αq E T 2 (6.1b) with gpq, L, αq " e α 1 L 1 p ξ 1 e p1´α 1 qL 1 p ξ 1 ¨¨¨e αn τ L N p ξ N e p1´αn τ qL N p ξ N (6.2a)

p ξi " { Bq i qi ; i P 1 ¨¨¨N (6.2b)
The block-diagonal generalized stiffness matrix K, the generalized actuation matrix H, and the vector of generalized position-dependent external forces F are now established as:

KpL, αq " diag ˜BT q 1 « ż X 1 1 Σ1dX, ż X 2 1 Σ1dX ff Bq 1 , ¨¨¨, B T q N « ż X 1 N ΣN dX, ż X 2 N ΣN dX ff Bq N ¸, (6.3a)
H pq, L, αq "

» -˜BT q 1 « ż X 1 1 Hτ dX, ż X 2 1 Hτ dX ff¸T , ¨¨¨, ˜BT q N « ż X 1 N Hτ dX, ż X 2 N Hτ dX ff¸T fi fl T (6.3b) F pq, L, αq " N ÿ i"1 « ż X 1 i J T i Fe i dX, ż X 2 i J T i Fe i dX ff (6.3c)
where X 1 i " r0, α k L i s and X 2 i " rα k L i , p1 ´αk qL i s, for i P 1 ¨¨¨N , and k P 1 ¨¨¨n τ . The configuration matrix g i of a body i and the strain twist p ξ i of a body i can now be defined as follows:

g i pX, L, αq " # g s1 i " e α k Li p ξi , X P X 1 i g s2
i " e p1´α k qLi p ξi , X P X 2 i (6.4a)

p ξ i pX, L, αq " g ´1 i g 1 i (6.4b)
The geometric Jacobian is now formulated as follows:

J i pq, L, α, Xq " i ÿ l"0
Ad ´1 g plq ¨¨¨gi T g plq B q plq (6.5)

Remark 8. It is important to mention that, in addition to optimizing the segment length parameter L and the actuators length parameter α, it is possible to optimize the parameter representing the distance of the actuator to the mid-line d and also the segment radius parameter r by following the same procedure that will be explained for the parameters L and α.

Optimization-based Method

In this section, we present an optimization approach for the purpose of optimizing the design of the investigated slender soft robot in order to reach a single specific point and multiple target points, as depicted in Figs. [START_REF] Eugene | Introduction to numerical continuation methods[END_REF].2a and 6.2b.

To achieve one single point (as shown in Fig. 6.2a), the procedure consists of finding the optimal parameters L and α in order that the soft robot design is optimized enough to the level that its workspace will contain this point. However, when trying to reach multiple target points (as shown in Figs. 6.2b), the parameters must be optimized in order that the workspace of the designed soft robot contains all the desired points, which is clearly more difficult since the cost function will be more complex than the first case, because these objective points may be competing with each other.

Optimal Design to Achieve One Single Point

The proposed method consists of determining the optimal values of both the length parameter L of each segment denoted as L ˚, and the actuators length parameter α designated by α ˚, in order that, given these optimal values, the end-effector's workspace of the optimally designed soft robot will contain the position to be accessed (as illustrated in Fig. 6.2a). Denote u and u d as the end-effector's position and the desired position to be reached, according to the workspace Definition 4, we can then define the distance between u and u d as f pq, L, αq " 1 2 upq, L, αq ´ud

2 2
Hence, an optimal design of a soft robot in order to reach u d can be achieved through solving the following nonlinear constrained optimization problem:

rL ˚, α ˚, τ ˚s " arg min where τ ˚is the optimal actuators vector to reach the desired position u d . The vectors L ˚and α ˚represent the optimal length of segments and the optimal actuators placement, respectively (since they satisfy the bounded constraints and the nonlinear constraints) to minimize the nonlinear cost function f pq, L, αq.

To solve the above optimization problem, gradient-based method [START_REF] Petros | Robust adaptive control[END_REF] can be used which depends on the knowledge of the gradient of f pq, L, αq with respect to the vectors L, α, and τ , noted as ∇ rL,α,τ s f . The gradient method to solve (6.6) can be established by the following [START_REF] Petros | Robust adaptive control[END_REF]: Y s`1 " rL, α, τ s s ´λ∇ rL,α,τ s f prL, α, τ s s q s.t. rL, α, τ s s`1 " arg min rL,α,τ s Y s`1 ´rL, α, τ s τ P T , L P L, α P A KpL, αqq " Hpq, L, αqτ `F pq, L, αq upq, L, αq " E 1 gpq, L, αqE T 2 (6.7) with: rL, α, τ s 0 is known, and λ is an infinitesimal small step size.

Or in an alternative form; The gradient method to solve (6.6) can be established by the following [START_REF] Petros | Robust adaptive control[END_REF]:

Y s`1 " rL, α, τ s s ´λ∇ rL,α,τ s f prL, α, τ s s q rL, α, τ s s`1 " P rpY s`1 q, Project Y s`1 onto the feasible set (6.8) where P r is the projection operator.

Within the PCS framework, the following subsections are devoted to establishing the analytical form:

∇ rL,α,τ s f "

" ∇ L f, ∇ α f, ∇ τ f ı (6.9) 6.2.3.1 Calculation of ∇ τ f
First, the analytical form of the gradient of the cost function with respect to the actuator vector ∇ τ f has already been demonstrated and established in (3.7) and can be written in the following form:

∇ τ f " " Bf Bτ ı T " ∇ τ u rupq, L, αq ´ud s (6.10)
with ∇ τ u represents the gradient of the end-effector with respect to τ , and is defined according to Theorem 1 by: ∇ τ u " ∇ τ q ∇ u q " ∇ τ q " r0 3 Rpq, L, αqsJpq, L, αq ı T (6.11) where Jpq, L, αq is the geometric Jacobian defined in (6.5), the orientation matrix Rpq, L, αq is computed via the definition of the configuration matrix gpq, L, αq through: Rpq, L, αq " E 1 gpq, L, αq E T 1 , and ∇ τ q is the gradient of the generalized strain vector with respect to the actuators force vector and formulated according to Assumption 1 as: ∇ τ q " H T pq, L, αq " B ´1pτ, q, L, αq ‰ T (6.12) with Bpτ, q, L, αq " KpL, αq ´Hq pq, L, αqpI n b τ q ´Fq pq, L, αq

6.2.3.2 Calculation of ∇ L f
The analytical form of the gradient of the cost function with respect to the segments length parameter ∇ L f can be calculated using the exponential map properties, but first let us write the formula of ∇ L f using the principle of variable separation as:

∇ L f " " Bf BL ı T " " Bf Bu Bu BL ı T (6.13)
By definition, the partial derivative of the cost function with respect to the end-effector position vector is established by: Bf Bu " rupq, L, αq ´ud s T (6.14)

Next, the partial derivative of the end-effector position vector with respect to the segment length parameter can be decomposed using (6.1b) as follows: where the value of q corresponding to the value of τ can be calculated by solving (6.1a) using Proposition 1.

Bu BL " E 1 Bg BL E N `Bu Bq Bq BL with E N " I N b E T 2 ,
Next, we calculate Bg BL using the exponential map properties. First, we use (6.4a) in order to reformulate the exponential map gpq, L, αq from (6.2a) in the following form: gpq, L, αq " pg s1 1 g s2 1 q ¨¨¨pg s1 N g s2 N q (6.17) with g s1 i and g s2 i being defined in (6.4a). Using the above exponential map formula (6.17), Bg BL can be established as follows [START_REF] Kowal | A note on matrix differentiation[END_REF]: Last, with ∇ τ f from (6.10), and ∇ L f being calculated from (6.13), the remaining gradient to be computed in order to complete (6.9) is ∇ α f , which will be treated hereafter.

Bg BL " Bg s 1 1 BL rIN b pg s 2 1 ¨¨¨g s 1 N g s 2 N qs `gs 1 1 Bg s 2 1 BL rIN b pg s 1 2 ¨¨¨g s 1 N g s 2 N qs `¨¨g s 1 1 g s 2 1 ¨¨¨B g s 1 N BL rIN b g s 2 N s `gs 1 1 g s 2 1 ¨¨¨g s 1 N Bg s 2 N BL (
The analytical form of the gradient of the cost function with respect to the actuators length parameter ∇ α f can be calculated by following the same steps as those for ∇ L f : Finally, with ∇ τ f being calculated from (6.10), ∇ L f in (6.13), and ∇ α f using (6.19), we deduce the value of ∇ rL,α,τ s f in (6.9), which is necessary for solving (6.6) via (6.7).

∇ α f " " Bf Bα ı T " " Bf Bu Bu Bα ı T ( 6

Optimal Design to Achieve Multiple Points

In a general sense, it is important to optimize the segments length parameter L and the actuators length parameter α so that the end-effector of the studied soft robot can reach multiple target points, which can be useful for various applications such as pick-and-place and trajectory planning. Equivalently, doing this will result in optimizing the end-effector's workspace in order to contain the target points (as depicted by Fig. 6.2b).

Given a set of target points, noted as S d " tu p1q d , ¨¨¨, u pn S q d u P R 3ˆn S , for the obtained optimized soft robot with optimized structure, it will exist j different configurations q pjq with the corresponding end-effector position u pjq for j " 1 ¨¨¨n S , such that u pjq ´upjq d 2 2 " 0. Therefore, let us define the total distance of the end-effector's position with respect to all target points in S d as:

hpq, L, αq " 1 2 n S ÿ j"1 u pjq pq pjq , L, αq ´upjq d 2 2
then the optimal design of the investigated soft robot for the purpose of reaching all target points in S d can be obtained by solving the following nonlinear constrained optimization problem:

rL ˚,α ˚, T ˚s " arg min ı T P R nτ , for j P 1, ¨¨¨, n S , is the actuators vector associated with the desired position u dj . The vectors L ˚and α ˚represent the optimal length of each segment of the robot and the optimal actuators length, respectively.

Similar to the case of a single point, the gradient-based method [START_REF] Petros | Robust adaptive control[END_REF] can be used to solve the above optimization problem which needs the calculation of hpq, L, αq with respect to the vectors L, α, and τ , and formulated as:

∇ rL,α,T s h " " ∇ L h, ∇ α h, ∇ T h ı (6.24)
which will be detailed hereafter.

Calculation of ∇ T h

According to the principle of variable separation, we have:

∇ T h " " Bh BT ı T " " Bh Bu Bu BT ı T " " Bh Bu p1q
Bu p1q Bτ p1q , ¨¨¨, Bh Bu pn S q Bu pn S q Bτ pn S q ı T (6.25)

Then by definition, the partial derivative of h with respect to u pjq can be written as: Bh Bu pjq " " u pjq pq pjq , L, αq ´udj ı T (6.26)

Moreover, the partial derivative of the end-effector position u pjq with respect to its associated actuators vector τ pjq can be deduced from (6.11) Bq pjq Bτ pjq " rB pjq pτ pjq , q pjq , L, αqs ´1H pjq pq pjq , L, αq (6.28) where B pjq pτ pjq , q pjq , L, αq " KpL, αq ´Hq pjq pq pjq , L, αqpI n b τ pjq q ´Fq pjq pq pjq q By substituting Bq pjq Bτ pjq into (6.27), then (6.26) and (6.27) into (6.25), we finally obtain the analytical form of ∇ T h.

Calculation of ∇ T L

Similarly, using the principle of variable separation gives us the following formula:

∇ L h " " Bh BL ı T " " Bh Bu Bu BL ı T " j"n S ÿ j"1 " Bh Bu pjq Bu pjq BL T (6.29)
where the partial derivative of h with respect to u pjq has already been established in (6.26).

Therefore, we need only to find the analytical form of Bu pjq BL . For this, using (6.15), we can write the following:

Bu pjq BL " E 1 Bg pjq BL E N `r∇ u q pjq s T Bq pjq BL (6.30) with ∇ u q pjq " " r0 3 Rpq pjq , L, αqsJpq pjq , L, αq ı T being deduced from (6.27). Moreover, using (6.16), we can get the following:

Bq pjq BL " K ´1 ˆBH pjq BL τ pjq `BF pjq BL ´BK BL pI N b q pjq q ˙(6.31)
where the value of q pjq corresponding to the value of τ pjq can be calculated by solving (6.1a) using Proposition 1.

Next, we need to find the analytical form of Bg BL using the exponential map properties. For this, we reformulate the exponential map g pjq pq pjq , L, αq from (6.2a) as follows:

g pjq pq pjq , L, αq " ´gpjq Next, we have to find the analytical form of Bg Bα using the properties of the exponential map.

For this, using (6.32), Bg pjq Bα can be written as follows [START_REF] Kowal | A note on matrix differentiation[END_REF]: Finally, with ∇ T h being calculated from (6.25), ∇ L h in (6.29), and ∇ α h using (6.34), we deduce the value of ∇ rL,α,T s h in (6.24), which is necessary for solving (6.23).

Bg pjq Bα " Bg pjq 1,s 1 Bα " In τ b ´gpjq 1,s 2 ¨¨¨g pjq N,s 1 g pjq N,s 2 ¯ı `gpjq 1,s 1 Bg pjq 1,s 2 Bα " In τ b ´gpjq 2,s 1 ¨¨¨g pjq N,s 1 g pjq N,s 2 ¯ı `¨¨g pjq 1,s 1 g pjq 1,s 2 ¨¨¨B g pjq N,s 1 Bα " In τ b g pjq N,s 2 ı `gpjq

Validation and Analysis

In this section, we illustrate the presented results by implementing the proposed optimization approach on tendon-driven soft robots (as depicted by Figs. 6.1a). Table 3.1 presents the numerical values of the investigated soft robot.

In the examined scenarios, W E is estimated using the forward approach (See Section 2.5.2) and illustrated by the gray area. The workspace boundary BW E is estimated using the continuation approach (See Section 5.2) and depicted by red curves and green bifurcation points. The mounted tendons are depicted by black lines, and the desired points to reach are blue-colored. Also, the presented simulation scenarios were implemented using MATLAB (the function "fmincon" was used to solve the optimization problems formulated in this chapter) on an Intel Xeon(R) with a 16-GB RAM and a 3.50 GHz processor.

Scenario 1: 2-tendons driven planar soft robot

In the first scenario, we investigate a soft robot consisting of two segments and driven by two tendons installed on the Ozx plan (as shown in Fig. 6.3a).

In the initial state, the segments initial length vector is given L 0 " r0.2, 0.1s T m, and the initial actuators length vector is α 0 " r0.15, 0.15s T . We select the bounded tension magnitude T 1,2 " r0, 100sN allowing us to apply a couple magnitude within C T 1 " r0, 2.14sN ¨m, and C T 2 " r0, 1.42sN ¨m. The segments length bounds are L " r0, 0.3s ˆr0, 0.2sm, and the actuators length bounds are A " r0, 1s 2 .

The obtained results of W E and BW E for the initial state of the designed soft robot with the given initial values of the parameters L 0 and α 0 are depicted by Fig. 6.3b where we can clearly observe that the workspace does not contain the desired points. Firstly, we want to show that there are some cases where the optimal solution can only reach a desired point due to the length constraint. For this, we impose the segments length bounds in L " r0, 0.24s ˆr0, 0.16sm.

Using the proposed optimization approach, we apply (6.23) on the designed soft robot in order to achieve the optimal values of L ˚and α ˚in order to reach the desired points. By doing this, we find the optimal segments length parameter L ˚" r0.225, 0.15s T , and the optimal actuators length parameter α ˚" r0.89, 0.41s T , which gives the optimal design shown by Fig. 6.4a. As illustrated by Fig. 6.4b, we can observe that the optimal solution allows only the reachability of one single point due to the restrictions imposed on the segments length parameter L. To achieve the accessibility to multiple points, we can select a larger possibility for the length bounds of each segment as L " r0, 0.3s ˆr0, 0.2sm.

Using the proposed optimization approach, we apply (6.23) on the designed soft robot. By doing this, we find the optimal segments length parameter L ˚" r0.28, 0.19s T , and the optimal Figure 6.5: Soft robot's optimal design -achieving multiple points.

The computational time necessary to obtain the optimal results took approximately 0.445 seconds for this scenario.

Scenario 2: 3-tendons driven planar soft robot

In the second scenario, we investigate a soft robot comprised of three segments and controlled by three tendons installed on the Ozx plan (as shown in Fig. 6.6a).

In the initial state, the segments initial length vector is given L 0 " r0.15, 0.1, 0.05s T m, and the initial actuators length vector is α 0 " r0.15, 0.15, 0.15s T . We choose the bounded tension magnitude T 1,2,3 " r0, 50sN allowing us to apply a couple magnitude within C T 1 " r0, 1.07sN ¨m, and C T 2 " r0, 0.71sN ¨m, and C T 3 " r0, 0.35s. The segments length bounds are L " r0, 0.3s ˆr0, 0.2s ˆr0, 0.1sm, and the actuators length bounds are A " r0, 1s 3 .

The obtained results of W E and BW E for the initial state of the designed soft robot with the given initial values of the parameters L 0 and α 0 are depicted by Fig. 6.6b where we can clearly observe that the workspace does not contain the desired points. optimally designed soft robot now encompass the desired points as illustrated by Fig. 6.8b, allowing the end-effector to reach these target locations. The computational time necessary to obtain the optimal results took approximately 0.612 seconds for this scenario.

Conclusion

In this chapter, we have proposed an optimization approach based on the adopted PCS model for the purpose of optimizing the design of a soft robot with a slender shape. The approach was then validated on different tendon-driven soft robot configurations.

However, it is important to mention that we can extend the proposed approach to also optimize the radius of the segments and the distance of the actuator to the mid-line (which will be detailed in Chapter 7). Moreover, the optimization approach can be extended to treat different cases of actuation since the PCS method was already applied to different methods of actuation, e.g., fluidic actuation [START_REF] Renda | A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation[END_REF].

On the other hand, to extend the applicability of the optimization approach in order to treat the FEM case, we need to consider what parameters should be optimized in addition to adapting the optimal problem for the FEM case (which will be explained in Chapter 7).

Chapter 7

Conclusions and Perspectives

Conclusions

In this thesis, we have investigated the workspace estimation and design optimization of soft robots. On the one hand, the workspace estimation of soft robots provides crucial information about the space where they can operate, which is useful for control applications (e.g., Pick-and-Place, trajectory planning, etc). Moreover, evaluating the workspace also helps in determining the singularities and impediments of motion to the robot which can be used as valuable information for the design optimization of soft robots. On the other hand, the design optimization of soft robots is a useful strategy to avoid the procedure of trial-and-error when designing soft robots. Given specific performance objectives, it is more logical and interesting for both economic and scientific reasons to determine the optimal design of a soft robot in order to achieve those objectives in a virtual environment before designing the final prototype of the robot.

In Chapter 2, we have adopted two mathematical methods, the first is PCS which focuses on slender-shaped soft robots, and the second is FEM which treats soft robots with general form.

To answer the direct problem of this thesis, we have firstly presented a classic approach (forward approach), which consists of discretizing the actuators space, in order to estimate the workspace of soft robots. However, due to its high complexity, alternative approaches were proposed to efficiently estimate the workspace of soft robots.

Accordingly, we have proposed a first approach in Chapter 3 (optimization-based) that consists of estimating the workspace of soft robots by mapping the exterior boundaries of the workspace. This method successfully reduced the complexity of the workspace estimation compared to the forward approach.

To address the challenge of estimating the interior/exterior points in the workspace by discretizing the end-effector space, we have proposed a new approach based on interval analysis (forward-backward) in Chapter 4 that consists of exploring all the feasible configurations in the workspace of soft robots. Compared to the forward approach, this method successfully decreased the complexity of the workspace estimation.

In a final approach to the workspace estimation problem, we have proposed in Chapter 5 a continuation approach, based on the PCS model, that consists of mapping the interior and exterior boundaries of the workspace of soft robots. Notably, it provides information about the singularities of soft robots.

To treat the inverse problem of this thesis, we have proposed an optimization method based on the PCS model in Chapter 6 that consists of optimizing the design of soft robots in order to achieve specific desired points.

In the following, we present a brief comparison of the proposed workspace estimation approaches for both the PCS and FEM models.

PCS-based Workspace Estimation Approaches Comparison

Based on the PCS modeling method, three different approaches have been proposed to evaluate the workspace of soft robots and were validated through the forward approach. Each approach has its advantages and disadvantages, and the following summarizes their characteristics.

For the sake of clarity and simplicity, we consider one specific scenario (e.g., the planar soft robot composed of three segments and actuated by three tendons). The structure of the robot and its corresponding results are given in Fig. 7.1. Then, we summarize the advantages, disadvantages, operations complexity, and computational time of each investigated workspace approach in Table 7.1. From the first point of view, Taylor expansions can indeed be used to treat all cases. In fact, the 1st order Taylor expansion has already been used to calculate tangent vectors to regular boundary points (see Theorem 4), and the 2nd order Taylor expansion is also used in Theorem 5 to calculate tangent vectors to simple bifurcation points. However, to treat multiple bifurcation points, this methodology requires the calculation of partial derivatives of G up to the n b th-order, with n b being the number of branches passing through the multiple bifurcation point must be previously known. Thus, the main difficulty of this technique lies first in identifying n b , and secondly in calculating partial derivatives of G up to the n b th-order.

In an alternative way, the eigenvalue problem to calculate tangent vectors of multiple bifurcation points (see Section 5.2.3.3) can be extended to treat the simple bifurcation as well as regular boundary point cases. However, the extension of this approach cannot be directly implemented since it requires the existence and the knowledge of two orthonormal solutions γ 1 and γ 2 of Ψ Z . However, this is difficult to obtain because for both simple bifurcation points and regular boundary points becaause the row-rank deficiency of Ψ Z is 1, i.e., we can only have one orthonormal solution γ 1 that verifies the property Ψ Z γ 1 " 0. Therefore, the main challenge using this proposed technique would be to determine the second solution γ 2 such that Ψ Z γ 2 " 0 for multiple bifurcation points as well as regular boundary points.

Continuation Approach: an Extension to treat FEM

As shown in Section 5.2, the continuation approach requires the analytical representation of the parameters describing the static/kinematic model of soft robots. This is mainly due to the fact that tangent vector calculation in Section (5.2.3) requires the computation of the Jacobians Ψ Z and G Y , and also the Hessian matrix G Y Y . With the PCS model, this is possible because we have the analytical form of these mathematical parameters.

However, the FEM static model uses numerical approximations to establish its parameters, therefore, in the FEM framework we don't have the analytical form of those parameters. Hence, extending the applicability of the continuation approach to treat the FEM model will not be easy, and we think that the only possible way to perform this would be to calculate the Jacobians Ψ Z and G Y , and also the Hessian matrix G Y Y , via numerical techniques.

Design Optimization: Including more Parameters in PCS

In Chapter 6, Section 6.2, given specific performance objectives (e.g., achieving specific desired points), we have explained how to optimize the segments length parameter L and the actuators length parameter α, in order to achieve those objectives. What has been established for those parameters can also be extended in order to optimize the segments radius parameter r of the robot and the actuators distance to the mid-line d, as illustrated by Fig. with X " rL, α, r, ds. One additional step that will be required for the implementation of the methodology is to calculate the gradient of the objective function with respect to r and d, i.e., ∇ r h and ∇ d h.

Finally, solving (7.1) will yield the optimal segments length parameter L ˚, the optimal actuators length parameter α ˚, the optimal segments radius parameter r ˚, and the optimal actuators distance to the mid-line d ˚, allowing us to achieve a full design optimization of the investigated soft robot based on the PCS model.

Design Optimization: From PCS to FEM

The optimization approach proposed to achieve the optimal design of soft robots modeled by the PCS (see Chapter 6, section 6.2) can be extended to treat the FEM case as well.

For this, the same procedure followed for the PCS model might be applied to the FEM framework with necessary adaptation. The first step is to establish the static/kinematic model

Résumé Substantiel Introduction

Les robots rigides présentent de nombreux inconvénients lorsqu'ils fonctionnent dans des environnements dynamiques et fragiles et par conséquent, les robots déformables ont été un engin émergent qui a été graduellement étudié par les chercheurs afin de surmonter ces limitations et de faire face à de nouvelles applications robotiques. Les robots déformables sont fabriqués à partir de matériaux souples et flexibles, ce qui leur permet d'avoir plusieurs caractéristiques telles qu'une grande dextérité, des collisions prudentes et sans danger ainsi qu'une flexibilité importante. Toutes ces fonctionnalités offrent de nombreux avantages pour différentes applications, notamment l'exploration de l'environnement [START_REF] Aracri | Soft robots for ocean exploration and offshore operations: A perspective[END_REF][START_REF] Della | Soft robots. Encyclopedia of Robotics[END_REF][START_REF] Giorgio-Serchi | Hybrid parameter identification of a multi-modal underwater soft robot[END_REF][START_REF] Mintchev | A soft robot for random exploration of terrestrial environments[END_REF][START_REF] Shan | Untethered soft robots for future planetary explorations? Advanced Intelligent Systems[END_REF][START_REF] Zhu | Tuna robotics: A high-frequency experimental platform exploring the performance space of swimming fishes[END_REF] et les opérations médicales [START_REF] Burgner | A telerobotic system for transnasal surgery[END_REF][START_REF] Burgner-Kahrs | Continuum robots for medical applications: A survey[END_REF][START_REF] Cianchetti | Biomedical applications of soft robotics[END_REF][START_REF] Cianchetti | Soft robotics technologies to address shortcomings in today's minimally invasive surgery: the stiff-flop approach[END_REF][START_REF] Haouas | 4-dof spherical parallel wrist with embedded grasping capability for minimally invasive surgery[END_REF][START_REF] Mattheis | Flex robotic system in transoral robotic surgery: the first 40 patients[END_REF][START_REF] Mitros | Design and modelling of a continuum robot for distal lung sampling in mechanically ventilated patients in critical care[END_REF][START_REF] Simaan | A dexterous system for laryngeal surgery[END_REF][START_REF] Tsz | Soft robotics in medical applications[END_REF].

Cependant, en raison de leur conformité naturelle, la modélisation des robots souples est plus complexe que celle des robots rigides, car les robots déformables comportent un nombre élevé de degrés de liberté, leur déformation est non linéaire et ils sont caractérisés par des lois mécaniques différentes de celles des robots rigides. Par conséquent, des problèmes scientifiques tels que la détermination de l'espace de travail et l'optimisation de la conception des robots déformables émergent et avec eux les possibilités de nouvelles contributions dans le domaine de la robotique souple.

Même si les robots mous ont fait des progrès intéressants au cours des dernières années, le processus de conception d'un robot mou s'inspire toujours principalement d'une série de systèmes biologiques [START_REF] Rus | Design, fabrication and control of soft robots[END_REF], notamment la structure et le comportement d'espèces animales telles que la trompe des éléphants [START_REF] Michael | Kinematics and the implementation of an elephant's trunk manipulator and other continuum style robots[END_REF], et les bras d'une pieuvre [START_REF] Laschi | Soft robot arm inspired by the octopus[END_REF]. Une telle procédure de conception peut être utile pour les tests initiaux et l'expérimentation du robot souple conçu. Cependant, lorsqu'il est confronté à des objectifs de performance spécifiques tels que la planification et le contrôle de la trajectoire [START_REF] Delettre | Robust control of a planar manipulator for flexible and contactless handling[END_REF][START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF], il est possible que la conception d'un tel robot souple ne soit pas en mesure d'atteindre ses objectifs, notamment parce que sa plage d'accessibilité (c'est-à-dire son espace de travail) peut être restreinte.

Par conséquent, il est utile pour la communauté de la robotique douce de proposer des méthodologies génériques afin d'évaluer l'espace de travail des robots souples, et également d'aider et de guider la conception systématique de robots souples dans le but d'optimiser des objectifs de performance spécifiques.

L'estimation de l'espace de travail en robotique douce reste un sujet ouvert, et son importance est due aux nombreux avantages qu'elle peut apporter pour résoudre différents défis scientifiques de la robotique douce, principalement liés à la conception mécanique du robot et à la synthèse du contrôleur. D'une part, le résultat de l'évaluation de l'espace de travail peut fournir des informations sur l'accessibilité de l'objet à contrôler [START_REF] Delettre | Robust control of a planar manipulator for flexible and contactless handling[END_REF][START_REF] Della | Dynamic control of soft robots interacting with the environment[END_REF][START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF], c'est-à-dire en identifiant si la position de l'objet appartient à l'espace de travail du robot souple, épargnant ainsi au contrôleur le temps de conception qui peut être perdu à essayer de contrôler un objet inaccessible en dehors de l'espace de travail du robot. D'autre part, les informations sur l'espace de travail sont également cruciales lors de la planification d'une trajectoire réalisable à suivre par des robots mous, car les positions de départ et d'arrivée d'une trajectoire réalisable doivent appartenir à l'espace de travail d'un robot [START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF]. En outre, l'analyse de l'espace de travail est également une étape nécessaire pour contribuer à une conception efficace, ciblée et optimale des robots mous [START_REF] Hiller | Automatic design and manufacture of soft robots[END_REF], en améliorant leur stabilité et en augmentant la portée de leur manipulabilité et de leur accessibilité.

Pour atteindre cet objectif, la thèse actuelle adopte deux modèles mathématiques différents, la méthode de déformation constante par morceaux (PCS) et la méthode des éléments finis (FEM) pour décrire la déformation des robots mous de forme élancée et de forme générale, respectivement. Sur la base de ces modèles mathématiques adoptés, plusieurs approches sont proposées pour estimer l'espace de travail des robots mous. Enfin, une méthode efficace d'optimisation de la conception des robots mous en vue d'atteindre des objectifs de performance spécifiques est présentée.

Modélisation

Dans ce chapitre, nous avons présenté les modèles mathématiques des méthodes PCS et FEM dans le but de modéliser les robots mous ayant une forme élancée et une forme générale, respectivement.

Ensuite, nous avons établi la définition de l'espace de travail d'un robot souple. Enfin, nous avons proposé une approche directe pour estimer l'espace de travail des robots mous. Cette approche a ensuite été validée en utilisant un robot souple de type tronc [START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF] pour les cas PCS et FEM.

Approche d'optimisation pour l'estimation de l'espace de travail Approche d'analyse intervalle pour l'estimation de l'espace de travail En raison des limitations rencontrées par l'approche d'optimisation qui ne peut pas fournir de connaissances sur les configurations intérieures de l'espace de travail, il est donc nécessaire de fournir une méthode qui peut surmonter cette limitation, mais qui est également basée sur la deuxième stratégie, c'est-à-dire la discrétisation de l'espace de l'effecteur.

Une telle méthode [START_REF] Amehri | Fem based workspace estimation for soft robots: a forward-backward interval analysis approach[END_REF] est basée sur des techniques d'analyse d'intervalles [START_REF] Gouttefarde | Interval-analysis-based determination of the wrench-feasible workspace of parallel cable-driven robots[END_REF][START_REF] Jaulin | Interval analysis[END_REF][START_REF] Merlet | Determination of 6d workspaces of gough-type parallel manipulator and comparison between different geometries[END_REF][START_REF] Merlet | Parallel robots[END_REF][START_REF] Merlet | On the workspace of suspended cable-driven parallel robots[END_REF], et consiste à discrétiser l'espace de l'effecteur en partant d'une configuration initiale atteignable puis en explorant l'ensemble de l'espace atteignable possible pour finalement estimer l'espace de travail d'un robot souple.

L'approche proposée a été appliquée avec succès aux modèles PCS et FEM, où nous avons montré son efficacité à réduire la complexité et le temps de calcul nécessaires pour estimer l'espace de travail des robots mous, contrairement à l'approche directe qui explose exponentiellement lorsque la dimension des actionneurs augmente. D'autre part, comme cette approche consiste à déterminer toutes les configurations possibles qui sont réalisables, elle est utile pour identifier les informations intérieures et extérieures de l'espace de travail. Cependant, cette approche est exhaustive dans le sens où elle explore toutes les configurations réalisables de l'espace de travail au lieu de cartographier uniquement ses limites intérieures et extérieures.

Approche de continuité pour l'estimation de l'espace de travail

En raison des limites rencontrées par l'approche d'optimisation qui ne peut pas fournir la connaissance des configurations intérieures et des limites de l'espace de travail des robots mous, et également en raison des limites de l'approche d'analyse d'intervalle qui est exhaustive en explorant l'espace de travail entier, il est logique de fournir une méthodologie alternative nous permettant de déterminer seulement mais toutes les limites des espaces de travail des robots mous, à la fois intérieurs et extérieurs.

Une telle méthode est appelée l'approche de continuation [START_REF] Amehri | Workspace boundary estimation for soft manipulators using a continuation approach[END_REF] qui est basée sur la théorie de la bifurcation. Cette méthode est fondée sur la théorie mathématique de la bifurcation qui a été introduite à l'origine par le mathématicien français Henri Poincaré [START_REF] Poincaré | Sur l'équilibre d'une masse fluide animée d'un mouvement de rotation[END_REF], et elle consiste à étudier les solutions d'équations différentielles non linéaires, en plus de l'étude des variations possibles de la structure topologique d'une famille de champs de vecteurs. Par la suite, Henri Poincaré a classé différents types de points de bifurcation [START_REF] Poincaré | New methods of celestial mechanics[END_REF]. En conséquence, cette approche prend en compte le comportement possible de bifurcation dans le calcul de continuation tout en cartographiant les limites intérieures et extérieures des espaces de travail des robots logiciels.

L'approche de continuation permet de cartographier les frontières intérieures et extérieures des robots mous aux formes élancées en utilisant le modèle PCS. Cependant, il est important de préciser que cette approche ne peut pas fournir d'informations sur tous les points de l'espace de travail, mais seulement sur les points limites.

Otpimisation de la conception des robots mous

Le processus actuel de conception d'un robot souple fait toujours appel à l'intuition et à la procédure d'essai-erreur. Si l'on considère des objectifs de performance spécifiques, l'approche

Estimation de l'Espace de Travail et l'Optimisation de la Conception des Robots Déformables Résumé

Les robots rigides présentent de nombreux inconvénients lorsqu'ils fonctionnent dans des environnements dynamiques et fragiles et par conséquent, les robots déformables ont été un engin émergent qui a été graduellement étudié par les chercheurs afin de surmonter ces limitations et de faire face à de nouvelles applications robotiques. Les robots déformables sont fabriqués à partir de matériaux souples et flexibles, ce qui leur permet d'avoir plusieurs caractéristiques telles qu'une grande dextérité, des collisions prudentes et sans danger ainsi qu'une flexibilité importante. Toutes ces fonctionnalités offrent de nombreux avantages pour différentes applications, notamment l'exploration de l'environnement et les opérations médicales. Cependant, en raison de leur conformité naturelle, la modélisation des robots souples est plus complexe que celle des robots rigides, car les robots déformables comportent un nombre élevé de degrés de liberté, leur déformation est non linéaire et ils sont caractérisés par des lois mécaniques différentes de celles des robots rigides. Par conséquent, des problèmes scientifiques tels que la détermination de l'espace de travail et l'optimisation de la conception des robots déformables émergent et avec eux les possibilités de nouvelles contributions dans le domaine de la robotique souple. L'évaluation de l'espace de travail offre de nombreux avantages pour différentes applications de la robotique souple principalement liées à leur conception et à leur contrôle. En conséquence, cette thèse étudie l'estimation de l'espace de travail et l'optimisation de la conception des robots souples. Afin d'accomplir cette tâche, on propose deux méthodes différentes pour la modélisation des robots déformables, la première est la méthode de déformation constante par morceaux (PCS), qui est utilisée pour la modélisation des robots déformables avec une géométrie continue et la deuxième est la méthode des éléments finis (FEM), qui est utilisé pour la modélisation des robots déformables avec une géométrie générale. Ensuite, basée sur ces modèles mathématiques, cette thèse propose différentes méthodologies pour estimer l'espace de travail des robots souples. Deux stratégies ont été proposées pour l'estimation de l'espace de travail, la première consiste à discrétiser l'espace des entrées (actionneurs), et la deuxième consiste à discrétiser l'espace des sorties (l'effecteur du robot). Cependant, la première stratégie est inefficace, car elle dépend de la dimension des actionneurs et celle-ci varie en fonction de la configuration du robot déformable étudiée. En revanche, la seconde stratégie présente une méthodologie stable et efficace pour l'estimation de l'espace de travail puisque l'espace de l'effecteur est toujours constant (inférieur ou égal à 3, si nous nous concentrons sur l'aspect position de l'espace de travail), quelles que soient la configuration du robot souple étudié et la dimension des actionneurs. Les approches proposées pour l'estimation de l'espace de travail ont ensuite été appliquées aux deux modèles mathématiques adoptés et validées à l'aide de différentes configurations de robots déformables. Enfin, cette thèse propose une approche d'optimisation basée sur les modèles mathématiques adoptés pour optimiser la conception des robots souples afin d'atteindre certains objectifs spécifiques.

Mots-clefs :

Robot Mou, Modèle Cinématique/Dynamique, Espace de Travail, Optimisation.

Workspace Estimation and Design Optimization of Soft Robots Abstract

Soft robots are an emergent instrument that has gradually been investigated by researchers in the recent years to overcome limitations of traditional rigid robots as well as to propose novel robotic applications. Rigid robots are challenged when operating in restricted and dynamic environments. Being made from a soft and flexible material, soft robots provide many benefits such as high dexterity, safe interactions, and increased adaptability, which are very useful for various applications, especially the manipulation of fragile objects, environment exploration, and medical operations. However, owing to their natural conformity and compliance, soft robots consist of a large number of degrees of freedom, their deformation is highly nonlinear, and they are characterized by different mechanical laws compared to that of rigid robots. All these aspects make their modeling more complex. Consequently, scientific challenges such as workspace evaluation and design optimization of soft robots arise and with them the opportunities of new contributions in the field of soft robotics. The workspace evaluation provides many benefits for different soft robotic applications mainly related to their design and control. Accordingly, the present thesis investigates the workspace evaluation and design optimization of soft robots. For this, two different methods are adopted for the mathematical modeling of soft robots. The first is the Piece-wise Constant Strain (PCS), which focuses on slender-shaped soft robots, and the second is the Finite Element Method (FEM), which concerns soft robots with a general shape. Based on the adopted mathematical models, this thesis proposes different methodologies to estimate the workspace of soft robots. In fact, two strategies can be followed for the workspace estimation of soft robots: the first consists of discretizing the inputs (actuators) space, and the second consists of discretizing the outputs (end-effector) space. However, the first strategy is not efficient as it depends on the dimension of the actuators which are used to control the investigated soft robot, and varies corresponding to the particular studied configuration. Conversely, the second strategy presents a stable and efficient way for the workspace estimation since the end-effector's space is always constant (smaller or equal to 3, if we focus on the position access) regardless of the studied soft robot's configuration and the dimension of actuators. Each proposed workspace estimation approach was then applied to both the PCS and the FEM models and validated via a variety of soft robots' configurations. Finally, this thesis proposes a model-based optimization approach in order to optimize the design of soft robots for the purpose of achieving specific performance objectives. Keywords: Soft Robots, Kinematic/Dynamic Modeling, Workspace Estimation, Optimization.

Figure 1 . 1 :

 11 Figure 1.1: Classification of robots [44, 52, 121].

Figure 1 . 2 :

 12 Figure 1.2: Rigid Robots Evolution: A-Mechanical Knight. B-Digesting Duck. C-Teleautomaton. D-Elektro. E-Elsie and Elmer[START_REF] Sabbatini | Imitation of life: A history of the first robots[END_REF]. F-Unimate[START_REF] Waurzyniak | Masters of manufacturing: Richard morley[END_REF]. G-WABOT-1[START_REF] Kuo | Robotics and Mechatronics: Proceedings of the 6th IFToMM International Symposium on Robotics and Mechatronics (ISRM 2019)[END_REF]. H-Famulus[START_REF] Singh | Evolution of industrial robots and their applications[END_REF]. I-PUMA. J-Delta and ABB Flex-Pitcher robots[START_REF] Douglas R Ewing | Robots in the operating room-the history[END_REF][START_REF] Gasparetto | A brief history of industrial robotics in the 20th century[END_REF][START_REF] Hägele | Industrial robotics[END_REF].

Figure 1 . 11 :

 111 Figure 1.11: Hybrid approaches [94].

Figure 2 . 1 :

 21 Figure 2.1: Schematic of the PCS kinematics.

Figure 2 . 2 :

 22 Figure 2.2: Schematic of the tendon-driven actuation of one body.

Figure 2 . 3 :

 23 Figure 2.3: Linear tetrahedron element. Take the trunk-like soft robot as an example, Fig. 2.4 shows the corresponding discretized version in the domain.

Figure 2 . 4 :

 24 Figure 2.4: Discretized mesh model for trunk-like soft robot using the linear tetrahedron element.

c 11 c 12 c 13 11 "

 1311 c 22 " c 33 " E p1 ´νq p1 `νq p1 ´2νq c 12 " c 13 " c 23 " E ν p1 `νq p1 ´2νq c 44 " c 55 " c 66 " E 2 p1 `νq

Figure 2 . 6 :

 26 Figure 2.6: Trunk-like soft robot [129].

  (a) W E : Real (gray) and PCS-simulation (red). (b) W E : Real (gray) and FEM-simulation (red).

Figure 2 . 7 :

 27 Figure 2.7: Forward approach -W E obtained from the real experimentation (gray) and W E achieved by the simulations of the adopted mathematical models (red).

s| « 7 . 3 | « 4 . 8

 7348 ) and FEM (noted V F EM s ) simulation results. The volume calculated from the set of points obtained by the real robot is V r « 0.0021m 3 , and the volumes obtained the PCS and FEM simulation results are V P CS s « 0.0028m 3 and V F EM s « 0.0025m 3 , respectively. Then, by comparing both results, we find |V r ´V P CS s ˆ10 ´4m 3 for the accuracy of the PCS result and |V r ´V F EM s ˆ10 ´4m 3

Fig 3 .

 3 1).

Figure 3 . 1 :

 31 Figure 3.1: For a radiating vector v determine its associated boundary point u b .

3. 2 . 2

 22 Finding a point on BW EIn order to map the surface enveloping the workspace W E , we need first to figure out how to determine a boundary point u b . For concrete applications, it is a nontrivial task to find a point on the workspace boundary BW E (as shown in Fig.3.2). To clearly visualize this problem, we present the following simplified illustration (see Fig.3.2).

Figure 3 . 2 :

 32 Figure 3.2: For a radiating vector v, calculate its corresponding actuators τ b via (3.6), and deduce its associated boundary point u b via (2.48).

Figure 3 . 3 :

 33 Figure 3.3: Numerical map of BW E -for each radiating vector v pjq , compute its corresponding actuators τ bj via (3.8), and deduce its associated boundary point u bj via (2.48).

Figure 3 . 4 :

 34 Figure 3.4: Non-convexity problem of BW E -sudden leap between the deduced boundary points u bn 1 and u bn 2 .

  3.5), we present a simplified figure describing two iterations of the proposed strategy for solving the non-convexity problem.

Figure 3 . 5 :

 35 Figure 3.5: Non-convexity solution -for each point v l , calculate its related actuators via (3.6), and deduce its associated boundary point via (2.48).

  (a) Soft robot structure. (b) Forward and Optimization approach.

Figure 3 . 6 :

 36 Figure 3.6: Scenario 1 -workspace (gray) W E , boundary points (green), and workspace boundary (red) BW E .

  (a) Soft robot structure. (b) Forward and Optimization approach.

Figure 3 . 7 :

 37 Figure 3.7: Scenario 2 -Workspace (gray) W E , boundary points (green), and workspace boundary (red) BW E .

  (a) Soft robot structure. (b) Forward and Optimization approach.

Figure 3 . 8 :

 38 Figure 3.8: Scenario 3 -Workspace (gray) W E , boundary points (green), and workspace boundary (red) BW E .

  (a) Soft robot structure. (b) Forward and Optimization approach.

Figure 3 . 9 :

 39 Figure 3.9: Scenario 4 -Workspace (gray) W E , boundary points (green), and workspace boundary (red) BW E .

  (a) Soft robot structure. (b) Oxz View. (c) Oyx View. (d) Oyz View.

Figure 3 . 10 :

 310 Figure 3.10: Scenario 5 -Workspace (gray) W E , boundary points (green), and workspace boundary (red) BW E .

Figure 3 . 11 :

 311 Figure 3.11: Operations complexity and computational time of the forward and optimization approach for the workspace estimation of soft robots using the PCS model.

  (a) Soft Robot Structure nτ " 2. (b) Oxy plan view. (c) Ozy plan view. (d) Ozx plan view.

Figure 3 . 12 :

 312 Figure3.12: Scenario 1 -W E estimation via forward approach (blue points) and BW E estimated via optimization approach (green points).

  (a) Soft Robot Structure nτ " 3. (b) Oxy plan view. (c) Ozy plan view. (d) Ozx plan view.

Figure 3 . 13 :

 313 Figure 3.13: Scenario 2 -W E estimation via forward approach (blue points) and BW E estimated via optimization approach (green points).

  (a) Soft Robot Structure nτ " 4. (b) Oxy plan view. (c) Ozy plan view. (d) Ozx plan view.

Figure 3 . 14 :

 314 Figure 3.14: Scenario 3 -W E estimation via forward approach (blue points) and BW E estimated via optimization approach (green points).

Figure 3 . 15 :

 315 Figure 3.15: Operations complexity and computational time of the forward and optimization approach for the workspace estimation of soft robots using the FEM model.

  p0, 2q T belongs to AX , whereas it does not belong to the actual value set B " tAx | x P X u, as shown by Fig. 4.1.

Figure 4 . 1 :

 41 Figure 4.1: Over-estimation introduced by the wrapping effect.

EĂ R 3

 3 with a prescribed precision and obtain, as shown by Step 2 in Figs. 4.2 and 4.3, the associated discretized space W pjq E d (subscript "d" stands for discretization). Next, based on the discretized space W pjq E d

  small neighborhoods (called the η-neighborhood in Step 2) around the current configuration u pj´1q (which is illustrated by Fig. 4.3, and will be explained in detail in Step 2, Section 4.3.

  Step 5 and depicted by Fig.4.7) is verified.

Figure 4 . 2 :

 42 Figure 4.2: Forward-Backward approach.

pjqs

  are determined by: δq pjq " A `δτ pj´1q `A´δ τ pj´1q `K´1 F pq pj´1q q and δq pjq " A `δτ pj´1q

  .

  neighborhood space W pjq Es are illustrated by Fig. 4.3.

Figure 4 . 3 :

 43 Figure 4.3: Feasible small neighborhood W pjq Es .

4

 4

  .4.

Figure 4 . 4 :

 44 Figure 4.4: Backward estimation of τ pjq .

Figure 4 . 5 :

 45 Figure 4.5: New actuators bound T pjq .

E

  calculated via Step 1, its related discretized space W pjq E d computed via Step 2, and its associated feasible small neighborhood W pjq Es deduced from (4.4

Figure 4 . 6 :

 46 Figure 4.6: Discard pre-explored feasible neighborhoods.

Figure 4 . 7 :

 47 Figure 4.7: Stop condition.

E

  as stated in Step. 1 of Fig. 4.2, based on which a feasible small neighborhood W p1q Es is determined as mentioned in Step. 2 of Fig. 4.2. Next, as indicated in Step. 3 of Fig. 4.2, we backwardly compute the feasible actuators vector τ p1q of each configuration in the feasible small neighborhood W p1q Es . Finally, we find the new actuators bound T p1q as demonstrated in Step. 4 of Fig. 4.2. The same process is pursued in the next iterations until we meet the stop condition as stated in Step. 5 and shown by Fig. 4.7.

Figure 4 . 8 :

 48 Figure 4.8: Scenario 1 -W E estimation via forward approach (black points) and forward-backward approach (gray area).

  (a) Soft robot structure. (b) Oxz plan view of W E .

Figure 4 . 9 :

 49 Figure 4.9: Scenario 2 -W E estimation via forward approach (black points) and forward-backward approach (gray area).

  (a) Soft robot structure. (b) Oxz plan view of W E .

Figure 4 . 10 :

 410 Figure 4.10: Scenario 3 -W E estimation via forward approach (black points) and forwardbackward approach (gray area).

  (a) Soft robot structure. (b) Oxz plan view of W E .

Figure 4 . 11 :

 411 Figure 4.11: Scenario 4 -W E estimation via forward approach (black points) and forwardbackward approach (gray area).

  (a) Soft robot structure. (b) Oxz plan view of W E . (c) Oyx plan view of W E . (d) Oyz plan view of W E .

Figure 4

 4 Figure 4.12: Scenario 5 -W E estimation via forward approach (black points) and forwardbackward approach (gray area).

  13a and 4.13b, we summarize the operations' complexity and computational time of each scenario.

Figure 4 . 13 :

 413 Figure 4.13: Operations complexity and computational time of the forward and forward-backward approaches for the workspace estimation of the investigated scenarios using the PCS model.

  (a) Soft Robot Structure nτ " 2. (b) Oxy plan view. (c) Ozy plan view. (d) Ozx plan view.

Figure 4 . 14 :

 414 Figure 4.14: Scenario 1 -W E estimation via forward approach (blue points) and forwardbackward (gray area) approach.

4

 4

  .6).

  (a) Soft Robot Structure nτ " 3. (b) Oxy plan view. (c) Ozy plan view. (d) Ozx plan view.

Figure 4 . 15 :

 415 Figure 4.15: Scenario 2 -W E estimation via forward approach (blue points) and forwardbackward (gray area) approach.

4

 4

  .6).

  (a) Soft Robot Structure nτ " 4. (b) Oxy plan view. (c) Ozy plan view. (d) Ozx plan view.

Figure 4 . 16 :

 416 Figure 4.16: Scenario 3 -W E estimation via forward approach (blue points) and forwardbackward (gray area) approach.

Figure 4

 4 Figure 4.17: Operations complexity and computational time of the forward and forward-backward approaches for the workspace estimation of the investigated scenarios using the FEM model.

Figure 5 . 1 :

 51 Figure 5.1: Overall view of the continuation method.

Figure 5 . 2 :

 52 Figure 5.2: A boundary point u b0 .

Figure 5 . 3 :

 53 Figure 5.3: Regular boundary point -tangent vector.

Figure 5 . 5 :

 55 Figure 5.5: Multiple bifurcation point -tangent vector.

  p (as depicted by Fig. 5.6) of the appropriate solution is calculated as follows: Y p " Y b0 ˘εhpY b0 q (5.41) Next, based on the predicted state Y p , its associated solution Y b (as shown in Fig. 5.6) must verify the nonlinear constraint condition (5.7) of the workspace boundary solution set.

Figure 5 . 6 :

 56 Figure 5.6: Mapping one-dimensional solution curves.

Figure 5 . 10 :

 510 Figure 5.10: Scenario 4 -W E estimation via forward approach (gray points) and workspace boundaries BW E computed via the continuation approach (red curves and green bifurcation points).

  (a) Soft robot structure. (b) Oxz plan view of W E and BW E . (c) Oyx plan view of W E and BW E . (d) Oyz plan view of W E and BW E .

Figure 5 . 11 :

 511 Figure 5.11: Scenario 5 -W E estimation via forward approach (gray points) and workspace boundaries BW E computed via the continuation approach (red curves and green bifurcation points).

  (a) Soft robot structure. (b) Oxz plan view of W E and BW E . (c) Oyx plan view of W E and BW E . (d) Oyz plan view of W E and BW E .

Figure 5 . 12 :

 512 Figure 5.12: Scenario 6 -W E estimation via forward approach (gray points) and workspace boundaries BW E computed via the continuation approach (red curves and green bifurcation points).

  (a) Soft robot structure. (b) Oxz plan view of W E and BW E . (c) Oyx plan view of W E and BW E . (d) Oyz plan view of W E and BW E .

Figure 5 . 13 :

 513 Figure 5.13: Scenario 7 -W E estimation via forward approach (gray points) and workspace boundaries BW E computed via the continuation approach (red curves and green bifurcation points).

  (a) Soft robot structure. (b) Oxz plan view of W E and BW E . (c) Oyx plan view of W E and BW E . (d) Oyz plan view of W E and BW E .

Figure 5 . 14 :

 514 Figure 5.14: Scenario 8 -W E estimation via forward approach (gray points) and workspace boundaries BW E computed via the continuation approach (red curves and green bifurcation points).

Figure 5 . 15 :

 515 Figure 5.15: Operations complexity and computational time of the forward and continuation approaches for the workspace and workspace boundaries estimation of the investigated scenarios using the PCS model.

Figure 5 . 16 :

 516 Figure 5.16: Multi-component workspaces (a); Holes/Voids (b).

  5.16 as an example. Consider, as depicted by Fig. 5.16b, a boundary point Y b E lying on the exterior boundary BW p1q

Figure 6 . 1 :

 61 Figure 6.1: Soft robot's design.

  (a) Optimal Design for a single point. (b) Optimal Design for multiple points.

Figure 6 . 2 :

 62 Figure 6.2: Soft robot's design to achieve a single point (a), and multiple points (b).

  (a) Initial design. (b) W E (gray) and BW E (red).

Figure 6 . 3 :

 63 Figure 6.3: Soft robot's initial design.

  (a) Optimal design. (b) W E (gray) and BW E (red).

Figure 6 . 4 :

 64 Figure 6.4: Soft robot's optimal design -achieving only one point.

  (a) Optimal design. (b) W E (gray) and BW E (red).

Figure 6 . 8 :

 68 Figure 6.8: Soft robot's optimal design -achieving multiple points.

  (a) Soft robot structure. (b) Optimization-based approach. (c) Interval Analysis approach. (d) Continuation Approach.

Figure 7 . 1 :

 71 Figure 7.1: Different workspace methodologies applied to the investigated soft robot modeled through the PCS model.

  7.3.

Figure 7 . 3 :u 2 2 2 ( 7

 732227 Figure 7.3: Design Optimization -Including more Parameters in PCS.

Table 1 -

 1 continued from previous page

	δq	Generalized strain vector with respect to prior configuration	Displacement vector with respect to prior configuration
	τ ‰ W E " " τ 1 ¨¨¨τ nτ P R nu Workspace of a soft robot
	BW E	P R nu´1 Soft robot workspace boundary
	W E d	Discretized space of the current configuration
	W Es	Feasible small displacement neighborhood space around the current configuration
	B	" rB, Bs Lower B and upper B bounds of a domain B
	R	P SOp3q Orientation matrix	-
	X	P r0, Ls Ă R Abscissa along the robot arm -
			Continued on next page

T Minimum actuators bound vector τ " rτ 1 ¨¨¨τ nτ s T Maximum actuators bound vector T " rτ 1 , τ 1 s ˆ¨¨¨ˆ"τ nτ , τ nτ ‰ Actuators bounds u P R nu End-effector position vector in the inertial frame.

Table 1 -

 1 continued from previous page

	Symbol Model	Piece-wise Constant Strain	Finite Element Method
	L	" : Maximum segments length bounds " ‰ T L 1 ¨¨¨L N	-
	L	" : Segments length bounds " L 1 , L 1 ‰ ˆ¨¨¨ˆ"L N , L N	‰	-

Table 1 :

 1 Acronyms and Glossary for the adopted mathematical modeling methods.

	Chapter 1
	Introduction
	1.1 General Introduction

1.1.1 Robotics: Definitions and Classes

  1. Walid, Amehri, Gang Zheng, Alexandre Kruszewski, and Federico Renda. "Discrete Cosserat Method for Soft Manipulators Workspace Estimation: An Optimization-Based Approach." Journal of Mechanisms and Robotics 14, no. 1 (2021): 011012. Summary: In this article, we presented an optimization-based approach in order to estimate the exterior workspace boundary of slender-shaped soft robots using the PCS modeling method. In this letter, we provided a continuation approach based on the theory of bifurcation in order to estimate the interior and exterior boundaries slender-shaped soft robots using the PCS modeling method. 3. Amehri, Walid, Gang Zheng, and Alexandre Kruszewski. "FEM-based Exterior Workspace Boundary Estimation for Soft Robots via Optimization." IEEE Robotics and Automation Letters (2022).

	2. Amehri, Walid, Gang Zheng, and Alexandre Kruszewski. "Workspace Boundary Estimation
	for Soft Manipulators Using a Continuation Approach." IEEE Robotics and Automation
	Letters 6, no. 4 (2021): 7169-7176.
	Summary:

Summary: In this paper, we extended the applicability of the optimization-based approach by applying it to the FEM modeling method in order to estimate the exterior workspace boundary of soft robots. 4. Amehri, Walid, Gang Zheng, and Alexandre Kruszewski. Soft Robotics SORO, FEM-based Reachable Workspace Estimation of Soft Robots using an Interval Analysis approach.

Table 3 .

 3

	Physical parameters	First	Second	Third
	L i (m)	0.3	0.2	0.1
	r i (m)	0.06	0.04	0.02
	J xi (m 4 )	2.03 ˆ10 ´5 4.02 ˆ10 ´6 2.51 ˆ10
				´7
	A i (m 2 )	0.0113	0.0050	0.0013
	E (Pa)		18 ˆ10 4
	ν		0.45
	G (Pa)		62.069 ˆ10 3

´7

J py,zq i (m 4 ) 1.07 ˆ10 ´5 2.01 ˆ10 ´6 1.25 ˆ10 1: Numerical values of the soft robot physical parameters.

Table 3 .

 3 2: Operations Complexity and computational time: forward vs optimization.

			Operations	Time (seconds)
	Inputs Approach Forward Optimization	Forward	Optimization
	n τ " 2	200 2	168	66.35	1.82
	n τ " 3	100 3	220	2962	3.52
	n τ " 4	100 4	168	« 1.65 ˆ10 5	3.92
	n τ " 6	100 6	220	« 1.65 ˆ10 9	13.58
	n τ " 9	100 9	8640	« 1.65 ˆ10 12	510

Table 3 .

 3 3: Operations and time computational complexity for the FEM case: forward vs optimization-based.

			Operations	Time (seconds)
	Inputs Approach Forward Optimization	Forward	Optimization
	n τ " 2	100 2	864	1463	984
	n τ " 3	100 3	864	« 1.4 ˆ10 5	1034
	n τ " 4	100 4	864	« 1.4 ˆ10 7	1039

  δq pjq ď δq

pjq with: δq pjq " A `δτ pj´1q `A´δ τ pj´1q `K´1 F pq pj´1q q and δq pjq " A `δτ pj´1q

Table 4 .

 4 1: Operations Complexity and computational time: forward vs forward-backward.

		Operations	Time (seconds)
	Inputs Approach Forward	Forward-Backward	Forward	Forward-Backward
	n τ " 2	200 2	1820	66.35	17
	n τ " 3	100 3	2275	2962	72
	n τ " 4	100 4	3640	« 1.65 ˆ10 5	187
	n τ " 6	100 6	4550	« 1.65 ˆ10 9	2436
	n τ " 9	100 9	13350	« 1.65 ˆ10 12	4679

(a) Operations. (b) Time

Table 4 .

 4 2: Operations and time computational complexity: forward vs forward-backward.

		Operations	Time (seconds)
	Inputs Approach Forward	Forward -Backward	Forward	Forward -Backward
	n τ " 2	100 2	635	1463	757
	n τ " 3	100 3	1941	« 1.4 ˆ10 5	2063
	n τ " 4	100 4	3667	« 1.4 ˆ10 7	3733

Table 5 .

 5 2: Operations Complexity and computational time: forward vs continuation.

			Operations	Time (seconds)
	Inputs Approach	Forward	Continuation	Forward	Continuation
	n τ " 2 (Sc1)	200 2	87	66.35	9.54
	n τ " 3 (Sc3)	100 3	327	2962.92	40.62
	n τ " 4 (Sc2, Sc5)	100 4	p687, 688q	« 1.65 ˆ10 5	p96, 107q
	n τ " 6 (Sc4, Sc6, Sc7)	100 6	p4126, 4107, 4122q	« 1.65 ˆ10 9 p1096, 980, 1392q
	n τ " 9 (Sc8)	100 9	7621	« 1.65 ˆ10 12	2674

  and formulated by:

		Bu pjq Bτ pjq " "	Bu pjq Bq pjq " r0 3 Rpq pjq , L, αqsJpq pjq , L, αq Bq pjq Bτ pjq ı Bq pjq Bτ pjq	(6.27)
	with	Bq pjq Bτ pjq deduced from (2.50) as:

  ¨¨¨N , k " 1 ¨¨¨n τ , and j " 1 ¨¨¨n S .

	Using the above exponential map formula (6.32),	Bg pjq BL	can be calculated as follows [66]:
	Bg pjq BL	"	Bg 1,s 1 pjq BL pjq 1,s 1 g pjq " IN b 1,s 2 ¨¨¨B ´gpjq 1,s 2 ¨¨¨g N,s 1 g pjq g pjq N,s 1 BL " IN b g pjq pjq N,s 2 ¯ı N,s 2 ı `gpjq `gpjq 1,s 1 1,s 1 g pjq Bg BL pjq 1,s 2 1,s 2 ¨¨¨g N,s 1 " IN b BL ´gpjq 2,s 1 ¨¨¨g N,s 1 g pjq N,s 2 ¯ı `¨¨g pjq pjq N,s 2 Bg pjq (6.33)
	where
						Bg i,s1 pjq BL	"	" 0 4 ¨¨¨α k g i pjq	p ξ i ¨¨¨0 4 pjq	ı
	and
				Bg i,s2 pjq BL	"	" 0 4 ¨¨¨p1 ´αk qg i pjq	p ξ i ¨¨¨0 4 pjq	ı
	With into (6.30) to obtain Bg pjq BL being calculated using (6.33), and Bu pjq BL , then, using (6.26), we finally deduce ∇ L h in (6.29). Bq pjq from (6.31), we implement both of them BL
	6.2.4.3 Calculation of ∇ T α
	As before, ∇ α h can be formulated as:
			∇ α h "	" Bh Bα		ı T	"	" Bh Bu	Bu Bα	ı T	"	j"n S ÿ j"1	"	Bh Bu pjq	Bu pjq Bα	T	(6.34)
	with this, using (6.20), we can write the following: Bh Bu being calculated from (6.26), thus, we only need to find the analytical form of	Bu Bα	. For
				Bu pjq Bα	" E 1	Bg pjq Bα	E nτ `r∇ u q pjq s T Bq pjq Bα	(6.35)
	Moreover, using (6.21), we can then have:
			Bq pjq Bα	" K	´1 ˆBH pjq Bα	τ pjq `BF pjq Bα	´BK Bα	pI nτ b q pjq q ˙(6.36)
												1,s1 g 1,s2 ¯¨¨¨´g pjq N,s1 g pjq N,s2 ¯(6.32) pjq
	with g i,s1 " e α k Li p pjq ξ pjq i , and g i,s2 " e p1´α k qLi p pjq ξ pjq i , for i " 1 123

Spτ, q ˚q " 2Φ T pτ, q ˚q 9 Φpτ, q ˚q

It eliminates wrong sub-spaces by the interval estimation technique, and discards preexplored configurations, yielding a significant decrease of the necessary iteration number to explore the whole workspace, which results in reducing the required computational time to achieve W E .

(a) Soft robot structure. (b) Oxz plan view of W E .

(a) Soft robot structure. (b) Oxz plan view of W E and BW E .

En raison de la complexité de calcul élevée de l'approche directe lors de l'estimation de l'espace de travail des robots mous, il est logique de proposer une approche qui peut discrétiser l'espace de l'effecteur final puisque sa dimension est plus petite et constante par rapport à l'espace des actionneurs pour les robots hyper-redondants, tels que les robots mous. De plus, il est plus efficace de cartographier uniquement la limite extérieure de l'espace de travail et d'éviter le calcul lourd de ses points intérieurs.Une telle méthode[START_REF] Amehri | Fem-based exterior workspace boundary estimation for soft robots via optimization[END_REF][START_REF] Walid | Discrete cosserat method for soft manipulators workspace estimation: An optimization-based approach[END_REF] consiste à discrétiser l'espace de l'effecteur final afin de cartographier la limite extérieure de l'espace de travail, et à déduire que la zone/volume délimitée par la limite de l'espace de travail est l'espace de travail accessible.L'approche proposée a été appliquée avec succès aux modèles PCS et FEM, où nous avons montré son efficacité à réduire la complexité et le temps de calcul nécessaires pour estimer l'espace de travail des robots mous, contrairement à l'approche directe qui explose exponentiellement lorsque la dimension des actionneurs augmente. D'autre part, comme cette approche ne consiste qu'à cartographier la limite extérieure de l'espace de travail en déterminant les points limites situés sur la surface entourant l'espace de travail, elle ne sera certainement pas en mesure d'identifier les informations internes de l'espace de travail telles que les limites intérieures.

Acknowledgments

First of all, I would like to express my thanks to my two supervisors, Gang Zheng and Alexandre Kruszewski for presenting me with the opportunity to work on such an interesting subject, and also for their guidance, help, and sharing their experience with me. I am grateful.

Thanks to the Region Hauts-de-France, the Project Inventor (I-SITE ULNE, le programme d'Investissements d'Avenir, Métropole Européenne de Lille), and the Project ROBOCOP [ANR PRCE 19 CE19] for the financial supporting of this thesis.

Many thanks to all members of DEFROST team, and especially Ke Wu for the help with the platform. Thanks to Shijie Li for his help with the validation part.

Many thanks to the reviewers of my manuscript and jury members. I would like also to express my gratitude to my family and friends for their care and support.

Firstly, we want to show that there are some cases where the optimal solution can only reach a desired point due to the length constraint. For this, we impose the segments length bounds in L " r0, 0.24s ˆr0, 0.16s ˆr0, 0.08sm.

Using the proposed optimization approach, we find the optimal segments length parameter L ˚" r0.225, 0.15, 0.075s T , and the optimal actuators length parameter α ˚" r0.89, 0.76, 0.10s T , which gives the optimal design shown by Fig. 6.7a. As illustrated by Fig. 6.7b, we can observe that the optimal solution allows only the reachability of one single point due to the restrictions imposed on the segments length parameter L. To achieve the accessibility to multiple points, we can select a larger possibility for the length bounds of each segment as L " r0, 0.3s ˆr0, 0.2sm.

Using the proposed optimization approach, we find the optimal segments length parameter L ˚" r0.29, 0.19, 0.09s T , and the optimal actuators length parameter α ˚" r0.88, 0.71, 0.17s T , which gives the optimal design as shown by Fig. 6.8a. Accordingly, the workspace of the In this specific scenario, the optimization-based approach is clearly the most efficient way to estimate the workspace of soft robots since it consists only in mapping the exterior boundary and avoiding the heavy computation of interior points, which also comes with the drawback that it cannot provide interior information to the workspace.

On the other hand, the interval analysis approach consists of exploring the whole workspace which allows us to have interior and exterior knowledge about the workspace, however, as a result, this approach is exhaustive.

Finally, the continuation approach allows the mapping of both interior and exterior boundaries to the workspace, this method allows us to identify the singularities of the robot which is useful for the design. However, this approach does not seek to find all achievable configurations of the workspace, but only the singularities of soft robots.

FEM-based Workspace Estimation Approaches Comparison

Based on the FEM modeling method, two different approaches have been proposed to evaluate the workspace of soft robots and were validated through the forward approach. Each approach has its advantages and disadvantages, and the following summarizes their characteristics.

For the sake of clarity and simplicity, we consider one specific scenario (e.g., the trunk-like soft robot actuated by four symmetric tendons). The structure of the robot and its corresponding results are given in Fig. 7.2. Then, we summarize the advantages, disadvantages, operations complexity, and computational time of each investigated workspace approach in Table 7.2. The interval analysis approach is computationally expensive compared to the optimizationbased approach. However, the interval analysis approach provides valuable interior and exterior information of the workspace.

On the other hand, the optimization-based approach is more efficient. However, this approach only maps the exterior boundary of the workspace ignoring its interior configurations.

In the following, several possible extensions of the obtained results will be discussed.

Perspectives

From Position-Access to Orientation-Access Workspace

This thesis has proposed different methodologies to estimate the position-access workspace of soft robots. However, the proposed approaches can also be extended to treat the orientation aspect of the workspace. Such an extension will come with necessary adaptations of the formulation of the workspace definition and the workspace set.

For the PCS model, we can treat the orientation aspect of the workspace by taking into consideration the orientation matrix R defined in (3.2) and then adapting the algorithms established in Chapters 3, 4, and 5 in a proper way for the purpose of studying the orientation aspect.

In the FEM model, we can treat the orientation aspect of the workspace by modifying the selection matrix C in (2.53). By considering 3 independent nodal positions around the end-effector q E 1 , q E 2 and q E 3 , which can be used to define a unique plan of the end-effector, we can always properly choose the matrix C P R 6ˆ3np such that q E now is a function of those 3 points and contains both position and orientation information of the defined plan.

From Tendons-driven Soft Robots to other Methods of Actuation

In this thesis, we have focused on tendon-driven soft robots in both the direct problem (i.e., workspace estimation) and inverse problem (i.e., design optimization). However, since the PCS method was applied to different methods of actuation (e.g., fluidic [START_REF] Renda | A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation[END_REF]), and that the FEM approach was also applied to different methods of actuation (e.g., pneumatic and hydraulic in [START_REF] Morales Bieze | Finite element method-based kinematics and closed-loop control of soft, continuum manipulators[END_REF][START_REF] Rodríguez | Real-time simulation of hydraulic components for interactive control of soft robots[END_REF]), then, we can also use different methods of actuation by adapting the actuation matrix H in both modeling methods to the adopted actuation manner.

Continuation Approach: Toward a Generalization of Tangent Vectors Calculations

The continuation approach was explained in detail in Section 5.2 of Chapter 5, where tangent vectors to boundary points met along the solution curve (5.7) are calculated in Section 5.2.3 based on a classification of those boundary points into three categories (see Table 5.1): regular boundary points (see Fig. 5.3), simple bifurcation points (see Fig. 5.4), and multiple bifurcation points (see Fig. 5.5).

The tangent vectors to each boundary point class were calculated using a different approach. On the one hand, the tangent vector of regular boundary points were calculated using Theorem 4 in Section 5.2.3.1 which consists of calculating the nullspace of the workspace boundary set Jacobian G Y . On the other hand, tangent vectors to simple bifurcation points was calculated in function of the parameters to be optimized, e.g., the robots' length L, the actuators length α, the radius of the robot r, and the distance of the actuators to the mid-line d, as follows:

Kpq, Xq q " Hpq, Xqτ `F pq, Xq (7.2a) q E " Cq (7.2b)

u " q E `up0q (7.2c) with X " rL, α, r, ds. The second step is to formulate the optimal problem in order to optimize the objective parameters X: Kpq pjq , Xq q pjq " Hpq pjq , Xqτ pjq `F ´qpjq , X ūpjq " Cq pjq `up0q

Finally, solving (7.3) yields the optimal parameters to achieve the optimal design of the investigated soft robot based on the FEM model.