The most commonly used approach to modeling CRs and HRRs is the constant curvature model. CRs and HRRs with multiple segments can be modeled using the piecewise constant curvature model (Burgner-Kahrs et al., 2015). Using the Jacobian matrix to establish differential kinematics is the standard approach for CR/HRR inverse kinematic models (IKMs) (Whitney, 1972) and was also employed by Wolf et al. (2005), Garg et al. (2014), Guardiani et al. (2022), Lee et al. (2014), Li et al. (2017), Rone and Ben-Tzvi (2014), Simaan et al. (2009), Jones and Walker (2006), and Yeshmukhametov et al. (2019).

The generalized coordinates q are updated by the multiplication of the pseudo-inverse of the Jacobian (Chirikjian and Burdick, 1994) by the difference in parameterized end-effector coordinates (Klein and Huang, 1983). Alternative approaches comprise the augmented (Seraji, 1989) and extended (John, 1983) Jacobian.

Bajo et al. (2012) applied a Jacobian-based IKM to their IREP platform, which consists of two independent, two-segment CRs. They computed the Jacobian matrices for each segment individually.

Giorelli et al. (2013) proposed an alternative approach tailored to non-constant curvature CRs, where the static model is governed by nonlinear differential equations, making exact solutions challenging. As analytical solutions are unavailable, the elements of the Jacobian matrix cannot be computed directly. Instead, a feed-forward neural network is trained to learn the IKM for a non-constant curvature cable-driven manipulator. Similar neural network-based methods were presented by Lai et al. (2019) and Thuruthel et al. (2016).

Wei et al. (2023) presented an IKM for a two-segment HRR under the constant curvature assumption. They simplified the complex nonlinear system by reducing it to a single nonlinear equation, resulting in faster solutions.

Xu et al. (2018) adopted an iterative approach to compute the IKM utilizing a least squares method to determine the joint angles from specified cable lengths. Their model incorporated multilevel mappings that describe the relationships between motors and cables, cables and joints, and joints and the robot’s end-effector.

Almanzor et al. (2023) introduced a method that combines a local inverse kinematics formulation in image space with a deep convolutional neural network trained on synthetic data. A 2D hand-drawn image of the desired shape is used as the input, and the network predicts motor commands to drive the robot toward that shape. Closed-loop control is achieved using visual feedback from a webcam.

Wild et al. (2023) developed a novel piecewise dual quaternion algorithm to model multi-section CRs, aiming to improve the traditional pseudo-inverse Jacobian method. The piecewise dual quaternion approach offers reduced computational complexity, faster convergence, and smoother backbone configurations, especially as the number of robot sections increases.

Kim and Wi (2025) introduced a tendon-driven discrete CR using ball–socket joints. Between one and three robot units connected in series were tested by applying proximal tendon tension, while distal tension was gradually increased to induce bending. The resulting bending curves were interpolated using third-to-fifth-order polynomials.

Li et al. (2025) proposed a novel multilevel motion control method for HRRs based on joint angle–cable force cooperative optimization. It comprises a finite-time optimization approach using macro–micro decomposition to solve inverse kinematics. The high-degree-of-freedom (DoF) robot is divided into smaller manipulators and solved using a dual neural network model.

Using the space vector method, the manipulator’s kinematic model is developed to dynamically determine its endpoint position (Wang et al., 2024). The workspace is then generated through the Monte Carlo method. The original search approach is enhanced by introducing an angle-decoupling mechanism between adjacent links to calculate each joint’s rotation angles.

Zhan et al. (2024) introduced a general approach for solving the real-time optimized IKM of redundant robots, while strictly enforcing hard limits in both joint and Cartesian spaces that cannot be violated. Instead of quadratic programming, the method employs constrained linear programming to address the IKM problem. Hard constraints—including joint limits, velocity, and acceleration bounds—are explicitly managed as inequality constraints.

Sheng et al. (2024) presented an IKM method that combines the end-following approach with a segmented, Jacobian-based iterative solver for CRs with redundant DoF and a moving base. The end-following method is first used to generate a smooth, constraint-compliant initial joint configuration by having joints sequentially follow the target points along a planned path. This initial guess is then refined using Jacobian-based inverse kinematics applied to robot segments, thus improving accuracy while reducing the computation time.

The following subsection provides a brief overview of the Jacobian-based IKM.

The state-of-the-art IKMs for CRs and HRRs are based on the Jacobian matrix J . q denotes the generalized (joint) coordinates, and X ( q ) denotes the end-effector pose corresponding to q . ( q ) 0 represents the initial guess of q . X t denotes the target pose. J + represents the Moore–Penrose inverse of the Jacobian as shown in Algorithm 1.

Algorithm 1

Jacobian-based IKM.

1: Initialize t o l ← 0.01 and q ← ( q ) 0

It should be noted that, in general, calculating the orientation difference requires careful handling beyond the simple subtraction of the target and achieved end-effector orientations.

The Jacobian matrix, as defined in Equation 1, is commonly used in the kinematics and dynamics of robotic systems. It relates variations in joint space to corresponding changes in configuration space. J q = ∂ X 1 ∂ q 1 ⋯ ∂ X 1 ∂ q n j ⋮ ⋱ ⋮ ∂ X m ∂ q 1 ⋯ ∂ X m ∂ q n j . (1)

The forward kinematic model of the robot describes the transformation matrices from the inertial frame I to the end-effector frame E , as provided in Equation 2. T I E q = T I 0 ⋅ ∏ k = 1 m T k − 1 , k q k ⋅ T m E . (2)

The Moore–Penrose inverse is computed using Equation 3. J + = J T J − 1 J T . (3)

Due to the iterative nature of the algorithm, the error converges against a predefined threshold that is considered acceptable by the user. The gain k increases the step size per update and influences the convergence speed.

This study introduces a novel IKM that departs from the traditional Jacobian-based approach. Instead, it relies on the numerical solution of a set of equations representing geometric constraints, referred to hereafter as NSGC. These equations are formulated using the unknown parameters of the arc representing the curvature of the bent CR/HRR. Solving this system enables the calculation of the robot’s length and bending parameters, which are subsequently converted into generalized coordinates q .

The target pose X t = [ x t , y t , z t , ψ t ] is given in 3D space. x t , y t and z t are the position coordinates defined in the global CS, while ψ t describes the end-effector orientation about the y ε -axis expressed in the local end-effector CS. The target pose must be projected onto a 2D plane, i.e., the global x–z-plane with proj ( X t ) = X ̂ t = [ x ̂ t , z t , ψ t ] . The projection x ̂ t is computed using Equation 4 and can be either positive or negative. x ̂ t = ± x t 2 + y t 2 . (4)

As shown in Figure 2, the angle α , about which X t is rotated, is defined by Equation 5. α = atan 2 y t , x t . (5)

The projection does not affect ψ t because ψ t is defined as the rotation about the y ε -axis in the local end-effector frame. In other words, ψ t represents the same local orientation in both the original and projected poses. Similarly, the z-axis also remains unchanged by the projection, i.e., the z-value of the original target pose X t in 3D is identical to the z-value of the projected target pose X ̂ t in 2D.

The projection of different robot configurations onto the global x–z-plane is displayed in Figure 2. X t denotes the original target pose, as defined by the user, and X ̂ t denotes the target pose projected onto the global x–z-plane. The angle of rotation of the pose is denoted by α . For example, the X t , 1 pose is rotated about α 1 = − 60 ° ; the X t , 2 pose is rotated about α 2 = − 155 ° ; and the X t , 3 pose is rotated about α 3 = 90 ° . A priori, it is unknown whether the projection x ̂ t should be positive (e.g., X t , 1 or X t , 3 ) or negative ( X t , 2 ) . A detailed treatment of the projection ambiguity is provided in Section 2.1.

Given a CR or HRR with an elbow joint and a length-adjustable bending part, the robot can achieve target poses X ̂ t = [ x ̂ t , z t , ψ t ] in a 2D plane, where x ̂ t and z t denote the end-effector position and ψ t denotes the end-effector orientation. The elbow joint is capable of deflection angles θ min ≤ θ ≤ θ max , where θ min and θ max are the (mechanical) limits. Figure 3 illustrates when ψ is positive or negative.

Generally, there are three possible robot configurations: downward bending, upward bending, and straight configuration. The distinction between the different cases is relevant as the equations representing the geometric constraints are slightly different for each case. ψ t describes the end-effector orientation of the projected target pose X ̂ t = [ x ̂ t , z t , ψ t ] and must be compared with ψ S , which is computed using Equation 6. ψ S = atan 2 x ̂ t , z t . (6)

In the straight configuration, the distal part l d , the bending part l b , and the end-effector l e are aligned as displayed by robot 3 in Figure 3. It should be noted that in this case, the bending part is straight and not bent, and θ = ψ S = ψ t . The case distinction is carried out according to Equation 7. ψ t < ψ S , downward bending = ψ S , straight configuration  > ψ S , upward bending . (7)

For both downward and upward bending, the set of equations representing the geometric constraints is provided in Equations 8–11. r 2 = a 2 + b 2 , (8) t a n θ = b a , (9) r 2 = c 2 + d 2 , (10) t a n ψ = c d . (11)

In Sections 2.4, 2.5, the placeholders a, b, c, and d are substituted by the arc parameters j = [ x ̂ 0 , z 0 , θ , r ] . x ̂ 0 and z 0 denote the center coordinates of the projected arc in the global x–z-plane, θ describes the deflection angle of the elbow joint, and r is the radius of the arc. After substitution of the placeholders, the set of equations will only depend on known hardware design parameters (i.e., l b and l e ) and the arc parameters j . This allows for the numerical solution of the arc parameters j , which in turn can be used to compute the generalized coordinates q.

For ψ t < ψ S , the robot performs downward bending. Notably, in downward bending, ψ t may be either positive ( 0 ° < ψ t < ψ S ) or negative ( ψ t < 0 ° ) . In Figure 4, the arc parameters j = [ x ̂ 0 , z 0 , θ , r ] are shown for downward bending.

Equations 8–11 can be trivially derived from the Pythagorean and angular relationships, as displayed in Figure 4. For downward bending, the placeholders a, b, c, and d, defined in the set of equations representing the geometric constraints (Equations 8–11), must be substituted as follows: a = − x ̂ 0 + s i n θ ⋅ l d , (12) b = z 0 − c o s θ ⋅ l d , (13) c = z 0 − z t + c o s ψ t ⋅ l e , (14) d = − x ̂ 0 + x ̂ t − s i n ψ t ⋅ l e . (15)

In the substituted set of equations, only the arc parameters j = [ x ̂ 0 , z 0 , θ , r ] are unknown.

It is evident that the same straightforward geometric relationships described in the set of equations also apply to upward bending. For ψ t > ψ S , the robot performs upward bending as shown in Figure 5. For upward bending, the placeholder variables defined in the set of equations representing the geometric constraints (Equations 8–11) must be substituted as follows: a = x ̂ 0 − s i n θ ⋅ l d , (16) b = c o s θ ⋅ l d − z 0 , (17) c = z t − c o s ψ t ⋅ l e − z 0 , (18) d = x ̂ 0 − x ̂ t + s i n ψ t ⋅ l e . (19)

The steps so far can be summarized as follows:

The given target pose X t in 3D is projected to 2D, i.e., X ̂ t . Then, it needs to be determined in which configuration the robot can reach X ̂ t , i.e., downward or upward bending. Depending on the configuration, the placeholders in the set of equations are substituted accordingly. In the substituted set of equations, only the arc parameters j = [ x ̂ 0 , z 0 , θ , r ] are unknown.

As ψ t approaches ψ S , the bending becomes less pronounced, and the radius r increases significantly, as described by Equation 20. ψ t → ψ S ⇒ r → ∞ . (20)

Finding a numerical solution where the variables approach infinity is unfeasible. Thus, to avoid divergence issues, in cases where the robot is straight or almost straight (i.e., | ψ t − ψ S | < 1°), the robot is regarded as being in a straight configuration. For ψ t = ψ S = θ , the robot performs no bending, as visualized in Figure 6; hence, there is no need to find arc parameters j using Powell’s hybrid method.

Instead, the generalized coordinate l b is computed analytically according to Equation 21. l b = x ̂ t − s i n ψ t ⋅ l e + l d 2 + z t − c o s ψ t ⋅ l e + l d 2 . (21)

Thus, in case of a straight configuration, the generalized coordinates are q = [ l b , Δ l b = 0   m m , α , θ = ψ t ] .

Reverting to a linear equation for the (almost) straight configuration ensures an extremely low runtime while maintaining acceptable accuracy. For the upward and downward bending scenarios, the set of equations is solved using Powell’s hybrid method, as described in the next subsection.

Powell’s hybrid method (Powell, 1970), also known as Powell’s dog-leg method, is an iterative optimization algorithm used to solve systems of nonlinear equations. It is designed to find the roots of systems of N nonlinear functions in N variables. In this article, it is used to numerically solve the set of nonlinear equations defined in Equations 8–11, where the placeholder variables a, b, c, and d are defined in Equations 12–15 for downward bending and in Equations 16–19 for upward bending. The objective function f X ̂ t is stated in Equation 22. f X ̂ t j = f X ̂ t x ̂ 0 , z 0 , θ , r = 0 , (22)

where f X ̂ t represents the substituted set of equations and x ̂ 0 , z 0 , θ , r represent the unknown arc parameters j . Now, Powell’s hybrid method is briefly outlined and applied to the equations representing the geometric constraints.

Initialization:

Set the initial guess j 0 = [ ( x ̂ 0 ) 0 , ( z 0 ) 0 , ( θ ) 0 , ( r ) 0 ] .

Iteration: 1)

Evaluate the function f X ̂ t j k and the Jacobian matrix J j k or its approximation B (Chen and Stadtherr, 1981).

The Cauchy point is then found by minimizing the quadratic function along the steepest descent direction as defined in Equation 27: Q μ g = 1 2 D 1 f X ̂ t j k + μ B g 2 . (27)

The multiplier μ is chosen to minimize Q ( μ g ) , which is a quadratic function of μ as stated in Equation 28:: μ = − f X ̂ t j k T B g g T B T B g . (28)

The Cauchy point is then given by p S = μ g . 5)

Determine the dog-leg step:

where α is found such that ‖ p ‖ = Δ . 6)

Update the solution: j k + 1 = j k + p .

This method ensures global convergence by combining the Gauss–Newton direction with the steepest descent direction within a trust region, making it robust for the set of equations representing the geometric constraints.

To ensure that NSGC finds a solution, it is necessary to provide multiple initial guesses x 0 for the algorithm to iterate through. In practice, it is expedient to establish different initial guesses for slight bending (e.g., | ψ t − ψ S | ≤ 15°) and strong bending (e.g., | ψ t − ψ S | > 15°). This approach allows the NSGC algorithm to avoid iterating through all the initial guesses, bypassing those that are potentially ill-suited for the current bending. For very slight bending (e.g., | ψ t − ψ S | < 1°), the robot is considered to be in a straight configuration, as detailed in Section 2.6.

If the projection of the target pose is always performed positively, such that proj ( X t ) = X ̂ t = [ x ̂ t , z t , ψ t ] with x ̂ t = x t 2 + y t 2 according to Equation 4, certain poses become unreachable due to mechanical joint limits. This issue is illustrated in Figure 2 for robot 2; if the projection X ̂ t were carried out positively (i.e., with a rotation of α 2 = 25 ° ) and assuming constant curvature of the bending section, the elbow joint would need to bend downward beyond its lower bound θ < θ min , thus violating the joint constraint. Consequently, the projection must instead be performed negatively with x ̂ t = − x t 2 + y t 2 and α 2 = − 155 ° .

As mentioned in Section 2.1, it is not known a priori whether a given target pose X t should be projected positively or negatively. The limit of the feasible end-effector orientation ψ ̄ t depends on both the target position ( x ̂ t , z t ) and the elbow joint constraint. For the case θ min = 0 ° , the boundary of the reachable orientation is defined using the following equations (Equations 30–32): r = x ̂ 0 , (30) r 2 = r − x ̂ t + s i n ψ ̄ ⋅ l w 2 + z t − l d − c o s ψ ̄ ⋅ l w 2 , (31) ψ ̄ = a t a n 2 z t − l d − c o s ψ ̄ ⋅ l w , r − x ̂ t + s i n ψ ̄ ⋅ l w . (32)

Here, the variables x ̂ 0 , r , and ψ ̄ are unknowns. After solving this system, ψ ̄ is compared to the target end-effector orientation ψ t to determine whether a positive or negative projection is appropriate. If ψ t ≤ ψ ̄ , a positive projection is feasible; otherwise, the projection must be carried out negatively. In the latter case, the robot must rotate about its main axis and approach the target from the − x ̂ direction, thereby staying within the joint constraint θ ≥ θ min = 0 ° and achieving X ̂ t = [ x ̂ t , z t , ψ t ] with ψ t > ψ ̄ , as illustrated in Figure 7a).

This rotation about the main axis significantly enlarges the workspace, as shown in Figure 7b. When a target pose is only reachable via this axial rotation, the projected pose must be transformed as X ̂ t = [ x ̂ t , z t , ψ t ] → X ̃ t = [ − x ̂ t , z t , − ψ t ] . This corresponds to the rotated configuration. In that case, α [Equation 5] has to be increased by 180°. The pose X ̃ t is then used as the target for inverse kinematics and solved using Powell’s hybrid method.

NSGC can provide false positives that are mathematically correct but do not represent the physical robot. The predominant cases of false positives are shown in Figure 8.

It must be ensured that the algorithm does not accept a false positive as a numerical solution. Thus, for each solution is checked against the following three convergence conditions:

0° ! ≤ θ ! ≤ 60°.

It should be noted that for upward bending, condition III is stronger than the statement “if ψ t > ψ S , then θ < ! ψ S .” As shown in Figure 8b, there exists a case where θ < ψ S , but the solution is false. Such a strong formulation is not possible for downward bending, i.e., condition II.

w is a design parameter of the robot that describes the distance between the backbone and the imaginary centerline. This distance is identical to the distance between the connectors and the imaginary centerline. The radius r denotes the radius of the imaginary centerline, as shown in Figure 1.

Using the numerically computed deflection angle θ and radius r , the length of the bending part (i.e., l b ) and the difference between l b and the length of the backbone (i.e., Δ l b ) are computed according to Equations 33, 34 for downward bending. l b = r − w ⋅ | ψ | + θ , (33) Δ l b = r + w ⋅ | ψ | + θ − l b . (34)

l b and Δ l b for upward bending are computed according to Equations 35, 36. l b = r + w ⋅ ψ − θ , (35) Δ l b = r − w ⋅ ψ − θ − l b . (36)

To derive the resulting end-effector pose based on the NSGC solution for the generalized joint coordinates q = [ l b , Δ l b , α , θ ] , q has to be plugged into the forward kinematic model defined by Equation 2. The resulting end-effector pose can then be compared to the target pose, with the target pose X t .

In summary, the final NSGC algorithm is expressed in Algorithm 2. Here, X t denotes the original target pose in 3D, X ̂ t denotes the 2D projection, and X ̃ t denotes a rotation thereof.

Algorithm 2

NSGC.

1: Receive X t ← [ x t , y t , z t , ψ t ] ⊳ target pose in 3D space

The flowchart in Figure 9 visualizes the individual steps of the NSGC algorithm.

The convergence times for all five IKMs are visualized in Figure 10 and summarized in Table 1. The mean convergence times of basic Jacobians and the adaptive damping LM IKMs are very similar to each other (in the range of 21   m s 29   m s ), which is significantly faster than 100   m s needed for real-time applications. NSGC achieves a mean convergence time of 16   m s , which is approximately 24% faster than the fastest competing IKM (i.e., the adaptive damping IKM). As shown in Figure 10, the quadratic programming approach converges with a mean time of 143   m s , which is notably slower than that of the four other models. Therefore, the quadratic programming IKM is not usable for real-time applications.

It has to be noted that the spread of convergence time of NSGC ( 39   m s ) is notably higher than that of all other IKMs, ranging from 6   m s − 28   m s . This large spread of 39   m s can be explained by the necessary iteration through many initial guesses in case the previous initial guess did not lead to a convergence. However, even the longest measured convergence time was only ≈ 40   m s . In contrast, in some cases, the very first initial guess led to convergence, which resulted in an extremely fast convergence time of less than 5   m s , outperforming all other models. IQR stands for interquartile range and measures statistical dispersion, specifically representing the spread of the middle 50% of a dataset.

The convergence behavior of NSGC is displayed in Figure 11. As observed, NSGC iterates through different initial guesses until it finds the correct solution for the set of equations defined in Section 2.3. It is evident that NSGC’s convergence behavior is fundamentally different from that of all other IKM approaches (Figure 11) as there is no convergence in the sense that the error progressively approaches 0. Rather, the convergence of an iteration is completely independent of the previous iteration’s convergence, as it depends solely on the initial guess. Note that the convergence time is not only influenced by the number of iterations until convergence but also by the time needed per iteration. Generally, NSGC and the quadratic programming approach need more time per iteration.

The NSGC algorithm, implemented in Python, exhibits a moderate memory footprint and computational complexity.

The NSGC algorithm demonstrates high robustness in practical inverse kinematic scenarios but retains certain sensitivities that merit discussion.

As CRs/HRRs are used in endoscopy (Sun and Chen, 2024) and biopsy procedures (Gao et al., 2020), this study uses endoscopic examination of the gallbladder as a case study. During a biopsy, the physician traverses the surface of the organ, and once a suitable location is identified, the biopsy needle is inserted into the tissue.

According to Can et al. (2012), the average gallbladder size is h × w × l = 50   m m × 50   m m × 30   m m . Hence, in this study, a cuboid was used to represent the size of the gallbladder. For this simulated case study, NSGC was used to solve the trajectory of the end-effector traversing the cuboid’s surface. The full simulation video is available in the Supplementary Material.

The following section describes the robotic prototype used for the experimental validation.

HyNiTi is an HRR that relies on backbone actuation. The backbone is made of highly elastic metal, such as superelastic NiTi or spring steel. Unlike tendons, the backbone is resistant to stretching and backlash, resulting in more stable and predictable robotic behavior.

Another advantage of this design is the reduced number of channels within the robotic elements. Tendon-driven systems typically require separate channels for each bending direction (e.g., one tendon for upward bending and another for downward bending). In contrast, the pull–push backbone actuation requires only one channel to facilitate movement in two opposing directions, enabling further miniaturization of the robot.

The prototype utilizes an elbow joint that is based on a rack–and-pinion mechanism (Fritsch et al., 2023). The rack is pushed or pulled by a motor (in this case, by motor 3), causing the pinion to rotate about the axis of rotation. The pinion is structurally integrated into the distal part. To mitigate the risk of tilting or jamming of the distal part, the rack and pinion mechanism is designed with redundancy, incorporating a two-sided configuration.

The elbow joint’s deflection angles are in the range of 0° = θ min ≤ θ ≤ θ max = 60°, where θ min and θ max are the mechanical limits.

A total of 100 target poses X t were solved using NSGC, the generalized coordinates q = [ l b , Δ l b , α , θ ] were calculated and sent to the motors, and the achieved pose was measured and compared to the target pose. A picture including the depth information of each pose was captured using the stereo camera and then analyzed. The analysis was based on two points: the proximal end of the wrist and the distal end of the end-effector. From the coordinates of both points, the achieved orientation about the local y ε -axis was derived. The differences between the target and measured poses showed good agreement, as shown in Figure 14 and Table 3. RMSE stands for root mean square error.

As shown in Figure 14, the range of error was the highest for the end-effector angle ψ t , with Δ ψ ranging from − 4 ° to 6° with an RMSE of 2.09°. The RMSE for all three positional coordinates ranged from 1.7   m m along the global z-axis to 3.53   m m along the global x-axis. The mean positional error ranged from 0.42   m m to 2.82   m m . A leading factor contributing to the pose error is the weight of the robot itself, which induces downward bending. The lack of rigidity can be alleviated to some degree by appropriate robot calibration and by increasing the bending section’s stiffness.

For the application of CR/HRR technology in the real world, good trajectory tracking performance is crucial. In this experiment, three trajectories were tracked: a linear 200   m m movement in the x-direction, a linear 200   m m movement in the y-direction, and a movement along a 3D helix. The helix was characterized by a radius of 50   m m and a length of 200   m m along the z-axis, as shown in Figure 15. The error metrics are given in Table 4.

The inverse kinematics for each point on the helix was computed using the proposed NSGC algorithm. The orientation at each point was kept constant; in other words, the end-effector reached each pose on the trajectory with the same orientation. In all three cases, the target trajectory was compared against the measured trajectory. During trajectory tracking, an image was captured and analyzed every 0.5   s . The trajectory tracking validation confirms that NSGC can be used for the control of length-extensible CRs/HRRs with an elbow joint since all poses along the path were successfully reached, as shown in Figure 15.

The average positional error during trajectory tracking was 4.05   m m along the x-axis and 4.22   m m along the y-axis. These errors are slightly higher than those observed in the measurements of individual points (Section 5.1). This discrepancy is most likely caused by the continuous movement of the robot’s end-effector during trajectory tracking, compared to the stationary end-effector used when measuring individual points. With an average positional error of 5.31   m m , trajectory tracking when following the helical path was considerably less accurate than the linear trajectories along the x- and y-axes. This is most likely due to the added dynamics of the end-effector’s movement along the z-direction. The superposition of movements along all three axes decreases the overall accuracy.

The inaccuracies are due to the mechanical structure of the robot: play in the joints connecting the rigid links in the bending section and a relatively heavy end-effector, which causes low-frequency mechanical vibrations.

The HyNiTi prototype, equipped with a biopsy needle and an endoscopic camera (as described in Section 4), was used in this case study to simulate a biopsy on a piece of meat representing the gallbladder. The objective of this case study is to demonstrate NSGC’s effectiveness in real-life applications.

Three distinct target locations were selected on the tissue sample. Position 1 was aligned along the main axis of the proximal region, while positions 2 and 3 were located on the top and bottom surfaces of the tissue, respectively. Reaching these off-axis targets required actuation of the robot’s rotational degree of freedom, denoted as α . As shown in Figure 16, the robot was able to reach all three points on the tissue without any issue.