CONTINUUM DEVICES AND CONTROL METHODS THEREOF

[0001] This application claims the benefit of Provisional Application Ser. No. 61/879,727, entitled “FLEXIBLE INSTRUMENT WITH PRE-CURVED ELEMENTS” and filed Sep. 19, 2013. This application is also a continuation-in-part of U.S. Non-Provisional application Ser. No. 13/808,994, entitled “CONTINUUM ROBOTS AND CONTROL THEREOF” and filed Feb. 5, 2013, which is a National Stage of International Application Serial No. PCT/US2011/38539, entitled “CONTINUUM ROBOTS AND CONTROL THEREOF” and filed May 31, 2011, which claims the benefit of Provisional Application Ser. No. 61/362,353, entitled “CONTINUUM ROBOTS AND CONTROL THEREOF” and filed Jul. 8, 2010. The above-identified documents are each herein incorporated by reference in their entirety.

[0002] This invention was made with government support under contract/grant no. IIS-I054331 awarded by National Science Foundation. The government has certain rights in the invention

[0003] The present invention relates to flexible instruments, and more specifically to apparatus and methods for configuring and operating such flexible instruments.

[0004] Continuum robots offer a number of potential advantages over traditional rigid link robots in certain applications, particularly those involving reaching through complex trajectories in cluttered environments or where the robot must compliantly contact the environment along its length. The inherent flexibility of continuum robots makes them gentle to the environment, able to achieve whole arm manipulation, and gives rise to a unique form of dexterity the shape of the robot is a product of both actuator and externally applied forces and moments. However, existing designs for continuum robots provide a limited range of motion and thus limit their usefulness.

[0035] The present invention is described with reference to the attached figures, wherein like reference numerals are used throughout the figures to designate similar or equivalent elements. The figures are not drawn to scale and they are provided merely to illustrate the instant invention. Several aspects of the invention are described below with reference to example applications for illustration. It should be understood that numerous specific details, relationships, and methods are set forth to provide a full understanding of the invention. One having ordinary skill in the relevant art, however, will readily recognize that the invention can be practiced without one or more of the specific details or with other methods. In other instances, well-known structures or operations are not shown in detail to avoid obscuring the invention. The present invention is not limited by the illustrated ordering of acts or events, as some acts may occur in different orders and/or concurrently with other acts or events. Furthermore, not all illustrated acts or events are required to implement a methodology in accordance with the present invention.

[0036] The various embodiments of the invention provide systems and methods for the control of continuum robots and similar flexible devices and instruments. Although, the discussion of the various embodiments will be in terms of wire-actuated (“tendons”) or rod-actuated continuum robots, this is solely for ease of illustration. Accordingly, the various aspects of the invention can be utilized with any type of wire-actuated or rod-actuated device or instrument.

[0037] One aspect of the various embodiments is directed to deformation or manipulation of some type of elastic structure using wires or string under tension (i.e., tendons or other elements with a lower stiffness than the elastic structure. These embodiments of the invention extend previous work on the Cosserat rod-based approach by taking into account not only the attachment point moment, but also the attachment point force and the distributed wrench that the tendon applies along the length of the elastic member. This approach couples the classical Cosserat string and rod models to express tendon loads in terms of the rod's kinematic variables.

[0038] The difference between this new coupled model and the point moment model for out of plane loads is shown in FIG. 1, and provide an experimental comparison of the two models described below. FIG. 1 is an illustration of the results of simulations of a continuum robot with single, straight, tensioned tendons with in-plane and out-of-plane forces applied at the tip. These illustrate the difference between the model proposed in this paper which includes distributed tendon wrenches, and the commonly used point moment approximation. For planar deformations and loads, the two models differ only by axial compression (which is small in most cases). However, for out of plane loads, the results differ significantly and including distributed wrenches enhances model accuracy (see Sec. V).

[0039] The various embodiments thus provide two new innovations over conventional methods. First, a new Cosserat rod-based model is provided for the spatial deformation of tendon actuated continuum robots under general external point and distributed wrench loads. This model is the first to treat the full effects of all of the tendon loads in a geometrically exact way for large 3D deflections. Second, the new model is the first to describe the mechanics of general tendon routing paths that need not run straight (along the undeformed robot configuration), as has been the case in prior prototypes. Thus, by providing a general model that can address most, if not all, types of tendon routing, this expands the design space and the set of shapes achievable for tendon-actuated robots.

[0040] In view of the foregoing, the various embodiments provide systems and methods for controlling continuum robots using exact models for the forward kinematics, statics, and dynamics and with general tendon routing experiencing external point and distributed loads. The models account for lance deformations due to bending, torsion, shear, and elongation. The static model is formulated as a set of nonlinear differential equations in standard form, and the dynamic model consists of a system of hyperbolic partial differential equations.

[0041] Using this approach, one can accurately predict the shape of a physical prototype with both straight and non-straight tendon routing paths and with external loading. With calibrated parameters, the mean tip error with respect to the total robot length can be significantly reduced as compared to conventional methods.

[0042] As illustrated in FIGS. 2A and 2B, the design space of achievable robot shapes can be expanded by considering alternative tendon paths. For example, FIG. 2A is an illustration of an exemplary of robot shape/workspace modification for a robot with four straight tendons spaced at equal angles around its periphery. In contrast, FIG. 2B is an illustration of an exemplary of robot shape/workspace modification for a robot with four helical tendons that each make one full revolution around the shaft, in accordance with an embodiment of the invention. The two designs differ significantly in tip orientation capability, and the helical design may be better suited to some types of tasks, e.g. a planar industrial pick and place tasks or surgical tasks. Thus, such continuum robots could be used to enhance the capabilities of medical instruments introduced using an orifice or minimal incision and allow additional control during procedures. Some types of devices can include devices for procedures in the throat and airways (introduced via the mouth), in the colon (introduced via the anus), in the stomach (introduced via the mouth and traveling through the esophogus), in the abdomen (either via a transgastric natural orifice transluminal endoscopic surgery approach, or via an incision in the abdomen similar to normal laparoscopic surgery), in skull base surgery (entering via the nose), in the brain and subarachnoid (entering an area around the brain stem and center of the spine) spaces (entering via a craniotomy), or in the bladder and kidneys (entering via the urethra).

[0043] The models in accordance with the various embodiments of the invention therefore allow new quasi-static and/or dynamic control techniques for tendon-actuated continuum robots in the future. Furthermore, the inclusion of general external loads in tendon actuated continuum robot models is an important step forward for future practical applications, given their significant sag under self-weight and when carrying payloads. Additionally, such models can be used to address the issue of modeling static friction, and real-time computation of static and dynamic robot shape.

[0044] A. Elastic Structure Kinematics

[0045] In Cosserat-rod theory, a rod or other elastic structure is characterized by its centerline curve in space p(s)ε and its material orientation, R(s)εSO(3) as functions of a reference parameter sε[0 L]|. Thus a homogeneous transformation can be used to describe the entire rod:

$g\left(s\right)=\left[\begin{array}{cc}R\left(s\right)& p\left(s\right)\\ 0& 1\end{array}\right]$

[0046] Kinematic variables v(s) and u(s) represent the linear and angular rates of change of g(s) with respect to s expressed in coordinates of the “body frame” g(s). Thus, the evolution of g(s) along s is defined by the following relationships;

{dot over (r)}(s)=r(s){circumflex over (u)}(s),

{dot over (p)}(s)=r(s)v(s)  (1)

where, the dot denotes a derivative with respect to s, and the and operators are as defined by R. M. Murray, Z. Li, and S. S. Sastry in “A Mathematical Introduction to Robotic Manipulation.” Boca Raton, Fla.: CRC Press, 1994. See also the Derivation Appendix for an explanation of these operators.

[0047] Letting the undeformed reference configuration of the rod be g*(s), where the z axis of R*(s) is chosen to be tangent to the curve p*(s). One could use the Frenet-Serret or Bishop's convention to define the x and y axes of R*(s), or, if the elastic structure has a cross section which is not radially symmetric, it is convenient to make the x and y axes align with the principal axes. The reference kinematic variables v* and u* can then be obtained by

[v*^Tu*^T]^T=(g*⁻¹(s)ĝ*(s)){hacek over ( )},

[0048] If the reference configuration happens to be a straight cylindrical elastic structure with s as the arc length along it, then v* =[0 0 1]^T and u*(s)=[0 0 0]^T.

[0049] B. Equilibrium Equations

[0050] One can the write the equations of static equilibrium for an arbitrary section of rod as shown in FIG. 3. The internal force and moment vectors (in global frame coordinates) are denoted by n and m, the applied force distribution per unit of s is f, and the applied moment distribution per unit of s is l. Taking the derivative of the static equilibrium conditions with respect to s, one arrives at the classic forms of the equilibrium differential equations for a special Cosserat rod,

{dot over (n)}(s)+f(s)=0,  (2)

{dot over (m)}(s)+{dot over (p)}(s)×n(s)+l(s)=0.  (3)

[0051] C. Constitutive Laws

[0052] The difference between the kinematic variables in the elastic structure's reference state and those in the deformed state can he directly related to various mechanical strains. For instance, transverse shear strains in the body-frame x and y directions correspond to v_x−v_x* and v_y−v_y* respectively, while axial elongation or stretch in the body-frame z direction corresponds to v_z−v_z*. Similarly, bending strains about the local x and y axes are related to u_x−u_x* and u_y−u_y* respectively, while torsional strain about the local z axis is related to u_z−u_z*.

[0053] One can use linear constitutive laws to map these strain variables to the internal forces and moments. Assuming that the x and y axes of g* are aligned with the principal axes of the cross section, one obtains

n(s)=R(s)D(s)(v(s)−v*(s)),

m(s)=R(s)C(s)(u(s)−u*(s)),  (4)

where

D(s)=diag(GA(s),GA(s),EA(s)), and

C(s)=diag(EI_xx(s),EI_yy(s),EI_xx(s)+EI_yy(s)),

where A(s) is the area of the cross section, E(s) is Young's modulus, G(s) is the shear modulus, and I_xx(s) and I_yy(s) are the second moments of area of the tube cross section about the principal axes. (Note that I_xx(s)+I_yy(s) is the polar moment of inertia about the centroid.) One can use these linear relationships here because they are notationally convenient and accurate for many continuum robots, but the Cosserat rod approach does not require it.

[0054] D. Explicit Model Equations

[0055] Equations (2) and (3) can then be written in terms of the kinematic variables using equation (4), their derivatives, and equation (1). This leads to the full set of differential equations shown below.

{dot over (p)}=Rv

{dot over (R)}=Rû

{dot over (v)}={dot over (v)}*−D⁻¹((ûD+{dot over (D)})(v−v*)+R^Tf)

{dot over (u)}={dot over (u)}*−C⁻¹((ûC+Ċ)(u−u*)+{circumflex over (v)}D(v−v*)+R^Tl)  (5)

[0056] Alternatively, an equivalent system can be obtained using m and n as state variables rather than v and u.

{dot over (p)}=R(D⁻¹R^Tn+v*)

{dot over (R)}=R(C⁻¹R^Tm+u*){tilde over ( )}

{dot over (n)}=−f

{dot over (m)}=−{dot over (p)}×n−l  (6)

[0057] Boundary conditions for a rod which is clamped at s=0 and subject to an applied force F_l and moment L_l at s=l would be R(0)=R₀, p(0)=p₀, m(l)=L_l, and n(l)=F_l.

[0058] Having reviewed the classic Cosserat-rod model, the derivation a new model for tendon driven continuum manipulators in accordance with the various embodiments of the invention will now be presented. The derivation uses the Cosserat model of Section II to describe the elastic member and the classic Cosserat model for extensible strings to describe the tendons. For purposes of the model, the string and elastic structure models are coupled together by deriving the distributed loads that the tendons apply to the elastic structure in terms of the elastic structure's kinematic variables, and then incorporating these loads into the model.

[0059] A. Assumptions

[0060] Two standard assumptions are employed in the derivation. First, an assumption of frictionless interaction between the tendons and the channel through which they travel. This implies that the tension is constant along the length of the tendon. Frictional forces are expected to increase as the curvature of the robot increases due to larger normal forces, but the assumption of zero friction is valid if low friction materials are used, which is the case for the experimental prototype discussed below. Second, the locations of the tendons within the cross section of the robot are assumed not to change during the deformation. This assumption is valid for designs which use embedded sleeves or channels with tight tolerances, as well as designs which use closely spaced tendon guide portions.

[0061] B. Tendon Kinematics

[0062] One can separate the terms f and l in the equations in (5) into truly external distributed loads, f_e and l_e, and distributed loads due to tendon tension, f_t and l_t.

f=f_e+f_t

l=l_e+l_t.  (7)

[0063] In order to derive f_t and l_t, one starts by defining the path in which the tendon is routed along the robot length. Note that this path can be defined by channels or tubes within the elastic structure. The elastic structure can be a homogeneous elastic structure or a series of support disks on an elastic member—both of which afford considerable flexibility in choosing tendon routing. In the experimental prototype, many holes are drilled around the periphery of each support disk, allowing easy reconfiguration of tendon path as desired.

[0064] A convenient way to mathematically describe the tendon routing path is to define the tendon location within the robot cross section as a function of the reference parameter s. Thus, the routing path of the i^th tendon is defined by two functions x_i(s) and y_i(s) that give the body-frame coordinates of the tendon as it crosses the x-y plane of the attached elastic structure frame at s. As shown in FIG. 4, a vector from the origin of the attached frame to the tendon location is then given in attached frame coordinates by

r_i(s)=[x_i(s)y_i(s)0]^T,  (8)

The parametric space curve defining the tendon path in the global frame when the robot is in its undeformed reference configuration is then given by

p_i*(s)=R*(s)r_i(s)+p*(s).

Similarly, when the robot is deformed due to tendon tension or external loads, the new tendon space curve will be

p_i(s)=R(s)r_i(s)+p(s).  (9)

[0065] C. Distributed Forces on Tendons

[0066] The governing differential equations for an extensible string can be derived by taking the derivative of the static equilibrium conditions for a finite section. This results in the same equation for the internal force derivative as in equation (2).

{dot over (n)}_i(s)+f_i(s)=0,  (10)

where f_i(s) is the distributed force applied to the i^th tendon per unit of s, and n_i(s) is the internal force in the tendon. In contrast to a Cosserat rod, an ideal string has the defining constitutive property of being perfectly flexible, meaning it cannot support internal moments or shear forces, but only tension which is denoted by τ_i. This requires that the internal force be always tangent to the curve p_i(s). Thus, one can write

$\begin{array}{cc}{n}_i\left(s\right)={\tau }_i\frac{{\stackrel{.}{p}}_i\left(s\right)}{{\stackrel{.}{p}}_i\left(s\right)}.& \left(11\right)\end{array}$

If friction were present, τ_i would vary with s, but under the frictionless assumption, it is constant along the length of the tendon. Using (10) and (11) one can derive the following expression for the distributed force on the tendon (see Appendix for Derivation):

$\begin{array}{cc}{f}_i\left(s\right)=-{\stackrel{.}{n}}_i={\tau }_i\frac{{\widehat{\stackrel{.}{p}}}_i^2}{{{\stackrel{.}{p}}_i}^2}{\stackrel{\_}{p}}_i.& \left(12\right)\end{array}$

[0067] D. Tendon Loads on Elastic Structure

[0068] One can now write the collective distributed loads f_t and l_t that the tendons apply to the elastic member, in terms of the individual forces on the tendons and their locations in the elastic structure cross-section (or in the case of guide discs, the cross-section of the elastic member supporting the discs). The total distributed force is equal and opposite to the sum of the individual force distributions on the tendons shown in equation (12), namely,

$f_i=-\sum _{i=1}^nf_i.$

The distributed moment at the elastic member centroid is the sum of the cross products of each moment arm with each force. Thus,

$l_i=-\sum _{i=1}^n{\left(p_i-p\right)}^{^}f_i=-\sum _{i=1}^n{\left({\mathrm{Rr}}_i\right)}^{^}f_i.$

Substituting equation (12), yields

$\begin{array}{cc}{f}_i=-\sum _{i=1}^n{\tau }_i\frac{{\stackrel{\widehat{.}}{p}}_i^2}{{{\hat{p}}_i}^3}{\hat{p}}_i,l_i=-\sum _{i=1}^n{{\tau }_i\left({\mathrm{Rr}}_i\right)}^{^}\frac{{\stackrel{\widehat{.}}{p}}_i^2}{{{\hat{p}}_i}^3}{\hat{p}}_i.& \left(13\right)\end{array}$

One can then express these total force and moment distributions in terms of the kinematic variables u, v, R and p so that one can substitute them into equations (7) and (5). To do this, one expands {dot over (p)} and {umlaut over (p)}. Differentiating equation (9) twice yields,

{dot over (p)}_i=R(ûr_i+{dot over (r)}_i+v),

{umlaut over (p)}_i=R({circumflex over (u)}(ûr_i+{dot over (r)}_i+v)+{dot over ({circumflex over (u)}r_i+û{dot over (r)}_i+{umlaut over (r)}_i+{dot over (v)}).  (14)

[0069] It is noted that {umlaut over (p)} is a function of {dot over (u)} and {dot over (v)}. Therefore, substituting these results into equation (13), and equation (13) into the rod model equation (5) via equation (7), one can obtain an implicitly defined set of differential equations. Fortunately, the resulting equations are linear in u and v, and it is therefore possible to manipulate them into an explicit form. Rewriting them in this way (such that they are amenable to standard numerical methods) is the topic of the following subsection.

[0070] E. Explicit Decoupled Model Equations

[0071] The coupled elastic structure & tendon model is given in implicit form by equations (5), (7), (13). and (14). In this subsection, these implicit equations are manipulated into explicit, first-order, state-vector form. To express the result concisely, some intermediate matrix and vector quantities are defined, starting with equation (14) expressed in body-frame coordinates, i.e.

{dot over (p)}_i^b=ûr_i+{dot over (r)}_i+v,

{umlaut over (p)}_i^b=û{dot over (p)}_i^b+{dot over ({circumflex over (u)}r_i+û{dot over (r)}_i+{umlaut over (r)}_i+{dot over (v)}.

[0072] Now define Matrices A_i, A, B_i, and B, as well as vectors a_i, a, b_i, and b, as follows:

$A_i=-r_i\frac{{\left({\stackrel{\widehat{.}}{p}}_i^b\right)}^2}{{{\stackrel{.}{p}}_i^b}^3},B_i={\hat{r}}_iA_i$ $A=\sum _{i=1}^nA_i,B=\sum _{i=1}^nB_i,a_i=A_i\left(\hat{u}{\stackrel{.}{p}}_i^b+\hat{u}{\stackrel{.}{r}}_i+{\hat{r}}_i\right),b_i={\hat{r}}_ia_i$ $a=\sum _{i=1}^na_i,b=\sum _{i=1}^nb_i,$

to find that f_t and l_t can now be expressed as

$\begin{array}{cc}{f}_i=R\left(a+A\stackrel{.}{\upsilon }+\sum _{i=1}^nA_i\stackrel{\widehat{.}}{u}r_i\right),l_i=R\left(b+B\stackrel{.}{\upsilon }+\sum _{i=1}^nB_i\stackrel{\widehat{.}}{u}r_i\right).& \left(15\right)\end{array}$

[0073] The vector terms Σ_i=1ⁿA_i{dot over (û)}r_i and Σ_i=1ⁿB_i{dot over (û)}r_i are both linear in the elements of {dot over (u)}. Therefore, it is possible to express them both by equivalent linear operations on {dot over (u)}. That is, one can define matrices G and H as

$G=\sum _{i=1}^n\left[\begin{array}{ccc}{A}_i{\hat{e}}_1r_1& A_i{\hat{e}}_2r_i& A_i{\hat{e}}_3r_i\end{array}\right]$ $H=\sum _{i=1}^n\left[\begin{array}{ccc}{B}_i{\hat{e}}_1r_1& B_i{\hat{e}}_2r_i& B_i{\hat{e}}_3r_i\end{array}\right]$

where e₁, e₂, and e₃ are the standard basis vectors [1 0 0], [0 1 0], and [0 0 1]. Then, equation (15) becomes

f_i=R(a+A{dot over (v)}+G{dot over (u)}),

ι_t=R(b+B{dot over (v)}+H{dot over (u)}),

[0074] Substituting tendon load expressions into the last two equations in (5) and rearranging them provides

(D+A){dot over (v)}+G{dot over (u)}=d

B{dot over (v)}+(C+H){dot over (u)}=c

where the vectors c and d are functions of the state variables as shown below.

d=D{dot over (v)}*−(ûD+{circumflex over (D)})(v−v*)−R^Tf_e−a

c=C{dot over (u)}*−(ûC+Ĉ)(u−u*)−{circumflex over (v)}D(v−v*)−R^Tl_e−b.

One can now easily write the governing equations as

$\begin{array}{cc}\stackrel{.}{p}=R\upsilon \stackrel{.}{R}=R\hat{u}\left[\begin{array}{c}\stackrel{.}{\upsilon }\\ \stackrel{.}{u}\end{array}\right]={\left[\begin{array}{cc}D+A& G\\ B& C+H\end{array}\right]}^{-1}\left[\begin{array}{c}d\\ c\end{array}\right].& \left(17\right)\end{array}$

Noting that the quantities on the right hand side of equation (17) are merely functions of the state variables and system inputs (u, R, τ_n, f_e and l_e) one arrives at a system of differential equations in standard explicit form, describing the shape of a continuum robot with any number of generally routed tendons and with general external loads applied.

[0075] This system can be solved by any standard numerical integration routine for systems of the form {dot over (y)}=f(s,y). The required matrix inverse may be calculated (either numerically or by obtaining a closed form inverse) at every integration step, or one could alternatively rewrite the equations as a system with a state dependent mass matrix on the left hand side and use any standard numerical method for solving M(y,s){dot over (y)}′=f(s,y). For purposes of the simulations and experiments in accordance with the various embodiments of the invention, numerically inversion is used.

[0076] F. Boundary Conditions

[0077] When tendon i terminates at s=l_i along the length of the robot, it applies a point force to its attachment point equal and opposite to the internal force in the tendon given by equation (11). Thus, the point force vector is given by

$\begin{array}{cc}{F}_i=-n_i\left({}_i\right)=-{\tau }_i\frac{{\stackrel{.}{p}}_i\left({}_i\right)}{{\stackrel{.}{p}}_i\left({}_i\right)}.& \left(18\right)\end{array}$

[0078] With a moment arm of p_i(l_i)−p(l_i), this force creates a point moment L_i at the elastic member centroid of,

$\begin{array}{cc}{L}_i=-{{\tau }_i\left(R\left({}_i\right)r_i\left({}_i\right)\right)}^{^}\frac{{\stackrel{.}{p}}_i\left({}_i\right)}{{\stackrel{.}{p}}_i\left({}_i\right)}.& \left(19\right)\end{array}$

If at some location s=σ, point loads F(σ) and L(σ) (resulting from tendon terminations or external loads) are applied to the elastic structure, the internal force and moment change across the boundary s=σ by,

n(σ⁻)=n(σ⁺)+F(σ),

m(σ⁻)=m(σ⁺)+L(σ),  (20)

where σ⁻ and σ⁺ denote locations just before and just after s=σ. Any combination of external point loads and tendon termination loads can be accommodated in this way.

[0079] G. Point Moment Model

[0080] In prior tendon robot models, tendon actuation has often been modeled by simply applying the pure point moment in equation (19) to an elastic structure model at the location where each tendon is attached, without considering the point force at the attachment point and the distributed tendon loads along the length (see FIG. 6). This approach is convenient because it allows one to use the classical Cosserat rod equations by simply applying boundary conditions that take into account the tendon termination moments.

[0081] This approximation for planar robots is justified since the effects of the point force and the distributed loads effectively “cancel” each other, leaving only the point moment. Thus, as shown in FIG. 1 this approach yields almost exactly the same final shape as the full coupled model when the robot deformation occurs in a plane.

[0082] However, as shown in FIG. 1, the two approaches diverge as the robot shape becomes increasingly non-planar due to a transverse load at the tip. In Section V, an investigation of the accuracy of both approaches a set of experiments on a prototype robot is provided.

[0083] Based on the coupled elastic structure and tendon model presented above for static continuum robot deformations, a model for the dynamics of a continuum robot with general tendon routing is derived. Such a model will be useful for analyzing the characteristics of specific designs as well as the development of control algorithms similar to those derived for planar robots with straight tendons. As shown below, adding the necessary dynamic terms and equations results in a hyperbolic system of partial differential equations, which can be expressed in the standard form

y_t=f(s,t,y,y_s),  (21)

where a subscript s or t is used in this section to denote partial derivatives with respect to the reference parameter s and time t respectively.

[0084] Two new vector variables are introduced, q and ω, which are the body frame linear and angular velocity of the rod at s. These are analogous to u and v respectively, but are defined with respect to time instead of arc length. Thus,

p_t=RqR_t=R{circumflex over (ω)}.  (22)

Recalling from equations in (5) that

p_s=RvR_s=Rû,  (23)

and using the fact that p_st=p_ts and R_st=R_ts one can derive the following compatibility equations,

u_t=w_s+û_wv_t=q_s+ûq−ŵv,  (24)

[0085] Equations (2) and (3) describe the static equilibrium of the rod. To describe dynamics, one can add the time derivatives of the linear and angular momentum per unit length in place of the zero on the right hand side, such that they become,

{dot over (n)}+f=ρAp_tt,  (25)

{dot over (m)}+{dot over (p)}×n+l=θ_t(RρJw),  (26)

where ρ is the mass density of the elastic structure, A is the cross sectional area of the elastic structure, and J is the matrix of second area moments of the cross section. Expanding these and applying the equations in (24) one can obtain a complete system in the form of equation (21),

$\begin{array}{cc}p_t=\mathrm{Rq}R_t=R\hat{\omega }{\upsilon }_i=q_s+\hat{u}q-\hat{\omega }\upsilon u_i={\omega }_s+\hat{u}\omega q_t=\frac{1}{\rho A}\left(D\left({\upsilon }_s-{\upsilon }_s^{*}\right)+\left(\hat{u}D+D_s\right)\left(\upsilon -{\upsilon }^{*}\right)+R^T\left(f_c+f_t\right)-\rho A\hat{\omega }q\right){\omega }_t={\left(\rho J\right)}^{-1}\left(C\left(u_s-u_s^{*}\right)+\left(\hat{u}C+C_s\right)\left(u-u^{*}\right)+\hat{\upsilon }D\left(\upsilon -{\upsilon }^{*}\right)+R^T\left(l_s+l_t\right)-\hat{\omega }\rho J\omega \right)& \left(27\right)\end{array}$

where f_t and l_t can be computed using the equations in (16). Typically, conditions at t=0 are given for all variables along the length of the robot, and the boundary conditions of Subsection IMF apply for all times.

[0086] A. Dynamic Simulation

[0087] To illustrate the capability of the equations in (27) to describe the time evolution of the shape of a continuum robot with general tendon routing, the following dynamic simulation of a robot whose elastic structure is identical to that of the experimental prototype described below is provided. The robot contains a single tendon routed in a helical where the tendon makes one complete revolution around the shaft as it passes from the base to the tip. This routing path is the same as the one for tendon 5 in the prototype, which is specified in Table 1.

TABLE I
TENDON ROUTING PATHS USED IN EXPERIMENTS
Tendon (i)	1	2	3	4	5	6
xi(s) (mm)	8	0	−8	0	8cos(2πs/l)	refer to (28)
yi(s) (mm)	0	8	0	−8	8sin(2πs/l)	refer to (28)

[0088] FIG. 7 shows snapshots of the robot backbone shape at millisecond intervals after a step input of 5 Newtons of tendon tension was applied. For the numerical simulation Richtmyer's two-step variant of the Lax-Wendroff finite-difference scheme was implemented.

[0089] The maximum length of the time step for any explicit time-marching algorithm for hyperbolic partial differential equations is limited by the Courant-Friedriechs-Lewy condition for stability. This is a fairly restrictive condition for dynamic rod problems because the shear, extension, and torsional vibrations are so fast that a very small is required in order to capture them without the simulation becoming unstable. An active research field in mechanics and computer graphics simulation is to find reduced—order models of elastic structures that are physically accurate and yet capable of being simulated in real-time. This simulation confirms the intuition that the elastic structure should move towards a helical very when the helical tendon undergoes a step in tension.

[0090] Below are described several different experiments conducted using a continuum robot prototype with a variety of tendon paths and external loading conditions applied

[0091] A. Prototype Constructions

[0092] A prototype robot in accordance with the various embodiments is shown in FIG. 8. The elastic structure consists of a central elastic member 802 is a spring steel rod (ASTM A228) of length l=242 mm and diameter d=0.8 mm with tendon guide portions 804 consisting of 12 stand-off disks, 20 mm in diameter, spaced 20 mm apart along its length. The disks were laser cut from 1.57 um thick PTFE filled Delrin plastic to minimize friction with the tendons. As shown in the inset of FIG. 8, 24 small pass-through holes 806 were laser cut in a circular pattern at a radius of 8 mm from the center of each disk. The elastic member rod 802 was passed through the center holes 808 of the disks and each was fixed to it using Loctite 401. For tendons 810, 0.36 mm PTFE coated fiberglass thread were used. Each tendon 810 was run through various pass-through holes along the robot and knotted at the end, after passing through the final support disk. The optimal ratio of tendon support spacing to offset distance from the elastic member was found to be 0.4, and the prototype was designed to exactly match this ratio.

[0093] Although the exemplary robot configuration utilizes standoff disks to provide the tendon guide portions, the various embodiments are not limited in this regard. Rather, any of means of coupling the tendons to the elastic member to cause deformation of the elastic member can be used in the various embodiments. Further, a particular combination of materials, spacing of guide portions, and openings in the guide portions is provided, the various embodiments are not limited in this regard. Rather, any variations on the combination recited above can be used with the various embodiments. Additionally, the methods above can be used with any number of tendons. In such embodiments, the tendons can extend along a same portion of the length of the elastic member or the tendons can extend over different portions of the length of the elastic member, including overlapping portions.

[0094] The tendon routing paths can be reconfigured on this robot by “re-threading” the tendons through a different set of holes in the various support disks. The robot's self-weight distribution was measured to be 0.47 N/m, which is enough to cause significant deformation, producing 44 mm of downward deflection at the tip (18% of the total arc length) for zero tendon tension. This weight was incorporated into all model calculations as a distributed force.

[0095] B. Experimental Procedure

[0096] In each of the following experiments, known tensions were applied to tendons behind the base of the robot by passing the tendons over approximately frictionless pulleys and attaching them to hanging calibration weights. In those cases with applied point loads, weights 812 were also hung from the tip of the robot, as shown in FIG. 8.

[0097] In each experiment, a set of 3D elastic member points was collected by manually touching the elastic member with the tip of an optically tracked stylus as shown in FIG. 8. A MicronTracker 2 H3-60 (Claron Technology, Inc.) was used to track the stylus, which has a specified fiducial measurement accuracy of 0.20 mm.

[0098] C. Calibration

[0099] The base frame position of the robot can be determined accurately using the optically tracked stylus. The angular orientation of the robot elastic member as it leaves the base support plate is more challenging to measure (Note that the elastic member cannot be assumed to exit exactly normal to the plate due to the tolerance between the elastic member and the hole drilled in the plate, and a 2° angular error in base frame corresponds to an approximately 8 mm tip error when the robot is straight). Also, the effective stiffness of the elastic member was increased due to the constraints of the standoff disks and Loctite adhesive at regular intervals. To account for these uncertainties the effective Young's modulus and the set of XYZ Euler angles (α, B, and γ) describing the orientation of the base frame were calibrated.

[0100] The calibration process was accomplished by sorting a preconstrained nonlinear optimization problem to find the set of parameters which minimizes the sum of the positional errors at the tip of the device for the set of 25 experiments with straight tendon paths described in Sec. V-D and Table II.

TABLE II
EXPERIMENTAL TENSIONS AND TIP LOADS
Experimental with Tendons 1-4 (Straight)
Tension (N)	0	0.98	1.96	2.94	2.94	2.94	4.91
Tip Load (N)	0	0	0	0	0.098	0.196	0
Experiments with Tendon 5 (Helical)
Tension (N)	0.98	1.96	2.94	4.91	4.91	4.91	6.87
Tip Load (N)	0	0	0	0	0.098	0.196	0
Experiments with Tendon 6 (Polynomial)
	Tension (N)	1.50	2.46	3.66	4.91	4.91
	Tip Load (N)	0	0	0	0	0.0196

[0101] In other words, for the parameter set P={E, α, B, γ}:

$P_{\mathrm{ca}l}=\underset{P}{\mathrm{argmin}}\left(\sum _{i=1}^{25}e_k\right)$ $\mathrm{where}e_k={p_{\mathrm{model}}\left(\right)-p_{\mathrm{date}}\left(\right)}_k$

is the Euclidean distance between the model tip prediction and the data in experiment k. To implement this minimization, the Nelder-Meade simplex algorithm was used.

[0102] To ensure fair comparison of the coupled model and the point moment model, the calibration procedure was performed separately for each model. Results are shown in Table III.

TABLE III
NOMINAL AND CALIBRATED PARAMETERS
	Nominal	Calibrated Value	Calibrated Value
Parameter	Value	(Point Moment Model)	(Coupled Model)
E (GPa)	210	227.9	229.6
α (deg)	180	177.7	177.9
β (deg)	0	2.2	2.2
γ (deg)	−90	−89.6	−89.7

Note that the similarity in calibrated Euler angles and their low deviation from nominal provides confidence that the correct base frame was obtained for both models. It is also important to note that the models contain the same number of parameters, so a fair comparison can he made, As expected, the calibrated values for Young's modulus are higher than the nominal value of 210 GPa for spring steel, due to the increased stiffness provided by the disks and glue. Poisson's ratio was held constant at v=0.3125 during calibration so that the shear modulus was correctly scaled relative to Young's modulus.

[0103] D. Straight Tendon Results and Model Comparison

[0104] Table I details the location of the tendon routing paths used in the experiments in terms of x_i(s) and y_i(s) as defined in (8). Twenty-five (25) experiments were performed (detailed in Table II) with straight tendon paths in order to compare the accuracy of the new coupled model with that of the point moment model. The tip error statistics for both models with calibrated parameters is detailed in Table IV.

TABLE IV
MODEL TIP ERRORS FOR STRAIGHT TENDON EXPERIMENTS
	Tip Error Statistic (mm)	mean	std. dev.	min	max
13 Cases with In-Plane Loads
	Point Moment Model	3.5	1.4	1.2	5.6
	Coupled Model	3.1	1.3	0.3	5.3
12 Cases with Out-of-Plane Loads
	Point Moment Model	9.8	5.5	1.7	16.2
	Coupled Model	4.1	2.1	0.6	7.9

The results for in-plane loading are accurate for both models, as shown in FIG. 9. FIG. 9 is an x-y-z plot of the actual and simulated result of in-plane loading of an exemplary continuum robot using a straight tendon. Shown in FIG. 9 are the 13 experimental cases with in-plane external loads. The tendons on the top and bottom of the robot (tendons 1 and 3) were tensioned and vertical tip loads were applied in four of the cases. Distributed gravitational loading is present in every case. As detailed in Table IV, both the coupled model and the point moment model are accurate and nearly identical for in-plane loads. In contrast, for out-of-plane loads, the coupled model provides more accurate predictions, as shown in FIG. 10. FIG. 10 is an x-y-z plot of the actual and simulated result of out-of-plane loading of an exemplary continuum robot using a straight tendon. Pictured in FIG. 10 are the twelve experimental cases with out-of-plane external loads. The tendons on the left and right of the robot (tendons 2 and 4) were tensioned. (a) Distributed loading (robot self-weight) applied, (b) additional tip loads applied. As detailed in Table IV, the data agrees with the coupled model prediction, but the point moment model becomes inaccurate as the out-of-plane load increases, and as the curvature increases.

[0105] With calibrated parameters, the mean tip error over all 25 straight tendon experiments was 3.6 mm for the coupled model. This corresponds to 1.5% of the total arc length of the robot. Note that experimental data points lie close to the model prediction along the entire robot length, and the error increases gradually along the robot length, so that tip error normalized by the robot length is a reasonable metric for the accuracy of the model.

[0106] E. A High-Tension, Large-Load, Straight Tendon Experiment

[0107] Also performed was one additional straight tendon experiment to see how the two approaches compare for a case of large tension and large out-of-plane load, similar to the case which is simulated in FIG. 1. Tendon 4 was tensioned to 6.38 N and a downward tip force of 0.196 N was applied. The resulting data and model predictions are shown in FIG. 11. FIG. 11 is an x-y-z plot of the actual and simulated result of out-of-plane loading for an exemplary continuum robot using a straight tendon and high tension. As illustrated in FIG. 11, the two models produce very different results. FIG. 11 shows that the coupled model prediction lies much closer to the data. Here, the tip error of the point moment model is 57 mm (23.5% of robot length), while the coupled model tip error is 12.8 mm (5.3% of robot length).

[0108] F. Experiments with Helical Tendon Routing

[0109] To explore more complex tendon routing, helical routing paths were also evaluated. As given in Table I, the helical routing path winds through one complete revolution as it traverses the robot from base to tip. The tensions and tip loads for these experiments are detailed in Table II. Using the parameters calibrated from the previous straight tendon dataset, the resulting data and model predictions are plotted in FIGS. 12A and 12B. FIG. 12A is an x-y-z plot of the actual and simulated results of operation of an exemplary continuum robot using a helical tendon without a load. FIG. 12B is an x-y-z plot of the actual and simulated results of operation of an exemplary continuum robot using a helical tendon with a tip load. As seen from Table V, the model agrees with the data with a mean tip error of 5.5 mm. The small increase in error over the straight tendon cases may be due to increased frictional forces since the tension for the helical cases was higher.

TABLE V
COUPLED MODEL TIP ERRORS FOR NON-STRAIGHT TENDON
EXPERIMENTS
	mean	std. dev.	min	max
	Tendon 5 (Helical)	5.5	2.7	1.9	10.0
	Tendon 6 (Polynomial)	4.6	1.9	2.7	7.2

[0110] G. Experiments with Polynomial Tendon Routing

[0111] In order to further illustrate the model's generality, an additional experiment with a general curved tendon routing choice was performed. In particular, the routing path variables were parameterized by two trigonometric functions whose arguments are defined by a polynomial function of degree 4 as follows:

x₆(s)=8 cos(5887s⁴−2849s³+320s²+6s)

y₆(s)=8 sin(5887s⁴−2849s³+320s²+6s),  (28)

where s is in meters and x₆ and y₆ are in millimeters. This routing path starts at the top of the robot, wraps around to the right side for most of the length, and then returns to the top at the end of the robot. The tensions and loads are given in Table II, and the results are detailed in Table V and illustrated in FIG. 13. FIG. 13 is an x-y-z plot of the actual and simulated results of operation of an exemplary continuum robot using a polynomial tendon with tip loads according to Table II. The coupled model's predictions agree with the data, with a mean tip error of 4.6 mm. This set of experiments confirms the coupled model's ability to handle an arbitrary tendon routing choices.

[0112] H. Sources of Error

[0113] The largest source of measurement is likely the procedure of manually placing the tip of the on the robot during data capture. It is estimated that this uncertainty is at most 2 mm. In general, the largest model errors occurred when the tendons were under the greatest tension. This agrees with the intuition that effects of static friction should become more significant as the tension and curvature increase, However, the low overall errors suggest that neglecting static friction is justifiable for this prototype.

[0114] Accordingly, in the view of the foregoing, the equations above can be integrated into a continuum robot system, as shown in FIG. 14. FIG. 14 is a schematic illustration of a continuum robot system 1400 in accordance with the various embodiments. The system 1400 can include a continuum robot 1402 similar to that illustrated in FIG. 8. That is, the continuum robot 1402 can include an elastic member or backbone 1404 having one end or portion coupled to a support or base (not shown), a number of guide portions 1406, and at least one tendon 1408 extending through the guide portions and defining a tendon path, as described above. The system 1402 can further include an actuator/sensor 1410 for applying a force or tension and for ascertaining a current tension on the tendon 1408 to the tendon 1408. The system 1400 can also include a control system 1412 for operating the system 1400, which can include a computing device.

[0115] The system can have at least two modes of operation. In a first mode of operation, the actuator/sensor 1410 can generate signals indicative of a current tension on the tendon 1408. This signal can be recited by the control system 1412. The control system 1412 can then use the equations described above, particularly the governing equations at (17), to estimate a current or resulting shape of the member 1404. In particular, the governing equations at (17) can be solved to extract the shape of the member 1404. Additional sensors 1414, such as video sensors, can also be coupled to the control system 1412 to allow verification of this estimated shape. In a second mode of operation, the control system 1412 can also use the equations described above, particularly the governing equations at (17), to determine an amount of tension required for the member 1404 to achieve a desired shape. Thereafter, the control system 1412 can cause the actuator/sensor 1410 to adjust the tension on the tendon 1408. The additional sensors 141 can then be used to verify that the target shape has been achieved.

[0116] Using these two modes of operation, it is then possible to control the robot 1402 to perform various types of tasks, as the equations above allow one to detect and adjust the configuration of the robot 1402 in real-time based on measurement and adjustment of the tension of the tendon 1408. That is, a robot with increased dexterity can be provided. Such a robot can be useful for various applications. In particular, such robots would be useful for carry out procedures in confined spaces, as the increased dexterity would allow the user to maneuver the tip around obstructions in such spaces. For example, such robots could be used to reduce the invasiveness of some existing surgical procedures which currently cannot be performed using conventional robotic tools. Such procedures in transnasal skull base surgery, lung interventions, cochlear implantation procedures, to name a few. However, the various embodiments are not limited in this regard and the various methods and systems described herein can be used for any other procedure in which increased dexterity of the robot is desired or required.

[0117] In the discussion above, it is generally assumed that the stiffness of the tendons will be relatively lower than the stiffness of the elastic structure. However, in other embodiments, rod-type actuators, with stiffnesses greater than that of the elastic structure can be utilized. Note that the rods need not be completely stiff, but rather have a sufficient stiffness to apply a deformation force to the corresponding elastic structure. The advantage of such arrangements is that maintaining tension on the tendons is not required. Further, a different control methodology can be implemented by utilizing rod displacement to effect deformation. However, of particular advantage are several configurations.

[0118] In some configurations, a pre-curved rod can be utilized. As used herein, the term “pre-curved rod” refers to a rod that changes direction in two- or three-dimensional space at least once (i.e., non-linear or have at least one bend or change in direction in space) in the absence of an external force. Thus, in the absence of any external forces being applied to a pre-curved rod, the pre-curved rod is configured to revert to its shape with changes in direction.

[0119] Pre-curved rods are useful for providing extended capabilities to a wide variety of surgical tools, including, but not limited to standard endoscopes, robotic endoscopes, continuum robots. The extended capabilities include: (1) enabling such devices to reach around “tighter corners” (i.e., achieve higher curvatures), (2) achieve larger workspaces, and (3) achieve novel motions (e.g., spinning the tip in place without moving the tip's position). Given these properties, the various embodiments of the invention can be used for a wide variety of surgical applications throughout the body. One example is retroflexed endoscopy, in which the endoscope enters a lumen or cavity and must turn through an angle between 90 and 360 degrees to operate “backwards”, i.e., in a retroflexed manner, which is useful for both identifying and surgically removing polyps hidden by intestinal edges. Another example is Peroral Endoscopic Myotomy in which the endoscope must enter the stomach and then operate on the muscle of the sphincter through which the endoscope entered the stomach cavity.

[0120] Another advantage of pre-curved rods is that they provide greater dexterity in retroflexed configurations. Often, in a given procedure, an instrument may be capable of reaching a particular anatomical target, but just barely (i.e. the location is at the edge of the mechanism's workspace). This makes it very difficult for the surgeon to dexterously and effectively control the tip of the device. It is why retroflexed endoscopy, ERCP, and POEM are possible for expert surgeons, but highly challenging, and beyond the capabilities of all but the best surgeons. With the pre-curved strip-actuated design, the overall workspace of the mechanism can be greater (with or without precurvature), and thus sites that previously lay near the edge of the workspace are now comfortably within the workspace, resulting in easier control and increased dexterity of the tip around those previously difficult-to-reach locations.

[0121] In one embodiment of the invention, a device with pre-curved rods can be configured similar to that shown in FIG. 8. As previously discussed, an elastic structure can be provided and can be configured for holding guide discs at a desired spacing between discs. The rods can then be used to push and pull on one or more discs of the device, causing the overall shape to change. Each rod is affixed to one disc and is free to slide within pass-through holes in other discs of the device.

[0122] As noted above, the rods can be pre-curved. For example, pre-curved, superelastic nitinol wires can be used. Thus the resulting shape of the device as a whole can be, for example, a circular arc when no forces are applied to the rods (i.e., when the device is in its “home” or undeflected position). FIG. 15 shows a straight rod device (top) and a pre-curved rod device (bottom) in the home position.

[0123] It should be noted that the home position of devices in accordance with the various embodiments is not limited solely to a circular arc configuration. Rather, the home position can be can any shape in three dimensional space. Thus, in a home position, a device in accordance with the various embodiments can include any number of bends and such bends can be circular and non-circular. In certain embodiments, the rod is not pre-curved and can be configured to have a low stiffness relative to the rods so as to conform to the home position of the rods, allowing the rod to sustain higher curvature changes than the rods.

[0124] In some embodiments, a rod can be configured have two degrees of freedom: (1) translated in/out (i.e., provide push/pull with respect to a disc) and (2) a rod can be rotated about its axis. Flexible surgical devices using straight wires or rods do not typically include a rotation feature. In particular, since rotation of a straight wire or rod would not result in a change in shape of the device, rotation would have no effect and thus, no use. In contrast, in the case of a pre-curved rod, the rotation of the rod alone can be utilized to change the shape of the device. Therefore, in the various embodiments, different combinations of rotation and translation for each pre-curved rod can permit a device with pre-curved rods in accordance with the various embodiments to produce shapes that would not be possible with straight rods or tendons. For example, non-circular curves. A mechanics-based model of the device, described in greater detail below, can then be utilized in selecting amounts of push/pull and rotation to produce a desired shape of the device.

[0125] One particular useful shape permitted by a device including pre-curved rods in accordance with the various embodiments is a circular are with high curvature, which can be obtained by pushing the rods on the outside of the curve and pulling those on the inside of the curve. These tighter possible curvatures result in a greater workspace when compared to the tendon configuration. In a tendon configuration, as shown in FIG. 16A, curvature is limited on either side at the point where the strain limit of one of the rods is reached. In the pre-curved wire case, the workspace is “biased” so that the strain limit is reached near a straight configuration on one side and at a very tight curvature on the other, as shown in FIG. 16B. Rotation is used to provide such biasing and therefore enables the device to reach this maximum curvature in any direction. Another interesting feature is that it is possible to rotate the tip of the tool in place without rotating the entire device.

This type of motion is accomplished by rotating the rods and the elastic structure in an equal and opposite fashion, so that shape is maintained in space while the tip of the tool rotates with the elastic structure.

[0126] As discussed above, devices with pre-curved rods in accordance with the various embodiments can be characterized and controlled using a mechanics-based mathematical model which can be used to predict the device's shape as a function of rod translation, rotation, and pre-curvatures as well as the geometric and material properties of the device. The model parameterizes the shape of the device by assuming that the corresponding elastic structure (e.g., an elastic support member and discs) is composed of multiple helical segments (one between each pair of spacing discs), where the shape of each helical segment is defined by the amount of bending and torsion within the elastic structure. The shape of the rods in the device can then be determined by the shape of the elastic structure and the locations of the pass-through holes relative to the elastic structure. The model also allows for the rods to twist axially within the pass-through holes.

[0127] The variables defining the shape are thus: (a) the amount of elastic structure bending in two directions and torsion within a segment (three variables per segment) and (b) the amount of torsion per segment in each rod (one variable per rod per segment). The elastic spring energy stored in the device and the total lengths of all rods can both be computed using this set of variables. To find the device shape one can then calculate the values of the variables which minimize the stored energy, subject to the constraints that the rod lengths are equal to their specified values. An example of the model's calculation of the shape of the device with pre-curved rods for a certain combination of tube translations and rotations is shown in FIG. 17.

[0128] The variables defining the shape may be determined by an energy minimization problem, wherein the values of the variables are those that provide a minimum of an energy. For example, define according to the preceding discussion the following variables: u_iS, which represents the curvature of the elastic structure over the ith segment; u_ij, which represents the curvature of the jth rod over the ith segment; K_S, which represents the stiffness matrix of the elastic structure, and is typically a diagonal matrix; K_j, which represents the stiffness of the jth rod; and the preformed shapes u_iS* and u_ij* which represent the preformed curvatures of the elastic structure and the rods, respectively. Then, the energy function may be written as

$E_{\mathrm{total}}=\sum _{i=1}^ML_i\left({\left(u_{\mathrm{iS}}-u_{\mathrm{iS}}^{*}\right)}^TK_S\left(u_{\mathrm{iS}}-u_{\mathrm{iS}}^{*}\right)+\sum _{j=1}^N{\left(u_{\mathrm{ij}}-u_{\mathrm{ij}}^{*}\right)}^TK_S\left(u_{\mathrm{ij}}-u_{\mathrm{ij}}^{*}\right)\right)$

Note that in this set of variables, the torsions are included implicitly through one of the components of the curvature variables. These variables are constrained by a number r functions, the form of which depends on the exact geometry of the device, which may be written implicitly as

g_k(u_1S, . . . ,u_MS,u₁₁, . . . ,u_MN)=0

for k ranging from 1 to r. Then, the values of the variables which solve the energy minimization problem may be found with a standard technique such as the method of Lagrange multipliers.

[0129] In the previous embodiments, it is generally assumed that the cross-section of the rods is substantially circular or substantially a regular polygon. However, the various embodiments are not limited in this regard. In some embodiments, the rods can have an elongated or oblong cross-section so as to form a strip, i.e., a cross-section other than a circle or a regular polygon. In such embodiments, the cross-section is substantially longer in a first direction than in a second direction perpendicular to the first direction. For example, in certain embodiments the strips can be formed by using a rectangular profile.

[0130] Such a configuration offers a number of potential advantages and can be used with straight or pre-curved rods. First, strip-type rods allow rods to be thin in one cross sectional dimension. Such a configuration allows strip-type rods to achieve higher curvatures than wire-shaped rods before reaching their strain limit. Second, by being thicker in the other cross sectional direction, strip-type rods can resist loads in the other direction extremely well. This becomes important in the context of continuum robots and endoscopes, which are typically not good at resisting loads normal to their bending plane. Thirds, strip-type rods allow for scalability. In particular, such rods offer a greater surface area for fastening tips to an end disc. Such a configuration is advantageous since the first point of failure in the devices of the various embodiments is likely to be the attachment point between an end of a rod and the disc it is attached to.

[0131] An exemplary device 1800 illustrating the use of strip-type rods is shown in FIG. 18. As shown in FIG. 18, the device 1800 is a continuously flexible mechanism, consisting of two elastic strip-type rods 1802, 1804. In the exemplary device 1800 of FIG. 8, these rods are made from nitinol. However, in other embodiments, other elastic materials could be used as well. In device 1800, the rods 1802, 1804 are constrained to remain a fixed distance from one another. As shown in FIG. 18, this is accomplished by passing the rods 1802, 1804, through a series of small spacer discs 1806, attached to an elastic support member 1808 so as to maintain even spacing between the support discs 1806. In the exemplary device 1800, the discs 1806 are 10 mm in diameter. The rods 1802, 1804, are formed using strips that are 0.33 mm thick and 3.3 mm wide, and are pre-curved into circular arcs. The device is shown in its undeflected, “home” position.

[0132] Note that other mechanisms could be used to maintain the spacing of the support discs. That is, mechanisms other than the elastic support member. For example, the strip-type rods 1802, 1804 themselves can not only serve to actuate bending in the device 1800, but can also be used to provide structural support within the device 1800. Bending actuation is achieved by pushing one of rods 1802, 1804 in by a fixed distance while pulling the other one of rods 1802, 1804 out by an equal distance. In such embodiments, the rods 1802, 1804 would not be configured to allow them to be individually axially rotated while leaving the rest of the device 1800 in place. However, axial rotation of the device 1800 can be achieved by axially rotating the device 1800 at its base.

[0133] In the various embodiments, the strips for rods may be straight or may be pre-curved. Thus, device 1800 can achieve extremely high curvatures by biasing the workspace. Note that the strip-type rods can be pre-curved into circular or non-circular profiles.

[0134] The use of thin flexible strips (i.e. elements with a rectangular cross-section in which one dimension is significantly greater than the other) offers several advantages over other flexible elements, which are outlined below.

[0135] First, the use of strips over other types of wires or other flexible shapes significantly increases the achievable curvature of the mechanism. This is because, in long beams, the achievable radius of curvature is a function of the thickness of the cross-section of the short dimension. The curvature is limited by a maximum allowable amount of strain in the material (typically the maximum strain before permanent deformation occurs, or some fraction thereof if a factor of safety is incorporated). The maximum strain occurring in a bent beam of a symmetrical cross-section is found by:

$\begin{array}{cc}\max \mathrm{bending}\mathrm{strain}=\frac{\mathrm{thickness}}{2\left(\mathrm{radius}\mathrm{of}\mathrm{curvature}\right)}& \left(29\right)\end{array}$

Based on this, the achievable radius of curvature for a given maximum strain is:

$\begin{array}{cc}\mathrm{achievable}\mathrm{radius}\mathrm{of}\mathrm{curvature}=\frac{\mathrm{thickness}}{2\left(\mathrm{allowable}\mathrm{strain}\right)}& \left(30\right)\end{array}$

Thus, as the thickness of the beam decreases, a smaller radius of curvature becomes achievable.

[0136] For a rod (or other non-strip-type member) of the same bending stiffness in the direction of “preferred” bending of the strip (see the following section), the diameter of the rod is approximately 2 times the thickness of the strip (assuming a width-to-thickness ratio of 10:1). Thus, for the equivalent rod, the smallest achievable radius of curvature would be double that of the strip (in other words, the strip would be able to reach twice the curvature the rod would be able to reach, where curvature is the inverse of radius of curvature). For example, FIG. 19 shows a photograph of straight strip 1900 of superelastic nitinol (0.33 mm thick and 3.3 mm wide), easily bending to a tight curvature. The high recoverable strain of superelastic nitinol makes it a useful material in the design of continuously flexible mechanisms, yet in principle any material that can achieve sufficient recoverable strain to satisfy application requirements may be used.

[0137] Another advantage of strip-type rods is their ability to achieve different bending stiffnesses in different directions. Using the same beam bending theory employed in the previous section on achievable curvatures, one can observe that a flat strip cross-section rod offers a substantially larger transverse stiffness (i.e. resistance to deflection by loads applied normal to the thin edge of the strip) than a circular cross-section with a bending stiffness equivalent to that of the strip in its direction of “preferred” bending. The stiffness of the rod is related to the modulus of elasticity and the moment of area by

K=EI  (31)

where K is the bending stiffness, E is the modulus of elasticity and I is the moment of area. For a rectangular cross-section of width w and thickness t, the moment of area for bending in each direction is:

I_preferred= 1/12t³w  (32)

I_transverse= 1/12t³w  (33)

Now take for example a strip-type member, such as a rod described above, for which the width is ten times the thickness (w=10t), as in the case of the rods shown in FIG. 18. The moment of area expressions can then be written as:

I_preferred=0.833t⁴  (34)

I_transverse=83.3t⁴  (33)

Thus, a strip-type member with a 10:1 width-to-thickness ratio is 100 times stiffer in the transverse direction than it is in the direction of preferred bending. In contrast, for a circular cross-section rod there is no way to achieve high transverse bending stiffness relative to the bending stiffness in the preferred bending direction, since bending stiffness will necessarily be symmetric for a rod with circular cross section.
Ability to increase torsional stiffness of the device

[0138] Use of strip-type rods also provides a way to increase the overall device's torsional stiffness, which is another noteworthy benefit. The wide cross-section of each strip provides a moment arm to resist torque as the strip-type rods press against the edges of the slits in the spacer discs. This means that in comparison with a wire-type rod mechanism, the strip-type design will be much more resistant to twisting under external loading. This is especially important as the lack of torsional stiffness is one of the major challenges associated with existing backbone-based designs.

[0139] Yet another advantage of strip-type rods is improved miniaturizability. The use of strip-type rods results in increased surface area of the rods compared to wire-type rods. Often in manufacturing, it is necessary to fabricate continuum robots and similar mechanisms at an increasingly small scale. However, a key limiting factor is simply the amount of surface area for attachment of rods to other structures. For example, the amount of surface area glue can adhere in order to attach a rod to a spacing disc. A circular cross section, by definition, has a minimal perimeter for a given cross sectional area. In contrast, so a strip-type rod will always have more surface area to which glue can be applied, and potentially much more if the rod's cross section has a large aspect ratio.

[0140] Still another advantage is simplicity in actuation and control. Control of of devices including strip-type rods will not require a complex robotic actuation system or a complex method of control, as in many continuum mechanisms. Rather, a simple geometrically derived model effectively describes the actuation of the mechanism. This model is based on the property that the mechanism remains in constant curvature, i.e. that it remains in a circular arc as it is actuated.

[0141] In operation, the pose of the mechanism can be described by two variables: bending of the mechanism (which can be equivalently described either by a radius of curvature, r, or the corresponding angle subtended by the arc of the mechanism, θ), and the axial rotation of the robot (described by an angle measured from some nominal or “home” position of the mechanism, φ). These variables constitute the configuration space of the mechanism. The variables used to define how the mechanism is actuated are the displacements of the two strip-type rods (δ_o for the outer rod and δ_i for the inner rod), and the rotation of the mechanism's base (the same φ previously defined). These variables constitute the actuation space of the mechanism.

[0142] The following relationship between the length of the centerline of the mechanism L and the radius of curvature r and 9 holds true as the mechanism moves: L=rθ. The length L remains constant during actuation. The lengths of the inner and outer rods are found using the same relationship to be:

L_o=(r+D)θ=L+Dθ  (36)

L_i=(r−D)θ=L−Dθ  (37)

where D is the distance from each strip to the centerline, i.e. half the distance between the strips. When the mechanism is in its nominal configuration, the lengths are functions of the nominal radius of curvature and angle, such that

L_o,nom=(r_nom+D)θ_nom=L+Dθ_nom  (38)

L_o,nom=(r_nom−D)θ_nom=L−Dθ_nom  (39)

The variables of the extension and retraction of the rods (δ_i and δ_o) are found to be the difference between each length for a given curvature and the length of the strip in its nominal curvature, i.e.

δ_o=ΔL_o=L_o−L_o,nom=(θ−θ_nom)D=+DΔθ  (40)

δ_i=ΔL_i=L_i−L_i,nom=−(θ−_nom)D=−DΔθ  (41)

implying that the required changes in length of the inner and outer rods are simply equal but opposite. Thus, a single variable δ can be used to describe both δ_i and δ_o by:

δ=δ_o=−δ_i  (42)

This one-to-one coupling of the extension and retraction of the two rods implies that the relationship can easily be enforced in software for robotically controlled systems, or could be enforced through a mechanical coupling of the two translations, making it possible to use the mechanism in a manual tool.

[0143] The rotational degree-of-freedom is fully described by the angle φ, which can be used interchangeably between the actuation space and configuration space variable. Coordinate frames used in modeling kinematics of the device can be defined as in FIG. 20. Given the coordinate frames defined in FIG. 20, the homogeneous transformation describing the tip frame with respect to the base frame for bending only (φ=0) is:

$\begin{array}{cc}{T}_{\theta }=\left[\begin{array}{cccc}\cos \theta & 0& \sin \theta & \frac{L}{\theta }\left(1-\cos \theta \right)\\ 0& 1& 0& 0\\ -\sin \theta & 0& \cos \theta & \frac{L}{\theta }\sin \theta \\ 0& 0& 0& 1\end{array}\right]& \left(43\right)\end{array}$

Similarly, the following transformation describes the motion for rotation only (θ=θ_nom):

$\begin{array}{cc}{T}_{\phi }=\left[\begin{array}{cccc}\cos \phi & -\sin \phi & 0& 0\\ \sin \phi & \cos \phi & 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]& \left(44\right)\end{array}$

The overall transformation is a combination of these two transformations:

T_T^B=T_φT_θ  (45)

Multiplying this expression out fully and substituting in the relationships between the configuration variables and the actuation variables gives the overall transformation in terms of the actuation variables as follows:

$\begin{array}{cc}{T}_T^B=\left[\begin{array}{cccc}\cos \mathrm{\phi cos}\left({\theta }_{\mathrm{nom}}+\frac{\delta }{D}\right)& -\sin \phi & \cos \mathrm{\phi sin}\left({\theta }_{\mathrm{nom}}-\frac{\delta }{D}\right)& \frac{-\mathrm{DL}\cos \phi }{\delta +D{\theta }_{\mathrm{nom}}}\left(\cos \left(\frac{\delta }{D}\right)-1\right)\\ \sin \mathrm{\phi cos}\left({\theta }_{\mathrm{nom}}+\frac{\delta }{D}\right)& \cos \phi & \sin \mathrm{\phi sin}\left({\theta }_{\mathrm{nom}}+\frac{\delta }{D}\right)& \frac{2\mathrm{DL}\sin \phi }{\delta +D{\theta }_{\mathrm{nom}}}\left({\sin }^2\left(\frac{\delta +D{\theta }_{\mathrm{nom}}}{2D}\right)\right)\\ -\sin \left({\theta }_{\mathrm{nom}}+\frac{\delta }{D}\right)& 0& \cos \left({\theta }_{\mathrm{nom}}+\frac{\delta }{D}\right)& \frac{\mathrm{DL}}{\delta +D{\theta }_{\mathrm{nom}}}\sin \left({\theta }_{\mathrm{nom}}+\frac{\delta }{D}\right)\\ 0& 0& 0& 1\end{array}\right]& \left(46\right)\end{array}$

This transformation can be used to easily compute the hybrid Jacobian for differential inverse kinematics.

[0144] The preceding discussion assumes the presence of elastic structure to support a platform or disc at the distal end of a continuum robot or other device. However, the various embodiments are not limited in this regard. In particular, multiple rods, attached in parallel to a platform or end disc, can provide all the support needed by such a platform. Thus, a parallel continuum manipulator can be defined. Such parallel continuum manipulators have the potential to inherit some of the compactness and compliance of continuum robots while retaining some of the precision, stability, and strength of rigid-link parallel robots, yet they represent a relatively unexplored area of the broad manipulator design space.

[0145] A. Concept and Motivation

[0146] The basic design paradigm of a parallel continuum robot is illustrated in FIG. 21. In particular, FIG. 21 schematically illustrates the parallel continuum design paradigm utilizes the translation of flexible rods or legs between constrained points to create large deflections in each of the flexible legs. In some embodiments, some or all of the legs can be strip-type rods, as discussed above. The distal ends of the multiple compliant legs are connected in some pattern to an end piece to effect a platform, and the base of each leg is independently translated by an actuator at its proximal end. There may also be some additional constraints on the relative positions of the rods at a certain point (e.g. linear bearings where the rods pass through a base platform). One can conceive of several different variations of parallel continuum robot designs. As with rigid-link parallel manipulators, kinematic behavior is defined by the number of legs and the types of actuation and constraints at each end, which could vary from spherical to clamped joints (our prototype manipulator in FIG. 3 uses clamped joints). The elastic legs themselves could be initially straight or have generally pre-curved stress free states, which could alter the kinematic properties of the device. Using tubes for the legs would allow wires or other infrastructure to be passed through their hollow center channels to the distal platform in order to actuate a gripper or other end effector. Using tubes would also enable hybrid “stacked” designs to be constructed (similar to multi-segment, snake-like designs, except with more degrees of freedom per segment) where the legs of a secondary parallel continuum mechanism extend out the ends the tubes through the distal platform of the primary mechanism.

[0147] In comparison to serial continuum robots, parallel continuum designs are likely to have higher payload capacity and accuracy/repeatability. Compared to rigid-link parallel manipulators, they may exhibit greater compliance, larger workspace, and easier miniaturization to the scale of a few millimeters in diameter or smaller. Their inherent mechanical compliance and low mass due to off-loading of the actuators can provide an important safety feature wherever parallel robots need to interact with humans. In the field of endoscopic robotic surgery, parallel continuum manipulators have the potential to provide precise, multi-DOF motion in a simple, compact, and short mechanism at the tip of an endoscope. This capability may be useful for manipulating objects within highly confined spaces where a long slender body would be constrained by anatomical structures. In these cases, access to the site could be gained through flexible endoscopy, and then a parallel continuum device could provide fine, multi-DOF motion at the tip location (actuated by rods passing through hollow channels).

[0148] B. Prototype Design And Construction

[0149] In this section described is the construction of a parallel continuum manipulator prototype structure designed to show proof-of-concept and study parallel continuum robot behavior. This is shown in FIG. 22. In FIG. 22, the middle picture shows our prototype parallel continuum manipulator constructed for initial proof of concept and evaluation of our modeling approach. Close-up views of the base and end platforms are shown to the left in FIG. 22. The dimensions and patterns of the base plate pass-through holes and the top platform connection points are detailed to the right in FIG. 22. The lengths of the compliant legs between the two platforms are actuated by manually translating linear slides connected to the rod ends below the base platform.

[0150] The prototype has six compliant legs connected in a similar arrangement to the legs of a 6-DOF rigid-link Stewart-Gough platform, and each leg can be manually actuated. The legs consist of 1.3 mm diameter spring steel music wire (ASTM A228) with an estimated Young's modulus of 207 GPa and Poisson's ratio of 0.305. As shown in FIG. 22, these rods are connected to an end-effector plate of clear 1=16 inch acrylic via 3=8 inch OD shaft collars constrained in channels within small blocks attached to the plate. The rods are routed through holes in the base plate, and the proximal ends are connected to linear slide carriages in the same fashion. The carriages translate along T-slotted aluminum rails (80/20® Inc.) that are bolted in a hexagonal pattern to the base plate. The linear slide carriages can be manually repositioned and locked in place with a brake so that the length of each rod between the base plate and the end plate can be actuated independently.

[0151] The connection locations of the 6 flexible legs are arranged in a conventional radial hexapod pattern of 3 pairs of rods spaced 120° apart at a radius of 87 mm. As depicted by the numbers in FIG. 22, the proximal holes for rods 1 and 2 are paired together (with a total separation of 20°), while at the end plate the connections for rods 2 and 3 are paired together. The same pattern follows for the other pairs of wires. Rods 3 and 4, and 5 and 6 are paired at the base plate, while 4 and 5, and 6 and 1 are paired at the distal end. In the neutral configuration shown in FIG. 4 with all leg lengths equal, this connection pattern causes the top plate to be rotated by 60° with respect to the hole pattern in the base plate and causes all of the rods to bend from their naturally straight state.

[0152] C. Modeling Approach

[0153] Adopted is a kinematic modeling approach for parallel continuum robots that is based on classical Cosserat rod mechanics. The framework outlined in the following subsections addresses the forward and inverse kinematics and statics problems.

[0154] i. Cosserat Rod Equations

[0155] The shape of each component rod in the robot is defined by its position p_i(s_i)ε and material orientation R_i(s_i)εSO(3)(, forming a material-attached reference frame

$g_i\left(s_i\right)=\left[\begin{array}{cc}{R}_i\left(s_i\right)& p_i\left(s_i\right)\\ 0& 1\end{array}\right]\in \mathrm{SE}\left(3\right)$

as a function of arc length s_iε measured from the proximal platform. The position and orientation evolve along the length of the rod according to kinematic variables v_i(s)ε and ui(s) 2 R3, which describe the linear and angular rates of change expressed in local or body frame coordinates of the material frame as follows:

p_i′=R_iv_i,

R_i′=R_iû_i,  (47)

where ′ denotes a derivative with respect to s_i, and ̂ denotes mapping from r to so(3) as follows,

$\begin{array}{cc}\hat{a}=\left[\begin{array}{ccc}0& -a_3& a_2\\ a_3& 0& -a_1\\ -a_2& a_1& 0\end{array}\right].& \left(48\right)\end{array}$

{hacek over ( )} is used to denote the inverse mapping of ̂, i.e., (û){hacek over ( )}=u.

[0156] The rates of change of the internal force vector n and internal moment vector m with respect to the arc length s_i are described by the classical Cosserat rod differential equations of static equilibrium:

n_i′=−f_i

m_i′=−p_i′×n_i−l_i,  (49)

where all vectors are assumed expressed in global coordinates, and f_i and l_i are distributed force and moment vectors respectively applied per unit length to rod i. Distributed self weight and any other external forces are straightforward to include within f_i and l_i.

[0157] The kinematic variables v_i and u_i are related to material strain (shear, extension, bending, and torsion) and can be used to calculate the internal force and moment vectors (denoted by vectors n and m respectively and expressed in global coordinates) via a material constitutive law. For the present, we use a linear constitutive relationship of the following form,

$\begin{array}{cc}{n}_i=R_iK_{\mathrm{se},i}\left(v_i-v_i^{*}\right),K_{\mathrm{se},i}=\left[\begin{array}{ccc}{A}_iG_i& 0& 0\\ 0& A_iG_i& 0\\ 0& 0& A_iE_i\end{array}\right]m_i=R_iK_{\mathrm{bt},i}\left(u_i-u_i^{*}\right),K_{\mathrm{bt},i}=\left[\begin{array}{ccc}{E}_iI_i& 0& 0\\ 0& E_iI_i& 0\\ 0& 0& J_iG_i\end{array}\right]& \left(50\right)\end{array}$

[0158] Where v_i* and u_i* are the kinematic variables of the rod in an assigned stress free reference state. For an initially straight rod, appropriate reference state variables are v_i*=[0 0 1]^T and u_i*=[0 0 0]^T. The matrices K_se,i and K_bt,i here contain the stiffness terms for a radially symmetric rod cross-section which could vary with arc length, involving the area A_i, Young's modulus, E_i, the shear modulus G_i, the second area moment I_i (about the local x and y axes), and the polar area moment J_i about the local z axis.

[0159] Thus, for each rod, the equations discussed above can be used to form a system of differential equations that describes the evolution of the state variables p_i, R_i, m_i, and n_i with respect to S_i.

[0160] ii. Boundary Conditions for Forward Kinematics

[0161] Each rod or tube in a parallel continuum robot is independently described by the system of differential equations above. However, the boundary conditions of each system are coupled because of the physical constraints inherent to the robot structure.

[0162] For the prototype robot design which describe herein, the proximal end of each rod is clamped with a set-screw shaft collar constrained within a groove in acrylic block, and each rod subsequently passes through a cylindrical hole in the base platform as shown in FIG. 3. This design constrains the position pi and the rod tangent vector at the base platform while allowing rotational freedom about the tangent axis so that no torsional moment can be supported. Thus, the torsional moment is set to zero at the base of each rod and R_i(0) is expressed as a rotation about the global z-axis by some angle θ_i as follows:

$\begin{array}{cc}{m}_{\mathrm{iz}}\left(0\right)=0R_{i}\left(0\right)=\left[\begin{array}{ccc}\cos {\theta }_i& -\sin {\theta }_i& 0\\ \sin {\theta }_i& \cos {\theta }_i& 0\\ 0& 0& 1\end{array}\right]& \left(51\right)\end{array}$

[0163] At the distal end of each rod (s_i=L_i), the following conditions of static equilibrium must hold for the top plate:

$\begin{array}{cc}\sum _{i=1}^n\left[n_{i}\left(L_i\right)\right]-F=0\sum _{i=1}^n\left[p_i\left(L_i\right)\times n_i\left(L_i\right)+m_i\left(L_i\right)\right]-p_c\times F-M=0,& \left(52\right)\end{array}$

where F and M are external force and moment vectors applied at the centroid of the top plate, p_c. The clamped connections of the distal end of the rods to the top plate are modeled by setting the material orientations of each rod at L_i equal to each other:

[log(R_i^T(L_i)R₁(L₁))]{hacek over ( )}=0 for i=2 . . . n,  (53)

where log ( ) is the matrix natural logarithm, which maps SO(3) to so(3), and the {hacek over ( )} operator subsequently maps so(3) to . This forces a common material orientation for all the distal rod ends. Finally, one can write the following equations for the rod end positions.

p₁(L₁)−p_i(L_i)−R₁(L₁)(r₁−r_i)=0 for i=2 . . . n,  (54)

where r_i is the connection position for rod i expressed with respect to the top plate coordinate system. These are analogous to loop closure equations because they are only satisfied when the positions of the rod ends have the same relative positions as the connection pattern in the top plate.

[0164] iii. Boundary Conditions for Inverse Kinematics

[0165] For the inverse kinematics problem, the boundary conditions given above in (5) and (6) still apply, but the geometric coupling of the rods at the distal platform is simplified by the fact that a desired position and orientation is known. This can be expressed by a constraint on the position and orientation of each distal leg end as follows:

[log(R_i^T(L_i)R_d)]{hacek over ( )}=0 for i=1 . . . n,

p_d+R_dτ_i−p_i=0 for i=1 . . . n.  (55)

where p_d and R_d are the desired position (of the centroid) and orientation of the distal platform. We note that these equations are applied for i=1 . . . n, so there are a total of six more scalar constraint equations to be solved than in the forward kinematics case, and the leg lengths constitute n additional unknowns.

[0166] iv. Numerical Computation

[0167] A simple shooting method provides an effective way to solve the systems of rod equations subject to these coupled boundary conditions. In this procedure, the unknown boundary conditions at the base of each rod (n_i(0), m_ix(0), m_iy(0), and θ_i) are guessed (in the inverse kinematics case, the leg lengths Li are included in this set of unknowns). Each system of rod equations is numerically integrated from s_i=0 to L_i as an initial value problem using a standard numerical routine such as a Runge-Kutta method, and the boundary condition equations are subsequently evaluated. This process is nested within an optimization loop which iteratively updates the guessed values for the unknown proximal boundary conditions until the distal boundary conditions are satisfied within an appropriate tolerance. This is an efficient and relatively stable method of computation. For the examples shown in the experimental section, convergence was achieved in every case, starting with an uninformed initial guess of zero for all unknown proximal conditions. Computation time was on the order of 10 seconds using unoptimized MATLAB code executed on a typical laptop computer. Computation time decreases to about 1 second if a good initial guess is available, which enables simulation of quasistatic motion since a previous solution can be used as the initial guess for the next iteration.

[0168] v. Multiple Solutions and Buckling Behavior

[0169] In contrast to rigid-link robots, the forward and inverse kinematics problems for parallel continuum manipulators are both complex (with the potential for multiple solutions possible) due to the large-deflection elastic behavior involved. For the numerical computation strategy above, the set of values for the unknown initial conditions (and leg lengths in the case of inverse kinematics) which satisfies the distal end boundary conditions may not always be unique. In this context, the phenomena of “buckling” describes the sudden transition of one valid solution to another one, which may occur when a particular equilibrium solution becomes unstable (or ceases to exist) as the result of a small change in actuation or external loading. The model equations are still valid in these situations, and they contain the ability to accurately describe all possible “buckled” states if the correct solution for the initial conditions is selected. The quasistatic simulation process described above tends to produce the correct solution until a point of instability is reached. The potential for buckling and instability exists for all continuum robots when external loads are considered, and for some designs even without external loading (e.g. concentric-tube robots).

[0170] D. Kinematic Simulation and Analysis

[0171] FIG. 23 demonstrates the kinematic degrees of freedom and range of motion of the prototype manipulator structure described herein by showing the manipulator shape in five different configurations. In particular, FIG. 23 shows photographs of five experimental configurations, demonstrating the ability of our prototype parallel continuum manipulator to execute axial twist, transverse tilt, translation, and bending (simultaneous translation and tilt). To the right of each experimental photograph, a MATLAB rendering of the rod-mechanics-based kinematic model prediction for the same actuator configurations is shown.

[0172] With all leg lengths equal, the robot is in a straight, neutral state, which can be raised or lowered by equal translation of the legs. By extending legs 2, 4, and 6, the distal platform twists about the z-axis. Translating legs 4 and 5 causes the platform to tilt, and translating legs 3 and 6 causes translation. A combination of rotation and translation (bending) is achieved by extending legs 4 and 5 while retracting legs 1 and 2. The collection of legs maintains a compact form throughout all these motions, which is encouraging for potential applications in confined spaces. FIG. 23 also shows a rendering of our forward kinematics model solution for these five cases, indicating the feasibility and qualitative accuracy of the modeling approach over a wide range of motion.

[0173] To further illustrate the kinematics of the prototype design, the body-frame manipulator is given a Jacobian, J_b, which was computed numerically for the prototype manipulator in the neutral configuration (where all leg lengths are 406 mm, and the rotation is aligned with the global reference frame.)

$J^b=\left[\begin{array}{cccccc}-1.62& -1.62& 1.83& -0.21& -0.21& 1.83\\ -1.18& 1.18& -0.82& -2.00& 2.00& 0.82\\ 0.17& 0.17& 0.17& 0.17& 0.17& 0.17\\ -0.12& 0.12& 0.24& 0.12& -0.12& -0.24\\ -0.20& -0.20& 0.00& 0.20& 0.20& 0.00\\ -0.65& 0.65& -0.65& 0.65& -0.65& 0.65\end{array}\right]$

where the top 3 rows are dimensionless, and the bottom three rows have units of degrees/mm. The matrix is full rank and well-conditioned, indicating that in the neutral configuration, actuators can easily move the top platform in any direction in the 6 DOF space of rigid body motion. In FIG. 24, depicted is the inverse kinematic mapping over 3 two-dimensional slices of the workspace. The figure shows the required leg length and axial tension at the base of each leg as a function of the desired end-effector pose, which is specified by position p_d=[x y z]^T and orientation R_d=exp({circumflex over (θ)}), θ=[θ_x θ_y θ_z]^T. For each slice, two of the six pose variables were varied over a 9×9 grid of values for which the inverse kinematics computation easily converged, while all other pose variables were held constant at their nominal values of x=y=θ_x=θ_y=θ_Z=0, and z=400 mm.

[0174] The first case (FIG. 24: (a) and (d)) shows an approximately linear kinematic mapping for desired rotation about the x and y axes. The second case (FIG. 24: (b) and (e)) shows an approximately quadratic leg length mapping and an approximately linear tension mapping for desired translation in the x-y plane. Note that in this case, the required lengths and tensions are identical for legs 1 and 4, 2 and 5, and 3 and 6 respectively. The final case (FIG. 24: (c) and (f)) shows an approximately linear leg length mapping and an approximately cubic tension mapping for desired motion which both rotates about and translates along the global z axis. In this case, the lengths and tensions are identical for legs 1, 3, and 5, and 2, 4, and 6 respectively.

[0175] The axial tension at the base of each leg corresponds to the actuator force which would be required to hold the robot in a particular configuration. Over these ranges of motion, moderate forces on the order of 5 N will be required to actuate the structure. A simulation of this kind can be used in the design process to size motors for a particular set of manipulator structural parameters, and to limit the length that the legs extend below the base platform to avoid buckling. One can also compute the output stiffness matrix at the top platform, which maps displacement of the end effector centroid to applied force dF=Kdp_c. For the neutral configuration in FIG. 4, the computed stiffness matrix is

$K=\left[\begin{array}{ccc}17& 0& 0\\ 0& 17& 0\\ 0& 0& 122\end{array}\right]N\text{/}\mathrm{mm},$

while for the bending configuration in FIG. 23, the stiffness matrix is computed to be

$K=\left[\begin{array}{ccc}12.7& 0& 12.8\\ 0& 0.5& 0\\ 12.8& 0& 14.5\end{array}\right]N\text{/}\mathrm{mm},$

[0176] These computations show that the stiffness in the y direction (out of the page for the bending case in FIG. 23), is only about 3% of its value in the straight case, which is consistent with this prototype. Care should be taken when designing and controlling a parallel continuum manipulator to ensure that output stiffness is sufficient for the desired tasks.

[0177] A set of model validation experiments was performed by photographing the shape of the prototype manipulator in front of a graph poster in 14 different actuator configurations, which are listed in Table I. FIG. 25 shows the basic elements of each image that were measured, the top platform centroid and orientation in the global x-z plane. The three main measurements taken from the images are the position of the centroid of the end effector plate in the x and z, which is denoted by the marker at the top, and the orientation of the end effector in the x-z plane. The orientation was measured by drawing a line along the end effector edge and then measuring the slope of the line with the grid poster.

[0178] These measurements were then compared to our forward kinematics model prediction for each actuator configuration. The photos of these planar cases were taken perpendicular to the graph plane approximately 30 feet away from the robot so that perspective error was minimized.

[0179] The resulting differences between the data and the model prediction are presented in Table I. The positional error was calculated as the total Cartesian error in the global x-z plane, and the percent error was calculated as the position error divided by the average leg length for each case. The maximum positional error was 11.74 mm with a 2.89% associated percent error. For the configurations resulting in a change in top platform orientation (7-14 in Table I), angular displacement was measured graphically as shown in FIG. 6, and the maximum angular difference between model prediction and experiment was 3.14 degrees.

TABLE I
EXPERIMENTAL CONFIGURATIONS AND
MODEL ERROR (MM)
#	L1	L2	L3	L4	L5	L6	Error	% Error
1	406	406	406	406	406	406	3.4	0.8
2	386	406	386	406	386	406	6.1	1.6
3	426	406	426	406	426	406	2.9	0.7
4	406	406	426	406	406	426	8.8	2.1
5	406	406	386	406	406	386	7.1	1.8
6	406	406	366	406	406	366	5.4	1.4
7	406	406	406	386	386	406	7.5	1.8
8	406	406	406	366	366	406	10.4	2.7
9	406	406	406	426	426	406	8.9	2.1
10	406	406	406	446	446	406	7.6	1.8
11	426	426	406	386	386	406	7.0	1.7
12	446	446	406	366	366	406	10.4	2.6
13	386	386	406	426	426	406	9.3	2.3
14	366	366	406	446	446	406	11.7	2.9

[0180] Observed was a small amount of flex in the top acrylic plate during the experiments. This unmodeled effect is a source of error between model prediction and actual manipulator shape.

[0181] A. Buckling Experiments

[0182] Also experimentally determined was a critical buckling load for the prototype structure in three different actuator configurations, (1, 4, and 5 in Table I, corresponding to a straight neutral case and two translating cases, respectively). Applied was an incrementally increasing force in the negative z direction (down) through the centroid of the top plate until at least one leg in the structure began to buckle. The resulting experimental buckling loads are shown by the black dashed lines in FIGS. 26A-C. FIGS. 26A-C show displacements of the top platform centroid in the x, y, and z directions are plotted by using our forward kinematics model to perform quasistatic simulation with an incrementally increasing vertical load for three different manipulator configurations (1, 4, and 5 in Table I). The dashed lines depict the experimentally determined critical buckling loads for the same cases. In FIG. 26A, the simulation converged to a buckled solution at an applied load very close to the experimentally determined buckling load. In FIGS. 26B and 26C the simulation failed to converge to any solution just prior to when the simulated load would have reached the experimentally determined buckling load.

[0183] When running a forward kinematics model simulation with the same incrementally increasing force, the simulation converged to a buckled solution (shown by the sudden large centroid displacement in FIG. 26A at a load close to the experimentally determined buckling load. For configurations 4 and 5 (FIGS. 26B and 26C), the simulation failed to converge to a valid forward kinematics solution just prior to reaching the experimentally determined buckling load, indicating that the buckled mode is not relatively close to the unbuckled state.

[0184] Turning now to FIG. 27, this is shown an exemplary system 2700 includes a general-purpose computing device 2700 for performing one or more of the various methods and processes described above. System 2700 includes a processing unit (CPU or processor) 2720 and a system bus 2710 that couples various system components including the system memory 2730 such as read only memory (ROM) 2740 and random access memory (RAM) 2750 to the processor 2720. The system 2700 can include a cache 2722 of high speed memory connected directly with, in close proximity to, or integrated as part of the processor 2720. The system 2700 copies data from the memory 2730 and/or the storage device 2760 to the cache 2722 for quick access by the processor 2720. In this way, the cache 2722 provides a performance boost that avoids processor 2720 delays while waiting for data. These and other modules can be configured to control the processor 2720 to perform various actions. Other system memory 2730 may be available for use as well. The memory 2730 can include multiple different types of memory with different performance characteristics. It can be appreciated that the disclosure may operate on a computing device 2700 with more than one processor 2720 or on a group or cluster of computing devices networked together to provide greater processing capability. The processor 2720 can include any general purpose processor and a hardware module or software module, such as module 1 2762, module 2 2764, and module 3 2766 stored in storage device 2760, configured to control the processor 2720 as well as a special-purpose processor where software instructions are incorporated into the actual processor design. The processor 2720 may essentially be a completely self-contained computing system, containing multiple cores or processors, a bus, memory controller, cache, etc. A multi-core processor may be symmetric or asymmetric.

[0185] The system bus 2710 may be any of several types of bus structures including a memory bus or memory controller, a peripheral bus, and a local bus using any of a variety of bus architectures. A basic input/output (BIOS) stored in ROM 2740 or the like, may provide the basic routine that helps to transfer information between elements within the computing device 2700, such as during start-up. The computing device 2700 further includes storage devices 2760 such as a hard disc drive, a magnetic disc drive, an optical disc drive, tape drive or the like. The storage device 2760 can include software modules 2762, 2764, 2766 for controlling the processor 2720. Other hardware or software modules are contemplated. The storage device 2760 is connected to the system bus 2710 by a drive interface. The drives and the associated computer readable storage media provide nonvolatile storage of computer readable instructions, data structures, program modules and other data for the computing device 2700. In one aspect, a hardware module that performs a particular function includes the software component stored in a tangible and/or intangible computer-readable medium in connection with the necessary hardware components, such as the processor 2720, bus 2710, display 2770, and so forth, to carry out the function. The basic components are known to those of skill in the art and appropriate variations are contemplated depending on the type of device, such as whether the device 2700 is a small, handheld computing device, a desktop computer, or a computer server.

[0186] Although the exemplary embodiment described herein employs the hard disc 2760, it should be appreciated by those skilled in the art that other types of computer readable media which can store data that are accessible by a computer, such as magnetic cassettes, flash memory cards, digital versatile discs, cartridges, random access memories (RAMs) 2750, read only memory (ROM) 2740, a cable or wireless signal containing a bit stream and the like, may also be used in the exemplary operating environment. Tangible, non-transitory computer-readable storage media expressly exclude media such as energy, carrier signals, electromagnetic waves, and signals per se.

[0187] To enable user interaction with the computing device 2700, an input device 2790 represents any number of input mechanisms, such as a microphone for speech, a touch-sensitive screen for gesture or graphical input, keyboard, mouse, motion input, speech and so forth. An output device 2770 can also be one or more of a number of output mechanisms known to those of skill in the art. In some instances, multimodal systems enable a user to provide multiple types of input to communicate with the computing device 2700. The communications interface 2780 generally governs and manages the user input and system output. There is no restriction on operating on any particular hardware arrangement and therefore the basic features here may easily be substituted for improved hardware or firmware arrangements as they are developed.

[0188] For clarity of explanation, the illustrative system embodiment is presented as including individual functional blocks including functional blocks labeled as a “processor” or processor 2720. The functions these blocks represent may be provided through the use of either shared or dedicated hardware, including, but not limited to, hardware capable of executing software and hardware, such as a processor 2720, that is purpose-built to operate as an equivalent to software executing on a general purpose processor. For example the functions of one or more processors presented in FIG. 27 may be provided by a single shared processor or multiple processors. (Use of the term “processor” should not be construed to refer exclusively to hardware capable of executing software.) Illustrative embodiments may include microprocessor and/or digital signal processor (DSP) hardware, read-only memory (ROM) 2740 for storing software performing the operations discussed below, and random access memory (RAM) 2750 for storing results. Very large scale integration (VLSI) hardware embodiments, as well as custom VLSI circuitry in combination with a general purpose DSP circuit, may also be provided.

[0189] The logical operations of the various embodiments are implemented as: (1) a sequence of computer implemented steps, operations, or procedures running on a programmable circuit within a general use computer, (2) a sequence of computer implemented steps, operations, or procedures running on a specific-use programmable circuit; and/or (3) interconnected machine modules or program engines within the programmable circuits. The system 2700 shown in FIG. 27 can practice all or part of the recited methods, can be a part of the recited systems, and/or can operate according to instructions in the recited tangible computer-readable storage media. Such logical operations can be implemented as modules configured to control the processor 2720 to perform particular functions according to the programming of the module. For example, FIG. 27 illustrates three modules Mod1 2762, Mod2 2764 and Mod3 2766 which are modules configured to control the processor 2720. These modules may be stored on the storage device 2760 and loaded into RAM 2750 or memory 2730 at runtime or may be stored as would be known in the art in other computer-readable memory locations.

[0190] A. Nomenclature

• • [0191] *: Denotes a variable in the reference state.
  • [0192] ′: Denotes a derivative with respect to s.
  • [0193] ̂: Coverts to so(3) and to o(3):

$u=\left[\begin{array}{ccc}0& -u_z& u_y\\ u_z& 0& -u_x\\ -u_y& u_x& 0\end{array}\right],\widehat{\left[\begin{array}{c}v\\ u\end{array}\right]}=\left[\begin{array}{cccc}0& -u_z& u_y& u_x\\ u_z& 0& -u_x& u_y\\ -u_y& u_x& 0& u_z\\ 0& 0& 0& 0\end{array}\right]$

• • [0194] {hacek over ( )}: Inverse of the ̂operation. (û){hacek over ( )}=u.

[0195] s: ε —Reference length parameter.

[0196] p(s): ε— Position of the robot backbone centroid in global frame coordinates.

[0197] R(s): εSO(3)—Orientation of the robot backbone material with respect to the global frame.

[0198] g(s): εSE(3)—Homogeneous transformation containing R(s) and p(s). (The “body frame”)

[0199] r_i(s) ε: Position of the i^th tendon with respect to the body frame. r_i(s)=[x_i(s)y_i(s)0]^T

[0200] p_i(s): Ε—Position of the i^th tendon in global frame coordinates. p_i(s)=Rr_i(s)+p(s).

[0201] u(s): ε—Angular rate of change of g with respect to s in body-frame coordinates. u=(R_T{dot over (R)}){hacek over ( )}

[0202] v(s): ε—Linear rate of change of g with respect to s expressed in body-frame coordinates. v=R^T{dot over (p)}

[0203] n(s): ε—Internal force in the backbone expressed in global frame coordinates.

[0204] m(s): ε—Internal moment in the backbone expressed in global frame coordinates.

[0205] f_e(s): ε—External force per unit s on the backbone expressed in global frame coordinates.

[0206] l_e(s): ε—External moment per unit s on the backbone expressed in global frame coordinates.

[0207] f_t(s): ε—Sum of all forces per unit s applied to the backbone by tendons, expressed in global frame coordinates.

[0208] l_t(s): ε—Sum of all moments per unit s applied to the backbone by the tendons, expressed in global frame coordinates.

f_i(s): ε—Force per unit s applied to the i^th tendon by it surrounds.

[0209] n_i(s): ε—Tension in the i^th tendon.

[0210] τ_i: ε—Tension in the i^th tendon. It is constant along s under the frictionless assumption.

[0211] B. Derivation of f_i(s)

[0212] Beginning with (11),

$n_i={\tau }_i\frac{{\stackrel{.}{p}}_i}{{\stackrel{.}{p}}_i},$

one can re-arrange and differentiate to obtain

${\stackrel{.}{p}}_i=\frac{1}{{\tau }_i}{\stackrel{.}{p}}_in_i,{\ddot{p}}_i=\frac{1}{{\tau }_i}\left(\frac{}{s}\left({\stackrel{.}{p}}_i\right)n_i+\left({\stackrel{.}{p}}_i\right){\stackrel{.}{n}}_i\right).$

Noting that n_i×n_i=0, one can take a cross product of the two results above to find,

${\ddot{p}}_i\times {\stackrel{.}{p}}_i=\frac{{{\stackrel{.}{p}}_i}^2}{{\tau }_i^2}\left({\stackrel{.}{n}}_i\times n_i\right)$

and so

${\stackrel{.}{p}}_i\times \left({\ddot{p}}_i\times {\stackrel{.}{p}}_i\right)=\frac{{{\stackrel{.}{p}}_i}^3}{{\tau }_i^3}\left(n_i\times \left({\stackrel{.}{n}}_i\times n_i\right)\right).$

Applying the vector triple product identity, a×(b×c)=b(a·c)−c(a·b), one can expand the right-hand side of this equation. Since τ_i (the magnitude of n_i) is constant with respect to s, then n_i·{dot over (n)}_t=0, and this results in

$f_i=-{\stackrel{.}{n}}_i=-{\tau }_i\frac{{\stackrel{.}{p}}_i\times \left({\ddot{p}}_i\times {\stackrel{.}{p}}_i\right)}{{{\stackrel{.}{p}}_i}^3}.$

Using the fact that a×b=−b×a, and writing the cross products in skew-symmetric matrix notation (a×b=âb), one arrives at

$f_i={\tau }_i\frac{{\stackrel{.}{p}}_i\times \left({\stackrel{.}{p}}_i\times {\ddot{p}}_i\right)}{{{\stackrel{.}{p}}_i}^3}={\tau }_i\frac{{\widehat{\stackrel{.}{p}}}_i^2}{{{\stackrel{.}{p}}_i}^3}{\ddot{p}}_i.$

[0213] While various embodiments of the present invention have been described above, it should be understood that they have been presented by way of example only, and not limitation. Numerous changes to the disclosed embodiments can be made in accordance with the disclosure herein without departing from the spirit or scope of the invention. Thus, the breadth and scope of the present invention should not be limited by any of the above described embodiments. Rather, the scope of the invention should be defined in accordance with the following claims and their equivalents.

[0214] Although the invention has been illustrated and described with respect to one or more implementations, equivalent alterations and modifications will occur to others skilled in the art upon the reading and understanding of this specification and the annexed drawings. In addition, while a particular feature of the invention may have been disclosed with respect to only one of several implementations, such feature may be combined with one or more other features of the other implementations as may be desired and advantageous for any given or particular application.

[0215] For example, in some embodiments, the method described above can also be used to determine the load or external forces and moments being applied to the elastic support member. In such embodiments, the 3D shape of the elastic support member can be determined via some kind of sensing method (i.e., with cameras, or optical fibers, or magnetic tracking coils, or ultrasound, or fluoroscopy, etc.). Thereafter, using a known tension on the rod and routing path for the rod, the iterative model equations described above can be iteratively solved to determine the external forces and moments (f_e and l_e) which result in the model-predicted shape that is close to the actual sensed shape. The resulting loads based on the model can then be used as an estimate of the loads acting on the elastic support member. Accordingly, these loads can be used to provide useful information to one who is operating the continuum robot. Alternatively, a similar method can be used to compute the required rod tension necessary to achieve forces and moments for the continuum robot to exert on its surroundings. In such embodiments, the 3D shape of the elastic support member is also determined via some kind of sensing method. Thereafter, the external loads are estimated using the above-mentioned procedure. Finally, the adjustment in tension needed to achieve a desired load or shape can be determined iteratively using the system of equations.

[0216] The terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. As used herein, the singular forms “a”, “an” and “the” are intended to include the plural forms as well, unless the context clearly indicates otherwise. Furthermore, to the extent that the terms “including”, “includes”, “having”, “has”, “with”, or variants thereof are used in either the detailed description and/or the claims, such terms are intended to be inclusive in a manner similar to the term “comprising.”

[0217] Unless otherwise defined, all terms (including technical and scientific terms) used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs. It will be further understood that terms, such as those defined in commonly used dictionaries, should be interpreted as having a meaning that is consistent with their meaning in the context of the relevant art and will not be interpreted in an idealized or overly formal sense unless expressly so defined herein.