Over the years, several strategies have been implemented to avoid adverse events within the tissue, such as needle buckling and target displacements (Rossa and Tavakoli, 2017). For example, to prevent needle buckling events as a result of increases in the needle–tissue interaction forces, novel designs for steerable needles were developed. The designs increase the critical buckling force of the needle shaft and eventually reduce insertion forces on the needle during insertions (Webster et al., 2006; Van De Berg et al., 2017). However, these methods do not provide real-time buckling control. It is more important to track the difference between insertion force and critical buckling force relative to the needle position in the tissue, rather than just the critical buckling force (Van De Berg et al., 2017; Hulburt et al., 2019). For this purpose, needle–tissue biomechanics models are used to derive a relationship between insertion forces and needle displacements (Tang et al., 2007; Misra et al., 2010b; Sakes et al., 2016). However, biomechanics models require prior information of tissue properties, which is challenging for intra-operative procedures when tissue environments are unknown. Similarly, imaging techniques (Abolhassani et al., 2006; Hadjerci et al., 2016; Li et al., 2017b), mechanical devices and probes as tissue constraints (Brown et al., 2010; Kobayashi et al., 2012), and path planners (Vrooijink and Abayazid, 2014; Li et al., 2017b; de Baere et al., 2022) have been implemented to control the deformations and displacements of the tissue. Though these methods have shown success in achieving accurate target placement of the needle tip, they rely on finite element and biomehanics models. Accurate characterization of needle–tissue dynamics can be time consuming, and thus, it can be challenging to implement real-time control (Rossa and Tavakoli, 2017; Barnoy et al., 2022).

Recently, research groups have taken inspiration from parasitic invertebrates such as wasp ovipositor and mosquito to design needles that can penetrate solid substrates easily. Inspired by the mechanism of the ovipositor, a thin needle composed of several wires sliding alongside the other was developed to maneuver needles with high curvatures while minimizing impact on the surrounding tissues (Sprang et al., 2016; Scali et al., 2017). Similarly, Li et al. (2020) developed a needle-cannula system that mimics the mosquito proboscis mechanism of incremental motions to reduce tissue deformations and displacements during insertion. They also used harpoon notches at the needle tip to reduce needle–tissue frictional forces (Li et al., 2020). These bio-inspired mechanisms of needle insertions can generally improve maneuverability, while reducing critical needle–tissue interaction events such as needle buckling, tissue deformations, target displacement, and needle deflections. However, they have complicated designs and require several fixtures to prevent buckling of individual wires. Furthermore, a systematic study to select design parameters is required, which usually change due to varying heterogeneous tissue substrate (Leibinger et al., 2016; Matheson and Rodriguez y Baena, 2020).

These limitations call for alternative solutions that do not depend on needle designs and unknown tissue structures, such as model-free techniques. Thus, optimizing on needle insertions while minimizing adverse events in real time using data-driven control techniques would be the focus of this work.

In the last decade, much attention has been given to the application of data-driven control techniques for modeling complex nonlinear processes (Abouaïssa and Chouraqui, 2019). The common classical technique being the proportional, integral, and derivative (PID) controller has shown to be generalizable to many systems and finds applications in most regulatory-based process control industries (Åström and Hägglund, 1995; McMillan, 2012). However, tuning of the gains is performed offline and is affected by process disturbances and the time-varying behavior of the system. This makes PID less robust and challenging to implement in real time. To address these issues, intelligent PID was proposed to predict dynamic system behavior through estimation techniques and online numerical differentiator, eliminating offline tuning and prior identification procedure (Fliess and Join, 2009). Yet, its performance is limited by issues of sensor noise, sampling rates, and computation power as the size of data increases (Hou and Jin, 2011).

There has been increasing interest in unconventional control strategies such as neural networks, fuzzy logic, and genetic algorithms (Thomas and Poongodi, 2009; Heidari et al., 2013; He et al., 2015; Bing et al., 2018; Glida et al., 2021). Neural network control is based on learning from the mapping between system input–output data, whereas fuzzy control is based on learning from past experience and expert knowledge to predict and control system behavior (Abouaïssa and Chouraqui, 2019). Genetic algorithms tune control models based on stochastic optimization techniques that mimic biological evolution (Rodríguez-Abreo et al., 2020). The advantage is that these methods do not require any mathematical model of the system; however, determining parameters that accurately satisfy system behavior can be computationally expensive. Moreover, methods to identify system stability with learning models are still yet to be developed and fully analyzed (Singh and Sukavanam, 2012). Efforts have been made to develop learning-based control techniques with low computation complexity, to optimize needle insertions in tissue while minimizing undesired needle deflections and tissue displacements (Buzurovic et al., 2010; Rossa et al., 2017; Pratt and Petruska, 2022). Yet, their methods do not entail prevention of these events before occurrence.

An alternative technique to control any complex nonlinear system without using learning models is MFAC (Hou and Jin, 2011). This technique uses the dynamic linearization principle to simplify the system nonlinearities and requires only system input–output data. Due to the ability to perform stability analysis online, simplified mathematical models, and robustness against external disturbances, MFAC finds applications in unmanned aerial vehicles, automated manufacturing, and mobile robots (Hou and Jin, 2019). Recently, Li et al. (2017a) derived a model-free adaptive control from Kalman filter techniques to predict needle buckling and control positions of a flexible robot. However, they have not provided a mechanism yet to minimize an adverse event when it occurs. To the best of our knowledge, model-free adaptive control has not yet been implemented on robotic systems that could control needle–tissue interactions while preventing or minimizing adverse events in real time.

If a discrete-time nonlinear system with output y and input u _c satisfies the Lipschitz condition | y k i − y k j | ≤ b ̃ | u c k i − u c k j | (1) for timestamps k _i ≠ k _j and i, j > 0, then the nonlinear system can be approximated as a dynamic linearization model, y k + 1 = y k + Θ c k Δ u c k , (2) where Θ _c is a slow time-varying parameter called PPD that compresses system nonlinearity. This approximation is useful in predicting the output of a nonlinear system without using dynamic models. PPD is a scalar for single-input single-output (SISO) systems and a vector for multi-input multi-output (MIMO) systems. For tracking and regulation-based control, the following quadratic cost function is employed to obtain control effort, Δu _c (k), as a function of prediction errors, y* − y(k + 1): J c u c k = y * − y k + 1 2 + λ c u c k − u c k − 1 , (3) where y* is the desired output of the system. The weighting factor λ _c > 0 restricts Δu _c (k) and, thus, is an important tuning parameter for output convergence and system stability. Substituting y(k + 1) from Eq. 2 and minimizing Eq. 3 with respect to u _c , the optimal input for tracking is computed as u c k = u c k − 1 + Θ c k λ c + Θ c k . Θ c k T y * − y k , (4) where Θ _c (k) is determined using numerical approximation methods due to unknown system dynamics. The following cost function is used to obtain an approximate Θ _c (k) by minimizing the error between the modeled output from Eq. 2 and observed output, y(k): J c Θ c k = y k − y k − 1 − Θ c k Δ u c k − 1 2 + Θ c k − Θ k − 1 2 . (5)

The second term in the cost function ensures robustness against system nonlinearity. Minimizing this cost with respect to Θ _c (k) yields Δ Θ c k = Δ u c k − 1 T . Δ y k − Θ c k − 1 . Δ u c k − 1 1 + Δ u c k − 1 T . Δ u c k − 1 , (6) where Δu _c (k) = u _c (k) − u _c (k − 1). However, Eq. 6 is not robust to process errors such as arbitrary initialization of system parameters, which affect output convergence. This issue is resolved when Θ _c (k) is assigned as a rapid time-varying parameter If  | Δ y k + 1 | > b ̃ | Δ u c k | Θ c k = | Δ u c k | | Δ y k + 1 | . (7)

The Lipschitz constant, b ̃ , is chosen from the theoretical analysis of system stability and output convergence. Finally, Eqs 2–7 are computed iteratively. This technique is similar to the existing framework (Hou and Jin, 2011); however, the novelty lies in the introduction of the rapid varying PPD and the methods of initializing control parameters, (λ _c , Θ _c (0)), and system parameters, (u _c (0), y(0)), for continuous real-time stability discussed in Section 3.2.

In the traditional MFAC scheme, Θ _c (0), y(0), u _c (0), and λ _c are hand-tuned using trial and error methods which can be tedious. This issue can be resolved if a systematic approach to parametric initialization is followed under certain assumptions of the system stated below.

Assumption 1

The nonlinear system (1) satisfies a special case of Lipschitz condition for all k > 0: y k + 1 − y k ≤ b ̃ u c k − u c k − 1 , (8) where u _c (k) ≠ u _c (k − 1) and b ̃ ∈ ( 0,1 ) . Therefore, |Θ _c | ∈ (0, 1).

This assumption holds for a class of systems that cannot move beyond the commanded external force, unless affected by disturbance factors. For example, in robotic needle steering, a tissue cannot displace beyond the needle insertion velocity unless affected by other factors such as patient breathing. This condition resembles the bounded-input bounded-output stability where energy in the system output cannot increase beyond the energy levels in the system inputs according to energy conservation law. Hence, the absolute values of Θ _c from Eq. 2 are also bounded between 0 and 1 according to Eq. 8.

Next, the sign of Θ _c depends on the system response to inputs, which is stated in the following assumption.

Assumption 2

If system output changes in the same direction as the inputs, then Θ _c ∈ (0, 1). Otherwise, Θ _c ∈ (−1, 0).

This assumption holds for systems whose outputs are affected by multiple inputs. For example, if two inputs (say, u _c1 and u _c2) influence an output such that an increase in u _c1 and decrease in u _c2 drive the system to desired outputs, then 0 < Θ _c (0) < 1 for u _c1 and −1 < Θ _c (0) < 0 for u _c2.

Third, if initial values of the system input and output are assigned close to the operating points, then output convergence is guaranteed.

Assumption 3

Output convergence, lim k → ∞ | y ( k + 1 ) | → y * , is guaranteed if −a < y(0) < a and a ̃ < u c ( 0 ) < a ̃ such that |y* − y(0)| → 0 for some a > 0, a ̃ > 0 .

The operating points ( a , a ̃ ) are usually identified during prior simulations. For physical systems, they are easily identified as safe operating limits provided by the system manufacturer.

Finally, output convergence also depends on the range of the input weighting factor, λ _c , stated as follows.

Assumption 4

If system (1) is controlled for a regulation problem, then there exists a minimum value, λ c min = 1 4 such that |y* − y(k + 1)| → 0 for all k and λ c > λ c min . There also exists a maximum value, λ c max , above which the system settles but with steady-state errors.

Proof. The proof is straightforward where the boundedness property of b ̃ from Assumption 1 is used to maintain both system stability and convergence. If the steady-state error is defined as e k + 1 = y * − y k + 1 , (9) then substituting y(k + 1) in Eq. 9 from Eq. 1, we have e k + 1 = y * − y k − b ̃ Δ u c k . (10)

Substituting Δu _c (k) from Eq. 4 in Eq. 10 and from the absolute value theorem, | e k + 1 | ≤ 1 − b ̃ Θ c T λ c + Θ c Θ c T | e k | . (11)

Now, using the boundedness of Θ _c from Assumption 2, we know that 0 < | 1 − b ̃ Θ c T λ c + Θ c Θ c T | < 1 . (12)

Since b ̃ cannot exceed 1 from Eqs 11, 12, we assume | 1 . Θ c T λ c + Θ c Θ c T | < 1 as the maximum limit to find λ c min . Then, using the general concept of quadratic equation, λ c 2 + Θ c Θ c T ≥ 2 λ c Θ c , we have | e k + 1 | ≤ 1 − Θ c T λ c + Θ c Θ c T | e k | ≤ 1 − 1 2 λ c < 1 . (13)

Hence, λ c min = 0.25 using the relation, 1 2 λ c = 1 in Eq. 13. Finding an exact maximum value, λ c max , is not important, as λ _c > 1 leads to lower rising times and larger steady-state errors for many systems (Hou and Jin, 2011). Since many real-time applications require fast rising time and low steady-state errors, λ c max = 1 is considered. Furthermore, Δu _c (k) needs to be constrained as most actuators have finite excitation. Thus, λ _c can be chosen arbitrarily in the range ( λ c min , 1 ) to balance between constrained inputs and tracking errors.

The improved framework (IMFAC) is validated on two test nonlinear systems to compare performance with MFAC (Hou and Jin, 2011). All equations are coded in the MATLAB^® 2018Rb environment and tested on simulated SISO and MIMO nonlinear systems. The algorithms are compared using the same initial conditions and control parameters.

For a SISO system with some nonlinear input–output model, y k + 1 = y k 1 + y k 2 + u c 3 k ; k ≤ 500 y k y k − 1 y k − 2 u c k − 1 y k − 2 − 1 + β k u c k 1 + y k − 1 2 + y k − 2 2 ; k > 500 , β k = round k / 500 . (14)

Let the initial conditions of the system be u _c (1) = u _c (2) = 0 and y(1) = y(2) = 0 with a constant reference output as y* = 0.5. The output of the system is expected to change at k = 500, 750 according to Eq. 14, and thus, some disturbances due to these transitions are expected. Simulation results of SISO are shown in Figures 2A, B, where solid green plots correspond to IMFAC, dashed red plots correspond to original MFAC, and the dash–dot blue plot is the reference signal. For λ c = 0.1 < λ c min , IMFAC converged to the desired value within k = 790 after disturbance at k = 750, while MFAC showed lack of output convergence due to sustained oscillations (upper plot of Figure 2A). When |Θ _c | > 1, IMFAC showed faster convergence at the disturbances, whereas MFAC showed slower convergence at k = 830 after disturbance at k = 750 (upper plot of Figure 2B). For the same control parameters and initial conditions, IMFAC showed better stability and faster convergence due to the corrective action of the rapid time-varying PPD, compared to MFAC with slow time-varying PPD (lower plots of Figures 2A, B).

We consider a two-input (u _c1, u _c2) one-output (y) MIMO system with some nonlinear model in Eq. 15. y k + 1 = 5 y k + 2 u c 1 k − 3 u c 2 k 2 + 2 u c 1 k 2 5 + u c 1 k + 5 u c 2 k , (15) where u c = [ u c 1 , u c 2 ] T ∈ R 2 × 1 . Arbitrary initial conditions of the system are u _c (1) = u _c (2) = [0.5,0.6] ^T , y(1) = 0, y(2) = 0.9, and Θ _c (1) = [0.5, −0.2]. Let λ _c1 = λ _c2 = λ _c = 0.5 for each input. We can have λ _c1 ≠ λ _c2; however, for the ease of validating Assumption 2, we choose λ _c1 = λ _c2. A constant reference trajectory close to zero (y* = 0.01) is considered since ideal y* = 0 would cause control inputs to be zero which cannot drive the system to desired outputs. Therefore, near zero reference is used to simulate stationary or minimum movement of the system. Simulation results show that negative Θ _c2 = −0.2 decreases u _c2 and positive Θ _c1 = 0.5 increases u _c1, which drive the system to the desired output (y* = 0.01), see Figure 2C. It should be noted that initial Θ _c can be assigned to control the system inputs in the desired direction only if there is basic knowledge about system response to the inputs.

Finally, the performance of IMFAC is compared with original MFAC in terms of frequency of unstable events computed as the number of occurrences of |Δy(k + 1)| > |Δu _c (k)| per simulation time, for different values of control parameters (λ _c and Θ _c ) as shown in Figure 2D. The blue diamond plot (IMFAC) for λ = 0.1 < λ c min shows that system stability was not affected by arbitrary initialization of λ _c and other disturbances. Otherwise, these results show the comparable performance of IMFAC with MFAC using a systematic approach to parameter initialization. On average, the convergence time of IMFAC was 20% faster than that of MFAC. This shows that irrespective of the control parameters and initial conditions, system stability and fast output convergence are always achieved with IMFAC.

The MIMO system in Section 3.3 is similar to the needle intervention applications where tissue motion, y = l, is controlled by two inputs: a) insertion velocity ν and b) duty-cycle factor as a measure of the number of needle shaft rotations per insertion, D. Tissue motion can be minimized with a decrease in insertion velocity and increase in the duty-cycle factor. Therefore, stationary tissue is desired (l* → 0.01) during needle insertion procedures. Through appropriate selection of parameters (λ _c , Θ _c ), IMFAC can minimize tissue motions during needle insertions as shown in Figure 3. The input–output equations of the SISO and MIMO systems from Section 3.3 are used only for simulation purpose to validate the control algorithms and are not available in real time. For real-time needle insertion experiments in unknown tissues, l is obtained from position sensor measurements, and inputs are the user-defined steering control inputs, ν and D. Based on the previous MIMO simulation results, initial values of Θ _c = [0.5,−0.2] ^T and same λ _c = 0.5 will be used for both the inputs (ν and D) for all needle insertions. Initial values of the inputs are commanded by a user. The range or operating points of these inputs is decided by the maximum rating of the motors used to drive the needle following Assumption 3. Since the output (tissue displacement) is zero when excitation (insertion inputs) is zero for needle–tissue systems, we will assume zero initial values, l(0) = 0, for all needle insertions. The next section will discuss application of IMFAC to prevent adverse events during needle insertions in a tissue.

Let the time series data obtained from a sensor at current time k be S(k). The steps to generate sensor data forecasts, S ̂ N ( k + 1 ) , for some future prediction horizon, N, using CFDL-MFP are S ̂ N k + 1 = S ̂ N k + Φ ̂ N Δ U N ∈ R N . (16)

The data vector S represents sensor data. For example, previous N samples of force data measured from a force sensor at time instant k are represented as S _N (k), and S ̂ N ( k + 1 ) represents N-step ahead forecasts of force data at the next time instant k + 1. The forecast parameters ΔU _N (k) and Φ ̂ N ( k ) are computed following the procedure in our previous work (Narayan and Fey, 2020): Φ ̂ N k = Δ U N k − 1 + 1 ⃗ ◦ − 1 ◦ Φ ̂ N k − 1 + Δ U N k − 1 ◦ Δ S ̂ N k Δ U N k = Φ ̂ N k λ . 1 ⃗ + 2 Φ ̂ N 2 k ◦ S N * k + S N k − 2 S ̂ N k , (17) where ◦ is array element-wise operation, 1 ⃗ is an N-dimension vector of ones, and Δ S ̂ ( k ) = S ̂ ( k ) − S ̂ ( k − 1 ) . The sensor data S _N (k) contain previous N observed data points, and the reference vector S N * ( k ) is the running average of the observed data points.

Forecasts of force sensor data are generated when a slight increase in force is detected, and if force forecasts show patterns of rapid increase, then needle buckling is predicted (Narayan and Fey, 2020). Similarly, tissue displacements can be predicted using forecasts of tissue position data. Tissue displacement is characterized based on the intensity of motion such as mild tissue displacement (MTD) and severe tissue displacement (STD). STD occurs when the needle tip stops moving relative to the tissue frame (Narayan et al., 2018a). To predict STD before occurrence, N-step forecasts of the tissue position relative to fixed frame l ̂ twN and the needle-tip position relative to fixed frame l ̂ nwN are generated using Eqs 16, 17. It should be noted that l _tw and l _nw are the Euclidean norm of 3D positions measured from position sensors p ⃗ t w and p ⃗ n w , respectively. The norm of the needle-tip position relative to the tissue frame is obtained as l _nt = l _nw − l _tw , and ( l ̂ twN , l ̂ nwN , l ̂ ntN ) are the forecasts of the corresponding position norms. STD events are then predicted using these norm position forecasts with the sigmoid slope detection methods in Eq. 18. m ̂ t k + i = 180 π . Δ l ̂ N k + i Δ s ̂ k + i ; i = 1,2 , . , N ; l ̂ N = l ̂ twN , l ̂ nwN , l ̂ ntN g ̂ m ̂ t k + i = 1 1 + e − m ̂ t k + i ; g ̂ = g ̂ n w , g ̂ t w , g ̂ n t , (18) If  g ̂ t w k + i > m 2  OR  g ̂ t w k + i < m 1  Then  If  m 1 < g ̂ n t k + i < m 2 ;  STD predicted Else  MTD predicted .

The needle insertion distance s ̂ is linearly projected into the prediction horizon N. The thresholds m ₁ and m ₂ are experimentally determined (Narayan et al., 2018b).

To prevent critical events such as needle buckling and STD before they occur, it is necessary for the robot to respond immediately to the prediction algorithms. For this purpose, reactive control is implemented where the robot stops and retracts the needle when buckling is predicted, and the robot stops insertion when STD is predicted to avoid potential damage to the tissue. However, when a mild event such as MTD is detected, adaptive control (IMFAC) computes optimal steering inputs (insertion velocity and duty cycle) to steer needles such that the occurrence of mild events is minimized and severe events are prevented. This strategy ensures continuity of the insertion procedure which otherwise is disrupted by reactive control. A complete event prevention algorithm with a combination of reactive and adaptive control is shown in Figure 4.

Algorithm 1

Adverse event prevention.

Require: Axial force sensor f _a (k), 3D needle-tip sensor p ⃗ n w ( k ) , 3D tissue sensor p ⃗ t w ( k ) , prediction horizon N, total insertion distance s, commanded insertion velocity ν, and commanded duty-cycle factor D.

The pseudo code for adverse event prevention is presented in Algorithm 1. The code is an integration of the proposed methods with the previous event detection, sensor forecasting (CFDL-MFP), and prediction algorithms (Narayan et al., 2018a; Narayan and Fey, 2020). The inputs to the algorithm are the filtered sensor data at current time k, CFDL-MFP prediction horizon N, user-commanded initial needle insertion velocity ν, and duty-cycle factor D. Before iteration, the IMFAC control parameters are initialized according to Section 3.2, and the initial event prediction flag is set to zero, see steps 3–4. The total insertion distance s is user-defined prior to the insertion procedure. Encoder readings s(k) are used to check if the measured insertion distance reached user-defined s, see step 5. All sensor data are filtered in real time using the Kalman filtering scheme from the work of Narayan et al. (2018b). For every iteration, the Euclidean norm of 3D position sensors and optimal mean forces f p * ( k ) are computed, since CFDL-MFP and detection algorithms work on univariate time series data, see steps 6–7. The Optimal Mean function is the running average of the observed sensor data (Narayan et al., 2018a). Checks for needle buckling detection are performed for every k, see step 8. If a force increase is alerted by the detection algorithm, N-step forecasts of force data f ̂ a N ( k ) are generated with CFDL-MFP and the buckling prediction algorithm checks for any buckling event using the force sensor forecasts (Narayan and Fey, 2020). If buckling is predicted, the robot stops insertions and needle is retracted, see steps 9–16. Similarly, checks for tissue displacement detection are performed, see step 17. If MTD is detected (k = k _m ), the displacement prediction algorithm checks for STD events. When STD is predicted, the robot stops insertions, see steps 18–24. Otherwise, steering control inputs for the needle are continuously updated until the end of insertion distance s to minimize or prevent further displacements of the tissue, see steps 25–29. The algorithm is executed in a multi-threaded framework to ensure simultaneous prediction and prevention of adverse events in real time.

The needle steering robot for needle insertion experiments is shown in Figure 5A. The robot is controlled with a Dell Inspiron 15 Gaming Laptop through custom C++ multi-threaded application (QT Creator 4.6.0) on Ubuntu Linux 16.04. The code manages low-level PID motor control and the execution of duty-cycled flipping algorithms to enable different needle steering directions (Majewicz et al., 2014). The robot is equipped with a 6-axis force–torque sensor (Nano-17, ATI Inc.) at the needle base to track axial forces due to needle–tissue interaction and an electromagnetic (EM) tracker (trakSTAR, NDI Inc.) to track needle-tip position. Data from these sensors are streamed at 50 Hz similar to (Narayan et al., 2018b). Algorithm 1 is coded as modules in Python 3.5 and interfaced with a low-level motor control code in QT Creator C++ through pybind11, an open-source Python-C++ wrapper that calls Python modules from C++ directly enabling multi-threaded operations. Real-time plotting of sensor data and prediction results of Algorithm 1 during needle insertions are displayed on QT GUI using the QCustomPlot library.

A homogeneous gelatin tissue of stiffness 60 kPa (1:3 gelatin to sugar in 500 mL water) and dimensions 14 × 8 × 7.5 cm was created to validate tissue displacement prevention methods from Algorithm 1. The tissue is placed on a platform smeared with soap water to simulate a slippery surface for the sliding motion of the tissue, as shown in Figure 5B. For short-hand notation, let this sample be denoted as T ₁. Another gelatin sample of varying stiffness 60 and 20 kPa (1:3 gelatin to sugar in 1000 mL water) was created by embedding an obstacle to validate needle buckling prevention methods in heterogeneous tissue with a known environment. To prevent damage to the needle and the EM needle position sensor, a standard rectangular Scotch-Brite dish scrubber of thickness 0.7 cm is used as an obstacle. The 20 kPa region is on one side of the obstacle at 7.1 cm from a short edge, and the stiffer 60 kPa region is on the other side of the obstacle at 6 cm from the opposite short edge as shown in Figure 5C. This tissue denoted as T ₂ is also placed on the slippery platform to validate Algorithm 1 during buckling and tissue displacement events. For more details on gelatin preparation methods, refer to the work of Narayan et al. (2018a). Both T ₁ and T ₂ are used for experiments within 12 h of refrigeration. A trakSTAR sensor was taped on the surface of these tissues to track displacements.

Two steerable needles with different dimensions and tips are used to validate the proposed prevention methods for robustness to needle types. The first steerable needle (N ₁) shown in Figures 5D is made from a hollow nitinol tube with dimensions 250 mm length, 0.8 mm OD, and 0.6 mm ID and a 30° prebent tip of 4 mm length (Majewicz et al., 2014). For integration with the robot, an EM needle tracker was inserted into the hollow end of the needle and secured with Luer locks. The second needle (N ₂) is a clinically available prostrate seeding needle with a stiff hollow tubing of length 175 mm and thickness 1.27 mm. A stainless steel wire of thickness 1.2 mm was cut into a small piece of length 10 mm and inserted into the tip of the hollow seeding tube such that 5 mm was jutting out of the tube and the rest 5 mm was inside the tube. The tip was ground to create a bevel angle of 60° to enable cutting through the tissue. Another EM needle tracker was inserted into the hollow end of the needle and secured with Luer locks for integration with the robot similar to N ₁, as shown in Figure 5E.

To specifically evaluate IMFAC control and tissue displacement prevention methods against different needle types, a total of 30 insertions with prebent needle N ₁ and six insertions with straight needle N ₂ were performed in tissue T ₁. To validate both buckling prevention and tissue displacement prevention algorithms against various needle and tissue types, a total of 10 insertions with needle N ₁ and four insertions with N ₂ were performed in tissue T ₂. These insertions were performed using four initial user-commanded insertion velocities (0.25,0.5,1,1.5)cm/s, and three initial user-commanded duty cycle (0,50,100)%. Videos of the insertions were recorded using an iPhone 11 camera focused on the top surface of the tissue and a digital camera (Sony W800/B20MP 5X Optical Zoom) focused on the overall robot-tissue setup. The videos were used as ground truth for verification of the prevention algorithms. For viewing and analysis, motion images were generated in MATLAB at 30 frames per second and overlaid using Adobe Photoshop CC 2020.

Around 16 needle insertions were performed with the ex vivo goat liver tissue labeled as T ₃ in Figure 5F at the University of Texas Southwestern Medical Center. The liver tissue was placed on a clean platform. Live fluoroscopic videos were collected using a 9-in XRII OEC series 9600 fluoroscope (X-ray source and CCD camera) and saved with a Canopus ADVC 110 frame grabber, as shown in Figure 5G. The fluoroscopic videos were used as ground truth to check if adverse events were detected or prevented using Algorithm 1.