A Sliding-Rod Variable-Strain Model for Concentric Tube Robots

In this work, the Piecewise Variable-strain (PVS) approach is applied to the case of Concentric Tube Robots (CTRs) and extended to include the tubes’ sliding motion. In particular, the currently accepted continuous Cosserat rod model is discretized onto a finite set of strain basis functions. At the same time, the insertion and rotation motions of the tubes are included as generalized coordinates instead of boundary kinematic conditions. Doing so, we obtain a minimum set of closed-form algebraic equations that can be solved not only for the shape variables but also for the actuation forces and torques for the first time. This new approach opens the way to torque-controlled CTRs, which is poised to enhance elastic stability and improve interaction forces’ control at the end-effector.


I. INTRODUCTION
C ONTINUUM robots can be considered a class of soft manipulators, particularly suited for Minimally Invasive Surgery (MIS) [1].One of the most promising continuum robotic systems developed so far is the Concentric Tube Robot (CTR), a collection of nested millimeters tubes whose elastic interaction is used to control the system's overall shape and the iteration force at the end-effector for surgical intervention [2].CTRs have demonstrated promising results in a variety of MIS applications.For a thorough survey of the clinical applications of CTRs, see [3] and the references therein.
The development of CTRs has been made possible thanks to rapid advancement in modeling these kinds of complex non-linear systems.These research efforts have converged on a select type of Cosserat rods model with a particular kinematic structure representing the internal tubes' additional rotational motion [4], [5].Although this Cosserat approach has been improved over the years to include tubes' clearance, friction [6], and inertial dynamics [7], it has been essentially maintained as initially proposed [8].Here, we refer to it as the standard CTR model, and we briefly recall it in section II.
On the other hand, a general soft robot is composed of flexible and rigid elements arranged in a parallel or serial fashion and capable of locomoting in the environment.A novel coordinate system has been proposed to model these kinds of systems, which discretizes the continuous Cosserat model of the flexible components onto a finite set of strain basis functions [9], [10].This Piecewise Variable-strain (PVS) approach is a generalization of traditional robotics' geometric model [11] to the case of highly flexible or soft robots [12].Thus, it provides the theoretical framework for applying traditional control strategies to the field of soft robots.One of the aims of the present work is to apply the PVS approach to CTR systems.In this way, we further extend the general PVS approach to this essential soft robotic technology.At the same time, we open new prospects for the control design of CTRs.In particular, the PVS model provides the equilibrium equations as a minimum set of closed-form algebraic equations, which may be easier to handle for control, design optimization, and stability assessment.
Furthermore, in this work, we tackle one of the main limitations of the standard CTR model concerned with the sliding structure problem, also known as the spaghetti problem [13].In the standard model, the tube rotations and insertions dictate the boundary conditions and locations of the section domains.While the rotational actuator torques are easily computed from the model solution, the insertion forces are not present in the model equations.The recent extension to dynamics also does not model the insertion forces and assumes negligible insertion velocity and acceleration [7].Here we relax this assumption, although still in a static setting, which allows calculation of the actuation insertion forces required for the equilibrium for the first time.The spaghetti problem appears in the planar models of animal locomotion [14], structural stability [15], and finite element analysis [16].Thus, to the authors' knowledge, it is also the first time the spaghetti problem is tackled in three dimensions and with an additional relative rotation.

II. MODEL PRELIMINARIES
Before applying the PVS approach to solve the CTR equilibrium, we revise the standard CTR model [4], [5] using geometric notation.For simplicity, we will present the case of fully overlapping tubes here.Note that additional details about the CTR kinematics will be provided for the more general cases analyzed in section III.Furthermore, a new kinematic model that includes the sliding and rotation of the the bases of the tubes as generalized coordinates of the system will be introduced in section IV.

A. Kinematics
Each tube of a concentric tube system can be modeled as a Cosserat rod, a continuous set of rigid cross sections of infinitesimal thickness along a material curvilinear abscissa X ∈ [0, L] where L is the total length of the rod.Identifying each rigid cross-section with the moving frame rigidly attached to it, the configuration of the tube j is completely defined by a curve g j (•) : X → g j (X) ∈ SE(3) represented by the homogeneous matrix: where r j (X) ∈ R 3 is the position vector of the origin of the moving frame and R j (X) ∈ SO(3) is the rotation matrix representing its orientation with respect to the spatial frame.
The high aspect ratio and material property of the conventional CTR allow assuming inextensibility and shearless deformation (Kirchhoff-Love kinematics).Furthermore, neglecting tube clearance, the rods are perfectly concentric.Thus, their centerlines must be the same at any configuration.Then, it follows that an inner tube j differs from its outer tube j − 1 only by a rotation around the tangent vector to the centerline [5], defined by: where θ j (X) is the rotation angle and we have assumed a local coordinate frame as shown in Figure 1.Thus, for an ordered set of tubes we obtain: where N is the total number of tubes and the outer tube is numbered with 1.A pictorial representation of this kinematics is shown in Figure 1.Note that in this model, the concentric tubes share the same material abscissa X. Let's obtain the time (∂/∂t = ˙) and space (∂/∂X = ) derivative of the tubes configuration g with the Lie algebra se(3).Starting with tube 1 we get: where v(X), u(X) ∈ R 3 represent the linear strains and velocity respectively, while w(X), k(X) ∈ R 3 are the angular strains and velocity respectively.All the quantities are defined in the local frame at X. To indicate the angular strain in the reference stress-free configuration, we use the notation k * .Finally, the superscript ∨ indicates the isomorphism between the Lie algebra se(3) and R 6 [11] (∧ will be used in the opposite direction).
Equating the mixed partial derivatives of the tube configuration, we obtain a compatibility equation between the velocity and strain twists.
As shown in [9], integrating (6) with respect to space yields the following useful relation: where the operator Ad is the Adjoint map in SE(3) defined in Appendix A.
For what concerns the inner tubes, first we define the derivatives of the relative rotation g θj .
The same procedure that led to equation (7), in this case yields: Finally, the chain rule together with (4), ( 5), (8), and (9) yields, for a general tube j:

B. Statics
In this work, we focus our attention on the static equilibrium of concentric tubes with no external applied force (except for the actuation in section IV) to get more neat results.However, it should be noted that there are no particular impediments that would prevent extending the approach to the general dynamic case.We have also neglected friction forces and torques interaction between the tubes as usually done with the standard model.
The static equilibrium of a Cosserat rod in a concentric tube setting can be derived from Hamilton's principle (see [17] for a derivation) with the addition of the constraints forces due to the concentricity constraint.Thus, for a tube j we obtain: where T ∈ R 6 is the vector of internal moment and force, F λj (X) ∈ R 6 is the vector of distributed constraint force, and ad * is the co-adjoint map in SE(3) defined in Appendix A. Note that the value of the internal force N ij is unknown due to the inextensibility and shearless constraints.On the other hand, the internal moment's value can be computed from a constitutive law that we assume linear for simplicity (Hooke's law).
where Σ = diag(GJ, EI, EI) ∈ R 3×3 is the elasticity matrix, E is the Young modulus, G the shear modulus, and J, I are, respectively, the polar and bending second moment of area of the circular cross-section.For what concerns the constraint contact force, although its value is also unknown, its basis can be obtained from the virtual constraints equations that enforce the concentricity constraint.These equations can be written in the form where j j and k j point in the y−, and z−axis directions of tube j at X (see Fig. 1) , and δ is a virtual displacement.Equation (15) says that the distributed constraint force is equal and opposite between two consecutive tubes and that it takes the form: where λ y (X) and λ z (X) are Lagrange multipliers, and Ad * is the co-Adjoint map in SE(3) defined in Appendix A.
The static equilibrium of a collection of tubes imposes the following boundary conditions on equation ( 13) [5].
(17) Note that ( 17) also holds at discontinuity points X where the tubes stop to overlap (in that case, the upper limit of the sum on the RHS will be smaller).Furthermore, since the concentric tubes can not transmit axial torque to one another, the following condition holds for each tube [5].
where i points in the x−axis direction, tangent to the midline (see Fig. 1).

III. VARIABLE-STRAIN CTR MODEL
In [9], [10], a novel variable-strain approach has been presented to model soft manipulators driven by tendons and pneumatic chambers.This section aims to apply this technique to CTRs to obtain a minimum set of closed-form equations describing the system equilibrium.
According to the PVS approach, the configuration g j of the collection of tubes is represented by the strain fields ξ 1 and ξ θj .In particular, the infinite dimensional strain fields are discretized on a finite set of basis functions as follows.
where B 1 (X) ∈ R 6×n1 and B θj (X) ∈ R 6×nθj are matrix functions whose columns form the basis for the strain field ξ 1 and ξ θj , respectively, while p 1 ∈ R n1 and p θj ∈ R nθj are the vectors of coordinates.Note that, due to the assumptions made in section II, ξ is equal to [0 0 0 1 0 0], while the last three rows of B 1 (X) and the last five rows of B θj (X) are all equal to zero.Then, the set of generalized coordinates become q = (p 1 , p θ2 , . . ., p θ N ), and the configuration g j (q) can be reconstructed through the integration of equations ( 5), ( 9) and the recursive formula given by (3).The differential relation between configuration and generalized coordinates is obtained by replacing equations ( 19), (20) in the velocity equations ( 7), (10), which yields and, in turn, replacing the result in the velocity equation (11).Finally, we obtain: Note that J θi (X) is an analytical function that can be computed offline given the choice of basis B θi .The differential equation (23) provides the required Jacobians to project the static equilibrium ( 13) by d'Alembert's principle.In particular, equation ( 13) is projected with dX for all j and i.The results are then summed over all j from 1 to N to obtain as much algebraic equilibrium equation as the dimension of q.This procedure is shown in detail for the case of two fully overlapping tubes in the next section.The general case of multiple non-overlapping tubes with straight actuation, as shown in Figure 3, will be illustrated in the following.

A. Single overlapping sections
Let us consider the case of two fully overlapping tubes.Applying d'Alembert's principle as described above yields: Considering the form of the constraint force ( 16) and the identity Ad * gj F ij , we get (note that Ad T g = Ad * g −1 ): Finally, integrating by part we obtain: where we have used the strain equation ( 12), the identity F i2 = 0 and the boundary conditions ( 17), (18) at L. Note that, since the generalized coordinates q are independent, the unknown constraint forces F λj cancel out as expected.
Equation (26) can be solved numerically for the unknowns p 1 and p θ2 that appear in F ij , ξ 1 , and g θ2 through equations ( 19), ( 20), (14), and the integration of ( 9).This integration becomes simple in this case, since the rotation angle θ of (2) can be obtained by: where () x extracts the first rotational element of a vector and α 2 is the rotation of the tube imposed at the base by the actuators.Note that the integrals of equation ( 26) can be expressed analytically with the use of hypergeometric functions, which makes their evaluation quick for numerical root finding.

Comparison with analytic torsionless case
As a first test to check the validity of the proposed model, we consider here the case of two torsionless, fully overlapping tubes with three-dimensional deformation.The torsionless approximation was used at the beginning of the CTR development, and analytic solutions are available for the constant curvature approximation [18], [19], [20] as well as the general variable curvature case [21].
where k * j is the reference strain of tube j.Thus, n 1 = 4 and the generalized coordinates T represent the amount of reference curvatures making up the equilibrium configuration.Rearranging (28), we can write: A similar equation can be obtained for q 3 and q 4 .Solving for q, we get: which correctly corresponds to the analytic solution k

B. Multiple non-overlapping sections
Let us now move to the case of multiple non-overlapping tubes.To address this problem, we divide the system domain into piecewise variable-strain sections corresponding to the curvature discontinuity located at the end of every overlapping portion.Thus, the kinematic equation ( 3) is applicable locally at each overlapping portion and can be generalized as follows (see Figure 1): where g js indicates the configuration of tube j with respect to the spatial frame, while g j represents the configuration with respect to the local overlapping portion.L i is the length of tube i, which corresponds to the end of an overlapping portion.g θj s (X) is the total rotation between tube j and tube j − 1 at X, which is obtained through the product of the partial rotation of each overlapping portion preceding X.For example, consider the total rotation between tube j and tube j − 1 at X located on the third overlapping section, then g θj s (X) = g θj 1 (L 1 )g θj 2 (L 2 )g θj 3 (X).
A procedure similar to the single overlapping sections case leads to the equilibrium equations for general multiple nonoverlapping tubes in the generalized variable-strain coordinates.For example, consider the case of three non-overlapping tubes, then we obtain: Each dashed line separates an overlapping section.As with (26), equation (33) can be solved numerically for the unknown q = p 1 , p θ2 1 , p θ3 1 , p 2 , p θ3 2 , p 3 .We note that the model grows linearly with the number of tubes, and number of degrees of freedom.

Comparison with simulation in literature
The multiple overlapping tubes model is tested by comparing the equilibrium configuration obtained for a two tubes robot with a literature benchmark reported in [21].Consider two tubes of length 140 mm (outer) and 200 mm (inner) with constant pre-curvature and physical properties reported in Table I.We seek for the equilibrium configuration when the inner tube is rotated by α 2 = 180 • .The system domain is divided into two variable-strain sections.The generalized coordinates are q = (p 1 , p θ2 , p 2 ), and the equilibrium equations are a subset of (33).The resulting model takes the form of a nonlinear algebraic system of equations that can be solved numerically with the MATLAB © fsolve function.We chose a linear nonhomogeneous base for the strains involved in the first section, i.e.: k 1 , and θ 2 .Then, we have 1 X 0 0 0 0 0 0 1 X 0 0 0 0 0 0 1 X 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ) For the second section, since it is unloaded, we use the reference strain k * 2 as a basis for the strain k 2 , which gives p 2 = 1.As illustrated in [21], the CTR has three equilibrium configurations.The solution to which the solver converges depends on the initial guess.Figure 2 illustrates the equilibrium configurations obtained with the variable-strain model.Comparing the Cartesian tip positions with respect to the standard model, we find a maximum error of 1.17mm, and a mean error of 7.8 * 10 −1 mm, which correspond to 0.58% and 0.39% of the backbone length, respectively.Furthermore, the inner tube's rotation with respect to the outer tube at the end of the overlapping section is equal to θ 2 (L 1 ) = ±85.5 • , which is only 1.3% bigger than the value reported in [21] (±84.4 • ).

C. Complete CTR system
The CTR systems considered so far only include the tubes' portions that follow the insertion orifice, assuming that the rotation input α was applied there (see equation ( 27)).However, a complete CTR system is also composed of a straight portion that precedes the insertion point and terminates with an actuation system responsible for applying the input angle and insertion motion.To include this part, we place the 0 of the system domain X at the innermost tube base.Then, we divide the system domain into several sections corresponding to the different discontinuities.Those include a change in the number of overlapping tubes (as before), the insertion orifice, and, eventually, a reference strain jump.For instance, the straight portion is usually obtained by letting the tube be originally straight at the proximal end and curved toward the distal end.Finally, we constrain the straight portion by choosing only a torsional strain basis for these proximal sections.Note that this method applies to fully curved tubes as well, where external constraints enforce the straightness condition.
For example, let us consider the two tubes CTR system shown in Figure 3, where each tube has a straight proximal reference strain of length L js as indicated.There are five sections in this example indicated with an additional subscript in the following.Define D 1 , D 2 to be the distance between the insertion orifice and the outer and inner tube base, respectively, and ∆D = D 2 − D 1 .Then, the equilibrium equations for the generalized coordinates q = p θ2 1 , p 12 , p θ2 2 , p 13 , p θ2 3 , p 14 , p θ2 4 , p 25 are given by: (35) In developing the equilibrium equation ( 35), we have used the identity J T θ ad * ξ F i = 0 for the straight portion that precedes the insertion orifice.Furthermore, note that B 12 has only torsional components, and B θ2 2 (X), B θ2 3 (X) are equal to zero at X < ∆D and X < D 2 , respectively.
Considering the CTR systems of equation ( 33) and ( 35), we notice that each section follows the same structure of the single overlapping case (26).

Comparison with experiments in literature
The results of the complete system model ( 35) are compared with experimentally validated simulations available in the literature [21].Consider two tubes with a straight portion followed by a constant pre-curvature, as shown in Figure 3, and physical properties reported in Table I.We seek the equilibrium configuration in two sets of insertion and rotation conditions of the inner tube.In the first set, the full overlap case, D 2 = 208.5 mm, and α 2 goes from 0 • to 280 • in 20 • increments.In the second set, the partial overlap case, D 2 = 170.5 mm, and α 2 goes from 0 • to 200 • in 20 • increments.D 1 and α 1 are always equal to 93.5 mm and 0, respectively.We chose a quadratic basis for all the strains involved, except for the last section, where we use k * 2 .Figure 4 illustrates the equilibrium configurations obtained with the variable-strain model.Comparing the Cartesian tip positions with respect to the standard model, we find a maximum error of 4.8 * 10 −1 mm, 8.9 * 10 −1 mm and a mean error of 4.4 * 10 −1 mm, 7.2 * 10 −1 mm, respectively, for the full overlap and partial overlap case.The mean errors correspond to 0.14% and 0.23% of the backbone length.

D. Benefits of the PVS model for CTR
The proposed piecewise variable-strain approach for concentric tube robots presents several benefits.It provides the equilibrium equations as a minimum set of closed-form algebraic equations.The system is easy to handle for control and design optimization since all the boundary and continuity conditions are already intrinsically embedded in the equilibrium equations.Furthermore, extension to dynamics can be done without any change to the number and type of generalized coordinates.Finally, the PVS model is well-posed to accommodate additional DOFs and incorporate design variation of CTRs as well as generalization to other rod-driven soft robots.In the next section, we include the tubes' insertion and rotation motion as generalized coordinates of the system for the first time.

IV. SLIDING CTR MODEL
In the standard model approach, the tube rotations and insertions are actuation variables that dictate the boundary conditions of the governing ODEs.Actuator motion is thus simulated by simply varying the boundary conditions and solving the static model at sampled points in time [4], [5], thus producing what is known as a "quasi-static" simulation.While the Jacobian can be numerically computed along the robot [22], the tube velocities are not automatically included in the standard model.The recent extension of the standard model to dynamics does describe velocities, but still assumes negligible actuator insertion velocity and acceleration [7].Further, while the torque required to rotate the tubes is computable from the solution boundary values and the constitutive law, the standard model does not include information about actuator insertion forces.This section will extend the proposed PVS model for CTR to include the input motions as generalized coordinates and overcome the standard model's above limitations.
Before that, we establish a general result for the differential kinematics of rod's sections with a variable domain.Consider a rod section with material abscissa X ∈ [a(t), b(t)] and g(a) = I.The tip of the section g(b) can be formally represented by the integration of the linear time-varying differential equation ( 4) through its state-transition matrix g(b) = Φ(b, a).Then, the derivative of g(b) with respect to time has to include the variation of the domain boundaries in addition to the variation of the shape.This can be calculated as follows.
where the first can be considered a "growing from the base" term and the second a "growing from the tip" term.Note that internal cross-sections g(X) will vary only due to the proximal boundary a, while they are indifferent to the distal boundary b.Finally, the total velocity twist becomes: Observe the differences with equation ( 7) and ( 10) used so far.
To keep track of the potentially varying boundaries, in the following, we replace the notation g(b) with g([a, b]).

A. Complete CTR system with actuation forces and torques
Let us consider, without loss of generality, two nonoverlapping concentric tubes constrained to be straight before the insertion orifice, as shown in Figure 3.We removed the reference strain discontinuity here for simplicity, reducing the number of required sections to four.First, we define two material abscissas respectively for the outer tube X 1 = [0, L 1 ], and the inner tube X 2 = [0, L 2 ].Overlapping cross-sections are then related by: where D 1 and D 2 are generalized coordinates now.In this example, the full set of generalized coordinates is q = D 1 , D 2 , α 1 , α 2 , p θ2 1 , p 12 , p θ2 2 , p 13 , p θ2 3 , p 24 .Contrary to section III-C, we fix the spatial frame on the insertion orifice.Then, the tubes' sections' kinematics can be expressed with equation (32) after being pre-multiplied by a translational homogeneous matrix g t , which depends on the insertions D 1 and D 2 .Let us focus on the overlapping section immediately after the insertion orifice.The kinematics there can be written as: where, for the inner tube, X 1 is given by (38).The last term g θ2 3 ([D 2 , X 2 ], X 1 (∆D)) indicates the relative rotation between X 1 (∆D) of the outer tube and X 2 of the inner tube accumulated in the last section starting from D 2 .Using the general formula (37) and the discretizations (19), and (20), we obtain the differential kinematics of the section.
where B α = [1 0 0 0 0 0] T , and ξk = [0 (k) y (k) z 1 0 0] T .Note that the variable upper boundary of g θ2 3 is due to X 1 (∆D).Thus, its variation has to be taken with X 2 fixed.Equation (40) provides the additional Jacobians corresponding to the input motions.Similar equations can be obtained for the other three sections.Projecting the differential equation ( 13) by d'Alembert's principle using these additional Jacobians yields the equilibrium equations for the insertion and rotation input forces τ D1 , τ D2 , and torques τ α1 , τ α2 .
(41) where ξ k = [0 (k) y (k) z 0 0 0] T .To obtain (41) we have also used the following boundary conditions. ξT According to the Newton law of conservation of momentum, the input forces' and torques' sums have to be zero.It is an excellent verification to check if equation (41) satisfies these conditions.For the input torques, the third equation of (41) ensures that the actuation torques are equal and opposite.For what concerns the input forces, we have: where the last identity, which states that the internal stress is null on the cross-section immediately after the orifice, is justified by the fact that no external forces act on the system from that point onward.Thus, also the input forces are equal and opposite as required.

B. Simulation tests
Some simulation tests are presented in this section to explore the behavior of the actuation forces and torques.We consider six configurations (three planar and three out-of-plane) of two concentric tubes with physical properties as for Table I.In all the cases, we consider linear non-homogeneous (or constant) pre-curvatures.The linear and non-homogeneous pre-curvature coefficients are varied to study the effect on the actuation forces and torques.In the out-of-plane cases, the inner tube is rotated by 90 • with respect to the outer tube.Table II reports the exact values of the non-homogeneous and linear coefficients and the calculated values of the actuation inputs.
The following observations can be made.For case a), the insertion motions D 1 and D 2 do not influence the actuation inputs value.In case b), while the insertion of the inner tube D 2  , and it is symmetric with respect to clockwise or counterclockwise rotations.In particular, the linear coefficients of e) and f) have an opposite contribution to the actuation force τ D2 .In case e), τ D2 decreases with the increment of the linear coefficient until it becomes negative, as reported in Table II.In case f), the opposite applies.

C. Benefits of the sliding-rod PVS model for CTR
The sliding-rod PVS model presented here proposes a growing (non-material) approach, which allows extending the PVS model to include the tubes' sliding motion without the need of calculating the unknown interaction forces.This new model may be used to control the CTR motion through the actuation force and torques instead of the insertion and rotation kinematics.Torque-controlled CTR can provide a new way to enhance elastic stability and improve interaction forces' control at the end-effector, currently a major concern in the design of CTR [23].In the present form, the base force equations (41) can also be used to improve the fast kinematic controllers of [24], which is based on actuation load sensing.

V. CONCLUSIONS
In this paper, the recently proposed piecewise variablestrain approach for modeling highly deformable rods has been adjusted and applied to the case of concentric tube robots.The performances of the new approach have been compared with analytic, simulated, and experimental data available in the literature.Furthermore, the PVS approach has been further extended to include the tubes' insertion motion for the first time, which opens new unexplored possibilities for controlling these kinds of systems.
Future works include experimental validations of the slidingrod PVS model, the addition of external forces, and the extension to dynamics.

Fig. 2 :
Fig. 2: Equilibrium configurations computed with the variable-strain model for two concentric tubes with constant pre-curvature.The inner tube is rotated by 180 • .The result of the variable-strain model match well with the one off the standard CTR model [21].

Fig. 4 :
Fig. 4: Equilibrium configurations computed with the variable-strain model for two concentric tubes with straight portion and constant pre-curvature.The inner tube is rotated from 0 • to 280 • in 20 • increments (right) and from 0 • to 200 • in 20 • increments (left).The result of the variable-strain model match well with the one of the standard CTR model [21].

TABLE I :
Physical Properties of Tubes

TABLE II :
Sliding-rod test results , the force input absolute value |τ D | increase with the outer tube's insertion and curvature.In c), the force input absolute value |τ D | increase with the insertion of the outer tube and the retraction of the inner tube.For the out-of-plane cases, in general, the absolute value |τ D | increases with the rotation angle α 2 up to 90