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 ∈ D ⊂ ℝ at time t ∈ ℝ as a position vector p b ( s , t ) ∈ E with respect to a specific frame of reference ℱ 0 in E . The differentiability requirement on p b is always at least that the first derivative of 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.

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 p b is placed on this reference configuration. The curve is then “framed” by a set of unit vectors d 1 ( s , t ) ,   d 2 ( s , t ) , and d 3 ( s , t ) 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 d 3 ( s , t ) = d 1 ( s , t ) × d 2 ( s , t ) . In terms of classical differential geometry, d 3 is the tangent vector, and d 1 and 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 ℱ ( s , t ) for each point along the curve. The unit vectors equivalently define a spatiotemporal field of rotation operators R b ( s , t ) ∈ S O ( 3 ) . 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 g b ( s , t ) ∈ S E ( 3 ) is defined by. g b ( s , t ) = { R b ( s ,   t ) ,   p b ( s , t ) } .

The vectors u ( s , t ) = u k d k and v ( s , t ) = v k d k are termed the “strain variables.” They describe deformation of the body and are invariant under rigid transformations. The vector 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 ∂ s p 2 = v ⋅ v . The strain variables are related to the framed curve by the following relationships. ∂ s p b ( s , t ) = v ( s , t ) ,           ∂ s d k ( s , t ) = u ( s , t ) × d k ( s , t ) . (1)

Finally, the vectors ∂ t p b ( s , t ) and ω ( s , t ) represent the linear and angular velocity of the backbone curve and director vectors. The angular velocity satisfies the equation ∂ t d k ( s , t ) = ω ( s , t ) × d k ( s , t ) . The four functions u , v , ∂ t p b , and ω are not independent; they must satisfy ∂ s ω = ∂ t u + u × ω . In the reference configuration, the flexure strains have non-zero values u 0 ( s ) if the backbone is not a straight. Generally, v 0 ( s ) = d 3 ( s ) , but other choices are possible.

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 ℱ 0 . Three coordinates X ∈ ℝ 3 identify the position of each point in the body via a one-to-one, differentiable vector-valued function P ( X ) . If X is chosen as the Cartesian coordinates with respect to ℱ 0 , then this function and its inverse are trivial. The final locations of the points are described by p ( X , t ) . In some cases, it is useful to define a displacement field U as follows. p ( X ,   t ) = P ( X ) + U ( X ,   t ) .

The amount of stretching can be quantified by the deformation gradient, defined by F ( X , t ) = ∂ p ∂ P | 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 X = ( s , X 2 , X 3 ) and p ( X , t ) = p b ( s ,   t ) + X 2 d 2 ( s , t ) + X 3 d 3 ( s , t ) (Antman, 2005).

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 q ( t ) ∈ ℝ N by assuming that some kinematic quantity is described by a truncated sum of “modal” shape functions ϕ i ( s ) in a manner analogous to a Fourier series. The general form is to have a scalar quantity S j represented as S j ( s , t ) = ∑ i = 1 N a i j ( q ( t ) ) ϕ i ( s ) .

The function a i j may be simply an index into the vector 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 ( k i π s ) and cos ( k i π s ) for a series of values k i ∈ ℝ (directly analogous to a truncated Fourier series), the standard monomials { 1 , s , s 2 , … } , 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 d 3 can be expressed in spherical angles θ ( s , t ) and ϕ ( s ,   t ) , and then S 1 = θ and S 2 = ϕ , and v = ( 0,0,1 ) 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 ℱ 0 . Parameterizations also possible which represent the strain variables u ( s , t ) and v ( s , t ) measured by an observer in the local frame ℱ ( s , t ) (Boyer et al., 2020). When coupled with a collocation method used to determine 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.

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 Γ e = [ s e − 1 ,   s e ] with D = ∪ e = 1 m   Γ e and an approximation of the flexural-torsional strain as the following sum: u ( s , t ) = ∑ e = 1 m u e ( t )   χ Γ e ( s )   ,     χ Γ e ( s ) = { 1 s ∈ Γ e 0 s ∉ Γ e v ( s , t ) = d 3 ( s ) .

A similar approach is called the piecewise constant strain (PCS) method and extends the definition to include an approximating sum representing v (Renda et al., 2016). The vectors are resolved in Cartesian coordinates of the frame ℱ ( s , t ) 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 p b and d k (equivalently R b ) are easily computed using an explicit recursion presuming that a single boundary condition on the pose of the curve, g ( s 0 ,   t ) = g 0 ( t ) , is known. On element Γ i , the transformation from ℱ 0 to ℱ ( s , t ) is given by g ( s ,   t ) = g ( s e − 1 ,   t ) exp ( ( s − s e − 1 ) ξ e ( t ) ) , ξ e ( t ) = [ u e ( t ) × v e ( t ) 0 0 ] ℱ ( s , t ) .

The mapping exp : s e ( 3 ) → S E ( 3 ) is continuous and differentiable, including at the element 0 ∈ s e ( 3 ) , 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 ( [ ω × v 0 0 ] ) = [ R ( ω ) p ( ω ,   v ) 0 1 ]   R ( ω ) = I + C ( ω ) ω × + A ( ω ) ω × 2     p ( ω ,   v ) = C ( ω ) v + A ( ω ) ( ω × v ) + B ( ω ) ( ω ⋅ v ) ω .

The coefficient functions are as follows, with sinc ( x ) the un-normalized sinc function, which is available in the standard libraries of many programming languages. A ( ω ) = 1 2 sinc 2 ( ‖ ω ‖ 2 ) ,     B ( ω ) = ‖ ω ‖ − sin ( ‖ ω ‖ ) ‖ ω ‖ 3 ,     C ( ω ) = sinc ( ‖ ω ‖ ) .

In IEEE-754 double-precision arithmetic, the author has found that B ( ω ) may be accurately computed by a truncated Taylor series B ( ω ) ≈ 1 / 6 − ‖ ω ‖ 2 / 120 + ‖ ω ‖ 4 / 5040 if ‖ ω ‖ < 1 × 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 d k and in terms of 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 ℱ ( s , t ) 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.

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 g ( s 0 ,   t ) 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). [ u ( s , t ) v ( s , t ) ] ℱ ( s , t ) = ∪ e = 1 m ∑ j = 1 N e ϕ j ( z e ) q e j ( t ) ,           z e = 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 ( p + 1 ) 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 p and/or 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 S O ( 3 ) , 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.

More general finite-element descriptions have also been used to model soft and continuum robots. In this case, the degrees of freedom q ( t ) directly interpolate the position field p ( 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).

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).

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). u ℱ ( s , t ) ( s , t ) = ∑ i = 1 m q i ( t ) δ ( s − s i ) n i .

A universal joint is the result if two orthogonal axes n i and 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).

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. ∂ s y = f ( s , y ) ,     G a ( y ( 0 , t ) ,   t j ) = 0 ,     G b ( y ( L , t ) ,   t j ) = 0.

Solutions can then be provided by numerical codes which automatically determine the degrees of freedom used to approximate the function y ( s ,   t j ) 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 y ( s ,   t j ) 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).

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. [ ∂ t p ω ] ( s , t ) = J q ( s , t ) ∂ t q ( t ) = [ J p J ω ] ∂ t q ( t ) . (2)

Letting J i be the column multiplied by ∂ t q i ( t ) , then the column can be calculated from the following differential relationships: ∂ s J p i = ∂ q i v ∂ s J ω i = ∂ q i u + u × J ωi .

One must take care when the interpolation is carried out on the strain variables in coordinates of the local frame ℱ ( s , t ) . If desired, the coordinates in the body frame representations J p i ℱ ( s , t ) and J ω i ℱ ( s , t ) may be calculated from the appropriate representation of these equations in the moving frame (Rucker and Webster, 2011). ∂ s   ( J p i ℱ ( s , t ) ) = − u ℱ ( s , t ) × J p i ℱ ( s , t ) + ∂ q i ( v ℱ ( s , t ) ) + J ω i ℱ ( s , t ) × v ℱ ( s , t )   ∂ s   ( J ω i ℱ ( s , t ) ) = − u ℱ ( s , t ) × J ω i ℱ ( s , t ) + ∂ q i ( u ℱ ( s , t ) ) .

From a known boundary condition where J p i ( 0 , t ) = 0 and J ω i ( 0 , t ) = 0 , the solution to these equations can be expressed in closed form as the solution to a linear time varying system. [ J p i J ω i ] ( s ) = ∫ 0 s exp ( − ∫ 0 r ad ξ ( τ ) d τ )   ∂ q [ v ℱ ( s , t ) u ℱ ( s , t ) ] d r ,     ad ξ ( τ ) = [ u × v × 0 u ×   ] ℱ ( s , t ) .

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). ∂ s n + f = ρ A ∂ t t p + ρ I k ∂ t t d k (3) ∂ s m + ∂ s p × n + ℓ = ρ I k d k × ∂ t t p + ∂ t ( ρ J ⋅ ω ) . (4)

The sum from k = 1 to 3 is implied over the terms involving I k and d k . The variables n ( s , t ) and m ( s , t ) 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). f and ℓ are externally applied forces and moments. Applied concentrated forces and moments require the Dirac δ distribution to express in this formulation.

In the case of a model which allows freedom in all the strain variables, m and 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, n is a basic unknown and is equivalent to a Lagrange multiplier which enforces the constraint v ( s , t ) = v 0 ( s ) .

The parameter ρ A is the mass density (expressed per unit length) of the cross-section. ρ J is the mass moment of inertia (per unit length) of the section, which makes ρ J ⋅ ω the angular momentum (per unit length) calculated about the mass center of the section. The three parameters ρ 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 ρ 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.

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). M ( q ) q ¨ + C ( q , q ˙ ) q ˙ + N ( q , q ˙ ) = B ( τ ) .

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 τ 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 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). M ( q ) q ¨ + F ( q , q ˙ ) + G ( q ) = H T λ

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

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 q ∈ ℝ 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. ∫ 0 L [ ( F ( s , t ) + F * ( s , t ) ) ⋅ γ j ( s , t ) + ( M ( s , t ) + M * ( s , t ) ) ⋅ β j ( s , t ) ]   d s = Q j , n c j = 1 , … , N F ( s , t ) = − ∂ s n ( s , t ) , F * ( s , t ) = ρ A ∂ t t p ( s , t ) + ρ I k ∂ t t d k ( s , t ) M ( s , t ) = − ∂ s m + ∂ s p × n ,    M * ( s , t ) = ρ I k d k × ∂ t t p + ∂ t ( ρ J ⋅ ω ) Q j , n c = ∫ 0 L f ( s , t ) ⋅ γ j ( s , t ) + ℓ ( s , t ) ⋅ β j ( s , t )   d τ . (5)

The velocity coefficient function and angular velocity coefficient function are defined as γ j ( s , t ) = ∂ q j p ( s , t ) = ∂ q ˙ j ∂ t p ,    β j ( s , t ) = ∂ q ˙ j ω .

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

Since the time derivatives of the momentum density and angular momentum density, ∂ t t p and ∂ t ( ρ J ⋅ ω ) , can be written as linear functions of the ∂ t t 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 ∂ s n and ∂ s m , if resolved in ℱ ( s , t ) , are zero except at the element boundaries. The equations may be integrated by parts into a form which shows the conjugacy of n and v and the conjugacy of m and u . ∫ 0 L ∂ s n ⋅ γ j   d s = [ n ⋅ γ j ] 0 L − ∫ 0 L n ⋅ ∂ q j v   d s ∫ 0 L ∂ s m ⋅ β j d s = [ m ⋅ β j ] 0 L − ∫ 0 L m ⋅ ( ∂ q j u + u × β j ) d s .

In the local frame, the equations take the following forms. ∫ 0 L ∂ s n ⋅ γ j   d s = [ n ℱ ( s , t ) ⋅ γ j ℱ ( s , t ) ] 0 L − ∫ 0 L n ℱ ( s , t ) ⋅ ( ∂ q j ( v ℱ ( s , t ) ) + β j ℱ ( s , t ) × v ℱ ( s , t ) ) d s ∫ 0 L ∂ s m ⋅ β j ds = [ m ℱ ( s , t ) ⋅ β j ℱ ( s , t ) ] 0 L − ∫ 0 L m ℱ ( s , t ) ⋅ ∂ q j ( u ℱ ( s , t ) )   d s .

Note also that if n ( s , t ) 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, 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 ∂ t τ to the time derivatives ∂ t p ( L , t ) and ∂ t d k ( L , t ) , 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 x = F ( q ) provided by sensors, where it is presumed that the dimension of x and q are the same and that F is invertible, the dynamics of the system can be posed as a one-to-one mapping ( τ , x , ∂ t x ) → ∂ t t 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 ( τ , x , 0 ) → 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.

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 ∂ t q . G i ( τ i ,   q , ∂ t q ,   Q n c ) = 0.

However, the nature of the model may change depending on the exact form of G i . If G i involves only τ i and 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 Q n c = G ^ i ( τ i ,   q ,   ∂ t q ) , 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, τ = P is the fluid pressure, and 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 δ W = τ   δ V = τ   J v ( q ) δ q = Q i , n c δ q i Q n c = J V T ( q ) τ ,    J V ( q ) = ∂ 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 p t ( s , t ) = p ( s , t ) + a ( s ) with a ( s ) = a 1 ( s ) d 1 ( s , t ) + a 2 ( s ) d 2 ( s , t ) (Rucker and Webster, 2011). For the sake of simplicity, restrict the tendon to a planar path with a 2 ( s ) = 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: ℓ t ( q ) = ∫ 0 L d ℓ ,       d ℓ 2 = d s 2 ( d 3 + ∂ s a ) F t T F t ( d 3 + ∂ s a ) ,         F t = ∂ p t ∂ P t .   (6)

If the tendons are not fully constrained, other models for ℓ t ( q ) 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 Q n c corresponding to the Lagrange multiplier. δ W = τ δ ℓ t = τ   J ℓ ( q )   δ q = Q j , n c δ q j ,     J ℓ ( q ) = ∂ q ℓ t ( q ) .

Therefore, the effect of the tendon alone (not considering any frictional forces) must be Q n c = J ℓ T ( q ) τ .

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 ℓ t ( q ) 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.

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: m ℱ ( s , t ) = K b t ℱ ( s , t ) ( u ℱ ( s , t ) − u 0 ℱ ( s , t ) )   n ℱ ( s , t ) = K se ℱ ( s , t ) ( v ℱ ( s , t ) − v 0 ℱ ( s , t ) ) .

The matrices K b t and K s e 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 α ( s , t ) = ( E J α β ) ( s ) [ u β ( s , t ) − u β 0 ( s ) ] ,     m 3 = D ( s ) u 3 ( t ) .

Note that bending about d 1 and 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 α = ( G A ) ( s )   v α ,     n 3 = ( E A ) ( s )   [ v 3 − 1 ] .

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 ( u ℱ ( s , t ) , v ℱ ( s , t ) ) m ℱ ( s , t ) = ∂ u W ,    n ℱ ( s , t ) = ∂ 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.

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): m ℱ ( s , t ) = K b t ℱ ( s , t ) ( u ℱ ( s , t ) − u 0 ℱ ( s , t ) ) + B b t ℱ ( s , t )   ∂ t ( u ℱ ( s , t ) ) n ℱ ( s , t ) = K se ℱ ( s , t ) ( v ℱ ( s , t ) − v 0 ℱ ( s , t ) ) + B se ℱ ( s , t )   ∂ t ( v ℱ ( s , t ) ) .

The matrices K b t and B b t are related by time constants referred to as the extensional retardation time constant, τ e = η E / E , and the shear retardation time constant, τ s = η / G , with η E the “extensional viscosity” and η the shear viscosity. B b t ℱ ( s , t ) = K b t ℱ ( s , t ) ⋅ diag ( τ e ,   τ e ,   τ s ) B se ℱ ( s , t ) = K se ℱ ( s , t ) ⋅ diag ( τ s ,   τ s ,   τ e ) .

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).

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 d 3 ( 0 , t ) and d 1 ( 0 , t ) , let the boundary condition g ( 0 , t ) be fixed, and let the load be concentrated at s = L and modeled by F = P δ ( s − L ) d 1 ( 0 , t ) 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 θ , which represents the tangent vector: d 3 ( s , t ) = cos ( θ ) d 3 ( 0 , t ) + sin ( θ ) d 1 ( 0 , t ) ∂ ξ ξ θ + λ ⁡ cos ( θ ) = 0 ,     ξ = x L ,     λ = P L 2 E I     θ ( 0 , t ) = 0 ,     ∂ ξ θ ( 1 , t ) = 0.

The boundary value problem has a known solution: θ ( ξ , t ) = 2 ⁡ sin − 1 ( k sn ( K ( k ) − ( 1 − ξ ) λ ; k 2 ) ) − π 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 θ ( 0 ) = 0 k sn ( K ( k ) − λ ; k 2 ) = 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 E I = 1 .

Piecewise constant curvature: u PCC = ∑ e = 1 3 u e   χ Γ e ( s )   d 2 ,     χ Γ e ( s ) = { 1 s ∈ Γ e 0 s ∉ Γ e ,     D = [ 0 , 1 3 ]   ∪   [ 1 3 , 2 3 ]   ∪   [ 2 3 , 1 ]   m ( s ) = E I   u ( s ) .

PRB: u PRB = ∑ i = 1 3 q i δ ( s − s i )   d 2 ,     s 1 = 0.125 ,     s 2 = 0.475 ,   s 3 = 0.863 ,     D = [ 0,1 ]   m ( s 1 ) q 1 = 3.25 ,     m ( s 2 ) q 2 = 2.84 ,     m ( s 3 ) q 3 = 2.95.

Spectral: u S = ∑ i = 1 3 q i s i − 1   d 2 ,     D = [ 0,1 ] . m ( s ) = E I   u ( s )

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 λ . For cantilevered loads, a single spectral element which is quadratic in 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.

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 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.

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 ( s ) is a linear function, say for example 1 + a s , then in the simplest planar case with a constant internal moment one would be tasked to find another polynomial function κ ( s ) such that κ ( s ) K ( s ) = C for some constant C . But this is obviously impossible, because κ = 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 ⋅ ∑ ( − 1 ) i a i s i . The series does not converge unless | s | < | 1 / a | 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 | s | , ensuring convergence. For this reason, single-element, spectral-type methods are not recommended as a first choice if non-uniform flexural rigidity is present.

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.

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). ( E I y , s s ) , s s = − k y .

The homogeneous solution to this equation has the following form. y h ( s ) = exp ( − β s ) [ C 1 ⁡ sin ( β s ) + C 2 ⁡ cos ( β s ) ] + exp ( β s ) [ C 3 ⁡ sin ( β s ) + C 4 ⁡ cos ( β s ) ] .

The constant β depends only on the properties of the beam and the foundation. β 2 = k 4 E I .

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 [ 0 , ℓ ] . Here we select the 2-norm and study the approximation error for 3^rd, 4^th, and 5^th order polynomials. Since β is related to the ratio of stiffnesses k and E I , and has dimension Length⁻¹, we restrict the range of the dimensionless group ℓ β to ( 0.1 ,   10 ) . This range is consistent with the idea of compliance matching as a form of embodied intelligence in biomedical continuum robots. Note that as ℓ β → 0 , the solution x h approaches a constant, which is easy to interpolate. As ℓ β → ∞ , 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 [ 0,1 ] for β ∈ ( 0.1 , 10 ) . x 1 = exp ( − β s ) sin ( β s )     x 2 = exp ( − β s ) cos ( β s ) .

See Figure 7B for examples with β = 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 = max s | p i ( s ) − x i ( s ) | max s | x i ( s ) |

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 ℓ β < 4.3 (Figure 8). As ℓ β increases beyond this value, the polynomial approximations to x 1 and x 2 begin to oscillate with increasing maximum error.

To put this in a practical perspective, a typical colonoscope has a linearized flexural rigidity of E I ≈ 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 β = 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.

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 q = f ( τ ) does not always exist because there may be two or more values of q which satisfy the equilibrium conditions given τ . This solution multiplicity is accompanied by a singular tangent stiffness matrix for some value of 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 F ( q , τ , Q n c ∗ ) = 0 where Q n c ∗ includes only those generalized forces which are not algebraically related with q and τ . Then a perturbation analysis yields the C-space Jacobians with respect to τ and Q n c ∗ δ F = ∂ q F  δ q + ∂ τ F  δ τ + ∂ Q n c F  δ Q n c ∗ = 0. δ q = − ( ∂ q F ) − 1 ( ∂ τ F  δ τ + ∂ Q n c F  δ Q n c ∗ ) .

It is evidently at configurations with singular ∂ q 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.

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.

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.

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).

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.

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 Q n c 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).

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 q ( t ) . 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.