Simplified Kinematics of Continuum Robot Equilibrium Modulation via Moment Coupling Effects and Model Calibration

TL;DR

This paper introduces a simplified kinematic model for continuum robots that captures micro and macro motion behaviors through moment coupling effects and calibration, validated by experiments with high accuracy.

Contribution

It presents a novel, calibration-based modeling approach that explains micro-scale motion phenomena in continuum robots, including turning point behavior.

Findings

Model accurately predicts micro and macro motions with RMS error of 5.82 micrometers.

Calibration parameters effectively reproduce turning point behavior.

Model fits experimental data well, especially when excluding motions past the turning point.

Abstract

Recently, a new concept for continuum robots capable of producing macro-scale and micro-scale motion has been presented. These robots achieve their multi-scale motion capabilities by coupling direct-actuation of push-pull back-bones for macro motion with indirect actuation whereby the equilibrium pose is altered to achieve micro-scale motion. This paper presents a first attempt at explaining the micro-motion capabilities of these robots from a modeling perspective. This paper presents the macro and micro motion kinematics of a single segment continuum robot by using statics coupling effects among its sub-segments. Experimental observations of the micro-scale motion demonstrate a turning point behavior which could not be explained well using the current modeling methods. We present a simplistic modeling approach that introduces two calibration parameters to calibrate the moment coupling effects among the sub segments of the robot. It is shown that these two parameters can reproduce the turning point behavior at the micro-scale. The instantaneous macro and micro scale kinematics Jacobians and the calibration parameters identification Jacobian are derived. The modeling approach is verified against experimental data showing that our simplistic modeling approach can capture the experimental motion data with RMS position error of 5.82 micrometers if one wishes to fit the entire motion profile with the turning point. If one chooses to exclude motions past the turning point, our model can fit the experimental data with an accuracy of 4.76 micrometers.

Tables

Table 1. Table 1: Nomenclature for Kinematic Modeling

Symbol	Description
Frame {F}	designates a right-handed frame with unit vectors ${\widehat{𝐱}}_f,{\widehat{𝐲}}_f,{\widehat{𝐳}}_f$ and point $𝐟$ as its origin.
\hdashline Frame {B}	the base disk frame with $𝐛$ located at the center of the base disk, ${\widehat{𝐱}}_b$ passing through the first secondary backbone and ${\widehat{𝐳}}_b$ perpendicular to the base disk.
\hdashline Frame {1}	Frame of the bending plane having ${\widehat{𝐳}}_1={\widehat{𝐳}}_b$ and ${\widehat{𝐱}}_1$ passing through with the project point of the end disk center. The angle $\delta$ is defined as from ${\widehat{𝐱}}_1$ to ${\widehat{𝐱}}_b$ about ${\widehat{𝐳}}_b$ according to right hand rule.
\hdashline Frames {E} & {G}	Frame {E} is defined with ${\widehat{𝐳}}_e$ as the normal to the end disk and ${\widehat{𝐱}}_e$ is the intersection of the bending plane and the end disk top surface. Frame {G} is the gripper frame that has the same $\widehat{𝐳}$ as {E}, i.e. ${\widehat{𝐳}}_g={\widehat{𝐳}}_e$, but ${\widehat{𝐱}}_g$ is passing through the first secondary backbone. It can be obtained by a rotation angle $\left(-{\sigma }_{1e}\right)$ about ${\widehat{𝐳}}_e$.
\hdashline Frames {P} & {C}	These frames are defined in a manner similar to the definition of frames {E} and {G} but for a specific arc insertion length $q_s$ as opposed to the full length of the robot segment $L$. The ${\widehat{𝐱}}_c-{\widehat{𝐲}}_c$ plane is the insertion plane as in shown in the planar case in Fig. 2.
\hdashline Frame {I}	designates the microscope image frame having the origin at the corner of the image and having its $x$-$y$ axes aligned with the width and height directions (Fig. 7(a, c)).
\hdashline Frame {M}	designates the marker frame that is determined by segmentation of the microscope image (Fig. 7(c)).
\hdashline Vector ${}_{}{}^x𝐩_{a/b}^{}$	designates the position of point $a$ relative to point $b$ that is expressed in frame {X}.

Table 2. Table 2: Robot Parameters Used in Simulations and Experiments

L	r	$E_p,E_i,E_s$	$d_p,d_i$	$d_s$	$I_p,I_i$	$I_s$
44.3mm	3mm	41 GPa	0.90 mm	0.38 mm	0.0312 ${\text{mm}}^4$	0.0010 ${\text{mm}}^4$

Table 3. Table 3: Calibration using full micro motion trajectory

Iteration	0	5	10	20	30	45	46
$\frac{k_{{\lambda }_0}}{100}$	0	4.22	6.72	9.06	9.88	10.23	10.24
$\frac{k_{{\lambda }_q}}{1000}$	0	2.7	4.3	5.7	6.3	6.5	6.5
RMSE [$\mu$m]	42.27	15.71	25.40	7.72	6.07	5.82	5.82

Table 4. Table 4: Calibration using partial micro motion trajectory

Iteration	0	5	10	20	30	58	59
$\frac{k_{{\lambda }_0}}{100}$	0	4.22	6.72	9.06	9.88	10.30	10.30
$\frac{k_{{\lambda }_q}}{1000}$	0	2.7	4.3	5.7	6.3	6.5	6.5
RMSE [$\mu$m]	19.51	13.01	9.30	6.12	5.17	4.76	4.76

Equations

$$m_1=E_p{I}_p\frac{\theta -{\theta }_0}{L}+{\sum }_i{E}_i{I}_i\frac{\theta -{\theta }_0}{L_i}$$

$${m_1}^{\prime }=E_p{I}_p\frac{{\theta }^{\prime }-{\theta }_s}{L-q_s}+{\sum }_i{E}_i{I}_i\frac{{\theta }^{\prime }-{\theta }_s}{L_{{\epsilon }_i}}$$

$${\Delta }_i=r\cos ({\sigma }_i),{\sigma }_i=\delta +(i-1)\beta$$

$$L_i=L+{\Delta }_i(\theta -{\theta }_0)$$

$$L_{s_i}=q_s+{\Delta }_i({\theta }_s-{\theta }_0)$$

$$L_{{\epsilon }_i}=(L-q_s)+{\Delta }_i({\theta }^{\prime }-{\theta }_s)=L_i-L_{s_i}$$

$$m_1={m_1}^{\prime }$$

$$m_2=$$

$$m_s=$$

$${m_1}^{\prime }+m_2+m_s=0$$

$$m_1=k_{{\theta }_0}(\theta -{\theta }_0),$$

$${m_1}^{\prime }=k_{{\theta }_1}({\theta }^{\prime }-{\theta }_s),$$

$$m_2=-k_{{\theta }_2}({\theta }_s-{\theta }_0),$$

$$m_s=-k_{{\theta }_s}({\theta }_s-{\theta }_0),$$

$${\theta }_s=\frac{k_{{\theta }_1}{\theta }^{\prime }+k_{{\theta }_2}{\theta }_0+k_{{\theta }_s}{\theta }_0}{k_{{\theta }_1}+k_{{\theta }_2}+k_{{\theta }_s}}$$

$${\theta }^{\prime }=\frac{k_0}{k_1}(\theta -{\theta }_0)+{\theta }_s=f_{{\theta }^{\prime }}({\theta }_s)$$

$${m_1}^{\prime }+m_2+m_s=\lambda (q_s,\theta ,{\mathbf{k}}_{\lambda })$$

$$\lambda (q_s,\theta ,{\mathbf{k}}_{\lambda })=k_{{\lambda }_0}+k_{{\lambda }_{\theta }}\theta +k_{{\lambda }_q}{q}_s;{\mathbf{k}}_{\lambda }\triangleq [k_{{\lambda }_0},k_{{\lambda }_q},k_{{\lambda }_{\theta }}{]}^{\mathrm{T}}$$

$${\theta }_s=\frac{k_1{\theta }^{\prime }+k_2{\theta }_0+k_s{\theta }_0-\lambda }{k_1+k_2+k_s}$$

$$\mathbf{\varphi }\triangleq [{\theta }_s,{\theta }_{\epsilon }{]}^{\mathrm{T}},{\theta }_{\epsilon }\triangleq {\theta }^{\prime }+\left(\frac{\pi }{2}-{\theta }_s\right)$$

$$\mathbf{\psi }\triangleq [\theta ,\delta {]}^{\mathrm{T}}$$

$$\mathbf{\varphi }={\mathcal{F}}_{\text{eqm}}(\mathbf{\psi },q_s,{\mathbf{k}}_{\lambda }),\mathbf{\varphi }\in I{R}^2,\mathbf{\psi }\in I{R}^2,{\mathbf{k}}_{\lambda }\in I{R}^3$$

$${\delta }_s={\delta }_{\epsilon }=\delta$$

$${}^{\text{base}}{\mathbf{p}}_{\text{end-disk / base}}=\frac{L_x}{{\theta }_x-\pi /2}\left[\begin{array}{c}\cos {\delta }_x(\sin {\theta }_x-1)\\ -\sin {\delta }_x(\sin {\theta }_x-1)\\ -\cos {\theta }_x\end{array}\right]$$

$${}^{\text{base}}{\mathbf{R}}_{\text{end-disk}}=e^{-{\delta }_x\left[{\mathbf{z}}^{\wedge }\right]}{e}^{(\frac{\pi }{2}-{\theta }_x)\left[{\mathbf{y}}^{\wedge }\right]}{e}^{{\delta }_x\left[{\mathbf{z}}^{\wedge }\right]}$$

$${}^b{\mathbf{p}}_{g/b}={}^b{\mathbf{p}}_{c/b}+{}^b{\mathbf{R}}_c{}^c{\mathbf{p}}_{g/c}\triangleq {\mathcal{F}}_g(\mathbf{\varphi },\delta ,q_s)$$

$${}^b{\mathbf{R}}_g={}^b{\mathbf{R}}_c{}^c{\mathbf{R}}_g={\mathbf{R}}_g(\mathbf{\varphi },\delta )$$

$${}^{b}\mathbf{T}_{g}=\left[\begin{array}[]{cc}{}^{b}\mathbf{R}_{g}&{}^{b}\mathbf{p}_{g/b}\\ \mathbf{0}&1\end{array}\right]\triangleq\boldsymbol{\mathscr{F}}_{T}(\bm{\phi},\delta,q_{s})$$

$${}^b{\mathbf{T}}_g\triangleq {\mathcal{F}}_T(\mathbf{\psi },q_s,{\mathbf{k}}_{\lambda })$$

$$\text{d}\mathbf{\xi }\triangleq {\left[(\text{d}\mathbf{x}{)}^{\mathrm{T}},(\text{d}\mathbf{\mu }{)}^{\mathrm{T}}\right]}^{\mathrm{T}},\text{d}\mathbf{\xi }\in I{R}^{6\times 1}$$

$$\text{d}\mathbf{x}\triangleq \text{d}({}^b{\mathbf{p}}_{g/b}),\text{d}\mathbf{\mu }\triangleq [\text{d}{\mu }_x,\text{d}{\mu }_y,\text{d}{\mu }_z{]}^{\mathrm{T}}$$

$$\text{d}\mathbf{\xi }=\frac{\partial \mathbf{\xi }}{\partial \mathbf{\varphi }}\text{d}\mathbf{\varphi }+\frac{\partial \mathbf{\xi }}{\partial \delta }\text{d}\delta +\frac{\partial \mathbf{\xi }}{\partial q_s}\text{d}{q}_s\text{vspace}1mm$$

$$\frac{\partial \mathbf{\xi }}{\partial \mathbf{\varphi }}\triangleq {\mathbf{J}}_{\mathbf{\xi }\mathbf{\varphi }}\in I{R}^{6\times 2},\frac{\partial \mathbf{\xi }}{\partial \delta }\triangleq {\mathbf{J}}_{\mathbf{\xi }\delta }\in I{R}^6,\frac{\partial \mathbf{\xi }}{\partial q_s}\triangleq {\mathbf{J}}_{\mathbf{\xi }{q}_s}\in I{R}^6$$

$$\mathbf{J}_{\bm{\xi}q_{s}}=\left[\begin{array}[]{c}\dfrac{\partial\;{}^{b}\mathbf{p}_{c/b}}{\partial q_{s}}+{}^{b}\mathbf{R}_{c}\;\dfrac{\partial\;{}^{c}\mathbf{p}_{g/c}}{\partial q_{s}}\vspace{1mm}\\ \hdashline[2pt/2pt]\mathbf{0}_{3\times 1}\end{array}\right]$$

Full text

1 Introduction

Current robotic manipulators for minimally invasive surgery (MIS) are capable of dexterous motion for surgical tasks requiring large workspace and position accuracy ranging from 0.5 to 1.5 mm. For example, the root mean square (RMS) localization accuracy of the da-Vinci Classic and da-Vinci S was evaluated experimentally as 1.02 mm and 1.05 mm respectively by Kwartowitz et al. [1, 2]. Despite recent increases in precision, current commercial surgical systems are unable to support micro-surgical precision (less than 0.1mm precision), and such precision can benefit micro-surgical tasks (e.g. micro-anastomosis and micro-vascular reconstruction [3, 4, 5, 6]).

This paper is motivated by a need for increased motion resolution at a micro-surgical scale and during deep surgical access minimally invasive surgery. In addition, the paper is equally motivated by the potential benefits of a new class of surgical devices capable of multi-scale motion. Such devices promise to provide a large workspace for traversal of deep passageways and for gross surgical manipulation while offering micro-scale motion suitable for cellular-level surgery. We refer to devices capable of macro and micro-scale motion as Multi Scale Motion (MSM) devices. With the advent of new devices with integrated optical coherence tomography imaging and confocal endo-microscopy (e.g. [7, 8]), the use of MSM can allow future image-based biopsy with imaging resolutions at cellular size [9, 10]. Such devices can in the future support surgical decisions on continued tumor excision, which can minimize the need for repeat follow-up surgeries due to incomplete resection of tumors.

To achieve MSM capabilities, this paper adopts the new design concept for Continuum Robot Equilibrium Modulation (CREM), which was first presented in [11]. CREM robots use a continuum structure that is primarily based on flexible elements to achieve large scale manipulation (i.e. robots without hinges and pin joints [12]). They also use fine adjustments to their static equilibrium pose in order to achieve micro-scale motion. The design concept for these robots is presented in Figure 1. This design is modified from a multi-backbone continuum robot presented in [13]. Each segment of a multi-segment continuum robot (MBCR) includes superelastic NiTi backbones. A single central backbone is surrounded by secondary backbones that are radially constrained by spacer disks and equidistantly distributed circumferentially around the central backbone. Macro motion of the robot tip is achieved by pushing and pulling on the secondary backbones (designated by the thick arrows in Fig. 1), which causes a deformation of the continuum robot body. We call this method of actuation direct actuation where the robot actuators directly change the length of the secondary backbones. In addition, CREM robots use indirect actuation whereby the equilibrium pose of the end effector is indirectly altered through a change of internal force distribution or by a change in material distribution altering cross sectional stiffness. For example, by inserting superelastic Ni-Ti beams (henceforth referred to as the Equilibrium Modulation Backbones (EMBs)) inside the secondary backbones (see thin arrows in Fig. 1), the static equilibrium of the robot is altered (modulated) by minute amounts.

Compared to prior designs, CREM robots possess a unique capability to allow MSM using a single design. Most prior works in the area of MSM rely on serial stacking of a macro-scale motion robot and a micro-scale robot, and such examples include Egeland’s pioneering work [14] and followed by several other works such as [15, 16, 17, 18, 19, 20, 21]. Other researchers investigated a variety of actuation methods and mechanisms to achieve micro motion capabilities, including a Steward/Gough parallel robot driven by hydraulic micro-actuators [22], twisted wire actuators for a planar parallel robot [23], a micromanipulation tool using shape memory alloy [24], and a piezo-electrically actuated parallel platform [25]. Although these works achieved micro-scale motion resolution, they are not suitable for MSM in confined spaces, in which continuum robots in general have advantages.

Within the context of continuum robots, the most relevant modeling works are [26] where a solution framework based on constraints of geometric compatibility and static equilibrium was derived using elliptic integrals for multi-backbone continuum robots and [27, 28] where Cosserat rod theory was used for dynamics modeling of wire-actuated continuum robots. One could use these methods to model the tip micro motion, however, due to the formulation complexity and the solutions of equilibrium direct kinematics based on energy minimization or boundary value problem solution it is hard to obtain an updated differential kinematics model that accounts for the exact bending shape curvatures. Also, due to uncertainties in material properties and friction, using an exact modeling method does not add value since a model calibration step has to be carried out anyhow when attempting to control a physical robot.

Another work is Li et al. [29], where the authors presented a constrained wire-driven flexible mechanism which used a constraint rod to selectively adjust its workspace by altering the length of its distal bending portion. The work showed that the constraint rod can change the workspace. The design however does not lend itself to easily allowing MSM and the work did not consider methods for achieving or modeling CREM.

Finally, our proposed design differs substantially from concentric tube robots [30, 31] in that the equilibrium modulation that creates the tip micro-motion is still governed by the strong geometric constraints employed by the parallel-backbone structure. Concentric tube continuum robots achieve their workspace through antagonistic bending of tube pairs and therefore they can theoretically be used for micro-scale equilibrium modulation. However to achieve micro-scale motion the designers are forced to use stiffness matched tube pairs with a very small difference in free curvature. The attainment of micro-scale motion by concentric tube robots therefore comes at the expense of sacrificing the macro-scale motion capabilities.

In contrast to the above-mentioned works, this paper takes a different approach. Instead of focusing on a high fidelity model, we present a simplified model that can be readily used to obtain the differential micro-motion kinematics Jacobian and is readily amenable to formulating a model calibration problem. This micro motion Jacobian is essential for control purpose, and an associated identification Jacobian is needed for calibration purposes. Therefore, the paper focuses on the derivation of the micro-motion kinematics and its associated identification Jacobian for calibration and error prorogation.

This work is built upon our previous work [11], in which we presented the concept of CREM and provided a visual measurement solution to observe micro-motion. The work in [11] lacked a modeling approach that can explain the experimental observations and that can be used for control and identification purposes. The contribution of this work is in presenting a simplified kinematic modeling framework that captures the micro-motion achieved by the equilibrium modulation of continuum robots, and in developing a calibration approach to capture the model parameters. We put forth the concept of moment coupling effect as a simplified approach to describe the equilibrium modulation behavior, and thereby, both direct kinematics and instantaneous kinematics are formulated for control purposes. To account for errors potentially caused by the simplistic modeling assumptions, a modeling uncertainty term is introduced, and the identification Jacobian along with a calibration framework to capture the parameterization is developed. Using the multi-backbone continuum robot design in [11] as a validation platform, we validate the kinematic model and model calibration experimentally while augmenting these results with additional simulation validations.

2 Equilibrium Modulation Backbone Insertion to Create Micro Scale Motion

This section presents the bending shape equilibrium modeling in the case where the Equilibrium Modulation Backbone (EMB) insertion is at a given depth $q_{s}$. To motivate the modeling approach taken we will first refer the readers to [13] where the simplified kinematics of multi-backbone continuum robots was presented. When the EMBs are not inserted and for proper design parameters (e.g. small spacing between the spacer disks) the continuum segment bends in a constant curvature [32]. We use this assumption to create a simplistic equilibrium model which lends itself to fast real-time computation. Since we have to account for modeling uncertainties due to friction and material parameter uncertainties, we later lump the error of the simplified model in an uncertainty term $\lambda$ that will be used to produce an updated CREM model.

2.1 Simplistic Equilibrium Model

Figure 2 shows the free body diagram of a continuum segment with and without an inserted EMB. In Fig. 2(b), a separation plane is defined at the insertion depth $q_{s}$, dividing the segment into two subsegment - Inserted and Empty. Though not accurate, the two subsegments are both assumed to have constant but different curvatures. The angles $\theta^{\prime}$ and $\theta_{s}$ denote the bending angles of the end-disk and at the insertion depth, respectively, when the EMB is inserted. The angle $\theta$ denotes the nominal bending angle when the EMB is not inserted. The angle $\theta_{0}=\pi/2$ denotes the angle at the base of the segment.

We first consider the resultant moment $m_{1}$ that the backbones apply on any imaginary cross section of the continuum segment when no EMB is inserted (Fig. 2(a)):

[TABLE]

Where $E_{p}$, $E_{i}$ and $I_{p}$, $I_{i}$ denote the Young’s moduli and cross-sectional bending moments of inertia of the central backbone and the $i^{\textit{th}}$ secondary backbone, respectively. Also $L$ and $L_{i}$ denote the lengths of the central backbone and the $i^{\textit{th}}$ secondary backbone.

We also consider the moment ${m_{1}}^{\prime}$ along the empty subsegment in the case of EMB being inserted (Fig. 2(b)):

[TABLE]

Where $L_{\varepsilon_{i}}$ denotes the $i^{\textit{th}}$ backbone length portion that belongs to the empty subsegment (this is the arc-length from the separation plane until the end-disk).

We next use key definitions from [33]. The radial distance between the secondary backbones and the primary backbone is denoted $r$. When $r$ is projected onto the plane in which a segment bends, we obtain the projected radial distance $\Delta_{i}$:

[TABLE]

where $\sigma_{i}$ designates the angular coordinate of the $i^{\,th}$ backbone relative to the bending plane. The angular coordinate of the first backbone relative to the bending plane is $\delta$ and the angular separation between secondary backbones is $\beta=\frac{2\pi}{n}$ where $n$ is the number of secondary backbones.

The length of the $i^{\,th}$ backbone, $L_{i}$ is derived using the fixed radial offset between the backbones:

[TABLE]

Using similar rationale, we calculate the empty length portion $L_{\varepsilon_{i}}$ and the inserted length portion of the $i^{\,th}$ secondary backbone $L_{s_{i}}$:

[TABLE]

In both Fig. 2(a) and (b), the static equilibrium at the end-disk is determined by the geometric constraints and the backbone loading forces at the end-disk. For example, coordinated pulling and pushing on all secondary backbones are assumed to form a force couple that generates only a moment at the end-disk.

We next use a simplifying assumption that the effect of EMB wire insertion on changes in the bending curvatures of the un-inserted subsegment backbones is negligible, hence:

[TABLE]

Next, we consider $m_{2}$ and $m_{s}$, the moments that the secondary backbones and the EMB apply on the separation plane as shown in Fig. 2(b):

[TABLE]

Where $E_{s}$ and $I_{s}$ denote the Young’s modulus and cross-sectional bending moment of inertia of the EMB.

Substituting equations (1, 2) into (7), results in one equation having two unknowns $\theta^{\prime}$ and $\theta_{s}$ as illustrated in Fig. 2(b). To obtain the second equation necessary for solving for these two unknowns, we use the moment balance on the separation plane:

[TABLE]

To solve equations (10) and (7) for the unknowns $\theta^{\prime}$ and $\theta_{s}$ we explicitly express the backbone moments using the beam equation $m=EI\kappa$ where $\kappa$ designates the radius of curvature and $EI$ designates the bending cross sectional stiffness of a beam. In doing so, we note that the curvature of a beam bent in a circular shape satisfies $\kappa=\frac{\theta}{L}$ where $\theta$ is the deflection angle and $L$ is the beam length. Since the backbone lengths are a function of the unknowns, we will rewrite the moment equation for a beam as $m=\frac{EI}{L}\theta$ and, by defining the beam angular deflection stiffness $k_{\theta}\triangleq\frac{EI}{L}$ we obtain a simple equation for the moment $m=k_{\theta}\theta$.

Using the above definition for beam angular deflection stiffness, we rewrite the moment equations for each beam as:

[TABLE]

Substituting Eq. (12)-(14) in Eq. (10) results in:

[TABLE]

Substituting equations (11) and (12) in Eq. (7) results in:

[TABLE]

As a final step in the solution, we substitute the result in Eq. (16) in Eq. (15), thereby obtaining $\theta_{s}$ and subsequently $\theta^{\prime}$.

2.2 Updated CREM Model

Equations (15) and (16) present the solution to the simplistic modeling approach that is fundamentally based on Eq. (7) and Eq. (10). In addition to the simplified assumption, the current model also neglects modeling uncertainties due to frictional effects and material property uncertainties. Prior works in [34, 35] show that these uncertainties include friction and strain along the actuation lines, non-uniformly distributed load on backbones that causes shape deviations from constant-curvature bending, deviations in the cross section of the backbones during bending, and uncertainties in the properties of the NiTi backbones.

To account for the modeling uncertainties caused by friction, material uncertainty111Manufacturer-specified Young’s modulus for superelastic NiTi is typically provided with a wide range of 40-70 GPa, and our simplistic model, we introduce an uncertainty term $\lambda$ to Eq. (10):

[TABLE]

The uncertainty term $\lambda$ captures effects of EMB insertion offset, bending angle uncertainty and a fixed offset:

[TABLE]

The solution in Eq. (15) is also updated to:

[TABLE]

Having obtained the solutions to the equilibrium tip bending angle $\theta^{\prime}$ and the bending angle at the separation plane $\theta_{s}$, we introduce an equilibrium configuration space variable vector $\bm{\phi}$ to combine them. With the purpose of preparing for kinematic derivations in later sections when we break a single continuum segment down to two subsegments, the vector $\bm{\phi}$ is defined as:

[TABLE]

Where $\theta_{\varepsilon}$ represents the bending angle of the empty subsegment.

We define the configuration space variable $\bm{\psi}$ as the nominal bending angle $\theta$ and the bending plane angle $\delta$:

[TABLE]

Finally, the solution to equilibrium modeling problem is presented as a mapping $\boldsymbol{\mathscr{F}}_{\text{eqm}}$ which is used in deriving the Jacobian matrices in the following sections:

[TABLE]

Equation (22) provides the end disk equilibrium angle for a combination of any given EMB insertion length $q_{s}$, nominal bending angle $\theta$, and bending plane angle $\delta$.

3 Kinematic Modeling

Kinematic modeling of CREM includes the mapping of configuration space to task space and its differential kinematics. The differential kinematics include the instantaneous kinematics and the error propagation.

The instantaneous kinematics is derived for control purpose, and it includes two motion Jacobian matrices that both relate actuation speeds to the robot tip velocity. The macro motion Jacobian $\mathbf{J}_{M}$ is associated with the joint velocities of push/pull on backbones (direct actuation) while the micro motion Jacobian $\mathbf{J}_{\mu}$ is related to the EMB insertion velocity (indirect actuation).

The kinematic error propagation investigates how errors in parameters contribute to errors in kinematic measurements of task space (e.g. measured positions). In this work, we focus on the vector $\mathbf{k}_{\lambda}$ that parameterizes the modeling uncertainty. Other robot geometric kinematic parameters can be calibrated following [36]. An identification Jacobian $\mathbf{J}_{\mathbf{k}}$ is derived and used in section 4 to estimate $\mathbf{k}_{\lambda}$ with experimental measurements.

3.1 Kinematic Modeling Using Mapping $\boldsymbol{\mathscr{F}}_{\text{eqm}}$

With the mapping $\boldsymbol{\mathscr{F}}_{\text{eqm}}$ in Eq. (22) derived as the result of static equilibrium, the kinematic mapping can be formulated by considering a single continuum segment as two concatenated subsegments - the inserted and the empty, divided at the insertion depth $q_{s}$. Figure 3 illustrates our approach to analyzing the two concatenated subsegments. The bending angles of both subsegments were introduced in Eq. (20), denoted as $\theta_{s}$ and $\theta_{\varepsilon}$, for the inserted and the empty subsegment. Since the whole segment is assumed to bend in plane, both subsegments have the equal bending plane angles:

[TABLE]

The kinematic nomenclature used in this paper refers to Table 1 (shown in Fig. 3).

Recalling the direct kinematics of a single segment [34] having length $L_{x}$ and an end disk angle $\theta_{x}$, the end disk pose (position and orientation) relative to the base are given by:

[TABLE]

Where $\delta_{x}$ designates the angle of the bending plane (analogous to $\delta$ in Fig. 3), $\left[\mathbf{v}^{\wedge}\right]$ represents the cross-product matrix of vector $\mathbf{v}$ and the matrix exponential $e^{\;\alpha\left[\mathbf{v}^{\wedge}\right]}$ represents a rotation matrix about the axis $\mathbf{v}$ by an angle $\alpha$.

To obtain the pose of the end disk of the inserted segment is given by ${}^{b}\mathbf{p}_{c/b}$ and ${}^{b}\mathbf{R}_{c}$ we substitute $L_{x}=q_{s},\theta_{x}=\theta_{s},\delta_{x}=\delta$ in Eqs. (24, 25). Similarly, the pose of the end disk of the empty segment relative to its base is obtained by substituting $L_{x}=L-q_{s},\theta_{x}=\theta_{\varepsilon},\delta_{x}=\delta$ in Eqs. (24, 25) to result in ${}^{c}\mathbf{p}_{g/c}$ and ${}^{c}\mathbf{R}_{g}$.

The pose of the free subsegment end disk relative to the segment base is given by:

[TABLE]

Casting the above two equations in a homogeneous transform format yields:

[TABLE]

With $\bm{\phi}$ expressed using mapping $\boldsymbol{\mathscr{F}}_{\text{eqm}}(\bm{\psi},q_{s},\mathbf{k}_{\lambda})$, the forward kinematics is determined, which can be also written as:

[TABLE]

3.2 Differential Kinematics

The total differential of a homogenous transformation $\mathbf{T}\in\text{SE}(3)$, may be represented as:

[TABLE]

Where $\textrm{d}\mathbf{x}$ and $\textrm{d}\bm{\mu}$ represent translational and rotational differentials in the base frame222A differential rotation is a sequence of rotations of small angles.. The vector $\bm{\mu}\in\hbox{I \kern-3.99994ptR}^{3\times 1}$ represents a chosen orientation parametrization (e.g. Euler angles).

The total differential of ${}^{b}\mathbf{T}_{g}$ is obtained by considering differentials on all variables, i.e. $\textrm{d}\bm{\phi}$, $\textrm{d}\delta$, and $\textrm{d}q_{s}$:

[TABLE]

Using the nomenclature of a Jacobian $\mathbf{J}_{ab}$ such that $\delta a=\mathbf{J}_{ab}\delta b$, we define the following Jacobian matrices:

[TABLE]

The Jacobian matrices $\mathbf{J}_{\bm{\xi}\bm{\phi}}$, $\mathbf{J}_{\bm{\xi}\delta}$, and $\mathbf{J}_{\bm{\xi}q_{s}}$, respectively, relate the differential on equilibrium configuration space variable $\textrm{d}\bm{\phi}=[\textrm{d}\theta_{s},\textrm{d}\theta_{\varepsilon}]^{\mathrm{T}}$, the differential on bending plane angle $\textrm{d}\delta$, and differential on EMB insertion depth $\textrm{d}q_{s}$, to the corresponding differential contributions on the pose vector $\textrm{d}\bm{\xi}$.

Both $\mathbf{J}_{\bm{\xi}\bm{\phi}}$ and $\mathbf{J}_{\bm{\xi}\delta}$ can be obtained by treating the inserted and empty subsegments as a concatenated two-segment robot, which is explained in section 3.3.

The third Jacobian, $\mathbf{J}_{\bm{\xi}q_{s}}$, defined as the partial derivative, ${\partial\bm{\xi}}/{\partial q_{s}}$, is derived with the other variables ($\bm{\phi}$ and $\delta$) held constant. The end-effector orientation, given by ${}^{b}\mathbf{R}_{g}$ in Eq. (27) is not a function of $q_{s}$. Therefore, by considering only the translational differential due to $\textrm{d}q_{s}$, we have:

[TABLE]

Where $\tfrac{\partial\;{}^{b}\mathbf{p}_{c/b}}{\partial q_{s}}$ and $\tfrac{\partial\;{}^{c}\mathbf{p}_{g/c}}{\partial q_{s}}$ are derived from Eq. (24). It is important to note that $\mathbf{J}_{\bm{\xi}q_{s}}$ differs from the micro motion Jacobian $\mathbf{J}_{\mu}$ derived later in that $\mathbf{J}_{\bm{\xi}q_{s}}$ is a contributing part of $\mathbf{J}_{\mu}$ - the length ‘tangential’ contribution, while $\textrm{d}q_{s}$ also propagates to $\textrm{d}\bm{\phi}$ that also causes change on $\textrm{d}\bm{\xi}$.

Having derived the above three Jacobian matrices, $\mathbf{J}_{\bm{\xi}\bm{\phi}}$, $\mathbf{J}_{\bm{\xi}\delta}$, and $\mathbf{J}_{\bm{\xi}q_{s}}$, we obtain the pose total differential $\textrm{d}\bm{\xi}$ expressed using differentials $\textrm{d}\bm{\phi}$, $\textrm{d}\delta$, and $\textrm{d}q_{s}$. Further, the differential $\textrm{d}\bm{\phi}$ is a result of multiple other differentials, which can be seen from mapping $\boldsymbol{\mathscr{F}}_{\text{eqm}}(\bm{\psi},q_{s},\mathbf{k}_{\lambda})$. To fully investigate and decouple the contributions of direct (macro) and indirect (micro) actuation, we express $\textrm{d}\bm{\phi}$ using differentials on $(\bm{\psi},q_{s},\mathbf{k}_{\lambda})$. This differentiation is also motivated by Eq. (29), where the variables are decoupled as $\bm{\psi}$ for macro motion, $q_{s}$ for micro motion, and $\mathbf{k}_{\lambda}$ for micro motion parameters. Such differentiation is derived as:

[TABLE]

Where the gradient terms are derived in Appendix A as:

[TABLE]

Matrices $\mathbf{A}^{\prime}_{x},\mathbf{B}^{\prime}_{x}$ are partial derivative matrices with respect to variable ‘$x$’, and $\mathbf{A},\mathbf{S}_{0},\mathbf{B},\mathbf{C}_{0},\mathbf{S}_{1}$ are defined as:

[TABLE]

Using Eq. (37-52), $\textrm{d}\bm{\phi}$ is fully expressed as Eq. (35). Substituting $\textrm{d}\bm{\phi}$ into the original differentiation in Eq. (32), we obtain the full differential kinematics that relates differentials on $\{\bm{\psi},q_{s},\mathbf{k}_{\lambda}\}$ to the pose total differential $\textrm{d}\bm{\xi}$:

[TABLE]

Rewriting Eq. (54) using the Jacobian definitions:

[TABLE]

Collecting like terms of $\textrm{d}\bm{\psi}$, $\textrm{d}q_{s}$, and $\textrm{d}\mathbf{k}_{\lambda}$, we obtain a differentiation that decouples differentials of the macro motion, the micro motion, and the parameters:

[TABLE]

The above result completes the mapping from configuration to task space. It clearly delineates the effects of EDM insertion and direct actuation to achieving macro and micro motion. For control purposes, a complete mapping from joint to task space is needed. We therefore consider next the mapping from direct (macro) actuation joint space $\mathbf{q}$ to task space $\bm{\xi}$. Since three secondary backbones are used in our experimental setup as illustrated in Figure 3, we will define $\mathbf{q}\triangleq[q_{1},q_{2},q_{3}]^{\mathrm{T}}$ where:

[TABLE]

When obtaining this mapping, we consider the nominal segment kinematics for multi-backbone continuum robots as in [13].

The Jacobian that relates the differential $\textrm{d}\mathbf{q}$ to the differential $\textrm{d}\bm{\psi}$ was reported in [32] as:

[TABLE]

Where $r$ denotes the constant radial distance between the central and surrounding backbones, and $\beta=2\pi/3$ denotes the backbone separation angle. Using Eq. (57), we substitute $\textrm{d}\bm{\psi}$ into Eq. (55), arriving at the final differential kinematics:

[TABLE]

Equation (58) fully decouples the end-effector pose differential to contributions of the direct (macro) actuation $\textrm{d}\mathbf{q}$, the indirect (micro) actuation $\textrm{d}q_{s}$, and the modeling uncertainty $\textrm{d}\mathbf{k}_{\lambda}$. The three Jacobian matrices are obtained from Eq. (58) as an important finding of this paper: $\mathbf{J}_{M}$ defined as the Macro motion Jacobian, $\mathbf{J}_{\mu}$ defined as the Micro motion Jacobian, and $\mathbf{J}_{\mathbf{k}}$ defined as the Identification Jacobian.

[TABLE]

Where $(\cdot)^{\dagger}$ is the Moore-Penrose pseudo inverse.

3.3 Deriving $\mathbf{J}_{{\xi}{\phi}}$ and $\mathbf{J}_{{\xi}\delta}$

The result in Eq. (58) builds on knowing the Jacobians $\mathbf{J}_{\bm{\xi}\bm{\phi}}\in\hbox{I \kern-3.99994ptR}^{6\times 2}$ and $\mathbf{J}_{\bm{\xi}\delta}\in\hbox{I \kern-3.99994ptR}^{6\times 1}$ as mentioned in Eqs. (32, 33). We now provide a derivation to these two Jacobians. Considering a single-segment CREM as two concatenated subsegments (inserted and empty), we apply the Jacobian formulation for a two-segment multi-backbone continuum robot (MBCR) while assuming that both subsegments share the bending plane angle $\delta$. For the ease of adapting formulations from [32], we introduce a vector notation:

[TABLE]

We next use the notation of ${}^{i-1}\bm{\xi}_{i/i-1}$ to denote the pose of the $i^{\text{th}}$ subsegment relative to the $(i-1)^{\text{th}}$ subsegment where $i\in\{s,\varepsilon\}$. Using $\mathbf{v}$ and $\bm{\omega}$ to denote linear and angular velocities, we define the corresponding four Jacobian matrices corresponding with the contributions of $\textrm{d}\theta_{i},\textrm{d}\delta_{i}$ where $i\in\{s,\varepsilon\}$ to the end-effector twist:

[TABLE]

Details of the derivations of the Jacobians are provided in Appendix B.

Following [37], the serial composition of two subsegments using twist transformation results in the end effector twist:

[TABLE]

These definitions of $\mathbf{J}_{\bm{\xi}\bm{\phi}}$ and $\mathbf{J}_{\bm{\xi}\bm{\delta}_{v}}$ complete the two missing terms needed in Eq. (58), but with a slight formulation modification. The Jacobian matrix $\mathbf{J}_{\bm{\xi}\delta}$ is slightly different from $\mathbf{J}_{\bm{\xi}\bm{\delta}_{v}}$ in Eq. (70), and using the differentiation chain rule it becomes:

[TABLE]

4 Calibration of Micro Motion Parameters

To calibrate the model uncertainty parameters $\mathbf{k}_{\lambda}$, we extract from Eq. (58) the following relation:

[TABLE]

Using this error propagation model, we construct an estimation method to estimate $\mathbf{k}_{\lambda}$. Let $\bm{\xi}_{j}\leftrightarrow\left[\bar{\mathbf{x}}_{j},\bar{\mathbf{R}}_{j}\right]$ designate the measured end-effector pose at the $j^{\text{th}}$ robot configuration (insertion depth) where $\bar{\mathbf{x}}_{j}$ and $\bar{\mathbf{R}}_{j}$ designate the measured position and orientation. Let $\mathbf{x}_{j}$ and $\mathbf{R}_{j}$ denote the modeled pose using the direct kinematics as presented in section 3.1 for a given $\mathbf{k}_{\lambda}$. The error between the measured and modeled poses are then defined as:

[TABLE]

where $\alpha_{e_{j}}$ and $\hat{\mathbf{m}}_{e_{j}}$ are the angle and axis parameterizing the orientation error $\mathbf{R}_{e_{j}}$. These parameters are given by:

[TABLE]

where the operator $(\cdot)^{\vee}$ designates the vector form of a skew-symmetric matrix.

An aggregated error vector $\widetilde{\mathbf{c}}_{\lambda}$ is defined to include errors of all $N$ robot configurations:

[TABLE]

The optimization objective function $M_{\lambda}$ is then defined as:

[TABLE]

Where $\mathbf{W}$ is a weight matrix encoding confidence in the measurements and the measurement unit scaling factors.

The first-order Taylor series approximation of $M_{\lambda}$ is given:

[TABLE]

where the aggregated Jacobian $\mathbf{J}_{M_{\lambda}}\in\hbox{I \kern-3.99994ptR}^{1\times 5}$ is given by:

[TABLE]

Equation (80) shows that minimizing $M_{\lambda}$ entails following the gradient descent direction along $({\partial\widetilde{\mathbf{c}}_{\lambda}}/{\partial\mathbf{k}_{\lambda}})$. The parameter $\mathbf{k}_{\lambda}$ is then obtained using an iterative nonlinear least squares estimation shown in Algorithm 1.

In the algorithm, $\mathbf{H}$ is the parameter scaling matrix and the task space variable scaling is achieved by adjusting $\mathbf{W}$, both of witch are discussed in details in [38].

5 Simulation Study of Direct Kinematics and Differential Kinematics

In this section, we present simulations to demonstrate the direct kinematics and differential kinematics. We also verify the differential kinematics through finite-difference simulations. We also carry out simulations to verify the differential kinematics model. In all simulations, we assumed the robot points vertically down at its home (straight) configuration.

5.1 Position Analysis of Micro Motion

Using the model in section 3.1, we present the simulated position analysis of the micro motion created by the EMB insertion. In both simulations and the experimental model validations, we use the parameters as in Table 2. They include the Young’s modulus of the superelastic NiTi material used for backbone tubes and EMBs ($E_{p}$, $E_{i}$, $E_{s}$), the diameters of backbones ($d_{p}$, $d_{i}$, $d_{s}$), and the cross-sectional moment of inertia ($I_{p}$, $I_{i}$, $I_{s}$).

Figure 4 shows the simulation results of the micro motion created by EMB insertion. Figure 4(a) shows the continuum segment at its initial bending angle $\theta=30^{\circ}$ . During simulation, the equilibrium bending angles $\{\theta^{\prime},\theta_{s}\}$ were computed at different EMB insertion depths. The resulting tip micro-motion is shown in Fig. 4(b) for the naïve kinematic model (i.e. $\lambda=0$). Figure  4(c) shows the tip motion for an updated model assuming $\lambda=0.2+0.025q_{s}$. We note that, as expected, in both cases the robot straightens with EMB insertion since the robot body straightens. However, the updated model exhibits a turning point behavior which relates to the combined effect of straightening and change in the end effector angle $\theta^{\prime}$. This same phenomenon was observed experimentally in section 6. The particular values used for the parameters defined in Eq. (18), $k_{\lambda_{0}}=0.2,k_{\lambda_{q}}=0.025$, are manually selected to illustrate the turning point behavior that is similar to what is observed in experiments. In practise, they are calibrated according to section 6.2.

5.2 Instantaneous Kinematics and Error Propagation

To verify the derivations of instantaneous kinematics and error propagation, we compute Jacobians following section 3.2. Since the simulation case represents the robot motion within a bending plane, the columns of the Jacobians represent $2\times 1$ vectors of induced velocities for unit change in the variables associated with each Jacobian. The following simulations verify the macro motion Jacobian $\mathbf{J}_{M}$, the micro motion Jacobian $\mathbf{J}_{\mu}$, and the identification Jacobian $\mathbf{J}_{\mathbf{k}}$ by plotting the Jacobian columns. The verification is carried out visually by verifying that the Jacobian columns induce tip velocity tangent to the trajectory generated by direct kinematics. In addition, each Jacobian has been also verified numerically via finite difference computations.

To verify $\mathbf{J}_{M}$, the EMB insertion depth $q_{s}$ was fixed and direct actuation of backbones was assumed. Sample tip positions along the trajectory were obtained via direct kinematics and the corresponding Jacobian $\mathbf{J}_{M}$ was computed. Figure 5 shows the simulation results. These results verify that the computed $\mathbf{J}_{M}$ is tangent to the macro motion trajectory.

To verify $\mathbf{J}_{\mu}$, the secondary backbones were assumed locked and the EMB insertion depth $q_{s}$ was varied. The Jacobian $\mathbf{J}_{\mu}$ was computed and plotted for each EMB depth. Two different scenarios of modeling uncertainty were considered: $\lambda=0$ and $\lambda=0.2+0.025q_{s}$. The results in Fig. 6(a)(b) verify that $\mathbf{J}_{\mu}$ is tangent to the micro scale trajectory generated by direct kinematics.

Figure 6(c) shows the plots of the identification Jacobian $\mathbf{J}_{\mathbf{k}}$ for the simulation scenario where $\lambda\neq 0$, revealing how the parameter errors of modeling uncertainty affect the tip positions and hence the shape of the trajectory.

6 Experimental Validations

In [11] the feasibility of micro motion through equilibrium modulation was demonstrated. The following experiments evaluate the ability of our simplified kinematic model to capture the micro-motion behavior, validate the calibration framework in section 4, and assess the accuracy of the updated kinematic model in reflecting the experimental data.

6.1 Experimental Setup & Ground Truth Data

A single-segment continuum robot with EMB insertion actuation was used as the experimental platform, Figure 7. The platform was presented in [11], and it was modified from an earlier multi-backbone continuum robot design [32]. To observe the robot tip motion at different scales, one HD camera (FLIR Dragonfly II®) was used to capture the macro motion and the bending shape while an identical camera mounted on a $22.5\times$ microscope lens to capture the micro motion. Custom “multi-circled” marker was used to track the tip motion under microscope while multiple custom “X” markers were attached to the continuum robot body to observe the bending shape. The vision measurement methods used were presented in [11] with the micro motion tracking accuracy being reported better than 2 $\mu$m.

Fig. 7 shows the frames used and also previously referred to in Fig. 4. The microscope is fixed at a known offset relative to the robot base, and such offset is represented as a constant transformation from the image frame {I} to the robot base frame {B}. The tracked marker frame {M} is placed at a known offset relative to the end disk (gripper frame {G}), and the transformation is represented as a constant transformation between {M} and {G}. The marker position and orientation in the image frame is obtained by the segmentation of the three circles that construct an asymmetric pattern, as shown in Fig. 7(c).

The multimedia extension and Fig. 8 show a sample marker frame trajectory during EMB insertion where the macro motion bending angle was chosen to be $\theta=45^{\circ}$. The marker positions were segmented from microscope images collected at 30 frames per second. Applying a butterworth infinite impulse response filter with the 3-dB frequency as 30 Hz, provided a smooth trajectory for calibration. Results of more experiments are plotted in Fig. 9: five macro bending angle configurations, $\theta=15^{\circ},30^{\circ},45^{\circ},60^{\circ},75^{\circ}$, were chosen to sample the workspace, and for each value of $\theta$, ten repetitions of the EMB insertion experiment were conducted. Figure 9 shows that the turning point phenomenon exists over the entire workspace.

6.2 Model Calibration

Using calibration method in section 4, we calibrated the modeling uncertainty parameter vector $\mathbf{k}_{\lambda}$ . The parameter vector $\mathbf{k}_{\lambda}$ in Eq. (18) consists of three elements, a bias term $k_{\lambda_{0}}$, a coefficient gain $k_{\lambda_{\theta}}$ that is associated with the nominal bending angle $\theta$, and a coefficient gain $k_{\lambda_{q}}$ that relates to the EMB insertion depth $q_{s}$. As a preliminary study, in this paper, we focus on investigating and calibrating $k_{\lambda_{0}}$ and $k_{\lambda_{q}}$. To obtain the value of $k_{\lambda_{\theta}}$, there are multiple possible solutions. In one method, after collecting observation data of micro motion trajectories from sufficient groups varying bending configuration $\theta$, one could use the same calibration method with the identification Jacobian $\mathbf{J}_{\mathbf{k}}$ that includes the third column corresponding to $k_{\lambda_{\theta}}$. In another method, one could first calibrate $k_{\lambda_{0}}$ and $k_{\lambda_{q}}$ for each micro motion trajectory that is observed in a particular bending configuration $\theta$, followed by generating an interpolating lookup table in the format of $(k_{\lambda_{0}},k_{\lambda_{q}})=f_{\text{Lookup}}(\theta)$.

Algorithm 1 was initialized with $k_{\lambda_{0}}=0$, $k_{\lambda_{q}}=0$. In each iteration, the modeled positions were computed using the current estimates of the parameters. The aggregated error vector was then calculated between the modeled and experimental positions. For each iteration, both of the current estimates of the parameters and the position root-mean squared errors (RMSE) of all insertion samples (382 in total) were reported. A relative convergence threshold of $0.1\%$ was used to determine the convergence.

For the particular experimental data collection shown in Fig. 8, the parameter estimation (model calibration) went through 46 iterations before converging, where a step size of $\eta=0.1$ was selected ($\eta$ as introduced in Algorithm 1). Figure 10(a) shows selected iterations during the estimation, and the details of the iterations are reported in Table 3. The estimation started with an initial position RMSE of 44.27 $\mu$m, and after its convergence, the position RMSE was reduced to 5.82 $\mu$m, showing an improvement of 86.8% in model errors.

By dividing the tip trajectory into two segments, we observe that the current simplistic modeling approach produced bigger errors after the turning point: the RMSEs were reported as $4.87\mu m$ and $6.63\mu m$ for the two segments before and after turning point that had the lengths of $48.11\mu m$ and $38.82\mu m$, respectively. If one wishes to further improve the model accuracy, a model that only predicts the trajectory before the turning point may be considered. We therefore considered another calibration where only the partial micro motion trajectory before the turning point was used. With the same iteration step size and convergence criterion, the estimation went through 59 iterations to converge, and the updated results were reported in Table 4 and plotted in Fig. 10(b). The position RMSE was then further improved to $4.76\mu m$.

6.3 Limitations

This work focused on creating a simplistic, yet fast model for equilibrium modulation control implementation. The kinematic model traded accurate mechanics modeling (which typically leads to solving nonlinear boundary value problems) with simplicity and speed of computation. In this simplistic model, one EMB is assumed to be inserted through the central backbone, or equivalently, multiple EMBs are assumed to be inserted in coordination such that a shared insertion separation plane can be defined (shown as Fig. 2). Hence the current method limits the use of modeled micro motion to be coplanar with the macro bending plane, making the controllable micro motion have only one degree of freedom (DoF). A more sophisticated modeling method could potentially capture any arbitrary insertions of multiple EMBs, which enables the full capability of the micro motion mechanism that can generate spatial motions, for example, spiral motion.

Our experimental data showed an unexpected motion behavior manifested by a turning point along the micro-motion trajectory. The model presented in this paper does not offer a physical explanation to this behavior, but can capture this behavior for a given robot. The model calibration results indicate that there is still a potential to improve the model performance by further investigating alternative modeling assumptions and different descriptions of modeling uncertainties. One of the limitations of our approach can be inferred from the simulation shown in Fig. 6, where both columns of the identification Jacobian are almost aligned with the tangent to the direct kinematics trajectory. The attainable correction directions that the column-space of $\mathbf{J}_{\mathbf{k}}$ affords is therefore limited in reshaping the model trajectory. This was also observed from the experimental validation. Shown from the iterations in Fig. 10, it is difficult to reshape the modeled tip trajectory in the direction that is perpendicular to the trajectory. The other limitation is potentially caused by the choice of linearity in expressing the modeling uncertainty, which may not be descriptive enough.

7 Conclusion

This work presented the first modeling attempt for a new class of continuum robots capable of multi-scale motion. These robots achieve macro-scale and micro-scale motions through direct and indirect actuation (equilibrium modulation). Instead of focusing on a high-fidelity mechanics-based model, which typically leads to non-linear boundary value problems not easy to adopt for real-time control or parameter identification, this paper presented a simplified mechanics-based model utilizing moment coupling effects between sub-segments of the continuum robot. This approach generates a differential kinematics model that covers both macro and micro-motion. As a result of unavoidable parameter uncertainty, we presented a model-calibration approach that can compensate for parameter inaccuracy, friction effects and modeling inaccuracies due to the simplistic modeling assumptions. The modeling approach along with the calibration framework was validated experimentally on a multi-backbone continuum robot. The calibrated model reported a positional root-mean-squared error as 5.83 $\mu$m if one wishes to use the model for the entire motion profile with the turning point. If one chooses to exclude motions past the turning point, the calibrated model fit the experimental data with an accuracy of 4.76 $\mu m$. Future work will focus on investigations on a more sophisticated models capable of incorporating geometric constraints as well as minimizing mechanical energy for improved model accuracy. In addition, effects of direction reversal of EMB insertion can manifest in hysteresis, which has not been explored in this work, but still remains the topic of ongoing research.

8 Multimedia Extension

The video shows the micro-motion trajectory during EMB insertion. The following is a link that would ideally be placed and linked on AMSE website if the paper is accepted. If ASME has no way of hosting these videos, then we can maintain them on our website.

Multimedia https://youtu.be/4GlLQwXUJpI. {acknowledgment} This work was supported by NSF grant #CMMI-1537659. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

Appendix A Deriving $\tfrac{\partial{\phi}}{\partial\theta}$, $\tfrac{\partial{\phi}}{\partial\delta}$, $\tfrac{\partial{\phi}}{\partial q_{s}}$ & $\tfrac{\partial{\phi}}{\partial\mathbf{k}_{\lambda}}$

Rewriting Eq. (19) and Eq. (16) in a matrix form yields:

[TABLE]

where $\mathbf{A}$ and $\mathbf{B}$ are defined as functions of $\{\bm{\psi},q_{s},\theta_{s}\}$ and $\{\bm{\psi},q_{s},\mathbf{k}_{\lambda},\theta_{s}\}$ respectively. Using the definition of $\bm{\phi}$, yields:

[TABLE]

By introducing two constant matrices in the above equation, $\mathbf{S}_{0}$ and $\mathbf{C}_{0}$, we have obtained the equation to differentiate:

[TABLE]

The full differentiation may be expressed as:

[TABLE]

Using $\mathbf{X}^{\prime}_{a}$ to denote the partial derivative of matrix $\mathbf{X}$ w.r.t the scalar variable $a$, i.e. $\mathbf{X}^{\prime}_{a}\triangleq\tfrac{\partial\mathbf{X}}{\partial a}$, then $\textrm{d}\mathbf{A}$ and $\textrm{d}\mathbf{B}$ may be written as:

[TABLE]

Let us define $\mathbf{C}_{\bm{\phi}}$ and $\bm{\Gamma}_{a}$ to provide ease in the derivations:

[TABLE]

where the letter $a\,\in\{\theta,\delta,q_{s},\theta_{s}\}$.

By substituting Eq. (A.5) and Eq. (A.6) into Eq. (A.4), and by using the definitions of $\mathbf{C}_{\bm{\phi}}$ and $\bm{\Gamma}_{a}$, we have:

[TABLE]

This equation shows the full differentiation of Eq. (35) and all the Jacobians can be obtained directly by their definitions, i.e., the expressions of $\left\{\tfrac{\partial\bm{\phi}}{\partial\theta},\tfrac{\partial\bm{\phi}}{\partial\delta},\tfrac{\partial\bm{\phi}}{\partial q_{s}},\tfrac{\partial\bm{\phi}}{\partial k_{\lambda_{i}}}\in\hbox{I \kern-3.99994ptR}^{2\times 1}\right\}$ may be written as:

[TABLE]

where $\mathbf{S}_{1}=[1,0]$ is just a selection matrix.

Appendix B Derivation of The Jacobian Partitions for the Multi-segment case

Equations (64), (67) and (70) refer to the Jacobian matrix partitions for the two-segment case where the first segment is the inserted subsegment (indicated vby superscript $s$) and the second segment is the empty subsegment (indicated by subscript $\varepsilon$). The Jacobian partitions represent the effects of differentials on $\theta_{i}$ and $\delta_{i}$ that contribute to the end-effector’s translational and rotational differential, labeled by subscripts ‘$\mathbf{v}$’ and ‘$\bm{\omega}$’, indicating ‘velocity’ and ‘angular velocity’, respectively. The expressions of {$\mathbf{J}_{\mathbf{v}\theta_{i}}$, $\mathbf{J}_{\bm{\omega}\theta_{i}}$, $\mathbf{J}_{\mathbf{v}\delta_{i}}$, $\mathbf{J}_{\bm{\omega}\delta_{i}}$} are extracted from [37] as:

[TABLE]

Where $c(\cdot)$ and $s(\cdot)$ denote the cosine and sine functions, and $D_{i}$ represents the length of the subsegment. For the inserted subsegment, $D_{s}=q_{s}$; and for the empty subsegment, $D_{r}=L-q_{s}$. In addition, the following shorthanded notations are used:

[TABLE]