A continuum robot backbone model is developed, based on the classical nonlinear Cosserat rod theory. The backbone can be described effectively by an elastic rod in three dimensions, which accounts for the nonlinearity of rod bending [see Alqumsan et al. (2019), Robinson and Davies (1999)]. This approach discounts the previously common modeling assumption of constant curvature of the backbone.

The centroid position along the backbone arc length s ∈ [ 0 , l ] together with its local orientation with respect to the global origin and coordinate system are given by r ( s ) and rotation matrix R ( s ) , respectively, as shown in Figure 1 with a Cartesian coordinate system. The backbone centroid is as depicted in Figure 1 and the local linear and angular rates of change of the backbone are denoted by u ( s ) and Ω ( s ) , respectively. The corresponding equations are given as r ( s ) = ( r x r y r z ) ,       u ( s ) = ( u x u y u z ) ,       Ω ⌢ ( s ) = ( 0 − Ω z Ω y Ω z 0 − Ω x − Ω y Ω x 0 ) , R ( s )   = ( cos ⁡ ϕ ⁡ cos ⁡ θ cos ⁡ ψ ⁡ sin ⁡ θ ⁡ sin ⁡ ϕ cos ⁡ ψ ⁡ sin ⁡ θ ⁡ cos ⁡ ϕ   +   sin ⁡ ψ ⁡ sin ⁡ ϕ sin ⁡ ψ ⁡ cos ⁡ θ sin ⁡ ψ ⁡ sin ⁡ θ ⁡ sin ⁡ ϕ   +   cos ⁡ ψ ⁡ cos ⁡ ϕ sin ⁡ ϕ ⁡ sin ⁡ θ ⁡ cos ⁡ ϕ   −   cos ⁡ ψ ⁡ sin ⁡ ϕ − sin θ cos ⁡ θ     sin ⁡ ϕ cos ⁡ θ     cos ⁡ ϕ ) , (1) where ψ, θ, and ϕ refer to the angles through which the local axes have turned about the global x, y, and z axes, respectively, and the hat denotes the skew matrix form. The transformation of the backbone along its arc length s can therefore be defined by r ˙ ( s ) = R ( s ) u ( s ) ,                   R ˙ ( s ) = R ( s ) Ω ⌢ ( s ) , (2) where R ( s ) acts to transform the points from a local to global coordinate system and R ˙ denotes the derivative of R with respect to the arc length s.

The backbone internal force and moments in global coordinates are given by n ( s ) = ( n x n y n z ) ,                     m ( s ) = ( m x m y m z ) . (3) The distributed load applied to the backbone f ( s ) is given by the negative of the derivative of n ( s ) with respect to s, which is taken to be the weight of the backbone. Similarly, due to the internal force n ( s ) acting at a given global coordinate r ( s ) , as in Figure 1, the derivative of m ( s ) with respect to s can be equated to the negative of the cross product of r ˙ ( s ) and n ( s ) . This results in the set of equations for the rate of change of internal forces and moments n ˙ ( s ) = − f ( s ) ,                     m ˙ ( s ) = − r ˙ ( s ) × n ( s ) . (4) The constitutive law dictates that the backbone internal forces and moments can be described by the difference between the deformed and undeformed states at arc length s along the backbone. The undeformed state is given by subscript * and is assumed to be a straight cylindrical rod extending in the local z axis, giving u * = ( 0,0,1 ) T and Ω * = ( 0,0,0 ) T . The internal forces and moments depend on the backbone stiffness in shear and extension for the forces and bending and torsion for the moments, characterized by matrices K s e and K b t , respectively; n ( s ) = R ( s ) K s e ( u ( s ) − u * ) ,                     m ( s ) = R ( s ) K b t ( Ω ( s ) − Ω * ) , (5) where K s e ( s ) = ( G A 0 0 0 G A 0 0 0 E A ) , K b t ( s ) = ( E I x 0 0 0 E I y 0 0 0 E I z ) , u * = ( 0 0 1 ) , Ω * = ( 0 0 0 ) . (6) G denotes the shear modulus, A the cross sectional area of the backbone, E the Young’s Modulus and I the second moment of area about the subscript coordinate.

Therefore, the set of governing equations for the global coordinates, formed from Eqs. 2–5, will be r ˙ ( s ) = R ( s ) u ( s ) ,                     R ˙ ( s ) = R ( s ) Ω ^ ( s ) ,                     n ˙ ( s ) = − f ( s ) ,                     m ( s ) = − r ˙ ( s ) × n ( s ) , (7) with boundary conditions r ( 0 ) = ( 0,0,0 ) ,                     R ( 0 ) = ( 0 0 0 0 0 0 0 0 0 ) , n ( l ) = F ,                     m ( l ) = L , (8) where F and L are the forces and moments applied at the backbone tip, respectively. The local coordinates are computed using u ( s ) = K s e − 1 R T ( s ) n ( s ) + u * ,                     Ω ( s ) = K b t − 1 R T ( s ) m ( s ) + Ω * . (9) The set of governing Eq. 7 with corresponding boundary conditions 8) are solved simultaneously using a boundary value problem solver based on a finite difference code that implements the four-stage Lobatto IIIa formula which has high mathematical efficiency and robustness (see Shampine and Kierzenka (2008) for details).

The experimental set-up for evaluating the model can be seen in Figure 2. The backbone is made of ASTMA228 spring steel and has a length of 240 mm, chosen such that large curvature can be achieved without plastic deformations. Four 0.28 mm 4-strands Dyneema (Dingbear, China) tendons are chosen to minimize elongation of the tendon itself while providing high fracture stress. The tendons are numbered 1–4 in a clockwise direction, when looking end on, starting from the middle right tendon, giving tendons 1 and 3 on the horizontal plane, and tendons 2 and 4 on the vertical plane, see Figure 6.

Acrylic circular support discs were used to minimize friction between the tendon and support discs; additionally Teflon spray was used to further reduce friction. A total of 13 support discs of 20 mm diameter and 1.5 mm thickness were used and spaced 20 mm apart along the backbone, such that an optimum ratio of support disc spacing to tendon offset from backbone of 2.5 mm was achieved as proposed by Rucker and Webster (2011). The tendons passed through the support discs at a distance of 8 mm from the backbone central axis and are equally spaced around the circle circumference.

Each tendon is actuated independently, through an assembly comprising of a lead screw mechanism driven by a stepper motor. The stepper motors have a resolution of 0.05625 ∘ and the lead screw provides 0.008 m of translation per rotation, resulting in a linear resolution of 1.25 × 10 − 5   mm per step. To stabilize the system and ensure linear translation of the lead crew nut, a v-slot linear rail and carriage system is used, together with its ability to translate with minimal friction which is necessary for the high tendon tension case.

A cantilever load cell is fixed to the carriage, which in turn is connected to the tendon, causing it to move through a set of forces and a moment offering a loading feedback for each of the tendons. A 1 kg load cell (Phidgets, Canada) was chosen to identify the tension as a compromise between the maximum expected tendon tension and resolution.

The position of the backbone was recorded using a contactless sensor, namely OptiTrack (Target3D, United Kingdom) motion capture system; this eliminated any error that may be associated with contact sensors. Reflective markers of 6 mm diameter were incorporated on the face of every second support disc. Seven cameras were used, to ensure uninterrupted tracking of the markers, with a maximum root mean square positional error of 0.105 mm achieved. The initial position of the backbone, before any tendon loading, was identified once a stable tension reading of the desired value ± 0.01 N had been achieved.

The backbone model in Section 2.1 is modified to represent the experimental set-up with the tendons coupled to the backbone as shown in Figure 3 using support discs to provide actuation of the backbone; moments on the end support disk of the backbone must be included, which are denoted by L i , where i = 1 − 4 for tendons 1 − 4 , respectively. The tendons also produce an axial point force at the tip location, however, these are assumed negligible as the distributed force that the tendons apply to the support discs are at least an order of magnitude larger Gravagne et al. (2003). Nevertheless, the boundary conditions can be adjusted for any externally applied end load through the internal force n ( l ) . Therefore, by neglecting the distributed forces from the tendon onto the backbone, and only considering the moment they impose on the backbone tip, the number of numerical computations that are needed at each iteration.

When the backbone is undeformed, it is assumed moments L 1 and L 3 act around the global z axis and L 2 and L 4 about the y axis as shown by Figure 1. Each applied tendon moment is calculated using L i = T i h , where T i is the tendon tension on the i th tendon and h is the tendon offset from the backbone. This leads to the boundary condition for the internal moment at the tip location m ( l ) = L = ( 0 L 4 − L 2 L 1 − L 3 ) . (10) Therefore, the boundary conditions in Eq. 8 are modified to accommodate Eq. 10.

Based on the boundary conditions in Eq. 8 and representation model of experimental facility in Section 2.3, the control variable is selected to be the tension of the tendons, T i . This firstly offers a simpler control approach and secondly ensures that any frictional effects are directly countered by the motor. The closed loop control of the tendon tension was implemented as a PID controller on an Arduino Uno (Arduino, Italy), see Figure 4. The Arduino Uno receives demand tension values from the model via serial communication from an i7-7500U, with 8 Gb memory running Matlab (Mathworks, United States) on Windows 10. The demand values are held via a zero-order hold and the control loop is run at 100 Hz. The tension feedback was received from the load cell via an integrated amplifier/ADC NAU7802 (Nuvoton, Taiwan) via an i2c bus interface, sampled at the control loop frequency.

A PID controller architecture was selected given the overall behavior of the system. Although PD is the preferred method in continuum robots George Thuruthel et al. (2018) it can be observed that PD control guarantees stability only when the PD gains tend to infinity Slotine and Li (1991). Moreover the steady state error does not tend to zero unless there is some form of model-based compensation. Since the presented model does not involve dynamic behavior and the tension controller, inclusions of the Integral term was done to ensure minimizing the steady state error by compensating tendon friction.

Initially, the tuning procedure for the PID controller followed the classic Ziegler-Nichols method based on a step response of a single tendon being tensioned from 0 to 3 N. This procedure resulted in zero steady state error and minimal oscillations with the calculated PID gains for the proportional, integral and derivative terms of K p = 0.15 , K i = 3 , and K d = 0.15 , respectively. However, as shown in Figure 5, it resulted in a rise time of 70 s was which was deemed insufficiently long.

Therefore, to reduce the system rise time to 25 s, a two-stage proportional controller was implemented. This is similar to a divide and conquer gain-scheduling controller Rugh and Shamma (2000) with the only evaluation criterion being the desired rise time and tendon tension. In stage one an initial high value of K p = 1 was used. This value then switched in the second-stage to the value identified by the Ziegler-Nichols method ( K p = 0.15 ), to maintain the desired steady state characteristics. The threshold at which this switch occurred was determined to be relative to the desired tendon tension ( T d e s ) at a value of T d e s ± 0.2 N. The integral, and derivative gains remained the same for the two control stages.

In order to evaluate the validity of the results produced by the experimental set-up, the effect of tendon-disk interface friction on the system repeatability was investigated. The tension in tendons 1 and 3 were kept at 0 N while the tension in tendons 2 and 4 were changed, but the tension in only one tendon was non-zero at any given time, to a value of either 1.5 N or 3 N. Experiments were carried out where these values were achieved by increasing the tension from 0 N as well as by over tensioning the tendon to 5 N before decreasing it to the desired value which allowed for the removal of tendon slack from within the support discs. The relative distance of the backbone tip from the origin was recorded and analyzed for each repeat.

The results of the relative distance of the backbone tip from the origin for different tensions in the tendons in the vertical plane are shown in Figure 7 when the desired value is achieved by increasing the tension from 0 N (red) or decreasing from 5 N (blue). Table 1 summarizes the standard deviation for each test case. From these results it can be seen that when the tension is 0 N in both tendon 2 and 4 the experimental facility is consistent in achieving the same position when either increasing or decreasing the tension to 0 N. However, for all other cases examined, when the desired tension is achieved by decreasing it from 5 N the backbone tip position is further from the fixed end at the origin than when the tension is increasing from 0 N. This indicates that the backbone has a greater amount of bending when the tension is achieved by decreasing it from 5 N for all loading cases. The difference between the two methods of tensioning the tendons increases with the tendon tension magnitude and is attributed to stiction/friction occurring between the tendon and the support discs, so that when the tension is increased from 0 N there is some slack left in the system. Given this observation the subsequent experiments will be tensioned to the desired value, then the backbone will be returned to the neutral position before performing the next force combination.

The impact of the support disks and the OptiTracker markers on the modulus of elasticity of the rod is investigated, to enable model validation. Based on these measurements the updated modulus is used to perform the calculation for the model and used subsequently. To identify the modified Young’s Modulus, values of E were swept in the model from 130 to 230 GPa. For each value, the error (distance) between the predicted backbone tip position and actual tip position was found for the different tendon tensions combinations. The errors for each of the tendon tension combinations were summed, giving a total error value for a given E. The total errors are shown in Figure 8 over the examined range, with the value E = 168 GPa minimizing the total error across the tendon tension range to 0.1422 m, thus this is the value used within the model for the rest of experimentation.

In this fist set of experiments the accuracy of the predictions from the derived model is investigated with no end effector. The updated Young’s Modulus from Section 3.2 and the friction countering approach from Section 3.1 are used. As stated above, it is assumed that the effect of the horizontal plane (tendons 1 and 3) would be symmetrical about 0 N. Each tension combination was measured three times and an average value was taken. Results showing the predicted and realized position of the backbone are given in Figure 9 for the three different regions defined within the tendon combinations; In-Plane, Semi-Out-of-Plane and Out-of-Plane test points. Figure 10 shows the mean positional error against the backbone arc length for the three different regions, at the sensor locations. In general there is an increase in the mean positional error as the backbone arc length increases. The In-Plane region has a significantly lower mean positional error with the Semi-Out-of-Plane region having a smaller error than the Out-of-Plane; this is expected as the backbone diverge from the fixed origin. This is due to the assumption that even for the Out-of-Plane bending, tendons are assumed to apply in plane moments to the end support disc.

Table 2 summarizes the mean positional tip location error together with the standard deviation for the three different regions. For the In-Plane region the error is 7.14 mm or 2.98% of the backbone length with a standard deviation of 3.21 mm. The Semi-Out-of-Plane and Out-of-Plane regions have an increased error of 10.9 and 12.1 mm or 4.56 and 5.05% of the backbone length with a standard deviation of 7.73 and 3.30 mm, respectively. This leads to an overall mean positional tip location error of 9.48 mm or 3.95% across all regions with a standard deviation of 4.99 mm.

In practice a continuum robot will have an end effector attached at the backbone tip location, therefore, investigations are carried out in which a vertical end load is acting on the tip of backbone. In this case the load was attached to the end support disk, where tendon 4 passes through it. The examined tensions in the different tendons are similar to those previously studied but with a reduce number of combinations of only 0 and 3 N for tendons 1, 2 and 4 with tendon 3 remaining slack. Four different loads were applied to the backbone tip of [ 0,96.14 , 189.92 , 293.32 ]   N , with three repeats carried out and averaged.

Figure 11 shows the predicted and realized backbone position for each load of the first three loads. The mean positional error between the predicted and realized positions along the backbone are given in Figure 12 for an increasing load on the backbone tip. Increasing the point load gives an increase in the mean positional error, with the tip error increasing significantly for the horizontal tensioned tendon cases; when tendon 1 has a non-zero tension. Table 3 gives the mean positional tip location error, together with the standard deviation, for increasing end load. Load 293.32 N results in the highest average positional tip error of 100 mm or 41.8% with a standard deviation of 73.3 mm, which is an order of magnitude larger than the smaller examined loads. The overall mean tip error for all cases is 42.3 mm or 17.63%, with a standard deviation of 73.33 mm.