On the Mathematical Modeling of Slender Biomedical Continuum Robots

The passive, mechanical adaptation of slender, deformable robots to their environment, whether the robot be made of hard materials or soft ones, makes them desirable as tools for medical procedures. Their reduced physical compliance can provide a form of embodied intelligence that allows the natural dynamics of interaction between the robot and its environment to guide the evolution of the combined robot-environment system. To design these systems, the problems of analysis, design optimization, control, and motion planning remain of great importance because, in general, the advantages afforded by increased mechanical compliance must be balanced against penalties such as slower dynamics, increased difficulty in the design of control systems, and greater kinematic uncertainty. The models that form the basis of these problems should be reasonably accurate yet not prohibitively expensive to formulate and solve. In this article, the state-of-the-art modeling techniques for continuum robots are reviewed and cast in a common language. Classical theories of mechanics are used to outline formal guidelines for the selection of appropriate degrees of freedom in models of continuum robots, both in terms of number and of quality, for geometrically nonlinear models built from the general family of one-dimensional rod models of continuum mechanics. Consideration is also given to the variety of actuators found in existing designs, the types of interaction that occur between continuum robots and their biomedical environments, the imposition of constraints on degrees of freedom, and to the numerical solution of the family of models under study. Finally, some open problems of modeling are discussed and future challenges are identified.

Introduction

Continuum robots use material deformation to move instead of joints. They may offer a technological solution to some of the difficult challenges of locomotion, perception, and manipulation found in a variety of unstructured and uncertain environments (Robinson and Davies, 1999). Biomedical applications have been a great motivator in the development of a wide variety of continuum and soft robots, ranging from surgery to therapy and other applications involving physical human-robot interaction. The great recent interest in these design paradigms stems from the observation that success in whatever form it is needed may be achieved without having complete control over the motion of a robot or its forces of interaction with the environment. In some cases, this is advantageous simply for reducing the complexity of engineered systems, and in other cases, performance may be increased beyond what is possible with rigid machines. Several excellent examples of this general principle come from tools of modern medicine. A flexible endoscope can navigate the intestines without a great degree of control over its own shape. The same is true for an intravascular catheter. In these examples, it is the particular combination of geometry and just the right amount of mechanical “softness” that facilitates the completion of the task. Beyond this snake-in-a-pipe approach to navigation, recent research has argued that physical compliance is advantageous in grasping, underwater swimming, robustness to collision, and locomotion on soft terrains where low ground pressure is required. The interested reader is referred to several review articles for a survey of the benefits, applications, challenges, and history of soft and continuum robots (Kim et al., 2013; Burgner-Kahrs et al., 2015; Walker et al., 2016; Cianchetti et al., 2018). Figure 1 shows four examples of continuum robot architectures which range from fully hard materials to fully soft and with composite structures in between these extremes.

FIGURE 1

FIGURE 1. (A) Concentric tube robots are comprised of hard (metallic) tubes which are precurved and nested inside one another. Rotating and translating the tubes results in motion. (B) Tendon-driven robots use one or more tendons or cables to provide internal actuation forces that bend a flexible, slender rod. (C) Pneumatic soft continuum robots use soft air muscles, which extend or contract with internal air pressure, to create bending in a composite structure. The supports could be hard or soft materials. (D) A fully soft pneumatic gripper uses asymmetry introduced by an inextensible fabric layer and an asymmetric air volume to create four slender fingers which bend to wrap around objects.

Though there is not universal agreement on definitions, the term continuum robot is generally used to imply that motion is generated without identifiable kinematic pairs, while the term soft robot implies at least a greater degree of mechanical compliance, defined as the ratio of displacement to force, exhibited in response to environmental forces than traditional approaches to robotic interaction. Many soft robots are made of soft materials, which may be characterized in terms of a material parameter such as the modulus of elasticity (Majidi, 2014). Continuum robots made of harder materials can be designed to exhibit high or low mechanical stiffness to external forces depending on the design details.

Continuum robots are classified as under-actuated mechanisms (Spong, 1998). This statement is taken to mean that in a practical sense, and within the context of a pre-defined scope of possible robot-environment interactions, more information than can be collected by a finite set of actuator-collocated sensors is needed to describe the shape and motion of the robot to the degree of accuracy demanded by engineering specifications or by the roboticist’s preference. The practical sense of the definition is emphasized since even rigid robots with revolute or prismatic pairs must deform to a small degree when interacting with their environment via forces. All mechanical systems are underactuated when there exist flexible modes that are not actuated but which should be controlled (Spong and Laurent, 1997). It is well known that the analysis of dynamics of underactuated robots is significantly more complex than for regular, fully actuated robots (Jain and Rodriguez, 1993).

Beyond being under-actuated, the modeler of a continuum manipulator also frequently faces other challenges. Designs are often difficult to separate into “components” since the structure and the actuator may be the same physical body. Actuators based on pneumatics, hydraulics, and composite structures may not be as easily characterized as electric motors. Friction and hysteresis models may be needed to explain observed mechanics, and environments rich with expected contacts may require the solution of contact models based on theories of nonlinear complementarity. Additionally, the standard kinematic descriptions based on the rigid transformations in the special Euclidean group $SE\left(3\right)$ are neither the most common approach to solid mechanics nor (necessarily) the most expedient approach to the description of solid continua undergoing deformation. With these considerations, one appreciates why the mathematical modeling of continuum and soft robots can be challenging.

This paper first reviews the state-of-the-art in the mathematical modeling of continuum manipulators having at least one “long” aspect in terms of its shape, which are termed slender in agreement with the mechanics literature. The goal of these models is to describe the dynamics (or statics) to relate actuator variables, other boundary conditions, and sensor measurements to the motion of the robot. The models are generally not concerned with other important aspects of robot design and analysis, such as repeatability, wear, safety, and other factors. For designs made of slender components, the motion of the robot is dominated by bending or beam-like deformations. This classification can be thought of as “arms,” “snakes,” or the individual “fingers” of a multi-fingered hand. Designs composed of individual components having this property are a natural extension, such as concentric tube robots (Mahoney et al., 2018) or multi-backbone continuum robots (Ding et al., 2013). For robots made of softer materials, such as the STIFF-FLOP designs, localized deformations may be complex, yet the dominant behavior is still beam-like (Cianchetti et al., 2014; Fraś et al., 2014). One of the goals of the work is to express the variety of methods encountered in the literature with a common notation. The review motivates a theoretical discussion rooted in the classical theories of solid mechanics. An analysis of the mechanics is used to support recommendations for future modeling efforts, with the conclusion that some choices for the model structure may result in better absolute model accuracy and efficiency (as quantified by the relationship between accuracy and dimensionality).

Review of the State of the Art

Table 1 presents the unified nomenclature that will be used throughout this paper. In the discussion of other works, the original nomenclature has been changed to match what is shown. There are three primary considerations in any physics-based approach to modeling of solid continua: the adoption of kinematic hypotheses and coordinates describing the configuration of the body, the application of the laws of mechanics, and the selection of mathematical models that describe the behavior of materials (Sadati et al., 2019). Kinematic hypotheses alone allow the modeler to describe the geometry of the robot, but this alone is insufficient for most purposes because it does not reveal which configurations are possible or likely. The mechanics, which are formulated naturally as partial differential equations, provide the relationships between the kinematic degrees of freedom that indicate which path of configurations will be taken if particular conditions (actuation, environments, etc.) are imposed. Finally, the material models are needed to close the relationship between the kinematic degrees of freedom and the kinetic quantities related by the mechanics.

TABLE 1

TABLE 1. Nomenclature used in this article.

Kinematic Descriptions

The forebears of continuum manipulators are the hyper-redundant robots, defined as those having a large (or infinite in the case of continuum robots) relative degree of redundancy (Chirikjian and Burdick, 1994). In any robot with material deformation which is substantial with regard to the kinematics or dynamics, both the relative degree of redundancy and the degree of under-actuation are theoretically infinite since the configuration space is infinite-dimensional. Here the usual definition of a robot configuration is used: “a complete specification of the location of every point on the robot” (Spong et al., 2006). There have been two primary methods to date of describing the configuration of continuum and soft robots: the curve-based description and the general continuum description.

The Curve-Based Description

The state of the art curve-based description is that of the special Cosserat rod (Antman, 2005). Figure 2 depicts the curve, its relationship to a solid body, and the quantities that are associated with the curve and the boundary conditions of a mechanical model. The elongated form of many continuum manipulators leads naturally to the concept of the “backbone curve,” which is typically defined to be a time-varying, piecewise differentiable curve in the standard three-dimensional affine Euclidean space $E$ with associated vector space E. A parametric representation gives the position of a point identified by a spatial parameter $s\in D\subset \mathrm{ℝ}$ at time $t\in \mathrm{ℝ}$ as a position vector ${\mathit{p}}_b\left(s\mathrm{,}t\right)\mathrm{\in }\mathbf{E}$ with respect to a specific frame of reference ${\mathrm{ℱ}}_0$ in $E$. The differentiability requirement on ${\mathit{p}}_b$ is always at least that the first derivative of ${\mathit{p}}_b$ with respect to $s$ exists, is piecewise continuous, and is nowhere equal to zero. This condition guarantees that the curve is rectifiable, or in other words has a measurable arc length (Kreyszig, 1991). The curve changes over time, modeling the motion of the robot, and it is presumed to describe the dominant features of the motion of the robot. Since there is no finite set of coordinates that describes every possible curve meeting these requirements, the description of the shape is infinite-dimensional.

FIGURE 2

FIGURE 2. Mathematical setup of the curve-based kinematic description of slender continuum robots.

The usual type of modeling hypothesis for slender bodies is that other points, which are not located on the backbone, are described by some auxiliary relationship that describes their positions relative to the positions on the backbone. The standard theories from beam mechanics may be adopted for this purpose, in which case the backbone curve may be affixed to the body at the neutral axis of bending¹. One example is the Euler-Bernoulli hypothesis, which states that sections normal to the backbone remain normal for all deformations. Another is the hypothesis due to Timoshenko stating that normal sections rotate relative to the backbone but remain planar. Standard “warping” theories can be used to couple motion of the points normal to the sections with twisting about the backbone if the sections are not circular. Regardless of these additional hypotheses, the curve is of fundamental importance to the kinematic description.

Explicitly, the body of the robot is identified by the curve through the consideration of a reference configuration $c_0$ of the robot. The backbone curve ${\mathit{p}}_b$ is placed on this reference configuration. The curve is then “framed” by a set of unit vectors ${\mathit{d}}_1\left(s,t\right), {\mathit{d}}_2\left(s,t\right),$ and ${\mathit{d}}_3\left(s,t\right)$ termed the director vectors. The first two are chosen to be orthogonal and to span the section of the body at $s$ which is normal to the curve. The third is taken to complete a right-handed, orthonormal coordinate frame as ${\mathit{d}}_3\left(s,t\right)={\mathit{d}}_1\left(s,t\right)\times {\mathit{d}}_2\left(s,t\right)$. In terms of classical differential geometry, ${\mathit{d}}_3$ is the tangent vector, and ${\mathit{d}}_1$ and ${\mathit{d}}_2$ could be selected as the normal and bi-normal vectors from Frenet’s formulas (Kreyszig, 1991). This procedure is problematic for general curves since torsion may be undefined, but many other alternative framings of the curve are possible which do not suffer this problem (Bishop, 1975). The backbone position and unit vectors together describe a local reference frame $\mathrm{ℱ}\left(s,t\right)$ for each point along the curve. The unit vectors equivalently define a spatiotemporal field of rotation operators ${\mathit{R}}_b\left(s,t\right)\in SO\left(3\right)$. The rotation field can be represented by matrices (Rucker and Webster, 2011), quaternions (Boyer et al., 2020), or any other suitable representation. Together with the position vector, a spatiotemporal field of transformations ${\mathit{g}}_b\left(s,t\right)\in SE\left(3\right)$ is defined by. ${\mathit{g}}_b\left(s,t\right)=\left\{{\mathit{R}}_b\left(s, t\right), {\mathit{p}}_b\left(s,t\right)\right\}.$

The vectors $\mathit{u}\left(s,t\right)=u_k{\mathit{d}}_k$ and $\mathit{v}\left(s,t\right)=v_k{\mathit{d}}_k$ are termed the “strain variables.” They describe deformation of the body and are invariant under rigid transformations. The vector $\mathit{u}$ has been widely called the “curvature” vector in the robotics literature, but this may be misleading since it is not generally the curvature of the deformed backbone curve. The term “flexural strain” is preferred for $u_1$ and $u_2$, and “torsional strain” for $u_3$. The variables $v_1$ and $v_2$ are called the shear strains, and $v_3$ is the dilation. The change in length or “extension” of the backbone curve is characterized by ${\partial }_s{\mathit{p}}^2=\mathit{v}\cdot \mathit{v}$. The strain variables are related to the framed curve by the following relationships.

$\begin{array}{c}{\partial }_s{\mathit{p}}_b\left(s,t\right)=v\left(s,t\right), {\partial }_s{\mathit{d}}_k\left(s,t\right)=\mathit{u}\left(s,t\right)\times {\mathit{d}}_k\left(s,t\right)\end{array}.(1)$

Finally, the vectors ${\partial }_t{\mathit{p}}_b\left(s,t\right)$ and $\mathit{\omega }\left(s,t\right)$ represent the linear and angular velocity of the backbone curve and director vectors. The angular velocity satisfies the equation ${\partial }_t{\mathit{d}}_k\left(s,t\right)=\mathit{\omega }\left(s,t\right)\times {\mathit{d}}_k\left(s,t\right)$. The four functions $\mathit{u},\mathit{v},{\partial }_t{\mathit{p}}_b,$ and $\mathit{\omega }$ are not independent; they must satisfy ${\partial }_s\mathit{\omega }={\partial }_t\mathit{u}+\mathit{u}\times \mathit{\omega }$. In the reference configuration, the flexure strains have non-zero values ${\mathit{u}}_0\left(s\right)$ if the backbone is not a straight. Generally, ${\mathit{v}}_0\left(s\right)={\mathit{d}}_3\left(s\right)$, but other choices are possible.

The General Continuum Description

The second approach to describing the configuration of continuum robots is to make as few prior kinematic hypotheses on the configuration as possible. The traditional description of a three-dimensional continuum in solid mechanics is used in this case. In this approach a reference configuration $c_0$ is identified by their position vector relative to a frame of reference ${\mathrm{ℱ}}_0$. Three coordinates $\mathit{X}\in {\mathrm{ℝ}}^3$ identify the position of each point in the body via a one-to-one, differentiable vector-valued function $\mathit{P}\left(\mathit{X}\right)$. If $\mathit{X}$ is chosen as the Cartesian coordinates with respect to ${\mathrm{ℱ}}_0$, then this function and its inverse are trivial. The final locations of the points are described by $\mathit{p}\left(\mathit{X},t\right)$. In some cases, it is useful to define a displacement field $\mathit{U}$ as follows.

$\mathit{p}\left(\mathit{X},\mathit{ }t\right)=\mathit{P}\left(\mathit{X}\right)+\mathit{U}\left(\mathit{X},\mathit{ }t\right).$

The amount of stretching can be quantified by the deformation gradient, defined by

$\mathit{F}\left(\mathit{X},t\right)={\frac{\partial \mathit{p}}{\partial \mathit{P}}|}_{\mathit{X},t}.$

The deformation gradient straightforwardly describes the local changes in length (amount of stretching) and therefore plays a major role in the definition of strain measures. Note also that the curve-based description of the configuration, together with the classical Euler-Bernoulli hypothesis, can be placed into this more general framework using $\mathit{X}=\left(s,X_2,X_3\right)$ and $\mathit{p}\left(\mathit{X},t\right)={\mathit{p}}_b\left(s, t\right)+X_2{\mathit{d}}_2\left(s,t\right)+X_3{\mathit{d}}_3\left(s,t\right)$ (Antman, 2005).

Perspective on Discretization and Configuration Spaces

There are two perspectives that one might take when describing the kinematics or mechanics of continua. In the first perspective, the model consists of a (possibly nonlinear) PDE, a domain on which the PDE applies, and boundary conditions in the form of constraints or measurements. The robot’s state space consists of the dependent variables related by the PDE. The state space is therefore a particular Cartesian product space that might involve, in general, both finite-dimensional spaces and infinite-dimensional function spaces. In the process of computing a numerical solution to a model, any part of the state that belongs to an infinite-dimensional space must be approximated by a finite set of coordinates in ${\mathrm{ℝ}}^n$, but the choice of coordinates does not need to be of great concern to the modeler. This perspective has been taken by numerous authors for general continuum manipulators (Trivedi et al., 2008; Till et al., 2019), concentric tube robots (Dupont et al., 2010; Rucker et al., 2010; Gilbert et al., 2016), parallel continuum robots (Black et al., 2018), and bioinspired locomotion by snakes and worms (Boyer et al., 2012). The modeler hopes that any approximation error is small enough to be ignored, and error-controlled numerical methods may provide some assurances. This first perspective is the natural one if, for example, the modeler selects an error-controlled, automatic step-size numerical integrator like the Dormand-Prince Runge Kutta pair to approximate the solution to a differential equation with a spatially distributed independent variable. The benefit to this perspective is that questions of convergence may generally be avoided. However, there are two main disadvantages: first, there is a relative paucity of tools available if the problem is not expressed with respect to a single independent spatial variable; second, the degrees of freedom chosen by automatic numerical methods may be unknowable in advance and may vary between model solutions, making it difficult to apply algorithms built on spaces like ${\mathrm{ℝ}}^n$ or on manifolds where coordinate charts are available.

In a second perspective, the equations of an infinite-dimensional model are explicitly discretized through a suitable method such as the finite element method or a finite difference method (Renda et al., 2014; Back et al., 2015; Gilbert and Godage, 2019) or via a spectral method involving a “modal” decomposition (Chirikjian and Burdick, 1994; Godage et al., 2015; Chen et al., 2020). In this perspective, the modeler takes control over the discretization and fixes the dimensionality of the resulting model. One is free to take the perspective that a new model has been created that is not necessarily subordinate in any way to the infinite-dimensional model. In other words, the infinite dimensional dependent variables, ODEs, and/or PDEs, were only a steppingstone to the finite-dimensional model. The dimension may be varied according to a model hyper-parameter $N$, and often one wishes that as $N\to \infty$, the solutions to the sequence of fixed-dimensional models approach the solution to a corresponding infinite-dimensional model.

The second perspective is the standard one in generally accepted theories of robot kinematics and dynamics, in which the goal is to find a suitable coordinate set that describes the displacement field $\mathit{u}\left(\mathit{X},t\right)$ that takes a material point located at initial position $\mathit{P}$ to its final position $\mathit{p}=\mathit{P}+\mathit{u}$. With rigid link manipulators, the space is partitioned into non-intersecting bodies indexed by number $i\in {\mathrm{\mathbb{Z}}}^{+}$ and equipped with local coordinate frames, and then the machinery of $SE\left(3\right)$ is used to associate each body with its own displacement field expressed in terms of one of the coordinate transformations ${}^{\mathrm{0}}{{\displaystyle \mathit{T}}}_i\in SE\left(3\right)$ representing the transformation between the base frame $0$ and the frame of the $i^{th}$ body. For serial, rigid-link robots, the choice of finite dimensional coordinates parameterizing the displacement field is usually one of two conventions, the Denavit-Hartenberg convention (Denavit and Hartenberg, 1955) or the twist coordinate system and product-of-exponentials formula (Brockett, 1984).

For continuum and soft robots, neither the perspective (finite vs. infinite-dimensional) nor the approach to discretization (choice of coordinates) appears to be standardized. In some cases, restrictive assumptions do allow a set of finite coordinates that uniquely specify the configuration of a continuum robot. For example, Bretl and McCarthy showed that for the Kirchhoff rod with no external loading, a configuration space isomorphic to ${\mathrm{ℝ}}^6$ can be selected, corresponding physically to a known internal force and moment at the same location in space as the known orientation of the rod (Bretl and McCarthy, 2014). A similar result is known for coordinates of the configuration space of concentric tube robots without any external loads (Gilbert et al., 2016). The general principle is a basic result on initial value problems. If the mechanics of the system can be modeled by a system of $n$ first-order initial value problems,

${\partial }_s\mathit{y}=\mathit{F}\left(s,\mathit{y}\right), \mathit{y}\left(s_0,t\right)={\mathit{y}}_0\left(t\right),$

with $F$ uniformly Lipschitz in $\mathit{y}$ and continuous in $s$, then the solutions are uniquely determined by ${\mathit{y}}_0$ (Schaeffer and Cain, 2016). Therefore, if all state information of the robot is contained in the functions $y_i\left(s, t\right)$, then it is clear that ${\mathit{y}}_0$ is a suitable set of coordinates for the configuration space of the robot. For curve-based models ${\mathit{y}}_0$ usually belongs to a space of the form $SE{\left(3\right)}^r\times {\mathrm{ℝ}}^n$.

However, with less restrictive assumptions, low-dimensional configuration spaces are not generally found. Such is the case for parallel continuum robots (Black et al., 2018), for growing robots (Greer et al., 2019), or soft robotic hands (Schlagenhauf et al., 2018). It is in general impossible to find a “minimal” set of coordinates for the C-space of any continuum manipulator when the locations and nature of external loads or contacts are a-priori unknown and when these loads cause substantial changes in the robot shape. The subsections that follow describe a variety of methods that have been used to mathematically represent the configurations of continuum robots.

Spectral Methods

Spectral methods were some of the earliest described methods for the kinematic modeling of backbone curves. In this method, the configuration is represented by a finite number of coordinates $\mathit{q}\left(t\right)\in {\mathrm{ℝ}}^N$ by assuming that some kinematic quantity is described by a truncated sum of “modal” shape functions ${\varphi }_i\left(s\right)$ in a manner analogous to a Fourier series. The general form is to have a scalar quantity $S_j$ represented as

$S_j\left(s,t\right)={\displaystyle \sum }_{i=1}^N{a}_{ij}\left(\mathit{q}\left(t\right)\right){\varphi }_i\left(s\right).$

The function $a_{ij}$ may be simply an index into the vector $\mathit{q}$ pulling out one of the components, or it may be a more complicated relationship. The mode shapes are generally selected among one of the standard families such as trigonometric functions $\sin \left(k_i\pi s\right)$ and $\cos \left(k_i\pi s\right)$ for a series of values $k_i\in \mathrm{ℝ}$ (directly analogous to a truncated Fourier series), the standard monomials $\left\{1,s,s^2,\dots \right\}$, the Legendre polynomials, Chebyshev polynomials, etc. (Chirikjian and Burdick, 1994; Zhang and Simaan, 2013; Chen et al., 2020). In general, to be classified as a spectral method, the mode functions should have global support rather than local support, which leads to the element-based methods described below.

There is a great deal of freedom within this approach. For example, the tangent vector ${\mathit{d}}_3$ can be expressed in spherical angles $\theta \left(s,t\right)$ and $\varphi \left(s, t\right)$, and then $S_1=\theta$ and $S_2=\varphi$, and $\mathit{v}=\left(\mathrm{0,0,1}\right)$ completes the kinematic description (Chirikjian and Burdick, 1994). $S_j$ could also be chosen directly as a component of the displacement field of the backbone curve (Godage et al., 2015). These methods are extrinsic because they seek to approximate kinematic quantities as measured by the observer in the inertial frame ${\mathrm{ℱ}}_0$. Parameterizations also possible which represent the strain variables $\mathit{u}\left(s,t\right)$ and $\mathit{v}\left(s,t\right)$ measured by an observer in the local frame $\mathrm{ℱ}\left(s,t\right)$ (Boyer et al., 2020). When coupled with a collocation method used to determine $\mathit{u}$, it was shown that the Magnus expansion can be used to efficiently recover the position and orientation field (Orekhov and Simaan, 2020).

In the context of continuum robots, to the best of the author’s knowledge, the spectral methods have only been applied in conjunction with the curve-based descriptions discussed in The Curve-Based Description and not for more general continuum descriptions.

Element-Based Methods: PCC

The element-based methods, in contrast to the spectral methods, break up the problem spatially into adjacent sub-domains and attempt to model the kinematics on each sub-domain using a simpler hypothesis. This procedure can be carried out for both the curve-based description and the general continuum description. Many authors have adopted the kinematic hypothesis that the backbone curve is a sequence of circular arcs which are concatenated by imposing tangency conditions. There is a natural extension of this idea to piecewise helical curves. This approximation is termed the “piecewise constant curvature” (PCC) method, and many continuum robots have even been designed to exhibit deformation of this kind, at least in the absence of external loads (Webster and Jones, 2010). For example, multi-backbone robots and tendon-driven robots will adopt, with actuation, shapes very close to circular arcs with appropriate design decisions (Camarillo et al., 2008; Xu and Simaan 2008). On the other hand, even gravitational loading may cause more flexible robots to adopt shapes more complex than a single circular arc (Trivedi et al., 2008).

Within the curve-based framework described above, the “standard” PCC hypothesis including inextensibility and shear-lessness is equivalent to a partitioning of the domain into $m$ elements ${\text{Γ}}_e=\left[s_{e-1}, s_e\right]$ with $D={\cup }_{e=1}^m{\text{Γ}}_e$ and an approximation of the flexural-torsional strain as the following sum:

$\mathit{u}\left(s,t\right)={\displaystyle \sum }_{e=1}^m{\mathit{u}}_e\left(t\right) {\text{χ}}_{{\text{Γ}}_e}\left(s\right) , {\chi }_{{\text{Γ}}_e}\left(s\right)=\{\begin{array}{cc}1& s\in {\text{Γ}}_e\\ 0& s\notin {\text{Γ}}_e\end{array}\mathit{v}\left(s,t\right)={\mathit{d}}_3\left(s\right).$

A similar approach is called the piecewise constant strain (PCS) method and extends the definition to include an approximating sum representing $\mathit{v}$ (Renda et al., 2016). The vectors are resolved in Cartesian coordinates of the frame $\mathrm{ℱ}\left(s,t\right)$ and it is these body-frame coordinates that are presumed constant over the element.

The major advantage of the PCC method is that the extrinsic variables ${\mathit{p}}_b$ and ${\mathit{d}}_k$ (equivalently ${\mathit{R}}_b$) are easily computed using an explicit recursion presuming that a single boundary condition on the pose of the curve, $\mathit{g}\left(s_0, t\right)={\mathit{g}}_0\left(t\right)$, is known. On element ${\text{Γ}}_i$, the transformation from ${\mathrm{ℱ}}_0$ to $\mathrm{ℱ}\left(s,t\right)$ is given by

$\mathit{g}\left(s, t\right)=\mathit{g}\left(s_{e-1}, t\right)\exp \left(\left(s-s_{e-1}\right){\mathit{\xi }}_e\left(t\right)\right),{\mathit{\xi }}_e\left(t\right)={}^{\mathrm{ℱ}\left(s,t\right)}\left[\begin{array}{cc}{\mathit{u}}_e{\left(t\right)}_{\times }& {\mathit{v}}_e\left(t\right)\\ 0& 0\end{array}\right].$

The mapping $\text{exp}:se\left(3\right)\to SE\left(3\right)$ is continuous and differentiable, including at the element $0\in se\left(3\right)$, and it may be computed in closed form, below in vector form but originally discovered by Euler and Rodrigues (Rodrigues, 1840; Cheng and Gupta, 1989).

$\exp \left(\left[\begin{array}{cc}{\mathit{\omega }}_{\times }& \mathit{v}\\ 0& 0\end{array}\right]\right)=\left[\begin{array}{cc}\mathit{R}\left(\mathit{\omega }\right)& \mathit{p}\left(\mathit{\omega }, \mathit{v}\right)\\ 0& 1\end{array}\right]\mathit{R}\left(\mathit{\omega }\right)=I+C\left(\mathit{\omega }\right){\mathit{\omega }}_{\times }+A\left(\mathit{\omega }\right){\mathit{\omega }}_{\times }^2 \mathit{p}\left(\mathit{\omega }, \mathit{v}\right)=C\left(\mathit{\omega }\right)\mathit{v}+A\left(\mathit{\omega }\right)\left(\mathit{\omega }\times \mathit{v}\right)+B\left(\mathit{\omega }\right)\left(\mathit{\omega }\cdot \mathit{v}\right)\mathit{\omega }\mathit{.}$

The coefficient functions are as follows, with $\text{sinc}\left(x\right)$ the un-normalized sinc function, which is available in the standard libraries of many programming languages.

$A\left(\omega \right)=\frac{1}{2}{\text{sinc}}^2\left(\frac{\Vert \omega \Vert }{2}\right), B\left(\omega \right)=\frac{\Vert \omega \Vert -\sin (\Vert \omega \Vert )}{{\Vert \omega \Vert }^3}, C\left(\omega \right)=\text{sinc}(\Vert \omega \Vert ).$

In IEEE-754 double-precision arithmetic, the author has found that $B\left(\omega \right)$ may be accurately computed by a truncated Taylor series $B\left(\omega \right)\approx 1/6-{\Vert \omega \Vert }^2/120+{\Vert \omega \Vert }^4/5040$ if $\Vert \omega \Vert <1\times {10}^{-4}$.

The PCC and/or PCS methods are the simplest explicit and consistent discretization methods for a framed curve which interpolate the intrinsic (strain) variables rather than the extrinsic (position and orientation) variables. In other words, given a curve with bounded flexural, torsional, shear, and extensional strains, the error between the curve and its PCC or PCS approximation (both in terms of ${\mathit{d}}_k$ and in terms of ${\mathit{p}}_b)$ shrinks as the number of elements increases. This property is important because it means that the PCC/PCS framework can describe robots with practically any backbone curve if the domain is broken into enough elements. Figure 3 depicts a single isotropic rod under a combined twisting and bending moment, resulting in a helical shape. The results of a PCC approximation under a linear elastic material law and subject to the virtual work principle discussed in Projection via D’Alembert’s Principle below are shown in Figure 4. The flexural strain variables of the PCC model approach those of the exact solution as the number of elements increases. The example also demonstrates an important distinction between material flexural strain and the usual notion of curvature of a shape. The exact solution under the end load depicted in Figure 3 is indeed a helix, which is a “constant curvature” shape, yet the flexural strain variables, when resolved in components of the material frame $\mathrm{ℱ}\left(s,t\right)$ or any fixed frame of reference, are not constant functions of the arc length. For this reason and others discussed in Considerations for Kinematic Hypotheses, the simplicity of the PCC approach relative to others may not outweigh the drawbacks.

FIGURE 3

FIGURE 3. Simulation of a cantilevered rod under a combined bending and twisting concentrated moment, forming a helix.

FIGURE 4

FIGURE 4. Convergence of the PCC discretization to the exact flexural strains of the helical rod shape depicted in Figure 3. Note that the exact flexural strain components are not constant functions of arc length.

Element-Based Methods: Higher-Order One-Dimensional

The discretization of the curve into elements can also be accomplished with higher-order schemes than the piecewise constant strain approach. In general, as with the PCC approach if a boundary condition $\mathit{g}\left(s_0, t\right)$ is provided at a single location along the backbone curve, then the C-space can be defined by approximating the strain variables on each element with an expansion having more terms than the PCS approach (Boyer et al., 2020).

${}^{\mathrm{ℱ}\left(s,t\right)}\left[\begin{array}{c}\mathit{u}\left(s,t\right)\\ \mathit{v}\left(s,t\right)\end{array}\right]=\underset{e=1}{\overset{m}{\cup }}{\displaystyle \sum _{j=1}^{N_e}{\varphi }_j\left(z_e\right)q_{ej}\left(t\right)}, z_e=\frac{s-s_e}{s_e-s_{e-1}}.$

The union operator is abused to mean here that the element-local terms contained within each argument of the union are “assigned” to provide the evaluation on that element. The coordinate $z_e$ is a normalized element-local length coordinate identifying the cross section within the element. If the model is to contain flexure, torsion, shear, and extension, then a polynomial expansion of order $p$ in each component entails $N_e=6\left(p+1\right)$ degrees of freedom within each element. Continuity conditions at the element boundaries may and/or boundary conditions may reduce the number of total degrees of freedom from $m{N}_e$ as implied by the formula.

Other representations for the C-space based on segmentation into elements are possible but have not been widely pursued within the robotics community. Certainly, the more widely adopted approach within the finite element and computational mechanics communities has been to make the primary kinematic hypothesis at the level of $\mathit{p}$ and/or ${\mathit{d}}_k$ (Simo and Vu-Quoc, 1986), and this has also been studied in the context of continuum robots (Sadati et al., 2019). This is the approach taken, for example, by the ANSYS simulation software in handling nonlinear beam elements via the BEAM188 and BEAM189 elements. These methods require interpolation on the rotation group $SO\left(3\right)$, and some care is required to ensure that the formulation is invariant to rigid-body transformations (Crisfield and Jelenić, 1999). One major advantage of strain interpolation is that it is frame-invariant directly by construction; the major disadvantage is that the calculation of inertial forces is greatly complicated by the spatial coupling of the degrees of freedom.

Element-Based Methods: General Continuum

More general finite-element descriptions have also been used to model soft and continuum robots. In this case, the degrees of freedom $\mathit{q}\left(t\right)$ directly interpolate the position field $\mathit{p}(\mathit{X}, t)$ over the three-dimensional domain of the body. Using typical first-order (linear) interpolants, the degrees of freedom are the Cartesian coordinates of the nodes of the mesh that breaks the body into discrete volumes. Direct nodal position discretization using finite elements can be used for closed-loop control using a dimensionality reduction scheme based on projection (Bieze et al., 2018). It has also been shown that high-order FEM models with an order reduction method involving fitting to PCC kinematics is effective (Runge et al., 2017). Finally, it has been demonstrated that general 3D FEM with model order reduction based on the Proper Orthogonal Decomposition can produce models amenable to dynamic closed-loop control (Katzschmann et al., 2019b).

Direct Nodal Discretization

Closely related to the element-based methods are those based on direct discretization of the variables. Differential operators in the mechanics can be replaced by their equivalent finite-difference operators to form algebraic equations directly, operating on the values of field variables specified at discrete spatial locations $s_i$ for the curve-based approach. A finite difference scheme applied directly to the geometrically exact Cosserat equations has been demonstrated for both the planar case and the spatial case (Hasanzadeh and Janabi-Sharifi, 2014; Gilbert and Godage, 2019; Wang et al., 2021) and described for concentric tube robots (Webster and Rucker, 2009). Finite-difference methods were also used with direct spatial discretization to model a soft underwater arm driven by cables (Renda et al., 2014).

Pseudo-Rigid Body Methods

The pseudo-rigid body methods replace the continuum with an approximating rigid linkage. If the curve is broken into a sequence of chords with rotational joints at the nodes joining the chords, then this is equivalent to a spatial “lumping” of the flexural strains into a discrete point via the use of the Dirac delta distribution (Chirikjian and Burdick, 1991; Greigarn et al., 2019).

${}^{\mathrm{ℱ}\left(s,t\right)}\mathit{u}\left(s,t\right)={\displaystyle \sum _{i=1}^m{q}_i\left(t\right)\delta \left(s-s_i\right){\mathit{n}}_i}.$

A universal joint is the result if two orthogonal axes ${\mathit{n}}_i$ and ${\mathit{n}}_{i+1}$ are placed with $s_i=s_{i+1}$ with both axes normal to the backbone curve. Three orthogonal axes create a spherical joint.

It has been shown that the kinematics of tip-loaded cantilever beams can be modeled adequately by a serial 3R mechanism (Su, 2009). Other PRB models have been created for modeling of catheters (Ganji and Janabi-Sharifi, 2009), tendon-driven continuum manipulators for minimally invasive surgery (Penning and Zinn, 2014), and MRI-actuated catheters (Greigarn et al., 2017). A 6-DOF PRB segment model has also been proposed (Venkiteswaran et al., 2019). An equivalence has also been shown between the coordinates of a PCC model and a suitably defined pseudo-rigid body model, indicating that PRB model segments with RPPR kinematics can be used to describe the same configuration space as PCC models (Katzschmann et al., 2019a).

Initial Value Problem Concepts

There are additionally a variety of other methods of analysis and computation which do not explicitly select the degrees of freedom in the kinematic description. In these methods, the unknowns are conceptually left as unknown functions, and numerical methods are used which automatically select the degrees of freedom used to represent the unknown functions, usually via an error estimation and control algorithm.

These methods have been used when the problem is re-cast as a one-dimensional boundary value problem with split boundary conditions.

${\partial }_s\mathit{y}=f\left(s,\mathit{y}\right), G_a\left(\mathit{y}\left(0,t\right), t_j\right)=0, G_b\left(\mathit{y}\left(L,t\right), t_j\right)=0.$

Solutions can then be provided by numerical codes which automatically determine the degrees of freedom used to approximate the function $\mathit{y}\left(s, t_j\right)$ for each discrete value of $t_j$. For continuum robots these methods have been demonstrated via collocation (Webster and Rucker, 2009) and shooting methods (Till et al., 2015; Mauzé et al., 2020) using numerical tools that approximate $\mathit{y}\left(s, t_j\right)$ via piecewise polynomials. It has also been shown recently that the dynamics problem for a wide variety of architectures based on single or multiple Cosserat rod sub-models can be cast as a shooting problem on an ODE once the time derivatives have been discretized using finite differences (Till et al., 2019).

Differential Kinematics for Strain-Variable Hypotheses

It is often necessary to calculate a manipulator “Jacobian field” based on the curve parametrization, and if the generalized coordinates are defined to interpolate the strain variables, this field is not trivial to calculate.

$\begin{array}{c}\left[\begin{array}{c}{\partial }_t\mathit{p}\\ \mathit{\omega }\end{array}\right]\left(s,t\right)={\mathit{J}}_{\mathit{q}}\left(s,t\right){\partial }_t\mathit{q}\left(\mathit{t}\right)=\left[\begin{array}{c}{\mathit{J}}_{\mathit{p}}\\ {\mathit{J}}_{\mathit{\omega }}\end{array}\right]{\partial }_t\mathit{q}\left(t\right)\end{array}.(2)$

Letting ${\mathit{J}}_i$ be the column multiplied by ${\partial }_t{q}_i\left(t\right)$, then the column can be calculated from the following differential relationships:

$\begin{array}{l}{\partial }_s{\mathit{J}}_{pi}={\partial }_{q_i}\mathit{v}\\ {\partial }_s{\mathit{J}}_{\omega i}={\partial }_{q_i}\mathit{u}+\mathit{u}\times {\mathit{J}}_{\mathit{\omega i}}.\end{array}$

One must take care when the interpolation is carried out on the strain variables in coordinates of the local frame $\mathrm{ℱ}\left(s,t\right)$. If desired, the coordinates in the body frame representations ${}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{J}}}_{pi}$ and ${}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{J}}}_{\omega i}$ may be calculated from the appropriate representation of these equations in the moving frame (Rucker and Webster, 2011).

${\partial }_s \left({}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{J}}}_{pi}\right)=-{}^{\mathrm{ℱ}\left(s,t\right)}\mathit{u}\times {}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{J}}}_{pi}+{\partial }_{q_i}\left({}^{\mathrm{ℱ}\left(s,t\right)}\mathit{v}\right)+{}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{J}}}_{\omega i}\times {}^{\mathrm{ℱ}\left(s,t\right)}\mathit{v}{\partial }_s \left({}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{J}}}_{\omega i}\right)=-{}^{\mathrm{ℱ}\left(s,t\right)}\mathit{u}\times {}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{J}}}_{\omega i}+{\partial }_{q_i}\left({}^{\mathrm{ℱ}\left(s,t\right)}\mathit{u}\right).$

From a known boundary condition where ${\mathit{J}}_{pi}\left(0,t\right)=0$ and ${\mathit{J}}_{\omega i}\left(0,t\right)=0$, the solution to these equations can be expressed in closed form as the solution to a linear time varying system.

$\left[\begin{array}{c}{\mathit{J}}_{pi}\\ {\mathit{J}}_{\omega i}\end{array}\right]\left(s\right)={\displaystyle \underset{0}{\overset{s}{\int }}\exp \left(-{\displaystyle \underset{0}{\overset{r}{\int }}{\text{ad}}_{\xi \left(\tau \right)}\text{d}\tau }\right)}{\partial }_{\mathit{q}}\left[\begin{array}{c}{}^{\mathrm{ℱ}\left(s,t\right)}\mathit{v}\\ {}^{\mathrm{ℱ}\left(s,t\right)}\mathit{u}\end{array}\right]\text{d}r, {\text{ad}}_{\xi \left(\tau \right)}={}^{\mathrm{ℱ}\left(s,t\right)}\left[\begin{array}{cc}{\mathit{u}}_{\times }& {\mathit{v}}_{\times }\\ 0& {\mathit{u}}_{\times }\end{array} \right].$

Mechanics

Regardless of how the shape of a robot is described, the principles of classical mechanics are frequently used to describe the relationships between the model’s degrees of freedom, the internal stresses, and any imposed boundary conditions which may include external forces, imposed positions or orientations of parts of the robot, contact conditions. The robot’s actuators may generally be modeled in one of two ways: either they are described as constraints (a form of boundary condition) or as sources of internal stress.

The Equations of Motion for the Special Theory of Cosserat Rods

In the curve-based description, the equations of motion of the special theory of Cosserat rods serve as the strong form differential equations governing the mechanics (Antman, 2005).

$\begin{array}{c}{\partial }_s\mathit{n}+\mathit{f}=\rho A{\partial }_{tt}\mathit{p}+\rho I_k{\partial }_{tt}{\mathit{d}}_k\end{array}(3)$

$\begin{array}{c}{\partial }_s\mathit{m}+{\partial }_s\mathit{p}\times \mathit{n}+\mathbf{\ell }=\rho I_k{\mathit{d}}_k\times {\partial }_{tt}\mathit{p}+{\partial }_t\left(\rho \mathit{J}\cdot \mathit{\omega }\right)\end{array}.(4)$

The sum from $k=1$ to 3 is implied over the terms involving $I_k$ and ${\mathit{d}}_k$. The variables $\mathit{n}\left(s,t\right)$ and $\mathit{m}\left(s,t\right)$ are the internal force and the internal moment, which are interpreted as the resultant force and resultant moment of the stress acting on section $s$. In the case of a slowly accelerating body, which is typical in many biomedical applications, a quasistatic approximation may be used, in which all terms on the right-hand side are neglected (Burgner-Kahrs et al., 2015). $\mathit{f}$ and $\ell$ are externally applied forces and moments. Applied concentrated forces and moments require the Dirac $\delta$ distribution to express in this formulation.

In the case of a model which allows freedom in all the strain variables, $\mathit{m}$ and $\mathit{n}$ are algebraically related to the kinematic variables through a suitable material constitutive law. On the other hand, in the shear-less and extension-less model, $\mathit{n}$ is a basic unknown and is equivalent to a Lagrange multiplier which enforces the constraint $\mathit{v}\left(s,t\right)={\mathit{v}}_0\left(s\right)$.

The parameter $\rho A$ is the mass density (expressed per unit length) of the cross-section. $\rho \mathit{J}$ is the mass moment of inertia (per unit length) of the section, which makes $\rho \mathit{J}\cdot \mathit{\omega }$ the angular momentum (per unit length) calculated about the mass center of the section. The three parameters $\rho I_k$ account for linear momentum density of the cross section caused by angular velocity of the curve. The author is not aware of any works in the robotics literature for which this term has been nonzero; if the backbone curve is chosen to pass through the mass centers of the cross sections, then $\rho I_k=0$ and the equations simplify considerably. However, it is noteworthy that this may in general result in the curve failing to pass through the cross-section centroids (if multiple materials are used) or it may be impossible to satisfy this requirement exactly if a single curve is used to model a body with complex geometric features.

The Equations of Motion for Pseudo-Rigid Body Models

With the PRB-type models, the equations of motion are exactly those of a classical multibody dynamical system with scleronomous, holonomic constraints. These equations are commonly given as follows (Murray et al., 2017).

$\mathit{M}\left(\mathit{q}\right)\ddot{\mathit{q}}+\mathit{C}\left(\mathit{q},\dot{\mathit{q}}\right)\dot{\mathit{q}}+\mathit{N}\left(\mathit{q},\dot{\mathit{q}}\right)=B\left(\mathit{\tau }\right).$

The right-hand side contains the non-conservative generalized forces associated with actuation and any other forces; since the robots are underactuated there are generally many more rows in this equation than actuator variables ${\tau }_i$. Also, it is noteworthy that the inertial forces are not trivial to calculate since the motion of the continuum body is not the same as the motion of the rigid-link approximation. Some assumptions about how the continuum “tracks” the rigid-link approximation as it moves is needed. One approach is to match the centers of mass of chords along the curves of a PCC model with centers of mass of the links in the rigid link model (Della Santina et al., 2018).

The Equations of Motion for General Deformable Bodies

The dynamic equilibrium conditions of classical continuum mechanics serve as the defining relationship for general three-dimensional finite element models of soft and continuum robots. Rarely are these equations encountered explicitly in the literature on continuum robots, with most authors preferring to state the result after the strong form equations have been converted to the weak form and integrated. The resulting equations, incorporating constraint forces, are of the following form (Goury et al., 2021).

$\mathit{M}\left(\mathit{q}\right)\ddot{\mathit{q}}+\mathit{F}\left(\mathit{q},\dot{\mathit{q}}\right)+\mathit{G}\left(\mathit{q}\right)={\mathit{H}}^T\mathit{\lambda }$

The form of this equation is directly analogous to the classical form of the dynamical equations for rigid multibody systems. $\mathit{M}\left(\mathit{q}\right)\ddot{\mathit{q}}$ accounts for the inertial forces, $\mathit{F}\left(\mathit{q},\dot{\mathit{q}}\right)$ accounts for the internal forces produced by deformation of the material, and $\mathit{G}\left(\mathit{q}\right)$ accounts for gravitational forces. The matrix $\mathit{H}$ is associated with the constraints and boundary conditions and encodes the effect of the boundary and actuation forces contained in the vector $\mathit{\lambda }$. The details of the construction procedure for this equation are out of the scope of this paper.

Projection via D’Alembert’s Principle

In the case of the curve-based models using either the PCS or higher-order models, the equations can be projected onto the degrees of freedom of the model using Galerkin’s principle, probably better known among mechanical engineers as the principle of virtual work (Greenwood, 1988). The method is also equivalent in results to Kane’s method of virtual power (Kane and Levinson, 1983; Rone and Ben-Tzvi, 2014). Because the backbone curve descriptions for the PCC, PCS, and higher order strain variable interpolants are described by independent degrees of freedom $\mathit{q}\in {\mathrm{ℝ}}^N$, a direct projection of the equilibrium equations is possible via D’Alembert’s principle, which amounts to an integration over the equations of motion.

$\begin{array}{l}\begin{array}{c}\begin{array}{c}\underset{0}{\overset{L}{{\displaystyle \int }}}\left[\left(\mathit{F}\left(s,t\right)+{\mathit{F}}^{\mathit{*}}\left(s,t\right)\right)\cdot {\mathit{\gamma }}_j\left(s,t\right)+\left(\mathit{M}\left(s,t\right)+{\mathit{M}}^{\mathit{*}}\left(s,t\right)\right)\cdot {\mathit{\beta }}_j\left(s,t\right)\right] ds=Q_{j,nc}\end{array}\\ j=1,\dots ,N\end{array}\\ \begin{array}{cc}\mathit{F}\left(s,t\right)=-{\partial }_s\mathit{n}\left(s,t\right),& {\mathit{F}}^{\mathit{*}}\left(s,t\right)=\rho A{\partial }_{tt}\mathit{p}\left(s,t\right)+\rho I_k{\partial }_{tt}{\mathit{d}}_k\left(s,t\right)\\ \mathit{M}\left(s,t\right)=-{\partial }_s\mathit{m}+{\partial }_s\mathit{p}\times \mathit{n},\mathit{ }& {\mathit{M}}^{\mathit{*}}\left(s,t\right)=\rho I_k{\mathit{d}}_k\times {\partial }_{tt}\mathit{p}+{\partial }_t\left(\rho \mathit{J}\cdot \mathit{\omega }\right)\end{array}\\ Q_{j,nc}=\underset{0}{\overset{L}{{\displaystyle \int }}}\mathit{f}\left(s,t\right)\cdot {\mathit{\gamma }}_j\left(s,t\right)+\ell \left(s,t\right)\cdot {\mathit{\beta }}_j\left(s,t\right) \text{d}\tau .\end{array}(5)$

The velocity coefficient function and angular velocity coefficient function are defined as

${\mathit{\gamma }}_j\left(s,t\right)={\partial }_{q_j}\mathit{p}\left(s,t\right)={\partial }_{{\dot{q}}_j}{\partial }_t\mathit{p},\mathit{ }{\mathit{\beta }}_j\left(s,t\right)={\partial }_{{\dot{q}}_{\mathit{j}}}\mathit{\omega }\mathit{.}$

The velocity coefficients are the “Jacobian field” satisfying the relation Eq. 2.

Since the time derivatives of the momentum density and angular momentum density, ${\partial }_{tt}\mathit{p}$ and ${\partial }_t\left(\rho \mathit{J}\cdot \mathit{\omega }\right)$, can be written as linear functions of the ${\partial }_{tt}\mathit{q}$, the equations of motion are linear in the accelerations of the generalized coordinates, as expected. In the case of the PCC/PCS kinematic description, the derivatives ${\partial }_s\mathit{n}$ and ${\partial }_s\mathit{m}$, if resolved in $\mathrm{ℱ}\left(s,t\right)$, are zero except at the element boundaries. The equations may be integrated by parts into a form which shows the conjugacy of $\mathit{n}$ and $\mathit{v}$ and the conjugacy of $\mathit{m}$ and $\mathit{u}$.

$\begin{array}{l}\underset{0}{\overset{L}{{\displaystyle \int }}}{\partial }_s\mathit{n}\cdot {\mathit{\gamma }}_j ds={\left[\mathit{n}\cdot {\mathit{\gamma }}_j\right]}_0^L-\underset{0}{\overset{L}{{\displaystyle \int }}}\mathit{n}\cdot {\partial }_{q_j}\mathit{v} \text{d}s\\ \underset{0}{\overset{L}{{\displaystyle \int }}}{\partial }_s\mathit{m}\cdot {\mathit{\beta }}_j\text{d}s={\left[\mathit{m}\cdot {\mathit{\beta }}_j\right]}_0^L-\underset{0}{\overset{L}{{\displaystyle \int }}}\mathit{m}\cdot \left({\partial }_{q_j}\mathit{u}+\mathit{u}\times {\mathit{\beta }}_j\right)\text{d}s.\end{array}$

In the local frame, the equations take the following forms.

$\begin{array}{l}{\displaystyle \underset{0}{\overset{L}{\int }}{\partial }_s\mathit{n}\cdot {\mathit{\gamma }}_j ds}={\left[{}^{\mathbf{ℱ}\left(s,t\right)}\mathit{n}\cdot {}^{\mathbf{ℱ}\left(s,t\right)}{{\displaystyle \mathit{\gamma }}}_{\mathit{j}}\right]}_0^L-{\displaystyle \underset{0}{\overset{L}{\int }}{}^{\mathbf{ℱ}\left(s,t\right)}\mathit{n}\cdot \left({\partial }_{q_j}\left({}^{\mathbf{ℱ}\left(s,t\right)}\mathit{v}\right)+{}^{\mathbf{ℱ}\left(s,t\right)}{{\displaystyle \mathit{\beta }}}_j\times {}^{\mathbf{ℱ}\left(s,t\right)}\mathit{v}\right)}ds\\ {\displaystyle \underset{0}{\overset{L}{\int }}{\partial }_s\mathit{m}\cdot {\mathit{\beta }}_j\mathit{ds}}={\left[{}^{\mathrm{ℱ}\left(s,t\right)}\mathit{m}\cdot {}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{\beta }}}_j\right]}_0^L-{\displaystyle \underset{0}{\overset{L}{\int }}{}^{\mathrm{ℱ}\left(s,t\right)}\mathit{m}\cdot {\partial }_{q_j}\left({}^{\mathrm{ℱ}\left(s,t\right)}\mathit{u}\right)} ds.\end{array}$

Note also that if $\mathit{n}\left(s,t\right)$ is constant over $s$, the first integral is trivially zero.

Finally, note that if more than one rod-like body is present, then a sum over the bodies takes place in Eq. 5. Explicit constraints between the bodies may be handled via the method of Lagrange multipliers.

Learning-Based Approaches

Learning-based approaches, which are also sometimes referred to as “model-free” approaches, may be able to describe the relationships between the actuator inputs and observable outputs such as the end-effector motion without recourse to physical parameters and the laws of mechanics. These models usually serve a complementary purpose to those based on physical first principles. Since they require training data from a real robot or from another simulation model, they may be used for on-line control, inverse and forward kinematics, or for off-line analysis and testing of other algorithms such as for navigation and control. The a-priori prediction of behaviors from only design data is generally not possible to date using only learning-based methods.

A variety of purely kinematic approaches have been proposed. One learning approach uses an on-line estimation of the Jacobian matrix relating the time derivatives of the actuation variables ${\partial }_t\mathit{\tau }$ to the time derivatives ${\partial }_t\mathit{p}\left(L,t\right)$ and ${\partial }_t{\mathit{d}}_k\left(L,t\right)$, and it has been shown that this approach works for both position control and hybrid position/force control when appropriate sensing is available in hardware (Yip and Camarillo, 2014; Yip and Camarillo, 2016). Since the method requires no information about the robot or the environment a-priori, it enables control in complex scenarios, where highly complex physics-based models may have poorly observable parameters or states. It has also been shown that inverse kinematics for continuum robots may be approximated by a multilayer perceptron network (George et al., 2017; Grassmann et al., 2018; Lai et al., 2019), with multi-agent reinforcement learning (Ansari et al., 2016), with K-nearest neighbors and Gaussian mixture regression (Chen and Lau, 2016), and with deep reinforcement learning (Satheeshbabu et al., 2019). For reconfigurable robots subject to varying loads, it has been shown that classification of the load state using long short-term memory networks can substantially improve open-loop kinematic control (Nicolai et al., 2020). For flexible catheters, a combination of a support vector machine classifier and a fully connect neural network regressor were demonstrated achieving sub-millimeter trajectory tracking errors (Jolaei et al., 2020). Learning may also play an important role in proprioception for continuum robots with distributed sensing, where additional sensors beyond actuator-collocated ones are available for measuring the robot shape (Truby et al.,2020).

It has also been shown that dynamic models may be learned. Under a state observation of the form $\mathit{x}=F\left(\mathit{q}\right)$ provided by sensors, where it is presumed that the dimension of $\mathit{x}$ and $\mathit{q}$ are the same and that $F$ is invertible, the dynamics of the system can be posed as a one-to-one mapping $\left(\mathit{\tau },\mathit{x},{\partial }_t\mathit{x}\right)\to {\partial }_{tt}\mathit{x}$. On a real or simulated robot, this map representing the dynamics of the observables of the system can be approximated in discrete time via recurrent neural network (Thuruthel et al., 2017). Note that the kinematic relationships under static conditions are obviously also contained in this trained map for all points satisfying $\left(\mathit{\tau },\mathit{x},0\right)\to 0$. A similar approach using deep neural networks has also been demonstrated (Gillespie et al., 2016). Data-driven system identification based on the Koopman operator theory has led to control-oriented dynamic models amenable to model-predictive control (Bruder et al., 2020). Autoregressive with exogenous input (ARX) and nonlinear autoregressive with exogenous input (NARX) models have been studied for a single-section tendon-driven continuum robot, with the conclusion that NARX provides advantages in modeled end-effector position accuracy (Parvaresh and Moosavian, 2019). For closed-loop dynamic control, the direct reinforcement learning of a control policy which learns the mapping from current robot states, previous robot states, and desired end-effector position to the appropriate control action is possible (Thuruthel et al., 2019).

There are also learning-based approaches to control which do not explicitly construct kinematic or dynamic models. One such approach is based on direct learning from demonstration in the actuator space, which was successfully demonstrated on a tendon-driven continuum manipulator (Xu et al., 2016). Learning can form a part of a.

Actuator Models

Actuators in continuum and soft robots have been classified as either extrinsic, in which case the actuators are not a part of the deformable body, or intrinsic, in which case the actuators are an integral part of the deformable body. Examples of the former include tendons, the boundary conditions placed on concentric tube robots. Examples of the latter include soft pneumatic muscles (Walker et al., 2016).

The actuators may be modeled (very generally) as relationships between the actuation variables, generalized forces, and the dynamic state of the robot consisting of q and ${\partial }_t\mathit{q}$.

${\mathit{G}}_i\left({\tau }_i, \mathit{q},{\partial }_t\mathit{q},\mathit{ }{\mathit{Q}}_{nc}\right)=0.$

However, the nature of the model may change depending on the exact form of $G_i$. If $G_i$ involves only ${\tau }_i$ and $\mathit{q}$, then it is exactly in the form of a holonomic constraint. In general, it may not be simple to find a reduced set of independent coordinates satisfying the constraint, and a Lagrange multiplier technique may be required to enforce it. On the other hand, if $G_i$ can be inverted to find ${\mathit{Q}}_{nc}={\widehat{G}}_i\left({\tau }_i, \mathit{q}, {\partial }_t\mathit{q}\right)$, then the actuation can be directly coupled to the equations of motion. Which of these two views of actuation is the more natural one depends on the characteristics of the particular actuator(s) and sensor(s) chosen.

A first example is the model of a fiber-reinforced elastic actuator, in which $V$ is the enclosed fluid volume of the actuator, $\tau =P$ is the fluid pressure, and ${\mathit{J}}_V$ is the Jacobian matrix relating the changes in the generalized coordinates to the change in volume of the fluid (Sedal et al., 2021). Then, the principle of virtual work indicates that

$\delta W=\tau \delta V=\tau {\mathit{J}}_v\left(q\right)\delta \mathit{q}=Q_{i,nc}\delta q_i$

${\mathit{Q}}_{nc}={\mathit{J}}_V^T\left(q\right)\tau ,\mathit{ }{\mathit{J}}_{\mathit{V}}\left(q\right)={\partial }_{\mathit{q}}V.$

Another explicit example is found in the case of a tendon-driven robot. If enough support for the tendon is provided, a reasonable model for the points occupied by the tendon is a continuous curve described by ${\mathit{p}}_t\left(s,t\right)=\mathit{p}\left(s,t\right)+\mathit{a}\left(s\right)$ with $\mathit{a}\left(s\right)=a_1\left(s\right){\mathit{d}}_1\left(s,t\right)+a_2\left(s\right){\mathit{d}}_2\left(s,t\right)$ (Rucker and Webster, 2011). For the sake of simplicity, restrict the tendon to a planar path with $a_2\left(s\right)=0$. The tendon length can then be calculated as an integral functional involving the deformation gradient evaluated along the tendon path using the curve-based kinematic hypotheses:

$\begin{array}{c}{\ell }_t\left(\mathit{q}\right)={\displaystyle \underset{0}{\overset{L}{\int }}d\ell , d{\ell }^2}=d{s}^2\left({\mathit{d}}_3+{\partial }_s\mathit{a}\right){\mathit{F}}_t^T{\mathit{F}}_t\left({\mathit{d}}_3+{\partial }_s\mathit{a}\right), {\mathit{F}}_t=\frac{\partial {\mathit{p}}_t}{\partial {\mathit{P}}_t}. \end{array}(6)$

If the tendons are not fully constrained, other models for ${\ell }_t\left(\mathit{q}\right)$ may be more appropriate (Rao et al., 2021). What is noteworthy about either length formulation is that the nonconservative generalized forces do not naturally appear. If the tendon lengths are a known quantity, the actuator model is a holonomic constraint on the generalized coordinates. The problem can be treated via the method of Lagrange multipliers. The Lagrange multiplier will be exactly the tendon tension, and the principle of virtual work can be used to reveal the exact form of the terms in ${\mathit{Q}}_{nc}$ corresponding to the Lagrange multiplier.

$\delta W=\tau \delta {\ell }_t=\tau {\mathit{J}}_{\ell }\left(q\right) \delta \mathit{q}=Q_{j,nc}\delta q_j,{\mathit{J}}_{\ell }\left(q\right)={\partial }_{\mathit{q}}{\ell }_t\left(\mathit{q}\right).$

Therefore, the effect of the tendon alone (not considering any frictional forces) must be

${\mathit{Q}}_{nc}={\mathit{J}}_{\ell }^T\left(q\right)\tau .$

Note that the causal form in which the tendon tensions are known is “easier” to handle since no additional equations must be added. The causal form involving known tendon lengths requires the addition of the nonlinear length constraints Eq. 6 to the equation set and the tension becomes an algebraic unknown along with the accelerations, forming a nonlinear differential-algebraic system in the dynamic case or a nonlinear algebraic system in the quasistatic case. The need to solve a DAE system disappears if the tendon is considered a spring element, since then the force is determined as a function of the difference between ${\ell }_t\left(\mathit{q}\right)$ and the tendon displacement input $d$.

The resulting model form as a set of ordinary differential equations or differential algebraic equations is shown for a variety of common continuum robot actuators in Table 2.

TABLE 2

TABLE 2. Model form as an ODE or DAE system based on actuator type, assuming a single rod model architecture for the model.

Materials

The kinematic hypotheses and mechanics models must be augmented by constitutive laws (material models) to complete the model of a continuum robot. For quasistatic models, the choice is usually between linear elasticity and other hyperelastic material models. For dynamic models, an additional choice of damping or friction laws is generally required to produce realistic responses.

Linear Elasticity

In the case of quasistatic models, a common assumption in the literature has been to assume a Hookean (linear) material response. In this case, if one assumes that the backbone curve passes through the neutral axis of bending, the following constitutive laws apply:

${}^{\mathrm{ℱ}\left(s,t\right)}\mathit{m}={}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{K}}}_{bt}\left({}^{\mathrm{ℱ}\left(s,t\right)}\mathit{u}-{}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{u}}}_0\right){}^{\mathrm{ℱ}\left(s,t\right)}\mathit{n}={}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{K}}}_{\mathit{se}}\left({}^{\mathrm{ℱ}\left(s,t\right)}\mathit{v}-{}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{v}}}_0\right).$

The matrices ${\mathit{K}}_{bt}$ and ${\mathit{K}}_{se}$ may be calculated from the classical Euler-Bernoulli or Timoshenko beam theories and the entries are the flexural and torsional rigidities and shear and extension rigidities, respectively. The explicit relationships follow below (Antman, 2005).

$m_{\alpha }\left(s,t\right)=\left(E{J}_{\alpha \beta }\right)\left(s\right)\left[u_{\beta }\left(s,t\right)-u_{\beta 0}\left(s\right)\right], m_3=D\left(s\right)u_3\left(t\right).$

Note that bending about ${\mathit{d}}_1$ and ${\mathit{d}}_2$ are characterized by the second moments of area and the Young’s modulus $E$, while the torsional moment is related to the torsional strain variable by a coefficient $D$ solving the St. Venant torsion problem. Only in the case of isotropic rods with circular cross section is this equal to the usual shear modulus $G$ times the polar moment of area $J_{33}$. Formulas for a wide variety of cross sections that are uniform over $s$ have been tabulated (Roark et al., 2002). The Timoshenko model for shear and elongation adds the following relationships.

$n_{\alpha }=\left(GA\right)\left(s\right) v_{\alpha }, n_3=\left(EA\right)\left(s\right) \left[v_3-1\right].$

Hyperelastic Material Models

Many other hyperelastic models are possible choices, such as Yeoh, neo-Hookian, Gent, Ogden, and Mooney-Rivlin (He et al., 2018; Shiva et al., 2019; Zhang et al., 2019; Antonelli et al., 2020; Bacciocchi and Tarantino, 2021; Zhao et al., 2021). Although in general one may expect that these more complex material models should offer improved model accuracy, it has been shown recently that, at least for some robot designs, a linear stress-strain response may be more than adequate (Shiva et al., 2019). Any hyperelastic law can be represented within the Cosserat rod framework as a strain energy density function $W.$

$W=W\left({}^{\mathbf{ℱ}\left(s,t\right)}\mathit{u},{}^{\mathbf{ℱ}\left(s,t\right)}\mathit{v}\right)$

${}^{\mathbf{ℱ}\left(s,t\right)}\mathit{m}={\partial }_{\mathit{u}}W,\mathit{ }{}^{\mathbf{ℱ}\left(s,t\right)}\mathit{n}={\partial }_{\mathit{v}}W.$

The details of these calculations for each of the respective hyperelastic models is omitted for the sake of brevity and can be found in the cited references.

Damping and Friction

The introduction of dissipative mechanisms is generally necessary to encourage numerical stability in dynamic models and to produce realistic dynamic responses. Additionally, in some cases static friction plays a significant role in determining the quasistatic solutions, such as in tendon-driven catheters (Jung et al., 2014). Viscous damping may be introduced via the Kelvin-Voigt material model, which extends the linear elastic models to include rate-dependence in the stress-strain relationship (Gilbert and Godage, 2019; Mustaza et al., 2019).

In the curve-based framework, the Kelvin-Voigt law takes the following form (Linn et al., 2013):

$\begin{array}{l}{}^{\mathbf{ℱ}\left(s,t\right)}\mathit{m}={}^{\mathbf{ℱ}\left(s,t\right)}{{\displaystyle \mathit{K}}}_{bt}\left({}^{\mathbf{ℱ}\left(s,t\right)}\mathit{u}-{}^{\mathbf{ℱ}\left(s,t\right)}{{\displaystyle \mathit{u}}}_0\right)+{}^{\mathbf{ℱ}\left(s,t\right)}{{\displaystyle \mathit{B}}}_{bt}\mathit{ }{\partial }_t\left({}^{\mathbf{ℱ}\left(s,t\right)}\mathit{u}\right)\\ {}^{\mathbf{ℱ}\left(s,t\right)}\mathit{n}={}^{\mathbf{ℱ}\left(s,t\right)}{{\displaystyle \mathit{K}}}_{\mathit{se}}\left({}^{\mathbf{ℱ}\left(s,t\right)}\mathit{v}-{}^{\mathbf{ℱ}\left(s,t\right)}{{\displaystyle \mathit{v}}}_0\right)+{}^{\mathbf{ℱ}\left(s,t\right)}{{\displaystyle \mathit{B}}}_{\mathit{se}}\mathit{ }{\partial }_t\left({}^{\mathbf{ℱ}\left(s,t\right)}\mathit{v}\right).\end{array}$

The matrices $K_{bt}$ and $B_{bt}$ are related by time constants referred to as the extensional retardation time constant, ${\tau }_e={\eta }_E/E$, and the shear retardation time constant, ${\tau }_s=\eta /G$, with ${\eta }_E$ the “extensional viscosity” and $\eta$ the shear viscosity.

$\begin{array}{l}{}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{B}}}_{bt}={}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{K}}}_{bt}\cdot \text{diag}\left({\tau }_e, {\tau }_e, {\tau }_s\right)\\ {}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{B}}}_{\mathit{se}}={}^{\mathrm{ℱ}\left(s,t\right)}{{\displaystyle \mathit{K}}}_{\mathit{se}}\cdot \text{diag}\left({\tau }_s, {\tau }_s, {\tau }_e\right).\end{array}$

Static friction models have also been considered for concentric tube robots (Lock and Dupont, 2011), tendon-driven continuum robots (Li et al., 2020), and continuum robots having sheathed tendons or multiple actuated backbones (Roy et al., 2017).

Discussion

The wide variety of modeling choices that have been described offer the modeler an almost paralyzing array of choices. In the subsections that follow, several questions are posed. The available evidence from the literature as well as analyses guided by classical theories of mechanics are used to discuss these questions and to provide guidance during the initial stages of selecting modeling approaches.

Considerations for Kinematic Hypotheses

The literature on modeling of continuum and soft robots suggests that errors in kinematic models, quantified by the absolute tip positioning error as a percentage of the overall root length, are typically on the order of a few percent. Therefore, there may be little benefit to increasing the order of a spectral method or to further subdividing the domain in an element-based method once the absolute accuracy with respect to the true solution reaches this point. In the sections that follow, analysis and recommendations for kinematic hypotheses which are derived from consideration of the mechanics of bending are offered. Table 3 provides a summary of the recommendations in terms of increasing either the number of elements or the order of the interpolation (assuming that $\mathit{u}$ is the interpolated variable). Figure 5 depicts the decisions leading to the various types of models that have been discussed.

TABLE 3

TABLE 3. Summary of recommendations to increase either the number of elements or the order of interpolants based on model assumptions and robot-environment conditions.

FIGURE 5

FIGURE 5. Flowchart depicting the modeling decisions to be made when selecting a model type for a biomedical continuum robot.

Considerations for Cantilevered Concentrated Loadings

For continuum robots which are soft enough to exhibit substantial compliance to environmental loads (for example those that may be presented by contact with human anatomy), one of the first considerations for modeling should be consistency with the requirements for accurately modeling cantilevered, concentrated loads.

Let the Cosserat rod equations be recast in terms of the angle of the tangent vector and the load and deformation fixed to the plane defined by ${\mathit{d}}_3\left(0,t\right)$ and ${\mathit{d}}_1\left(0,t\right)$, let the boundary condition $\mathit{g}\left(0,t\right)$ be fixed, and let the load be concentrated at $s=L$ and modeled by $\mathit{F}=P\delta \left(s-L\right){\mathit{d}}_1\left(0,t\right)$ for the scalar force magnitude $P$. Furthermore, assume the material response is linear and that the robot is inextensible. These restrictions simplify the problem and result in the following dimensionless boundary value problem governing the angle $\theta$, which represents the tangent vector:

${\mathit{d}}_3\left(s,t\right)=\cos \left(\theta \right){\mathit{d}}_3\left(0,t\right)+\sin \left(\theta \right){\mathit{d}}_1\left(0,t\right)$

${\partial }_{\xi \xi }\theta +\lambda \cos \left(\theta \right)=0, \xi =\frac{x}{L}, \lambda =\frac{P{L}^2}{EI}\theta \left(0,t\right)=0, {\partial }_{\xi }\theta \left(1,t\right)=0.$

The boundary value problem has a known solution:

$\theta \left(\xi ,t\right)=2{\sin }^{-1}\left(k\text{sn}\left(K\left(k\right)-\left(1-\xi \right)\sqrt{\lambda };k^2\right)\right)-\frac{\pi }{2}.$

The quantity $k$ is a scalar that may be found by Newton-Raphson iteration on the following equation, which is implied by the boundary condition $\theta \left(0\right)=0$

$k\text{sn}\left(K\left(k\right)-\sqrt{\lambda };k^2\right)=\sqrt{2}/2.$

The PRB models have the attractive property that they map the problem back into the domain of traditional robotic manipulators, with the obvious advantage that all the tools and knowledge that have been developed in that context (in general, restricted to underactuated mechanisms) now apply to the continuum robot. In the traditional PRB models, the inertia properties are lumped into the links formed by the model, and the stiffness and damping properties are lumped into the joints between links. This lumping introduces error, but it has been shown that optimization of the parameters of the rigid body model can lead to accurate mechanical responses for both cantilevered transverse loads and for applied or internal moments (Chen et al., 2011). Given that the optimal 3R planar PRB model has three degrees of freedom, it is a fair comparison to place the model in competition with other three-DOF models. Here we consider the following three sets of potential kinematic hypotheses and matching constitutive laws and compare them with each other and with the exact solution. Without loss of generality, let $L=1$ and $EI=1$.

Piecewise constant curvature:

${\mathit{u}}_{\text{PCC}}={\displaystyle \sum _{e=1}^3{u}_e {\text{χ}}_{{\text{Γ}}_e}\left(s\right) {\mathit{d}}_2, {\chi }_{{\text{Γ}}_e}\left(s\right)}=\{\begin{array}{cc}1& s\in {\text{Γ}}_e\\ 0& s\notin {\text{Γ}}_e\end{array}, D=\left[0,\frac{1}{3}\right]\cup \left[\frac{1}{3},\frac{2}{3}\right]\cup \left[\frac{2}{3},1\right]m\left(s\right)=EIu\left(s\right).$

PRB:

${\mathit{u}}_{\text{PRB}}={\displaystyle \sum _{i=1}^3{q}_i\delta \left(s-s_i\right) {\mathit{d}}_2}, s_1=0.125, s_2=0.475, s_3=0.863, D=\left[\mathrm{0,1}\right]\frac{m\left(s_1\right)}{q_1}=3.25, \frac{m\left(s_2\right)}{q_2}=2.84, \frac{m\left(s_3\right)}{q_3}=\mathrm{2.95.}$

Spectral:

${\mathit{u}}_{\text{S}}={\displaystyle \sum _{i=1}^3{q}_i{s}^{i-1} {\mathit{d}}_2}, D=\left[\mathrm{0,1}\right].$

$m\left(s\right)=EI u\left(s\right)$

How well does each of the strategies perform when given 3-DOF to capture the deformation? The answer is depicted in Figure 6, showing tip position error in percent of robot length versus the dimensionless cantilevered load index $\lambda$. For cantilevered loads, a single spectral element which is quadratic in $\mathit{u}$ is a far better choice than either a 3-element PCC model or a 3R PRB model. If nonzero shear forces are expected, the PCC model seems to have little in its favor; the Jacobian for the PRB model is simpler to calculate, meaning that the statics equations in Eq. 5 are easier to formulate, and the tip position is predicted more accurately, which also implies that the overall structural stiffness is more accurate for the PRB model than for the PCC model. The Jacobian for the spectral model, unlike the other two, cannot be obtained in an exact closed form.

FIGURE 6

FIGURE 6. Error in reproducing the correct behavior under cantilevered loading conditions for three-DOF kinematic hypotheses of the PCC, PRB, and spectral types.

The results imply that the typical piecewise constant-curvature assumption used in the development of geometrically nonlinear models for robots is a poor choice from the perspective of mechanics whenever a concentrated external load is present and is expected to produce internal shear forces which are transverse to $\mathit{v}$. In summary, if point loads are present on the robot, a linear interpolant of internal moment (equivalently curvature) is necessary to accurately capture the static equilibrium configurations for unrefined elements even in the small deflection case, and degrees of freedom are better spent on increasing the order of the interpolants than on increasing the number of elements.

Considerations for Non-Uniform Flexural Rigidity

In the case of non-uniform flexural rigidity, element refinement is more effective at reducing approximation error than increases in order. This conclusion is easily justified by the observation that if $K\left(s\right)$ is a linear function, say for example $1+as$, then in the simplest planar case with a constant internal moment one would be tasked to find another polynomial function $\kappa \left(s\right)$ such that $\kappa \left(s\right)K\left(s\right)=C$ for some constant $C$. But this is obviously impossible, because $\kappa =C/K$ is a rational function, not a polynomial, and the Maclaurin series at $s=0$ has a finite region of convergence. In the example case, the expansion is $C\cdot {\displaystyle \sum {\left(-1\right)}^i{a}^i{s}^i}$. The series does not converge unless $\left|s\right|<\left|1/a\right|$ and as $s$ approaches this upper bound, the number of terms in the series required to obtain convergence to a fixed tolerance increases without bound. Element refinement, on the other hand, has exactly the effect of reducing $\left|s\right|$, ensuring convergence. For this reason, single-element, spectral-type methods are not recommended as a first choice if non-uniform flexural rigidity is present.

Considerations for Uniformly Distributed Loads

Distributed loads may act on biomedical continuum robots. The most obvious of these loads is a gravitational force distributed along the length of the robot. Other common forces may include buoyancy forces, electric forces, magnetic forces, and aerodynamic and hydrodynamic forces. The simplest possible model of a distributed load is a uniform one that is applied normal to the body of a robot which is initially in a straight configuration. In this case, the solution to the linearized Euler-Bernoulli equation is in general a fourth-order polynomial in position. The shear force is a linear function of arc length and the internal moment (and hence curvature in the linear elastic case) is quadratic. If the shape is discretized at the level of angle, the discretization should be cubic to accommodate a uniform load.

Considerations for Elastic Environmental Contact

For continuum robots in contact with soft bodies such as the soft tissues of the human anatomy, the contact might be well-described, at least in the region of contact, by a model like the linear elastic foundation model. For small deflections, the linearized Euler-Bernoulli model with a linear elastic foundation is modeled by the following differential equation (see Figure 7A).

$\left(EIy{,}_{ss}\right){,}_{ss}=-ky.$

FIGURE 7

FIGURE 7. (A) Schematic diagram for the beam on an elastic foundation as a model for a continuum robot interacting with soft tissue. (B) Example with $\beta =3$, showing the shape of the displacement that must be approximated.

The homogeneous solution to this equation has the following form.

$y_h\left(s\right)=\exp \left(-\beta s\right)\left[C_1\sin \left(\beta s\right)+C_2\cos \left(\beta s\right)\right]+\exp \left(\beta s\right)\left[C_3\sin \left(\beta s\right)+C_4\cos \left(\beta s\right)\right].$

The constant $\beta$ depends only on the properties of the beam and the foundation.

${\beta }^2=\sqrt{\frac{k}{4EI}}.$

To what degree of accuracy does a polynomial shape function (assuming the small-deflection case) approximate $y_h$?

To answer this question, one should find the best approximation of $y_h$ under a particular norm on $L^2\left[0,\ell \right]$. Here we select the 2-norm and study the approximation error for 3^rd, 4^th, and 5^th order polynomials. Since $\beta$ is related to the ratio of stiffnesses $k$ and $EI$, and has dimension Length⁻¹, we restrict the range of the dimensionless group $\ell \beta$ to $\left(0.1, 10\right)$. This range is consistent with the idea of compliance matching as a form of embodied intelligence in biomedical continuum robots. Note that as $\ell \beta \to 0$, the solution $x_h$ approaches a constant, which is easy to interpolate. As $\ell \beta \to \infty$, the elastic foundation is becoming infinitely stiff relative to the body of the robot, modeling a hard contact. In this case, the internal forces and moments and the resulting deformations become strongly localized, and a point load may be a more appropriate model for the contact than an elastic foundation.

The physical solutions to the equation decay away from the application of a point load. Therefore, we restrict the approximation problem to the consideration of the two functional forms that follow on a domain $\left[\mathrm{0,1}\right]$ for $\beta \in \left(0.1,10\right)$.

$x_1=\exp \left(-\beta s\right)\sin \left(\beta s\right)x_2=\exp \left(-\beta s\right)\cos \left(\beta s\right).$

See Figure 7B for examples with $\beta =3$. Errors for polynomials $p_1$ and $p_2$ approximating $x_1$ and $x_2$, respectively, are shown in Figure 4 as the maximum absolute errors.

$e_i=\frac{\underset{s}{\max }\left|p_i\left(s\right)-x_i\left(s\right)\right|}{\underset{s}{\max }\left|x_i\left(s\right)\right|}$

For a single 5^th order polynomial in shape, the maximum absolute error in approximating either $x_1$ or $x_2$ remains below 1% if $\ell \beta <4.3$ (Figure 8). As $\ell \beta$ increases beyond this value, the polynomial approximations to $x_1$ and $x_2$ begin to oscillate with increasing maximum error.

FIGURE 8

FIGURE 8. Approximation errors for best polynomial fits in the L² norm to the solution for the linear beam on an elastic foundation problem. Higher-order polynomials permit greater elastic foundation stiffnesses.

To put this in a practical perspective, a typical colonoscope has a linearized flexural rigidity of $EI\approx 0.02 N m^2$ (Wehrmeyer et al., 1998). Soft tissues may have an elastic foundation stiffness of approximately 4 kN/m² (Asadian et al., 2011). This results in $\beta =15$ and therefore a hypothesis which is 5^th order in position (3^rd order in strain variables) should not have elements longer than approximately 0.28 m. Note that for a hypothesis that is linear in $s$ for the strain variables, the length requirement would drop to approximately 96 mm, and for PCC elements, the length would drop to only 38 mm. For a spatial robot model that is inextensible and un-shearable and is 1 m long, this would result in a PCC model with approximately 81 degrees of freedom (27 elements at 3 DOF/element), a linear strain variable model with approximately 66 degrees of freedom (11 elements at 6 DOF/element), or a cubic strain model with approximately 48 degrees of freedom (4 elements at 12 DOF/element). Therefore, if environmental contacts are soft and distributed over a long length, there is a strong incentive to develop models with higher-order strain variable hypotheses.

Considerations for Numerical Methods

Solution Multiplicity

In general, the problem defined by Eq. 5 together with any constraints is a nonlinear algebraic problem, even if linear material models are used. This is either a consequence of the nonlinear geometry, which shows up in any finite-strain relationship between the strain variables and the position and orientation of the body, or a consequence of nonlinear material behavior, or both. In special cases, the problem may become linear; for example, if the actuators and generalized forces are related linearly, linear constitutive laws are used, and no external loads are present. For nonlinear static problems, the Newton-Raphson method and trust-region methods like the Levenberg-Marquardt method generally work well, but the modeler must be cautious of the possibility of solution multiplicity.

In other words, a function $\mathit{q}=f\left(\mathit{\tau }\right)$ does not always exist because there may be two or more values of $\mathit{q}$ which satisfy the equilibrium conditions given $\mathit{\tau }$. This solution multiplicity is accompanied by a singular tangent stiffness matrix for some value of $\mathit{q}$ and possibly a tangent stiffness matrix with negative eigenvalues, as is the case for so-called “negative-stiffness mechanisms” (Platus, 1999). The coupling between kinematics and mechanics means that it is not always safe to assume the existence of a “forward kinematic mapping” which computes the C-space coordinates from the actuator variables and then the task-space variables from the C-space coordinates. Consider the case in which Eq. 5 is of the form $\mathit{F}\left(\mathit{q},\mathit{\tau },{\mathit{Q}}_{nc}^{\ast }\right)=0$ where ${\mathit{Q}}_{nc}^{\ast }$ includes only those generalized forces which are not algebraically related with $\mathit{q}$ and $\mathit{\tau }$. Then a perturbation analysis yields the C-space Jacobians with respect to $\mathit{\tau }$ and ${\mathit{Q}}_{nc}^{\ast }$

$\delta F={\partial }_{\mathit{q}}\mathit{F }\delta \mathit{q}+{\partial }_{\mathit{\tau }}\mathit{F }\delta \mathit{\tau }+{\partial }_{{\mathit{Q}}_{nc}}\mathit{F }\delta {\mathit{Q}}_{nc}^{\ast }=0.$

$\delta \mathit{q}=-{\left({\partial }_{\mathit{q}}\mathit{F}\right)}^{-1}\left({\partial }_{\mathit{\tau }}\mathit{F }\delta \mathit{\tau }+{\partial }_{{\mathit{Q}}_{nc}}\mathit{F }\delta {\mathit{Q}}_{nc}^{\ast }\right).$

It is evidently at configurations with singular ${\partial }_{\mathit{q}}\mathit{F}$ where multiple solutions may arise. This is one reason that quasistatic “resolved-rate” or continuation-type methods may fail to converge; dynamic models do not suffer this problem since the accelerations are resolved.

Time Stepping

For time stepping, explicit ode integrators can become prohibitively computationally expensive. This is a consequence of the fact that unresolved vibrational modes (as defined for linear test problems) become unstable using explicit methods. Implicit integrators and those designed for solving stiff ODEs and DAEs, such as the trapezoidal method or the backwards difference formulae, are preferable. Energy-preserving integrators have the benefit that the damping behavior is caused entirely by the material model, ensuring repeatable dynamic behavior with different time steps.

Current and Future Challenges in Modeling

Generalizability and Re-Usability

Despite the growing body of evidence that models built on the foundation of the Cosserat rod equations are an adequate description of many continuum robots, one challenge that still faces practitioners is a lack of standardized tools to build new model simulation codes. For rigid robots, a wide variety of domain-specific modeling languages are available and permit concise descriptions within an easy-to-use interface to build new models. One example of this is the Universal Robot Description Format and Gazebo simulator within the Robot Operating System, but there are many others presently available including Simulink/Simscape, Dymola, and other Modelica-language based toolsets such as OpenModelica (Brück et al., 2002; Fritzson et al., 2006; Miller and Wendlandt, 2011; Sucan and Kay, 2019). To enable the widespread re-use of validated modeling components, a library of reusable “model building blocks” for continuum robots should be designed. Some important capabilities of such a library would be the following:

• Coupling of curve-based models to rigid multibody models.

• Coupling of curve-based models and general finite element models.

• Incorporation of common actuator models.

• Incorporation of common constraints (length, concentricity, no-penetration, selective inextensibility/strong anisotropy, revolute and prismatic joints, etc.)

• User-selected switching between dynamic and quasi-static model generation.

For biomedical continuum robots in particular, models which couple to mechanical models of human anatomy are needed. Coupling of state-of-the-art models for continuum robots or their direct incorporation with real-time finite element codes using GPU acceleration is a promising approach (Allard et al., 2007; Duriez and Bieze, 2017).

Novel Kinematic Hypotheses

There is a great deal of freedom in element-based kinematic hypotheses which has yet to be explored. One interesting avenue is the use of a shared or constrained DOF between elements. The motivation for this idea is that for dynamic models, time stepping is sometimes restricted or difficult for “stiff” problems having many eigenvalues. The equations of motion for solid continua are wave equations, which means that if many elements are stacked end-to-end, acoustic waves (axial compression and tension) and twist waves (torsional waves) through the structure may be resolved by the model. For most robotics applications, these modes are likely to be irrelevant, and constraining the problem so that they do not exist in the model may improve computational performance. The elimination of twist waves in elastic rod models was previously considered by an energy minimization argument (Bergou et al., 2008).

Furthermore, adaptive kinematic hypotheses based on pre-defined, switchable degrees of freedom that permit local, automatic refinement of the model may allow greatly improved computational efficiency in problems involving a-priori unknown environmental interactions or constraints. This will permit, for example, a single high-order element to describe the deformation in free-space, while local refinement can take place where a catheter contacts a vessel wall, a robotic endoscopic system contacts the colon, or where multi-fingered hands contact an object to manipulate it.

Learning

Within the context of continuum and soft robotics, data-driven methods have begun to demonstrate strong utility. For example, Long Short Term Memory networks can capture hysteresis in pneumatically actuated catheters (Wu et al., 2021), and offline simulation of first-principles models can be used to learn reduced-order models using the snapshot-based proper orthogonal decomposition, resulting in new models suitable for real-time control and other applications requiring fast computation (Goury and Duriez, 2018; Katzschmann et al., 2019b). The continued development of learning methods enabling low-DoF representations will be an important future area of research.

There are also interesting opportunities for learning that amalgamate first-principles models with data-driven model “correctors,” or which use constrained learning techniques to identify models which are topologically like a curve-based model. One possibility is to use a low-DOF curve-based model capturing some of the behavior and to introduce a nonconservative generalized force ${\mathit{Q}}_{nc}$ which is learned from observed data to close the gap between simulation and reality. Learning-based methods which are constrained to obey fundamental principles are another emerging area of research, such as learning the Lagrangian or Hamiltonian function of systems directly from data (Ahmadi et al., 2018; Lutter et al., 2018).

Dynamic Model Validation

Although many dynamic models have been proposed, the validation of these models is currently lacking. There are many opportunities for rigorous evaluation and comparison of models with experimentally obtained data. The best and strongest form of model validation would be to instrument real robots with enough sensors to measure all the quantities appearing in Eq. 5 or the equivalent formulations for PRB and general continuum models, and to calculate the model residuals over conditions ranging over static, low-acceleration, and high-acceleration (e.g. sudden contact) regimes. This is clearly a challenging experimental task that may require state reconstruction and many sensors just to measure the configuration trajectory $\mathit{q}\left(t\right)$. Other options for validation may include comparison of standard test signal response characteristics (e.g. rise time, percent overshoot, settling time, steady-state error, and oscillation period) in response to both actuator inputs and environmental perturbations.

There are also many other interesting questions that can be asked and answered which are quantitative in a different sense, but which may be even more aligned with the spirit of soft and continuum robotics theory. For example, a model and simulated controller could be used to predict the success or failure of the navigation of a robotic catheter through tortuous vasculature parameterized by some measure of “tortuosity,” and then the classification error could be assessed via experiment matching the simulations.

Conclusion

Continuum robots offer solutions to problems in biomedical applications which may not be solvable by traditional robotics technologies. With these new robots came the need for new models. A wide variety of physics-based and learning-based approaches to the modeling of continuum manipulators—both those made of hard materials and soft materials—are now available to the roboticist who needs them. This can lead to a dizzying array of choices for the uninitiated. This manuscript has reviewed the state-of-the-art approaches using a common language, discussed considerations which can guide the modeler when selecting which methods to use and some numerical difficulties to be aware of, and offered a view of the current and future challenges in the modeling of continuum robots. As modeling techniques continue to improve in terms of predictive power, as techniques begin to standardize, and as system identification techniques for soft and continuum robots mature, there is every reason to expect that the field will continue to expand, find new applications, and ultimately lead to transformative robotic solutions for human problems.

Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author Contributions

HG is solely responsible for all aspects of this manuscript.

Funding

This material is based upon work supported by the National Science foundation under Grant No. 2024795. Any opinions, findings, and conclusion or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

Conflict of Interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

¹There are additional considerations for this placement in the case of dynamical models, discussed below.

References

Ahmadi, M., Topcu, U., and Rowley, C. (2018). “Control-Oriented Learning of Lagrangian and Hamiltonian Systems,” in 2018 Annual American Control Conference (ACC), Milwaukee, WI, June 27–29, 2018, 520–525. doi:10.23919/ACC.2018.8431726

CrossRef Full Text | Google Scholar

Allard, J., Cotin, S., Faure, F., Bensoussan, P.-J., Poyer, F., Duriez, C., et al. (2007). “SOFA - an Open Source Framework for Medical Simulation,” in MMVR 15 - Medicine Meets Virtual Reality, Palm Beach, 13–18.

Google Scholar

Ansari, Y., Falotico, E., Mollard, Y., Busch, B., Cianchetti, M., and Laschi, C. (2016). “A Multiagent Reinforcement Learning Approach for Inverse Kinematics of High Dimensional Manipulators with Precision Positioning,” in 2016 6th IEEE International Conference on Biomedical Robotics and Biomechatronics (BioRob), Singapore, June 26–29, 2016 (IEEE), 457–463. doi:10.1109/BIOROB.2016.7523669

CrossRef Full Text | Google Scholar

Antman, S. S. (2005). Nonlinear Problems of Elasticity. 2nd ed. Applied Mathematical Sciences, (Springer).

Google Scholar

Antonelli, M. G., Beomonte Zobel, P., D’Ambrogio, W., and Durante, F. (2020). Design Methodology for a Novel Bending Pneumatic Soft Actuator for Kinematically Mirroring the Shape of Objects. Actuators 9 (4), 113. doi:10.3390/act9040113

CrossRef Full Text | Google Scholar

Asadian, A., Kermani, M. R., and Patel, R. V. (2011). “An Analytical Model for Deflection of Flexible Needles during Needle Insertion,” in 2011 IEEE/RSJ International Conference on Intelligent Robots and Systems, San Francisco, CA, September 25–30, 2011, 2551–2556. doi:10.1109/IROS.2011.6094959

CrossRef Full Text | Google Scholar

Bacciocchi, M., and Tarantino, A. M. (2021). Finite Bending of Hyperelastic Beams with Transverse Isotropy Generated by Longitudinal Porosity. Eur. J. Mech. - A/Solids 85 (1), 104131. doi:10.1016/j.euromechsol.2020.104131

CrossRef Full Text | Google Scholar

Back, J., Manwell, T., Karim, R., Rhode, K., Althoefer, K., and Liu, H. (2015). “Catheter Contact Force Estimation from Shape Detection Using a Real-Time Cosserat Rod Model,” in 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Hamburg, Germany, September 28–October 2, 2015 (IEEE), 2037–2042. doi:10.1109/IROS.2015.7353647

CrossRef Full Text | Google Scholar

Bergou, M., Wardetzky, M., Robinson, S., Audoly, B., and Grinspun, E. (2008). Discrete Elastic Rods. ACM Trans. Graph. 27, 1–12. doi:10.1145/1360612.1360662

CrossRef Full Text | Google Scholar

Bieze, T. M., Largilliere, F., Kruszewski, A., Zhang, Z., Merzouki, R., and Duriez, C. (2018). Finite Element Method-Based Kinematics and Closed-Loop Control of Soft, Continuum Manipulators. Soft Robotics 5 (3), 348–364. doi:10.1089/soro.2017.0079

PubMed Abstract | CrossRef Full Text | Google Scholar

Bishop, R. L. (1975). There Is More Than One Way to Frame a Curve. The Am. Math. Monthly 82 (3), 246–251. doi:10.1080/00029890.1975.1199380710.2307/2319846

CrossRef Full Text | Google Scholar

Black, C. B., Till, J., and Rucker, D. C. (2018). Parallel Continuum Robots: Modeling, Analysis, and Actuation-Based Force Sensing. IEEE Trans. Robot. 34 (1), 29–47. doi:10.1109/TRO.2017.2753829

CrossRef Full Text | Google Scholar

Boyer, F., Ali, S., and Porez, M. (2012). Macrocontinuous Dynamics for Hyperredundant Robots: Application to Kinematic Locomotion Bioinspired by Elongated Body Animals. IEEE Trans. Robot. 28 (2), 303–317. doi:10.1109/TRO.2011.2171616

CrossRef Full Text | Google Scholar

Boyer, F., Lebastard, V., Candelier, F., and Renda, F. (2021). Dynamics of Continuum and Soft Robots: A Strain Parameterization Based Approach. IEEE Trans. Robot. 37, 847–863. doi:10.1109/TRO.2020.3036618

CrossRef Full Text | Google Scholar

Bretl, T., and McCarthy, Z. (2014). Quasi-Static Manipulation of a Kirchhoff Elastic Rod Based on a Geometric Analysis of Equilibrium Configurations. Int. J. Robotics Res. 33 (1), 48–68. doi:10.1177/0278364912473169

CrossRef Full Text | Google Scholar

Brockett, R. W. (1984). “Robotic Manipulators and the Product of Exponentials Formula,” in Mathematical Theory of Networks and Systems. Editor P. A. Fuhrmann (Berlin, Heidelberg: Springer), 120–129. Lecture Notes in Control and Information Sciences. doi:10.1007/BFb0031048

CrossRef Full Text | Google Scholar

Brück, D., Elmqvist, H., Olsson, H., and Mattsson, S. E. (2002). “Dymola for Multi-Engineering Modeling and Simulation,” in 2nd International Modelica Conference, Oberpfaffenhofen, Germany, March 18–19, 2002, 55-1–55–58.

Google Scholar

Bruder, D., Fu, X., Gillespie, R. B., Remy, C. D., and Vasudevan, R. (2021). Data-Driven Control of Soft Robots Using Koopman Operator Theory. IEEE Trans. Robot. 37, 948–961. doi:10.1109/TRO.2020.3038693

CrossRef Full Text | Google Scholar

Burgner-Kahrs, J., RuckerRucker, D. C., and Choset, H. (2015). Continuum Robots for Medical Applications: A Survey. IEEE Trans. Robot. 31 (6), 1261–1280. doi:10.1109/tro.2015.2489500

CrossRef Full Text | Google Scholar

Camarillo, D. B., Milne, C. F., Carlson, C. R., Zinn, M. R., and Salisbury, J. K. (2008). Mechanics Modeling of Tendon-Driven Continuum Manipulators. IEEE Trans. Robot. 24 (6), 1262–1273. doi:10.1109/TRO.2008.2002311

CrossRef Full Text | Google Scholar

Chen, G., Xiong, B., and Huang, X. (2011). Finding the Optimal Characteristic Parameters for 3R Pseudo-rigid-body Model Using an Improved Particle Swarm Optimizer. Precision Eng. 35 (3), 505–511. doi:10.1016/j.precisioneng.2011.02.006

CrossRef Full Text | Google Scholar

Chen, J., and Lau, H. Y. K. (2016). “Learning the Inverse Kinematics of Tendon-Driven Soft Manipulators with K-Nearest Neighbors Regression and Gaussian Mixture Regression,” in 2016 2nd International Conference on Control, Automation and Robotics (ICCAR) (Hong Kong: IEEE), 103–107. doi:10.1109/ICCAR.2016.7486707

CrossRef Full Text | Google Scholar

Chen, Y., Wang, L., Galloway, K., Godage, I., Simaan, N., and Barth, E. (2021). Modal-Based Kinematics and Contact Detection of Soft Robots. Soft Robotics 8, 298–309. doi:10.1089/soro.2019.0095

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, H., and Gupta, K. C. (1989). An Historical Note on Finite Rotations. J. Appl. Mech. 56 (1), 139–145. doi:10.1115/1.3176034

CrossRef Full Text | Google Scholar

Chirikjian, G. S., and Burdick, J. W. (1994). A Modal Approach to Hyper-Redundant Manipulator Kinematics. IEEE Trans. Robot. Automat. 10 (3), 343–354. doi:10.1109/70.294209

CrossRef Full Text | Google Scholar

Chirikjian, G. S., and Burdick, J. W. (1991). “Hyper-Redundant Robot Mechanisms and Their Applications,” in IEEE/RSJ International Workshop on Intelligent Robots and Systems ’91, Osaka, Japan, November 3–5, 1991, 185–190. doi:10.1109/IROS.1991.174447

CrossRef Full Text | Google Scholar

Cianchetti, M., Laschi, C., Menciassi, A., and Dario, P. (2018). Biomedical Applications of Soft Robotics. Nat. Rev. Mater. 3 (6), 143–153. doi:10.1038/s41578-018-0022-y

CrossRef Full Text | Google Scholar

Cianchetti, M., Ranzani, T., Gerboni, G., Nanayakkara, T., Althoefer, K., Dasgupta, P., et al. (2014). Soft Robotics Technologies to Address Shortcomings in Today's Minimally Invasive Surgery: The STIFF-FLOP Approach. Soft Robotics 1 (2), 122–131. doi:10.1089/soro.2014.0001

CrossRef Full Text | Google Scholar

Crisfield, M. A., and Jelenić, G. (1999). Objectivity of Strain Measures in the Geometrically Exact Three-Dimensional Beam Theory and its Finite-Element Implementation. Proc. R. Soc. Lond. A. 455 (1983), 1125–1147. doi:10.1098/rspa.1999.0352

CrossRef Full Text | Google Scholar

Della Santina, C., Katzschmann, R. K., Biechi, A., and Rus, D. (2018). “Dynamic Control of Soft Robots Interacting with the Environment,” in 2018 IEEE International Conference on Soft Robotics (RoboSoft), Livorno, Italy, April 24–28, 2018 (IEEE), 46–53. doi:10.1109/ROBOSOFT.2018.8404895

CrossRef Full Text | Google Scholar

Denavit, J., and Hartenberg, R. S. (1955). A Kinematic Notation for Lower-Pair Mechanisms Based on Matrices. J. Appl. Mech. 22 (2), 215–221. doi:10.1115/1.4011045

CrossRef Full Text | Google Scholar

Ding, J., Goldman, R. E., Xu, K., Allen, P. K., Fowler, D. L., and Simaan, N. (2013). Design and Coordination Kinematics of an Insertable Robotic Effectors Platform for Single-Port Access Surgery. Ieee/asme Trans. Mechatron. 18 (5), 1612–1624. doi:10.1109/TMECH.2012.2209671

PubMed Abstract | CrossRef Full Text | Google Scholar

Dupont, P. E., Lock, J., Itkowitz, B., and Butler, E. (2010). Design and Control of Concentric-Tube Robots. IEEE Trans. Robot. 26 (2), 209–225. doi:10.1109/TRO.2009.2035740

PubMed Abstract | CrossRef Full Text | Google Scholar

Duriez, C., and Bieze, T. (2017). “Soft Robot Modeling, Simulation and Control in Real-Time,” in Soft Robotics: Trends, Applications and Challenges (Springer), 103–109. doi:10.1007/978-3-319-46460-2_13

CrossRef Full Text | Google Scholar

Fraś, J., Jan, C., Maciaś, M., and Główka, J. (2014). “Static Modeling of Multisection Soft Continuum Manipulator for Stiff-Flop Project,” in In Recent Advances in Automation, Robotics and Measuring Techniques. Editors R. Szewczyk, C. Zieliński, and M. Kaliczyńska (Cham: Springer International Publishing), 365–375.

Google Scholar

Fritzson, P., Aronsson, P., Pop, A., Lundvall, H., Nystrom, K., Saldamli, L., Broman, D., and Sandholm, A. (2006). “OpenModelica - A Free Open-Source Environment for System Modeling, Simulation, and Teaching,” in 2006 IEEE Conference on Computer Aided Control System Design, 2006 IEEE International Conference on Control Applications, Munich, Germany, October 4–6, 2006, 1588–1595. doi:10.1109/CACSD-CCA-ISIC.2006.4776878

CrossRef Full Text | Google Scholar

Ganji, Y., and Janabi-Sharifi, F. (2009). Catheter Kinematics for Intracardiac Navigation. IEEE Trans. Biomed. Eng. 56 (3), 621–632. doi:10.1109/TBME.2009.2013134

PubMed Abstract | CrossRef Full Text | Google Scholar

George Thuruthel, T., Falotico, E., Manti, M., Pratesi, A., Cianchetti, M., and Laschi, C. (2017). Learning Closed Loop Kinematic Controllers for Continuum Manipulators in Unstructured Environments. Soft Robotics 4 (3), 285–296. doi:10.1089/soro.2016.0051

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilbert, H. B., and Godage, I. S. (2019). “Validation of an Extensible Rod Model for Soft Continuum Manipulators,” in 2019 2nd IEEE International Conference on Soft Robotics (RoboSoft), Seoul, South Korea, April 14–18, 2019 (IEEE), 711–716. doi:10.1109/ROBOSOFT.2019.8722721

CrossRef Full Text | Google Scholar

Gilbert, H. B., Hendrick, R. J., and Webster, R. J. (2016). Elastic Stability of Concentric Tube Robots: A Stability Measure and Design Test. IEEE Trans. Robot. 32 (1), 20–35. doi:10.1109/TRO.2015.2500422

PubMed Abstract | CrossRef Full Text | Google Scholar

Gillespie, M. T., Best, C. M., and Killpack, M. D. (2016). “Simultaneous Position and Stiffness Control for an Inflatable Soft Robot,” in 2016 IEEE International Conference on Robotics and Automation (ICRA), Stockholm, Sweden, May 16–21, 2016 (IEEE), 1095–1101. doi:10.1109/ICRA.2016.7487240

CrossRef Full Text | Google Scholar

Godage, I. S., Medrano-Cerda, G. A., Branson, D. T., Guglielmino, E., and Caldwell, D. G. (2015). Modal Kinematics for Multisection Continuum Arms. Bioinspir. Biomim. 10 (3), 035002. doi:10.1088/1748-3190/10/3/035002

PubMed Abstract | CrossRef Full Text | Google Scholar

Goury, O., Carrez, B., and Duriez, C. (2021). Real-Time Simulation for Control of Soft Robots with Self-Collisions Using Model Order Reduction for Contact Forces. IEEE Robot. Autom. Lett. 6 (2), 3752–3759. doi:10.1109/LRA.2021.3064247

CrossRef Full Text | Google Scholar

Goury, O., and Duriez, C. (2018). Fast, Generic, and Reliable Control and Simulation of Soft Robots Using Model Order Reduction. IEEE Trans. Robot. 34 (6), 1565–1576. doi:10.1109/tro.2018.2861900

CrossRef Full Text | Google Scholar

Grassmann, R., Modes, V., and Burgner-Kahrs, J. (2018). “Learning the Forward and Inverse Kinematics of a 6-DOF Concentric Tube Continuum Robot in SE(3),” in 2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Madrid, Spain, October 1–5, 2018 (IEEE), 5125–5132. doi:10.1109/IROS.2018.8594451

CrossRef Full Text | Google Scholar

Greenwood, D. T. (1988). “Principles of Dynamics,” in Prentice-hall International Series in Dynamics. 2nd ed. (Upper Saddle River: Prentice-Hall).

Google Scholar

Greer, J. D., Morimoto, T. K., Morimoto, A. M., and Okamura, E. W. (2019). A Soft, Steerable Continuum Robot that Grows via Tip Extension. Soft Robotics 6 (1), 95–108. doi:10.1089/soro.2018.0034

PubMed Abstract | CrossRef Full Text | Google Scholar

Greigarn, T., Jackson, R., Liu, T., and Cavusoglu, M. C. (2017). “Experimental Validation of the Pseudo-rigid-body Model of the MRI-Actuated Catheter,” in 2017 IEEE International Conference on Robotics and Automation (ICRA), 3600–3605. doi:10.1109/ICRA.2017.7989414

PubMed Abstract | CrossRef Full Text | Google Scholar

Greigarn, T., Poirot, N. L., Xu, X., and Cavusoglu, M. C. (2019). Jacobian-Based Task-Space Motion Planning for MRI-Actuated Continuum Robots. IEEE Robot. Autom. Lett. 4 (1), 145–152. doi:10.1109/LRA.2018.2881987

PubMed Abstract | CrossRef Full Text | Google Scholar

Hasanzadeh, S., and Janabi-Sharifi, F. (2014). An Efficient Static Analysis of Continuum Robots. J. Mech. Robotics 6 (3), 031011. doi:10.1115/1.4027305

CrossRef Full Text | Google Scholar

He, L., Lou, J., Dong, Y., Kitipornchai, S., and Yang, J. (2018). Variational Modeling of Plane-Strain Hyperelastic Thin Beams with Thickness-Stretching Effect. Acta Mech. 229 (12), 4845–4861. doi:10.1007/s00707-018-2258-4

CrossRef Full Text | Google Scholar

Jain, A., and Rodriguez, G. (1993). An Analysis of the Kinematics and Dynamics of Underactuated Manipulators. IEEE Trans. Robot. Automat. 9 (4), 411–422. doi:10.1109/70.246052

CrossRef Full Text | Google Scholar

Jolaei, M., Hooshiar, A., Dargahi, J., and Packirisamy, M. (2021). Toward Task Autonomy in Robotic Cardiac Ablation: Learning-Based Kinematic Control of Soft Tendon-Driven Catheters. Soft Robotics 8, 340–351. doi:10.1089/soro.2020.0006

PubMed Abstract | CrossRef Full Text | Google Scholar

Jung, J., Penning, R. S., and Zinn, M. R. (2014). A Modeling Approach for Robotic Catheters: Effects of Nonlinear Internal Device Friction. Adv. Robotics 28 (8), 557–572. doi:10.1080/01691864.2013.879371

CrossRef Full Text | Google Scholar

Kane, T. R., and Levinson, D. A. (1983). The Use of Kane's Dynamical Equations in Robotics. Int. J. Robotics Res. 2 (3), 3–21. doi:10.1177/027836498300200301

CrossRef Full Text | Google Scholar

Katzschmann, R. K., Santina, C. D., Toshimitsu, Y., Bicchi, A., and Rus, D. (2019a). “Dynamic Motion Control of Multi-Segment Soft Robots Using Piecewise Constant Curvature Matched with an Augmented Rigid Body Model,” in 2019 2nd IEEE International Conference on Soft Robotics (RoboSoft), Seoul, South Korea, April 14–18, 2019 (IEEE), 454–461. doi:10.1109/ROBOSOFT.2019.8722799

CrossRef Full Text | Google Scholar

Katzschmann, R. K., Thieffry, M., Goury, O., Kruszewski, A., Guerra, T.-M., Duriez, C., and Rus, D. (2019b). “Dynamically Closed-Loop Controlled Soft Robotic Arm Using a Reduced Order Finite Element Model with State Observer,” in 2019 2nd IEEE International Conference on Soft Robotics (RoboSoft), Seoul, South Korea, April 14–18, 2019 (IEEE), 717–724. doi:10.1109/ROBOSOFT.2019.8722804

CrossRef Full Text | Google Scholar

Kim, S., Laschi, C., and Trimmer, B. (2013). Soft Robotics: A Bioinspired Evolution in Robotics. Trends Biotechnol. 31 (5), 287–294. doi:10.1016/j.tibtech.2013.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Kreyszig, E. (1991). Differential Geometry. New York: Dover.

Google Scholar

Lai, J., Huang, K., and Chu, H. K. (2019). “A Learning-Based Inverse Kinematics Solver for a Multi-Segment Continuum Robot in Robot-independent Mapping,” in 2019 IEEE International Conference on Robotics and Biomimetics (ROBIO), Dali, China, December 6–8, 2019 (IEEE), 576–582. doi:10.1109/ROBIO49542.2019.8961669

CrossRef Full Text | Google Scholar

Li, J., Zhou, Y., Wang, C., Wang, Z., and Liu, H. (2020). A Model-Based Method for Predicting the Shapes of Planar Single-Segment Continuum Manipulators with Consideration of Friction and External Force. J. Mech. Robotics 12 (4), 041013. doi:10.1115/1.4046035

CrossRef Full Text | Google Scholar

Linn, J., Lang, H., and Tuganov, A. (2013). Geometrically Exact Cosserat Rods with Kelvin-Voigt Type Viscous Damping. Mech. Sci. 4 (1), 79–96. doi:10.5194/ms-4-79-2013

CrossRef Full Text | Google Scholar

Lock, J., and Dupont, P. E. (2011). “Friction Modeling in Concentric Tube Robots,” in 2011 IEEE International Conference on Robotics and Automation, Shanghai, China, May 9–13, 2011 (IEEE), 1139–1146. doi:10.1109/ICRA.2011.5980347

PubMed Abstract | CrossRef Full Text | Google Scholar

Lutter, M., Ritter, C., and Peters, Jan. (2018). “Deep Lagrangian Networks: Using Physics as Model Prior for Deep Learning,” in International Conference on Learning Representations, New Orleans, LA, May 6–9, 2019.

Google Scholar

Mahoney, A. W., Gilbert, H. B., and Webster, R. J. (2018). “A Review of Concentric Tube Robots: Modeling, Control, Design, Planning, and Sensing,” in The Encyclopedia of Medical Robotics (World Scientific), 181–202. doi:10.1142/9789813232266_0007

CrossRef Full Text | Google Scholar

Majidi, C. (2014). Soft Robotics: A Perspective—Current Trends and Prospects for the Future. Soft Robotics 1 (1), 5–11. doi:10.1089/soro.2013.0001

CrossRef Full Text | Google Scholar

Mauzé, B., Dahmouche, R., Laurent, G. J., André, A. N., Rougeot, P., Sandoz, P., et al. (2020). Nanometer Precision with a Planar Parallel Continuum Robot. IEEE Robot. Autom. Lett. 5 (3), 3806–3813. doi:10.1109/LRA.2020.2982360

CrossRef Full Text | Google Scholar

Miller, S., and Wendlandt, J. (2011). Real-time Simulation of Physical Systems Using Simscape. Available at: https://www.mathworks.com/company/newsletters/articles/real-time-simulation-of-physical-systems-using-simscape.html (Accessed September 24, 2021).

Google Scholar

Murray, R. M., Li, Z., and Sastry, S. S. (2017). A Mathematical Introduction to Robotic Manipulation. Boca Raton: CRC Press. doi:10.1201/9781315136370

CrossRef Full Text | Google Scholar

Mustaza, S. M., Elsayed, Y., Lekakou, C., Saaj, C., and Fras, J. (2019). Dynamic Modeling of Fiber-Reinforced Soft Manipulator: A Visco-Hyperelastic Material-Based Continuum Mechanics Approach. Soft Robotics 6 (3), 305–317. doi:10.1089/soro.2018.0032

PubMed Abstract | CrossRef Full Text | Google Scholar

Nicolai, A., Olson, G., Menguc, Y., and Hollinger, G. A. (2020). “Learning to Control Reconfigurable Staged Soft Arms,” in 2020 IEEE International Conference on Robotics and Automation (ICRA), Paris, France, May 31–August 31, 2020 (IEEE), 5618–5624. doi:10.1109/ICRA40945.2020.9197516

CrossRef Full Text | Google Scholar

Orekhov, A. L., and Simaan, N. (2020). “Solving Cosserat Rod Models via Collocation and the Magnus Expansion,” in 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Las Vegas, NV, October 24–January 24, 2021 (IEEE), 8653–8660. doi:10.1109/IROS45743.2020.9340827

CrossRef Full Text | Google Scholar

Parvaresh, A., and Moosavian, S. A. A. (2019). “Linear vs. Nonlinear Modeling of Continuum Robotic Arms Using Data-Driven Method,” in 2019 7th International Conference on Robotics and Mechatronics (ICRoM), Tehran, Iran, November 20–21, 2019 (IEEE), 457–462. doi:10.1109/ICRoM48714.2019.9071914

CrossRef Full Text | Google Scholar

Penning, R. S., and Zinn, M. R. (2014). A Combined Modal-Joint Space Control Approach for Continuum Manipulators. Adv. Robotics 28 (16), 1091–1108. doi:10.1080/01691864.2014.913503

CrossRef Full Text | Google Scholar

Platus, D. L. (1999). “Negative-Stiffness-Mechanism Vibration Isolation Systems,” in Proc. SPIE 3786, Optomechanical Engineering and Vibration Control (Denver: SPIE). doi:10.1117/12.363841

CrossRef Full Text | Google Scholar

Rao, P., Peyron, Q., Lilge, S., and Burgner-Kahrs, J. (2021). How to Model Tendon-Driven Continuum Robots and Benchmark Modelling Performance. Front. Robot. AI 7 (February), 630245. doi:10.3389/frobt.2020.630245

PubMed Abstract | CrossRef Full Text | Google Scholar

Renda, F., Cacucciolo, V., Dias, J., and Seneviratne, L. (2016). “Discrete Cosserat Approach for Soft Robot Dynamics: A New Piece-Wise Constant Strain Model with Torsion and Shears,” in 2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Daejeon, South Korea, October 9–14, 2016 (IEEE), 5495–5502. doi:10.1109/IROS.2016.7759808

CrossRef Full Text | Google Scholar

Renda, F., Giorelli, M., Calisti, M., Cianchetti, M., and Laschi, C. (2014). Dynamic Model of a Multibending Soft Robot Arm Driven by Cables. IEEE Trans. Robot. 30 (5), 1109–1122. doi:10.1109/TRO.2014.2325992

CrossRef Full Text | Google Scholar

Roark, R. J., Young, W. C., and Budynas, R. G. (2002). Roark’s Formulas for Stress and Strain. 7th ed. McGraw-Hill.

Google Scholar

Robinson, G., and Davies, J. B. C. (1999). “Continuum Robots-A State of the Art.,” in Proceedings 1999 IEEE International Conference on Robotics and Automation (Cat. No. 99CH36288C), 4, Detroit, MI, May 10–15, 1999 (IEEE), 2849–2854.

Google Scholar

Rodrigues, O. (1840). “Des Lois Géométriques Qui Régissent Les Déplacements d’un Système Solide Dans l’espace, et de La Variation Des Coordonnées Provenant de Ces Déplacements Considérés Indépendamment Des Causes Qui Peuvent Les Produire. J. de Mathématiques Pures Appliquées 5, 380–440.

Google Scholar

Rone, W. S., and Ben-Tzvi, P. (2014). Continuum Robot Dynamics Utilizing the Principle of Virtual Power. IEEE Trans. Robot. 30 (1), 275–287. doi:10.1109/TRO.2013.2281564

CrossRef Full Text | Google Scholar

Roy, R., Wang, L., and Simaan, N. (2017). Modeling and Estimation of Friction, Extension, and Coupling Effects in Multisegment Continuum Robots. Ieee/asme Trans. Mechatron. 22 (2), 909–920. doi:10.1109/TMECH.2016.2643640

CrossRef Full Text | Google Scholar

Rucker, D. C., and Webster III, R. J. (2011). Statics and Dynamics of Continuum Robots with General Tendon Routing and External Loading. IEEE Trans. Robot. 27 (6), 1033–1044. doi:10.1109/TRO.2011.2160469

CrossRef Full Text | Google Scholar

Rucker, D. C., Webster, R. J., Chirikjian, G. S., and Cowan, N. J. (2010). Equilibrium Conformations of Concentric-Tube Continuum Robots. Int. J. Robotics Res. 29 (10), 1263–1280. doi:10.1177/0278364910367543

PubMed Abstract | CrossRef Full Text | Google Scholar

Rucker, D. C., and Webster, R. J. (2011). “Computing Jacobians and Compliance Matrices for Externally Loaded Continuum Robots,” in 2011 IEEE International Conference on Robotics and Automation, Shanghai, China, May 9–13, 2011 (IEEE), 945–950. doi:10.1109/ICRA.2011.5980351

CrossRef Full Text | Google Scholar

Runge, G., Wiese, M., Gunther, L., and Raatz, A. (2017). “A Framework for the Kinematic Modeling of Soft Material Robots Combining Finite Element Analysis and Piecewise Constant Curvature Kinematics,” in 2017 3rd International Conference on Control, Automation and Robotics (ICCAR) (Nagoya, Japan: IEEE), 7–14. doi:10.1109/ICCAR.2017.7942652

CrossRef Full Text | Google Scholar

Sadati, S., Ali, S., Naghibi, S. E., Rucker, D. C., Renson, L., Bergeles, C., et al. (2019). “Reduced Order vs. Discretized Lumped System Models with Absolute and Relative States for Continuum Manipulators,” in Robotics: Science and Systems, 10. doi:10.15607/rss.2019.xv.076

CrossRef Full Text | Google Scholar

Satheeshbabu, S., Uppalapati, N. K., Chowdhary, G., and Krishnan, G. (2019). “Open Loop Position Control of Soft Continuum Arm Using Deep Reinforcement Learning,” in 2019 International Conference on Robotics and Automation (ICRA) (Montreal, QC, Canada: IEEE), 5133–5139. doi:10.1109/ICRA.2019.8793653

CrossRef Full Text | Google Scholar

Schaeffer, D. G., and Cain, J. W. (2016). “Ordinary Differential Equations: Basics and beyond,” in Texts in Applied Mathematics (Springer), Availabe at: http://libezp.lib.lsu.edu/login?url=https://search.ebscohost.com/login.aspx?direct=true&db=cat00252a&AN=lalu.4596309&site=eds-live&scope=site&profile=eds-main. doi:10.1007/978-1-4939-6389-8

CrossRef Full Text | Google Scholar

Schlagenhauf, C., Bauer, D., Chang, K.-H., King, J. P., Moro, D., Coros, S., and Pollard, N. (2018). “Control of Tendon-Driven Soft Foam Robot Hands,” in 2018 IEEE-RAS 18th International Conference on Humanoid Robots (Humanoids), Beijing, China, November 6–9, 2018 (IEEE), 1–7. doi:10.1109/HUMANOIDS.2018.8624937

CrossRef Full Text | Google Scholar

Sedal, A., Wineman, A., Gillespie, R. B., and Remy, C. D. (2021). Comparison and Experimental Validation of Predictive Models for Soft, Fiber-Reinforced Actuators. Int. J. Robotics Res. 40 (1), 119–135. doi:10.1177/0278364919879493

CrossRef Full Text | Google Scholar

Shiva, A., Sadati, S. M. H., Noh, Y., Fraś, J., Ataka, A., Würdemann, H., et al. (2019). Elasticity versus Hyperelasticity Considerations in Quasistatic Modeling of a Soft Finger-Like Robotic Appendage for Real-Time Position and Force Estimation. Soft Robotics 6 (2), 228–249. doi:10.1089/soro.2018.0060

PubMed Abstract | CrossRef Full Text | Google Scholar

Simo, J. C., and Vu-Quoc, L. (1986). A Three-Dimensional Finite-Strain Rod Model. Part II: Computational Aspects. Comput. Methods Appl. Mech. Eng. 58 (1), 79–116. doi:10.1016/0045-7825(86)90079-4

CrossRef Full Text | Google Scholar

Spong, M. W., Hutchinson, S., and Vidyasagar, M. (2006). Robot Modeling and Control. Hoboken: John Wiley & Sons.

Google Scholar

Spong, M. W., and Laurent, P. (1997). “Control of Underactuated Mechanical Systems Using Switching and Saturation,” in Control Using Logic-Based Switching. Editor A. Stephen Morse (Berlin, Heidelberg: Springer Berlin Heidelberg), 162–172.

Google Scholar

Spong, M. W. (1998). “Underactuated Mechanical Systems,” in Control Problems in Robotics and Automation (Springer), 135–150.

Google Scholar

Su, H.-J. (2009). A Pseudorigid-Body 3R Model for Determining Large Deflection of Cantilever Beams Subject to Tip Loads. J. Mech. Robotics 1 (2), 021008. doi:10.1115/1.3046148

CrossRef Full Text | Google Scholar

Sucan, I., and Kay, J. (2019). Urdf - ROS Wiki. Available at: http://wiki.ros.org/urdf (Accessed September 24, 2021).

Google Scholar

Thuruthel, T. G., Falotico, E., Renda, F., Laschi, C., and Laschi, C. (2017). Learning Dynamic Models for Open Loop Predictive Control of Soft Robotic Manipulators. Bioinspir. Biomim. 12 (6), 066003. doi:10.1088/1748-3190/aa839f

PubMed Abstract | CrossRef Full Text | Google Scholar

Thuruthel, T. G., Falotico, E., Renda, F., and Laschi, C. (2019). Model-Based Reinforcement Learning for Closed-Loop Dynamic Control of Soft Robotic Manipulators. IEEE Trans. Robot. 35 (1), 124–134. doi:10.1109/TRO.2018.2878318

CrossRef Full Text | Google Scholar

Till, J., Aloi, V., and Rucker, C. (2019). Real-Time Dynamics of Soft and Continuum Robots Based on Cosserat Rod Models. Int. J. Robotics Res. 38 (6), 723–746. doi:10.1177/0278364919842269

CrossRef Full Text | Google Scholar

Till, J., Bryson, C. E., Chung, S., Orekhov, A., and Rucker, D. C. (2015). “Efficient Computation of Multiple Coupled Cosserat Rod Models for Real-Time Simulation and Control of Parallel Continuum Manipulators”, in 2015 IEEE International Conference on Robotics and Automation (ICRA), Seattle, WA, May 26–30, 2015 (IEEE), 5067–5074. doi:10.1109/ICRA.2015.7139904

CrossRef Full Text | Google Scholar

Trivedi, D., Lotfi, A., and Rahn, C. D. (2008). Geometrically Exact Models for Soft Robotic Manipulators. IEEE Trans. Robot. 24 (4), 773–780. doi:10.1109/TRO.2008.924923

CrossRef Full Text | Google Scholar

Truby, R. L., Santina, C. D., and Rus., D. (2020). Distributed Proprioception of 3D Configuration in Soft, Sensorized Robots via Deep Learning. IEEE Robot. Autom. Lett. 5 (2), 3299–3306. doi:10.1109/LRA.2020.2976320

CrossRef Full Text | Google Scholar

Venkiteswaran, V. K., Sikorski, J., and Misra, S. (2019). Shape and Contact Force Estimation of Continuum Manipulators Using Pseudo Rigid Body Models. Mechanism Machine Theor. 139 (9), 34–45. doi:10.1016/j.mechmachtheory.2019.04.008

CrossRef Full Text | Google Scholar

Walker, I. D., Choset, H., and Chirikjian, G. S. (2016). “Snake-Like and Continuum Robots”, in Springer Handbook of Robotics. Editors B. Siciliano, and O. Khatib (Cham: Springer International Publishing), 481–498. doi:10.1007/978-3-319-32552-1_20

CrossRef Full Text | Google Scholar

Wang, L., Zheng, D., Harker, P., Patel, A. B., Guo, C. F., and Zhao, X. (2021). Evolutionary Design of Magnetic Soft Continuum Robots. Proc. Natl. Acad. Sci. USA 118 (21), e2021922118. doi:10.1073/pnas.2021922118

PubMed Abstract | CrossRef Full Text | Google Scholar

Webster, R. J., and Jones., B. A. (2010). Design and Kinematic Modeling of Constant Curvature Continuum Robots: A Review. Int. J. Robotics Res. 29 (13), 1661–1683. doi:10.1177/0278364910368147

CrossRef Full Text | Google Scholar

Webster, R. J., and Rucker., D. C. (2009). Parsimonious Evaluation of Concentric-Tube Continuum Robot Equilibrium Conformation. IEEE Trans. Biomed. Eng. 56 (9), 2308–2311. doi:10.1109/TBME.2009.2025135

PubMed Abstract | CrossRef Full Text | Google Scholar

Wehrmeyer, J. A., Barthel, J. A., Roth, J. P., and Saifuddin, T. (1998). Colonoscope Flexural Rigidity Measurement. Med. Biol. Eng. Comput. 36 (4), 475–479. doi:10.1007/bf02523217

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, D., Zhang, Y., Ourak, M., Niu, K., Dankelman, J., and Poorten, E. V. (2021). Hysteresis Modeling of Robotic Catheters Based on Long Short-Term Memory Network for Improved Environment Reconstruction. IEEE Robot. Autom. Lett. 6 (2), 2106–2113. doi:10.1109/LRA.2021.3061069

CrossRef Full Text | Google Scholar

Xu, K., and Simaan, N. (2008). An Investigation of the Intrinsic Force Sensing Capabilities of Continuum Robots. IEEE Trans. Robot. 24 (3), 576–587. doi:10.1109/TRO.2008.924266

CrossRef Full Text | Google Scholar

Xu, W., Jie Chen, Jie., Lau, H. Y. K., and Ren., H. (2016). “Automate Surgical Tasks for a Flexible Serpentine Manipulator via Learning Actuation Space Trajectory from Demonstration,” in 2016 IEEE International Conference on Robotics and Automation (ICRA), Stockholm, Sweden, May 16–21, 2016 (IEEE), 4406–4413. doi:10.1109/ICRA.2016.7487640

CrossRef Full Text | Google Scholar

Yip, M. C., and Camarillo, D. B. (2014). Model-Less Feedback Control of Continuum Manipulators in Constrained Environments. IEEE Trans. Robot. 30 (4), 880–889. doi:10.1109/TRO.2014.2309194

CrossRef Full Text | Google Scholar

Yip, M. C., and Camarillo, D. B. (2016). Model-Less Hybrid Position/Force Control: A Minimalist Approach for Continuum Manipulators in Unknown, Constrained Environments. IEEE Robot. Autom. Lett. 1 (2), 844–851. doi:10.1109/LRA.2016.2526062

CrossRef Full Text | Google Scholar

Zhang, C., He, B., Ding, A., Xu, S., Wang, Z., and Zhou, Y. (2019). Motion Simulation of Ionic Liquid Gel Soft Actuators Based on CPG Control. Comput. Intelligence Neurosci. 2019 (2), 1–11. doi:10.1155/2019/8256723

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., and Simaan, N. (2013). Design of Underactuated Steerable Electrode Arrays for Optimal Insertions. J. Mech. Robotics 5 (1). doi:10.1115/1.4007005

CrossRef Full Text | Google Scholar

Zhao, W., Zhang, Y., and Wang, N. (2021). Development and Performance Analysis of Pneumatic Soft-Bodied Bionic Actuator. Appl. Bionics Biomech. 2021, 1–13. doi:10.1155/2021/6623059

CrossRef Full Text | Google Scholar