A modal-space formulation for momentum observer contact estimation and effects of uncertainty for continuum robots

Abstract

Contact detection for continuum and soft robots has been limited in past works to statics or kinematics-based methods with assumed circular bending curvature or known bending profiles. In this paper, we adapt the generalized momentum observer contact estimation method to continuum robots. This is made possible by leveraging recent results for real-time shape sensing of continuum robots along with a modal-space representation of the robot dynamics. In addition to presenting an approach for estimating the generalized forces due to contact via a momentum observer, we present a constrained optimization method to identify the wrench imparted on the robot during contact. We also present an approach for investigating the effects of unmodeled deviations in the robot’s dynamic state on the contact detection method and we validate our algorithm by simulations and experiments. We also compare the performance of the momentum observer to the joint force deviation method, a direct estimation approach using the robot’s full dynamic model. We also demonstrate a basic extension of the method to multisegment continuum robots. Results presented in this work extend dynamic contact detection to the domain of continuum and soft robots and can be used to improve the safety of large-scale continuum robots for human–robot collaboration.

1. Introduction

Human–robot collaboration presents a host of sensing and situational awareness challenges in order to ensure worker safety. These challenges are exacerbated for in-situ collaborative robots (ISCRs) where robots must safely coexist and interact with a human co-worker within a confined space. An example of such an ISCR is shown in Figure 1. In Abah et al. (2022), we presented our work on achieving active safety using sensory disks ① along the length of the continuum robot. While these sensory disks sense proximity, contact, and contact force along the length of the robot, in this paper, we consider an alternate approach for contact estimation using the robot’s dynamics because we realize that achieving enhanced active safety requires several independent methods to ensure fail-safe operation. The scope of this paper is therefore focused on identifying the challenges of contact detection and wrench estimation for continuum robots and on presenting what we believe is the first adaptation of the generalized momentum observer (GMO) method for continuum robots.

Figure 1.

A collaborative robot for cooperative manipulation in confined spaces using continuum segments. Contact detection and wrench estimation is a key component of ensuring worker safety during human–robot interaction. ① Robot sensory skin as described in Abah et al. (2022).

2. Related works

2.1. General contact estimation

In past works, to address the challenge of contact detection and estimation, other researchers have also incorporated sensors into robots in order to estimate the location and magnitude of the wrench applied during contact. Examples of this include sensory skins such as Yamada et al. (2005), Abah et al. (2019), and Abah et al. (2022), strain gauges placed onto the body of the robot as in García et al. (2003) and Malzahn and Bertram (2014), force/torque sensors in the robot’s wrist and/or base as in Lu et al. (2005) and Shimachi et al. (2008), vibration and acoustic analysis as in Richmond and Pai (2000), Min et al. (2019), and Fan et al. (2020), and joint torque sensors as in Popov et al. (2017).

Sensor-based methods can suffer from dead zones where sensors are not present and can be cost-ineffective. Therefore, researchers explored methods for contact detection that do not rely on adding extra sensors to the robot. The earliest technique is the joint force/torque deviation (JFD) method given in Takakura et al. (1989), Suita et al. (1995), Yamada et al. (1996), and Bajo and Simaan (2010). This method is also referred to as direct estimation by some authors (e.g., in Haddadin et al. (2017)). The JFD method consists of subtracting model-predicted generalized forces from the sensed generalized forces (e.g., using joint torque sensing). The model-predicted generalized forces may be from a full dynamic model as in Suita et al. (1995) or from a statics model as in Bajo and Simaan (2010). A model of the joint level friction is sometimes added to increase the accuracy of the method. The downside of this method is that it relies on joint-level current (or torque) sensing and is highly sensitive to the accuracy of the dynamic/static model. This can lead to false positives due to modeling errors and/or unmodeled friction. For methods that rely on the full dynamic model to predict the generalized forces such as Takakura et al. (1989), Suita et al. (1995), and Yamada et al. (1996), a numerical second-derivative of the joint positions is required which can lead to extremely noisy estimation signals. Additionally, JFD works best for slow motions where viscous friction and the effect of unmodeled dynamics are small.

Due to the aforementioned issues with estimation noise, it may be desirable to use methods that do not rely on the second derivative of the joint position and provide estimation noise filtering. These methods are based on the fault detection literature such as Wünnenberg and Frank (1990), Freyermuth (1991), Caccavale and Walker (1997), Hammouri et al. (1999), and De Persis and Isidori (2001) and typically rely on knowledge of the robot’s dynamic parameters. They work by observing an unexpected deviation in some aspect of the robot’s dynamic state such as its energy, joint velocity, or generalized momentum. A detailed review of these methods is given in Haddadin et al. (2017). The most common of these methods is the generalized momentum observer (GMO). This method was first developed in De Luca and Mattone (2003) and has been used in a wide range of applications such as flying robots in Tomić et al. (2017), humanoids in Vorndamme et al. (2017), redundant manipulators in De Luca and Ferrajoli (2008) and Cho et al. (2012), and robots with series-elastic actuators in De Luca et al. (2006) and Kim et al. (2015). The GMO has also been combined with Kalman filters in Wahrburg et al. (2015, 2018)bib_wahrburg_et_al_2018, additional sensory input in Buondonno and De Luca (2016), and friction models in Lee et al. (2015) and Lee and Song (2016) to improve the performance of the observer. The performance of the observer under characterized dynamic model uncertainty has also been studied in Briquet-Kerestedjian et al. (2017) and Li et al. (2021). It should also be noted that the GMO bares some similarities to the time-delay control paradigm introduced in Youcef-Toumi and Ito (1988) as it relies on previous measurements of the system state to estimate the contributions of external disturbances.

The GMO method for contact estimation has been used predominantly for rigid-link robots with/without series-elastic actuators. A detailed review of contact detection and localization methods for continuum robots is presented next in the related works section. Despite its successful use for various robot architectures, adapting the GMO method has not been previously applied to continuum robots. Furthermore, investigating the effects of state uncertainty on the performance of the contact estimator of a continuum robot has not been carried out.

The contribution of this paper is in presenting a modal space formulation of the GMO and verification of the GMO method to variable curvature continuum robots while taking into account uncertainty in the dynamic state of the robot. Unlike serial rigid-link robots, the shape of continuum robots is heavily dependent on external loading. To enable successful adaptation of the GMO method for continuum robots, one must overcome the challenge due to the load-induced shape changes. One, therefore, needs a modeling framework that takes this information into account. For these reasons, we leverage our recent result on real-time shape sensing in Orekhov et al. (2023) along with a dynamic formulation in the modal curvature space to present a formulation of the GMO method that can be used for continuum and soft robots. Relative to prior works in GMO-based contact estimation, we present a unique constrained optimization formulation for solving the problem of contact wrench estimation. We also present a model of the effects of dynamic state uncertainty on the output of the GMO and leverage this result to make recommendations for setting detection thresholds. Lastly, we compare the performance of the GMO to the performance of the JFD method in both simulation and experiment.

The next section presents related works while focusing on the most relevant works on contact detection/estimation for continuum robots. The following sections present the shape sensing approach and a dynamic formulation in the space of curvature modal factors. We then present an adaptation of the GMO method and an investigation of state parameter uncertainty on expected performance. Finally, sections 11 and 12 present a simulation-based study and experimental validation on the isolated distal continuum segment of the robot shown in Figure 1. Lastly, in section 13 we present a preliminary extension of the method to multisegment continuum robots.

2.2. Contact estimation for continuum robots

In addition to the seminal works on collision detection reviewed above, the most relevant work to this paper is De Luca and Mattone (2003), which we build on and extend in this paper. Within the scope of contact detection/estimation for continuum robots, the most relevant works are Xu and Simaan (2008), Bajo and Simaan (2010), Bajo and Simaan (2012), and Santina et al. (2020). In Bajo and Simaan (2010), the authors considered the JFD method and the effects of bounds of joint-force sensing uncertainty on the detectability of contact. In Bajo and Simaan (2012), a kinematics-based contact detection and localization approach based on motion tracking and detection in the shift of the instantaneous screws of motion was presented and shown to be quite effective for continuum robots that bend in circular shapes. Roy et al. (2016) simulated the effects of dynamics on this method assuming constant curvature bending profiles. In Rucker and Webster (2011), the authors combine Cosserat rod theory and extended Kalman filters to estimate the applied wrench to the tip of a continuum robot in quasi-static configurations. More recently, in Ashwin et al. (2021), the authors presented contact detection and estimation for wire-actuated continuum robots subject to quasi-static modeling assumptions. Lastly, in Leung et al. (2024), the authors present two tip-force estimation techniques for a two segment continuum robot using constant curvature assumptions. Their first technique uses a cantilever beam model to estimate the tip force and their second technique uses support vector regression to estimate the difference between no-load tendon forces and loaded tendon forces.

Table 1 summarizes the key modeling assumptions and capabilities supported by this work compared to the most relevant prior works. The limitations of these prior works on continuum robot contact estimation/detection are:

➣ The assumption of a known external load state (mostly unloaded) and known bending shape profiles of the continuum robot segments (e.g., most works in Table 1 assume constant curvature or use piece-wise constant curvature (PCC)). For example, Xu and Simaan (2008) assumed constant curvature and a force at the segment tip while Bajo and Simaan (2012) assumed constant curvature and detected contact along the length of a segment. Finally, Chen et al. (2021) assumed a planar curvature experimentally captured via a modal representation of the local tangent angle of a fiber-reinforced bellow actuator acting as a single-segment soft robot.

➣ The methods do not estimate the magnitude of the contact force and assume negligible dynamics. For example, Bajo and Simaan (2012) assumed negligible dynamics in their kinematics-based contact detection and localization approach. In Santina et al. (2020), the authors estimated contact forces in soft robots while assuming PCC, ignoring Coriolis/centrifugal and inertial effects, and assuming extrinsic sensing of the robot shape using optical tracking.

Table 1.

Summary of the continuum robot contact estimation literature.

Method	Dynamics included?	Contact localization?	Wrench estimation?	General curvature?	Intrinsic shape sensing?
Xu and Simaan (2008)			✓
Bajo and Simaan (2010)		✓
Bajo and Simaan (2012)		✓
Roy et al.(2016)	✓	✓			(Sim. only)
Rucker and Webster (2011)			✓	(PCC)	(Sim. only)
Santina et al. (2020)			✓	(PCC)
Chen et al. (2021)		✓		✓
Ashwin et al. (2021)		✓	✓
Zeng et al. (2023)			✓	✓	(Sim. only)
Leung et al. (2024)			✓
This work	✓		✓	✓	✓

While these approaches can work very well for small continuum robots with negligible dynamics and gravitation-induced deflections, large robots such as in Figure 1 require special consideration of their dynamics and the ensuing shape deflections.

The goal of this work is to address the limitations of prior works as shown in Table 1 via an extension/adaptation of the GMO method to continuum robots subject to dynamic loading and with general bending shape profiles. The robot shown in Figure 1 is very large and therefore its segments undergo significant passive deflections violating the constant curvature assumptions. This robot is also designed to work within confined spaces where external optical tracking may not be possible, therefore a combined approach achieving intrinsic shape sensing, contact detection, and contact wrench estimation is sought in this work.

Another set of key works that we build on and extend is the literature on continuum robot dynamics and kinematics. In this work, we build on Chirikjian’s (Chirikjian and Burdick, 1994, 1995; Chirikjian, 1994) early works on modal space representation of hyperredundant robots. We present a modal-space dynamics model for the purpose of formulating the GMO contact estimation approach and analyzing the effects of uncertainty.

3. Problem definition and modeling assumptions

This paper focuses on presenting an extension of the GMO contact estimation method to continuum robots. Because these robots are very different from traditional serial rigid-link robots, these robots have significant sources of geometric, dynamic, and mechanics uncertainty. This combination of challenges requires solutions for online shape sensing, updated kinematics, and updated dynamics in a way that allows reliable application of the GMO algorithm. Therefore, the two problems targeted by this paper are:

➣ Collision Estimation: Detect contact and estimate the wrench applied to the robot during contact.

➣ Effects of Uncertainty: Formulate the effects of uncertainty and provide an estimate of the effects of uncertainty on the GMO estimator, thereby delineating the limits of contact detectability due to uncertainty.

In dynamics model presented in this paper, we make the following assumptions:

➣ The friction in the actuation tendons can be modeled as concentrated forces at the contact points between the actuation tendons and the robot body.

➣ The friction can be modeled as a smooth Coulomb friction model

➣ Viscous effects that may arise due to gearhead losses can be neglected

➣ The robot is torsionally rigid. This assumption is justified since the robot shown in Figure 1 has bellows with very high torsional stiffness as reported in Orekhov et al. (2023).

➣ The central backbone of the robot includes a superelastic nickel-titanium (NiTi) rod. NiTi is known to exhibit a nonlinear, hysteretic stress-strain curve (Liu et al., 1999). Because we mostly operate in the second stress-strain curve plateau that occurs past approximately 2% strain (Liu et al., 1999), we neglect these effects and assume a linear elasticity model for simplicity.

➣ Although the method can handle more complex contact types, we evaluate the method in simulation and experiment using only point contacts.

4. Modal kinematics

The instantaneous shape of the continuum segment shown in Figure 2 will be represented using a modal representation following the model presented in Orekhov et al. (2023). For completeness, we briefly summarize the modal kinematics below while using the nomenclature presented in Table 2.

Figure 2.

(a) Actuation tendons and passive shape sensing strings whose lengths are used to calculate c. The figure is reproduced with permission from Orekhov et al. (2023). (b) Cross section of the robot showing the backbone arc length coordinate s, the local backbone frame {T(s)}, and the robot base frame {T(0)}.

Table 2.

Nomenclature.

Symbol	Description
{F}	Designates a right-handed frame with unit vectors ${\widehat{\mathbf{x}}}_f,{\widehat{\mathbf{y}}}_f,{\widehat{\mathbf{z}}}_f$ and f as its origin
{T(s)}	Local frame at segment arc length s. It has its ${\widehat{\mathbf{z}}}_t$ axis tangent to the central backbone of the continuum robot
${}^a\mathbf{R}_b$	A rotation matrix describing the orientation of frame {B} relative to frame {A}
(⋅)′	Is a shorthand notation for d(⋅)/ds
$\dot{(\cdot )}$	Is a shorthand notation for d(⋅)/dt, i.e., time derivative
x ∧	The exterior product matrix representation of a vector. ${\mathbf{x}}^{\wedge }=\left\{\begin{array}{cc}\left[\begin{array}{ccc}\hfill 0\hfill & \hfill -x_3\hfill & \hfill x_2\hfill \\ \hfill x_3\hfill & \hfill 0\hfill & \hfill -x_1\hfill \\ \hfill -x_2\hfill & \hfill x_1\hfill & \hfill 0\hfill \end{array}\right]\in so(3)\hfill & \mathrm{f}\mathrm{o}\mathrm{r}\mathbf{x}=\left[\begin{array}{c}\hfill x_1\hfill \\ \hfill x_2\hfill \\ \hfill x_3\hfill \end{array}\right]\hfill \\ \left[\begin{array}{cc}\hfill {\mathit{\omega }}^{\wedge }\hfill & \hfill \mathbf{v}\hfill \\ \hfill \mathbf{0}\hfill & \hfill 0\hfill \end{array}\right]\in se(3)\hfill & \mathbf{x}=\left[\begin{array}{c}\hfill \mathbf{v}\in {\text{IR}}^3\hfill \\ \hfill \mathit{\omega }\in {\text{IR}}^3\hfill \end{array}\right]\hfill \end{array}\right.$
${\left(\mathbf{X}\right)}^{\vee }$	An inverse operator of (⋅)∧ such that: ${\mathbf{X}}^{\vee }=\left\{\begin{array}{cc}{[x_1,x_2,x_3]}^{\text{T}},\hfill & \mathrm{f}\mathrm{o}\mathrm{r}\mathbf{X}\in so(3)\hfill \\ {[{\mathbf{v}}^{\text{T}},{\mathit{\omega }}^{\text{T}}]}^{\text{T}},\hfill & \mathrm{f}\mathrm{o}\mathrm{r}{\mathbf{X}}^{\wedge }=\left[\begin{array}{cc}\hfill {\mathit{\omega }}^{\wedge }\hfill & \hfill \mathbf{v}\hfill \\ \hfill \mathbf{0}\hfill & \hfill 0\hfill \end{array}\right]\in se(3)\hfill \end{array}\right.$
I n	An n × n identity matrix
Κ	A generalized force

The continuum segment of Figure 2 uses two capstans with two fixed-length wire loops to bend in any direction. The robot uses torsionally rigid metal bellows assembled on a superelastic NiTi rod. Since the bellows are torsionally very stiff, we assume the robot does not experience torsional strains and that the robot has a fixed arc length L. For a given arc length s ∈ [0, L] along a continuum segment, a local frame T(s) is assigned with its z-axis tangent to the backbone and following the increment direction of the arc-length (from the base to the segment tip):

${}^0\mathbf{T}_t(s)=\left[\begin{array}{cc}\hfill {}^0\mathbf{R}_t\left(s\right)\hfill & {}^0\mathbf{p}\left(s\right)\\ \hfill \mathbf{0}\hfill & \hfill 1\hfill \end{array}\right]\in SE(3)$

(1)

As the robot undergoes a deformation parameterized in the local frame as a twist ${}^t\mathit{\eta }(s)$ , the local frame described by ${}^0\mathbf{T}_t(s)$ can be found using the following differential equation:

${}^0\mathbf{T}_t^{\prime }(s)={}^0\mathbf{T}_t(s){}^t\mathit{\eta }^{\wedge }(s)$

(2)

The body twist ${}^t\mathit{\eta }(s)$ can be expressed using the local curvature of the backbone u(s):

${}^t\mathit{\eta }^{\wedge }(s)=\left[\begin{array}{cc}\hfill {\mathbf{u}}^{\wedge }(s)\hfill & \hfill {\mathbf{e}}_3\hfill \\ \hfill \mathbf{0}\hfill & \hfill 0\hfill \end{array}\right]\in se(3)$

(3)

where e₃ = [0, 0, 1]^T. This differential equation can be solved using either numerical integration or the Magnus expansion as in Orekhov and Simaan (2020). For the experiments and simulations presented herein, we used the 4th order Magnus expansion. For more background on the Magnus Expansion, see Magnus (1954) and Blanes et al. (2009).

In this work, we choose to represent u(s) using a weighted sum of polynomial basis functions. These polynomial functions ϕ _x (s) and ϕ _y (s) represent the curvature in the x and y directions, respectively. Because the robot is assumed to be perfectly rigid in torsion, there is zero curvature in the z direction (i.e. ϕ _z (s) = 0). Using these polynomial basis functions and the scalar weights c _x , c _y ∈ IR³, the robot’s curvature distribution u(s) can be written as:

$\mathbf{u}(s)=\left[\begin{array}{cc}\hfill {\mathit{\varphi }}_x^{\text{T}}(s)\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill {\mathit{\varphi }}_y^{\text{T}}(s)\hfill \\ \hfill \mathbf{0}\hfill & \hfill \mathbf{0}\hfill \end{array}\right]\left[\begin{array}{c}\hfill {\mathbf{c}}_x\hfill \\ \hfill {\mathbf{c}}_y\hfill \end{array}\right]=\mathbf{\Phi }(s)\mathbf{c}$

(4)

We chose these polynomial basis functions to be the first three Chebyshev polynomials of the first kind:

${\mathit{\varphi }}_x^{\text{T}}(s)={\mathit{\varphi }}_y^{\text{T}}(s)=\left[1,\frac{2s-L}{L},\frac{8s^2}{L^2}-\frac{8s}{L}+1\right]$

(5)

Chebyshev polynomials of the first kind were chosen over a simple monomial basis because Chebyshev polynomials are bounded to within ±1 and therefore improve scaling of the modal coefficients. Additionally, monomial bases are known to suffer from poor conditioning at higher powers as shown in Gautschi (1975). For a discussion on how the choice of Chebyshev polynomial order affects the accuracy of the kinematics model see Orekhov and Simaan (2020).

4.1. Shape sensing

As shown in Figure 2, the robot used in this work is endowed with four string encoders that measure the change in arc length between the end-disk and one of the spacer disks. As reported in Orekhov et al. (2023), these string encoders are custom-made using Renishaw AM4096 12 b incremental rotary encoders. These custom-made string encoders have an average measurement error of 0.1 mm (0.08% of the total stroke) and a maximum error of 0.27 mm (0.23% of total stroke). Using the string encoder measurements, as well as the change in backbone lengths, we can calculate the modal coefficients c that describe the curvature distribution u(s). To do this, we define the augmented space of tendon and sensory string extensions relative to the straight home configuration, that is, $\mathrm{\Delta }\mathit{\ell }\triangleq {\left[\mathrm{\Delta }{\ell }_1,\dots ,\mathrm{\Delta }{\ell }_4,\mathrm{\Delta }{\ell }_{q_1},\mathrm{\Delta }{\ell }_{q_2}\right]}^{\text{T}}$ . We then follow the derivation in Orekhov et al. (2023) to get the matrix J _ℓc ∈ IR^6×6 that maps the modal coefficients c to Δ ℓ :

$\mathrm{\Delta }\mathit{\ell }={\mathbf{J}}_{\ell c}\mathbf{c}$

(6)

Because the robot used in this paper is torsionally stiff and has constant pitch-radius string encoder/actuation tendon routing, J _ℓc is a constant, square, full-rank matrix (Orekhov et al., 2023) and can therefore be inverted to solve for c given measurements of Δ ℓ .

4.2. Modal body Jacobian

In this section, we will derive the Jacobian J _ξc (s) ∈ IR^6×6 that relates the temporal rates of the modal coefficients $\dot{\mathbf{c}}\in {\text{IR}}^6$ to the body twist coordinates such that ${}^t\mathit{\xi }(s)={\mathbf{J}}_{\xi c}(s)\dot{\mathbf{c}}$ .

The body twist of the end-effector can be found using:

${}^t\mathit{\xi }^{\wedge }(s)={}^0\mathbf{T}_t^{-1}(s){}^0\dot{\mathbf{T}}_t(s)\in se(3)$

(7)

After integration of (2), ${}^0\mathbf{T}_t$ becomes a function of curvature c, therefore, using the chain rule results in:

${}^t\mathit{\xi }^{\wedge }(s)=\sum _{i=1}^6\left[{}^0\mathbf{T}_t^{-1}(s)\left(\frac{\partial }{\partial c_i}{}^0\mathbf{T}_t(s){\dot{c}}_i\right)\right]$

(8)

where c _i is the ith modal coefficient. The ith column of J _ξc can be found using.

${\mathbf{J}}_{\xi c}^{[i]}(s)={[{}^0\mathbf{T}_t^{-1}(s)(\frac{\partial }{\partial c_i}{}^0\mathbf{T}_t(s))]}^{\vee }$

(9)

To calculate (9), we must compute the partial derivative of ${}^0\mathbf{T}_t(s)$ with respect to c _i . Since ${}^0\mathbf{T}_t(s)$ was computed using the Magnus expansion, $\partial /\partial c_i{}^0\mathbf{T}_t(s)$ is computed using the derivative of the matrix exponential (see Selig (2005) and Orekhov et al. (2023)).

4.3. Modal capstan Jacobian

In this section, we will derive the Jacobian J _cq (s) ∈ IR^6×2 that relates the angular velocity of the capstans $\dot{\mathbf{q}}\in {\text{IR}}^2$ to the velocity of the modal coefficients $\dot{\mathbf{c}}\in {\text{IR}}^6$ . The tendon speeds ${\dot{\mathit{\ell }}}_q={[{\dot{\ell }}_{q_1},{\dot{\ell }}_{q_2}]}^{\text{T}}$ can be calculated from $\dot{\mathbf{q}}$ using:

${\dot{\mathit{\ell }}}_q=\frac{1}{2\pi }\sqrt{{(2\pi r_c)}^2+{\gamma }^2}\dot{\mathbf{q}}$

(10)

where r _c is the capstan radius and γ is the capstan lead. Using J _ℓc , we can write the velocity of the modal coefficients that result from the angular velocity of the capstans:

$\dot{\mathbf{q}}=\underset{{\mathbf{J}}_{qc}}{\underbrace{2\pi {\left(\sqrt{{(2\pi r_c)}^2+{\gamma }^2}\right)}^{-1}\left[\begin{array}{cc}\hfill \mathbf{0}\hfill & \hfill {\mathbf{I}}_2\hfill \end{array}\right]{\mathbf{J}}_{\ell c}}}\dot{\mathbf{c}}$

(11)

where the matrix [0, I₂] accounts for the fact that only the last two values of $\dot{\mathit{\ell }}$ correspond to the control tendon lengths.

5. Continuum robot dynamics in the space of modal coefficients

Since the GMO method requires the robot’s dynamics, the dynamics of the continuum segment shown in Figure 2 is derived in the space of modal coefficients using the Euler–Lagrange method. The closest work to this section is the work of Roy et al. (2016). However, the dynamic model presented in Roy et al. (2016) assumes constant curvature (i.e., circular) bending. In this section, we will relax this assumption. The derivations provided in sections 5.1 and 5.2 are extended to general curvature from the constant-curvature models presented in Roy et al. (2016). The friction model presented in section 5.3 is unique to this paper.

5.1. Kinetic energy

The total kinetic energy of the robot can be modeled as the sum of the central backbone’s kinetic energy T _CB , the kinetic energy of the spacer disks T _SD , and the kinetic energy of the actuation capstans and motors T _A .

$T=T_{CB}+T_{SD}+T_A$

(12)

The models presented in sections 5.1.1 and 5.1.2 are adapted to general curvature from the constant curvature models presented in Roy et al. (2016). However, the derivations provided in section 5.1.3 are unique to this paper.

5.1.1. Central backbone

The kinetic energy of the central backbone can be found using the robot body twist coordinates ${}^t\mathit{\xi }(s)$ , the mass per unit length of the backbone ρ, and the backbone radius r:

$T_{CB}=\frac{1}{2}{\int }_0^L{}^t\mathit{\xi }^{\mathrm{T}}(s)\left[\begin{array}{cc}\hfill \rho \mathbf{I}\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill {\mathcal{I}}_{CB}\hfill \end{array}\right]{}^t\mathit{\xi }(s)ds$

(13)

In this equation, ${\mathcal{I}}_{CB}$ is the inertia tensor of a thin, circular cross-section of the backbone:

${\mathcal{I}}_{CB}=\left[\begin{array}{ccc}\hfill \frac{1}{4}\rho r^2\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \frac{1}{4}\rho r^2\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{1}{2}\rho r^2\hfill \end{array}\right]$

(14)

Using the modal Jacobian matrix J _ξc (s), we can write T _CB as a function of the modal coefficient velocities $\dot{\mathbf{c}}$ :

$T_{CB}=\frac{1}{2}{\dot{\mathbf{c}}}^{\text{T}}\underset{{\mathbf{M}}_{CB}}{\underbrace{{\int }_0^L{\mathbf{J}}_{\xi c}^{\text{T}}(s)\left[\begin{array}{cc}\hfill \rho \mathbf{I}\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill {\mathcal{I}}_{CB}\hfill \end{array}\right]{\mathbf{J}}_{\xi c}(s)ds}}\dot{\mathbf{c}}$

(15)

$T_{CB}=\frac{1}{2}{\dot{\mathbf{c}}}^{\text{T}}{\mathbf{M}}_{CB}\dot{\mathbf{c}}$

(16)

5.1.2. Spacer disks

In the following, we define for each spacer disk the local backbone frame {T $(s_{d_i})$ }, i ∈ [1, …, n _d ] where n _d denotes the number of spacer disks. Since each spacer disk may have its center of mass offset from {T $(s_{d_i})$ }, we define a center-of-mass body frame {U $(s_{d_i})$ } having its origin at the center of mass of the ith disk and parallel to {T $(s_{d_i})$ }. Using ${}^t\mathbf{p}_{c{m}_i}\in {\text{IR}}^3$ to denote the location of the disk center of mass in frame {T $(s_{d_i})$ }, we obtain the body twist of the ith spacer disk in a center of mass frame {U $(s_{d_i})$ } using the adjoint twist transformation:

${}^u\mathit{\xi }_{c{m}_i}=\underset{{\mathbf{S}}_{c{m}_i}}{\underbrace{\left[\begin{array}{cc}\hfill \mathbf{I}\hfill & \hfill -{}^t\mathbf{p}_{c{m}_i}^{\wedge }\hfill \\ \hfill \mathbf{0}\hfill & \hfill \mathbf{I}\hfill \end{array}\right]}}{}^t\mathit{\xi }(s_{d_i})$

(17)

Using $m_{d_i}$ and ${\mathcal{I}}_{d_i}$ to denote the mass and inertia tensor of the ith disk expressed in {U}, the kinetic energy of all disks is given by:

$T_{SD}=\frac{1}{2}\sum _{i=1}^{n_d}{}^u\mathit{\xi }_{c{m}_i}^{\mathrm{T}}\left[\begin{array}{cc}\hfill m_{d_i}\mathbf{I}\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill {\mathcal{I}}_{d_i}\hfill \end{array}\right]{}^u\mathit{\xi }_{c{m}_i}$

(18)

Using ${}^t\mathit{\xi }(s)={\mathbf{J}}_{\xi c}(s)\dot{\mathbf{c}}$ and (17), we write T _SD as a function of the modal coefficient rates:

$T_{SD}=\frac{1}{2}{\dot{\mathbf{c}}}^{\text{T}}{\mathbf{M}}_{SD}\dot{\mathbf{c}}$

(19)

where M _SD is given by:

${\mathbf{M}}_{SD}=\sum _{i=1}^{n_d}{\mathbf{J}}_{\xi c}^{\text{T}}(s_{d_i}){\mathbf{S}}_{c{m}_i}^{\text{T}}\left[\begin{array}{cc}\hfill m_{d_i}\mathbf{I}\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill {\mathcal{I}}_{d_i}\hfill \end{array}\right]{\mathbf{S}}_{c{m}_i}{\mathbf{J}}_{\xi c}(s_{d_i})$

(20)

5.1.3. Actuation capstans and motors

Each axis of the continuum robot actuation unit consists of a motor (Figure 3①) with rotor inertia J _m attached to gearhead (Figure 3②) with inertia J _gh as viewed by the gearhead input and reduction ratio R _gh . The gearhead is attached to a pinion (Figure 3③) with inertia J _p which meshes with a gear (Figure 3④) with inertia J _g . The gear is rigidly attached to a shaft and capstan (Figure 3⑤) with combined inertia J _c . Using this information, the inertia of the entire actuation chain as viewed by the capstan can be written as:

$J_A=R_g^2\left[R_{gh}^2(J_{gh}+J_m)+J_p\right]+J_g+J_c$

(21)

Figure 3.

Continuum robot actuation unit: ① Motor with rotor inertia J _m , ② gearhead with inertia J _gh and reduction ratio R _gh , ③ pinion with inertia J _p about the center of rotation, ④ gear with inertia J _g about the center of rotation, and ⑤ shaft and capstan with combined inertia J _c . The gear and pinion have a reduction ratio of R _gp .

The kinetic energy of the capstans T _A can be found using J _A and the angular velocity of the capstans $\dot{\mathbf{q}}\in {\text{IR}}^2$ :

$T_A=\frac{1}{2}{\dot{\mathbf{q}}}^{\text{T}}\left[\begin{array}{cc}\hfill J_A\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill J_A\hfill \end{array}\right]\dot{\mathbf{q}}$

(22)

T _A can be written in terms of the modal coefficient rates using:

$T_A=\frac{1}{2}{\dot{\mathbf{c}}}^{\text{T}}{\mathbf{M}}_A\dot{\mathbf{c}}$

(23)

where M _A is given by:

${\mathbf{M}}_A={\mathbf{J}}_{qc}^{\text{T}}\left[\begin{array}{cc}\hfill J_A\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill J_A\hfill \end{array}\right]{\mathbf{J}}_{qc}$

(24)

Using (16), (19), and (23), the total kinetic energy of the continuum robot can be written as:

$T=\frac{1}{2}{\dot{\mathbf{c}}}^{\text{T}}\underset{\mathbf{M}}{\underbrace{\left({\mathbf{M}}_{CB}+{\mathbf{M}}_{SD}+{\mathbf{M}}_A\right)}}\dot{\mathbf{c}}$

(25)

Note that, in reality, the capstans also slowly translate along their supporting shafts as they rotate. The kinetic energy from this effect is ignored in this model.

5.2. Potential energy

The total potential energy of the robot can be modeled as the sum of the bending energy V _B of the central backbone, the gravitational potential energy of the central backbone V _CB , and the gravitational potential energy of the spacer disks V _SD .

$V=V_B+V_{CB}+V_{SD}$

(26)

The derivations given in this section are extended to general curvature from the constant curvature models presented in Roy et al. (2016).

5.2.1. Elastic potential energy

The bending stiffness matrix of the central backbone can be written as:

$\mathbf{K}=\mathrm{diag}([E{I}_x,E{I}_y,JG])$

(27)

where EI _x is the flexural rigidity of the central backbone in the x direction, EI _y is the flexural rigidity in the y direction, and JG is the torsional rigidity of the backbone. Using the stiffness matrix, the potential energy due to bending of the central backbone can be written as:

$V_B=\frac{1}{2}{\int }_0^L{\mathbf{u}}^{\text{T}}(s)\mathbf{K}\mathbf{u}(s)ds$

(28)

Using (4), V _B can be written in terms of c as:

$V_B=\frac{1}{2}{\mathbf{c}}^{\text{T}}\underset{{\mathbf{K}}_c}{\underbrace{\left({\int }_0^L{\mathbf{\Phi }}^{\text{T}}(s)\mathbf{K}\mathbf{\Phi }(s)ds\right)}}\mathbf{c}$

(29)

5.2.2. Gravitational potential energy

The gravitational potential energy of the central backbone can be written as:

$V_{CB}=-\rho {\mathbf{g}}^{\text{T}}{\int }_0^L{}^0\mathbf{p}(s)ds$

(30)

where g = [0,0,−9.81]^T m/s² is the acceleration of gravity, ${}^0\mathbf{p}(s)$ is the position vector of the backbone at arc length s, and ρ is the mass per unit length of the backbone.

The gravitational potential energy of the spacer disks can be written as:

$V_{SD}=-\sum _{i=1}^{n_d}{m}_{d_i}{\mathbf{g}}^{\text{T}}\left[{}^0\mathbf{p}(s_{d_i})+{}^0\mathbf{R}_t(s_{d_i}){}^t\mathbf{p}_{{cm}_i}\right]$

(31)

where $m_{d_i}$ is the mass of the ith disk, ${}^0\mathbf{p}(s_{d_i})$ is the position of the central backbone at spacer disk anchoring arc length $s_{d_i}$ and ${}^0\mathbf{R}_t(s_{d_i})$ was defined in (1).

5.3. Approximate joint friction

In this section, we will approximate the effect of friction between the actuation tendons and the brass bushings that guide the actuation tendons through the spacer disks. Our approach, which is unique to this paper, uses a simplistic Coulomb friction model that ignores nonlinear, velocity-dependent effects such as Stribeck and viscous friction. As shown in Figures 2 and 4, the robot’s actuation tendons are wrapped around a motorized capstan and then routed through holes in the spacer disks. Once the tendon reaches the end-disk, they are rerouted back through the spacer disks using a pulley (see Figure 2) and terminated on a fixed shaft. This is repeated on the same capstan, but separated by 180° to control bending in both directions in the same plane. To ensure the tendons are always in tension, each side of the tendon are pre-tensioned to f _pl = 208 N.

Figure 4.

Continuum segment cross section showing the terms used to calculate the joint friction.

The pretensions cannot create a moment about the capstan because they act in equal and opposite directions on the capstan. ¹ The pretension does, however, affect the friction between the tendons and the bushings because it contributes to the normal force. When torque τ _j is applied to capstan j ∈ [1, 2], the tension in the tendon being pulled by the capstan is $f_{0,j}^{+}$ as shown in Figure 4. On the release side of the capstan j, the tendon tension is reduced to $f_{0,j}^{-}$ . The values of $f_{0,j}^{+}$ and $f_{0,j}^{-}$ can be found from τ _j and f _pl using:

$f_{0,j}^{+}=\frac{|{\tau }_j|}{r_c}+f_{pl},f_{0,j}^{-}=\max \left(f_{pl}-\frac{|{\tau }_j|}{r_c},0\right)$

(32)

As the tendon passes through the bushing in disk i, the tendon tension is reduced by the friction force between the bushing and the tendon ( ${\mathfrak{f}}_{i,j}^{+}$ or ${\mathfrak{f}}_{i,j}^{-}$ ):

$f_{i,j}^{+}=f_{i-1,j}^{+}-{\mathfrak{f}}_{i,j}^{+},f_{i,j}^{-}=f_{i-1,j}^{-}-{\mathfrak{f}}_{i,j}^{-}$

(33)

In this section, we will use the convention that i = 1 at the segment base as shown in Figure 4. The friction forces depend on the magnitude of normal force between the tendon and the bushing ( ${\mathfrak{n}}_{i,j}^{+}$ or ${\mathfrak{n}}_{i,j}^{-}$ ) which can be found by projecting the tendon forces onto the plane of the disk:

$\begin{array}{cc}\hfill {\mathfrak{n}}_{i,j}^{+}=& \Vert {\mathbf{P}}_i(f_{i-1,j}^{+}{\widehat{\mathbf{t}}}_{i-1}+f_{i,j}^{+}{\widehat{\mathbf{t}}}_i)\Vert ,\hfill \\ \hfill {\mathfrak{n}}_{i,j}^{-}=& \Vert {\mathbf{P}}_i(f_{i-1,j}^{-}{\widehat{\mathbf{t}}}_{i-1}+f_{i,j}^{-}{\widehat{\mathbf{t}}}_i)\Vert \hfill \end{array}$

(34)

where P _i is the matrix that projects vectors onto the plane of the disk $({\mathbf{P}}_i=\mathbf{I}-{\widehat{\mathbf{n}}}_i{\widehat{\mathbf{n}}}_i^{\text{T}})$ and ${\widehat{\mathbf{n}}}_i$ is a unit vector normal to the disk (3 ^rd column of ${}^0\mathbf{R}_t(s_{d,i})$ ). Also, ${\widehat{\mathbf{t}}}_i$ and ${\widehat{\mathbf{t}}}_{i-1}$ are unit vectors in the direction of the tendon segments after disks i and i − 1, respectively. These unit vectors can be found from the geometry of the disks and the local disk frames.

Given the tension $f_{i-1,j}^{+}$ or $f_{i-1,j}^{-}$ , the tension $f_i^{+}$ or $f_i^{-}$ can be calculated using a constrained minimization approach. The following optimization problem finds the minimum value the tension subject to the constraints that the wire stays in tension and the friction force is proportional to the normal force by the friction coefficient μ:

$\begin{array}{cc}\hfill \min & f_{i,j}^2\hfill \\ \hfill s.t.& 0\le f_{i,j},{\mathfrak{f}}_{i,j}=\mu {\mathfrak{n}}_{i,j}\hfill \end{array}$

(35)

This optimization problem must be solved starting from i = 1 (i.e., the segment base) and progressing at each section of tendon until the end-disk and back from the end-disk to the termination point. In our approach, the constrained minimization problem given in (35) is solved using the sequential quadratic programming algorithm implemented in MATLAB’s fmincon function. At each disk, after $f_i^{+}$ and $f_i^{-}$ is known, the friction force can be found using (33).

Once the friction force at each bushing is known, the total friction force experienced by capstan j is the sum of each value of ${\mathfrak{f}}_{i,j}^{+}$ and ${\mathfrak{f}}_{i,j}^{-}$ :

${\mathfrak{f}}_j=\sum _{i=1}^{2(n_d+1)}{\mathfrak{f}}_{i,j}^{+}+\sum _{i=1}^{2(n_d+1)}{\mathfrak{f}}_{i,j}^{-}$

(36)

The torque on the capstans due to friction can now be found using the capstan radius r _c :

${\mathit{\tau }}_F=\mathrm{sign}(\dot{\mathbf{q}})r_c{\left[\begin{array}{cc}\hfill {\mathfrak{f}}_1\hfill & \hfill {\mathfrak{f}}_2\hfill \end{array}\right]}^{\text{T}}$

(37)

For the purposes of simulation, the $\mathrm{sign}(\dot{\mathbf{q}})$ function was replaced with $\tanh (10\dot{\mathbf{q}})$ to not introduce discontinuities into the dynamic model which can cause issues with numerical differential equation solvers. Using J _qc , the generalized forces due to friction can be found from τ _F using:

${\mathit{\kappa }}_{\mathit{fric}}={\mathbf{J}}_{qc}^{\text{T}}{\mathit{\tau }}_F$

(38)

5.4. Equations of motion

Using the total kinetic energy of the segment from section 5.1 and the total potential energy of the segment from section 5.2, the Lagrangian $\mathcal{L}$ can be written as:

$\mathcal{L}=T-V$

(39)

Using $\mathcal{L}$ , the segment’s equations of motion can be found using the Euler–Lagrange equation:

$\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{\mathbf{c}}}-\frac{\partial \mathcal{L}}{\partial \mathbf{c}}=\mathit{\kappa }$

(40)

where κ are the non-conservative generalized forces:

$\mathit{\kappa }={\mathbf{J}}_{qc}^{\text{T}}\mathit{\tau }+{\mathbf{J}}_{\xi c}^{\text{T}}(s_c){\mathbf{w}}_c-{\mathit{\kappa }}_{\mathit{fric}}$

(41)

In this equation, τ ∈ IR² are capstan moments and w _c ∈ IR⁶ is the external wrench applied to the robot at arc length s _c . In this paper, we use the convention that a wrench consists of a force f ∈ IR³ proceeding a moment m ∈ IR³ (i.e., ${\mathbf{w}}_c={[{\mathbf{f}}_c^{\text{T}},{\mathbf{m}}_c^{\text{T}}]}^{\text{T}}$ ). Expanding (40) results in ² :

$\mathbf{M}\ddot{\mathbf{c}}+\mathbf{N}\dot{\mathbf{c}}+\frac{\partial V}{\partial \mathbf{c}}=\mathit{\kappa }$

(42)

where M is the generalized inertia matrix given in (25) and N is the centrifugal/Coriolis matrix. The derivation of N and ∂V/∂c is given in the appendix.

6. Modal-space momentum observer collision detection

6.1. Observer formulation

In this section, we adapt the momentum observer framework presented in De Luca and Mattone (2003) to our modal approach and present our implementation of the observer. The momentum observer approach relies on a detection of a change in the robot’s generalized momentum as a means to detect an incident of collision. The generalized momentum p is defined as:

$\mathbf{p}=\frac{\partial \mathcal{L}}{\partial \dot{\mathbf{c}}}=\mathbf{M}\dot{\mathbf{c}}$

(43)

Since the contact impulse is related to a temporal change in momentum, we calculate $\dot{\mathbf{p}}$ :

$\dot{\mathbf{p}}=\dot{\mathbf{M}}\dot{\mathbf{c}}+\mathbf{M}\ddot{\mathbf{c}}$

(44)

Substituting the solution for $\mathbf{M}\ddot{\mathbf{c}}$ from (42) results in:

$\dot{\mathbf{p}}=\dot{\mathbf{M}}\dot{\mathbf{c}}-\mathbf{N}\dot{\mathbf{c}}-\frac{\partial V}{\partial \mathbf{c}}+\mathit{\kappa }$

(45)

Since the matrix $\dot{\mathbf{M}}-2\mathbf{N}$ is skew-symmetric as shown in Koditschek (1984), we use $\dot{\mathbf{M}}-2\mathbf{N}=-{\left(\dot{\mathbf{M}}-2\mathbf{N}\right)}^{\text{T}}$ and solve for $\dot{\mathbf{M}}$ . This results in:

$\dot{\mathbf{M}}={\mathbf{N}}^{\text{T}}+\mathbf{N}$

(46)

Using this result along with (41), equation (45) can now be expressed as:

$\dot{\mathbf{p}}={\mathbf{N}}^{\text{T}}\dot{\mathbf{c}}-\frac{\partial V}{\partial \mathbf{c}}+{\mathbf{J}}_{qc}^{\text{T}}\mathit{\tau }-{\mathit{\kappa }}_{\mathit{fric}}+{\mathbf{J}}_{\xi c}^{\text{T}}(s_c){\mathbf{w}}_c$

(47)

Since we do not know the contact wrench w _c , we define a residual r representing the effect of an unknown contact impulse on the system. Since we plan to estimate r, the momentum temporal rate is rewritten as:

$\dot{\tilde{\mathbf{p}}}={\mathbf{N}}^{\text{T}}\dot{\mathbf{c}}-\frac{\partial V}{\partial \mathbf{c}}+{\mathbf{J}}_{qc}^{\text{T}}\mathit{\tau }-{\mathit{\kappa }}_{\mathit{fric}}+\mathbf{r}$

(48)

In the following derivation, we assume that a contact estimation algorithm is updated at every sensory acquisition cycle i, but the value of r is updated from the preceding cycle:

${\dot{\tilde{\mathbf{p}}}}_i={\left({\mathbf{N}}^{\text{T}}\dot{\mathbf{c}}-\frac{\partial V}{\partial \mathbf{c}}+{\mathbf{J}}_{qc}^{\text{T}}\mathit{\tau }-{\mathit{\kappa }}_{\mathit{fric}}\right)}_i+{\mathbf{r}}_{i-1}$

(49)

To reject sensory noise, we assume that r_i−1 includes estimate updates from the previous n sensory cycles. At any kth cycle (k ∈ [1, n _w ], n _w >> n), we define the momentum residual as:

${\mathbf{r}}_k={\mathbf{K}}_o\left[{\mathbf{p}}_k-{\mathbf{p}}_1-\sum _{j=k-n}^k{\dot{\tilde{\mathbf{p}}}}_j\text{d}t\right],k\in [1,n_w]$

(50)

In the above definition, dt is the sensory update cycle time, n _w denotes the total number of cycles where the estimator is allowed to run before it is reset. This is needed to avoid false positives due to cumulative summation/integration effects of measurement noise and model error. The matrix K _o is a diagonal positive definite gain matrix determining the dynamics of the momentum residual as discussed in section 6.2, p₁ is the initial generalized momentum at the beginning of a contact estimation span k = 1.

Substituting (49) in (50) results in the explicit form of the update law for the residual:

$\begin{array}{cc}\hfill & {\mathbf{r}}_k={\mathbf{K}}_o\left[{\mathbf{M}}_k{\dot{\mathbf{c}}}_k-{\mathbf{p}}_1-\right.\hfill \\ \hfill & \left.\sum _{j=k-n}^k\left({\left({\mathbf{N}}^{\text{T}}\dot{\mathbf{c}}-\frac{\partial V}{\partial \mathbf{c}}+{\mathbf{J}}_{qc}^{\text{T}}\mathit{\tau }-{\mathit{\kappa }}_{\mathit{fric}}\right)}_j+{\mathbf{r}}_{j-1}\right)\text{d}t\right]\hfill \end{array}$

(51)

6.2. Estimator dynamics and stability

In this section, we will repeat a result first introduced by De Luca and Mattone (2003) which proves that the GMO is a stable, first-order, filtered estimate of the generalized forces due to contact. To investigate the estimation dynamics, we will cast (50) into the continuous time domain and then take the derivative with respect to time:

$\dot{\mathbf{r}}={\mathbf{K}}_o\left[\dot{\mathbf{p}}(t)-\frac{d}{dt}{\int }_0^t\dot{\tilde{\mathbf{p}}}(\sigma )d\sigma \right]$

(52)

This equation simplifies to:

$\dot{\mathbf{r}}={\mathbf{K}}_o\left[\dot{\mathbf{p}}(t)-\dot{\tilde{\mathbf{p}}}(t)\right]$

(53)

Plugging (45) and (49) in for $\dot{\mathbf{p}}$ and $\dot{\tilde{\mathbf{p}}}$ , respectively, and simplifying yields:

$\dot{\mathbf{r}}(t)={\mathbf{K}}_o\left[{\mathit{\kappa }}_c(t)-\mathbf{r}(t)\right]$

(54)

where ${\mathit{\kappa }}_c(t)={\mathbf{J}}_{\xi c}^{\text{T}}(s_c){\mathbf{w}}_c$ . In the Laplace domain, the ith element of (54) can be written as:

$s{r}_i(s)=K_{o,i}\left[{\kappa }_{c,i}(s)-r_i(s)\right]$

(55)

Solving for κ_c,i yields:

$r_i(s)=\frac{K_{o,i}}{s+K_{o,i}}{\kappa }_{c,i}(s)$

(56)

This equation can be re-written to reveal the time constant of the filter:

$r_i(s)=\frac{1}{\frac{1}{K_{o,i}}s+1}{\kappa }_{c,i}(s)$

(57)

As shown in De Luca and Mattone (2003), for a positive estimation gain K_o,i > 0, the transfer function r _i (s)/κ_c,i(s) given in (57) has one stable (i.e., negative) pole at s = −K_o,i. Therefore, in ideal scenarios, the estimation signal r _i will approach the true generalized force κ_c,i with time constant 1/K_o,i. From this analysis, it can be seen that selecting a gain is a balancing act between responsiveness of the estimator and noise filtering. Selecting a higher K_o,i value will create a more responsive estimator, but the estimator will be more sensitive to noise. Conversely, a lower gain will more effectively filter noise, but will have a delay in the estimation.

7. Wrench estimation

The output of the momentum observer r is an estimate of the generalized forces due to contact. To determine the corresponding estimate of the applied wrench, ${\tilde{\mathbf{w}}}_c$ , the following equation must be solved for ${\tilde{\mathbf{w}}}_c$ :

$\mathbf{r}={\mathbf{J}}_{\xi c}^{\text{T}}(s_c){\tilde{\mathbf{w}}}_c$

(58)

In scenarios where the location of contact is known (e.g., from a sensor Abah et al. (2019; 2022)) and ${\mathbf{J}}_{\xi c}^{\text{T}}(s_c)\in {\text{IR}}^{6\times 6}$ is full rank, ${\tilde{\mathbf{w}}}_c$ can be found by inverting ${\mathbf{J}}_{\xi c}^{\text{T}}(s_c)$ .

7.1. Constrained minimization-based wrench estimation

However, the continuum segment shown in Figure 2 cannot instantaneously translate along the local $\widehat{\mathbf{z}}$ axis due to its nickel-titanium backbone. Therefore, for this robot, $rank\left[{\mathbf{J}}_{\xi c}(s_c)\right]\le 5$ . This has two implications (1) (58) can no longer be solved by simple inversion of ${\mathbf{J}}_{\xi c}^{\text{T}}(s_c)$ and (2) if ${\mathbf{w}}_c\in Null\left[{\mathbf{J}}_{\xi c}^{\text{T}}(s_c)\right]$ the wrench is fundamentally insensible. Therefore, we introduce a constrained minimization based approach to the wrench estimation problem:

$\begin{array}{cc}\hfill \underset{{\tilde{\mathbf{w}}}_c}{\min }& {\tilde{\mathbf{w}}}_c^{\text{T}}\mathbf{W}{\tilde{\mathbf{w}}}_c\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \mathbf{r}={\mathbf{J}}_{\xi c}^{\text{T}}(s_c){\tilde{\mathbf{w}}}_c\hfill \end{array}$

(59)

In this equation, W ∈ IR^6×6 is a diagonal weight matrix. This minimization problem will find the weighted minimum-norm solution to (58).

7.2. Incorporating knowledge of the contact type

In some scenarios, information about the type of external contact is known which imposes constraints on the components of the applied wrench. For instance, if the contact is assumed to be a point contact, it cannot impart a moment to the segment and the moment components of ${\tilde{\mathbf{w}}}_c$ can be assumed to be zero. For a comprehensive list of contact constraints, see Mason and Salisbury (1985). The approach presented here is a modification of the constrained optimization approach initially presented in Xu and Simaan (2008). Additionally, due to the assumption of infinite torsional rigidity and inextensibility of the central backbone, any applied moments/forces about/along the local z axis cannot alter the generalized momentum of the robot and will therefore be insensible by the momentum observer and can be assumed to be zero. These assumptions can be incorporated in the wrench estimation using the constraint matrix $\mathbf{A}\in {\text{IR}}^{n_c\times 6}$ where n _c is the number of components of ${\tilde{\mathbf{w}}}_c$ that can be assumed to be zero:

$\begin{array}{cc}\hfill \underset{{\tilde{\mathbf{w}}}_c}{\min }& {\tilde{\mathbf{w}}}_c^{\text{T}}\mathbf{W}{\tilde{\mathbf{w}}}_c\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \mathbf{r}={\mathbf{J}}_{\xi c}^{\text{T}}(s_c){\tilde{\mathbf{w}}}_c,\mathbf{A}{\tilde{\mathbf{w}}}_c=\mathbf{0}\hfill \end{array}$

(60)

For the simulations and experiments in this paper, the applied wrench is assumed to be a point contact and the force along the local z axis is assumed to be zero for the reasons described above. Therefore, we will use the following A matrix:

$\mathbf{A}=[{\mathbf{0}}_{4\times 2},{\mathbf{I}}_4]$

(61)

For the simulations and experiments in this paper, (60) was implemented in MATLAB 2023a and was compiled using MATLAB coder. During the experiments shown in Figure 9(a), (60) was solved with a median execution time of 0.0390 seconds and a maximum execution time of 0.1297 seconds. The computer used for this analysis has 32 GB of RAM and an Intel™ i7-13700 processor. The operating system was Windows 10.

8. Effect of state uncertainty on wrench estimation

In the following, we consider the effect of parameter and model uncertainty on the GMO method. Understanding how uncertainty in the segment’s dynamic state affects the output of GMO is critical to characterizing the expected performance of the method. We therefore investigate how an unmodeled error in the state vector will affect the GMO residual.

Since our formulation is in modal space, we define the segment state vector as:

$\mathbf{x}={[{\mathbf{c}}^{\text{T}},{\dot{\mathbf{c}}}^{\text{T}}]}^{\text{T}}\in {\text{IR}}^{12}$

(62)

The effect of an unmodeled change Δx on the ith element of r can be expressed using a Taylor series expansion:

$\begin{array}{ll}\hfill r_i(\mathbf{x}+\mathrm{\Delta }\mathbf{x})& =r_i(\mathbf{x})+\sum _{j=1}^{12}\frac{\partial r_i(\mathbf{x})}{\partial x_j}\mathrm{\Delta }{x}_j\hfill \\ \hfill & +\frac{1}{2}\sum _{k=1}^{12}\sum _{j=1}^{12}\frac{{\partial }^2{r}_i(\mathbf{x})}{\partial x_k\partial x_j}\mathrm{\Delta }{x}_k\mathrm{\Delta }{x}_j+\dots \hfill \end{array}$

(63)

In this work, we will approximate the effect of an unmodeled deviation in dynamic state Δx using a first-order approximation of the Taylor series:

$\mathbf{r}(\mathbf{x}+\mathrm{\Delta }\mathbf{x})\approx \mathbf{r}(\mathbf{x})+\frac{d\mathbf{r}(\mathbf{x})}{d\mathbf{x}}\mathrm{\Delta }\mathbf{x}$

(64)

The derivative of r with respect to x can be found by differentiating (50):

$\frac{d{\mathbf{r}}_k(\mathbf{x})}{d\mathbf{x}}={\mathbf{K}}_o\left[\frac{d{\mathbf{p}}_k}{d\mathbf{x}}-\sum _{j=k-n}^k\frac{d{\dot{\tilde{\mathbf{p}}}}_j}{d\mathbf{x}}\text{d}t\right]$

(65)

First, we will calculate dp _k /dx. The i^th element of p _k can be represented using the sum:

$p_i=\sum _{l=1}^6{M}_{il}{\dot{c}}_l$

(66)

Note that in this equation and the subsequent analysis, the subscript k, indicating the sample number, has been dropped for clarity. Using this summation representation of p, the ijth element of dp/dx can be written as:

$\frac{\partial p_i}{\partial x_j}=\sum _{l=1}^6\frac{\partial M_{il}}{\partial x_j}{\dot{c}}_l+M_{il}\frac{\partial {\dot{c}}_k}{\partial x_j}$

(67)

In our implementation, the partial derivatives of the elements of M are computed using a central finite differences approximation. The partial derivative of ${\dot{c}}_l$ with respect to x _j can be computed simply as:

$\frac{\partial {\dot{c}}_l}{\partial x_j}=\left\{\begin{array}{cc}0,\hfill & {\dot{c}}_l\ne x_j\hfill \\ 1,\hfill & {\dot{c}}_l=x_j\hfill \end{array}\right.$

(68)

Next, we will compute the derivative of ${\dot{\tilde{\mathbf{p}}}}_j$ with respect to the state vector by differentiating (49) with respect to x:

$\begin{array}{c}\hfill \frac{d{\dot{\tilde{\mathbf{p}}}}_j}{d\mathbf{x}}=\frac{d}{d\mathbf{x}}\left({\mathbf{N}}_j^{\text{T}}{\dot{\mathbf{c}}}_j\right)-\frac{d}{d\mathbf{x}}{\left[\frac{\partial V}{\partial \mathbf{c}}\right]}_j+\frac{d}{d\mathbf{x}}{\mathbf{J}}_{qc}^{\text{T}}\mathit{\tau }{0}_j\\ \hfill -\frac{d}{d\mathbf{x}}{\mathit{\kappa }}_{\mathit{fric,j}}+\frac{d}{d\mathbf{x}}{\mathbf{r}}_{j-1}\end{array}$

(69)

Because J _qc is not a function of c or $\dot{\mathbf{c}}$ , the generalized forces due to the capstan torques cancel to zero. The term d/dxr_j−1 is the value of d/dxr from the previous time step. In our implementation, the other derivatives in this equation are computed via central finite differences.

Given a characterized uncertainty in x, (64) can be used to predict the associated uncertainty in the GMO estimation residual r and subsequently the estimated applied wrench. This can be used to set thresholds on the elements of r where, if the value is less than the threshold, it cannot be determined if the residual is due to uncertainty in the dynamics or due to a small contact force and therefore should be ignored in practical applications. Conversely, if the value is above the threshold, the residual cannot be fully explained by dynamic model errors and thus the robot should be considered in a contact state. To establish these thresholds, we plug x = 0 into (64) and set Δx to the characterized uncertainty in the dynamic state. Because (64) is a first-order truncation of the Taylor series, it is inherently an approximation and there may be small differences between the predicted and actual observer residuals for a given error in x.

9. Joint force/torque deviation

One intuitive method for contact detection/estimation is the joint force/torque deviation (JFD) method. This method was proposed in some of the earliest works on robot contact estimation (e.g., Suita et al. (1995)) and has been applied to continuum robots for contact detection for static scenarios in Bajo and Simaan (2010). This method is also referred to in the literature as the direct estimation method (e.g., Haddadin et al. (2017)). The key idea of this method is to rearrange the robot’s full dynamic model (42) to directly solve for the generalized forces due to contact:

${\tilde{\mathit{\kappa }}}_c=\mathbf{M}\ddot{\mathbf{c}}+\mathbf{N}\dot{\mathbf{c}}+\frac{\partial V}{\partial \mathbf{c}}-{\mathbf{J}}_{qc}^{\text{T}}\mathit{\tau }+{\mathit{\kappa }}_{\mathit{fric}}$

(70)

Using ${\tilde{\mathit{\kappa }}}_c$ , the method described in (60) can be used to estimate the applied wrench during the contact. The major drawback of the JFD is the need to compute the second derivative of the modal coefficients. Unlike the noise filtering nature of the GMO, taking the numerical second derivative of c is a noise amplification process. Significant care needs to be taken to apply low-pass filters to the numerical derivatives of c and to the estimation signal ${\tilde{\mathit{\kappa }}}_c$ .

In the experimental evaluation section of this paper, the JFD method is used as a baseline of comparison for the performance of the GMO.

10. Dynamic model calibration

The dynamic model presented in section 5 contains numerous parameters (masses, inertias, etc.). Most of these parameters can be found from physical measurements, known from the robot geometry, or estimated from CAD software. These parameters are summarized in Table 3. Note that because of the assumption of zero torsional deflections implicit to the formulation of (4), the value used for JG does not affect the system dynamics. Some parameters, however, need to be calibrated from experimental data. These parameters were chosen to be the backbone flexural rigidities, EI _x and EI _y , and friction coefficients μ₁ and μ₂.

Table 3.

Dynamic model parameters.

L = 300.65 mm	r c = 15.255 mm	γ = 2.83 mm
ρ = 83.1532 g/m	r = 2 mm	n d = 6
JG = 1	J A = 14.323 gm2	g = [0, 0, −9.81]T m/s2
s d = [53.08, 102.62, 153.16, 203.70, 254.24, L] mm	md,i = 308.81 g, i ∈ [1, n d − 1]	md,i = 743.12 g, i = n d
$\begin{array}{cc}\hfill {\mathcal{I}}_{d,i}=& \left[\begin{array}{ccc}\hfill 0.5211\hfill & \hfill 0.0024\hfill & \hfill -0.0002\hfill \\ \hfill 0.0024\hfill & \hfill 0.5273\hfill & \hfill 0.0007\hfill \\ \hfill -0.0002\hfill & \hfill 0.0007\hfill & \hfill 0.9934\hfill \end{array}\right]g{m}^2\hfill \\ \hfill & i\in [1,n_d-1]\hfill \end{array}$		$\begin{array}{cc}\hfill {\mathcal{I}}_{d,i}& =\left[\begin{array}{ccc}\hfill 1.1580\hfill & \hfill -0.0357\hfill & \hfill 0.0001\hfill \\ \hfill -0.0357\hfill & \hfill 1.4871\hfill & \hfill 0\hfill \\ \hfill 0.0001\hfill & \hfill 0\hfill & \hfill 2.0564\hfill \end{array}\right]g{m}^2\hfill \\ \hfill & i=n_d\hfill \end{array}$
$\begin{array}{cc}\hfill {\mathbf{p}}_{cm,i}=& \left[\begin{array}{c}\hfill -0.0739\hfill \\ \hfill 0.1954\hfill \\ \hfill 5.7456\hfill \end{array}\right]mm,i\in [1,n_d-1]\hfill \end{array}$	$\begin{array}{cc}\hfill {\mathbf{p}}_{cm,i}=& \left[\begin{array}{c}\hfill -0.0133\hfill \\ \hfill 0\hfill \\ \hfill 20.8966\hfill \end{array}\right]mm,i=n_d\hfill \end{array}$	f pl = 208 N

To calibrate these parameters, the segment’s capstans were commanded to move with a chirp trajectory that starts from f₁ = 0.01 Hz and linearly increases to f₀ = 0.25 Hz in t _f = 6 seconds, has an amplitude of q_max = 360° (which corresponds to 42 degrees of segment bending), and a phase offset of ϕ₀ = 90°. The chirp trajectory is given by:

$q(t)=q_{\mathit{max}}\sin \left[2\pi \left(\left(\frac{f_1-f_0}{2t_f}\right)t^2+f_{0}t\right)+{\varphi }_0\right]$

(71)

During the experiment, the string encoders were read using a Teensy 4.1 microcontroller and sent to the central control computer using user datagram protocol (UDP) communication. For more information on the string encoders see Orekhov et al. (2020, 2023). The segment’s motors (Maxon™ DCX22L GB KL) were controlled using Maxon™ ESCON 50/5 motor drivers, Sensoray™ 526 encoder reading cards, and a Diamond Systems™ Aries PC/104 CPU running Ubuntu 14.04 with the PREEMPT-RT real-time patch. The segment motors were current controlled and the commanded current was used to estimate the capstan torque using the motor’s torque current constant, the rated efficiency of the motor and gearhead (Maxon™ GPX22HP 62:1), and the reduction ratios of the gearhead and the additional gearing (Figure 3③ and ④). The encoders and capstan torque readings are sent to the central control computer using UDP communication. The central control computer runs Ubuntu 18.04 with the Melodic distribution of the Robot Operating System (ROS). For the duration of the experiment, data was recorded from the motor encoders, motor torques, and string encoders using ROS topics at an average sampling rate of ∼100 Hz.

Using the recorded starting pose of the segment and the recorded capstan torques, the segment’s motion was simulated using MATLAB’s implementation of the Dormand–Prince version of the Runge-Kutta method (ode45) first described in Dormand and Prince (1980). Using the recorded modal coefficients from this experiment, the backbone flexural rigidities and friction coefficients were calibrated using the following constrained minimization problem:

$\begin{array}{cc}\hfill \mathit{\alpha }& =\arg \hspace{0.17em}\min \sum _{i=1}^n\frac{1}{2}{({\mathbf{c}}_i-{\tilde{\mathbf{c}}}_i)}^{\text{T}}({\mathbf{c}}_i-{\tilde{\mathbf{c}}}_i)\hfill \\ \hfill s.t.& \mathbf{0}<\mathit{\alpha }<{\mathit{\alpha }}_{\mathit{max}}\hfill \end{array}$

(72)

In this equation $\mathit{\alpha }={[E{I}_x,E{I}_y,{\mu }_1,{\mu }_2]}^{\text{T}}$ , c _i is the experimentally determined modal coefficients at sample i, and ${\tilde{\mathbf{c}}}_i$ is the simulated modal coefficients at sample i. Since we know the geometric parameters of the robot exactly, the uncertainty in the moments of inertia is very small (negligible). The primary source of uncertainty is due to the Young’s modulus of the combined superelastic NiTi and bellow backbone system. We used α _max = [2.5,2.5,1,1]^T during the calibration process after some experimental tuning. This minimization problem was solved using the interior point algorithm implemented in MATLAB’s fmincon function. The calibrated parameters are summarized in Table 4.

Table 4.

Calibrated dynamic model parameters.

EI x	EI y	μ 1	μ 2
1.1440	1.0373	0.0312	0.1637

Figure 5 compares the simulated and experimental x, y, and z components of the segment tip position ${}^0\mathbf{p}(L)$ using the calibrated parameters. Figure 5① shows the experimental tip x position and Figure 5② shows the simulated tip x position. The root-mean-square error (RMSE) in the x direction was 21.7 mm (7.2% of segment length). Figure 5③ shows the experimental tip y position and Figure 5④ shows the simulated tip y position. The RMSE in the y direction was 35.3 mm (11.7% of segment length). Figure 5⑤ shows the experimental tip z position and Figure 5⑥ shows the simulated tip z position. The RMSE in the z direction was 20.3 mm (6.8% of segment length).

Figure 5.

Simulated versus measured continuum segment tip positions after calibration: ① Experimental tip x position, ② simulated tip x position, ③ experimental tip y position, ④ simulated tip y position, ⑤ experimental tip z, and ⑥ simulated tip z position.

Other than the errors in the dynamic model, potential sources of error in this experiment include the use of commanded current to estimate the applied torque (the actual current supplied to the motor may differ from the command), the lack of exact data synchronization inherent in the decentralized structure of the control framework, errors due to numerical differentiation of the modal coefficients, backlash in the drive train, and sensor noise. Another major source of error was the shape sensing approach first presented in Orekhov et al. (2023). In this paper, the authors report an average end-disk positional error of 5.9 mm (1.9% of L) and a maximum error of 14.4 mm (4.8% of L).

11. Simulation studies

11.1. Effect of state-uncertainty

To investigate the effects of state uncertainty on the output of the momentum observer and the joint force/torque deviation method, we generated 21 body frame wrenches with x and y forces ranging between ±50 N. For each wrench, we also generated capstan torque values ranging between ±2 Nm. Using MATLAB’s implementation of the Dormand–Prince version of the Runge-Kutta method (ode45) (Dormand and Prince (1980)), the segment’s motion was simulated for 1 second for each randomly generated applied wrench/capstan torques. The simulations were started with the segment in the straight configuration and the applied wrench and capstan torques were increased from zero to the generated value over the 1 second simulation. The wrench was applied to the tip of the continuum robot (s _c = L). For each of the 21 generated wrenches, the applied wrench was estimated with constant errors in the state vector ranging from a 0% to 20% in steps of 2.22% for a total of 210 simulation conditions. For the GMO, the observer gain used in the simulation was K _o = 25I ∈ IR^6×6, the estimation weight matrix was the identity matrix W = I ∈ IR^6×6, and the estimation window was the entire 1 second trajectory. The constraint matrix used in these simulations is given in (61). For the JFD method, $\ddot{\mathbf{c}}$ was estimated using the numerical derivative of the output of dynamic simulation.

The results of these simulations are shown in Figure 6 and Table 5. Figure 6 shows polar scatter plots of the simulation results. Figure 6(a) shows the results for the GMO and Figure 6(b) shows the results for the JFD method. In both plots, the radial coordinate indicates the RMSE in the estimated force magnitude $w_e=\Vert {\mathbf{w}}_c-{\tilde{\mathbf{w}}}_c\Vert$ and the angular coordinate indicates the RMSE in the force direction ${\theta }_e=\arccos \left({\mathbf{w}}_c^{\text{T}}{\tilde{\mathbf{w}}}_c\right)/\left(\Vert {\mathbf{w}}_c\Vert \Vert {\tilde{\mathbf{w}}}_c\Vert \right)$ The color of the dots indicates the percent error in the state vector.

Figure 6.

Simulation results: (a) GMO and (b) JFD. In both plots, the radial axis shows the RMSE in the sensed force norm w _e , the angular axis shows the RMSE in the estimated wrench direction θ _e , and the color of the data points represent the percent error in the state vector.

Table 5.

Simulation results: Mean and max estimated force RMSE values at representative state vector percent errors.

Estimator	State error (%)	Mean(x)	Max(x)	Mean(y)	Max(y)
GMO	0	1.02 N	2.88 N	1.07 N	2.73 N
	4.44	1.50 N	3.48 N	1.28 N	3.14 N
	11.11	1.98 N	3.98 N	2.68 N	6.91 N
	15.56	2.52 N	6.77 N	3.95 N	11.60 N
	20	2.96 N	6.12 N	5.17 N	13.56 N
JFD	0	1.77 N	5.10 N	1.30 N	4.82 N
	4.44	4.66 N	16.11 N	4.10 N	20.29 N
	11.11	4.92 N	17.00 N	5.51 N	20.66 N
	15.56	5.57 N	17.55 N	6.64 N	21.22 N
	20	5.93 N	18.46 N	7.75 N	22.28 N

Table 5 shows the mean and max RMSE values at representative percent errors in the state vector. These results indicate that, assuming the applied forces are bounded to ±50 N and the error in the dynamic state is bounded to 20%, the average RMSE in the GMO force estimation is expected to be around 2.96 N in the x direction and 5.17 N in the y direction. The RMSE can be expected to be below 6.12 N in the x direction and 13.56 N in the y direction. For the JFD method, the average RMSE can be expected to be around 5.93 N in the x direction and 7.75 N in the y direction. The RMSE for the JFD method can be expected to be at or below 18.46 N in the x and 22.28 N in the y.

The largest reason for the difference between the errors in the GMO and JFD method is the need for a numerical second derivative of the modal coefficients which introduces additional errors into the system. In these simulations and the experiments in section 12, the estimation errors for forces in the x direction are lower than in the y direction. There are many sources that contribute to this difference. The biggest contributors are the difference between the calibrated flexural rigidities and friction coefficients in Table 4. Additionally, there is small difference in the elements of ${\mathcal{I}}_{d,i}$ as shown in Table 3.

11.2. Effect of sensor noise

As previously stated, a major advantage of the GMO is the ability to filter sensor noise. In this section, we present a simple simulation study showing how sensor noise affects the performance of the GMO and JFD methods. In this simulation study, we simulated the dynamics of the robot for a two second interval. During the first second, w _c was ramped up from zero to [10,−10,0,0,0,0]^T and then held steady for the rest of the simulation. The noiseless values of c, $\dot{\mathbf{c}}$ , and $\ddot{\mathbf{c}}$ were then plugged directly into the GMO and JFD estimators. In these simulations, the contact location was the robot tip (s _c = L), no joint torque was applied to the robot τ = [0,0]^T Nm, and the GMO estimation gain was set to K _O = 10I. No filtering was applied to the observer outputs or wrench estimates. The estimated wrenches for this noiseless simulation are given for the GMO and JFD in Figure 7(a) and (b), respectively.

Figure 7.

Effect of sensor noise simulation results: (a) GMO with no sensor noise. (b) JFD with no sensor noise. (c) GMO with 0.001 peak-to-peak amplitude random noise added to c. (d) JFD with 0.001 peak-to-peak amplitude random noise added to c. (e) GMO with 0.01 peak-to-peak amplitude random noise added to c. (f) JFD with 0.01 peak-to-peak amplitude random noise added to c.

Random noise with 0.001 peak-to-peak amplitude was introduced to c. Using this noisy version of c, $\dot{\mathbf{c}}$ and $\ddot{\mathbf{c}}$ were estimated using numerical differentiation. For this numerical differentiation, a 10-point Gaussian filter was used to smooth c and $\dot{\mathbf{c}}$ before differentiating. Figure 7(c) shows the result of the GMO estimation with this noisy state vector, and Figure 7(d) shows the results for the JFD. Lastly, this same process was repeated with 0.01 peak-to-peak amplitude random noise added to c. The results for the GMO and JFD are given in Figure 7(e) and (f), respectively. The RMSE for each noise level is given in Table 6.

Table 6.

Sensor noise simulation results for GMO and JFD at increasing noise levels.

Method	Direction	RMSE (N)
		Noiseless	0.001 Amp noise	0.01 Amp noise
GMO	x	0.67	1.35	3.20
	y	0.78	1.44	3.59
JFD	x	9.57 × 10−7	3.89	7.71
	y	1.26 × 10−6	5.37	9.43

As can be seen in Figure 7 and Table 6, when the state vector and dynamic parameters are known exactly, the JFD perfectly predicts the applied wrench to within numerical precision. On the other hand, the GMO has small differences in the estimation mostly due to the sensing delay and estimator dynamics as discussed in section 6.2. However, when noise is introduced into the simulation, the GMO becomes more accurate for two primary reasons: (1) as discussed in section 6.2 the GMO acts as a low-pass filter and (2) the GMO does not require the numerical second derivatives of c.

12. Experimental evaluation

In order to experimentally validate our method and investigate how the speed, location, and direction of contact affects the results, we mounted a Bota Systems™ Rokubi force/torque sensor (Figure 8①) in a known location relative to the continuum segment such that the positive z axis of the force/torque sensor aligned with the positive x axis of the segment base frame {T(0)}. The force/torque sensor was positioned such that the tip was barely touching the center of first spacer disk. The segment was then commanded to move away from the sensor with capstan velocities $\dot{\mathbf{q}}=[0,-180]°/s$ until q = [0, −180] °. The segment was then commanded to move towards the force/torque sensor with capstan velocities $\dot{\mathbf{q}}=[0,36]°/s$ until q = [0, 225] ° or the magnitude of the z axis force/torque sensor reading was $\ge 20$ N. Finally, the segment was commanded to move away from the sensor with capstan velocities $\dot{\mathbf{q}}=[0,-36]°/s$ . This was repeated for each disk and with the z axis of the sensor aligned with the positive y, negative y, positive x, and negative x axes of {T(0)} for a total of twenty experiments. The capstan angular velocities (and therefore contact point velocities) were varied among the test cases. Table 7 summarizes the conditions for each of the twenty test cases. The first row of the table shows the disk number contacted during the experiment, the second row shows the axis of {T(0)} that the z axis of the force/torque sensor was aligned with during the experiment, the third row shows contact arc lengths s _c , the fourth and fifth rows correspond with the speed of capstans 1 and 2 when the robot is moving towards the force/torque sensor, respectively. The sixth and seventh rows show the RMS errors of the estimated force in the x and y directions of the body frame at the contact point {T(s _c )}. The force/torque sensor was read using RS422 serial communication and published on a ROS topic. Videos of a few representative experiments are shown in Multimedia Extension I.

Figure 8.

Film strip showing an example of the experimental results given in Table 7: ① Bota Systems™ Rokubi force/torque sensor.

Table 7.

Experimental test cases and results. The highlighted column indicates the experiment that was plotted in Figure 9.

In the experiments, the modal coefficients c were estimated from the string-encoder and actuation tendon readings using (6). $\dot{\mathbf{c}}$ was estimated from c by applying a 30 point, backwards facing Gaussian filter to c and then taking a numerical derivative. Similarly, $\ddot{\mathbf{c}}$ was estimated from $\dot{\mathbf{c}}$ by applying a 30 point, backwards facing Gaussian filter to c and then taking a numerical derivative. For all experiments, the contact was assumed to be a point contact (i.e., using the A matrix given in (61)) and the estimation weight matrix was the identity matrix W = I ∈ IR^6×6.

For the GMO, the residual gain was K _o = diag([10, 25, 25, 10, 25, 0.02]). The value of K _o was chosen using trial and error. The lower gain values for the first, fourth, and last entries of K _o was selected due to higher levels of noise in their respective modal coefficients during the experiments. Because of noise introduced into the system from numerical differentiation of the modal coefficients, a 100 point, backwards-facing Gaussian filter was applied to both the estimation residual r and the estimated generalized forces ${\tilde{\mathbf{w}}}_c$ . The RMSE of the estimated forces in the disk frame x and y directions are given in Table 7. Figure 9 shows the estimated forces ${\tilde{\mathbf{w}}}_c$ and estimator residuals r for one of the twenty evaluation experiments. Figure 9(a) compares the estimated x and y forces in the local frame {T(s _c )} to the measured x and y forces in the {T(s _c )} frame. Figure 9(c) shows the estimated generalized forces r during the experiment.

Figure 9.

Plots showing an example of the experimental results given in Table 7. This is the same experiment shown in Figure 8. (a) The measured versus GMO-estimated applied force in the local disk frame. (b) The measured versus JFD estimated applied force in the local disk frame. (c) The GMO-estimated generalized forces r. (d) The JFD-estimated generalized forces ${\tilde{\mathit{\kappa }}}_c$ . (e) The absolute value of the error in the local x and y directions for the GMO. (f) The absolute value of the error in the local x and y directions for the JFD. ① Time span without contact, ② time span with contact, and ③ threshold on r above which the generalized force is assumed to be from external contact.

For the JFD, a 200 point, backwards-facing Gaussian filter was applied to both the estimated generalized forces ${\tilde{\mathit{\kappa }}}_c$ and the estimated wrench ${\tilde{\mathbf{w}}}_c$ . The RMSE of the estimated forces in the disk frame x and y directions are given in Table 7. Figure 9 shows the estimated forces ${\tilde{\mathbf{w}}}_c$ and estimator residuals ${\tilde{\mathit{\kappa }}}_c$ for one of the twenty evaluation experiments. Figure 9(b) compares the estimated x and y forces in the local frame {T(s _c )} to the measured x and y forces in the {T(s _c )} frame. Figure 9(d) shows the estimated generalized forces r during the experiment.

The shaded region of Figure 9(c) labeled ③ shows the threshold on the elements of r above which the generalized forces are assumed to be from external contact and not from uncertainty in the dynamic state as described in section 8. The value of Δx chosen to establish this threshold was 1.84 for all elements. This was 75% of the maximum error in the dynamic state during the calibration experiments described in section 10. The use of these thresholds helps guard against false positives, but comes at the cost of delayed detection response and sensitivity. Note that each element of r has its own threshold and Figure 9(c)③ only shows the threshold associated with r₄ for visual clarity. Figure 9(e) and (f) show the absolute value of the error in the x and y directions of the body frame at the contact point {T(s _c )} for the GMO and JFD, respectively.

Table 8 summarizes the mean and max RMSEs across all 20 experiments and for each intermediate disk. For the GMO, across all 20 experiments, the average RMSE was 2.59 N in the disk-frame x direction and 5.59 N in the disk-frame y direction. For the JFD method, across all 20 experiments, the average RMSE was 3.19 N in the disk-frame x direction and 5.57 N in the disk-frame y direction. When comparing the performance of the GMO and the JFD method, the GMO performed better on average than the JFD method in the x direction. The GMO and JFD method performed similarly on average in the y direction. It is worth noting that for the JFD method to have a similar performance to the GMO, the filter window on the estimation signals had to be twice as long as the GMO. This is because of the GMO only needs the numerical first derivative of c while the JFD method requires the numerical second derivative of c. The largest sources of error in these experiments were most likely the accuracy of the dynamic model and the numerical differentiation of the modal coefficients. Another major source of error is the accuracy of the shape sensing approach as detailed in Orekhov et al. (2023). Other sources of error include that the point of contact was assumed to be fixed relative to the spacer disks. However, during the experiments the contact point did move slightly relative to the disk as can be seen in Figure 8. Additionally, the segment’s NiTi backbone has a slight permanent bend in the yz-plane. This contributes to the larger errors in y direction forces. Neither method was able to detect contact on the first and second disks. This is because the forces were not able to make a significant enough change in the dynamic state to register as a contact. When comparing the simulations in Table 5 and the experimental result 8, it can be seen that the performance gain of the GMO over the JFD was less in experiment than in simulation. This is most likely due to the uncertainty in the dynamic parameters in real experimental scenarios.

Table 8.

Statistics that summarize the experimental results presented in Table 7.

Estimator	Statistic	Disk 1	Disk 2	Disk 3	Disk 4	Disk 5	All disks
GMO	Average RMSE x (N)	2.56	3.49	2.22	2.11	2.57	2.59
	Max. RMSE x (N)	4.62	6.63	3.97	4.37	4.68	6.63
	Average RMSE y (N)	6.89	5.17	6.06	5.51	4.31	5.59
	Max. RMSE y (N)	14.81	9.80	12.67	11.45	9.56	14.81
JFD	Average RMSE x (N)	2.55	2.95	1.80	4.54	4.10	3.19
	Max. RMSE x (N)	4.60	6.48	2.85	12.89	9.34	12.89
	Average RMSE y (N)	6.89	5.16	6.15	5.39	4.28	5.57
	Max. RMSE y (N)	14.80	9.60	12.35	11.44	9.57	14.80

The vast majority of prior works on continuum robot contact detection assume a constant curvature (i.e., circular) bending model. A key contribution of this paper is a method for robots that do not bend in constant curvature arcs. To demonstrate that the robot did not bend in constant curvature during the experiment shown in Figure 8, we introduce a simple measure of how non-circular the robot was bending at any given time. First we define, the maximum curvature along the robot’s arclength, u_max = max(u(s)), s = [0, L], and the minimum curvature along the robot’s arclength, u_min = min(u(s)), s = [0, L]. The measure β is simply:

${\beta }_x=u_{x,\max }-u_{x,\min },{\beta }_y=u_{y,\max }-u_{y,\min }$

(73)

The closer β is to zero, the better the robot’s bending can be approximated by a circular arc. Figure 10(a) shows the maximum curvature during the experiment shown in Figure 8. Figure 10(b) shows the circular-ness measure β during the experiment.

Figure 10.

(a) Maximum curvature in the local x and y direction during the experiment shown in Figure 8. (b) Circular-ness measure β in the local x and y direction during the experiment shown in Figure 8.

13. Preliminary experimental extension to multi-segment continuum robots

The results given in the previous section were for the distal segment in isolation. In this section, we will show preliminary results extending the method to the multisegment continuum robot shown in Figure 1. To do this, we must first discuss how the reaction moments from the distal segment’s actuation affects the dynamics of the proximal segment.

13.1. Reaction wrenches on the proximal segment

In these experiments, the mass and inertia of the distal segment was statically lumped into the end-disk. Additionally, the proximal segment must resist the reaction moments from the distal segment’s actuation motors as well as the reaction moment from the elbow motor placed between the two segments (see Figure 11④). To include these moments, we included the generalized forces due to the distal segment, κ _dist , into the non-conservative forces given in (41):

$\mathit{\kappa }={\mathit{\kappa }}_{\mathit{dist}}+{\mathbf{J}}_{qc}^{\text{T}}\mathit{\tau }+{\mathbf{J}}_{\xi c}^{\text{T}}(s_c){\mathbf{w}}_c-{\mathit{\kappa }}_{\mathit{fric}}$

(74)

Figure 11.

Multisegment experimental setup. (a) Robot at its starting configuration. (b) Close-up of the contact point. ① Proximal segment, ② distal segment, ③ Rokubi force/torque sensor, and ④ elbow motor between segments. Videos of these experiments are provided in multimedia extension II.

To calculate κ _dist , the wrench on the end-disk due to the elbow motor, w _elbow , and from the distal segment actuation motors, w_τ,dist are multiplied by the transpose of the proximal segment Jacobian:

${\mathit{\kappa }}_{\mathit{dist}}={\mathbf{J}}_{\xi c}^{\text{T}}(L)\left({\mathbf{w}}_{\mathit{elbow}}+{\mathbf{w}}_{\tau ,dist}\right)$

(75)

In these equations, ${\mathbf{w}}_{\mathit{elbow}}={[\mathrm{0,0,0},{\tau }_{\mathit{elbow}},\mathrm{0,0}]}^{\text{T}}$ and w_τ,dist = [0, 0, 0, τ _dist,1 , τ _dist,2 , 0].

13.2. Proximal segment calibration

Now that the effect of the distal segment’s actuation on the proximal segment is accounted for, we can calibrate the proximal segment’s dynamic parameters. To do this, a similar process to the one described in section 10 was completed. In the calibration experiment, the distal segment was attached to the proximal segment and the mass and inertia of the distal segment were lumped into the end-disk parameters of the proximal segment. The calibrated parameters are given in Table 9.

Table 9.

Proximal segment calibrated dynamic model parameters.

EI x	EI y	μ 1	μ 2
1.6003	2.1502	0.0900	0.2000

13.3. Experimental setup

In these experiments, we placed the continuum segments in a horizontal configuration as shown in Figure 11 (g = [0,9.81,0]^T m/s²). Again, the Bota Systems Rokubi force/torque sensor was placed in the workspace of the robot. In this configuration, similar to the single-segment experiments, both segments were commanded to move away from the sensor with capstan velocities $\dot{\mathbf{q}}=[0,-180]°/s$ until q = [0, −180] °. The segments were then commanded to move towards the force/torque sensor with capstan velocities $\dot{\mathbf{q}}=[0,36]°/s$ until q = [0, 225] ° or the magnitude of the z axis force/torque sensor reading was $\ge 30$ N. Finally, the segments were commanded to move away from the sensor with capstan velocities $\dot{\mathbf{q}}=[0,-36]°/s$ . Videos of these experiments are provided in multimedia extension II.

In the experiments, the robot came into contact with the 4th spacer disk of the distal segment s _c = 213.293 mm. The observer gains and filter windows were the same as the single-segment experiment.

13.4. Results

Figure 12 shows the experimental results for the proximal continuum segment and Figure 13 shows the results for the distal segment. In both figures, (a) shows the estimated wrench using the GMO and (b) shows the estimated wrench on the proximal segment using the JFD. Subfigures (c) and (d) show the estimated generalized forces due to contact for the GMO and JFD methods, respectively. Lastly, (e) shows the maximum curvature during the experiment and (f) shows the circular-ness measure during the experiment. Table 10 shows the estimation RMSE for the distal and proximal segments using both the JFD and GMO methods.

Figure 12.

Experimental results for the proximal segment. (a) Estimated forces in the x and y directions using the GMO. (b) x and y directions using the JFD. (c) GMO-estimated generalized forces r. (d) The JFD-estimated generalized forces ${\tilde{\mathit{\kappa }}}_c$ . (e) The maximum curvature in the local x and y directions. (f) The circular-ness measure β in the local x and y directions.

Figure 13.

Experimental results for the distal segment. (a) Estimated forces in the x and y directions using the GMO. (b) x and y directions using the JFD. (c) GMO-estimated generalized forces r. (d) The JFD-estimated generalized forces ${\tilde{\mathit{\kappa }}}_c$ . (e) The maximum curvature in the local x and y directions. (f) The circular-ness measure β in the local x and y directions.

Table 10.

RMSE values for the proximal and distal segments.

Segment	Method	Direction	RMSE (N)
Proximal	GMO	x	16.07
		y	10.51
	JFD	x	68.01
		y	10.74
Distal	GMO	x	10.45
		y	4.65
	JFD	x	15.96
		y	6.55

13.5. Discussion and limitations

As can be seen in Figure 11, the proximal segment is bent in a noticeable S-curve shape. This shape could not be reasonable approximated using a constant curvature assumption. In these experiments, both methods performed worse during the multisegment experiments than the single segment experiments. This is mostly because the segments GMO and JFD estimators were mostly run in isolation with minimal coupling effects. For example, the kinetic energy due to the motion in the proximal segment was not taken into account for the distal segments and the reaction forces due to motion of the distal segments were not taken into account for the proximal segment. Additionally, the distal segment was approximated as a fixed mass and inertia and lumped into the end-disk parameters. This approximation was a major source of error.

Also, although there was not contact with the proximal segment, the GMO detected contact as can be seen in Figure 12(a). This is due to the contact on the distal segment. The GMO was therefore estimating an equivalent wrench at its fourth disk that would cause the same change in dynamics as the true contact with the distal segment. In practice, the contact sensors in the robot’s sensor disk (see Abah et al. (2022)) could be used to determine which portion of the robot was in contact and the GMO could be used to determine the location of contact.

14. Conclusions

To safely operate alongside workers in a confined space, in-situ collaborative robots (ISCRs) must be endowed with a host of sensory information not typically present in traditional industrial robots. One of the most important of these sensory modalities is the ability to detect and estimate wrenches applied to the robot by the environment or by the human operator.

In this paper, we presented a formulation of the generalized momentum observer (GMO) in the space of modal coefficients of curvature. This formulation enables contact estimation for continuum robots with variable curvature. We presented a modal-space dynamic model for the continuum segment shown in Figure 8 and calibrated the central backbone flexural-rigidities and static friction coefficients that minimize errors in the dynamic model. Additionally, we presented a constrained minimization approach for estimating the wrench applied to the segment during the contact. This approach allows for incorporating knowledge of the contact type. The effect of dynamic state uncertainty on the estimation signal was then investigated and recommendations were made about establishing contact detection thresholds. A simulation study was presented that establishes limits on the expected performance of the method for up to 20% error in the segment’s dynamic state and compared the performance of the GMO to the joint force/torque deviation (JFD) method. The GMO and JFD methods were also compared using 20 different experimental conditions. In these experiments, the JFD was found to have similar performance in the local y direction (which had worse performance in general), but the GMO was found to have a better performance in the local x direction. We also demonstrated a basic extension of the method for multisegment continuum robots and showed preliminary experimental results.

A major limitation of this work is that the contact location cannot be determined. In order to achieve this, this method could be used in conjunction with the screw-deviation method presented in Bajo and Simaan (2012). This would have the added benefit of increasing sensitivity of the method at very low speeds. Another limitation of our work is that it does not work effectively for the first and second spacer disks.

Future work will include investigating contact localization, improving the dynamic model by including hysteretic elastic effects, and viscous friction effects. Although our method is the first to be specifically tailored to general curvature robots with significant inertial effects, we will work on a systematic comparison to additional existing contact detection methods to establish a baseline of comparison. Additionally, we will work on real-time non MATLAB-based implementation of this method to enable real-time contact estimation. Lastly, we will also investigate how the choice of Chebyshev polynomial order affects the performance of the method.

Dynamic matrices derivation

In this appendix, we expand upon the derivations of the terms in the dynamic model described in section 5. First, we will expand (40):

$\frac{d}{dt}\frac{\partial T}{\partial \dot{\mathbf{c}}}-\frac{d}{dt}\frac{\partial V}{\partial \dot{\mathbf{c}}}0-\frac{\partial T}{\partial \mathbf{c}}+\frac{\partial V}{\partial \mathbf{c}}=\mathit{\kappa }$

(76)

The potential energy does not depend $\dot{\mathbf{c}}$ so therefore the term $d/dt\partial V/\partial \dot{\mathbf{c}}$ goes to 0. The partial derivative of the kinetic energy can be written as:

$\frac{\partial T}{\partial \dot{\mathbf{c}}}=\mathbf{M}\dot{\mathbf{c}}$

(77)

The time derivative of $\partial T/\partial \dot{\mathbf{c}}$ can be written as:

$\frac{d}{dt}\frac{\partial T}{\partial \dot{\mathbf{c}}}=\dot{\mathbf{M}}\dot{\mathbf{c}}+\mathbf{M}\ddot{\mathbf{c}}$

(78)

The ith row of (78) can be written as:

$\frac{d}{dt}\frac{\partial T}{\partial {\dot{c}}_i}=\sum _{j=1}^6{M}_{ij}\ddot{c_j}+\sum _{j=1}^6\sum _{k=1}^6\frac{\partial M_{ij}}{\partial c_k}{\dot{c}}_j{\dot{c}}_k$

(79)

The partial derivative of M with respect to the ith modal coefficient can be found by differentiating the mass matrix as defined in (25):

$\frac{\partial \mathbf{M}}{\partial c_i}=\frac{\partial {\mathbf{M}}_{CB}}{\partial c_i}+\frac{\partial {\mathbf{M}}_{SD}}{\partial c_i}+\frac{\partial {\mathbf{M}}_A}{\partial c_i}0$

(80)

Because M _A is a constant matrix and not a function of the modal coefficients, ∂M _A /∂c _i goes to 0. In our implementation, the partial derivatives of M _CB and M _SD are found using a central finite difference approximation.

The third term in (40), the partial derivative of the kinetic energy with respect to the ith modal coefficient, can be written as:

$\frac{\partial T}{\partial c_i}=\sum _{j=1}^6\sum _{k=1}^6\frac{1}{2}\frac{\partial M_{jk}}{\partial c_i}{\dot{c}}_j{\dot{c}}_k$

(81)

Next, we will find the partial derivative of the potential energy V with respect to the modal coefficients:

$\frac{\partial V}{\partial \mathbf{c}}=\frac{d{V}_B}{d\mathbf{c}}+\frac{d{V}_{CB}}{d\mathbf{c}}+\frac{d{V}_{SD}}{d\mathbf{c}}$

(82)

The derivative of the bending potential energy:

$\frac{d{V}_B}{d\mathbf{c}}=\frac{1}{2}\left({\mathbf{K}}_c+{\mathbf{K}}_c^{\text{T}}\right)\mathbf{c}$

(83)

The derivative of the gravitational potential energy of the central backbone:

$\frac{d{V}_{CB}}{d\mathbf{c}}=-\rho {\int }_0^L{\left(\frac{\partial \mathbf{p}(s)}{\partial \mathbf{c}}\right)}^{\text{T}}\mathbf{g}ds$

(84)

The derivative of the gravitational potential energy of the spacer disks:

$\frac{d{V}_{SD}}{d\mathbf{c}}=-\sum _{i=1}^{n_d}{m}_{d_i}{(\frac{d\mathbf{p}(s_{d_i})}{d\mathbf{c}}+\frac{d}{d\mathbf{c}}({}^0\mathbf{R}_t(s_{d_i})){}^t\mathbf{p}_{c{m}_i})}^{\text{T}}\mathbf{g}$

(85)

Combining the above terms, the ith row of the segment dynamics can be written as:

$\sum _{j=1}^6{M}_{ij}\ddot{c_j}+\sum _{j=1}^6\sum _{k=1}^6\left(\frac{\partial M_{ij}}{\partial c_k}-\frac{1}{2}\frac{\partial M_{jk}}{\partial c_i}\right){\dot{c}}_j{\dot{c}}_k+\frac{\partial V}{\partial c_i}={\kappa }_i$

(86)

Using Christoffel symbols, the second term in the above equation can be written in terms of the centrifugal/Coriolis matrix, N. The ijth element of N is:

$N_{ij}=\frac{1}{2}\sum _{k=1}^6\left(\frac{\partial M_{ij}}{\partial c_k}+\frac{\partial M_{ik}}{\partial c_j}-\frac{\partial M_{kj}}{\partial c_i}\right){\dot{c}}_k$

(87)

Supplemental Material

Supplemental Material